Unbiased diffeomorphic atlas construction for computational anatomy

Construction of population atlases is a key issue in medical image analysis, and particularly in brain mapping. Large sets of images are mapped into a common coordinate system to study intra-population variability and inter-population differences, to provide voxel-wise mapping of functional sites, and help tissue and object segmentation via registration of anatomical labels. Common techniques often include the choice of a template image, which inherently introduces a bias. This paper describes a new method for unbiased construction of atlases in the large deformation diffeomorphic setting. A child neuroimaging autism study serves as a driving application. There is lack of normative data that explains average brain shape and variability at this early stage of development. We present work in progress toward constructing an unbiased MRI atlas of two year of children and the building of a probabilistic atlas of anatomical structures, here the caudate nucleus. Further, we demonstrate the segmentation of new subjects via atlas mapping. Validation of the methodology is performed by comparing the deformed probabilistic atlas with existing manual segmentations.

Unbiased Diffeomorphic Atlas Construction for Computational Anatomy S. Joshi a,b,, Brad Davis a,b Matthieu Jomier c Guido Gerig b,c aUniversity of North Carolina, Department of Radiation Oncology bUniversity of North Carolina, Department of Computer Science cUniversity of North Carolina, Department of Psychiatry Abstract Construction of population atlases is a key issue in medical image analysis, and par-ticularly in brain mapping. Large sets of images are mapped into a common coordi-nate system to study intra-population variability and inter-population differences, to provide voxel-wise mapping of functional sites, and help tissue and object seg-mentation via registration of anatomical labels. Common techniques often include the choice of a template image, which inherently introduces a bias. This paper de-scribes a new method for unbiased construction of atlases in the large deformation diffeomorphic setting. A child neuroimaging autism study serves as a driving application. There is lack of normative data that explains average brain shape and variability at this early stage of development. We present work in progress toward constructing an unbiased MRI atlas of two year of children and the building of a probabilistic atlas of anatomical structures, here the caudate nucleus. Further, we demonstrate the segmentation of new subjects via atlas mapping. Validation of the methodology is performed by comparing the deformed probabilistic atlas with existing manual segmentations. Preprint submitted to Neuroimage 21st June 2004 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript Key words: Computational Anatomy, Brain Atlases, Registration, Image Segmentation 1 Introduction Since Broadman,1909 the construction of brain atlases has been central to the understanding of the variabilities of brain anatomy. More recently, since the advent of modern computing and digital imaging techniques intense research has been directed towards the development of digital three dimensional atlases of the brain. Most digital brain atlases so far are based on a single subject's anatomy [29,16]. Although these atlases provide a standard coordinate sys-tem, they are limited because a single anatomy cannot faithfully represent the complex structural variability between individuals. A major focus of com-putational anatomy has been the development of image mapping algorithms [33,24,10,27] that can map and transform a single brain atlas on to a pop-ulation. In this paradigm the atlas serves as a deformable template[12]. The deformable template can project detailed atlas data such as structural, bio-chemical, functional as well as vascular information on to the individual or an entire population of brain images. The transformations encode the variability of the population under study. A statistical analysis of the transformations can also be used to characterize different populations [5,30,17]. For a detailed review of deformable atlas mapping and the general framework for computa-tional anatomy see [32,13]. One of the fundamental limitations of using a single anatomy as a template is the introduction of a bias based on the arbitrary Corresponding author. Email address: [e-mail] (S. Joshi). 2 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript choice of the template anatomy. Thompson and Toga [31] very elegantly address this bias in their work by mapping a new data set on to every scan in a brain image database. This approach addresses the bias by in effect forgoing the formal construction of a representative template image. Although this framework is mathematically elegant and power full it results in a computationally prohibitive approach in which each new scan has to be mapped independently to all the data sets in a database. This is analogous to comparing each subject under study to every previously analyzed image. As brain image databases grow the analysis problem grows combinatorially. In more recent and related work Avants and Gee [2] developed an algorithm in the large deformation diffeomorphic setting for template estimation by av-eraging velocity fields. Most other previous work [?,3] in atlas formation has focused on the small de-formation setting in which arithmetic averaging of displacement fields is well defined. Guimond et. al develop an iterative averaging algorithm to reduce the bias[?]. In the latest work of [3], explicit constraints requiring that the sum of the displacement fields add to zero is enforced in the proposed atlas con-struction methodology. These small deformation approaches are based on the assumption that a transformations of the form h(x) = x+u(x), parameterized via a displacement field, u(x), are close enough to the identity transformation such that composition of two transformations can be approximated via the addition of their displacement fields: h1 h2(x) x + u1(x) + u2(x) . (1) 3 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript The focus of this paper is a development of a methodology that simultane-ously estimates the transformations and an unbiased template, in the large deformation setting. The method developed herein does not assume the above approximation of Equation 1 and build atlases of populations with large geo-metric variability. 2 Methods Given a collection of anatomical images a natural problem is the construction of a statistical representative of the population. If the data associated with the population under study can be easily parameterized by a flat euclidean space, classical statistical methods of simple averaging can be applied to generate such a representative. The imaging data under the Gaussian assumption can be easily represented as member of a flat space. The image can be represented as a member of very large dimensional euclidean space (RN, where N is the number of voxels in the image). Alternatively, using appropriate interpolation assumptions the image can be assumed to be a square integrable function, that is a member of the (flat) Hilbert space L2( ) where is the underlaying coordinate space, usually a compact subset of R3. The geometric variability of the anatomy itself usually cannot be represented by elements of a flat space. If the geometry of the underlying anatomy can be adequately represented by a finite number of landmarks, representative tem-plate landmark configuration can be estimated using the Procrustes method pioneered by Kendal [20] and championed by Bookstein [4]. The study of anatomical shape is inherently related to the construction of transformations of the underlying coordinate space that map one anatomy to another. Various 4 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript transformation groups of R3 have been studied for understanding anatomical geometry. These groups vary in dimensionality from simple global transla-tions (R3) and rigid rotations (SO(3)) to the infinite dimensional group of diffeomorphisms (H) [23]. In this paper we address the problem of anatomical template construction as the joint estimation of the most representative image and the associated anatomical geometry given a database of brain images. 2.1 Averages in metric spaces In this paper, the problem of building an anatomical template is posed as a statistical estimation problem. For anatomical representations in which the underlying geometry is parametrized as a Euclidean vector space, training data can be represented as a set of vectors x1, · · · , xN in a vector space V . In the small deformation elastic image mapping setting this is assumed to be true, as the deformations are assumed to be close enough to the identity mapping. Under this assumption the displacement vector fields parameterizing the transformations can be assumed to be elements of the Hilbert space of square integrable functions L2( ). In a vector space, with addition and scalar multiplication well defined, an av-erage representation of the training set can be computed as the linear average μ = 1 N XN i=1 xi . Linear averaging cannot be directly applied to the large deformation setting as under the large deformation model the space of transformations is not a vector space but rather the infinite dimensional group H of diffeomorphisms 5 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript of the underlying coordinate system . In the group of diffeomorphisms, the addition of two diffeomorphisms is not generally a diffeomorphism and, hence, a template based on linear averaging of transformations is not well defined. In this paper we extend the notion of averaging to general metric spaces first proposed by Frechet [9]. For a general metric space M, with a distance d : M × M ! R, the intrinsic mean for a collection of data points xi can be defined as the minimizer of the sum-of-squared distances to each of the data points. That is μ = argmin x2M XN i=1 d(x, xi)2 . In our previous work [7,8], we have used these concepts to extend first and second order statistical analysis to finite dimensional Riemannian Manifolds for statistical analysis of medial representations of objects. In this paper, we apply this approach to the construction of large deformation diffeomorphic templates. This work builds heavily on the mathematical metric theory of diffeomorphisms developed by Miller and Younes [22]. 2.2 Rigid Template Estimation To exemplify the template estimation problem first, consider a collection of N images Ii(x) of the same anatomy acquired with unknown rigid registra-tion. For such a collection of images the variability is the noise introduced by the imaging modality and the rigid pose of the anatomy parameterized by rigid translations and rotations ((Ti, Si) 2 R3 × SO(3)). The optimal tem-plate is one that requires the ‘minimum amount of transformation' to each of the anatomical images and which is most likely given the imaging modality 6 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript noise model. If the imaging modality is MRI than we can assume an additive Gaussian noise. The template estimation problem can then be stated as the quadratic minimization problem ˆI(x) = argmin I(x),(Ti,Si) XN i=1 Z (Ii(Six + Ti) − I(x))2dx + XN i=1 D(Ti, Si)2 , (2) where D(Ti, Si) is the metric on the space of rigid transformations from the identity transformation given by D2(Ti, Si) = ||Ti||2 + || log(Si)||2 . In the above expression the log, is the matrix log which for a rotation ma-trix can be easily calculated using the Rodriguez formula. A straight forward extension of the well known Procrustes method minimizes Equation 2. above [34]. Although the above methodology has been extensively studied for rigid trans-formations, the concept can be readily extended to general transformations. Given a metric on a group of transformations, the template construction prob-lem can be stated as that of estimating an image ˆI that requires the minimum amount of deformation to transform into every population image Ii. More pre-cisely, given a transformation group S with associated metric D : S ×S ! R, along with an image dissimilarity metric E(I1, I2), we wish to find the image ˆI such that {ˆh i, ˆI} = argmin hi2S,I XN i=1 E(Ii hi, I)2 + D(e, hi)2 (3) where e is the identity transformation. 7 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript 2.3 Large deformation Diffeomorphic Template Estimation Now consider a collection of N MRI images of different anatomies acquired by the same imaging modality in rigid registration. For such a collection of im-ages the variability is a result of not only of the noise of the imaging modality but also the inherent biological variability of the geometry of the underlying anatomy. In sharp contrast to the above setting, to study the geometry of the underlying anatomy, the finite dimensional matrix group of rigid transfor-mations must be replaced by the infinite dimensional analogue, the group of diffeomorphisms H. Let R3 be the underlying co-ordinate system of the template image. Let i R3, i = 1, · · · ,N be the coordinate systems of the images Ii. The problem of estimating the most representative template image can be stated as the estimation of an image ˆI, an associated independent coordinate system , which requires the least deformation represented by diffeomorphisms hi(x), to best match each of the input images. This framework is depicted in Figure 2. We apply the theory of large deformation diffeomorphisms [19,6,24] to generate deformations hi 2 H that are solutions to the Lagrangian ODEs d dthi(x, t) = vi(hi(x, t), t), t 2 [0, 1]. The transformations hi are generated by integrating velocity fields vi forward in time. Inverse transformations h−1 i are generated by integrating the negative velocity fields ˜ vi backward in time. The relationship between vi and ˜ vi is simply given by vi(·, t) = −˜ vi(·, (1−t). For a single trans-formation, h, this relationship is shown in Figure 1. The location y = h(x, 1) is described in terms of the forward integration of the velocity field v starting from the location x. Similarly, x can be described in terms of integrating the 8 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript negative velocity field ˜v backward in time starting at y. We induce a metric on the space of diffeomorphisms by using a Sobolev norm via a partial differential operator L on the velocity field v. Let h be a dif-feomorphism isotopic to the identity transformation e. We define the squared distance D2(e, h) as D2(e, h) = min v Z 1 0 Z ||Lv(x, t)||2 dxdt subject to h(x) = Z 1 0 v(h(x, t), t) dt. The distance between any two diffeomorphisms is defined by D(h1, h2) = D(e, h−1 1 h2). This distance satisfies all of the properties of a metric [22]. Namely it is non-negative, symmetric, and satisfies the triangle inequality. D is trivially non-negative. Symmetry follows from the fact that h−1 is generated by integrating backward in time the negative of the velocity field that generates h. Hence the minimizer of Equation 2.3 below is the same for both h and h−1, implying that D(e, h) = D(e, h−1). Miller et al. [22] give a detailed discussion of D and show that it satisfies the triangle inequality. Having defined a metric on the space of diffeomorphisms, the minimum energy 9 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript Figure 1. Velocity Field Integration Figure 2. Template Construction Framework template estimation problem (Equation 3) is formulated as {ˆh i, ˆI} = argmin hi,I XN i=1 E(Ii hi, I)2+ Z 1 0 Z kLvi(x, t)k2 dxdt subject to: hi(x) = Z 1 0 vi(hi(x, t), t) dt. (4) 10 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript The above problem of estimating a template image ˆI that is the best rep-resentative for a population of N anatomical images {Ii}Ni =1 is depicted in Figure 2. Throughout this paper we use the squared error dissimilarity metric but other metrics such as the Kullback-Leibler divergence can also be used [21]. Un-der the squared error dissimilarity measure the template estimation problem becomes {ˆh i, ˆI} = argmin hi,I XN i=1 Z (Ii (hi (x)) − I(x))2 dx + Z 1 0 Z kLvi(x, t)k2 dxdt. (5) This minimization problem can be simplified by noticing that for fixed trans-formations hi, the ˆI that minimizes Equation 5 is given by ˆI (x) = 1 N XN i=1 Ii (hi (x)) . (6) That is, ˆI is the voxel-wise arithmetic mean of the deformed images Ii (hi (x)). Note that the method for computing ˆI from {Ii} is determined by the image dissimilarity metric used. Other image dissimilarity metrics would imply dif-ferent methods for computing ˆI. Combining Equations 5 and 6 results in ˆh i = argmin hi XN i=1 Z 0 @Ii (hi (x)) − 1 N XN j=1 Ij (hj (x)) 1 A 2 dx + Z 1 0 Z kLvi(x, t)k2 dxdt. (7) Note that the solution to this minimization problem is independent of the ordering of the N images. 11 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript This template construction framework produces transformations ˆh i such that ˆh i : ! i. Since each ˆh i is a diffeomorphism, its inverse ˆh −1 i : i ! exists and can be calculated by integrating the negative velocity fields backward in time (Figure 1). Image to image correspondences can be computed from these transformations using the composition rule ˆh i,j = ˆh j ˆh −1 i : i ! j . (8) 3 Inverse Consistent Image Registration When the template construction framework presented in the previous section is applied to two images the result is an inherently inverse consistent image registration algorithm-no correction penalty for consistency is required. A registration framework is inverse consistent if image ordering does not af-fect the registration result. Many image registration algorithms are not inverse consistent because their image dissimilarity metrics are computed in the co-ordinate system of one of the images being registered. The choice of such a reference image can bias the result of the registration. Inverse consistent reg-istration is desired when there is no a priori reason to choose one image over another as a reference image. Previous work ([14,1]) has introduced methods for computing approximate inverse consistent registrations by applying inverse consistency constraints on intermediate incremental transformations. 12 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript For two images I1 and I2, Equation 7 reduces to {ˆh 1,ˆh 2} = argmin h1,h2 1 2 Z (I1 (h1 (x)) − I2 (h2 (x)))2 dx + Z 1 0 Z kLv1(x, t)k2 dxdt + Z 1 0 Z kLv2(x, t)k2 dxdt. The transformations h1 and h2 map to 1 and 2 respectively. Using the composition rule (Equation 8), we define the transformations h1,2 : 1 ! 2 = h2 h−1 1 and h2,1 : 2 ! 2 = h1 h−1 2 . In other words, h1,2 is a transformation from I1 to I2 and h2,1 is a transformation from I2 to I1. This method is inverse consistent since h1,2 h2,1 = h2,1 h1,2 = e, the identity transformation. 4 Implementation In this paper we present results based on a greedy fluid flow algorithm and are currently working on implementing the full space time optimization based on the Euler-Lagrange equations derived in [22]. Following the greedy fluid algorithm of propagating templates described in [24], we approximate the solution to the minimization problem in Equation 7 using an iterative greedy method. At each iteration k, the updated transformation hk+1 i , for each image Ii, is computed using the update rule hk+1 i = hki x + "vk i (x) . hki and vk i are the current estimated transformation and velocity for the ith image, and " is the step size. In other words, each final transformation hi is built up from the composition of k transformations. The velocity vk i for each iteration k is computed as follows. First, compute the 13 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript updated template estimate ˆIk(x) = 1 N XN i=1 Ik i (x), where Ik i = Ii(hki (x)) is the ith image deformed by hki . Next, define force functions Fk i (x) = − h Ik i (x) − ˆIk(x) i rIk i (x). This is the variation of the image dissimilarity term in Equation 7 with respect to hi. The velocity field vk i is computed at each iteration by applying the inverse of the differential operator L to the force function, i.e. vk i (x) = L−1Fk i (x), where L = r2 + rr· + is the Navier-Stokes operator. This computation is carried out in the Fourier domain [18]. For each iteration the dominating computation is the Fast Fourier Transform. Thus, the order of the algorithm is MNn log n where M is the number of iterations, N is the number of images to be registered, and n is the number of voxels in each image. The complexity increases only linearly as images are added, making the algorithm extremely scalable. Satisfactory correspondence is typically achieved after 200 iterations. In practice, we use a multi scale approach that initializes the fine (voxel) scale registration with the up-sampled correspondence computed at a coarser scale level. The finer scale levels only need to account for residue from coarser scale levels and thus require far fewer iterations to converge. 5 Driving application: Autism neuroimaging study In the following, we describe the application of the unbiased atlas construction to an ongoing infant autism study directed by Joseph Piven at UNC Chapel 14 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript Hill. From the partnership with our Psychiatry department, we have access to a morphologic MRI study with a large set of autistic children (N=50), de-velopmentally delayed subjects (N=25) and control subjects (N=25), scanned at age two with follow-up at age four. This project employs structural (MRI) and functional (fMRI) neuroimaging techniques to examine the neural basis of social cognitive, affective, and executive functioning deficits in autism. The combined use of neuroimaging and neuropsychological data should further our understanding of the causes of autism, thus assisting identification, prevention, and potentially treatment of this disorder. This study of early development imposes several significant challenges to im-age analysis. Cross-sectional comparison between healthy and autistic subjects at ages two and four and the study of brain growth over the two year follow-up period require adequate image analysis methods. Normative imaging data of healthy and patient populations in this age range are not available, suggesting the development of atlases that describe the mean and the statistical variabil-ity in early development. We plan to build population atlases that help us to study differences between groups and differences in growth patterns. Methods have to cope with the relatively large variability of head size and shape in these age groups and with significant longitudinal changes due to growth. 5.1 Gray-level MRI data for atlas-building The atlas-building described here uses high-resolution T1-weighted MRI with voxel size 0.781252×1.5mm3. The standard segmentation pipeline of the UNC autism image analysis group rigidly aligns the raw MR images to a Talairach coordinate space by specifying anterior and posterior commissure (AC-PC) 15 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript and the interhemispheric plane. The transformation also interpolates the im-ages to a standard isotropic voxel-size of 1×1×1mm3. We randomly selected 8 images from the total of currently over 80 cases in our database. The images were cropped for efficiency reasons. 5.2 Construction of unbiased gray-level atlas To evaluate the performance of this method we applied our algorithm to a set of intensity adjusted 3D MR brain images taken from eight different subjects. As a preprocessing step, these images were aligned using an affine transforma-tion. Axial slices from five of these eight initial images are shown in the first column of Figure 3. There is noticeable large deformation variation between these anatomies, especially around the ventricles. The second column of Fig. 3 shows the deformed (3D) images after 500 iterations of our algorithm. The de-formed images look very similar, as they have been deformed into the common coordinate space of the evolving template . Column 3 of Fig. 3 shows the ab-solute error between each input image and the initial template estimate. After applying our algorithm the deformed neuroanatomies are very close, resulting in a sharp final template estimate. Column 4 of Fig. 3 shows the absolute error between each deformed image and the final template estimate. This image is amplified to a maximum range to show the residual error. The squared error as a function of the number of iterations is displayed in column 5 of Fig. 3. The plots demonstrate the monotonous decrease of the error with increasing number of iterations. Figure 4 shows the initial and final estimate of the template. The initial tem-plate estimate is blurry since it is an average of the varying individual neu- 16 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript roanatomies. Ghosting is evident around the lateral ventricles and near the boundary of the brain. The final template estimate Fig. 4 right is close to each individual final results (column 2 Fig. 3). The sharpening of the resulting at-las in comparison to the initial estimate is also shown in Fig. 5. Whereas the initial template obtained by affine deformation and averaging does not show any details in cortical regions, these regions appear much sharper in the final template. 5.3 Annotated Atlas: Example Caudate Nucleus Atlases might also carry anatomical labels, which can be used to explain loca-tion and variability of anatomical structures if transformed to the atlas space. Also, they find use in the segmentation of new subjects by atlas deformation. Rohlfing et al. [25,26] deformed a set of labeled atlas images to a new image and showed that the combination of these independent classifiers was sub-stantially better than an individual classifier. Combination of the deformed label images was done using an extension of the STAPLE [28,35] algorithm. This classification method requires multiple high-dimensional deformations to be applied to each new individual, which poses a computational problem problem if applied to large clinical studies. The caudate study allows us to perform similar experiments. Each infant MRI is segmented into brain tissue, fluid, ventricles, and subcortical structures. Here, we will use the left and right caudate nucleus. In the following, we will describe the construction of a probabilistic caudate atlas which, similarly to Rohlfing's work, is built by deformation of a set of segmented images into a template. We will further show how this caudate atlas can be applied to 17 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript a b c d e Figure 3. Template Construction shown for five of the eight initial images. a) Set of initial images. b) Deformed images after 500 iterations. c) Absolute error between the initial images and the initial template estimate. Notice the structure of the residual error around the ventricles and the cortex. the deformed images and the final template estimate. These errors are normalized to the maximum intensity range to show the locations of the residual errors. Notice the lack of structure in the residue after the deformation. e) Sum of squared errors as a function of the number of iterations. segment new subject's images. 18 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript Figure 4. Estimated templates initially and after 500 iterations. The figure shows an axial slice through the ventricles and basal ganglia. Atlas construction is based on eight infant MRI. Figure 5. Estimated templates initially (left) and after 500 iterations (middle). The right graph illustrates a profile along the dashed line for the initial (line) and final template (dotted). 5.3.1 Caudate nucleus segmentations A specific structure of interest in autism research is the caudate nucleus, as this structure is associated with controlling motor function and indirectly with the repetitive behavior as observed in autistic children. On a first sight, the caudate seems easy to segment since the largest fraction of its boundary is bounded by the lateral ventricles and white matter. Portions of the boundary of the caudate can be localized with standard edge detection (provided the appropriate scales are chosen). However, the caudate is also adjacent to the 19 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript nucleus accumbens and the putamen where there are no visible boundaries in the MRI (see Fig. 6). The caudate and nucleus accumbens are distinguishable on histological slides, but not on MRI of this resolution. Another "trouble-spot" for the caudate is where it borders the putamen; there are "fingers" of cell bridges which span the gap between the two. We have developed a highly reliable manually assisted caudate segmentation using SNAP, a tool based on 3D geodesic snakes [15]. This tool reduced segmentation from about 2 hours to approximately 10 minutes, still including manual experts' definition of nucleus accumbens and putamen boundaries. Here, we used the caudate segmentations of the 8 cases selected for atlas build-ing to construct a probabilistic caudate atlas. Further, we selected 5 new cases from our database which are not part of these 8 images for testing caudate seg-mentation by atlas deformation. These 5 cases were taken from the reliability test series, with the advantage that we have 6 manual segmentations (2 raters with 3 segmentations each) for each case, which represent a probabilistic gold standard. 5.3.2 Probabilistic caudate template using voxel voting The caudate segmentations of the eight atlas images were transformed into the atlas space by applying the individual deformations obtained during atlas building. The set of deformed caudate segmentations can be used to build a probabilistic caudate atlas. The segmentations were combined by a voxel-wise voting scheme counting the number of segmentations at each voxel. Normal-ization by the number of images finally results in a probabilistic caudate tem-plate. This segmentation template can be seen as a level set within the range 20 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript Figure 6. Two and three-dimensional views of the caudate nucleus. Coronal slice of the caudate: original T1-weighted MRI (left), and overlay of segmented structures (middle). Right and left caudate are shown shaded in green and red; left and right putamen are sketched in yellow, laterally exterior to the caudates. The nucleus accumbens is sketched in red outline. Note the lack of contrast at the boundary between the caudate and the nucleus accumbens, and the fine-scale cell bridges between the caudate and the putamen. At right is a 3D view of the caudate and putamen relative to the lateral ventricles. [0 · · · 1] where the level 0.5 defines the average shape. Figure 7 left illustrates a sagittal cut and the corresponding 3D surface of the average shape. 5.3.3 Probabilistic caudate template using STAPLE Warfield and Zou et al. [28,35] developed an algorithm to calculate the compos-ite gold standard estimate from multiple manual segmentations. The Simultaneous Truth and Performance Level Estimation (STAPLE) method is based on an expectation maximization (EM) algorithm. Given a set of binary segmenta-tions of the same object, STAPLE calculates the maximum likelihood estimate of the composite "gold standard" or the best estimate of the unknown gold standard. The algorithm calculates the specificity and sensitivity of each seg-mentation in an iterative way. The major difference over voxel-wise voting is its ability to assign weights for each individual segmentation proportional to the performance, i.e. segmentations which are closer to the estimated gold 21 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript standard get larger weights. The STAPLE algorithm is applied to the set of eight segmented left and right caudates after they were deformed to the atlas. It is important to notice that we limited the STAPLE calculations to disputed voxels only, so that regions with only background and agreement within re-gions were not taken into account. Figure 7 right illustrates a sagittal cut and the corresponding 3D average object. a b Figure 7. Probabilistic caudate atlas. The manual segmentations of the eight images used for atlas building are deformed using the individual deformation fields obtained during atlas building. a) The eight deformed segmentations are superimposed by voxel-voting and normalized to form a caudate probability map. b) The STAPLE algorithm is applied to represent a probabilistic best estimate (gold standard) of the true structure. The images show sagittal slice through this probabilistic images and 3D surfaces of the average structure (probability 0.5). 5.3.4 Segmentation of new unknown subjects The combined unbiased MRI and caudate atlases can be used to segment new subjects by atlas deformation. The use of the unbiased atlas constructed from a representative set of images eliminates the need to deform and combine multiple labeled atlases. The MRI atlas is deformed using fluid deformation. The same transformation is then applied to the caudate template to trans-fer the probabilistic caudate to the new image. Figure 8 illustrates manual segmentation and segmentation by deformation-segmentation. The compar- 22 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript ison demonstrates the potential of this segmentation technique not only to segment well-visible boundaries but also transition regions (like nucleus ac-cumbens, e.g.) which can only be segmented in the anatomical context of embedding neighboring structures. Manual Segmentation Segmentation by Atlas Deformation Figure 8. Coronal and sagittal view of the right caudate segmented manually (left two images) and by nonlinear deformation of the atlas (right two images). Atlas deformation allows to capture the inferior and lateral boundaries of the caudate with nucleus accumbens and putamen although there are no visible contours in the MRI image (see also Fig. 6 for anatomical reference). 5.3.5 Validation The deformation-segmentation is validated against the gold standard of hu-man expert segmentation. As discussed earlier, each of the 5 new subjects selected for segmentation also comes with a set of 6 expert segmentations (3 repeated segmentations by 2 raters). This allows to compare not only binary segmentations but also probabilistic segmentations. We use a previously developed validation package VALMET [11] that includes a probabilistic overlap measure between two fuzzy segmentations. This metric is derived from the normalized L1 distance between two probability distribu-tions 23 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript POV (A,B)=1 − R |PA − PB| 2 R PAB . (9) PA and PB are the probability distributions representing the two fuzzy seg-mentations and PAB is the joint probability distribution. In this study, PA and PB are calculated by integrating a set of binary segmentations with sub-sequent normalization to the range of [0 · · · 1], whereas PAB is calculated by integrating the set all binary segmentations and appropriate normalization. The numerator expresses the probabilities of non-intersecting regions. Table 1 lists the left and right caudate volumes for manual segmentation (user assisted geodesic snake) and deformation segmentation. The results are en-couraging but also show the limitations of high dimensional deformation with-out using landmarks. There is one case (5007 right) with very large volumetric difference, probably due to very thin ventricles creating local large scale de-formation. Table 5.3.5 lists overlap measures for pairs of binary segmentations (top) and pairs of probabilistic segmentations (bottom). Binary objects are extracted from the probabilistic segmentations by choosing level 4 close to the middle level. The overlap ratio is defined as the intersection divided by the average. The probabilistic overlap uses the probabilistic caudate atlas constructed from 8 cases and the manual experts segmentations (6 cases). Overlap results are in the range 0.85 to 0.90, which is encouraging given the small size of the objects. Manual raters still have a significantly better intra- and inter-rater reliability, however this only comes after several months of training with several reliability studies. 24 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript Case Manual Vol. Atlas Deformed Vol. Diff. Manual vs. Def. (L) (R) (L) (R) (L) (R) Case 5003 3583.83 3596.00 3852.13 3638.25 7% 1% Case 5004 4054.67 3910.83 3833.88 3881.13 5% 1% Case 5007 4045.17 4319.67 3672.25 3562.25 9% 18% Case 5011 3933.67 3899.67 4151.13 4089.25 6% 5% Case 5020 3997.67 4112.17 3833.88 3881.13 4% 6% Average 3922.99 3967.67 3868.65 3810.40 6.2% 6.2% Std Dev 195.57 269.62 173.97 211.53 Table 1 Comparison between manual segmentation and automatic segmentation by atlas deformation shown for five cases. The table illustrate the volumes for manual seg-mentation and for segmentation by atlas deformation (volumes in mm3). The last two columns list the absolute percentage differences. Volumes are represented by the average objects at probability level 0.5. 6 Conclusions In this paper a new concept for unbiased construction of atlases is presented based on Frechet means in metric spaces. This approach results in an itera-tive algorithm of simultaneous deformation of a population of subject images into a new average image that evolves iteratively. This technique avoids the systematic bias introduced by selecting a template but also the combinatorial problem of deformation of a large number of datasets into each new subject. 25 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript The new techniques produces a population average image which might serve as a template to represent the population group. Sharpness of structures indi-cates the quality of match and residual biological variability. Local variability of brain structures is encoded in the set of deformation maps. We plan to explore this information in our future work. Results demonstrate the application of the new technique to eight 3D MRI of children at age 2-years. A visual comparison of the resulting average atlas with each individual image suggests that the atlas represents the average while still being sharp. As each individual deformation is diffeomorphic, we can apply transformations in both directions from the individual into the atlas and back. We can also transform images into each other by cascading their transformation and inverse transformations. The caudate segmentation study with a set of over 80 segmented images and a reliability study of 5 subjects with sets of repeated segmentations by several experts form an excellent database to test and validate intermediate stages of our development. Moreover, the caudate is an excellent example of an anatom-ical structure that is not fully delineated by strong contrast boundaries but can only be segmented in the context of embedding structures or a geometric model. This supports use of deformable atlas registration where constraints are provided through the use of a volumetric, unbiased atlas. The caudate segmentation experiment clearly demonstrates the accuracy to be obtained by deformation segmentation without landmarking. It can be seen from the results that in one of the validation cases (subject 5007) lume of the right caudate was substantially underestimated. This was primarily as a result of a local extrema in the greedy optimization strategy used. We 26 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript are in the process of augmenting the current completely automated process with manual landmarking and a complete space-time optimization which will greatly improve the accuracy of the segmentations. Acknowledgments This research is supported by the NIH NIBIB grant P01 EB002779, the NIH Conte Center MH064065, DOD Prostate Cancer Research Program DAMD17- 03-1- 0134,and the UNC Neurodevelopmental Research Core NDRC, subcore Neuroimaging.. The MRI images of infants, caudate images and expert man-ual segmentations are funded by NIH RO1 MH61696 and NIMH MH 64580 (PI: Joseph Piven). Manual segmentations are by Michael Graves and Rachel Gimpel, with protocol development in collaboration with Cody Hazlett. We would like to especially thank Mark Foskey and Peter Lorenzen for help with the preparation of the manuscript. References [1] Magnotta V A, Bockholt H J, Johnson H F, Christensen G E, and Andreasen N C. Subcortical, cerebellar, and magnetic resonance based consistent brain image registration. NeuroImage, 19:233-245, June 2003. [2] B. Avants and J.C. Gee. Symmetric geodesic shape averaging and shape interpolation. 2004. in press. [3] K.K. Bhatia, J.V. Hajnal, B.K. Puri, A.D. Edwards, and D. Rueckert. Consistent groupwise non-rigid registration for atlas construction. In IEEE International Symposium on Biomedical Imaging. IEEE, 2004. 27 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript [4] F.L. Bookstein. Morphometric Tools for Landmark Data. Cambridge University Press, New York, 1991. [5] John C. Csernansky, Sarang Joshi, Lei Wang, Mokhtar Gado, J. Philip Miller, Ulf Grenander, and Michael I. Miller. Hippocampal morphometry in schizophrenia by high dimensional brain mapping. Proceedings of the National Academy of Science, 95:11406-11411, September 1998. [6] P. Dupuis, U. Grenander, and M.I. Miller. Variational problems on flows of diffeomorphisms for image matching. Quarterly of Applied Mathematics, 1997. [7] P. Thomas Fletcher, Conglin Lu, and Sarang Joshi. Statistics of shape via principal geodesic analysis on lie groups. In CVPR2003. CVPR, 2003. [8] PT Fletcher, S Joshi, C Lu, and S Pizer. Gaussian distributions on lie groups and their application to statistical shape analysis. In C Taylor and JA Noble, editors, Information Processing in Medical Imaging (IPMI), volume 2732, pages 450-462. Lecture Notes in Computer Science, Springer, July 2003. [9] Maurice Frechet. Les elements aleatoires de nature quelconque dans un espace distancie. Annales De L'Institut Henri Poincare, 10:215-310, 1948. [10] J. C. Gee, M. Reivich, and R. Bajcsy. Elastically deforming an atlas to match anatomical brain images. J. COmput. Assis. Tomogr., 17:225-236, 1993. [11] G. Gerig, M. Jomier, and M. Chakos. VALMET: a new validation tool for assessing and improving 3D object segmentation. In W. Niessen and M. Viergever, editors, Medical Image Computing and Computer-Assisted Intervention MICCAI, volume 2208, pages 516-523. Springer-Verlag, October 2001. [12] U. Grenander. General Pattern Theory. Oxford Univ. Press, 1994. [13] U. Grenander and M. I. Miller. Computational anatomy: An emerging discipline. Quarterly of Applied Mathematics, 56:617-694, 1998. 28 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript [14] A. Guimond, J. Meunier, and J.-P. Thirion. Average brain models: A convergence study. Computer Vision and Image Understanding, 77(2):192-210, 2000. [15] Jianchun He and Gary E. Christensen. Large deformation inverse consistent elastic image registration. In C. J. Taylor and J. A. Noble, editors, IPMI3003, LNCS 2732, pages 438-449. IPMI, Springer-Verlag, 2003. [16] S. Ho, E. Bullitt, and G. Gerig. Level set evolution with region competition: Automatic 3-D segmentation of brain tumors. In R. Katsuri, D. Laurendeau, and C. Suen, editors, Proc. 16th International Conference on Pattern Recognition, pages 532-535. IEEE Computer Society, August 2002. [17] K. H. Hohne, M. Bomans, M. Riemer, U. Tiede R. Shubert, and W. Lierse. A 3d anatomical atlas based on a volume model. IEEE Comp. Grpahics. Appl., pages 72-78, December 1992. [18] S. Joshi, U. Grenander, and M.I. Miller. On the geometry and shape of brain sub-manifolds. International Journal of Pattern Recognition and Artificial Intelligence: Special Issue on Processing of MR Images of the Human, 11(8):1317-1343, 1997. [19] S. Joshi, P. Lorenzen, G. Gerig, and E. Bullitt. Structural and radiometric asymmetry in brain images. Medical Image Analysis, 7(2):155-170, June 2003. [20] Sarang C. Joshi and Michael I. Miller. Landmark matching via large deformation diffeomorphisms. IEEE Transactions On Image Processing, 9(8):1357-1370, August 2000. [21] D.G. Kendall. Shape manifolds, procrustean metrics and complex projective spaces. Bulletin London Mathematical Society, 16:81-121, 1984. [22] Peter Lorenzen and Sarang Joshi. High-dimensional mult-modal image 29 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript registration. Workshop on Biomedical Image Registration (WBIR), LNCS- 2717:234-243, 2003. [23] M. I. Miller, A. Trouve, and L. Younes. On the metrics and euler-lagrange equations of computational anatomy. Annual Review of Biomedical Engineering, 4:375-405, 2002. [24] M.I. Miller and L. Younes. Group actions, homeomorphisms, and matching: A general framework. International Journal of Computer Vision, 41:61-84, January 2001. [25] Michael I. Miller, Sarang C. Joshi, and Gary E. Christensen. Large deformation fluid diffeomorphisms for landmark and image matching. In Arthur W. Toga, editor, Brain Warping, chapter 7. Academic Press, 1999. [26] Torstens Rohlfing, Daniel B. Russakoff, and Calvin.R. Maurer. Expectation maximization strategies for multi-atlas multi-label segmentation. In Ch. Taylor and A. Noble, editors, Proceedings of Information Processing in Medical Imaging IPMI, volume 2732 of Lecture Notes in Computer Science LNCS, pages 210-221. Springer-Verlag Berlin, 2003. [27] Torstens Rohlfing, Daniel B. Russakoff, and Calvin.R. Maurer. Extraction and application of expert priors to combine multiple segmentations of human brain tissue. In R.E. Ellis and T.M. Peters, editors, Proceedings Medical Image Computing and Computer Assisted Intervention MICCAI, volume 2879 of Lecture Notes in Computer Science LNCS, pages 578-585. Springer-Verlag Berlin, 2003. [28] Dinggang Shen and Christos Davatzikos. Hammer: Hierarchical attribute matching mechanism for elastic registration. IEE Transactions on Medical Imaging, 21(11):1421-1439, November 2002. [29] Warfield S.K., K.H. Zou, andW.M. III.Wells. Validation of image segmentation 30 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript and expert quality with an expectation-maximization algorithm. In T. Dohi and R. Kikinis, editors, Medical Image Computing and Computer Assisted Interventions, number 2488 in Lecture Notes in Computer Science LNCS, pages 298-306. Springer, Sept 2002. [30] J. Talairach and P. Tournoux. Co-planar Stereotaxis Atlas of the Human Brain. Thieme Medical Publishers, 1988. [31] P. M. Thompson, J. Moussai, S. Zohoori, A Goldkorn, A. A. Khan, M. S. Mega, G. W. Small, J. L. Cummings, and A. W. Toga. Cortical variability and asymmetry in normal aging and alzheimer's disease. Cerebral Cortex, 8(6):492- 509, September 1998. [32] Paul M. Thompson and Arthur W. Toga. Detection, visualization and animation of abnormal anatomic structure with a deformable probabilistic brain atlas based on random vector field transformations. Medical Image Analysis, 1:271-294, 1997. [33] Paul M. Thompson and Arthur W. Toga. A framework for computational anatomy. Computing and Visualization in Science, (5):13-34, 2002. [34] A. W. Toga. Brain Warping. Academic Press, 1999. [35] R. P. Woods, , S. T. Grafton, J. D. G. Watson, N. L. Sicotte, A. W. Toga, and J. C. Mazziotta. Automated image registration:ii intersubject validation of linear and nonlinear models. Journal of Computer Assited Tomogrophy, 22:155- 165, August 1998. [36] K.H. Zou, S.K. Warfield, J.R. Fielding, C.M.C. Tempany, W.M. Wells, M.R. Kaus, F.A. Jolesh, and R. Kikinis. Statistical validation based on parametric receiver operating characteristic analysis of continuous classification data. Academic Radiology, 10(12):1359-1368, Dec 2003. 31 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript Overlap Ratio Case Mean Level (L) Mean Level (R) STAPLE (L) STAPLE (R) Case 5003 0.88 0.86 0.88 0.87 Case 5004 0.91 0.92 0.90 0.91 Case 5007 0.88 0.86 0.88 0.87 Case 5011 0.88 0.85 0.88 0.85 Case 5020 0.90 0.86 0.90 0.86 Average 0.89 0.87 0.89 0.87 Probabilistic Overlap Case Fuzzy (L) Fuzzy (R) STAPLE (L) STAPLE (R) Case 5003 0.84 0.82 0.89 0.86 Case 5004 0.86 0.87 0.90 0.90 Case 5007 0.85 0.84 0.88 0.85 Case 5011 0.84 0.82 0.88 0.85 Case 5020 0.87 0.84 0.90 0.85 Average 0.85 0.84 0.89 0.86 Table 2 Comparison between manual segmentation and automatic segmentation by atlas deformation shown for five cases. The tables illustrate the differences for the prob-abilistic caudate atlas templates represented as a probability atlas (fuzzy caudate atlas) and derived by applying the STAPLE algorithm. Top: Overlap ratio between pairs of binary objects derived as the mean level of the probability atlas templates. Bottom: Probabilistic overlap calculated between pairs of probabilistic atlas tem-plates. 32 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript