A large deformation diffeomorphic metric mapping solution for diffusion spectrum imaging datasets

Abstract

Spatial transformation for diffusion spectrum imaging (DSI) is an important step for group analyses of DSI datasets. In this study, we developed a transformation method for DSI datasets under the framework of large deformation diffeomorphic metric mapping (LDDMM), which is termed LDDMM-DSI. The proposed method made use of the fact that a DSI dataset is 6D, and generalized the original 2D/3D LDDMM algorithm to the 6D case with some modifications made for the DSI datasets. In this manner, the conventional reorientation problem that arises from transforming diffusion-weighted datasets was avoided by making the DSI datasets capable of being freely deformed in the q-space. The algorithm treated the data-matching task as a variational problem under the LDDMM framework and sought optimal velocity fields from which the generated transformations were diffeomorphic and the transformation curve was a geodesic. The mathematical materials and numerical implementation are detailed in the paper, and experiments were performed to analyze the proposed method on real brain DSI datasets. The results showed that the method was capable of registering different DSI datasets in both global structural shapes and local diffusion profiles. In conclusion, the proposed method can facilitate group analyses of DSI datasets and the generation of a DSI template.

Introduction

Diffusion weighted (DW) magnetic resonance imaging (MRI), by applying diffusion-sensitive gradients, is capable of detecting the Brownian motion of water molecules, and it has been considered to be a promising tool for investigating brain white matter (WM) in vivo. To capture the microstructure properties in more detail, acquiring more DW volume images using multiple diffusion-sensitive gradients is usually desirable. Diffusion tensor imaging (DTI) (Basser et al., 1994), one of the most elegant methods, can represent microscopic architectures using a second-order tensor, where the principal eigenvector is usually assumed to coincide with the orientation of the fiber bundle in a voxel. Nevertheless, the inherent Gaussian probability density function (PDF) assumption limits DTI to capture at most single fiber orientation in a voxel, and this is insufficient because it has been shown that a significant portion of WM voxels contain two or more fiber bundles with different orientations (Behrens et al., 2007, Jeurissen et al., 2012). One of the methods to address the “crossing fiber” problem is diffusion spectrum imaging (DSI) (Wedeen et al., 2005), which attempts to explicitly measure the diffusion PDF by sampling signals in the q-space, which is the Fourier reciprocal space of the PDF space. The q-space coordinate is defined by x_q = γgδ, and the q-space field-of-view (FOV) is practically controlled by the maximum b-value, defined by b = (2π|x_q|)²(Δ − δ/3), where γ is gyromagnetic ratio, g the diffusion-sensitive gradient, δ the gradient pulse duration and ∆ the pulse separation. A DSI dataset usually consists of a number of DW volume images, which spread on uniformly spaced grid points in the q-space. One can reconstruct the ensemble diffusion PDF from a DSI dataset by applying the inverse Fourier transform (FT) to the q-space signals of each voxel. Therefore, the DSI is a 6D dataset (Wedeen et al., 2005, Wedeen et al., 2008): 3D for the image space, which defines the anatomical position and the other 3D for the q-space, which encodes the diffusion displacements. To estimate the fiber orientations of a voxel from a DSI dataset, the usual approach is to construct the orientation distribution function (ODF) (Tuch et al., 2003) from the PDF (Wedeen et al., 2005, Wedeen et al., 2008), and the local maxima of the ODF are assumed to coincide with the fiber orientations within the considered voxel.

While DW-MRI is an essential technique for characterizing WM architectures, how to perform a spatial transformation for the DW-MRI datasets becomes an important issue for group analyses. Because of the fact that DW-MRI datasets are images of vectors (Cao et al., 2005), tensors (Basser et al., 1994), ODFs (Tuch et al., 2003), or PDFs (Wedeen et al., 2005), it is not readily applicable to employ the transformation methods designed for 3D scalar images for the DW-MRI datasets. A complete registration of the DW-MRI datasets usually consists of two transformations; one is to spatially align the 3D anatomical shapes, and the other is to align the local diffusivity profiles. While the former is similar to the classical transformation of 3D structures, the latter usually involves reorientation.

