1 Introduction
Diffusionweighted MRI methods are promising tools for characterizing tissue microstructure. While diffusion tensor imaging (DTI) and high angular resolution diffusion imaging (HARDI) methods are widely used methods, they do not provide a complete description of the diffusion distribution. In order to more accurately reconstruct the ensemble average propagator (EAP), a thorough sampling of the
space sampling is needed, which requires multiple value diffusion weighted imaging (mDWI). The EAP estimation using mDWI better characterizes more complex neural fiber geometries and nonGaussian diffusion behavior when compared to single value techniques. Recently, new space imaging techniques, diffusion spectrum imaging (DSI) Wedeen et al. (2005) and hybrid diffusion imaging (HYDI) Wu and Alexander (2007) have been developed for estimating the EAP. HYDI is a mDWI technique that samples the diffusion signal along concentric spherical shells in thespace, with the number of encoding directions increased with each shell to increase the angular resolution with the level of diffusion weighting. Originally, HYDI employed the fast Fourier transform (FFT) to reconstruct the EAP. However, the recent advent of analytical EAP reconstruction schemes, which obtain closedform expressions of the EAP, obviate the use of the FFT in HYDI. One such technique successfully validated on HYDI datasets is Bessel Fourier orientation reconstruction (BFOR)
Hosseinbor et al. (2013). While mDWI techniques like HYDI have not been widely used, the new human connectome project Essen et al. (2013) and spinoff projects will likely significantly increase the application. However, there is a lack of fundamental image analysis tools for mDWI and EAP maps, such as registration and atlas generation, that can fully utilize their information.In the last decades, researchers have spent great efforts on developing registration algorithms to align diffusion tensors derived from DTI and orientation distribution functions (ODFs) derived from HARDI [e.g., Alexander et al. (2001); Raffelt et al. (2011); Du et al. (2012)]. However, registration algorithms directly based on DWIs are few. The direct alignment of DWIs in the space utilizes the full diffusion information, is independent of the choice of diffusion models and their reconstruction algorithms (e.g., tensor, ODF), and unifies the transformation to align the local diffusion profiles defined at each voxel of two brains Dhollander et al. (2010); Yap and Shen (2012); Zhang et al. (2012). Dhollander et al.Dhollander et al. (2010) developed an algorithm that transforms the diffusion signals on a single shell of the space and preserves anisotropic as well as isotropic volume fractions. Yap et.al Yap and Shen (2012) proposed to decompose the diffusion signals on a single shell of the space into a series of weighted diffusion basis functions, reorient these functions independently based on a local affine transformation, and then recompose the reoriented functions to obtain the final transformed diffusion signals. This approach provides the representation of the diffusion signal and also explicitly models the isotropic component of the diffusion signals to avoid undesirable artifacts during the local affine transformation. Zhang et al. Zhang et al. (2012) developed a diffeomorphic registration algorithm for aligning DW signals on a single shell of the space.
Only recently, Dhollander et al. Dhollander et al. (2011) aligned DWIs on multiple shells of the space by first estimating transformation using a multichannel diffeomorphic mapping algorithm, in which generalized fractional anisotrophy (GFA) images computed from each shell were used as mapping objects, and then applying the transformation to DWIs in each shell using the DWI reorientation method in Dhollander et al. (2010). This approach neglected possible influences of the DWI reorientation on the optimization of the spatial transformation. Hsu et al. Hsu et al. (2012) generalized the large deformation diffeomorphic metric image mapping algorithm Beg et al. (2005) to DWIs in multiple shells of the space and considered the image domain and space as the spatial domain where the diffeomorphic transformation is applied to. The authors claimed that the reorientation of DWIs is no longer needed as the transformation also incorporates the deformation due to the shape differences in the diffusion profiles in the space. It is a robust registration approach with the explicit consideration of the large deformation in both the image domain and the space. However, its computational complexity and memory requirement are high.
While limited research has been done for aligning the HYDI images, efforts on the white matter atlas from HYDI is even less. Only recently, Dhollander et al. Dhollander et al. (2011) employed their multichannel diffeomorphic matching algorithm for the atlas generation using HYDI datasets. To our best knowledge, there is no research on probablistic atlas genertion for HYDI.
In this paper, we propose a new large deformation diffeomorphic metric mapping (LDDMM) algorithm to align HYDI datasets, denoted as LDDMMHYDI, and then develop a Bayesian estimation framework for generating the brain atlas. In particular, we adopt the BFOR framework in representing the space signal Hosseinbor et al. (2013). Unlike the diffeomorphic mapping of mDWIs in Hsu et al. Hsu et al. (2012), the BFOR signal basis provides the representation of the space signal and thus reduces memory requirement. In addition, since the BFOR signal basis is orthonormal, the norm that quantifies the differences in the space signals can be easily computed as the sum of the squared differences in the BFOR expansion coefficients. In this work, we will show that the reorientation of the space signal due to spatial transformation can be easily defined on the BFOR signal basis. Unlike the work in Dhollander et al. (2011), we will incorporate the BFOR signal basis into the LDDMM framework and derive the gradient descent algorithm for solving the LDDMMHYDI variational problem with explicit orientation optimization. Using this registration approach, we will further estimate the brain white matter atlas from the space based on a Bayesian model. This probabilistic model is the extension of the previous Bayesian atlas estimation for scalarbased intensity images Ma et al. (2008). With the aids of the BFOR representation and reorientation of mDWIs introduced in this work, we show that it is feasible to adopt the previous Bayesian atlas estimation model for scalarvalued images Ma et al. (2008) to HYDI. As shown below, the main contributions of this paper are:

