|Home | About | Journals | Submit | Contact Us | Français|
In this paper, we describe ongoing work in the Image Processing and Analysis Group (IPAG) at Yale University specifically aimed at the analysis of structural information as represented within magnetic resonance images (MRI) of the human brain. Specifically, we will describe our applied mathematical approaches to the segmentation of cortical and subcortical structure, the analysis of white matter fiber tracks using diffusion tensor imaging (DTI), and the intersubject registration of neuroanatomical (aMRI) data sets. Many of our methods rally around the use of geometric constraints, statistical (MAP) estimation, and the use of level set evolution strategies. The analysis of gray matter structure and connecting white matter paths combined with the ability to bring all information into a common space via intersubject registration should provide us with a rich set of data to investigate structure and variation in the human brain in neuropsychiatric disorders, as well as provide a basis for current work in the development of integrated brain function–structure analysis.
Accurate and robust extraction and measurement of neuroanatomic structure in the human brain from magnetic resonance images (MRI) remain a challenging area of research. Within the Image Processing and Analysis Group (IPAG) at Yale University, we have been developing mathematical approaches to aspects of this problem, with focus on using appropriate geometrical and statistical constraints and decision-making strategies for particular subproblems in the domain. We have identified four key areas important to our efforts, each of which requires a unique approach: (1) the segmentation and measurement of cortical gray matter structure, (2) the segmentation and measurement of subcortical gray matter structure, (3) tracking and analysis of white matter fiber pathways, and (4) structurally focused intersubject image registration for developing multisubject measures. There are, of course, other processing issues in the analysis of human brain images, such as bias field correction, that are not discussed here. Methods for the analysis of subcortical gray matter, cortical gray matter, and white matter combined with intersubject registration techniques provide the foundation for detailed investigation of brain structure in neuropsychiatric disorders. These structural methods also provide the basis for our ongoing research of integrated brain function–structure analysis techniques. In this section, we present related work in each area and describe our own efforts in the ensuing sections.
Work in segmenting the cortex from three-dimensional MR images has fallen into two broad categories: voxel classification and deformable models. Classifying gray and white matter by voxel intensity can incorporate voxel continuity or homogeneity using, for example, Markov random fields (Geman and Geman, 1984; Leahy et al., 1991) to model probabilistic constraints on the image or fuzzy logic (Barra and Boire, 2001). The approach of Wells et al. (1994) estimates tissue classes (gray matter, white matter, cerebrospinal fluid (CSF)) while simultaneously estimating the bias field using an expectation-maximization (EM) strategy. Cline et al. (1990) use multispectral voxel classification in conjunction with connectivity to segment the brain into tissue types. Material mixture models (Liang et al., 1992) have also been used. Region-based methods of this type typically require further processing to group segmented regions into coherent structures.
Snakes or active contour models (ACMs) (Kass et al., 1988) are energy minimizing parametric contours with smoothness constraints. Unlike level set implementations (Osher and Paragios, 2003; Osher and Sethian, 1988), the direct implementation of this energy model is not capable of handling topological changes of the evolving contour without explicit discrete model parameter manipulations. Deformable models using level set methods (Malladi et al., 1995, 1996; Sethian, 1996) analyze and compute interface motion representing the boundary as a propagating wavefront. Level set methods are quite powerful and have useful properties of topologic flexibility. They need to be carefully constrained to avoid stopping short or flowing too freely beyond the desired boundaries. Methods have been devised to preserve topology when necessary (Han et al., 2003). See below for further discussion of level set methods.
Deformable models specifically tailored for segmenting the cortical gray matter have been devised. Davatzikos and Prince (1995) proposed an active contour algorithm for determining the spine of a ribbon modeling the outer cortex. Xu et al. (1998) use gradient vector flow fields in conjunction with tissue membership to control the deformation of active contour models for finding the central gray matter layer. Teo et al. (1997) designed an approach in which white matter and CSF regions were first segmented, followed by connectivity analysis and a constrained deformation from the white matter boundary. MacDonald et al. (1998) simultaneously deform multiple surfaces to segment the brain using multiple cost function constraint terms. This approach takes advantage of the information of the interrelation between the surfaces of interest, but is computationally expensive with many parameters to tune. Kapur et al. (1996) also use a snake approach in conjunction with EM segmentation and mathematical morphology. Unlike subcortical structure (discussed below), the cortex is highly variable and thus generally not well suited for specific prior shape models (Staib et al., 2000) unless limited to specific regions of the cortex; more generic constraints are necessary. In work we describe in more detail in the next section Zeng et al. (1999a,b) developed a coupled level set algorithm to segment the cortex using a geometric constraint based on the consistency of cortical thickness.
Once the cortical layers are segmented, a variety of useful neuroanatomical measurement parameters can be derived. Examples of these include sulcal depth, cortical thickness, or cortical shape (e.g., Vaillant and Davatzikos, 1997; Zeng et al., 1998). In the work of Vaillant and Davatzikos (1997), special “ribbon” operators were developed to characterize the depth of the sulcal grooves. A related example is the automated extraction of sulcal and gyral curves using crest lines from the differential geometric properties of a segmented brain surface (Thirion, 1996). The extraction of shape features from a segmented anatomical object in this way may be useful not only for measurement but also feature-based registration as well.
Statistical models can be powerful tools to directly capture the variability of structures being modeled. Such techniques are a necessity for the segmentation of subcortical structure that has consistent shape but is poorly defined by image features. Atlas registration for the purposes of segmentation (Collins et al., 1995; Declerck et al., 1995) is one way of using prior shape information. Collins et al. (1992), for example, segment the brain using an elastic registration to an average brain, based on a hierarchical local correlation. The average brain provides strong prior information about the expected image data and can be used to form probabilistic brain atlases (Collins et al., 1992; Thompson et al., 1997). Specific models for prior shape have been used successfully in our lab (Staib and Duncan, 1992, 1996; Wang and Staib, 1998; Yang et al., 2003, 2004a) and by other groups (Cootes et al., 1993; Szekély et al., 1995) for segmentation. The statistics of a sample of images can be used to guide the deformation in a way governed by the measured variation of individuals. Region-based information can also be combined with prior models (Chakraborty and Duncan, 1999; Chakraborty et al., 1996; Staib et al., 1997) to enhance robustness.
Level set methods and new energy terms have been reported to increase the capture range of deformable models and incorporate prior shape information. Chan and Vese (2001) proposed a level set method that can detect objects whose boundaries are not necessarily defined by gray level gradients. Leventon et al. (2000) extended Caselles et al.'s (1997) geodesic active contours by incorporating shape information into the evolution process. In our work, described in this paper, we adopt a level set approach using prior information of the shape of an object and its neighbors (Yang et al., 2003, 2004a). We note that Tsai et al. (2001, 2003) describe a similar approach that uses a global multishape model. We feel that our approach is more flexible in the modeling of joint priors and accommodates situations with limited interobject information and variation in contrast and discernibility among the objects.
Typically, the accurate and meaningful measurement of normal and abnormal subcortical structure includes volume, surface area, as well as shape measures. These analysis methods require processing beyond the initial steps of segmentation or registration of data sets (Gerig et al., 2001).
While the basic anatomy of white matter tracts in the human brain is generally known from anatomical dissection, much is unknown about its interconnections and its natural variations. The characterization and quantitative measurement of its connections is of fundamental importance in understanding brain function. Diffusion tensor magnetic resonance imaging (DT-MRI) has emerged as a noninvasive imaging modality capable of providing this information in vivo, enabling the detailed study of white matter structure in the human brain.
Brain white matter, because of the long and fibrous nature of axons, exhibits higher restriction to water diffusion across the fibers than along them. This directional variation is measured in diffusivity rates and can be captured by diffusion-weighted MRI. By acquiring diffusion-weighted data in at least six non-collinear directions, it is possible to estimate a 3 × 3 symmetric matrix (i.e., diffusion tensor) that characterizes diffusion in anisotropic systems (Basser and Pierpaoli, 1996). After tensor diagonalization, the eigenvector corresponding to the largest eigenvalue is considered to point along the direction of a fiber bundle.
Many connectivity studies relying on the straightforward integration of the principal tensor eigenvector have been described in the literature (Mori and van Zijl, 2002). Their reliability, however, is limited by acquisition noise and partial voluming due to fiber tracts that cross, branch, and merge. To account for these variations, level set methods (Osher and Fedkiw, 2003) have been employed (e.g., Lenglet et al., 2003; O'Donnell et al., 2002; Parker et al., 2002). These techniques model the evolution of an advancing front through the white matter tracts by following the local directionality provided by the diffusion tensor field. Such methods have been shown to be more robust to noise and singularities than classical streamlining methods. A tractography technique based on the fast marching method (FMM) was used by Parker et al. (2002). A front was evolved with a speed proportional to the collinearity between the front normal and the principal tensor eigenvector. A discrete approximation of front direction had to be used to drive the evolution through the eigenvector field, because the original FMM cannot handle propagation in oriented domains. Others have posed the connectivity problem in a Riemannian framework (e.g., Lenglet et al., 2003; O'Donnell et al., 2002).
As will be described later, in our group, we have developed a novel wavefront propagation method for estimating the connectivity in the white matter of the brain using DT-MRI.
To gain statistical power for comparing groups of subjects in either structural or functional brain image analysis, it can be useful to pool image information from many subjects. In these cases, a key step is the formation of multisubject composite activation maps where each subject is nonrigidly mapped to a common coordinate space. Group characterization using composite activation or structural variation maps can then be computed on a voxel-wise basis. However, nonrigid registration of brain images is a difficult task.
Previously, Talairach-type methods (Talairach and Tournoux, 1988) have been used but the limitations of these methods are well known (Mazziotta et al., 1995; Toga and Thompson, 1999). True nonrigid registration is necessary to bring the subjects into a common space. There have been many approaches recently to nonrigid registration, with a particular emphasis on applications to brain imaging (see the recent special journal issues Goshtasby et al., 2003; Pluim and Fitzpatrick, 2003). Most commonly, nonrigid registration methods use image intensities to compute the transformation (Christensen et al., 1996; Friston et al., 1995; Rueckert et al., 1999). The high anatomic variability of the cortex can often result in intensity based methods yielding inaccurate results due to local minima, as was recently demonstrated (Hellier et al., 2003). These errors are of particular concern in functional MRI (fMRI) analysis for evaluating cortical activations.
Feature-based methods have been developed to overcome the limitations of gray-level methods (Collins et al., 1998; Corouge et al., 2001; Davatzikos, 1997). None of these methods, however, is able to handle large variations in sulcal anatomy, as well as irregular sulcal branching and discontinuity. In our work, which will be described in a later section, we extend the robust point matching (RPM) algorithm (Chui et al., 2003) to better address brain registration. We use a point matching strategy to register structurally salient regions as part of our strategy to generate higher resolution composite functional maps tailored to specific regions.
We have developed a coupled surfaces approach (Zeng et al., 1999a,b) to segment cortical structure. In this work, a local gray level operator designed to obtain the likelihood of each voxel lying on the outer and inner cortical surfaces, respectively, is used instead of simple gradient information. A gradient operator will respond to all gray level transitions at a given scale. By focusing on the gray-level transitions of interest, we remove extraneous information that could impede the propagation of the level sets. Key to the approach is that the surfaces representing the gray matter–white matter cortical interface (inner) and gray–CSF interface (outer) are found in an integrated manner by evolving two coupled level sets. Thus, starting from inside the inner bounding surface (gray–white boundary), with an offset in between, the interfaces propagate along the normal direction stopping at the desired location, while maintaining a distance between them.
Embedding each surface as the zero level set in its own level function ψ, we have two equations:
where Fin and Fout are speed functions of the surface normal, image-derived information and distance between the two surfaces. The coupling is embedded in the design of Fin and Fout. Where the distance between the two surfaces is within the normal range for cortical thickness, the two surfaces propagate according to the image-based information; where the distance is out of the normal range, the distance constrains the propagation. The level set implementation provides an easy and natural way to evaluate the distance between the two surfaces because the value of the level function at any point is simply the distance from this point to the current front, which as in Sethian (1996) is calculated as the shortest distance from this point to the points on the front. In our case of two moving surfaces, for any point on the inner moving surface, the distance to the outer moving surface is the value ψout at this point and vice versa for the point on the outer moving surface, that is, its distance is ψin.
This approach has been tested on simulated data (with known underlying segmentations) created from the Montreal Neurological Institute's MR simulator (Brainweb) (using different noise conditions and 1 mm3 voxels) with very encouraging results: true positive rates greater than 92% and false positive rates less than 6% for segmentation of cortical gray matter (Zeng et al., 1999a,b). In addition, as reported in part in Zeng et al. (1999a,b), we have run our coupled-surface cortical layer segmentation strategy on 20 normal brains and compared our results to known manual segmentations from the Internet Brain Segmentation Repository (IBSR) (about 1 × 1 × 3 mm3 voxels). In this work, we were able to illustrate that our approach outperformed six other automated algorithms in the literature in terms of the average overlap with pooled-manual-expert segmentations (with 1.0 perfect and 0.0 the worst, our approach on cortical gray matter achieved 0.701, and the next best algorithm, 0.564). In addition, we have run this approach on over 60 subjects acquired on our 1.5-T GE Signa system at Yale (SPGR acquisition, (1.2 mm)3 voxels). A typical result from these data is shown in Figs. Figs.11 and and2.2. For one group of 14 human subjects, there was 87% true positive findings for each subject in comparison to manual tracings of the frontal cortical gray matter layer, although with false positives in the 20% range. Errors can occur where the outer cortical surface is obscured, such as due to susceptibility artifacts in the orbital frontal cortex.
We have developed a novel method for the segmentation of multiple objects from three-dimensional medical images using interobject constraints that we have employed as our subcortical structure segmentation strategy (Yang and Duncan, 2003; Yang et al., 2002, 2003, 2004a). Our method is motivated by the observation that neighboring structures have consistent locations and shapes that provide configurations and context that aid in segmentation. In this effort, we have defined a maximum a posteriori (MAP) estimation framework using the constraining information provided by neighboring objects to segment several objects simultaneously. Within the MAP strategy, we assume that the likelihood (data adherence) term and the prior term are Gaussian. The segmented surface can be found using discrete approximations to level set functions and computing the associated Euler–Lagrange equations. The contours evolve both according to prior information related to the shape of the object of interest and relational shape and position information, as well as the image gray level information. We have recently compared (Yang et al., 2004b) these level-set-based shape models with point-based models and have been able to show that the priors formed from either method are not statistically different for the test shapes examined, although the level set approach can accommodate varying topology, a notion that we intend to ultimately exploit in our work. These methods are useful in situations where there is limited interobject information as opposed to robust global atlases.
Relating the development from Yang et al. (2002), but recasting the notation slightly, we describe our approach as follows. Consider a structural image represented as a random field I that has i = 1,…,M structures of interest. Next assume that the boundary of each object i can be embedded as the zero level set of a three-dimensional (Euclidean) distance function ψi. These distance functions can be sampled yielding discrete three-dimensional versions for each structure of interest Si, i = 1,…,M. Furthermore, we can model the distributions of shape, size, and position of any structure i by assuming Si to be a random field (rf) and then measuring instances of each rf over several subjects, thereby constructing the probability density functions (pdfs) p(Si).
Each outcome of Si is represented by an N3-long column vector that is a row-by-row stacking of the matrix found by discretely sampling the corresponding level set distance function ψi in the image space. We have shown (Yang and Duncan, 2003) that the level set distance function of any one object Si that depends on its neighbors can be estimated using the following MAP framework:
for i = 1, 2,…,M and where p(I|S1, S2,…,SM) is the probability of producing an image I given S1, S2,…,SM. In three-dimensional, assuming gray level homogeneity within each object, we implemented the following image-data term (adapted from Chan and Vese, 2001):
where c1i and r1i are the average and variance of I inside the zero level set of Si. c2i and σ2i are the average and variance of I outside the zero level set of Si but within a domain Ωi that contains Si. In our work, we typically set the domain Ωi to be where the level set distance is no more than the average diameter of the object of interest.
Furthermore, p(S1,S2,...,SM) is the joint density function of all the M objects. It contains neighbor prior information such as the relative position and shape among the objects. A variety of relational assumptions can be used here, but in situations where the neighbors are sometimes difficult to locate, we initially assumed that each object is related to the key object (denoted object k) independently. In this case, the joint density function can be simplified to:
where Δjk = p(Sj|Sk)= Sj – Sk is the difference between the level sets of object j and k. The process of defining the joint density function p(S1,S2,…,SM) is simplified to building only the self prior p(Sk) and the local neighbor priors p(Δjk), j,k = 1,2,…,M; j ≠ k.
Consider a training set of n aligned images, with M objects or structures in each image. As described above, each object i in the training set is embedded as the zero level set of a higher dimensional level set ψi whose discretized distance function is the N3 column vector Si. As in our previous efforts in the estimation of the segmentation of a single object in an image (e.g., Chakraborty and Duncan, 1999; Staib and Duncan, 1992), we deem it important to first be able to model the range of plausible object self-shape information. Here, we assume that the shapes vary smoothly in a relatively compact portion of a high dimensional manifold such that we can model their variation using principal component analysis (PCA). Thus, the pdf of the level function of object i can be computed using PCA similar to what is done for point distribution models (Cootes et al., 1993). An estimate of Si can be represented by the mean level set i, p principal components Ui, and a p dimensional vector of coefficients (where p < n), αi: i = Uiαi + i. Under the assumption of a Gaussian distribution of shape represented by αi, we can compute the probability of object i: p(αi)= (0,Σ).
The level set representation of shape provides tolerance to slight misalignment of object shape in an attempt to avoid having to solve the general correspondence problem. In practice, the variations captured by the principal components in the level set distribution model (Ui) in this paper are based on a rigid alignment of the training data and may contain undesired residuals due to misalignment. We are looking to improve the alignment method to reduce such residuals and undesired topology changes.
Due to the assumptions in our initial work above, that is, that each neighbor is independently related to the object of interest, we can use the difference between the two level sets Sj – Sk as the representation of the neighbor difference Δjk, j = 1, 2, 3,…,M shown above where this implicitly represents P(Sj|Sk. We again assume that the distribution of Δjk's for any key object k form a compact portion of a high dimensional space such that the distribution can be parameterized using a linear PCA formulation. Thus, the range of neighbor-to-object variation for each object can be found from the mean neighbor difference jk, and p principal components Pjk and a p dimensional vector of coefficients, βjk: jk = Pjkβjk + jk. We assume the neighbor difference Δjk to be Gaussian distributed over βjk: p(βjk) = (0,Λjk).
Using the techniques described in Chan and Vese (2001),we compute the associated Euler–Lagrange equation for each unknown level set function (written here for convenience in terms of the continuous functions ψk):
To simplify the complexity of the segmentation system, we generally choose the parameters in our experiments as follows: Ωik = Ωkk = Ωk, λ1k = λ2k = λk = 1 − Ωk, μk = 0.00005 × 2552, νk = 0 (Chan and Vese, 2001). This leaves us only one free parameter (ωk) to balance the influence of two terms: the image data term and the neighbor prior term for each object. G denotes the conversion from a matrix to a vector by column scanning. g is the inverse of G. The trade-off between neighbor prior and image information depends on how much faith one has in the neighbor prior model and the image-derived information for a given application. We set these parameters empirically given the image quality and the neighbor prior information.
In extensive testing using simulated image data, we have found that using the shape prior alone reduced the final error in the presence of varying noise and surface seed initialization, but the addition of the neighbor prior always reduced the segmentation errors limiting them to a very small range even for large noise variance and for varying seed points (Yang et al., 2003, 2004a).
To illustrate the utility of our neighbor-constrained approach to segment a particular neuroanatomical structure of interest to us, we have applied it to the segmentation of the left and right amygdalae and the left and right hippocampi, structures of particular interest in autism. Note that here, contralateral structures provide context due to the strong bilateral symmetry in the brain. In this work, we compared the results of our fully automated, neighbor-constrained segmentations to slice-by-slice manual tracings of the same structure using a set of N = 12 rigidly prealigned normal male subjects with an age range of 14–43. Note that even a relatively small sample provides a valuable constraint for images of subjects in the same population group. These three-dimensional anatomical MR images have (1.2 mm)3 resolution. Prior distributions for the automated algorithm were created for each test subject using a leave-one-out approach, where both the self-shape–size priors on any one object (e.g., the left amygdala) and neighboring-objects (e.g., right amygdala, left–right hippocampi) priors are found from manual tracings of these structures on the other 11 test images. The initial training data were found from 12 three-dimensional MRI structural T1-weighted images. Thus, using the any one object to key on, we assumed it was independently related to the other three structures. Using PCA, we built the object self shape–size model of each object and the neighbor difference models between it and the other three structures. In this initial work, we used additional regularizing terms in the formulation of the objective function (different for each object), one to enforce smooth boundaries and one a scalar term related to the approximate volume of the object.
One of the results of running our neighbor-constrained MAP algorithm on a anatomical MR image using manually traced priors is shown in Fig. 4. We report the errors of the application of this algorithm on a set of 12 subjects in Table 1. The mean and standard deviation for each structure over the 12 subjects are also shown in Table 1. Virtually all the boundary points lie within 1 or 2 mm of the manual segmentation. Note that we have begun to extend this approach with the inclusion of coupled intensity-appearance–shape priors (Yang and Duncan, 2003).
As part of our overall approach to structural brain analysis, we have been actively investigating the use of wavefront-level set-based strategies to estimate white matter fiber connectivity in the human brain (Jackowski et al., 2004). In this work, we first employ an anisotropic version of the static Hamilton–Jacobi (HJ) equation, and solve it by a sweeping method to obtain accurate front arrival times and determine connectivity. We briefly describe each part of this strategy below and present some early results as well.
White matter connectivity can be viewed as an instance of the minimum-cost path problem in an oriented-weighed domain. One would like to find a fiber path that minimizes the cumulative travel cost from a starting point A to some destination point B in the white matter. Because of the directionality of the tensor field, the cost function τ or its reciprocal speed F = 1/τ, is anisotropic, as it is a function of both position P(s) and direction P′(s). The minimum cumulative cost at x is defined as:
where L is pathway length, and the starting and ending points are given by P(0) = A and P(L)= x. A solution to Eq. (6) also satisfies the wave propagation equation:
which describes a wavefront propagating with speed 1/τ. In continuous space, solutions to Eq. (7) are given by the Hamilton–Jacobi (HJ) equations. While a classical solution to Eq. (7) may not exist, the viscosity solution is commonly sought using numerical approximations (Kao et al., 2002, 2003; Sethian and Vladimirsky, 2001).
Once the evolution Eq. (7) is solved for all points in the domain, one can use the resulting arrival times and find a solution for Eq. (6). The minimum-cost path between point A and an arbitrary point B in the white matter then becomes a solution to:
given X(0) = B. This optimal path can be constructed by integrating Eq. (8) at point B back to the seed point A using standard techniques.
To trace connectivity in the white matter, we use the entire tensor in our propagation model to avoid the possible misclassification of the principal eigenvector in oblate tensor regions, which may lead to wrong assignment of front arrival times. We design our wavefront to evolve from a seed point A, T(A) = 0, at a speed governed by a function of the diffusivity magnitude d in the front normal direction :
where F is designed to be large for large diffusivity and to slow the front rapidly when the diffusivity decreases. In addition, F accounts for the local anisotropy to stop the front from advancing into areas such as the ventricles or gray matter. In this way, we let the speed vary locally according to the tensor profile, descriptive of the underlying tissue structure.
The propagation equation belongs to a family of static Hamilton–Jacobi equations described by:
where Ω is the domain in 3, V(x) = 1 and q(x) is a function prescribing boundary condition values T(A)= q(A) = 0. While the propagation equation can be reformulated as a time-dependent HJ equation and solved by recovering each zero-level set, it is more convenient and less computationally expensive to model it as a static problem and determine arrival times instead.
Hamiltonians such as these cannot be correctly solved by isotropic propagation methods. However, carefully crafted methods have been devised (Kao et al., 2002, 2003; Sethian and Vladimirsky, 2001) to construct accurate solutions for anisotropic equations. We use a Lax–Friedrichs (LF) discretization of our Hamiltonian and employ a nonlinear Gauss–Seidel updating scheme (Kao et al., 2002) to solve the propagation equation. No minimization is required when updating an arrival time, and thus it is very easy to implement. Details on the LF sweeping (LFS) scheme can be found in Kao et al. (2003).
Fig. 5a shows a synthetic model consisting of two fiber bundles, A and B, oriented along helical paths that cross each other at their middle section. The background was filled with nearly isotropic tensors and diffusion-weighted images were created with added Gaussian noise. The signal-to-noise ratio is roughly half what is encountered in real diffusion MR scans, considering a single acquisition average. Fig. 5b depicts a close-up of the fiber-crossing region where oblate tensors resulting from the crossing are found. Streamlining techniques fail in such regions because of the oblate tensors in the fiber-crossing region. Using our method, we reconstructed the pathways properly without deviating the fiber trajectories. One of the resulting tensor images is shown in Fig. 5c. A seed point A1 was fixed at the bottom of bundle A (Fig. 5c) and our wavefront (γ = 2) was propagated using the LFS method on the tensor images. Points A2,B1, and B2 were fixed at the extreme ends of each bundle (Fig. 5c) and corresponding pathways A2A1, B1A1, and B2A1 were traced on T images using a Runge–Kutta 4th-order integration. Fig. 5d illustrates the resulting connectivity pathways embedded in the arrival time images and corresponding arrival isocurves (up to time 100). Darker areas in the maps show earlier arrivals. The fiber-crossing region did not prevent pathways connecting different branches or the same branch (A2A1) from being recovered. To assess the variability of the extracted paths, we propagated the same front in the diffusion tensor image without added noise and then computed the mean distance between corresponding pathways. The mean distance for all paths under noise σ2 = 0.05 was 1.06 voxels; for σ2 = 0.15 was 1.22 voxels; and for σ2 = 0.20 was 1.65 voxels. At the intersection, the maximum distance to the original pathway under noise variance σ2 = 0.20 was an average of 1.34 pixels. The largest errors, an average of 2.67 pixels, occurred not at the intersection itself, but rather in the remaining trajectories after propagating through the crossing. We attribute this deviation to the added noise in the region of the intersection, thereby altering the original trajectory. Nonetheless, the recovered pathways remained close to their original trajectories in spite of the added noise. Variations in the seed point location (A1) result in different pathway solutions, depending on the local tensor field directionality.
We applied our method on real human data, fixing the seed point in the splenium of the corpus callosum (Fig. 6a). We then propagated our wavefront throughout the image using the LFS method. Fig. 6a depicts the resulting arrival time level sets between 0 and 500. All pathways between points on the white matter boundary and the point in the splenium were traced using the map of arrival times. Fig. 6b shows the resulting 20,817 pathways, colored by the fractional anisotropy at each point, where brighter points represent higher anisotropy. Thus, we can observe the main routes of connection between various brain regions and the splenium. Not all connections shown in Fig. 6 represent true anatomical pathways. In future work, we will incorporate a metric to rate their anatomical likelihood.
We have developed a method that builds on and extends the robust point matching framework (RPM) previously developed by our group (Chui et al., 2001, 2003). RPM was originally presented in the context of joint estimation of rigid transformations (affine and piecewise-affine) and correspondences using “softassign” and deterministic annealing (Rangarajan et al., 1997a,b). Developed in an optimization framework, this previous work models point matching as a linear assignment-least squares problem. The algorithm estimates fuzzy correspondences while enforcing the one-to-one constraint. A deterministic annealing optimization framework gradually reduces the fuzziness of the correspondence in a controlled manner. Progressively refined transformation parameters are estimated by the resulting alternating update algorithm. The framework was first extended to the nonlinear registration case by Chui and Rangarajan (2000) and was evaluated on a synthetic data set as has been reported (Chui et al., 2001, 2003). The algorithm was shown to map landmarks to an accuracy of the order of 1–2.5 mm on average, which is superior to the typical performance of intensity based methods in the cortex (Hellier et al., 2003). However, this algorithm was hampered by (a) the limitation to a relatively small number of points and (b) the lack of outlier rejection, as a result of a clustering scheme used to increase the actual number of points.
As we have reported (Papademetris et al., 2003), we have extended the robust point matching algorithm to address both of these issues. We explicitly use the measure of “outlierness” estimated in the correspondence stage of this algorithm as a weight in the transformation estimation step to account for outliers in the template. Advanced numerical and graphics techniques such as sparse matrices and proper search strategies are also employed to enable the algorithm to use many points, which is important in the context of nonrigid brain registration.
We present here a modified form of the standard RPM methodology (Chui et al., 2003). The registration procedure consists of two alternative steps: (i) the correspondence estimation step and (ii) the transformation estimation step. In the following discussion, we label the reference point set as X and the transform point set as Y. The goal of the registration is to estimate the transformation G : X → Y . We label Gk the estimate of G at the end of iteration k. G0 is the starting transformation that can be the identity transformation. The whole process is embedded in a deterministic annealing framework, where the temperature T is progressively lowered, gradually decreasing the fuzziness of the correspondence and the smoothness of the transformation.
Given the point sets X and Y, we estimate the match matrix M, where Mij is the distance metric between points Gk(Xi) and Yj. The standard distance metric is defined as:
where |Xi – Yj| is the Euclidean distance between points Xi and Yj, and T is the temperature term that controls the fuzziness of the correspondence. Further, we restrict the correspondences to resemble a linear assignment problem—this requires that the rows and columns of M must sum to 1. To handle outliers, we introduce an outlier column C and an outlier row R, where Ci is a measure of the degree of “outlierness” of a point in the reference point set Xi and Rj is the same for a point in the transform point set Yj. C and R are initialized with constant values and then, using an iterative procedure (Rangarajan et al., 1997a,b), the values of M, C, and R are normalized to satisfy Eqs. (12) and (13).
Once the normalization is completed, we can compute the corresponding points as follows. Let Vi be the corresponding point to Xi and wi, the confidence in the match. Then Vi is defined as a normalized weighed sum of the points Yj where the weights are the elements of the match matrix M.
Note that a point that has a high value in the outlier column C will have low confidence and vice versa.
The transformation is estimated simply using a regularized weighted least squares fit between Xi and Vi as follows:
where S(g) is a regularization functional (e.g., a bending energy function) weighted by a function of the temperature f(T). This last weighting term is used to decrease the regularization as we approach convergence.
As an initial evaluation of the methodology, we tested it on a set of magnetic resonance three-dimensional SPGR images, with isotropic (1.2 mm)3 voxels, acquired using a GE 1.5-T scanner on a sample of 39 normal male controls (age 23.5 ± 10.2, IQ in the average range). We segmented the gray matter using our coupled level set algorithm (Zeng et al., 1999a,b) that was described earlier. Next, the gray matter ribbon surrounding the fusiform gyri was isolated using manual tracing and points from the gray matter surfaces were extracted and used as an input to our extended RPM algorithm (Papademetris et al., 2003), together with points from the cortical surfaces. Nonrigid registration using RPM to a separate normal template brain image yielded a transformation parameterized as a dense displacement field. An analysis of the registration results showed that the average errors were (a) 1.1 mm for points on the outer cortical surface and (b) 1.4 mm for points on the fusiform surfaces. These errors were determined by computing the distance from points on the warped reference surface to the corresponding points on the target surface, thus ensuring one-to-one correspondences.
To illustrate the utility of the approach, we performed a tensor-based morphometric study of the structural differences between normal controls and Williams syndrome subjects. Jackowski and Schultz (2003, in press) previously reported that the distances between the ends of the central sulci and the midline were larger in Williams patients compared to normal controls once total brain size was factored out. In collaboration with the authors of this study, we attempted to reproduce this finding using our point-based registration method.
In this work, three-dimensional SPGR magnetic resonance images with isotropic (1.2 mm)3 voxels were acquired using a GE 1.5-T scanner on a sample of 27 Williams syndrome subjects and 20 controls (groups were matched in terms of gender and age). We manually traced the superior frontal sulci and the central sulci (as reported in Jackowski and Schultz, 2003, in press) that were used as inputs to our extended RPM algorithm (Papademetris et al., 2003), together with points from the cortical surfaces, registering to a normal template image.
A tensor-based deformation measure (the determinant of the Jacobian of the displacement field) was computed from the displacement field to provide a quantitative metric of the local expansion–contraction needed for each voxel to register to the template. The deformation measure was sampled every 1.5 mm and smoothed (FWHM = 3.1 mm), and groups were then compared by creating a t-map (adjusted for overall brain size scaling differences) keeping all values more significant than P < 0.001 (uncorrected). Differences shown on the t-map away from the sulci may not be reliable because sulcal points were not used there to define local anatomy. However, as shown in Fig. 7, and in accordance with the earlier reported finding, we observe a bilateral morphometric difference between the ends of the central sulci and the midline. We note that using an intensity based nonlinear registration yielded less significant differences only in the left hemisphere as a result of inaccuracies in the registration of the sulci.
While we believe the registration results near the regions of interest such as the central and superior frontal sulci are accurate, we feel that the results in other regions are more variable. To address this, we are currently working on an integrated feature–intensity nonrigid registration strategy (Papademetris et al., 2004).
The reliable and accurate derivation of parametric information from structural MR images useful for quantifying differences related to a range of neuropsychiatric disorders requires several stages of processing (some of which are not mentioned here). In this paper, we have described four related image analysis efforts aimed at the quantification of cortical and subcortical gray matter and white matter interconnections, along with a strategy for registering multiple subject structural image information. Our overall approach to these problems is to focus on different specific regions rather than attempting to analyze the entire brain at once. In this context, we have often found it advantageous to focus on methods that rely on Bayesian reasoning strategies that permit the incorporation of spatial relationship priors representing the geometrical–biological variability across a population. In addition, level set and wavefront propagation strategies, such as were used for subcortical and cortical segmentation, as well as for diffusion image analysis, have proven to be useful and generally robust methodologies for our purposes. Finally, we note that our registration strategy focuses on regionally segmented structure, an approach that we have found useful in a variety of analyses. Our future directions are aimed at integrating functional and structural information within this same basic framework.
We would like to acknowledge the help of collaborators Robert Schultz, Andrea Jackowski, and Todd Constable as well as the support of NIH-NINDS grant R01NS035193 and NIH-NIBIB grant R01EB000311.