Neighbor-Constrained Segmentation with 3D ... - Semantic Scholar
Neighbor-Constrained Segmentation with 3D Deformable Models Jing Yang1 , Lawrence H. Staib12 , and James S. Duncan12 Departments of Electrical Engineering1 and Diagnostic Radiology2 , Yale University, P.O. Box 208042, New Haven CT 06520-8042, USA, [email protected], [email protected], [email protected]
Abstract. A novel method for the segmentation of multiple objects from 3D medical images using inter-object constraints is presented. Our method is motivated by the observation that neighboring structures have consistent locations and shapes that provide configurations and context that aid in segmentation. We define a Maximum A Posteriori(MAP) estimation framework using the constraining information provided by neighboring objects to segment several objects simultaneously. We introduce a representation for the joint density function of the neighbor objects, and define joint probability distributions over the variations of the neighboring positions and shapes of a set of training images. By estimating the MAP shapes of the objects, we formulate the model in terms of level set functions, and compute the associated Euler-Lagrange equations. The contours evolve both according to the neighbor prior information and the image gray level information. We feel that this method is useful in situations where there is limited inter-object information as opposed to robust global atlases. Results and validation from various experiments on synthetic data and medical imagery in 2D and 3D are demonstrated.
1
Introduction
Segmentation and quantitative analysis of structures in an image has tremendous applications in medical imaging. In order to fully realize the value of medical imaging in both clinical and research settings, information about anatomical structures must be extracted and quantified with accuracy, efficiency, and reproducibility. Snakes or Active Contour Models(ACM)(Kass et al. (1987)) [1] have been widely used for segmenting non-rigid objects in a wide range of applications. ACMs are energy minimizing parametric contours with smoothness constraints deformed according to the image data. Unlike level set implementations[2], the direct implementation of this energy model is not capable of handling topological changes of the evolving contour without explicit discrete pixel manipulations. Usually, ACMs can detect only objects with edges defined by the gradient. Recently, methods using level set methods and new energy terms have been reported to increase the capture range of deformable models. For example, Chan and Vese [3] have proposed an active contour model that can detect objects whose boundaries are not necessarily defined by gray level gradients.
In addition to smoothness model, the incorporation of more specific prior information into deformable models has received a large amount of attention. Cootes et al. [4] find corresponding points across a set of training images and construct a statistical model of shape variation from the point positions. Staib and Duncan [5] incorporate global shape information into the segmentation process by using an elliptic Fourier decomposition of the boundary and placing a Gaussian prior on the Fourier coefficients. Zeng et al. [6] develop a coupled surfaces algorithm to segment the cortex by using a thickness prior constraint. Leventon et al. [7] extend Caselles’ [8] geodesic active contours by incorporating shape information into the evolution process. In many cases, objects to be detected have one or more neighboring structures whose location and shape provide information about the local geometry that can aid in the delineation. The relative positions or shapes among these neighbors can be modeled based on statistical information from a training set. Though applicable in many domains, these models are particularly useful for medical applications. In anatomy, neighboring structures provide a strong constraint on the relative position and shape of a structure. Without a prior model to constrain the segmentation, algorithms often fail due to the difficult challenges of poor image contrast, noise, and missing or diffuse boundaries. Segmentation can be made easier if suitable neighbor prior models are available. Our model is based on a MAP framework using the neighbor prior constraint. We introduce a representation for the joint density function of the neighbor objects and define the corresponding probability distributions. Formulating the segmentation as a MAP estimation of the shapes of the objects and modeling in terms of level set functions, we compute the associated Euler-Lagrange equations. The contours evolve both according to the neighbor prior information and the image gray level information. The neighboring objects can be automatically detected simultaneously.
2 2.1
Description of the Model MAP Framework with Neighbor Prior
A probabilistic formulation is a powerful approach to deformable models. Deformable models can be fit to the image data by finding the model shape parameters that maximize the posterior probability. Consider an image I that has M shapes of interest; a MAP framework can be used to realize image segmentation combining neighbor prior information and image information: Sˆi = arg maxSi p(S1 , S2 , ..., Si , ..., SM /I) ∝ arg maxSi p(I/S1 , S2 , ..., SM )p(S1 , S2 , ..., SM ), i = 1, 2, ..., M
(1)
where S1 , S2 , ..., SM are the evolving surfaces of all the shapes of interest. p(I/S1 , S2 , ..., SM ) is the probability of producing an image I given S1 , S2 , ... , SM . In 3D, assuming gray level homogeneity within each object, we use the
following imaging model[3]: QM
2 exp[−(I(p, q, r) − c1i )2 /(2σ1i )] 2 2 (p,q,r)outside(Si ),inside(Ωi ) exp[−(I(p, q, r) − c2i ) /(2σ2i )]} (2) where c1i and σ1i are the average and variance of I inside Si , c2i and σ2i are the average and variance of I outside Si but also inside a certain domain Ωi that contains Si . p(S1 , S2 , ..., SM ) is the joint density function of all the M objects. It contains the neighbor prior information such as the relative position and shape among the objects. By the chain rule, we have:
p(I/S1 , S2 , ..., SM ) = Q
i=1 {
Q
(p,q,r)inside(Si )
p(S1 , S2 , ..., SM ) = p(SM /SM −1 , SM −2 , ..., S1 )p(SM −1 /SM −2 , SM −3 , ..., S1 ) ...p(S3 /S2 , S1 )p(S2 /S1 )p(S1 ) (3) 2.2
Neighbor Priors
Let us define a binary matrix RM ×M where each element rij describes the independence of Si and Sj . rij has value zero when Si and Sj are independent and has value one otherwise. Obviously, the more ones there are in R, the more neighbor prior information is incorporated in the MAP segmentation model. When: 1 1 ... 1 1 1 ... 1 R= (4) ... 1 1 ... 1 all of the M objects are related to each other. The shape prior, as well as neighbor prior of all the neighbors, are included. In this case, equation (3) cannot be simplified and finding the joint density function of all the M objects is complicated. We incorporate the most neighbor prior information(full order prior) we can use but with the corresponding loss of computational efficiency. If all the M objects are independent to each other, i.e., 1 0 ... 0 0 1 ... 0 (5) R= ... 0 0 ... 1 then equation (3) can be simplified to: p(S1 , S2 , ..., SM ) = p(SM )p(SM −1 )...p(S2 )p(S1 )
(6)
No neighboring information is included here since all the objects are independent to each other. The only prior information contained in the MAP model for each object is the object’s self shape prior p(Si ), which we designate the first order
prior. In this case, equation (3) can be most simplified; we can achieve good computational efficiency but with no neighbor prior information. Each object in the image can be segmented independently according to its shape prior and image gray level information. This formulation corresponds to previous work[7][9]. In order to achieve a better segmentation result by using neighboring information, as well as good computational efficiency, we can consider second order prior information, i.e. the neighboring information from only one neighbor and the first order prior, i.e. the self shape information. Let us consider a simple case where each object is related to object 1 independently; the corresponding R is: 1 1 ... 1 1 1 ... 0 R= (7) ... 1 0 ... 1 The joint density function can be simplified to: p(S1 , S2 , ..., SM ) = p(SM /S1 )p(SM −1 /S1 )...p(S2 /S1 )p(S1 ) = p(∆M,1 )p(∆M −1,1 )...p(∆2,1 )p(S1 )
(8)
where ∆i,1 = Si − S1 is the difference between shape i and shape 1 (to be defined in the next section). The process of defining the joint density function p(S1 , S2 , ..., SM ) is simplified to building only the shape prior, p(S1 ), and the local neighbor prior p(∆i,1 ), i = 2, 3, ..., M . In our MAP model, we consider this case for the rest of the paper. 2.3
Neighbor Prior Model
To build a model for the neighbor prior and shape prior, we choose level sets as the representation of the shapes, and then define the joint probability density function in equation (8). Consider a training set of n aligned images, with M objects or structures in each image. Each shape in the training set is embedded as the zero level set of a higher dimensional level set Ψ . For object 1, the training set consists of a set of level set functions {Ψ1,1 , Ψ2,1 , ..., Ψn,1 }. We can use the difference between the two level sets, Ψi −Ψ1 , as the representation of the neighbor difference ∆i,1 , i = 2, 3, ..., M . Thus, the corresponding training set is {Ψ1,i − Ψ1,1 , Ψ2,i − Ψ2,1 , ..., Ψn,i − Ψn,1 }, i = 2, 3, ..., M . Our goal is to build the shape model and neighbor difference model over these distributions of the level set functions and level sets differences. The mean and variance of shape 1 can be computed Pn using Principal Component Analysis(PCA)[4]. The mean shape, Ψ¯1 = n1 l=1 Ψl,1 , is subtracted from each Ψl,1 to create the deviation from the mean. Each such deviation is placed as a column vector in a N d × n dimensional matrix Q where d is the number of spatial dimensions and N d is the number of samples of each level set function. Using Singular Value Decomposition(SVD), Q = U ΣV T . U is a matrix whose column vectors represent the set of orthogonal modes of shape variation and Σ
is a diagonal matrix of corresponding singular values. An estimate of the shape Ψ1 can be represented by k principal components and a k dimensional vector of coefficients(where k < n), α1 [7]: Ψ˜1 = Uk α1 + Ψ¯1
(9)
Under the assumption of a Gaussian distribution of shape represented by α1 , we can compute the probability of a certain shape: p(α1 ) = p
1 (2π)k |Σ
1 exp[− αT1 Σk−1 α1 ] 2 k|
(10)
Similarly, an estimate of the neighbor difference ∆i,1 can be represented ¯i,1 and k principal components Pik and a k from the mean neighbor difference ∆ dimensional vector of coefficients, βi,1 : ˜i,1 = Pik βi,1 + ∆ ¯i,1 ∆
(11)
The neighbor difference ∆i,1 can also be assumed to be Gaussian distributed over βi,1 : 1 T −1 1 Λi,1k βi,1 ] (12) exp[− βi,1 p(βi,1 ) = p k 2 (2π) |Λi,1k |
Fig. 1. Training set:outlines of 4 shapes in 12 3D MR brain images.
Figure 1 shows a training set of four sub-cortical structures:left and right amygdalas and hippocampuses, where we assume the left amygdala is related to all of the other three structures, independently. By using PCA, we can build the shape model of the left amygdala and the neighbor difference models of the other three structures. Figure 2 shows the three primary modes of variance of the left amygdala. Figure 3 shows the three primary modes of variance of the neighbor difference between the left hippocampus and the left amygdala.
Fig. 2. The three primary modes of variance of the left amygdala.
Fig. 3. The three primary modes of variance of the left hippocampus relative to the left amygdala.
In our active contour model, we also add some regularizing terms[10]: a genH QM −µi S ds i , and a eral boundary smoothness prior, pB (S1 , S2 , ..., SM ) = i=1 e QM −νi Ac i prior for the size of the region, pA (S1 , S2 , ..., SM ) = i=1 e , where Ai is the size of the region of shape i, c is a constant and µi and νi are scalar factors. Here we assume the boundary smoothness and the region size of all the objects are independent. Thus, the prior joint probability p(S1 , S2 , ..., SM ) in equation (8) can be approximated by a product of the following probabilities: M Y p(S1 , S2 , ..., SM ) = [ p(βi,1 )]·p(α1 )·pB (S1 , S2 , ..., SM )·pA (S1 , S2 , ..., SM ) (13) i=2
Since:
Sˆi = arg maxSi p(S1 , S2 , ..., Si , ..., SM /I) = arg minSi [− ln p(S1 , S2 , ..., Si , ..., SM /I)], i = 1, 2, ..., M
(14)
combining equation (1), (2), and (13), we introduce the energy functional E defined by E = − ln p(S1 , S2 , ..., Si , ..., SM /I) R P |I(p, q, r) − c1i |2 dpdqdr ∝ M · i=1 {λ R 1i (p,q,r)inside(Si ) +λ2i · (p,q,r)outside(Si ),inside(Ωi ) |I(p, q, r) − c2i |2 dpdqdr} H PM PM PM 1 T −1 T + i=1 µi Si ds + i=1 νi Aci + i=2 12 βi,1 Λ−1 i,1k βi,1 + 2 α1 Σk α1
(15)
The MAP estimation of the shapes in equation (1), Sˆi (i = 1, 2, ..., M ), is also the minimizer of the above energy functional E. This minimization problem can be formulated and solved using the level set method and we can realize the segmentation of multiple objects simultaneously. 2.4
Level Set Formulation of the Model
In the level set method, Si is the zero level set of a higher dimensional level set ψi corresponding to the i th object being segmented, i.e., Si = {(x, y, z)|ψi (x, y, z) =
0}. The evolution of surface Si is given by the zero-level surface at time t of the function ψi (t, x, y, z). We define ψi to be positive outside Si and negative inside Si . Each of the M objects being segmented in the image has its own Si and ψi . For the level set formulation of our model, we replace Si with ψi in the energy functional in equation (15) using regularized versions of the Heaviside function H and the Dirac function δ, denoted by Hε and δε [3](described below): Z δε (ψi (x, y, z))|∇ψi (x, y, z)|dxdydz E (c1i , c2i , ψi |i = 1, 2, ...M ) = µi Ω Z + νi (1 − Hε (ψi (x, y, z)))dxdydz ZΩ + λ1i |I(x, y, z) − c1i |2 (1 − Hε (ψi (x, y, z)))dxdydz Ω Z |I(x, y, z) − c2i |2 Hε (ψi (x, y, z))dxdydz + λ2i Ωi
+
M X 1 i=2
+
2
T ¯ [G(ψi − ψ1 ) − ∆¯i,1 ]T Pik Λ−1 i,1k Pik [G(ψi − ψ1 ) − ∆i,1 ]
1 [G(ψ1 − ψ¯1 )]T Uk Σk−1 UkT [G(ψ1 − ψ¯1 )] 2
(16)
where Ω denotes the image domain. G(·) is an operator to generate the vector representation(as in equation 9) of a matrix by column scanning. g(·) is the inverse operator of G(·). To compute the associated Euler-Lagrange equation for each unknown level set function ψi , we keep c1i and c2i fixed, and minimize E with respect to ψi (i = 1, 2, ...M ) respectively. Parameterizing the descent direction by artificial time t ≥ 0, the evolution equation in ψi (t, x, y, z) is: ∂ψi ∇ψi = δε (ψi )[µi · div[ ] + νi + λ1i |I − c1i |2 − λ2i |I − c2i |2 ] ∂t |∇ψi | T −H(i − 1.5) · g{Pik Λ−1 Pik [G(ψi − ψ1 ) − ∆¯i,1 )} i,1k
−[1 − H(i − 1.5)] · g{Uk Σk−1 UkT [G(ψ1 − ψ¯1 )]} 2.5
(17)
Evolving the Surface 2 z π arctan( ε )], δε (z) = I(x,y,z)·(1−H(ψi (x,y,z)))dxdydz
We approximate Hε and δε as follows [11]: Hε (z)R= 12 [1 + ε π(ε2 +z 2 ) .
c1i and c2i are defined by: c1i (ψi ) = R I(x,y,z)·H(ψi (x,y,z))dxdydz c2i (ψi ) = Ωi R . Ωi
Ω
R
Ω
(1−H(ψi (x,y,z)))dxdydz
,
H(ψi (x,y,z))dxdydz
Given the surfaces ψi (i = 1, 2, ...M ) at time t, we seek to compute the evolution steps that bring all the zero level set curves to the correct final segmentation based on the neighbor prior information and image information. We first set up p(α1 ) and p(βi,1 ), i = 2, 3, ..., M from the training set using PCA. At each stage
of the algorithm, we recompute the constants c1i (ψit ) and c2i (ψit ) and update ψit+1 . This is repeated until convergence. The parameters µi , νi , λ1i , and λ2i are used to balance the influence of the neighbor prior model and the image information model. The tradeoff between neighbor prior and image information depends on how much faith one has in the neighbor prior model and the imagery for a given application. We set these parameters empirically for particular segmentation tasks, given the general image quality and the neighbor prior information.
3
Experimental Results
We have used our model on various synthetic and real images, with at least two different types of contours and shapes. In Figure 4 top, we show the segmentation of the left and right ventricles using only image information, by which the curves cannot lock onto the shapes of the objects. In Figure 4 bottom, we show the results obtained using our model. The curves are able to converge on the desired boundaries even though some parts of the boundaries are too blurred to be detected using only gray level information. Both of the segmentations converged in several minutes on an SGI Octane with a 255MHz R10000 processor.
Fig. 4. Three steps in the segmentation of 2 shapes in a 2D cardiac MR image without(top) and with (bottom) neighbor prior. The right ventricle is the reference shape S1 . The training set consists of 16 images.
In Figure 5, we show that our model can detect multiple objects of different intensities and with blurred boundaries. Figure 5 top shows the results of using only gray level information. Only the lower (posterior) portions of the lateral ventricles can be segmented perfectly since they have clearer boundaries. Figure 5 bottom shows the results obtained using our neighbor prior model. Segmenting all eight subcortical structures took approximately twenty minutes.
Fig. 5. Detection of 8 sub-cortical structures(the lateral ventricles, heads of the caudate nucleus, and putamen) in a MR brain image. Top: results with no prior information. Bottom: results with neighbor prior. The left lateral ventricle is the reference shape S1 . The training set consists of 12 images.
Figure 6 shows the segmentation of the right amygdala and hippocampus in a 2D MR image. In Figure 6 top, we show results of using only gray level information. The segmentations are poor since both structures have very poorly defined boundaries. The middle row in Figure 6 shows the results of using the shape prior but with no neighbor prior. The results are much better, but the boundaries of the amygdala and the hippocampus are overlapped at the part where the two structures are connected. This is because the two structures are treated independently here without the constraint of the neighbor. In Figure 6 bottom, we show results of using our neighbor prior model, the two structures can be clearly segmented, and there is no overlap of the boundaries. We also tested our method in 3D images. We have generated a training set of 9 synthetic images of two uniform ellipsoids with added Gaussian noise. Figure 7 illustrates several steps in the segmentation of the two ellipsoids. Figure 8 shows initial, middle, and final steps in the segmentation of the left and right amygdalas and hippocampuses in a 3D MR brain image using training set model shown in Figures 1,2 and 3. Three orthogonal slices and the 3D surfaces are shown for each step. To validate the segmentation results, we compute the undirected distance between the boundary of the computed segmentation A(NA points) and the boundary of the P manual segmentation B: H(A, B) = max(h(A, B), h(B, A)), h(A, B) = N1A a∈A minb∈B ka − bk. Table 1 shows the computed results for the synthetic image, the heart image and the brain images. Virtually all the boundary points lie within one or two voxels of the manual segmentation.
Fig. 6. Four steps in the segmentation of right amygdala and hippocampus. Top: results with no prior information. Middle: results using individual shape priors. Bottom: results using our neighbor prior model. The right amygdala is the reference shape S1 . The training set consists of 12 brain images.
Fig. 7. Initial, middle, and final steps in the segmentation of 2 shapes in a synthetic image. Three orthogonal slices and the 3D surfaces are shown for each step.
Table 1. Distance between the computed boundary and the manual boundary Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Without neighbor prior With neighbor prior
4.2 1.8
9.6 1.9
6.7 0.8
11.2 1.2
7.8 1.7
Fig. 8. Initial, middle, and final steps in the segmentation of 4 shapes in a brain image. Three orthogonal slices and the 3D surfaces are shown for each step.
We also test the robustness of our algorithm to noise as well as to the location of the initial seeds. First, we add Gaussian noise with different variances to the synthetic image(as in Figure 7), and run our algorithm to segment the two ellipsoids, where we set the initial seeds at the centers of the objects. Figure 9 shows the segmentation error in three cases: with no prior, with shape prior, and with neighbor prior. When the variance of the noise is small, the errors are also small for all the three cases. As the variance of the noise goes up, the error for no prior increases rapidly since the objects are too noisy to be recognized using only gray level information. However, for the methods with the shape prior and with the neighbor prior, the segmentation errors are much lower and are locked in a very small range even when the variance of the noise is very large. We also notice that our neighbor prior model achieves the least error among all the cases. Next, we fix the standard deviation of the noise to be 40, but vary the location of the initial seed inside the right ellipsoid and run the segmentation for the same three cases again. The segmentation error for different seed locations with each method is shown in Figure 10. As the initial seed goes far away from the center of the ellipsoid, the errors are kept in a small range for all the cases because the models are based on level sets. Still, the method with the neighbor prior achieves the smallest error.
4
Conclusions
A new model for automated segmentation of images containing multiple objects by incorporating neighbor prior information in the segmentation process has been presented. We wanted to capture the constraining information that neighboring objects provided and use it for segmentation. We define a MAP estimation framework using the prior information provided by neighboring objects to segment several objects simultaneously. We introduce a representation for the joint density function of the neighbor objects, and define joint probability distributions over the variations of the neighboring positions and shapes in a
12
4.5
Without Prior With Shape Prior With Neighbor Prior
4
10 3.5
8 Without Prior With Shape Prior With Neighbor Prior
Error
Error
3
6
2.5
2
4 1.5
2 1
0
0
10
20
30 40 Standard Deviation of the Noise
50
60
70
Fig. 9. Segmentation errors with different variances of the noise.
0.5
0
5
10 15 20 25 Distance between the initial seed and the center of the object
30
Fig. 10. Segmentation errors with different locations of initial seed.
set of training images. We estimate the MAP shapes of the objects using evolving level sets based on the associated Euler-Lagrange equations. The contours evolve both according to the neighbor prior information and the image gray level information. Multiple objects in an image can be automatically detected simultaneously.
References 1. M. Kass, A. Witkin, D. Terzopoulos.: Snakes: Active contour models. Int’l Journal on Computer Vision, 1 (1987) 321–331 2. S. Osher and J. A. Sethian.: Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi Formulation. J. Comp. Phy. 79 (1988) 12–49 3. T. Chan, L. Vese.: Active Contours Without Edges. IEEE Transactions on Image Processing, vol.10 No. 2 (2001) 266–277 4. T.F. Cootes, A. Hill, C.J. Taylor, and J. Haslam.: Use of active shape models for locating structures in medical images. Image and Vision Computing, 12(6):355-365, July 1994. 5. L.Staib, J. Duncan.: Boundary finding with parametrically deformable models. PAMI, 14(11) (1992) 1061–1075 6. X. Zeng, L.H.Staib, R.T.Schultz and J.S.Duncan.: Volumetric Layer Segmentation Using Coupled Surfaces Propagation. IEEE Conf. on Comp. Vision and Patt. Recog. (1998). 7. M. Leventon, E. Grimson, and O. Faugeras.: Statistical shape influence in geodesic active contours. IEEE Conf. on Comp. Vision and Patt. Recog. 1 (2000) 316–323 8. V. Caselles, R. Kimmel, and G. Sapiro.: Geodesic active contours. Int. J. Comput. Vis. vol.22 No. 1 (1997) 61–79 9. A. Tsai, A. Yezzi, et al.: Model-based curve evolution technique for image segmentation. IEEE Conf. on Comp. Vision and Patt. Recog. vol.1 (2001) 463–468 10. Z. Tu, and S. Zhu: Image segmentation by data-driven Markov Chain Monte Carlo. IEEE Trans. on Patt. Ana. and Machine Intelligence. vol.24 No. 5 (2002) 657–673 11. J. Yang, L. Staib and J. Duncan: Statistical Neighbor Distance Influence in Active Contours. MICCAI. vol.1 (2002) 588–596
Abstract. A novel method for the segmentation of multiple objects from 3D medical images using inter-object constraints is presented. Our method is motivated by the observation that neighboring structures have consistent locations and shapes that provide configurations and context that aid in segmentation. We define a Maximum A Posteriori(MAP) estimation framework using the constraining information provided by neighboring objects to segment several objects simultaneously. We introduce a representation for the joint density function of the neighbor objects, and define joint probability distributions over the variations of the neighboring positions and shapes of a set of training images. By estimating the MAP shapes of the objects, we formulate the model in terms of level set functions, and compute the associated Euler-Lagrange equations. The contours evolve both according to the neighbor prior information and the image gray level information. We feel that this method is useful in situations where there is limited inter-object information as opposed to robust global atlases. Results and validation from various experiments on synthetic data and medical imagery in 2D and 3D are demonstrated.
1
Introduction
Segmentation and quantitative analysis of structures in an image has tremendous applications in medical imaging. In order to fully realize the value of medical imaging in both clinical and research settings, information about anatomical structures must be extracted and quantified with accuracy, efficiency, and reproducibility. Snakes or Active Contour Models(ACM)(Kass et al. (1987)) [1] have been widely used for segmenting non-rigid objects in a wide range of applications. ACMs are energy minimizing parametric contours with smoothness constraints deformed according to the image data. Unlike level set implementations[2], the direct implementation of this energy model is not capable of handling topological changes of the evolving contour without explicit discrete pixel manipulations. Usually, ACMs can detect only objects with edges defined by the gradient. Recently, methods using level set methods and new energy terms have been reported to increase the capture range of deformable models. For example, Chan and Vese [3] have proposed an active contour model that can detect objects whose boundaries are not necessarily defined by gray level gradients.
In addition to smoothness model, the incorporation of more specific prior information into deformable models has received a large amount of attention. Cootes et al. [4] find corresponding points across a set of training images and construct a statistical model of shape variation from the point positions. Staib and Duncan [5] incorporate global shape information into the segmentation process by using an elliptic Fourier decomposition of the boundary and placing a Gaussian prior on the Fourier coefficients. Zeng et al. [6] develop a coupled surfaces algorithm to segment the cortex by using a thickness prior constraint. Leventon et al. [7] extend Caselles’ [8] geodesic active contours by incorporating shape information into the evolution process. In many cases, objects to be detected have one or more neighboring structures whose location and shape provide information about the local geometry that can aid in the delineation. The relative positions or shapes among these neighbors can be modeled based on statistical information from a training set. Though applicable in many domains, these models are particularly useful for medical applications. In anatomy, neighboring structures provide a strong constraint on the relative position and shape of a structure. Without a prior model to constrain the segmentation, algorithms often fail due to the difficult challenges of poor image contrast, noise, and missing or diffuse boundaries. Segmentation can be made easier if suitable neighbor prior models are available. Our model is based on a MAP framework using the neighbor prior constraint. We introduce a representation for the joint density function of the neighbor objects and define the corresponding probability distributions. Formulating the segmentation as a MAP estimation of the shapes of the objects and modeling in terms of level set functions, we compute the associated Euler-Lagrange equations. The contours evolve both according to the neighbor prior information and the image gray level information. The neighboring objects can be automatically detected simultaneously.
2 2.1
Description of the Model MAP Framework with Neighbor Prior
A probabilistic formulation is a powerful approach to deformable models. Deformable models can be fit to the image data by finding the model shape parameters that maximize the posterior probability. Consider an image I that has M shapes of interest; a MAP framework can be used to realize image segmentation combining neighbor prior information and image information: Sˆi = arg maxSi p(S1 , S2 , ..., Si , ..., SM /I) ∝ arg maxSi p(I/S1 , S2 , ..., SM )p(S1 , S2 , ..., SM ), i = 1, 2, ..., M
(1)
where S1 , S2 , ..., SM are the evolving surfaces of all the shapes of interest. p(I/S1 , S2 , ..., SM ) is the probability of producing an image I given S1 , S2 , ... , SM . In 3D, assuming gray level homogeneity within each object, we use the
following imaging model[3]: QM
2 exp[−(I(p, q, r) − c1i )2 /(2σ1i )] 2 2 (p,q,r)outside(Si ),inside(Ωi ) exp[−(I(p, q, r) − c2i ) /(2σ2i )]} (2) where c1i and σ1i are the average and variance of I inside Si , c2i and σ2i are the average and variance of I outside Si but also inside a certain domain Ωi that contains Si . p(S1 , S2 , ..., SM ) is the joint density function of all the M objects. It contains the neighbor prior information such as the relative position and shape among the objects. By the chain rule, we have:
p(I/S1 , S2 , ..., SM ) = Q
i=1 {
Q
(p,q,r)inside(Si )
p(S1 , S2 , ..., SM ) = p(SM /SM −1 , SM −2 , ..., S1 )p(SM −1 /SM −2 , SM −3 , ..., S1 ) ...p(S3 /S2 , S1 )p(S2 /S1 )p(S1 ) (3) 2.2
Neighbor Priors
Let us define a binary matrix RM ×M where each element rij describes the independence of Si and Sj . rij has value zero when Si and Sj are independent and has value one otherwise. Obviously, the more ones there are in R, the more neighbor prior information is incorporated in the MAP segmentation model. When: 1 1 ... 1 1 1 ... 1 R= (4) ... 1 1 ... 1 all of the M objects are related to each other. The shape prior, as well as neighbor prior of all the neighbors, are included. In this case, equation (3) cannot be simplified and finding the joint density function of all the M objects is complicated. We incorporate the most neighbor prior information(full order prior) we can use but with the corresponding loss of computational efficiency. If all the M objects are independent to each other, i.e., 1 0 ... 0 0 1 ... 0 (5) R= ... 0 0 ... 1 then equation (3) can be simplified to: p(S1 , S2 , ..., SM ) = p(SM )p(SM −1 )...p(S2 )p(S1 )
(6)
No neighboring information is included here since all the objects are independent to each other. The only prior information contained in the MAP model for each object is the object’s self shape prior p(Si ), which we designate the first order
prior. In this case, equation (3) can be most simplified; we can achieve good computational efficiency but with no neighbor prior information. Each object in the image can be segmented independently according to its shape prior and image gray level information. This formulation corresponds to previous work[7][9]. In order to achieve a better segmentation result by using neighboring information, as well as good computational efficiency, we can consider second order prior information, i.e. the neighboring information from only one neighbor and the first order prior, i.e. the self shape information. Let us consider a simple case where each object is related to object 1 independently; the corresponding R is: 1 1 ... 1 1 1 ... 0 R= (7) ... 1 0 ... 1 The joint density function can be simplified to: p(S1 , S2 , ..., SM ) = p(SM /S1 )p(SM −1 /S1 )...p(S2 /S1 )p(S1 ) = p(∆M,1 )p(∆M −1,1 )...p(∆2,1 )p(S1 )
(8)
where ∆i,1 = Si − S1 is the difference between shape i and shape 1 (to be defined in the next section). The process of defining the joint density function p(S1 , S2 , ..., SM ) is simplified to building only the shape prior, p(S1 ), and the local neighbor prior p(∆i,1 ), i = 2, 3, ..., M . In our MAP model, we consider this case for the rest of the paper. 2.3
Neighbor Prior Model
To build a model for the neighbor prior and shape prior, we choose level sets as the representation of the shapes, and then define the joint probability density function in equation (8). Consider a training set of n aligned images, with M objects or structures in each image. Each shape in the training set is embedded as the zero level set of a higher dimensional level set Ψ . For object 1, the training set consists of a set of level set functions {Ψ1,1 , Ψ2,1 , ..., Ψn,1 }. We can use the difference between the two level sets, Ψi −Ψ1 , as the representation of the neighbor difference ∆i,1 , i = 2, 3, ..., M . Thus, the corresponding training set is {Ψ1,i − Ψ1,1 , Ψ2,i − Ψ2,1 , ..., Ψn,i − Ψn,1 }, i = 2, 3, ..., M . Our goal is to build the shape model and neighbor difference model over these distributions of the level set functions and level sets differences. The mean and variance of shape 1 can be computed Pn using Principal Component Analysis(PCA)[4]. The mean shape, Ψ¯1 = n1 l=1 Ψl,1 , is subtracted from each Ψl,1 to create the deviation from the mean. Each such deviation is placed as a column vector in a N d × n dimensional matrix Q where d is the number of spatial dimensions and N d is the number of samples of each level set function. Using Singular Value Decomposition(SVD), Q = U ΣV T . U is a matrix whose column vectors represent the set of orthogonal modes of shape variation and Σ
is a diagonal matrix of corresponding singular values. An estimate of the shape Ψ1 can be represented by k principal components and a k dimensional vector of coefficients(where k < n), α1 [7]: Ψ˜1 = Uk α1 + Ψ¯1
(9)
Under the assumption of a Gaussian distribution of shape represented by α1 , we can compute the probability of a certain shape: p(α1 ) = p
1 (2π)k |Σ
1 exp[− αT1 Σk−1 α1 ] 2 k|
(10)
Similarly, an estimate of the neighbor difference ∆i,1 can be represented ¯i,1 and k principal components Pik and a k from the mean neighbor difference ∆ dimensional vector of coefficients, βi,1 : ˜i,1 = Pik βi,1 + ∆ ¯i,1 ∆
(11)
The neighbor difference ∆i,1 can also be assumed to be Gaussian distributed over βi,1 : 1 T −1 1 Λi,1k βi,1 ] (12) exp[− βi,1 p(βi,1 ) = p k 2 (2π) |Λi,1k |
Fig. 1. Training set:outlines of 4 shapes in 12 3D MR brain images.
Figure 1 shows a training set of four sub-cortical structures:left and right amygdalas and hippocampuses, where we assume the left amygdala is related to all of the other three structures, independently. By using PCA, we can build the shape model of the left amygdala and the neighbor difference models of the other three structures. Figure 2 shows the three primary modes of variance of the left amygdala. Figure 3 shows the three primary modes of variance of the neighbor difference between the left hippocampus and the left amygdala.
Fig. 2. The three primary modes of variance of the left amygdala.
Fig. 3. The three primary modes of variance of the left hippocampus relative to the left amygdala.
In our active contour model, we also add some regularizing terms[10]: a genH QM −µi S ds i , and a eral boundary smoothness prior, pB (S1 , S2 , ..., SM ) = i=1 e QM −νi Ac i prior for the size of the region, pA (S1 , S2 , ..., SM ) = i=1 e , where Ai is the size of the region of shape i, c is a constant and µi and νi are scalar factors. Here we assume the boundary smoothness and the region size of all the objects are independent. Thus, the prior joint probability p(S1 , S2 , ..., SM ) in equation (8) can be approximated by a product of the following probabilities: M Y p(S1 , S2 , ..., SM ) = [ p(βi,1 )]·p(α1 )·pB (S1 , S2 , ..., SM )·pA (S1 , S2 , ..., SM ) (13) i=2
Since:
Sˆi = arg maxSi p(S1 , S2 , ..., Si , ..., SM /I) = arg minSi [− ln p(S1 , S2 , ..., Si , ..., SM /I)], i = 1, 2, ..., M
(14)
combining equation (1), (2), and (13), we introduce the energy functional E defined by E = − ln p(S1 , S2 , ..., Si , ..., SM /I) R P |I(p, q, r) − c1i |2 dpdqdr ∝ M · i=1 {λ R 1i (p,q,r)inside(Si ) +λ2i · (p,q,r)outside(Si ),inside(Ωi ) |I(p, q, r) − c2i |2 dpdqdr} H PM PM PM 1 T −1 T + i=1 µi Si ds + i=1 νi Aci + i=2 12 βi,1 Λ−1 i,1k βi,1 + 2 α1 Σk α1
(15)
The MAP estimation of the shapes in equation (1), Sˆi (i = 1, 2, ..., M ), is also the minimizer of the above energy functional E. This minimization problem can be formulated and solved using the level set method and we can realize the segmentation of multiple objects simultaneously. 2.4
Level Set Formulation of the Model
In the level set method, Si is the zero level set of a higher dimensional level set ψi corresponding to the i th object being segmented, i.e., Si = {(x, y, z)|ψi (x, y, z) =
0}. The evolution of surface Si is given by the zero-level surface at time t of the function ψi (t, x, y, z). We define ψi to be positive outside Si and negative inside Si . Each of the M objects being segmented in the image has its own Si and ψi . For the level set formulation of our model, we replace Si with ψi in the energy functional in equation (15) using regularized versions of the Heaviside function H and the Dirac function δ, denoted by Hε and δε [3](described below): Z δε (ψi (x, y, z))|∇ψi (x, y, z)|dxdydz E (c1i , c2i , ψi |i = 1, 2, ...M ) = µi Ω Z + νi (1 − Hε (ψi (x, y, z)))dxdydz ZΩ + λ1i |I(x, y, z) − c1i |2 (1 − Hε (ψi (x, y, z)))dxdydz Ω Z |I(x, y, z) − c2i |2 Hε (ψi (x, y, z))dxdydz + λ2i Ωi
+
M X 1 i=2
+
2
T ¯ [G(ψi − ψ1 ) − ∆¯i,1 ]T Pik Λ−1 i,1k Pik [G(ψi − ψ1 ) − ∆i,1 ]
1 [G(ψ1 − ψ¯1 )]T Uk Σk−1 UkT [G(ψ1 − ψ¯1 )] 2
(16)
where Ω denotes the image domain. G(·) is an operator to generate the vector representation(as in equation 9) of a matrix by column scanning. g(·) is the inverse operator of G(·). To compute the associated Euler-Lagrange equation for each unknown level set function ψi , we keep c1i and c2i fixed, and minimize E with respect to ψi (i = 1, 2, ...M ) respectively. Parameterizing the descent direction by artificial time t ≥ 0, the evolution equation in ψi (t, x, y, z) is: ∂ψi ∇ψi = δε (ψi )[µi · div[ ] + νi + λ1i |I − c1i |2 − λ2i |I − c2i |2 ] ∂t |∇ψi | T −H(i − 1.5) · g{Pik Λ−1 Pik [G(ψi − ψ1 ) − ∆¯i,1 )} i,1k
−[1 − H(i − 1.5)] · g{Uk Σk−1 UkT [G(ψ1 − ψ¯1 )]} 2.5
(17)
Evolving the Surface 2 z π arctan( ε )], δε (z) = I(x,y,z)·(1−H(ψi (x,y,z)))dxdydz
We approximate Hε and δε as follows [11]: Hε (z)R= 12 [1 + ε π(ε2 +z 2 ) .
c1i and c2i are defined by: c1i (ψi ) = R I(x,y,z)·H(ψi (x,y,z))dxdydz c2i (ψi ) = Ωi R . Ωi
Ω
R
Ω
(1−H(ψi (x,y,z)))dxdydz
,
H(ψi (x,y,z))dxdydz
Given the surfaces ψi (i = 1, 2, ...M ) at time t, we seek to compute the evolution steps that bring all the zero level set curves to the correct final segmentation based on the neighbor prior information and image information. We first set up p(α1 ) and p(βi,1 ), i = 2, 3, ..., M from the training set using PCA. At each stage
of the algorithm, we recompute the constants c1i (ψit ) and c2i (ψit ) and update ψit+1 . This is repeated until convergence. The parameters µi , νi , λ1i , and λ2i are used to balance the influence of the neighbor prior model and the image information model. The tradeoff between neighbor prior and image information depends on how much faith one has in the neighbor prior model and the imagery for a given application. We set these parameters empirically for particular segmentation tasks, given the general image quality and the neighbor prior information.
3
Experimental Results
We have used our model on various synthetic and real images, with at least two different types of contours and shapes. In Figure 4 top, we show the segmentation of the left and right ventricles using only image information, by which the curves cannot lock onto the shapes of the objects. In Figure 4 bottom, we show the results obtained using our model. The curves are able to converge on the desired boundaries even though some parts of the boundaries are too blurred to be detected using only gray level information. Both of the segmentations converged in several minutes on an SGI Octane with a 255MHz R10000 processor.
Fig. 4. Three steps in the segmentation of 2 shapes in a 2D cardiac MR image without(top) and with (bottom) neighbor prior. The right ventricle is the reference shape S1 . The training set consists of 16 images.
In Figure 5, we show that our model can detect multiple objects of different intensities and with blurred boundaries. Figure 5 top shows the results of using only gray level information. Only the lower (posterior) portions of the lateral ventricles can be segmented perfectly since they have clearer boundaries. Figure 5 bottom shows the results obtained using our neighbor prior model. Segmenting all eight subcortical structures took approximately twenty minutes.
Fig. 5. Detection of 8 sub-cortical structures(the lateral ventricles, heads of the caudate nucleus, and putamen) in a MR brain image. Top: results with no prior information. Bottom: results with neighbor prior. The left lateral ventricle is the reference shape S1 . The training set consists of 12 images.
Figure 6 shows the segmentation of the right amygdala and hippocampus in a 2D MR image. In Figure 6 top, we show results of using only gray level information. The segmentations are poor since both structures have very poorly defined boundaries. The middle row in Figure 6 shows the results of using the shape prior but with no neighbor prior. The results are much better, but the boundaries of the amygdala and the hippocampus are overlapped at the part where the two structures are connected. This is because the two structures are treated independently here without the constraint of the neighbor. In Figure 6 bottom, we show results of using our neighbor prior model, the two structures can be clearly segmented, and there is no overlap of the boundaries. We also tested our method in 3D images. We have generated a training set of 9 synthetic images of two uniform ellipsoids with added Gaussian noise. Figure 7 illustrates several steps in the segmentation of the two ellipsoids. Figure 8 shows initial, middle, and final steps in the segmentation of the left and right amygdalas and hippocampuses in a 3D MR brain image using training set model shown in Figures 1,2 and 3. Three orthogonal slices and the 3D surfaces are shown for each step. To validate the segmentation results, we compute the undirected distance between the boundary of the computed segmentation A(NA points) and the boundary of the P manual segmentation B: H(A, B) = max(h(A, B), h(B, A)), h(A, B) = N1A a∈A minb∈B ka − bk. Table 1 shows the computed results for the synthetic image, the heart image and the brain images. Virtually all the boundary points lie within one or two voxels of the manual segmentation.
Fig. 6. Four steps in the segmentation of right amygdala and hippocampus. Top: results with no prior information. Middle: results using individual shape priors. Bottom: results using our neighbor prior model. The right amygdala is the reference shape S1 . The training set consists of 12 brain images.
Fig. 7. Initial, middle, and final steps in the segmentation of 2 shapes in a synthetic image. Three orthogonal slices and the 3D surfaces are shown for each step.
Table 1. Distance between the computed boundary and the manual boundary Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Without neighbor prior With neighbor prior
4.2 1.8
9.6 1.9
6.7 0.8
11.2 1.2
7.8 1.7
Fig. 8. Initial, middle, and final steps in the segmentation of 4 shapes in a brain image. Three orthogonal slices and the 3D surfaces are shown for each step.
We also test the robustness of our algorithm to noise as well as to the location of the initial seeds. First, we add Gaussian noise with different variances to the synthetic image(as in Figure 7), and run our algorithm to segment the two ellipsoids, where we set the initial seeds at the centers of the objects. Figure 9 shows the segmentation error in three cases: with no prior, with shape prior, and with neighbor prior. When the variance of the noise is small, the errors are also small for all the three cases. As the variance of the noise goes up, the error for no prior increases rapidly since the objects are too noisy to be recognized using only gray level information. However, for the methods with the shape prior and with the neighbor prior, the segmentation errors are much lower and are locked in a very small range even when the variance of the noise is very large. We also notice that our neighbor prior model achieves the least error among all the cases. Next, we fix the standard deviation of the noise to be 40, but vary the location of the initial seed inside the right ellipsoid and run the segmentation for the same three cases again. The segmentation error for different seed locations with each method is shown in Figure 10. As the initial seed goes far away from the center of the ellipsoid, the errors are kept in a small range for all the cases because the models are based on level sets. Still, the method with the neighbor prior achieves the smallest error.
4
Conclusions
A new model for automated segmentation of images containing multiple objects by incorporating neighbor prior information in the segmentation process has been presented. We wanted to capture the constraining information that neighboring objects provided and use it for segmentation. We define a MAP estimation framework using the prior information provided by neighboring objects to segment several objects simultaneously. We introduce a representation for the joint density function of the neighbor objects, and define joint probability distributions over the variations of the neighboring positions and shapes in a
12
4.5
Without Prior With Shape Prior With Neighbor Prior
4
10 3.5
8 Without Prior With Shape Prior With Neighbor Prior
Error
Error
3
6
2.5
2
4 1.5
2 1
0
0
10
20
30 40 Standard Deviation of the Noise
50
60
70
Fig. 9. Segmentation errors with different variances of the noise.
0.5
0
5
10 15 20 25 Distance between the initial seed and the center of the object
30
Fig. 10. Segmentation errors with different locations of initial seed.
set of training images. We estimate the MAP shapes of the objects using evolving level sets based on the associated Euler-Lagrange equations. The contours evolve both according to the neighbor prior information and the image gray level information. Multiple objects in an image can be automatically detected simultaneously.
References 1. M. Kass, A. Witkin, D. Terzopoulos.: Snakes: Active contour models. Int’l Journal on Computer Vision, 1 (1987) 321–331 2. S. Osher and J. A. Sethian.: Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi Formulation. J. Comp. Phy. 79 (1988) 12–49 3. T. Chan, L. Vese.: Active Contours Without Edges. IEEE Transactions on Image Processing, vol.10 No. 2 (2001) 266–277 4. T.F. Cootes, A. Hill, C.J. Taylor, and J. Haslam.: Use of active shape models for locating structures in medical images. Image and Vision Computing, 12(6):355-365, July 1994. 5. L.Staib, J. Duncan.: Boundary finding with parametrically deformable models. PAMI, 14(11) (1992) 1061–1075 6. X. Zeng, L.H.Staib, R.T.Schultz and J.S.Duncan.: Volumetric Layer Segmentation Using Coupled Surfaces Propagation. IEEE Conf. on Comp. Vision and Patt. Recog. (1998). 7. M. Leventon, E. Grimson, and O. Faugeras.: Statistical shape influence in geodesic active contours. IEEE Conf. on Comp. Vision and Patt. Recog. 1 (2000) 316–323 8. V. Caselles, R. Kimmel, and G. Sapiro.: Geodesic active contours. Int. J. Comput. Vis. vol.22 No. 1 (1997) 61–79 9. A. Tsai, A. Yezzi, et al.: Model-based curve evolution technique for image segmentation. IEEE Conf. on Comp. Vision and Patt. Recog. vol.1 (2001) 463–468 10. Z. Tu, and S. Zhu: Image segmentation by data-driven Markov Chain Monte Carlo. IEEE Trans. on Patt. Ana. and Machine Intelligence. vol.24 No. 5 (2002) 657–673 11. J. Yang, L. Staib and J. Duncan: Statistical Neighbor Distance Influence in Active Contours. MICCAI. vol.1 (2002) 588–596
