|Home | About | Journals | Submit | Contact Us | Français|
Segmentation of anatomic structures from magnetic resonance brain scans can be a daunting task because of large inhomogeneities in image intensities across an image and possible lack of precisely defined shape boundaries for certain anatomical structures. One approach that has been quite popular in the recent past for these situations is the atlas-based segmentation. The atlas, once constructed, can be used as a template and can be registered nonrigidly to the image being segmented thereby achieving the desired segmentation. The goal of our study is to segment these structures with a registration assisted image segmentation technique.
We present a novel variational formulation of the registration assisted image segmentation problem which leads to solving a coupled set of nonlinear Partial Differential Equations (PDEs) that are solved using efficient numeric schemes. Our work is a departure from earlier methods in that we can simultaneously register and segment in three dimensions and easily cope with situations where the source (atlas) and target images have very distinct intensity distributions.
The proposed atlas-based segmentation technique is capable of simultaneously achieve the nonrigid registration and the segmentation; unlike previous methods of solution for this problem, our algorithm can accommodate for image pairs having very distinct intensity distributions.
In medical imaging applications, segmentation can be a daunting task because of possibly large inhomogeneities in image intensities across an image (eg, in magnetic resonance imaging [MRI]). These inhomogeneities combined with volume averaging during the imaging and possible lack of precisely defined shape boundaries for certain anatomical structures complicates the segmentation problem immensely. One possible solution for such situations is atlas-based segmentation. The atlas, once constructed, can be used as a template and can be registered nonrigidly to the image being segmented (henceforth called a target image), thereby achieving the desired segmentation. This approach has been quite popular in the recent past for segmenting cortical and subcortical structures in the human brain from MRI images. Many of the methods that achieve atlas-based segmentation are based on a two-stage process involving estimating the nonrigid deformation field between the atlas image and the target image and then applying the estimated deformation field to the desired shape/atlas to achieve the segmentation of the corresponding structures in the target image. In the following, we will briefly review some of these methods.
In Chen et al. (1), an atlas-based segmentation scheme is described that uses expert input to define the atlas and then warp it to the subject brain MRI to segment the subject brain, which is then followed by morphometrics. Their algorithm is almost identical to the work reported in Dawant et al. (2) and Thirion (3). All these methods are not strictly derived from a variational principle and hence can not be interpreted as minimizers of any cost functionals. Doing so will accrue the benefits of mathematical analysis of the solution (eg, it is possible to discuss the existence and uniqueness of the solution, which is otherwise not possible). Moreover, they do not achieve registration and segmentation simultaneously in a unified mathematical framework. Achieving this in a unified framework leads us to solve and analyze the errors in a single framework, which is much simpler and more elegant.
Among the image registration methods that can be put to use in achieving model-based segmentation, one popular approach is based on maximizing mutual information (MI) reported simultaneously in Viola and Wells (4) and Collignon et al. (5). Impressive experiments have been reported using this scheme for rigid and nonrigid registration (6-8). These methods have been used for constructing atlases (9). Wang et al. presented a novel measure of information called cumulative residual entropy (CRE) and developed a cost function called the cross-CRE to measure the similarity between two images (10). They presented results of registration under a variety of noise levels and showed significantly superior performance over MI-based methods. None of the methods that exist in literature use shape priors in an information theoretic cost functional for achieving atlas-based segmentation of the hippocampus or other neuroanatomical structures.
Atlas-based segmentation (11-13) was achieved using the so called fluid-flow model introduced by Christensen et al. (14). A more general and mathematically thorough treatment of the nonrigid registration, which subsumes the fluid-flow methods was presented in Trouve (15). The method in Trouve (15) has not been used in atlas-based segmentation applications, however. Once again, these methods do not use a unified framework for registration + segmentation and hence may be limited in their accuracy by the inaccuracies in each stage.
On the topic of hippocampal shape recovery from brain MRI (16,17), a level-set based image registration algorithm was introduced that was designed to nonrigidly register two three-dimensional (3D) volumes from the same modality of imaging. This algorithm was computationally efficient and was applied to achieve atlas-based segmentation. Although this scheme was mathematically well grounded, the task of atlas-based segmentation was achieved in two separate stages as in the other methods mentioned previously. Moreover, the technique was restricted in its ability to deal with image pairs (atlas image and target image) that were almost similar in intensity distributions. Yezzi et al. (18) and Noble et al. (19) presented a joint registration and segmentation scheme wherein the registration could only accommodate rigid motion. Also, in Fischl et al. (20), a Bayesian method is presented that simultaneously estimates a linear registration and the segmentation of a novel image. Note that linear registration does not involve nonrigid deformations and hence their method does not handle nonrigid registration in the joint registration + segmentation. The case of joint registration and segmentation with nonrigid registration has not been addressed adequately in literature with the exception of the recent work reported in Soatto et al. (21,22) and Vemuri et al. (23). However, all of these methods can only work with image pairs that are necessarily from the same modality or the intensity profiles are not too disparate.
In this article, we present a unified variational principle that will simultaneously register the atlas shape to the novel brain image and segment the desired shape in the novel image. The salient features of our algorithm are it accommodates for image pairs having distinct intensity distributions. To the best of our knowledge, it is the first algorithm for joint multimodal image registration and segmentation, the formulation of the joint (nonrigid) registration and segmentation problem is in a unified variational framework. In this work, the atlas serves in the segmentation process as a prior and the registration of this before the novel brain scan will assist in segmenting it.
We now present our formulation of joint registration & segmentation model. Let I1 be the atlas image containing the atlas shape C, I2 the novel image that needs to be segmented and v be the vector field, from I2 to I1 (i.e., the transformation is centered in I2), defining the nonrigid deformation between the two images. The variational principle describing our formulation of the registration assisted segmentation problem is given by:
where the first term denotes the segmentation functional. is the boundary contour (surface in 3D) of the desired anatomic shape in I2. The second term measures the distance between the transformed atlas v(C) and the current segmentation in the novel brain image (i.e., the target image and the third term denotes the nonrigid registration functional between the two images). Our joint registration and segmentation model is illustrated in Figure 1.
For the segmentation functional, we use a piecewise constant Mumford-Shah model, which is one of the well-known variational models for image segmentation, wherein it is assumed that the image to be segmented can be modeled by piecewise constant regions, as was done previously (24). This assumption simplifies our presentation but our model itself can be easily extended to the piecewise smooth regions case. Additionally, because we are only interested in segmenting a desired anatomical shape (eg, the hippocampus, the corpus callosum), we will only be concerned with a binary segmentation (i.e., two classes, namely voxels inside the desired shape and those that are outside it). These assumptions can be easily relaxed if necessary but at the cost of making the energy functional more complicated and hence computationally more challenging. The segmentation functional takes the following form:
where Ω is the image domain and α is a regularization parameter. u = ui if × in and u = uo if × out. in and out denote the regions inside and outside of the curve, representing the desired shape boundaries in I2.
For the nonrigid registration term in the energy function, we use the recently introduced information theoretic-based criteria (10) called the cross-cumulative residual entropy (CCRE). CCRE was shown to outperform mutual information-based registration in the context of noise immunity and convergence range (10), motivating us to pick this criterion over the MI-based cost function. The new registration functional is defined by
where cross-CRE C(I1, I2) is given by
with E(I1) = −∫+ (I1 > λ) log(I1 > λ) dλ and + = (x ; x ≥ 0). v(x) is before and μ is the regularization parameter and . denotes Frobenius norm. Using a B-spline representation of the nonrigid deformation, one need only compute this field at the control points of the B-splines and interpolate elsewhere, thus accruing computational advantages. Using this representation, we have derived analytic expressions for the gradient of the energy with respect to the registration parameters. This in turn makes our optimization more robust and efficient.
For the registration and the segmentation terms to “talk” to each other, we need a connection term and that is given by
where in is the region enclosed by , øv(C) (x) is the embedding signed distance function of the contour v(C), which can be used to measure the distance between v(C) and . The level set function ø: 2 → 1 is chosen so that its zero level set corresponds to the transformed template curve v(C). Let Edist :=dist(v(C),), one can show that where N is the normal to . The corresponding curve evolution equation given by gradient descent is then
Not only does the signed distance function representation make it easier for us to convert the curve evolution problem to the level set framework, it also facilitates the matching of the evolving curve and the transformed template curve v(C), and yet does not rely on a parametric specification of either or the transformed template curve. Note that because dist(v(C), ) is a function of the unknown registration v and the unknown segmentation , it plays the crucial role of connecting the registration and the segmentation terms.
Combining these three functionals together, we get the following variational principle for the simultaneous registration + segmentation problem:
αi are weights controlling the contribution of each term to the overall energy function and can be treated as unknown constants and either set empirically or estimated during the optimization process. This energy function is quite distinct from those used in methods existing in literature because it is achieving the Mumford-Shah type of segmentation in an active contour framework jointly with nonrigid registration and shape distance terms. We can solve for in a level set framework by representing it as a level set function ø. The unknown displacement field (represented as B-spline) and the level set function ø are updated iteratively until converges so that we can achieve the registration and segmentation simultaneously. The details of our algorithm are summarized in the following section.
Given the atlas image I1 and the unknown subject's brain scan I2, we want the segmentation result in I2. Initialize in I2 to C and set the initial displacement field to zero.
In this section, we present several examples results from an application of our algorithm. The results are presented for synthetic and real data. The first three experiments were performed in two dimensions, whereas the fourth one was performed in 3D. Note that the image pairs used in all these experiments have significantly different intensity profiles, which is unlike any of the previous methods, reported in literature, used for joint registration and segmentation. The synthetic data example contains a pair of MR T1- and T2-weighted images, which are from the Montreal Neurological Institute Brain Web site (25). They were originally aligned with each other. We use the MR T1 image as the source image and the target image was generated from the MR T2 image by applying a known nonrigid transformation that was procedurally generated using kernel-based spline representations (cubic B-spline). The possible values of each direction in deformation vary from −15 to 15 in pixels. In this case, we present the error in the estimated nonrigid deformation field, using our algorithm, as an indicator of the accuracy of estimated deformations.
Figure 2 depicts the results obtained for this image pair. With the MR T1 image as the source image, the target was obtained by applying a synthetically generated nonrigid deformation field to the MR T2 image. Notice the significant difference between the intensity profiles of the source and target images. Figure 2 is organized as follows, from left to right: the first row depicts the source image with the atlas-segmentation superposed in red, the registered source image that is obtained using our algorithm followed by the target image with the unregistered atlas-segmentation superposed to depict the amount of misalignment; the second row depicts ground truth deformation field, which we used to generate the target image from the MR T2 image, followed by the estimated non-rigid deformation field and finally the segmented target. As visually evident, the registration + segmentation are quite accurate from a visual inspection point of view. As a measure of accuracy of our method, we estimated the average, μ, and the standard deviation, δ, of the error in the estimated nonrigid deformation field. The error was estimated as the angle between the ground truth and estimated displacement vectors. The average and standard deviation are 1.5139° and 4.3211°, respectively, which is quite accurate.
Table 1 depicts statistics of the error in estimated non-rigid deformation when compared with the ground truth. For the mean ground truth deformation (magnitude of the displacement vector) in column 1 of each row, five distinct deformation fields with this mean are generated and applied to the target image of the given source-target pair to synthesize five pairs of distinct data sets. These pairs (one at a time) are input to our algorithm and the mean μ of the mean deformation error (MDE) is computed over the five pairs and reported in column 2 of Table 1.
MDE is defined as dm = ΣxiR v0(xi) / card(R), where v0(xi) and v(xi) are the ground truth and estimated displacements respectively at voxel xi. . denotes the Euclidean norm, and R is the volume of the region of interest. Column 3 depicts the standard deviation of the MDE for the five pairs of data in each row. As evident, the mean and the standard deviation of the error are reasonably small indicating the accuracy of our joint registration + segmentation algorithm. Note that this testing was done on a total of 20 image pairs (=40) because there are five pairs of images per row.
For the first real data experiment, we selected two image slices from two different modalities of brain scans. The two slices depict considerable changes in shape of the ventricles, the region of interest in the data sets. One of the two slices was arbitrarily selected as the source and segmentation of the ventricle in the source was achieved using an active contour model. The goal was then to automatically find the ventricle in the target image using our algorithm given the input data along with the segmented ventricles in the source image. Figure 3 is organized as follows, from left to right: the first row depicts the source image with the atlas segmentation superposed in black followed by the target image with the unregistered atlas segmentation superposed to depict the amount of misalignment; second row depicts the estimated nonrigid vector field and finally the segmented target. As evident from Figure 3, the accuracy of the achieved registration + segmentation is visually very good. Note that the nonrigid deformation between the two images in these two examples is quite large and our method was able to simultaneously register and segment the target data sets quite accurately.
The second real data example is obtained from two brain MRIs of different subjects and modalities, the segmentation of the cerebellum in the source image is given. We selected two “corresponding” slices from these volume data sets to conduct the experiment. Note that even though the number of slices for the two data sets is the same, the slices may not correspond to each other from an anatomical point of view. However, for the purposes of illustration of our algorithm, this is not very critical. We use the corresponding slice of the 3D segmentation of the source as our atlas-segmentation. The results of an application of our algorithm are organized as before in Figure 4. Once again, as evident, the visual quality of the segmentation and registration are very high.
Finally, we present a 3D real data experiment. In this experiment, the input is a pair of 3D brain scans with the segmentation of the hippocampus in one of the two images (labeled the atlas image) being obtained using the well-known Principle Component Analysis (PCA) on the several training data sets. Each data set contains 19 slices of size 256 × 256. The goal was then to automatically find the hippocampus in the target image given the input. Figure 5 depicts the results obtained for this image pair. From left to right, the first image shows the given (atlas) hippocampus surface followed by one cross-section of this surface overlaid on the source image slice; the third image shows the segmented hippocampus surface from the target image using our algorithm and finally the cross-section of the segmented surface overlaid on the target image slice. To validate the accuracy of the segmentation result, we randomly sampled 120 points from the segmented surface and computed the average distance from these points to the ground truth hand-segmented hippocampal surface in the target image. The hand segmentation was performed by an expert neuroanatomist. The average and standard deviation of the error in the aforementioned distance in estimated hippocampal shape are 0.8190 and 0.5121 (in voxels), respectively, which is very accurate.
In this article, we presented a novel variational formulation of the joint (nonrigid) registration and segmentation problem, which requires the solution to a coupled set of nonlinear Partial Differential Equation (PDEs) that are solved using efficient numerical schemes. Our work is a departure from earlier methods in that we presented a unified variational principle wherein nonrigid registration and segmentation are simultaneously achieved. Unlike earlier methods presented in literature, a key feature of our algorithm is that it can accommodate for image pairs having distinct intensity distributions. We presented 20 examples on synthetic and 3 real data sets along with quantitative accuracy estimates of the registration in the synthetic data case. The accuracy as evident in these experiments is quite satisfactory. Further evaluation using more hippocampal datasets will be investigated in our future work.
In our algorithm, the atlas played a fundamental role. We did not construct an atlas image, however, and the atlas hippocampal shape from a training set, but picked one arbitrarily. This would yield a biased atlas image and a biased atlas shape. Our future work will focus on creating an unbiased atlas image and shape, respectively, for use in the joint registration and segmentation algorithm. Figure 1.
Authors thank Dr. C. M. Leonard of the University of Florida Neuroscience Department for the hippocampal data sets.
1Funded in part by the NIH grant RO1 NS046812.