Alexander et al. (2001) proposed a transformation method for DTI datasets that consists of two steps: components of the DTI datasets are first spatially registered using the linear affine transformation, and then, the diffusion tensor at each voxel is reoriented, along with the affine transformation matrix, according to the “preservation of principal direction” (PPD) rule, which assumes that the transformation would alter the orientations of the tensors but not the eigenvalues of the tensors. The PPD reorientation strategy was later employed in nonlinear transformation methods for DTI (Cao et al., 2006, Ceritoglu et al., 2009, Chiang et al., 2008, Park et al., 2003, Van Hecke et al., 2007, Van Hecke et al., 2008, Van Hecke et al., 2009), in which the local affine transformation matrices were estimated from the local Jacobian matrices of the deformation fields. To reduce the error in estimating the amount of DTI reorientation, Xu et al. (2003, 2008) proposed another approach, which statistically estimated the rotation matrices for reorienting the tensors basing on the Procrustean method. While the reorientation issue is also addressed for vector images through local Jacobian matrices (Cao et al., 2005), it is more challenging to reorient ODFs because the ODFs essentially do not have a single fiber orientation. A simpler way to reorient ODFs is to employ the finite strain (FS) method (Alexander et al., 2001), in which the ODFs are rigidly rotated via the rotation matrices estimated from the local Jacobian matrices (Geng et al., 2011); however, the FS method does not account for the shearing or scaling effects of transformations. Raffelt et al. (2011) proposed a more complicated approach, which approximates the ODF as the sum of a number of weighted point spread functions (PSF) and reorients each PSF according to the local affine transformation, and then combines the reoriented PSFs using their original PSF weights. Hong et al. (2009) proposed a method that assumes that the volume fraction of water molecules diffusing along a given direction is preserved after applying a transformation, and starting from this assumption, an analytical closed form formula could be derived for the reoriented ODFs (Du et al., 2012).

Along with the reorientation strategies, several approaches have been proposed for transforming the vector (Cao et al., 2005), the DTI (Alexander et al., 2001, Cao et al., 2006, Ceritoglu et al., 2009, Chiang et al., 2008, Park et al., 2003, Van Hecke et al., 2007, Van Hecke et al., 2008, Van Hecke et al., 2009, Xu et al., 2003, Xu et al., 2008, Yeo et al., 2009), and the ODF (Du et al., 2012, Geng et al., 2011, Hong et al., 2009, Raffelt et al., 2011, Yap et al., 2011, Yeh and Tseng, 2011) datasets, either applying the reorientation after the transformation maps are estimated (Alexander et al., 2001, Ceritoglu et al., 2009, Hong et al., 2009, Park et al., 2003, Raffelt et al., 2011, Van Hecke et al., 2007, Van Hecke et al., 2008, Van Hecke et al., 2009, Xu et al., 2003, Xu et al., 2008, Yeh and Tseng, 2011), i.e., the estimated transformation maps are independent of the reorientation methods, or imposing the reorientation in the cost function (Cao et al., 2005, Cao et al., 2006, Chiang et al., 2008, Du et al., 2012, Geng et al., 2011, Yap et al., 2011, Yeo et al., 2009), i.e., engaging the reorientation during the estimation of the transformation maps. The former approach is relatively simple because it employs a two-step scenario similar to the method adopted in Alexander et al. (2001), thus providing the opportunity of combining different transformation methods and reorientation strategies (Yeh and Tseng, 2011). However, it has been shown that the latter approach gives better registration results than the former (Yeo et al., 2009). One of the difficulties in using the later approach is the derivation of the exact gradients of the cost functions (Geng et al., 2011, Yeo et al., 2009), and the difficulty of determining the exact gradients has been successfully addressed under the framework of large deformation diffeomorphic metric mapping (LDDMM).

