A statistical model is presented that combines the registration of an atlas with the segmentation of magnetic resonance images. We use an Expectation Maximization-based algorithm to find a solution within the model, which simultaneously estimates image artifacts, anatomical labelmaps, and a structure-dependent hierarchical mapping from the atlas to the image space. The algorithm produces segmentations for brain tissues as well as their substructures. We demonstrate the approach on a set of 22 magnetic resonance images. On this set of images, the new approach performs significantly better than similar methods which sequentially apply registration and segmentation. D 2005 Elsevier Inc. All rights reserved. Keywords: Registration; Segmentation; Subcortical segmentation; Bayesian modeling; Expectation -Maximization To better understand brain diseases, many neuroscience studies focus on the anatomical differences between control and diseased subjects. In order to find these differences, scientists often analyze medical images for brain structures which seem to be influenced by the disease. The analysis is frequently based on segmentations of the structures of interest that are mostly performed by human experts. However, this manual process is not only very expensive, but in addition, it increases risks related to inter-and intra-observer reliability (Kikinis et al., 1992) . Neuroscientists are keenly interested in automatic methods, which often rely on prior information, to perform this task (Collins et al., 1999; Leventon et al., 2000; Marroquin et al., 2003; Fischl et al., 2004; Pohl et al., 2004a; Ashburner and Friston, 2005) . With notable exceptions, these methods first register the prior information, i.e., an atlas, to the medical image and then segment the medical image into anatomical structures based on that aligned information. The goal of this work is to unify this process into a single Bayesian framework in order to overcome biases caused by commitment to the initial registration. When automatic segmentation methods are guided by prior information, they frequently are used to segment anatomical structures defined by weakly visible boundaries in medical images. For example, the intensity properties of the thalamus in T1-weighted magnetic resonance (MR) images are very similar to those of the neighboring white matter (Fig. 1) . Algorithms cannot rely on the MR images alone in order to distinguish these two structures. However, the ventricles, the dark structures above the thalamus, are more easily identified. In order for the ventricles to guide the detection of the boundary between the thalamus and the white matter, automatic segmentation algorithms use spatial priors (Mazziotta et al., 1995; Thompson et al., 1996) . These spatial priors capture the relationship between structures such as the fact that the ventricles are above the thalamus. As mentioned previously, most atlas-based algorithms perform registration and segmentation sequentially (Cocosco et al., 2003; Van Leemput et al., 1999; Fischl et al., 2002) . In this paper, we present an approach that exploits complementary aspects of both problems. In the example of the thalamus, segmentation simplifies the registration of the region defined by ambiguous intensity patterns between thalamus and white matter (Pohl et al., 2004a) . On the other hand, aligning an atlas aids the detection of thalamus as discussed previously. The idea of joining registration and segmentation has been utilized by boundary localization techniques using level set representation, such as Leventon et al. (2000) and Tsai et al. (2003) , which align an atlas to the subject and simultaneously estimate the shape of a structure. These methods relate both problems to each other by extending the definition of the shape to include its pose. This paper describes an integrated segmentation and registration approach related to voxel-based classification methods, which considers the anatomical structure associated with each voxel within a Bayesian framework (Wells et al., 1996; Van Leemput et al., 1999; Marroquin et al., 2003; Pohl et al., 2004a) . These voxel-based classification methods explicitly model image artifacts in order to segment large data sets without manual intervention. Voxel-based classification methods that have coupled registration and segmentation of misaligned spectral images include (Wyatt and Noble, 2002; Xiaohua et al., 2005 our purpose is to align an atlas to MR images and separate the images into anatomical structures, as also suggested by (Ashburner and Friston, 2005) . The proposed non-rigid registration framework of Ashburner and Friston (2005) performs a local deformation of the atlas space under the assumption that the atlas is globally aligned to the image space. This is an assumption which does not hold in practice and which we address in our approach. Section 2 describes the Bayesian framework that models the relationship between the unknown segmentation, the observed MR image data, and the registration parameters. The registration parameters capture the deformation between atlas and image space. The solution within our joint registration and segmentation framework is defined by a Maximum A posteriori Probability (MAP) estimation problem and is found by adopting an instance of the expectation maximization (EM) algorithm (Dempster et al., 1977; McLachlan and Krishnan, 1997) . In Section 3, we specify a hierarchical description of the registration parameters. Global registration parameters capture the correspondence between the entire atlas and image space. Structure-dependent registration parameters represent the local deformations within a region of an anatomical structure. This hierarchical model simplifies the correspondence problem between atlas and image space. We describe in Section 4 a robust implementation of our algorithm, where the registration parameters capture affine mappings. In Section 5, we perform an experiment comparing the robustness of this implementation to other EM segmenters. In this experiment, the automatic methods outline a set of 22 MR images into the major brain tissue classes as well as the thalamus, which is a structure with indistinct boundaries. When comparing the automatic results to manual segmentations, our new method significantly outperforms the other implementations. The accuracy of segmenting structures that have indistinct boundaries in MR images, using tissue classification methods, significantly depends on properly modeling image artifacts and accurately registering prior information to the subject. In this section, we develop a unified framework which performs registration, segmentation, and image artifact estimation simultaneously. In this paper, we are primarily interested in studying the dependencies between segmentation and registration. We therefore view image artifacts as nuisance parameters 1 , which, if known, would greatly simplify the segmentation of MR images. For example, nuisance parameters can represent the image inhomogeneity caused by the acquisition of the MR images or the subject specific intensity histogram of an anatomical structure. The proper modeling of these nuisance parameters is an ongoing discussion in neuroimaging (Wells et al., 1996; Leventon, 2000; Kapur, 1999; Marroquin et al., 2003; Ashburner and Friston, 2005) . Due to their complex dependencies, it is difficult to extract the nuisance parameters h and the registration parameters R from the MR images I without explicit knowledge of the unknown true segmentation. However, this problem is greatly simplified when the solution is determined within an EM framework. In this framework, the solution is formulated as an incomplete data problem marginalizing over all possible segmentations T . The remainder of this section derives an instance of the EM algorithm by modeling the dependencies between the nuisance parameters h, the registration parameters R, the hidden labelmap T , and the observed image data I . To determine h and R within this framework, we define the following MAP estimation problem: In general, this results in a system of equations for which there is no analytical solution. We simplify the problem by adding the labelmap T to the model. The unknown labelmap T assigns each voxel in the image to an anatomical structure. Combining T with the known intensity histograms of anatomical structures produces images, which then can be compared to the observed MR images I in order to determine the solution to Eq. (1). We can restate Eq. (1) by marginalizing with respect to all possible labelmaps T ðĥ h;R RÞ ¼ arg max Next, we incorporate the conditional labelmap probability P T j I ;h V ;RV ð Þ , where h V;RV ð Þare estimates of ðĥ h;R RÞ, into Eq. (2) ðĥ h;R RÞ ¼ arg max The purpose of these operations is to put Eq. (3) into a form such that we can exploit the following bound derived via Jensen's Inequality (Neal and Hinton, 1998) log E T j I;h V ;RV P h;R; The right side of Eq. (4) defines a lower bound on the objective function (the left side). The EM approach uses this lower bound to find a local maximum to the objective function. The Expectation . From the intensity pattern in panel (a), the boundary of the thalamus and the neighboring white matter is not clearly defined. 1 It is common in discussions of statistical estimation to refer to parameters, which are not of primary importance, as nuisance parameters in formulating an estimation problem. For the purposes of this work, the relevant nuisance parameters are the values of the intensity correction at every voxel. Step (E-Step) establishes a new objective function or the lower bound (McLachlan and Krishnan, 1997) The Maximization Step (M-Step) updates the approximations h V and RV with the maximum of the lower bound For many problems, Eq. (6) defines a simpler maximization problem than Eq. (3). Finding the solution to Eq. (6) is still complicated, as we simultaneously estimate the nuisance parameters h and registration parameters R. We will make certain simplifying assumptions. Before we do so, we rewrite Eq. (6) by applying Bayes' rule and dropping terms that do not depend on The optimization procedure decomposes nicely as a consequence of the following independence assumptions (see Fig. 2 ): First, we assume that I is independent of R conditioned on T and h. We therefore characterize each anatomical structure with a stationary intensity distribution which is not influenced by the mapping between atlas and image space such as described by the models (Wells et al., 1996; Van Leemput et al., 1999; Kapur, 1999; Ashburner and Friston, 2005) . Based on the same argument, we also assume R is independent of h conditioned T . Finally, we assume independence of h and T . 2 Note that in the graphical model shown in Fig. 2 , the link joining the registration parameters R with the labelmap T could be reversed by the use of Bayes' rule, which would yield an equivalent model (this conditioning of T on R is used below in the development of the EM iteration). Based on the three independence assumptions, Eq. (7) simplifies to We demonstrate the utility of the independence assumptions in the experiment of Section 5 when we segment the MR images of 22 subjects into white matter, gray matter, cerebrospinal fluid, and cerebrospinal structures. In its current form, Eq. (8) contains the expected value, which sums over all possible settings of T . To further simplify the problem, we assume spatial independence of T and independence between R and h on the basis of the previous discussion. Eq. (8), which characterizes the EM approach, is restated as the sum over all voxel locations x of the expected value with respect to indicator random vector T x a e 1 ; . . . ;e N f g . The vector e a is zero at every position but a, where its value is one. For example, if T x ¼ e a then voxel x is assigned to the anatomical structure a. If we now assume that the two conditional probabilities P I j T ;h ð Þand P Rj T ð Þ are defined by the product of the corresponding conditional probabilities over all the voxels in the image space, then we can rewrite Eq. (8) as The weights, which will be calculated in the E-Step instead of Eq. (5), capture the posterior probability of the structure a being present at voxel x and are now defined as As each factor in Eq. (9) only depends on either R or h, the M-Step updates the approximations of the inhomogeneities and the registration parameters separately The EM algorithm iterates between the E-Step (Eq. (10)) and the M-Steps (Eqs. (11) and (12)) until the lower bound function converges to a local extrema, which is guaranteed by the EM framework if the iteration sequence has an upper bound (McLachlan and Krishnan, 1997) . Based on Eq. (11), the algorithm adjusts the registration parameters so that the atlas matches the weights defined by the E-Step. However, the weights might mislead the algorithm in cases where the alignment between atlas and image coordinate system is difficult. For example, the MR image of Fig. 3 shows a subject with an unusual head position. In this case, the weights wrongly assert that the neck is part of the brain. This causes the algorithm to adjust the nuisance parameters so that intensity histogram of In this model, the observed image data I is directly linked to the labelmap T and the nuisance parameters h. In addition, the labelmap T is linked to the registration parameter R. 2 Note, if we incorporate a nuisance parameter model into our algorithm that does not fulfill these independence requirements, e.g., we use spatial varying intensity distributions as suggest by Fischl et al. (2004) , then we perform the maximization for the nuisance and registration parameters sequentially. The resulting algorithm belongs to the class of generalized-EM algorithms, which have the same convergence properties as EM. the neck is equivalent to the ones inside the brain. The registration parameters scale the atlas of the brain to cover both brain and neck. The algorithm, therefore, converges to a suboptimal solution, which does not separate the neck from the brain. In summary, this section defined the problem of mapping the atlas to the image space as a MAP estimation problem. Finding the solution to this problem is generally too difficult, so we iteratively solve the simpler lower bound optimization problem of Eq. (4) in an instance of an EM algorithm. Even computing the solution to the optimization problem is generally too complicated as we have to simultaneously maximize the underlying objective with respect to the nuisance parameters h and the registration parameters R. We therefore split Eq. (4) into two separate estimation problems for h and R by making certain simplifying independence assumptions. The resulting EM approach is composed of an E-Step determining the W of Eq. (10) and the M-Step, which solves the estimation problems of Eqs. (11) and (12). While the solution to Eq. (12) has been discussed in detail in the literature (Wells et al., 1996; Leventon, 2000; Kapur, 1999; Marroquin et al., 2003; Ashburner and Friston, 2005) , the next section focuses on determining the registration parameters of Eq. (11). To solve the estimation problem of Eq. (11), we define the registration parameters R, the conditional log probability log P T j R ð Þ, and the prior probability P R ð Þ. We model the parameters R as a hierarchical registration framework, which distinguishes between global-and structure-dependent deformations. We then apply the registration framework to the estimation problem and determine its solution with another optimization algorithm. The hierarchical registration parameters R capture the correspondence between atlas, brain, and structures within the brain. The mapping of the atlas to the image space is performed by an interpolation function r R; x ð Þ, which maps voxel x into the coordinate system defined by R. r R;I ð Þ can be rigid, affine, or more general non-rigid. As a rigid registration can be interpreted as a special case of an affine transformation, we focus our discussion on affine and non-rigid registration methods. The parameters R can be structure-dependent or -independent. The structure-independent parameters capture the correspondence between atlas and image space without knowledge of the underlying brain anatomy. If constrained to affine interpolation, the degrees of freedom of R are too low to capture the characteristics of individual brain structures (Pohl et al., 2002; Srivastava et al., 2004) . The alternative is a more general non-rigid framework, which often has problems aligning structures with weakly visible boundaries (Pohl et al., 2004b) . Instead, we favor an approach that registers each anatomical structure based on the specific requirements of that structure. Our structure-dependent registration parameters treat the relationship between the atlas and image space for each structure independently. However, this independence assumption does not hold in practice. Our registration problem is mostly characterized by the structure-independent misalignment between the coordinate system of the image space (represented by the position of the head) and the atlas coordinate system. Thus, the difference between different structure-dependent parameters should be relatively small. We therefore propose a registration framework that consists of structure-dependent and -independent parameters. We model dependency across structures with a hierarchical registration framework R ¼ R G ;R C ð Þas shown in Fig. 4 . R G are the global registration parameters, which describe the nonstructure-dependent deformations between atlas and image. The Fig. 3 . For this MR image, our approach achieved a suboptimal result as it failed to differentiate the neck from the brain (outlined by the line). We refer to the text for a detailed discussion. Fig. 4 . The graph represents the hierarchical affine registration framework developed throughout this section. In this approach, the registration parameter R G captures the global deformations of the atlas coordinate system to the image space. R a represents the structure-specific alignment of structure a. The atlas is aligned to the structure-specific coordinate system in the image space by applying the interpolation function r R G ;R a ; I ð Þ . structure-dependent parameters R C ç R 1 ; . . . ;R N ð Þare defined in relation to R G and capture the residual structure-specific deformations that are not adequately explained by R G . In the remainder of this work, we refer to R a , the ath entry of R C , as the registration parameters specific to structure a with a Z {1, . . . , N}. The similarity between structure-specific parameters is encoded in Eq. (11) through the prior probability P R ð Þ. Section 4 defines R explicitly and models P R ð Þ with a simple Gaussian distribution. This model enforces tight bounds on R C and weak constraints on R G . To further specify the model, we now define the conditional probability P T x ¼ e a j R ð Þ , which is the spatial prior probability of structure a in the image space. In our model, this prior is composed by an atlas of the spatial distributions across all structures. The spatial distribution for a single structure a is represented in our model by f a (I). Fig. 5 shows examples of f a (I) for the right and left thalamus as well as the remaining gray matter in the images space. These examples were generated by computing the overlap of between the segmentations of structures of interest according to (Warfield et al., 2001 ). Blue indicates low and red high probability of the structure at that voxel location. f a (I) is defined in the coordinate system of the atlas space, which is in general different from the image space. We align the atlas to the image space by making use of the interpolation function r R G ;R a ;I ð Þ . The conditional structure probability is then defined as which is the conditional probability of the presence of an individual structure normalized by the spatial prior probability of all structures. The normalization across all structures is necessary as the coordinate system of each structure is characterized by the structure-dependent registration parameters R a . Unlike global affine registration methods, this results in structure-dependent coordinate systems that are not aligned with each other. In other words, multiple voxels in the atlas space can be mapped to one location in the image space. Substituting Eq. (13) into Eq. (11) changes the MAP estimation problem to The spatial conditioned probabilities f a r R G ;R a ; x ð Þ ð Þare not sufficient to account for all occurrences within the model. For example, the product W x a ð Þ I log f a r R G ;R a ; x ð Þ ð Þ is not defined in the case where W x a ð Þ 6 ¼ 0 and f a r R G ;R a ; x ð Þ ð Þ¼0. To overcome this problem, one can simply add a uniform distribution to the spatial conditioned probabilities by replacing log f a r R G ;R a ; x ð Þ ð Þ with log f a r R G ;R a ; x ð Þ ð Þ þ Þ. Note, that our approach is relatively insensitive to the exact value of the . The objective function Q(I) measures the disagreement between f a (I) and W a ðÞ (see Fig. 6 ). We show in Appendix A that the MAP estimation is composed by the Kullback -Leibler (KL) divergence Fig. 5. Examples of space conditioned probabilities of anatomical structures. In our framework, these probabilities are represented by f a (I). Blue indicates low and red high probability with respect to presence of a structure at that voxel location. Fig. 6 . We applied the EM approach to the MR image (a) and the misaligned spatial atlas of white matter and gray matter in panel (b). The result of the analysis is shown in panel (c), which is the weights W of the E-Step, and (d), which is the initial objective function Q(I) of the M-Step. The impact of the misalignment between atlas and image space is especially apparent in the region of the ventricles, whose location are pointed out in the image space and in the atlas space by arrowheads. D(I||I) with respect to the labelmap T and the prior P R ð Þ of the registration parameters R: Ignoring the log prior log P R ð Þ, the solution of the estimation problem is the registration parameters R, which minimize the difference between the weights W T (see Appendix A) and the spatial prior P T j R ð Þ with respect to all labelmaps T , in the KL sense. The example of Fig. 6 shows how the misalignment between atlas and image space impacts the analysis of the EM approach. The intensity of the images corresponds to the value of the underlying function with white indicating high and black low values. In this example, we especially focus on the region around the ventricles. Arrowheads point out their location in the MR image in (a) and the atlas space in (b), which shows the combined spatial atlas of white matter and gray matter. The misalignment between (a) and (b) impacts the weights W of the brain in (c) as it detects the ventricles twice-once according to the location of the image and once according to the spatial prior. The initial objective function Q(I) defined by the M-Step in (d) detects this ambiguity as it measures the disagreement between spatial prior of the brain f a (I) in (b) and the weights W a I ð Þ in (c). The objective function shows large disagreement in the area around the location of the ventricles in the image space as (b) assigns this region to the brain and (c) does not. In order to resolve this disagreement, the EM approach aligns the atlas to the image space as we will observe later in the example of Fig. 7 . In summary, we developed an EM approach based on the newly defined registration parameters R. The E-Step calculates the weights W at voxel x of Eq. (10) based on the aligned spatial priors f a r R G V ;R a V ;I ð Þ ð Þ , intensity I, and the approximation of the nuisance parameters hV: The M-Step updates the approximation of the nuisance parameters hV and registration parameters RV based on the current weights W. These parameters will be described in more detail in the next section. The update rule for R is defined by Eq. (14) as the maximum of objective functions Q(I) The solution to this update rule depends on the underlying model with respect to R G and R a , which is discussed in detail in the next section. This section describes an implementation of the previously derived EM approach. We give an example for the interpolation function r(I, I, I), the corresponding registration parameters R, the prior of the registration parameter P R ð Þ, and the maximization algorithm to solve the MAP estimation problem defined in Eq. (14). The interpolation function r(I, I, I) of Section 3 can represent a variety of mapping approaches. For simplicity, we choose a restricted affine interpolation function. The global and structure-specific parameters are defined as . ., N}, where G is the index of the global parameters and {1,. . ., N} represents the list of indices for the structure-specific parameters. The parameters capture displacement Y t z a R 3 , rotation Y r z a R 3 , and scaling Y s z a R 3 . The mapping is defined by the function r I; I; I ð Þ: where A Rz is an affine transformation matrix based on the parameter setting R z . This framework makes no assumptions about the correspondence between the atlas and the image space by uniformly distributing the global registration parameter R G . As opposed to R G , the structure or class specific parameters R C ç R 1 ; . . . ;R N ð Þ describe the residual of structure-specific deformations that are not well explained by R G . In general, our model should penalize large deviations of R C from the expected mean, which is approximated by the average structure-specific registration parameters of the training data. We enforce this penalty by modeling the probability density function of R C as a Gaussian distribution N l RC ;Y RC À Á with structure-dependent mean l RC and variance Y RC based on the mapping parameters of the training data. We choose a Gaussian distribution, as small variance Y RC discourages large deformations from the mean l RC . In addition, Gaussian distributions simplify the calculations in the M-Step (McLachlan and Krishnan, 1997). Based on the previous modeling assumptions, the objective function Q(I) of Eq. (14) is defined as To find a solution to the MAP problem defined by Q(I), we decouple the search for R G and R S as their dependencies can cause instability. It is generally difficult to determine a closed form solution for this problem so that we estimate the solutions through the Powell's method (Press et al., 1992) . The second MAP estimation problem (Eq. (12)) of the M-Step updates the approximation of the nuisance parameters In our implementation, the nuisance parameters h are defined by the image inhomogeneity model of Wells et al. (1996) . We choose this model as it is quite simple but also flexible enough for the given domain. Note that other models, such as Leventon (2000); Kapur (1999); Marroquin et al. (2003); and Ashburner and Friston (2005) , are appropriate choices for the nuisance parameters. The choice of model influences the behavior of our algorithm. In this paper, however, we will not further investigate the impact of these different models on our algorithm as we are primarily interested in studying the dependencies between segmentation and registration. One of Wells' models defines the likelihood P I jT x ¼ e a ;h ð Þof Eq. (12) by the following Gaussian distribution: (l a , Y a ) represents the mean and variance of the intensity distribution of the anatomical structure a. As a consequence of this definition, Eq. (12) reduces to a system of linear equations. In Appendix B, we explain in further detail the module of the image inhomogeneities and derive the solution of this system of linear equations. The solution is the product of the weighted residual Þand a low pass filter represented by the matrix H. The E-Step, characterized by Eq. (15), calculates the weights The entire implementation is summarized by the pseudo code below: Algorithm 1. Segmentation and Registration repeat E-Step: Update soft assignment according to Eq. (17) The last row shows the objective function Q(I) corresponding to the segmentations above. Red indicates high and blue low disagreement between the atlas space and the weights of the M- Step. Initially, the area around the ventricles and cerebrospinal fluid close to the neck are highlighted in red, as they are correctly segmented regions that disagree with their location within the atlas space. As the method progresses, the overall accuracy of the registration as well as segmentation increases. The disagreement captured by the objective function in the final iteration is mostly caused by the smooth boundaries of the spatial priors. In summary, this section described an implementation of the joint registration and segmentation algorithm. Affine registration parameters represent the global and structure-specific deformations between atlas and image space. The nuisance parameters capture the inhomogeneity of the MR images caused by the acquisition of these images. The next section discusses the accuracy of our implementation in outlining the thalamus. The experiment described in this section compares the accuracy of four EM implementations which differ in how they map the atlas to the patient. The first approach (EM-NonRigid) maps the atlas to the patient using an intensity based non-rigid registration approach and then runs our EM implementation without registration parameters (Pohl et al., 2004a) . The second approach (EM-Affine) is similar to EM-NonRigid but uses the affine mapping method by (Warfield et al., 2001) . The remaining two approaches are based on our new algorithm, which solves the registration and segmentation problem simultaneously. EM-Sim-Affine uses the same global affine registration model as EM-Affine and EM-Sim-Hierarchical incorporates the hierarchical registration approach of Section 3. All four methods utilize the same atlas, which captures the spatial distributions of the anatomical structures. The training data of the atlas consist of semi-automatic as well as manual segmentations of up to 80 different subjects 3 . The atlas is generated through the approach suggested by Warfield et al. (2001) , which aligns the training subjects to a preselected training case and then measures the overlap between the corresponding segmentations. Examples of the atlas are shown in Fig. 5 . Each method segments 22 test cases (1.5-T MRI system (GE Medical Systems, Milwaukee), three-dimensional Fourier transformation spoiled gradient-recalled acquisition sequence, echo time (TE) = 5 ms, repetition time (TR) = 35 ms, repetition = 1, nutation angle = 45, field of view = 24 cm, acquisition matrix = 256 Â 256 Â 124, voxel dimension = 0.9375 Â 0.9375 Â 1.5 mm) into the three brain tissue classes -white matter, gray matter, and cerebrospinal fluid -and further parcellates gray matter into the right and left thalamus. We evaluate the accuracy of the approaches by measuring the agreement of the automatic segmentations of the thalamus to manual ones, which we view as ground-truth. Note that our EM approach segments anatomical structure by considering all neighboring structures. Thus, the accuracy of segmenting the thalamus greatly depends on the precise segmentation of white matter, the remaining gray matter, and the cerebrospinal fluid. We focus on the thalamus because it is a challenging structure for registration and segmentation. Many fiber tracks pass through the thalamus so that its intensity properties in MR images are very similar to the neighboring white matter (Fig. 8) . Intensity based registration methods have difficulties aligning this structure because of its weakly visible boundary. In our experience, EM without spatial priors will fail and EM with spatial priors heavily relies on these priors. Thus, registration errors greatly influence the segmentation quality of our three methods. To measure the quality of the automatic generated results, we compare them to the manual segmentations using the volume overlap measure DICE (Dice, 1945) . The results of the experiment are summarized in the Fig. 9 . The graph shows the mean and standard error of the DICE measure for the three algorithms in the 22 cases. For the thalamus, EM-NonRigid performed worst (82.35 T 1.2% À average DICE score T standard error) because the intensity based registration method is too unreliable for structures with smooth boundaries. The method generally overestimates white matter and underestimated the thalamus in this region, which can be also observed by the segmentation in Fig. 8 (EM-NonRigid) . EM-Affine (87.28 T 1.2%) performs much better than EM-NonRigid, but the method is sensitive towards initial misalignments. In Fig. 8 (EMAffine), this sensitivity causes a vertical offset in the segmentation generated by EM-Affine. Note that the average DICE score of EM-Affine for the thalamus is already excellent as we experienced large disagreement even between different expert segmentations. However, EM-Sim-Affine and EM-Sim-Hierarchical, the approaches based on our new framework, achieve a higher average score than the other two methods. EM-Sim-Affine is more robust (88.65 T 0.4%) than EM-Affine even though both methods use the same affine registration framework. This empirically demonstrates the utility of performing registration in an EM segmenter. EM-Sim-Affine is constrained by the dependency between segmentation and registration. In comparison to EM-Affine, these constraints further reduce the space of possible solutions, which simplifies the search for the optimal solution. In addition, the spatial priors of the structures are directly mapped to the segmentation model. In contrast, EM-Affine and EM-NonRigid align an MR image in the atlas space to the image of the patient. They use the resulting deformation map to align the spatial priors, which inherently increases the risk of systematic biases in the model. EM-Sim-Hierarchical (89.40 T 0.3%), however, significantly outperforms the other three approaches. The hierarchical registration framework has greater flexibility than the single global affine registration models of EM-Affine and EM-Sim-Affine. We believe that the structure-specific registration framework of EM-Sim-Hierarchical better incorporates the anatomical knowledge gained from the segmentation domain than the continuous model of EM-Sim-Affine. We end the discussion by further investigating the previously made observation with respect to the ventricles. As pointed out in Section 3, the ventricles are very important in guiding the joint registration and segmentation method. In order to validate this statement, we repeat the previous experiment. This time gray matter, white matter, and the ventricles are combined into one structure. This change in our model causes the registration cost function to ignore the misalignment of the location of the ventricles between atlas and image space (see also Fig. 6 ). As a result, the average DICE score for the thalamus reduces by 5.2% and the standard error increases by a factor of 2.5. We therefore conclude that the accuracy of our new approach greatly depends on anatomical structures that are easily identifiable in the segmentation problem. We have demonstrated that our new approach EM-SimHierarchical and EM-Sim-Affine achieve higher accuracy than the methods we denoted as EM-Affine and EM-NonRigid as a consequence of the fact that our approach directly maps the spatial priors of the structures to the segmentation model. In contrast, EMAffine and EM-NonRigid align an MR image in the atlas space to the image of the patient (Pohl et al., 2004a) , using the resulting deformation map to align the spatial priors. This inherently increases the risk of systematic biases in the model. Furthermore, EM-Sim-Hierarchical performs significantly better then the other three methods as the hierarchical registration model utilizes well the anatomical knowledge gained from the segmentation domain. We presented a statistical framework combining atlas registration and segmentation of MR images. Unlike other voxel-based classification methods, our framework models these problems as a single MAP estimation problem. We implemented the framework as an instance of an EM algorithm using a hierarchical affine mapping approach for anatomical structures. Our approach was validated by automatically segmenting 22 sets of MR images into the major brain tissue classes and the thalamus, a structure with indistinct boundaries. Using manual segmentations, we then compared our results to other EM implementations which sequentially register and segment. The study empirically demonstrated the utility of combining registration and segmentation in an EM framework. In addition, our hierarchical registration approach achieved a significantly higher score than the single continuous registration models.
Pohl et al. REF used an expectation maximization (EM) type algorithm again with shape priors to perform segmentation.
A Bayesian model for joint segmentation and registration
NEUROIMAGE
Computer Science Medicine