to seek large deformation for aligning HYDI datasets based on the BFOR representation of mDWI.

to derive the rotationbased reorientation of the space signal via the BFOR signal basis. This is equivalent to applying Wigner matrix to the BFOR expansion coefficients, where Wigner matrix can be easily constructed by the rotation matrix (see Section 3.1).

to derive the gradient descent algorithm for the LDDMMHYDI variational problem with the explicit orientation optimization. In particular, we provide a computationally efficient method for calculating the variation of Wigner matrix due to the small variation of the diffeomorphic transformation (see Section 3.4).

to show that the LDDMMHYDI gradient descent algorithm does not involve the calculation of the BFOR signal bases and hence avoids the discretization in the space.

to propose a Bayesian estimation model for the space signals represented via the BFOR signal basis and derive an expectationmaximization algorithm for solving it (see Section 4).
2 Review: BFOR Signal Basis
According to the work in Hosseinbor et al. (2013), the space diffusion signal, , can be represented as
(1) 
where and respectively denote the image domain and space. is the th BFOR signal basis with its corresponding coefficient, , at . is given as
(2) 
Here, is the root of the order spherical Bessel (SB) function of the first kind . is the radial distance in space at which the Bessel function goes to zero. are the modified real and symmetric spherical harmonics (SH) bases as given in Goh et al. (2009). is the Bessel function of the first kind. is the number of terms in the modified SH bases of truncation order , while is the truncation order of radial basis. Note that in the original publication on BFOR Hosseinbor et al. (2013), the BFOR basis was not normalized to unity, but we have rectified this in Eq. (2). We refer readers to Hosseinbor et al. (2013) for more details on the BFOR algorithm.
Using the fact that the BFOR signal basis is orthonormal, the norm of can be easily written as
(3) 
3 Large Deformation Diffeomorphic Metric Mapping for HYDI
3.1 RotationBased Reorientation of
We now discuss the reorientation of when rotation transformation is applied. We assume that the diffusion profile in each shell of the space remains in the same shell after the reorientation. However, its angular profile in each shell of the space is transformed according to the rotation transformation. Hence, we define
According to the BFOR representation of in Eq. (1), we thus have
This indicates that the rotation reorientation of mDWI is equivalent to applying the rotation transformation to the real spherical harmonics, . According to the work in Geng et al. (2011), the rotation of can be achieved by the rotation of their corresponding coefficients, yielding
(4) 
where is the th element of Wigner matrix constructed based on (see details in Geng et al. (2011)). We can see that the same Wigner matrix is applied to at a fixed . For the sake of simplicity, we rewrite Eq. (4) in the matrix form, i.e.,
where is a sparse matrix with diagonal blocks of .
is a vector that concatenates coefficients
in the order such that at a fixed , corresponds to . concatenates the BFOR signal basis.3.2 Diffeomorphic Group Action on
We define an action of diffeomorphisms on , which takes into consideration of the reorientation in the space as well as the transformation of the spatial volume in . Based on the rotation reorientation of in Eq. (4), for a given spatial location , the action of on can be defined as
where can be defined in a way similar to the finite strain scheme used in DTI registration Alexander et al. (2001). That is, , where is the Jacobian matrix of at . For the remainder of this paper, we denote this as
(5) 
where indicates as the composition of diffeomorphisms.
3.3 Large Deformation Diffeomorphic Metric Mapping for HYDIs
The previous sections equip us with an appropriate representation of HYDI mDWI and its diffeomorphic action. Now, we state a variational problem for mapping HYDIs from one subject to another. We define this problem in the “large deformation” setting of Grenander’s group action approach for modeling shapes, that is, HYDI volumes are modeled by assuming that they can be generated from one to another via flows of diffeomorphisms
, which are solutions of ordinary differential equations
starting from the identity map . They are therefore characterized by timedependent velocity vector fields . We define a metric distance between a HYDI volume of a subject and an atlas HYDI volume as the minimal length of curves in a shape space such that, at time , . Lengths of such curves are computed as the integrated norm of the vector field generating the transformation, where , where is a reproducing kernel Hilbert space with kernel and norm . To ensure solutions are diffeomorphic, must be a space of smooth vector fields. Using the duality isometry in Hilbert spaces, one can equivalently express the lengths in terms of , interpreted as momentum such that for each , , where we let denote the inner product between and , but also, with a slight abuse, the result of the natural pairing between and in cases where is singular (e.g., a measure). This identity is classically written as , where is referred to as the pullback operation on a vector measure, . Using the identity and the standard fact that energyminimizing curves coincide with constantspeed lengthminimizing curves, one can obtain the metric distance between the template and target volumes by minimizing such that at time . We associate this with the variational problem in the form of(6) 
where is a positive scalar. quantifies the difference between the deformed atlas and the subject . Based on Eq. (3) and (5), is expressed in the form of
(7) 
3.4 Gradient of with respect to
We now solve the optimization problem in Eq. (6) via a gradient descent method. The gradient of with respect to can be computed via studying a variation on such that the derivative of with respect to is expressed in function of . According to the general LDDMM framework derived in Du et al. (2011), we directly give the expression of the gradient of with respect to as
(8) 
where
(9) 
Eq. (9) can be solved backward given . is the partial derivative of with respect to .
In the following, we discuss the computation of . We consider a variation of as and denote the corresponding variation in as . Denote for the simplicity of notation. We have
(10)  
As the calculation of Term (A) is straightforward, we directly give its expression, i.e.,
(11) 
This term is similar to that in the scalar image mapping case. It seeks the optimal spatial transformation in the gradient direction of image weighted by the difference between the atlas and subject’s images.
The computation of Term (B) involves the differential of with respect to rotation matrix and the variation of with respect to the small variation of . Let’s first compute the derivative of with respect to rotation matrix . According to the work in Cetingul et al. (2012), the analytical form of this derivative can be solved using the Euler angle representation of but is relatively complex. Here, we consider Wigner matrix and the coefficients of the BFOR signal basis together, which leads to a simple numeric approach for computing the derivative of with respect to rotation matrix , i.e., . Assume , where
is a skewsymmetric matrix parameterized by
. From this construction, is the tangent vector at on the manifold of rotation matrices and is also a rotation matrix. Based on Taylor expansion, we have the first order approximation of asHence, we can compute as follows. Assume to be skewsymmetric matrices respectively constructed from . We have
(12) 
It is worth noting that this formulation significantly reduces the computational cost for . Since is independent of spatial location , , , and are only calculated once and applied to all .
We now compute the variation of with respect to the small variation of . This has been referred as exact finitestrain differential that was solved in Dorst (2005) and applied to the DTI tensorbased registration in Yeo et al. (2009). Here, we directly adopt the result from Yeo et al. (2009) and obtain
(13) 
where . denotes as the cross product of two vectors. denotes the th column of matrix . .
Given Eq. (12) and (13), we thus have
Term (B)  (14)  
where
(15) 
and is approximated as
where are the neighbors of in directions, respectively. is the distance of these neighbors to . Here, term (B) seeks the spatial transformation such that the local diffusion profiles of the atlas and subject’s HYDIs have to be aligned.
3.5 Numerical Implementation
We so far derive and its gradient in the continuous setting. In this section, we elaborate the numerical implementation of our algorithm under the discrete setting. Since HYDI DW signals were represented using the orthonormal BFOR signal bases, both the computation of in Eq. (6) and the gradient computation in Eq. (16) do not explicitly involve the calculation . Hence, we do not need to discretize the space. In the discretization of the image domain, we first represent the ambient space, , using a finite number of points on the image grid, . In this setting, we can assume to be the sum of Dirac measures, where is the momentum vector at and time . We use a conjugate gradient routine to perform the minimization of with respect to . We summarize steps required in each iteration during the minimization process below:

Use the forward Euler method to compute the trajectory based on the flow equation:
(17) 
Compute based on Eq. (16).

Solve in Eq. (9) using the backward Euler integration, where indices , with the initial condition .

Compute the gradient .

Evaluate when , where is the adaptive step size determined by a golden section search.
4 Bayesian HYDI Atlas Estimation
The above registration problem requires an atlas. In this section, we introduce a general framework for Bayesian HYDI atlas estimation. Given observed HYDI datasets for , we assume that each of them can be estimated through an unknown atlas and a diffeomorphic transformation such that
(18) 
The total variation of relative to is then denoted by . The goal here is to estimate the unknown atlas and the variation . To solve for the unknown atlas , we first introduce an ancillary “hyperatlas” , and assume that our atlas is generated from it via a diffeomorphic transformation of such that . We use the Bayesian strategy to estimate and from the set of observations by computing the maximum a posteriori (MAP) of . This can be achieved using the ExpectationMaximization algorithm by first computing the loglikelihood of the complete data () when are introduced as hidden variables. We denote this likelihood as . We consider that the paired information of individual observations, for , as independent and identically distributed. As a result, this loglikelihood can be written as
(19)  
where
is the shape prior (probability distribution) of the atlas given the hyperatlas,
. is the distribution of random diffeomorphisms given the estimated atlas (). is the conditional likelihood of the HYDI data given its corresponding hidden variable and the estimated atlas (). In the remainder of this section, we first adopt and introduced in (Ma et al., 2008; Qiu et al., 2010) under the framework of LDDMM and then describe how to calculate in §4.2 based on the BFOR representation of HYDI.4.1 The Shape Prior of the Atlas and the Distribution of Random Diffeomorphisms
In the framework of LDDMMHYDI, we adopt the previous work in Ma et al. (2008); Qiu et al. (2010) and briefly describe the construction of the shape prior (probability distribution) of the atlas, . Under the LDDMM framework, we can compute the prior via , i.e.,
(20) 
where is initial momentum defined in the coordinates of such that it uniquely determines diffeomorphic geodesic flows from to the estimated atlas. When remains fixed, the space of the initial momentum provides a linear representation of the nonlinear diffeomorphic shape space, , in which linear statistical analysis can be applied. Hence, assuming is random, we immediately obtain a stochastic model for diffeomorphic transformations of . More precisely, we follow the work in (Ma et al., 2008; Qiu et al., 2010) and assume to be a centered Gaussian random field (GRF) model. The distribution of is characterized by its covariance bilinear form, defined by
where are vector fields in the Hilbert space of with reproducing kernel .
We associate with . The “prior” of in this case is then
where is the normalizing Gaussian constant. This leads to formally define the “logprior” of to be
(21) 
where we ignore the normalizing constant term .
We can construct the distribution of random diffeomorphisms, , in the similar manner. We define via the corresponding initial momentum defined in the coordinates of . We also assume that is random, and therefore, we again obtain a stochastic model for diffeomorphic transformations of . is assumed to be a centered GRF model with its covariance as , where is the reproducing kernel of the smooth vector field in a Hilbert space . Hence, we can define the log distribution of random diffeomorphisms as
(22) 
where as before, we ignore the normalizing constant term .
4.2 The Conditional Likelihood of the HYDI Data
Given the representation of the diffusion signals in the space using BFOR bases, we construct the conditional likelihood of the HYDI data via the BFOR coeffcients. We assume that
has a multivariate Gaussian distribution with mean of
and covariance , where are the BFOR coefficients associated with .From the Gaussian assumption, we can thus write the conditional “loglikelihood” of given and as
(23)  
where as before, we ignore the normalizing Gaussian term.
4.3 ExpectationMaximization Algorithm
We have shown how to compute the loglikelihood shown in Eq. (19). In this section, we will show how we employ the ExpectationMaximization algorithm to estimate the atlas, , for , and . From the above discussion, we first rewrite the loglikelihood function of the complete data in Eq. (19) as
(24)  
where and can be computed based on Eq. (5).
The EStep. The Estep computes the expectation of the complete data loglikelihood given the previous atlas
and variance
. We denote this expectation as given in the equation below,(25)  
The MStep. The Mstep generates the new atlas by maximizing the function with respcet to and . The update equation is given as
(26)  
where we use the fact that the conditional expectation of is constant. We solve and by separating the procedure for updating using the current value of , and then optimizing using the updated value of .
We now derive how to update values of and from function in Eq. (26). It is straightforward to obtain by taking the derivative of with respect to and setting it to zero. Hence, we have
(27) 
For updating , let . Using the change of variables strategy, we have
Since the second item in the above equation is independent of and ,we have
Comments
There are no comments yet.