The LDDMM method is proposed by Beg et al. (2005) for registering 2D/3D structural images. Different from other registration approaches which usually estimate the transformation maps directly, the LDDMM method instead models the transformation process as a fluid flow and attempts to estimate the velocity fields, from which the generated transformations are diffeomorphic and the transformation path is a geodesic. The LDDMM framework has been extended to register various DW-MRI related data types, including multi-contrast images (LDDMM-MC) (Ceritoglu et al., 2009), vector images (LDDMM-vector) (Cao et al., 2005), DTI (LDDMM-DTI) (Cao et al., 2006) and ODF (LDDMM-ODF) (Du et al., 2012). All of the methods characterize their own data structures, which represent different diffusion information. The LDDMM-MC method estimates the transformation maps using the standard LDDMM algorithm, but the velocity fields are weighted by the contrast images. The conventional contrast images for DTI are the null images (i.e., the images without applying the diffusion-sensitive gradients) and the fractional anisotropy (FA) (Basser and Pierpaoli, 1996) maps. After applying the estimated transformation maps to each component of the DTI dataset, the tensor of each voxel is reoriented using the PPD algorithm (Ceritoglu et al., 2009). Therefore, the LDDMM-MC method takes only the reduced diffusion information (i.e., the FA), not the whole DTI dataset, into consideration when estimating the velocity fields. The LDDMM-vector method is more complicated because it explicitly estimates the velocity fields by accounting for the whole vector dataset (estimated from the principal eigenvector of the DTI dataset). In contrast to the LDDMM-MC method, in the LDDMM-vector method, the components of the vector dataset interact with each other when estimating the velocity fields. This is attributed to the definition of the action of the diffeomorphism on the vector dataset, which defines how the vectors are reoriented according to the local Jacobian matrices (Cao et al., 2005). The LDDMM-DTI method also employs the whole DTI dataset in the estimation procedure, and defines the action of the diffeomorphism on the DTI dataset by incorporating the PPD strategy for local tensor reorientation (Cao et al., 2006). The LDDMM-ODF method, by defining its own action of diffeomorphism on the ODF dataset, which is similar to that employed by Hong et al. (2009), could explicitly account for the whole ODF dataset when estimating the solutions (Du et al., 2012). All of the methods (except LDDMM-MC) impose the reorientation in the cost function, and the analytical forms of the gradients of the cost functions with respect to the velocity fields are successfully derived (Cao et al., 2005, Cao et al., 2006, Du et al., 2012), suggesting that these methods could effectively register their corresponding DW-MRI datasets, not only in terms of the global shapes but also the local diffusivity profiles.

To the best of our knowledge, neither the transformation nor the reorientation issues are studied for the DSI datasets. In this paper, we proposed a registration method for DSI datasets by making use of the fact that DSI is 6D (Wedeen et al., 2005, Wedeen et al., 2008) and generalized the 2D/3D LDDMM algorithm to the 6D case with some modifications made for DSI datasets. A fundamental difference of our proposed method from other registration methods for DW-MRI datasets is that we treat the reorientation issue for DSI datasets as a free deformation in the PDF space or in the q-space, rather than operations like rotation or affine transformation that is governed by the deformation fields. This choice is primarily motivated by the recent progresses of DW-MRI analysis methods as well as the unique properties of the LDDMM framework. Recent studies show that the q-space signals (especially those with high b-values) have the potential to provide more information about the WM structures than the conventional diffusion indices (e.g., FA of DTI or generalized fractional anisotropy (GFA) (Tuch, 2004) of ODF). For example, Raffelt et al. (2012) recently showed that the fiber ODF amplitude could be a measure of apparent fiber density (AFD), especially for the fiber ODFs reconstructed from DW volume images with b-value ≧ 3000 s/mm². They also demonstrate that AFD has the potential to distinguish differences between populations in terms of not only the image space location, but also the orientation along which the differences exist. Therefore, it would be beneficial to explore the behaviors of the q-space signals. Since we model the alignment of diffusivity profiles in the PDF space or in the q-space as a free deformation, the unique properties of the LDDMM framework provide the possibilities to explicitly extract the relevant information in the q-space. For example, because the estimated transformation path is a geodesic, it is possible to employ the amplitude of deformation (AOD) (Risser et al., 2011) as a measure of the amount of deformation for each q-space voxel. Another potential application is to make use of the fact that the initial momentum encodes the entire geodesic flow, such that not only the image space component but also the q-space component of the initial momenta of a group of DSI datasets could be analyzed on a linear space.