Image registration driven by combined probabilistic and geometric descriptors

Deformable image registration in the presence of considerable contrast dierences and large-scale size and shape changes represents a signicant challenge for image registration. A representative driving application is the study of early brain development in neuroimaging, which requires co-registration of images of the same subject across time or building 4-D population atlases. Growth during the first few years of development involves significant changes in size and shape of anatomical structures but also rapid changes in tissue properties due to myeli-nation and structuring that are reflected in the multi-modal Magnetic Resonance (MR) contrast measurements. We propose a new registration method that generates a mapping between brain anatomies represented as a multi-compartment model of tissue class posterior images and geometries. We transform intensity patterns into combined probabilistic and geometric descriptors that drive the matching in a diffeomorphic framework, where distances between geometries are represented using currents which does not require geometric correspondence. We show preliminary results on the registrations of neonatal brain MRIs to two-year old infant MRIs using class posteriors and surface boundaries of structures undergoing major changes. Quantitative validation demonstrates that our proposed method generates registrations that better preserve the consistency of anatomical structures over time.

Image Registration Driven by Combined Probabilistic and Geometric Descriptors ? Linh Ha1, Marcel Prastawa1, Guido Gerig1, John H. Gilmore2, Claudio T. Silva1, Sarang Joshi1 1 Scientic Computing and Imaging Institute, University of Utah 2 Department of Psychiatry, University of North Carolina Abstract. Deformable image registration in the presence of consider- able contrast dierences and large-scale size and shape changes repre- sents a signicant challenge for image registration. A representative driv- ing application is the study of early brain development in neuroimaging, which requires co-registration of images of the same subject across time or building 4-D population atlases. Growth during the rst few years of development involves signicant changes in size and shape of anatomi- cal structures but also rapid changes in tissue properties due to myeli- nation and structuring that are re ected in the multi-modal Magnetic Resonance (MR) contrast measurements. We propose a new registration method that generates a mapping between brain anatomies represented as a multi-compartment model of tissue class posterior images and ge- ometries. We transform intensity patterns into combined probabilistic and geometric descriptors that drive the matching in a dieomorphic framework, where distances between geometries are represented using currents which does not require geometric correspondence. We show pre- liminary results on the registrations of neonatal brain MRIs to two-year old infant MRIs using class posteriors and surface boundaries of struc- tures undergoing major changes. Quantitative validation demonstrates that our proposed method generates registrations that better preserve the consistency of anatomical structures over time. 1 Introduction Image registration is a basic task in dening a standard space for analyzing populations that change over time, which is essential to determine development in normal growth and neurological disorders. The growth process can involve large-scale size and shape changes, as well as changes in tissue properties and appearance. These factors pose signicant challenges in image registration, as image intensities need to be interpreted dierently at dierent stages. A strong example is the human brain at early development stages. A driving research question is to determine the process of white matter myeli- nation, which manifests as two distinct white matter appearance patterns pri- marily during the rst year of development. Other clinical research questions are ? Supported by NIH grants R01 EB000219-NIH-NIBIB and Conte Center MH064065, and NSF grant CNS-0751152. UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript related to nding a link between cognitive development and the rapid, locally varying growth of specic anatomical structures. To approach these questions, a robust registration method is necessary for mapping longitudinal brain MRI to a reference space so that we can perform reliable analysis of the tissue property changes re ected in the MR measurements. Knickmeyer et al.[1] showed that the total brain volume grows by 100% the rst year and 15% the second year, whereas the cerebellum shows 220% volume growth for the rst and another 15% for the second year, indicating the very dierent growth rates of dierent anatomical structures. Through regression on shape representations, Datar et al.[2] illustrated that the rapid volume changes are also paralleled by signi- cant shape changes which describe the dynamic pattern of localized, nonlinear growth. These challenges require a method that does not rely on raw intensity measurements, while also being capable of estimating large structural deforma- tions. Xue et al.[3] addressed these issues by proposing a registration scheme for neonatal brains by registering in ated cortical surfaces extracted from the MRI. We propose a new registration framework for longitudinal brain MRI that makes use of the underlying anatomies, which are represented by both class pos- teriors and boundary surfaces. This framework is able to match internal anatom- ical regions and simultaneously preserving a consistent mapping for the bound- aries of relevant anatomical objects.We show results of registering neonatal brain MRI to 2-year old brain MRI of the same subjects obtained in a longitudinal neu- roimaging study. Our method consistently provides transformations that better preserve time-varying structures than obtained by intensity-only registration. 2 Registration between Anatomies Fig. 1. Overview of the proposed registration method that can handle large defor- mations and dierent contrast properties, applied to mapping brain MRI of neonates to 2-year olds. We segment the brain MRIs and then extract equivalent anatomical descriptors by merging the two dierent white matter types present in neonates. The probabilistic and geometric anatomical descriptors are then used to compute the trans- formation h that minimizes the distance between the class posterior images, as well as the distance between surfaces represented as currents. UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript We propose a registration method that makes use of the underlying anatomy in the MR images. Fig. 1 shows an overview of the registration process. We begin by extracting anatomical descriptors from the images, followed by computing a transformation that minimizes the distance between the anatomical descriptors. 2.1 Anatomical Descriptors We represent brain anatomy as a multi-compartment model of tissue classes and manifolds. We associate each position x with a vector of tissue probability den- sities. In a given anatomy, we capture the underlying structures by estimating, for each image, the class posterior mass functions associated with each of the classes. Given as the underlying coordinate system of the brain anatomies, each anatomy Ai=1; ;N is represented as Ai = fpi;c=1(x); ; pi;c=Nc (x);Mi;j=1(2); ;Mi;j=Ns(2) g (1) where Nc is the number of probability images, Ns is the number of surfaces, pc(x) is the class posterior for tissue c at location x, and Mj(2) are 2-dimensional submanifolds of (surfaces). The classication of brain MR images with mature white matter structures into class posteriors are well studied. We extract the posteriors from 2-year old brain MR images using the segmentation method proposed by van Leemput et al.[4]. The method generates posterior probabilities for white matter (wm), gray matter (gm), and cerebrospinal uid csf. These probabilities can then be used to generate surfaces from the maximum a posteriori tissue label maps. The classication of neonatal brain MR images is challenging as the white matter structure undergoes myelination, where the bers are being covered in myelin sheathes. Several have proposed methods that make use of prior informa- tion from an atlas or template that takes into account the special white matter appearance due to myelination [5] [6]. We use the method described by Prastawa et al.[6] for extracting the tissue class posteriors of neonatal brain MRI which includes for myelinated wm, non-myelinated wm, gm, and csf. These can then be used to create an equivalent anatomy to the 2-year old brain by combining the two white matter class probabilities and surfaces. For the results in this paper, we use the probabilities fpwm(x); pgm(x); pcsf (x)g and we use the surfaces of white matter, gray matter, and cerebellum. The cere- bellum surfaces are generated from semi-automated segmentations that are ob- tained by anely registering a template image followed by a supervised level set segmentation. The cerebellum has a signicant role in motor function and it is explicitly modeled as it undergoes the most rapid volume change during the rst year of development and thus presents a localized large-scale deformation. 2.2 Registration Formulation Given two anatomies A1 and A2, the registration problem can be formulated as an estimation problem for the transformation h that minimizes ^h = arg min h E(h A1;A2)2 + D(h; e)2 (2) UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript where h A1 is the transformed anatomy, E(; ) is a metric between anatomies and D(; e) is a metric on a group of transformations that penalizes deviations from the identity transformation e. We dene distance between anatomies E by dening a norm on an anatomy as a combination of the L2 norm on the class posteriors and a Reproducing Kernel Hilbert space norm on the manifolds dened as \currents" through Glaunes [7]. The currents norm does not require geometric correspondence and thus can be used to register manifolds with dierent resolutions. For an oriented surface M(2) in R3 the norm [M(2)] is the vector valued Borel measure corresponding to the collection of unit-normal vectors to M(2), distributed with density equal to the element of surface area ds and can be written as (x)ds(x), where (x) is the unit normal and ds(x) is the surface measure at point x. When M(2) is a discrete triangular mesh with Nf faces, a good approximation of the norm can be computed by replacing [M(2)] by a sum of vector-valued Dirac masses: k[M(2)]k2 k = XNf f=1 XNf f0=1 h(f); (f0)i k(c(f); c(f0)); (3) where k(; ) is a shift-invariant kernel (e.g., Gaussian or Cauchy), Nf is the number of faces of the triangulation, and for any face f, c(f) is its center and (f) its normal vector with the length capturing the area of each triangle. Having dened the norm on probability images and surfaces, the dissimilarity metric between anatomies [A1] 􀀀 [A2] 2 k is given by: XNc c=1 Z jp1;c(x) 􀀀 p2;c(x)j2dx + XNs j=1 k[M1;j(2) [ (􀀀M2;j(2))]k2 k (4) where the distance between two surface currents k[M1;j(2) 􀀀M2;j(2)]kk = k[M1(2) [ (􀀀M2(2))]kk is computed as the norm of the union between surface M1(2) and surface M2(2) with negative measures. We use the large deformation framework [8] that generates dense defor- mation maps in Rd by integrating time-dependent velocity elds. The ow equation is given by @hv(t;x) @t = v(t; hv(t; x)); with h(0; x) = x, and we dene h(x) := hv(1; x), which is a one-to-one map in Rd (dieomorphism). We dene an energy functional that ensures the regularity of the transformations on the velocity elds: kv(t; )k2 V = R Rd hLv(t; x); Lv(t; x)i dx; where L is a dierential operator acting on vector elds. This energy also denes a distance in the group of dieomorphisms: D2(h; e) = inf v;pv(1;)=h Z 1 0 kLv(t)k2 V dt: (5) The registration optimizations in this paper are performed using a greedy ap- proach by iteratively performing gradient descent on velocity elds and updating the transformations via an Euler integration of the O.D.E. At each iteration of UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript the algorithm the velocity eld is calculated by solving the p.d.e Lv = F(h) where F(h) is the variation of [h A1] 􀀀 [A2] 2 k with respect to h. This varia- tion is a combination of the variation of the L2 norm on the class posteriors and of the currents norm; computed using the gradient @k[M(2)]k2 k @xr : X fjxr2f @(f) @xr XNf f0=1 k(c(f0); c(f))(f0) + 2 3 XNf f0=1 @k(c(f); c(f0)) @c(f) (f0)t(f) (6) given that points fxr; xs; xtg form the triangular face f and its center c(f) = xr+xs+xt 3 and its area-weighted normal (f) = 1 2 (xs 􀀀 xr) (xt 􀀀 xr). 2.3 Ecient Norm Computation using Particle Mesh Approximation on GPU The major challenge of implementing the currents norm (Eq. 3) for realistic brain surfaces is the high computational cost to compute the dissimilarity met- ric of all pairs of surface elements, which is O(N2 f ) where Nf is the number of faces (can be up to millions). For computational tractability, Durrleman et al. [9] used a sparse representation of the surface based on matching pursuit algorithm. An ecient framework based on the standard fast Gauss transform [10] requires the construction and maintenance of the kd-tree structure on the y, which is slow on the current CPU model. Our method, however, exploits the Particle-Mesh approximation to reduce the complexity to M logM where M is the volume size of the embedded grid. The grid size is chosen as the size of the input images to limit the approximation error to the same order of matching term for the class posteriors. This approximation have been extensively studied in the cosmological N-body simulation literature (see Hockney and Eastwood [11] for details). Particle mesh framework shares the same computational grid with the class posteriors, moving the computation to the grid level and enabling an ecient parallel computation. Even with the particle mesh approximation of the norm computation, the total complexity of the method is still very high. On a high-end workstation with 8-CPU cores, a highly optimized multi-threaded implementation in C++ takes several hours for one matching pair, hence can not be used for parame- ter exploration and real-time analysis. Fortunately, this computational cost can be amortized using the massive parallel power and bandwidth provided by the graphics processing unit (GPU). Based on the GPU framework by Ha et al.[12], we developed an implementation that runs entirely on the GPU. The main bene- t of a GPU implementation is the ability to exploit parallel eciency of regular grid presentation. Computing the currents u and gradient between a surface with 160535 trian- gular faces and another with 127043 faces takes approximately 504 seconds on an AMD Phenom II X4 955 CPU, while it takes 0.33 seconds on an NVIDIA GTX 260 GPU. The speed gain is in order of three magnitudes over the equiva- lent CPU implementation using particle mesh, while the computing time for the UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript (a) (b) (c) (d) Fig. 2. Registration results of neonates mapped to two-year olds. From left to right: (a) neonatal T1 image after ane registration, (b) reference T1 image at 2-years, followed by (c) neonatal T1 after deformable mutual information registration using B-splines, and (d) after combined probabilistic and geometric registration. From top to bottom: subject 0012, 0102, 0106. exact norm on CPU is dicult to measure since it takes signicantly longer. The proposed algorithm typically converges in 1000 iterations, so on average it takes less than eight minutes to register two anatomies. This allows us to perform pa- rameter exploration and real-time analysis on a single desktop with commodity GPU hardware. The eciency of the GPU method also provides an opportunity to apply the algorithm for high quality atlas formation using our framework on a GPU cluster, which gives us the ability to perform statistical tests that are previously impossible due to excessive time requirements. 3 Results We have applied the registration method for mapping neonatal MRI scans to two-year MRI scans of the same subjects in ten datasets. The datasets are taken from an ongoing longitudinal neuroimaging study with scans acquired at ap- proximately two weeks, one year, and two years of age. Due to rapid early brain development, each longitudinal MR scan shows signicant changes in brain size and in tissue properties. For comparison, we also applied the standard inten- sity based deformable registration using mutual information (MI) metric and B-spline transformation proposed by Rueckert et al.[13] which has been applied for registering 1-year old and 2-year old infants. The T1 weighted images before and after registration using the dierent approaches for the rst three subjects are shown in Fig. 2. A quantitative study of the performance of the registration method is per- formed by measuring the overlap between the transformed segmentation maps of neonates to the segmentation maps of two-year olds. Since we consider the UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript Subject 0012 0102 0106 0121 0130 0146 0156 0174 0177 0180 White matter MI 0.329 0.155 0.281 0.384 0.379 0.230 0.257 0.300 0.350 0.301 P+G 0.497 0.397 0.442 0.453 0.482 0.414 0.461 0.440 0.478 0.442 Cerebellum MI 0.755 0.212 0.588 0.515 0.732 0.820 0.713 0.569 0.631 0.777 P+G 0.881 0.821 0.875 0.878 0.858 0.899 0.907 0.885 0.896 0.892 Fig. 3. Overlap measures comparing the registered segmentation maps against the reference segmentation maps for the white matter and cerebellum structure, obtained through deformable mutual information registration (MI) and our proposed method (P+G). Subject 0012 0102 0106 0121 0130 0146 0156 0174 0177 0180 MI on CPU 92 63 103 92 101 112 106 99 91 96 P+G on GPU 9 8 8 8 8 7 9 8 7 7 Fig. 4. Time elapsed, in minutes, for registration using deformable mutual information on CPU (MI) and our proposed approach (P+G) on GPU with 1000 iterations of gradient descent. segmentation maps at two years of age to be the standard, we use the following overlap metric: Overlap(h S0; S2) = jh S0 \ S2j jS2j (7) where h S0 is the transformed neonate segmentation map, S2 is the reference two-year segmentation map, and j j indicates the volume of a binary map. Please note that this metric gives considerably lower values than the standard Dice coecient. Fig. 3 shows the quantitative analysis for the white matter and cerebellum segmentation maps. Registration using both probabilistic and geo- metric descriptors provides consistently better results and are generally more stable for the structures of interest. In particular, our method better preserves the shape of the cerebellum, which has weak intensity boundaries in regions where it touches the cerebrum and thus cannot be registered properly using only image based information. Another signicant challenge is that the cere- bellar growth is distinctly dierent from the growth of neighboring structures. Registrations using our approach on the GPU takes 8 minutes on average, while registration on the CPU using mutual information metric and B-spline transfor- mation takes 100 minutes on average. Detailed time measures are listed in Fig. 4. We have also performed validation using only the probabilistic descriptor, which generates results that are less accurate compared to our method (particularly for the cerebellum) while more accurate than image registration using MI. 4 Conclusions We have proposed a registration framework that makes use of the probabilistic and geometric structures of anatomies embedded in the images. This allows us to UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript enforce matching of important anatomical features represented as regional class posteriors and tissue boundaries. Our framework allows us to register images with dierent contrast properties by using equivalent anatomical representations, and we have demonstrated results for registering brain MRIs with dierent white matter appearances at early stages of growth. The overlap validation measures in Fig. 3 show that geometric constraints, particularly for the cerebellum, is crucial for registering structures undergoing signicant growth changes. In the future, we plan to apply this framework in early neurodevelopmental studies for analyzing the eects of neurological disorders such as autism and Fragile X syndrome. The proposed registration framework is generic and independent of the application domain, it can thus be applied to any registration where one encounters large-scale deformation and dierent appearance patterns. References 1. Knickmeyer, R.C., Gouttard, S., Kang, C., Evans, D., Wilber, K., Smith, J.K., Hamer, R.M., Lin, W., Gerig, G., Gilmore, J.H.: A structural MRI study of human brain development from birth to 2 years. J. Neurosci. 28 (Nov 2008) 12176{12182 2. Datar, M., Cates, J., Fletcher, P., Gouttard, S., Gerig, G., Whitaker, R.: Particle based shape regression of open surfaces with applications to developmental neu- roimaging. In: MICCAI. Number 5762 in LNCS, Springer Verlag (2009) 167{174 3. Xue, H., Srinivasan, L., Jiang, S., Rutherford, M.A., Edwards, A.D., Rueckert, D., Hajnal, J.V.: Longitudinal cortical registration for developing neonates. In: MICCAI. (2007) 127{135 4. Van Leemput, K., Maes, F., Vandermeulen, D., Suetens, P.: Automated model- based tissue classication of MR images of the brain. IEEE Trans. Medical Imaging 18 (October 1999) 897{908 5. Wareld, S.K., Kaus, M., Jolesz, F.A., Kikinis, R.: Adaptive, template moderated, spatially varying statistical classication. MedIA 4(1) (Mar 2000) 43{55 6. Prastawa, M., Gilmore, J.H., Lin, W., Gerig, G.: Automatic segmentation of mr images of the developing newborn brain. MedIA 9(5) (2005) 457{466 7. Glaunes, J., Trouve, A., Younes, L.: Dieomorphic matching of distributions: a new approach for unlabelled point-sets and sub-manifolds matching. In: CVPR. (2004) 8. Miller, M, I., Younes, L.: Group action, dieomorphism and matching: a general framework. Int. J. Comp. Vis 41 (2001) 61{84 9. Durrleman, S., Pennec, X., Trouve, A., Ayache, N.: Sparse approximations of currents for statistics on curves and surfaces. In: MICCAI 2008. (2008) 10. Greengard, L., Strain, J.: The fast gauss transform. SIAM J. Sci. Stat. Comput. 12(1) (1991) 79{94 11. Hockney, R.W., Eastwood, J.W.: Computer Simulation Using Particles. Taylor and Francis (1989) 12. Ha, L., Kruger, J., Fletcher, P.T., Joshi, S., Silva, C.T.: Fast parallel unbiased dieomorphic atlas construction on multi-graphics processing units. In: EGPGV. (2009) 13. Rueckert, D., Sonoda, L., Hayes, C., Hill, D., Leach, M., Hawkes, D.: Nonrigid registration using free-form deformations. IEEE Trans Med Imag 18(8) (1999) 712{721 UU IR Author Manuscript UU IR Author Manuscript University of Utah Institutional Repository Author Manuscript