Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2927818

Formats

Article sections

- Abstract
- 1 Introduction
- 2 Problem Definition and Probabilistic Model
- 3 Probabilistic View of the Level Set Framework
- 4 Alternating Minimization Algorithm
- 5 Experimental Results
- 6 Discussion and Conclusions
- References

Authors

Related links

Med Image Comput Comput Assist Interv. Author manuscript; available in PMC 2010 August 25.

Published in final edited form as:

Med Image Comput Comput Assist Interv. 2009; 12(Pt 1): 272–280.

PMCID: PMC2927818

NIHMSID: NIHMS179638

Tammy Riklin Raviv,^{1} Koen Van Leemput,^{1,}^{2,}^{3} William M. Wells, III,^{1,}^{4} and Polina Golland^{1}

See other articles in PMC that cite the published article.

Spatial priors, such as probabilistic atlases, play an important role in MRI segmentation. However, the availability of comprehensive, reliable and suitable manual segmentations for atlas construction is limited. We therefore propose a joint segmentation of corresponding, aligned structures in the entire population that does not require a probability atlas. Instead, a latent atlas, initialized by a single manual segmentation, is inferred from the evolving segmentations of the ensemble. The proposed method is based on probabilistic principles but is solved using partial differential equations (PDEs) and energy minimization criteria. We evaluate the method by segmenting 50 brain MR volumes. Segmentation accuracy for cortical and subcortical structures approaches the quality of state-of-the-art atlas-based segmentation results, suggesting that the *latent atlas* method is a reasonable alternative when existing atlases are not compatible with the data to be processed.

Probabilistic atlases are crucial for most MR segmentation methods due to the absence of well defined boundaries. Derived from comprehensive sets of manually labeled examples, atlases provide statistical priors for tissue classification and structure segmentation [1,9,13,14,18]. Although atlas-based segmentation methods often achieve accurate results, the need for spatial priors can be problematic. First, the availability of suitable atlases is limited since manual segmentation of a significant number of volumes requires expensive effort of an experienced physician. Second, the suitability of existing atlases for images from different populations is questionable. Examples include using normal adult brain atlas for brain parcellation of young children or patients with severe brain pathologies.

Recently, a few methods have been proposed to reduce or avoid the dependency on possibly incompatible atlases. In the atlas-based segmentation method in [3], topological constraints are used to avoid possible bias introduced by the atlas. In [19], manually labeled structures are used to support the automatic segmentation of neighboring structures within the same image. Tu *et al.* [17] propose a discriminative approach for the segmentation of adjacent brain structures using a set of discriminative features learned from training examples. Lord *et al.* [11] suggest a group-wise smoothing, segmentation and registration method for cross sectional MR scans. Bhatia *et al.* [5] update an initial atlas constructed from a different population using the evolving segmentations of multiple images.

Here we propose and demonstrate a method that does not use a set of training images or probabilistic atlases as priors. Instead we extract an ensemble of corresponding structures simultaneously. The evolving segmentation of the entire image set supports each of the individual segmentations. In practice, a subset of the model parameters, called the spatial parameters, is inferred as part of the joint segmentation processes. These latent spatial parameters, which can be viewed as a ‘dynamic atlas’, are estimated exclusively from the data at hand and a single manual segmentation. The latent atlas is used as a Markov Random Field (MRF) prior on the tissue labels. The main novelty of the method with respect to other group-wise segmentation methods such as [11,5] is the consistent statistically-driven variational framework for MR ensemble segmentation.

Our contribution is two-fold. We introduce a level set framework, that is based on probabilistic principles, in which segmentation uncertainty is expressed by the logistic function of the associated level set values, similar to [15]. We then use it for group-wise segmentation. We evaluate our method by segmenting the amygdala, temporal gyrus and hippocampus in each hemisphere in 50 MR brain scans. The dice scores achieved by our method approach the atlas-based segmentation results of [14].

