Abstract
Nonlinear registration of 3D surfaces is important in many medical imaging applications, including the mapping of longitudinal changes in anatomy, or of multi-subject functional MRI data to a canonical surface for comparison and integration. To register 3D surfaces, such as the cortical surface of the brain, one approach is to transform them first to planar or spherical objects. Internal landmarks can then be matched on these simpler parameter domains. Here we study the diffeomorphic matching of landmarks on the sphere. Our method builds on the level set technique of Leow et al. [1] for the plane. Both forward and backward matching terms are included, thus ensuring the invertibility of the representation. We demonstrate our technique on a pair of lines on the sphere. The overall approach improves on earlier work in cortical matching by allowing the matching energy to be relaxed along sulcal landmarks, minimizing distortion, and also enables point and curve landmarks to be aligned in the same general framework as densely-defined scalar fields, such as curvature or cortical thickness maps.
1. INTRODUCTION
Image registration to align anatomical structures is of fundamental importance in medical imaging. Such methods have enabled the pooling and comparison of data across subjects, yielding valuable information on brain anatomy and function in healthy and diseased individuals. These techniques have also been used to monitor changes within a single individual due to normal brain growth [2], as well as to map progressive changes due to aging and illness.
Although the cortical surface is the focus of most brain mapping studies, its anatomy varies widely between individuals. Because of this, the comparison of cortical imaging data is often performed via non-linear registration techniques. For many applications, sets of sulci or gyri that are common to all normal brains are used as landmarks to guide the transformation, and the cortices are registered with the constraint that the corresponding structures be matched [3].
Flat or spherical representations are often used to visualize the entire cortical surface [4], [5]. While cuts are required to obtain a planar map, the topology of the cortex can be preserved by using a spherical transformation. Furthermore, for shape comparisons, artificial boundary conditions need to be imposed at the edges of the planar region to which the cortex is flattened during the registration process. Cuts may also need to be introduced to avoid excessive distortion in some regions when flattening the cortex to a planar domain. No such cuts or boundary constraints are required for the mapping between spheres.
In [6], Christensen et al. proposed a viscous fluid nonlinear registration model. A template image T is evolved over time t into a study S using a fluid partial differential equation (PDE). The model allows for the stress of the transformation to relax over time, thereby permitting large deformations.
Within this framework, the transformation ϕ satisfies
| (1) |
where
The final solution is given by ϕ(x, 1). Here the coordinate x is taken in a fixed reference frame, while the image deforms in time. Thus v(ϕ(x, t)) defines a velocity for the particles of the fluid in which the image is considered to be embedded. Joshi and Miller [7] (see also [8] ) constructed diffeomorphic solutions of this transport equation for the landmark matching problem in Euclidean space ϕ : [0,1]3 → [0,1]3. Their map was chosen to minimize the smoothness constraint ∫||Lv(x⃗, t)||2dx⃗dt, where L is a differential operator, while minimizing the distance between landmarks at the final time-point.
The diffeomorphic landmark matching procedure was first extended to the sphere by Bakircioglu et al. [9], and the method was generalized to a coordinate independent framework in [10]. In the latter approach, the smoothness constraint is written at each time point in terms of a propagator K between points on the trajectories of the landmarks. For point-to-point matching, ϕ(xi, 1) = yi, where the yi are the study landmarks, and xi are corresponding landmarks in the template. The method was tested for point landmarks on the sphere.
For landmarks representing more general shapes (e.g. curves and subregions on the surface), we would like to be able to match the entire shape rather than enforce point-by-point matching. However, if the landmarks are not constrained to match exactly, the numerical estimation of the reverse transformation S → T will not necessarily yield the inverse of the original mapping from T → S. For elastic landmark matching in Euclidean space, Christensen and Johnson [11] (see also [12]) solved this problem by simultaneously calculating both the forward and backward transformations, while requiring them to be inverses of each other.
In the context of large deformations, Leow et al. [1] used the smoothness constraint described in [7] in the plane, while replacing their matching term by a sum of the backward and forward distance functions at the final positions of the template landmarks. This method was shown to significantly improve the results, relative to unidirectional matching. Furthermore, the technique described by those authors allows relaxation of the transformation energy along curved landmarks (and therefore less distortion), rather than using an explicit point-by-point correspondence.
In this work, we extend the method proposed by Leow et al. to the sphere. This paper is organized as follows. In the next section, we review the large deformation diffeomorphism method of [10] in our context, and the new landmark matching term is derived. Our method is first applied to a pair of lines on the sphere. All our results are compared to ones obtained using the method in [10].
2. METHOD
2.1. Landmark matching on the sphere
In this paper, we consider transformations between 2D spheres S2 embedded in ℜ3. The landmark problem for our case is defined as follows: consider a structure which exists on both images, which may be any combination of points, curves or surfaces on S2. This shape may be discretized as a set of landmark points xi in the template, and yi in the study. For the case of point-to-point landmark matching, we would like a registration procedure to match the points xi to the yi as closely as possible, while enforcing a smoothness constraint for the transformation.
Following Glaunes et al. [10], we use the regularization energy E(v) over the velocity field:
| (2) |
Here <.,.> is the usual dot product in ℜ3 for vectors in the tangent plane at x, while L is the Laplace-Beltrami operator dδ + δd on the vector field v (see for instance [13]). This energy is a good choice since
| (3) |
was shown by those authors to be a geodesic distance in the group of deformation maps generated by Eq. 1.
Point-to-point matching of the landmarks is enforced by constraining the deformation ϕ(xi, 1) to be equal to the yi’s. One problem with this definition is that it does not take into account the fact that some of the landmark points are discretizations of higher dimensional structures such as curves or surfaces. In that case, the matching should not be done point by point, but rather for the entire structure at once. Some relaxation of the matching energy is then permitted to occur along the curved landmarks, while still fitting all of one curve to the other. Thus, a more natural way to proceed consists of calculating the unsigned distance functions Ty(ϕ(., 1)) between the ϕ(xi, 1) and the associated structure in the study. This method was used in the plane by Leow et al. [1], and is extended to the sphere in this paper.
Furthermore, if we are to use distance functions, the matching needs to be enforced in both the backward and forward directions. The bi-directionality of the constraint ensures full matching of the geometry. To see why this is true, consider the choice an open curve as a landmark. In this case, it is possible to obtain a minimum for the cost function if the template curve is deformed to a subset of its counterpart in the study. This problem disappears if the mapping is done in both directions at once. Thus, denoting the inverse transformation by ϕ−1, the problem consists of finding solutions v such that
| (4) |
is minimal.
Glaunes et al. derived a coordinate independent method to calculate the regularization term Eq. 2. A simplified description of it is given here. The basic idea consists of computing the deformation in a basis of eigenvectors of the Laplacian operator. Initially, estimated trajectories xi(t) for the deformation between the initial and final landmark points are chosen. These may be for instance arcs of great circles between each of the associated initial and final landmark points. Given these trajectories, the energy E(v(t)) at time t can be written in a simple form in the eigenbasis of L, as described more explicitly below. The time component of the energy integral is then computed numerically. Eq. 4 is minimized via steepest descent until the optimal trajectories are found.
More explicitly, the central piece of the calculation is the value of the spatial part of the energy integral. If we denote the eigenvectors of L by Elm and its eigenvalues by λm, a simple calculation yields the energy integral at time t in the Elm basis as
| (5) |
This energy is to be minimized given the constraint
where vi(t) denotes the velocity vectors along the estimated trajectories xi(t) between the landmarks.
The velocity field vmin which solves this problem is found in [10] to be
| (6) |
where αi(t) is defined through the constraint values for the velocities
Furthermore, by replacing this value for vmin into Eq.5, the minimum of the energy at time t is
| (7) |
Thus, the problem reduces to finding the trajectories for the landmarks which minimize the time integral of Eq.7 (see also [7], for an analogous proof for landmark paths in 3D).
In our case, in order to obtain both the forward and backward distance functions, the quantity that we minimize is the following
| (8) |
where β is the scale of the distance function contribution (i.e. the weighting on the data driven term relative to the regularizer).
3. RESULTS
To understand the differences in performance between our distance function matching approach and the exact matching one on the sphere, we evaluated the geodesic distance for a pair of curved lines on the sphere. The same landmarks were matched by each method and the results compared. Twenty landmark points were selected on each line. Figures 1 display the trajectories for both point-to-point and distance function methods. The energy of the transformation for both cases is mapped in Figure 2 and the total geodesic distance for each transformation was found to be 7.35x104 and 2.24x104, respectively. Thus, the distance function mapping does significantly better than point-to-point matching. Notice also that for distance function matching, the final template points do not necessarily fall on the study landmarks.
Fig. 1.
Trajectories using point-to-point matching (left) and distance function matching (right). The leftmost curve was chosen as the template, while that on the right is the study. Initial template points (x), final template points (o) and study points (+) are also shown.
Fig. 2.
Energy of the velocity field vs time for the case of point-to-point matching (left) and distance function matching (right).
4. CONCLUSION
In this paper, we described a new landmark matching method on the sphere, within the large deformation diffeomorphism framework. The use of distance functions as part of the matching constraints allowed us to create mappings that match the whole landmark curves, allowing the transformation energy to be relaxed along them. This produces correspondence maps with less spatial distortion. Furthermore, consistent backward and forward matching constraints were imposed. We applied our results to match a pair of lines and compared them to the point-to-point matching case developed in [10]. Our method showed significant improvements over the one found in that paper, partly because it allows for the minimization of distortion when matching submanifolds lying in the sphere. We are currently validating this approach for matching cortical surfaces in 3D, in studies of brain growth and degeneration as well as multi-subject functional imaging studies, where cortical registration is especially important.
Acknowledgments
We thank Chui-Yen Kao and Brian Thorndyke for useful discussions. Algorithm development was supported by NIH grants R21 EB01651, R21 RR019771, P50 AG016570, U54 RR021813, and P41 RR13642 (to P.T.),
References
- 1.Leow A, Yu CL, Lee SJ, Huang C, Nicolson SR, Hayashi KM, Protas H, Toga AW, Thompson PM. Brain structural mapping using a novel hybrid implicit/explicit framework based on the level-set method. Neuroimage. 2005;24:910–927. doi: 10.1016/j.neuroimage.2004.09.022. [DOI] [PubMed] [Google Scholar]
- 2.Thompson PM, Giedd JN, Woods RP, MacDonald D, Evans AC, Toga AW. Growth pattern in the developing brain detected by using continuum mechanical tensor maps. Nature. 2000;404:190–193. doi: 10.1038/35004593. [DOI] [PubMed] [Google Scholar]
- 3.Thompson PM, Hayashi KM, Sowell ER, Gogtay N, Giedd JN, Rapoport JL, de Zubicaray GI, Janke AL, Rose SE, Semple J, Doddrell DM, Wang YL, van Erp TGM, Cannon TD, Toga AW. Mapping cortical change in alzheimer’s disease, brain development, and schizophrenia. NeuroImage. 2004 Sep;23(Suppl):2–18. doi: 10.1016/j.neuroimage.2004.07.071. [DOI] [PubMed] [Google Scholar]
- 4.Fischl B, Sereno MI, Tootell RBH, Dale AM. High-resolution inter-subject averaging and a coordinate system for the cortical surface. Hum Brain Mapp. 1999;8:272–284. doi: 10.1002/(SICI)1097-0193(1999)8:4<272::AID-HBM10>3.0.CO;2-4. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 5.Hurdal MK, Stephenson K. Cortical cartography using the discrete conformal approach of circle packings. NeuroImage, Special Issue on Mathematics in Brain Imaging. 2004 Sep; doi: 10.1016/j.neuroimage.2004.07.018. [DOI] [PubMed] [Google Scholar]
- 6.Christensen GE, Rabbitt RD, Miller MI. Deformable templates using large deformation kinematics. IEEE Trans Image Process. 1996;5:1435–1447. doi: 10.1109/83.536892. [DOI] [PubMed] [Google Scholar]
- 7.Joshi SC, Miller MI. Landmark matching via large deformation diffeomorphisms. IEEE Trans Image Process. 2000;9:1357–1370. doi: 10.1109/83.855431. [DOI] [PubMed] [Google Scholar]
- 8.Miller MI. The emerging field of computational anatomy. Neuroimage, Special Issue on Mathematics in Brain Imaging. 2004 Sep; doi: 10.1016/j.neuroimage.2008.10.033. [DOI] [PubMed] [Google Scholar]
- 9.Bakircioglu M, Joshi S, Miller MI. Landmark matching on brain surfaces via large deformation diffeomorphisms on the sphere. Proc SPIE Medical Imaging 1999: Image Processing. 1999;3661:710–715. [Google Scholar]
- 10.Glaunés J, Vaillant M, Miller MI. Landmark matching via large deformation diffeomorphisms on the sphere. J of Math Imaging and Vision. 1994;20:179–200. [Google Scholar]
- 11.Johnson HJ, Christensen GE. Consistent landmark and intensity based image registration. IEEE Trans Med Imaging. 2002;21:450–461. doi: 10.1109/TMI.2002.1009381. [DOI] [PubMed] [Google Scholar]
- 12.Christensen GE, Johnson HJ. Consistent image registration. IEEE Trans Image Process. 2001;20:568–582. doi: 10.1109/42.932742. [DOI] [PubMed] [Google Scholar]
- 13.Choquet-Bruhat Y, DeWitt-Morette C, Dillard-Bleick M. Analysis, manifolds and physics, part 1: basics. Elsevier Science Publishers B. V; Amsterdam, Netherlands: 1991. [Google Scholar]