Our objective is to segment a particular structure or region of interest in *N* aligned MR images. Specifically, we consider the 2-partition problem where each voxel in image *I _{n}* (

Let each image *I _{n}*:

Let {*I*_{1} … *I _{N}* } be the given set of aligned images that form the observed variable in our problem and let

We jointly optimize for the segmentations *Γ* and the parameters *Θ*, assuming that *I*_{1} … *I _{N}* are independent:

$$\{\widehat{\mathrm{\Theta}},\widehat{\mathrm{\Gamma}}\}=arg\underset{\{\mathrm{\Theta},\mathrm{\Gamma}\}}{max}logp({I}_{1}\dots {I}_{N},{\mathrm{\Gamma}}_{1}\dots {\mathrm{\Gamma}}_{N};\mathrm{\Theta})$$

(1)

$$=arg\underset{\{\mathrm{\Theta},\mathrm{\Gamma}\}}{max}\sum _{n=1}^{N}logp({I}_{n},{\mathrm{\Gamma}}_{n};\mathrm{\Theta})$$

(2)

$$=arg\underset{\{\mathrm{\Theta},\mathrm{\Gamma}\}}{max}\sum _{n=1}^{N}[logp({I}_{n}\mid {\mathrm{\Gamma}}_{n};{\theta}_{I,n})+logp({\mathrm{\Gamma}}_{n};{\theta}_{\mathrm{\Gamma}})].$$

(3)

We propose to alternate between estimating the maximum a posteriori (MAP) segmentations and updating the model parameters. For a given value of the model parameters , Equation (3) implies that the segmentations can be estimated by solving *N* separate MAP problems:

$${\widehat{\mathrm{\Gamma}}}_{n}=arg\underset{{\mathrm{\Gamma}}_{n}}{max}[logp({I}_{n}\mid {\mathrm{\Gamma}}_{n};{\theta}_{I,n})+logp({\mathrm{\Gamma}}_{n};{\theta}_{\mathrm{\Gamma}})].$$

(4)

We then fix and estimate the model parameters *Θ* = {*θ _{Γ}*,

$${\widehat{\theta}}_{I,n}=arg\underset{{\theta}_{I,n}}{max}logp({I}_{n}\mid {\mathrm{\Gamma}}_{n},{\theta}_{I,n})$$

(5)

$${\widehat{\theta}}_{\mathrm{\Gamma}}=arg\underset{{\theta}_{\mathrm{\Gamma}}}{max}\sum _{n=1}^{N}logp({\mathrm{\Gamma}}_{n}\mid {\theta}_{\mathrm{\Gamma}}).$$

(6)

In the following sections we present a level-set framework that is motivated by this probabilistic model. We reformulate the estimation problem stated in Eq. (4) such that the soft segmentations *p*(*Γ _{n}*) rather then the

Now we draw the connection between the probabilistic model presented above and the level set framework for segmentation. Let * _{n}*:

$${\stackrel{\sim}{H}}_{\epsilon}({\phi}_{n})=\frac{1}{2}\left(1+tanh\left(\frac{{\phi}_{n}}{2\epsilon}\right)\right)=\frac{1}{1+{e}^{-{\phi}_{n}/\epsilon}}.$$

(7)

For *ε* = 1, the function * _{ε}*(·) is the logistic function. Similar to [15], we define the level set function

$${\phi}_{n}(\mathbf{x})\triangleq \epsilon \phantom{\rule{0.16667em}{0ex}}\text{logit}(p)=\epsilon log\frac{p(\mathbf{x}\in w)}{1-p(\mathbf{x}\in \omega )}=\epsilon log\frac{p(\mathbf{x}\in \omega )}{p(\mathbf{x}\in \mathrm{\Omega}\backslash \omega )}.$$

(8)

The scalar *ε* determines the scaling of the level set function * _{n}* with respect to the ratio of the probabilities. Substituting this definition into Eq. (7) we obtain

$${\stackrel{\sim}{H}}_{\epsilon}({\phi}_{n}(\mathbf{x}))=\frac{1}{1+p(\mathbf{x}\in \mathrm{\Omega}\backslash \omega )/p(\mathbf{x}\in \omega )}=p(\mathbf{x}\in \omega ),$$

(9)

which implies that the function * _{ε}*(

In the following subsections we relate the terms in Eq. (3) to the energy terms in the classical level set functional.

Let us first consider the image likelihood term in Eq. (3):

$$logp({I}_{n}\mid {\mathrm{\Gamma}}_{n},{\widehat{\theta}}_{I,n})=\sum _{\{v\mid {\mathrm{\Gamma}}_{n}^{v}=1\}}log{p}_{\text{in}}({I}_{n}^{v};{\theta}_{I,n})+\sum _{\{v\mid {\mathrm{\Gamma}}_{n}^{v}=0\}}log{p}_{\text{out}}({I}_{n}^{v};{\theta}_{I,n}),$$

(10)

where *p*_{in} and *p*_{out} are the probability distributions of the foreground and back-ground image intensities, respectively.

Let *E _{I}* −log

$${E}_{I}({\phi}_{n},\mathrm{\Theta})=-{\int}_{\mathrm{\Omega}}\left[log{p}_{\text{in}}({I}_{n};{\theta}_{I,n})\stackrel{\sim}{H}({\phi}_{n}(\mathbf{x}))+log{p}_{\text{out}}({I}_{n};{\theta}_{I,n})\stackrel{\sim}{H}(-{\phi}_{n}(\mathbf{x}))\right]d\mathbf{x}.$$

(11)

If we use, for example, Gaussian densities for *p*_{in} and *p*_{out} we get the familiar minimal variance term [6,12]. Here, we use a Gaussian mixture to model the background, as described later in the paper.

We define the prior probability *p*(*Γ _{n}*|

$$p({\mathrm{\Gamma}}_{n}\mid {\theta}_{\mathrm{\Gamma}})=\frac{1}{Z({\theta}_{\mathrm{\Gamma}})}\prod _{v=1}^{V}{({\theta}_{\mathrm{\Gamma}}^{v})}^{{\mathrm{\Gamma}}_{n}^{v}}{(1-{\theta}_{\mathrm{\Gamma}}^{v})}^{(1-{\mathrm{\Gamma}}_{n}^{v})}{e}^{-f({\mathrm{\Gamma}}_{n}^{v},{\mathrm{\Gamma}}_{n}^{\mathcal{N}(v)})},$$

(12)

where *Z*(*θ _{Γ}*) is the partition function and (

Using the continuous level set formulation with soft segmentation, we define the spatial energy term. as follows:

$${E}_{S}({\phi}_{n},\mathrm{\Theta})=-{\int}_{\mathrm{\Omega}}\left[log{\theta}_{\mathrm{\Gamma}}(\mathbf{x})\stackrel{\sim}{H}({\phi}_{n}(\mathbf{x}))+log(1-{\theta}_{\mathrm{\Gamma}}(\mathbf{x}))\stackrel{\sim}{H}(-{\phi}_{n}(\mathbf{x}))\right]d\mathbf{x}$$

(13)

As in [4] we ignore the partition function, approximating the MRF model above. The logarithm of the pairwise clique potential term *f* (·) can be configured to act as a finite difference operator approximating the gradient of *Γ _{n}* at the voxel

$${E}_{\text{LEN}}({\phi}_{n})={\int}_{\mathrm{\Omega}}\mid \nabla \stackrel{\sim}{H}({\phi}_{n}(\mathbf{x}))\mid d\mathbf{x},$$

(14)

which is the commonly used smoothness constraint as reformulated in [6].

We now construct the cost functional for _{1} … * _{N}* and the parameters

$$E({\phi}_{1}\dots {\phi}_{N},\mathrm{\Theta})=\gamma {E}_{\text{LEN}}+\beta {E}_{I}+\alpha {E}_{S},$$

(15)

where *α* = 1 − *β* − *γ*. As in [16] we tune the weights such that the contributions of the energy terms *E*_{LEN}, *E _{I}* and

We optimize the unified functional (15) in a set of alternating steps. For fixed model parameters *Θ*, the evolution of each of the level set functions * _{n}* follows the gradient descent equation:

$$\frac{\partial {\phi}_{n}}{\partial t}=\stackrel{\sim}{\delta}({\phi}_{n})\left\{\gamma \text{div}(\frac{\nabla {\phi}_{n}}{\mid \nabla {\phi}_{n}\mid})+\beta \left[log{p}_{\text{in}}({I}_{n}(\mathbf{x});{\widehat{\theta}}_{I,n})-log{p}_{\text{out}}({I}_{n}(\mathbf{x});{\widehat{\theta}}_{I,n})\right]+\alpha \left[log{\widehat{\theta}}_{\mathrm{\Gamma}}-log(1-{\widehat{\theta}}_{\mathrm{\Gamma}})\right]\right\},$$

(16)

where $\stackrel{\sim}{\delta}({\phi}_{n})\triangleq {\stackrel{\sim}{\delta}}_{\epsilon}({\phi}_{n})={\scriptstyle \frac{d{\stackrel{\sim}{H}}_{\epsilon}({\phi}_{n})}{d{\phi}_{n}}}$ is the derivative of the Heaviside function, i.e.,

$${\stackrel{\sim}{\delta}}_{\epsilon}({\phi}_{n})=\frac{1}{2\epsilon}\text{sech}(\frac{{\phi}_{n}}{2\epsilon})=\frac{1}{\epsilon cosh({\scriptstyle \frac{{\phi}_{n}}{\epsilon}})}.$$

For fixed segmentations * _{n}* the model parameters are recovered by differentiating the cost functional (15) with respect to each parameter.

We assume that the intensities of the structure of interest are drawn from a normal distribution, i.e.,
${p}_{\text{in}}({I}_{n};{\theta}_{I,n})=\mathcal{N}({I}_{n};{\mu}_{n},{\sigma}_{n}^{2})$. The distribution mean and variance of the foreground region of *I _{n}* are updated at every iteration. The intensities of the background tissues are modeled as a K-Gaussian mixture:

$${p}_{\text{out}}({I}_{n};{\theta}_{I,n})=\text{GMM}({\mu}_{n}^{1}\cdots {\mu}_{n}^{K},{\sigma}_{n}^{1}\cdots {\sigma}_{n}^{K},{w}_{n}^{1}\cdots {w}_{n}^{K}),$$

where
${w}_{n}^{k}$ is the mixing proportion component *k* in the mixture. The Gaussian mixture model parameters are estimated using expectation maximization (EM) [7].

We estimate the spatial function *θ _{Γ}*(

$${\widehat{\theta}}_{\mathrm{\Gamma}}=arg\underset{{\theta}_{\mathrm{\Gamma}}}{max}\sum _{n=1}^{N}{\int}_{\mathrm{\Omega}}[\stackrel{\sim}{H}({\phi}_{n}(\mathbf{x}))log({\theta}_{\mathrm{\Gamma}}(\mathbf{x}))+(1-\stackrel{\sim}{H}({\phi}_{n}(\mathbf{x})))log(1-{\theta}_{\mathrm{\Gamma}}(\mathbf{x}))]d\mathbf{x},$$

yielding ${\widehat{\theta}}_{\mathrm{\Gamma}}(\mathbf{x})={\scriptstyle \frac{1}{N}}{\sum}_{n=1}^{N}\stackrel{\sim}{H}({\phi}_{n}(\mathbf{x}))$.

We test the proposed approach on 50 MR brain scans. Some of the subjects in this set are diagnosed with the first episode schizophrenia or affective disorder. The MR images (T1, 256 × 256 × 128 volume, 0.9375 × 0.9375 × 1.5 mm^{3} voxel size) were acquired by a 1.5-T General Electric Scanner. The data was originally acquired for brain morphometry study [8]. In addition to the MR volumes, manual segmentations of three structures (superior temporal gyrus, amygdala, and hippocampus) in each hemisphere were provided for each of the 50 individuals and used to evaluate the quality of the automatic segmentation results. MR images are preprocessed by skull stripping. The volumes were aligned using B-spline registration according to [2].

Assuming that the manual segmentation of a single instance is given, we initialize the latent atlas *θ _{Γ}* by the Heaviside function of a single manual segmentation smoothed with a Gaussian kernel of width

We used the Dice score to evaluate segmentations. The results are shown in Figure 1. In contrast to [14] we do not use spatial priors, neither do we use hierarchical multi-stage segmentation model. Nevertheless, the Dice scores obtained by our method approache the state-of-the-art atlas-based segmentation results reported there. To exemplify the significance of a latent atlas, generated concurrently with the segmentation, we also show a comparison to the segmentation results obtained by using an atlas constructed from a single manual segmentation smoothed by a Gaussian kernel, without the update procedure. Figure 2 shows segmentation examples of the three pairs of structures in representative individual brains. Coronal views of the resulting 3D atlases for each pair of structures are shown in the fourth column of Figure 2.

We presented a level set framework for segmentation of MR image ensembles, that is motivated by a generative probabilistic model. Unlike most previous methods, we do not use spatial priors in the form of a probabilistic atlas. Instead, spatial latent parameters, which form a ‘dynamic atlas’, are inferred from the data set through an alternating minimization procedure.

The quality of the segmentation results obtained for ensembles of brain structures shows that the proposed method presents a reasonable alternative to standard segmentation techniques when a compatible atlas is not available.

An on-going research is now conducted to demonstrate the ability of the proposed algorithm to handle pathological cases (e.g., in a longitudinal or a multimodal study) where the atlas-based approach still fails.

This work was supported in part by NIH NIBIB NAMIC U54-EB005149, NIH NCRR NAC P41-RR13218, NIH NINDS R01-NS051826, NIH NCRR mBIRN U24-RR021382 grants and NSF CAREER Award 0642971.

1. Ashburner J, Friston K. Unified segmentation. NeuroImage. 2005;26:839–851. [PubMed]

2. Balci SK, Golland P, Shenton M, Wells WM. Free-form B-spline deformation model for groupwise registration. Proc. of MICCAI 2007 Statistical Registration Workshop: Pair-wise and Group-wise Alignment and Atlas Formation; 2007. pp. 23–30. [PMC free article] [PubMed]

3. Bazin P-L, Pham DL. Statistical and topological atlas based brain image segmentation. In: Ayache N, Ourselin S, Maeder A, editors. MICCAI 2007, Part I. LNCS. Vol. 4791. Springer; Heidelberg: 2007. pp. 94–101. [PubMed]

4. Besag J. Statistical analysis of non-lattice data. The Statistician. 1975;24(3):179–195.

5. Bhatia KK, Aljabar P, Boardman JP, Srinivasan L, Murgasova M, Counsell SJ, Rutherford MA, Hajnal JV, Edwards AD, Rueckert D. Groupwise combined segmentation and registration for atlas construction. In: Ayache N, Ourselin S, Maeder A, editors. MICCAI 2007, Part I. LNCS. Vol. 4791. Springer; Heidelberg: 2007. pp. 532–540. [PubMed]

6. Chan T, Vese L. Active contours without edges. TIP. 2001;10(2):266–277. [PubMed]

7. Dempster A, Laird N, Rubin D. Maximal likelihood form incomplete data via the EM algorithm. Proc of the Royal Statistical Society. 1977;39:1–38.

8. Hirayasu Y, Shenton ME, Salisbury DF, Dickey CC, Fischer IA, Mazzoni P, Kisler T, Arakaki H, Kwon JS, Anderson JE, Yurgelun-Todd D, Tohen M, McCarley RW. Lower left temporal lobe MRI volumes in patients with first-episode schizophrenia compared with psychotic patients with first-episode affective disorder and normal subjects. Amer J Psychiatry. 1998;155(10):1384–1391. [PubMed]

9. Fischl B, Salat D, Busa E, Albert M, Dietrich M, Haselgrove C, Van Der Kouwe A, Killany R, Kennedy D, Klaveness S, Montillo A, Makris N, Rosen B, Dale A. Whole brain segmentation: automated labeling of neuroanatomical structures in the human brain. Neuron. 2002;33(3):341–355. [PubMed]

10. Li S. Markov Random Field Modeling in Computer Vision. Springer; Heidelberg: 1995.

11. Lord N, Ho J, Vemuri B. ICCV. 2007. Ussr: A unified framework for simultaneous smoothing, segmentation, and registration of multiple images; pp. 1–6.

12. Paragios N, Deriche R. Geodesic active regions: A new paradigm to deal with frame partition problems in computer vision. JVCIR. 2002;13:249–268.

13. Pohl KM, Fisher J, Grimson WEL, Kikinis R, Wells WM. A bayesian model for joint segmentation and registration. NeuroImage. 2006;31(1):228–239. [PubMed]

14. Pohl KM, Bouix S, Nakamura M, Rohlfing T, McCarley RW, Kikinis R, Grimson WEL, Shenton ME, Wells WM. A hierarchical algorithm for MR brain image parcellation. TMI. 2007;26(9):1201–1212. [PMC free article] [PubMed]

15. Pohl KM, Fisher J, Bouix S, Shenton M, McCarley R, Grimson WEL, Kikinis R, Wells WM. Using the logarithm of odds to define a vector space on probabilistic atlases. Medical Image Analysis. 2007;11(6):465–477. [PMC free article] [PubMed]

16. Riklin-Raviv T, Sochen N, Kiryati N. Shape-based mutual segmentation. International Journal of Computer Vision. 2008;79:231–245.

17. Tu Z, Narr KL, Dollár P, Dinov I, Thompson PM, Toga AW. Brain anatomical structure segmentation by hybrid discriminative/generative models. IEEE TMI. 2008;27(4):495–508. [PMC free article] [PubMed]

18. Van Leemput K, Maes F, Vandermeulen D, Suetens P. Automated model-based tissue classification of MR images of the brain. IEEE TMI. 1999;18(10):897–908. [PubMed]

19. Yang J, Duncan J. CVPR. 2004. Joint prior models of neighboring object for 3D image segmentation; pp. 314–319. [PMC free article] [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |