Fluid Registration of Diffusion Tensor Images Using Information Theory

Ming-Chang Chiang; Alex D Leow; Andrea D Klunder; Rebecca A Dutton; Marina Barysheva; Stephen E Rose; Katie L McMahon; Greig I de Zubicaray; Arthur W Toga; Paul M Thompson

doi:10.1109/TMI.2007.907326

. Author manuscript; available in PMC: 2009 Oct 29.

Published in final edited form as: IEEE Trans Med Imaging. 2008 Apr;27(4):442–456. doi: 10.1109/TMI.2007.907326

Fluid Registration of Diffusion Tensor Images Using Information Theory

Ming-Chang Chiang ¹, Alex D Leow ², Andrea D Klunder ³, Rebecca A Dutton ⁴, Marina Barysheva ⁵, Stephen E Rose ⁶, Katie L McMahon ⁷, Greig I de Zubicaray ⁸, Arthur W Toga ⁹, Paul M Thompson ^*

PMCID: PMC2770435 NIHMSID: NIHMS47057 PMID: 18390342

Abstract

We apply an information-theoretic cost metric, the symmetrized Kullback-Leibler (sKL) divergence, or J-divergence, to fluid registration of diffusion tensor images. The difference between diffusion tensors is quantified based on the sKL-divergence of their associated probability density functions (PDFs). Three-dimensional DTI data from 34 subjects were fluidly registered to an optimized target image. To allow large image deformations but preserve image topology, we regularized the flow with a large-deformation diffeomorphic mapping based on the kinematics of a Navier-Stokes fluid. A driving force was developed to minimize the J-divergence between the deforming source and target diffusion functions, while reorienting the flowing tensors to preserve fiber topography. In initial experiments, we showed that the sKL-divergence based on full diffusion PDFs is adaptable to higher-order diffusion models, such as high angular resolution diffusion imaging (HARDI). The sKL-divergence was sensitive to subtle differences between two diffusivity profiles, showing promise for nonlinear registration applications and multisubject statistical analysis of HARDI data.

Keywords: Diffusion tensor imaging (DTI), fluid registration, high angular resolution diffusion imaging (HARDI), Kullback-Leibler divergence

I. INTRODUCTION

Diffusion tensor imaging (DTI) is becoming increasingly popular for in vivo studies of fiber architecture in the brain (tractography), investigations of anatomical and functional connectivity, and clinical studies of white matter (WM) integrity [1], [2]. Because DTI is sensitive to subtle differences in WM fiber orientation and diffusion anisotropy, it provides a powerful new approach to reveal 3-D deficit patterns in population studies of neurological diseases, e.g., Alzheimer’s disease, multiple sclerosis, or schizophrenia [2].

Multisubject studies and group comparisons of DTI require a scheme to correlate or compare data across subjects, and to perform statistical analysis of DTI signals in a standard coordinate space. DTI data is inherently multidimensional (with effectively six independent parameters defining a positive—definite symmetric matrix at each voxel) and non-Euclidean (in that Euclidean interpolation for diffusion tensors may lead to loss of anisotropy, or tensor swelling [3]). DTI, therefore, presents unique challenges for image analysis. Extending work on the differential geometry of manifolds [4], researchers at INRIA in Sophia-Antipolis, France [3], [5] developed an extensive calculus for performing computations on tensors. They derived a log-Euclidean metric for comparing and averaging diffusion tensors on a Riemannian manifold by transforming each diffusion tensor into its matrix logarithm counterpart. The processing of diffusion tensors then becomes Euclidean in the logarithmic domain, and standard formulae can be applied to compute means and variances. Fletcher and Joshi [6] further developed the notion of principal geodesic analysis (PGA), which extends principal component analysis to manifolds, by using an affine—invariant metric on the space of symmetric positive—definite tensors. Work treating diffusion tensors as elements of a Lie group with a non-Euclidean metric is being pursued on many fronts, for DTI smoothing and anisotropic filtering [7], interpolation, segmentation and fiber tracking, and estimation of probabilistic connectivity matrices [8]–[11]. Khurd et al. [12], for example, use isometric mapping and manifold learning techniques (eigendecomposition of the distance matrix) to directly fit a manifold to the tensors, compute its dimensionality, and distinguish groups using Hotelling’s T²statistics. Taking a different conceptual approach, diffusion tensors, which are positive—definite, may instead be interpreted as first-order summaries of 3-D probability density functions (representing diffusion processes) and these distributions can also be processed using information theory [13]. Tsuda et al. [14] introduced the Kullback-Leibler (KL) divergence as the distance between two Gaussian probability density functions (PDFs), parameterizing the PDFs by using a positive—definite matrix or tensor as their covariance matrices. Wang and Vemuri [15], as well as Lenglet et al. [16], pioneered the use of the symmetrized KL-divergence in the context of DTI segmentation, in which the symmetrized KL-divergence served as a natural measure of dissimilarity among symmetric positive—definite (SPD) tensors, and used it in a region-based active contour model. The KL-divergence was part of a cost function for segmentation, which was minimized in a level-set evolution scheme. However, to our knowledge, the symmetrized KL-divergence has not to date been used for DTI registration. This is the focus of our efforts in this paper.

For nonlinear registration of DT images, Guimond et al. [17] proposed a multidimensional “Demons” method (building on work on demons by Thirion et al. [18]), which treats the DTI as multiple-channel data. Ruiz-Alzola et al. [19] proposed a multiresolution scheme in which areas with local detail in the tensor data were matched at successively refined resolution levels; in this approach, deformation fields were interpolated using a Kriging estimator. Alexander and Gee [20] compared the performance of different DTI matching approaches using various DT similarity measures in a multiresolution elastic registration scheme, pioneering the concept of tensor reorientation, which, given an affine transform A, correctly applies only the rotational component of A to the tensors during warping, preserving the underlying fiber connectivity. More recently, Zhang et al. [21] included a “finite strain” strategy for tensor reorientation [22] in the cost function, and approximated a nonlinear deformation using a piecewise affine algorithm. Cao et al. [23], [24] have worked on a problem that is perhaps most closely related to the topic of this paper. They registered two DT images by minimizing the Euclidean distance between the principal eigenvectors [23] or between diffusion tensors themselves (using the Frobenius norm) [24], with images deformed using diffeomorphic flows.

In this paper, we match diffusion tensors based on variationally minimizing the symmetrized KL-divergence between Gaussian PDFs whose covariance matrices are given by the diffusion tensors. The registration is computed in a fluid regularization scheme, in which diffeomorphic and inverse—consistent mappings are guaranteed. Moreover, we extend the symmetrized KL-divergence to model higher-order water diffusivity profiles based on the concept of spherical harmonics [25], [26]. Our work is most related to that presented by Cao et al. [23], [24].

The remainder of the paper is organized as follows. Section II explains the variational scheme to minimize the cost function based on the symmetrized KL-divergence, the extension of the symmetrized KL-divergence to higher-order diffusion models, the concept of reorientation of diffusivity profiles, and the theory of the inverse—consistent fluid model. Section III details the experimental procedures in which DT images are registered. Section IV provides the measures by which the registration performance is evaluated. Section V presents the results on the registration of DTI. We also compared the symmetrized KL-divergence to other cost functions in DTI, in synthetic two-fiber diffusivity models, and in real high angular resolution diffusion imaging (HARDI) data. Section VI compares our algorithm with other approaches for DTI registration.

II. REGISTRATION SCHEME FOR DIFFUSION TENSORS

A. DTI Registration Using the Symmetrized KL-Divergence

In a diffusion tensor—mangnetic resonance (DT-MR) image, the random diffusion of water molecules with respect to a specific image point × can be modeled, to the first order, using a 3-D Gaussian probability density function (PDF), whose spatial covariance can be approximated by a positive—definite second-order symmetric tensor T(x)[27]

p_{x} (r) = \frac{1}{\sqrt{{(2 π)}^{3} | 2 T (x) τ |}} exp (\frac{- r^{T} T^{- 1} (x) r}{4 τ})

(1)

where r is the displacement of water molecules relative to x, and τ ≥ 0 is the diffusion time for a molecule initially located at x at time τ = 0. The superscript T denotes transpose. Without loss of generality, we can set D(x) = 2T(x)τ to absorb the diffusion time. Equation (1) then becomes

p_{x} (r) = \frac{1}{\sqrt{{(2 π)}^{3} | D (x) |}} exp (\frac{- r^{T} D^{- 1} (x) r}{2}) .

(2)

Throughout this paper, we denote tensors and matrices by uppercase bold fonts, vectors by lowercase bold fonts, and scalars by nonbold fonts.

Given two continuous PDFs p(x)and q(x), the KL-divergence (also known as the information gain, or relative entropy) of q(x)from p(x) is defined as [15]

KL (p | | q) = \int p (x) log \frac{p (x)}{q (x)} d x .

(3)

Because this measure is not symmetric, instead we use the symmetrized KL-divergence, (also known as the J-divergence), as proposed in [15], which is given by

\begin{array}{l} sKL (p, q) & = \frac{1}{2} {KL (p | | q) + KL (q | | p)} \\ = \frac{1}{2} + {\int p (x) log \frac{p (x)}{q (x)} d x \\ + \int q (x) log \frac{q (x)}{p (x)} d x} . \end{array}

(4)

Let S and T denote tensor-valued images, with diffusion tensors S(x) and T(x), respectively, at each voxel. As we model S(x) and T(x) as the covariance matrix of the associated Gaussian PDFs, they were replaced by positive—definite counterparts [28], if any of their eigenvalues was negative due to noise [29].

Denoting by y=x − u the mapping of S at voxel y to T at voxel x with the deformation vector u, taken together with (2) and (4), the symmetrized KL-divergence of the Gaussian PDFs p_S,x−u(r) and p_T,x(r) has a closed form given by [15]

\begin{array}{l} sKL (p_{S, x - u,} p_{T, x) =} & \frac{1}{4} {t r (T (x) S^{- 1} (x - u)) \\ + t r (S (x - u) T^{- 1} (x)) - 6} \end{array}

(5)

where tr(·) represents the trace of a square matrix. sKL(p_{S, x−u}, p_T,x) has been proved to be affine—invariant [15]. Vemuri et al. [15] pioneered the use of this functional, and have used it extensively for DT segmentation, but not registration.

We adopted the variational framework proposed in [30] to derive a driving force for registering the tensor-valued image S(x) to T(x) by minimizing the symmetrized KL-divergence of the two images

\begin{array}{l} E (u, S, T) = & \frac{1}{4} \int_{Ω} {t r (T (x) S^{- 1} (x - u)) \\ + t r (S (x - u) T^{- 1} (x)) - 6} d x \end{array}

(6)

where Ω is a bounded domain in which the two images are defined. If w is a small perturbation of the deformation field u that registers the images, then taking the definition of the first variation of E(u,S,T), we have

\begin{array}{l} {\frac{\partial E (u + ε w, S, T)}{\partial ε} |}_{ε = 0} \\ = \frac{1}{4} \int_{Ω} {{\frac{\partial [t r (T (x) S^{- 1} (x - (u + ε w)))]}{\partial ε} |}_{ε = 0} \\ + {\frac{\partial [t r (S (x - (u + ε w)) T^{- 1} (x))]}{\partial ε} |}_{ε = 0}} d x \\ = \frac{1}{4} \int_{Ω} {[{(S^{- 1} (x - u) T (x) S^{- 1} (x - u) - T^{- 1} (x))}_{c}]}^{T} \\ \cdot \nabla S (x - u) w (x) d x \end{array}

(7)

where (X)_c denotes the pq × 1 column vector formed by concatenating the columns of a p × q tensor X, and ∇S(x-u) = [(∂(S(x))_c/∂x₁)(∂(S(x))_c/∂x₂)…(∂(Sx))_c/∂x_n)] taking the value at x − u. See Appendix A for details of this derivation. Therefore, the driving force F(x,u) can be defined as

\begin{array}{l} F (x, u) & = & \nabla_{u} E \\ = & \frac{1}{4} {{[{(S^{- 1} (x - u) T (x) S^{- 1} (x - u) - T^{- 1} (x))}_{c}]}^{T} \\ {\times \nabla S (x - u)}}^{T} \\ = & \frac{1}{4} {[\nabla S (x - u)]}^{T} \\ \times {(S^{- 1} (x - u) T (x) S^{- 1} (x - u) - T^{- 1} (x))}_{c} . \end{array}

(8)

B. Applicability of the Symmetrized KL-Divergence to Higher-Order Diffusion Models

The Gaussian PDF in (1) is simply a first-order approximation of the diffusion behavior of water molecules within a voxel, assuming that the water diffusion at each voxel can be fully described by the covariance matrix of the Gaussian PDF. This is not the case in regions where fiber tracts are crossing within a single voxel, and a higher-order approximation is required, as illustrated in Fig. 1. Below we show that the symmetrized KL-divergence is extensible beyond this first-order diffusion model, if higher-order spherical harmonic terms are available, e.g., using diffusion spectrum imaging (DSI) [31], q-space reconstruction methods [32], or high angular resolution diffusion imaging (HARDI) [33]. We give the definition of symmetrized KL-divergence for higher-order diffusion profiles and compare it to other standard cost functions in HARDI data. The registration framework for HARDI data using the symmetrized KL-divergence is very similar to that for DTI, which we will implement and validate more comprehensively in a separate paper.

Fig. 1 — (a) To a first-order approximation, water diffusion is modeled by a single Gaussian PDF. The covariance matrix or tensor, with eigenvalue—eigenvector pairs (λ_i, e_i), characterizes the diffusion properties, and can be visualized as an ellipsoid in 3-D space. (b) However, a single Gaussian model is unable to capture multidirectional water diffusion in brain regions where intravoxel fiber crossing is significant. In this case, the probability distribution for water diffusion, which may have a shape such as that displayed in (b), can be estimated from the MR signals, as in diffusion spectrum or q-ball imaging.

For higher-order diffusion profiles, the symmetrized KL-divergence can no longer be defined using a second-rank tensor model. Inspired by [34], we model the scalar diffusivity function (apparent diffusion coefficient) D as a probability density function. The diffusivity function D is given by the signal attenuation in a specific diffusion-encoded gradient direction g based on the Stejskal-Tanner equation [35]

S (g) = S_{0} exp (- b D (g))

(9)

where b is the diffusion weighting factor, and g = g(θ,ϕ) = [sinθcosϕ sinϕsinϕ cosθ]^T : θ and ϕ are the polar and azimuthal angles. The PDF is then defined by normalizing the integral of D over the spherical angle Ω to 1

\begin{array}{l} p (θ, ϕ) & = & D (θ, ϕ) / g t r (D) \\ g t r (D) & = & \int_{Ω} D (θ, ϕ) d Ω \end{array}

(10)

where gtr(D)is the generalized trace of D [34]. For two diffusivity functions D_p and D_q, similar to (4), we define the symmetrized KL-divergence on the unit sphere based on the corresponding PDFs p(θ,ϕ and q(θ,ϕ)

\begin{array}{l} sKL (p, q) & = \frac{1}{2} {\int_{Ω} p log (\frac{p}{q}) d Ω + \int_{Ω} q log (\frac{p}{q}) d Ω} \\ = \frac{1}{2} {\int_{Ω} p log (p) d Ω - \int_{Ω} p log (q) d Ω \\ + \int_{Ω} q log (q) d Ω - \int_{Ω} q log (p) d Ω} . \end{array}

(11)

Applying (10) to the integrals in (11), for example

\begin{array}{l} \int_{Ω} p log (q) d Ω & = & \int_{Ω} \frac{D_{p}}{g t r (D_{p})} [log (D_{q}) - log (g t r (D_{q}))] d_{Ω} \\ = & [\frac{1}{g t r (D_{p})} \int_{Ω} D_{p} log (D_{q}) d Ω] \\ - log (g t r (D_{q})) \end{array}

(12)

we have

\begin{array}{l} sKL (p, q) = & \frac{1}{2} {\frac{1}{g t r (D_{p})} [\int_{Ω} D_{p} \log (D_{p}) d Ω \\ - \int_{Ω} D_{p} \log (D_{q}) d Ω] \\ + \frac{1}{g t r (D_{q})} [\int_{Ω} D_{q} \log (D_{q}) d Ω \\ - \int_{Ω} D_{q} \log (D_{p}) d Ω]} . \end{array}

(13)

Direct estimation of sKL in (13) is computationally expensive, but is faster if we expand the diffusivity functions D(θ,ϕ)in terms of a spherical harmonic (SH) series. See Appendix B for derivations.

C. Reorientation of Diffusion Profiles

Unlike in the case of warping a scalar-valued image, adjustment of the direction of the diffusivity profiles is required when warping a tensor-valued image or HARDI, in order to keep the direction of diffusivity profiles oriented with the direction of the underlying fibers. For DTI, the tensor reorientation step has been well described in [22] and [36]. In this paper, tensor orientations were adjusted at each iteration of the fluid registration using the preservation of principal direction (PPD) method [22]. When image S is deformed to match image T with the mapping h : T → S, the eigenvectors of the tensor (ordered as e₁, e₂, and e₃) at the lattice point x of S are distorted according to the Jacobian matrix J = ∇h⁻¹(x). Using the PPD procedure, the rotation matrix by which e₁is rotated onto Je₁/|Je₁|and e₂ onto the plane spanned by Je₁/|Je₁| and Je₂/|Je₂|, was applied to the whole tensor.

In HARDI, however, a single diffusivity function may have multiple local maxima of diffusivity, with respect to spherical angle, that are computationally expensive to determine [37], and thus impractical to perform at every iteration during image warping. Here, we propose a fast algorithm to determine the principal direction of the diffusivity function, based on principal component analysis (PCA) of its shape [38]. At each direction sampled by the diffusion gradient, g(θ_i,ϕ_i), 0 ≤ i < n_s, where n_s is the number of gradient directions, we define a point with distance to the origin proportional to D(θ_i,ϕ_i), d_i = D(θ_i,ϕ_i) · g(θ_i,ϕ_i) = (d_i0, d_i1, d_i2), where the last term is the Cartesian coordinates of d_i, and the diffusivity function is assumed to be centered at the origin. When the diffusivity functions are sampled at radially symmetric angles, d_iat opposing locations cancel, and therefore the mean vector µ of the point set {d_i} is a zero vector. The covariance matrix ∑ of {d_i} is given by

\sum = \frac{1}{n_{s}} \sum_{i = o}^{n_{s} - 1} d_{i} d_{i}^{T} .

(14)

The principal direction of the diffusivity function is determined by the first eigenvector of ∑. The rotation matrix R, which adjusts the direction of the diffusivity function, is then obtained using the PPD procedure, and the diffusivity is sampled discretely at the new gradient directions $g (θ_{i}^{'}, ϕ_{i}^{'}) = R \cdot g (θ_{i}, ϕ_{i})$ .

It is advantageous to compute values of the reoriented diffusivity function in the original directions (θ_i,ϕ_i), rather than the new ones $(θ_{i}^{'}, ϕ_{i}^{'})$ for a couple of reasons. First, we do not need to keep track of new reoriented gradient directions at each iteration of image warping. Second, the symmetrized KL-divergence (or other cost functions) must be computed from diffusivity functions that are sampled at the same gradient directions for the target image (which is fixed) and the source (which is moving and in which the diffusivity functions are being reoriented). Because the gradient directions are sampled discretely, the reoriented diffusivity functions are constructed by “pushing” the values sampled at the original directions to the new directions, so the values of the reoriented function at the original directions are not known. These new values can be computed from the SH series, as the basis functions Y_lm are continuous and defined at all spherical angles. If $(Δ θ_{i}, Δ ϕ_{i}) = ({θ^{'}}_{i} - θ_{i}, {ϕ^{'}}_{i} - ϕ_{i})$ , then the values of the reoriented diffusivity functions D′ in the original directions are given by

\begin{array}{l} D^{'} (θ_{i}, ϕ_{i}) & = D (θ_{i} - Δ θ_{i}, ϕ_{i} - Δ ϕ_{i}) \\ = \sum_{l = 0, e v e n}^{\infty} \sum_{m = - 1}^{l} c_{l m} Y_{l m} (θ_{i} - Δ θ_{i}, ϕ_{i} - Δ ϕ_{i}) \end{array}

(15)

where c_lm are the SH coefficients of the original diffusivity function D [see (B2)].

D. Inverse—Consistent Fluid Regularization

We adopted the inverse—consistent regularization scheme proposed in [39], and extended it from an elastic to a fluid model to deal with large local deformations while maintaining one-to-one topology. We refer to [6], [40]—[42] for background on these fluid models, where essentially mappings regularized by linear elastic or Laplacian penalties are not guaranteed to be diffeomorphic (i.e., differentiable one-to-one mappings with smooth inverses). Because of this disadvantage of small-deformation models, smoothness can be guaranteed by moving to a fluid or other large-deformation models (e.g., LDDMM [43] or geodesic shooting [44]). Let h(x) = x-u₁(x) be the displacement field mapping image S to image T (in the forward direction), and h⁻¹(x) = x−u₂(x) be the mapping from T to S (in the backward direction). The fluid regularization scheme is then formulated by minimizing the following matching cost functional, regarding forward and backward velocities v₁(x) and v₂(x), respectively

E_{total} (S, T) = \underset{E_{1}}{\underset{︸}{E (u_{1}, S, T) + \int_{0}^{1} \int_{Ω} {‖ L v_{1} (x) ‖}^{2} d x d t}} + \underset{E_{2}}{\underset{︸}{E (u_{2}, T, S) + \int_{0}^{1} \int_{Ω} {‖ L v_{2} (x) ‖}^{2} d x d t}}

(16)

where L = a∇²+b∇∇·+cI; I is the identity linear operator. We set a = µ and b µ + λ, and c = 0, µ and λ, where and are the viscosity constants (c is an inertial term that vanishes for very low Reynolds number flow, see [40]); E(u₁,S,T) and E(u₂,T,S) are given in (6). This cost functional is invariant to the order of S and T, i.e., E_total(S,T) = E_total(T,S).

We adopted the “greedy descent” algorithm, which updates the mapping by perturbing the displacement fields from the previous time step. This type of approach does not update the whole time-dependent velocity field at each iteration, but accumulates velocity fields from the gradient at each iteration. The iterative process is an extension of that in [39] for elastic registration, and detailed in Appendix D.

III. EXPERIMENTAL METHODS

A. DTI Acquisition

DTI datasets from 34 elderly adults (14F/20M, age 73.6 ± 9.0 SD years) were analyzed. These images were acquired as part of a previous study where subject demographics are detailed [45]. DTI was performed with a 1.5 T Siemens Sonata MRI scanner using an optimized diffusion tensor sequence [46]. The imaging parameters were 43 axial slices, FOV = 23cm, TR/TE 6000/106 ms, 2.5 mm slice thickness with a 0.25-mm gap, with an acquisition matrix of 128 × 128 and 60 images acquired at each location consisting of 16 low (b = 0)and 44 high diffusion-weighted images in which the encoding gradient vectors are uniformly radially distributed in space (b = 1100s mm⁻²) using the electrostatic approach described elsewhere [46]. The reconstruction matrix was 256 × 256, resulting in an in-plane resolution of 0.89 × 0.89mm². The total scan time was 8 min. A fractional anisotropy (FA) map was constructed from each DT image using the equation shown at the bottom of the page, where λ₁, λ₂, λ₃ are eigenvalues [1].

B. Preprocessing of DTI

Each subject’s FA image was manually stripped of nonbrain tissues, yielding a binary brain-extraction mask (cerebellum included) for the DT images. The masked images were then coregistered with global scaling (affine nine-parameter transformation) to the upsampled FA image (dimension: 256 × 256 × 132 voxels, resolution: 0.89 × 0.89 × 0.89 mm³) from a randomly selected control subject. The resulting transformation parameters were used to rotationally reorient the tensor at each voxel [21], [22], and then affine align the DT images based on trilinear interpolation of the log-transformed tensors. As noted in [3], [5], [6], and in other Lie group work, transforming diffusion tensors using the matrix logarithm transformation and performing subsequent arithmetic operations in the logarithmic domain is beneficial as it makes the processing Euclidean. The interpolation of log-transformed tensors causes less “swelling,” or loss of anisotropy, than direct interpolation of the original tensors. We used this interpolation method in the following fluid registration experiments.

FA = \sqrt{({(λ_{1} - λ_{2})}^{2} + {(λ_{2} - λ_{3})}^{2} + {(λ_{1} - λ_{3})}^{2}) / (2 (λ_{1}^{2} + λ_{2}^{2} + λ_{3}^{2}))}

C. Target Selection

As in several other studies, e.g., in [47] and [48], we elected to register each subject’s data to a typical single-subject target image included in the data set. This choice was motivated by the fact that, unlike a multisubject average intensity atlas target, a single-subject target typically has comparable contrast-to-noise, and spatial frequency content to the other images. However, to reduce the bias toward the particular geometry of a single subject, we adopted the “best individual target” (BIT) method proposed by Kochunov et al. [49]. (We note that arguments can be made in favor of using other templates, such as minimal deformation templates [50], or geodesic means [48], and the best approach may depend on the application). For each FA image serving as the target, a target quality score (TQS) [49] was computed, defined as the average voxel-wise product of the mean and the standard deviation of the deformation fields of all other subjects fluidly mapped onto the target by maximizing the Jensen-Rényi divergence [51]. Although our proposed fluid registration algorithm, which minimizes the symmetrized KL-divergence between two DT images, can also be used for target selection here, we decided to use an established scalar method for this step. The subject with the lowest TQS was then defined as the best individual target.

D. Fluid Registration of Diffusion Tensor Images

The DT images were registered respectively by minimizing the symmetrized KL-divergence under inverse—consistent fluid regularization. The time step ε was set to 0.02. In all of the following experiments, µ and λ were set to 0.9 and 6.0 based on our earlier experiences fluidly registering scalar-valued MR images [51]. We have found that a relatively high setting for λ is useful to penalize extreme volume distortions in the deforming template, and to keep any such distortions spatially smooth. To facilitate convergence, we modified the hierarchical multiresolution scheme used in [52] to solve [in (D2)] with a 1/8 sampling of data along each axis (32 × 32 × 32 voxels) for 2000 iterations, 1/4 subsampling (64×64×64 voxels) for 1000 iterations, and 1/2 resolution (128×128×128 voxels) for 50 iterations. At each resolution level, initial displacement vectors were obtained from the previous level using a cubic B-spline interpolation [53]. Iterative registration was stopped (or proceeded to the next resolution level) when the value of the symmetrized KL-divergence (either forward or backward direction of registration) was no longer monotonically decreasing. Fig. 2 shows the symmetrized KL-divergence plotted against the iteration number. The values of the symmetrized KL-divergence are consistent in both the forward and backward directions, and decrease smoothly. The final displacement vector fields were interpolated to the full resolution of the image. The registration processes were executed in parallel by submitting one process to one compute node of a 306-node, dual-processor SUN Microsystems V20z cluster. Each compute node has a dual 64-bit 2.4 GHz AMD Opteron CPU. The average computation time is about 5 h by a single processor, excluding preprocessing steps such as target selection.

Fig. 2 — The sKL-divergence for both the forward and backward fluid registration mappings decreases with increased iterations. The abrupt increases in the sKL-divergence are due to upsampling of the DTI to a higher resolution level in the multiresolution scheme.

IV. EVALUATION OF THE REGISTRATION ALGORITHM

A. Evaluation of Registration Accuracy

There is no widely-accepted gold standard for validation of nonlinear registration algorithms, so we evaluated the performance of our fluid mapping method on the basis of two criteria: 1) the similarity between the source and the target image is maximized (in terms of data agreement and the regularity of the solution), and 2) the variability of a population of deformed source images is minimized, when a mapping is optimized. To assess the accuracy of the mapping, we defined several overlap measures that were not explicitly optimized in the cost function, such as the inverse consistency (IC) on a test sample and several measures of agreement for the eigenstructures of the DTI before and after alignment. Clearly, each of these represents different desirable characteristics in the mappings that are arguably necessary but not sufficient for good registration. We used the tensor overlap index (OVL) [54] as a measure for the first criterion and compared the dyadic coherence of principal eigenvectors [21], [54] and the variance of FA, before and after fluid registration, to evaluate the second criterion.

Measurement of Tensor Overlap

We adopted the measure of overlap between eigenvalue—eigenvector pairs proposed by Basser and Pajevic [54] to define the matching of two diffusion tensors

OVL = \frac{\sum_{i = 1}^{3} λ_{i} λ_{i}^{'} {(e_{i} \cdot e_{i}^{'})}^{2}}{\sum_{i = 1}^{3} λ_{i} λ_{i}^{'}}

(17)

where (λ_i,e_i) is the eigenvalue—eigenvector pairs. For each voxel, OVL = 0 indicates no overlap and 1 indicates complete overlap of the three principal axes of diffusion tensors. The average of the OVL (aOVL) across all the voxels in the registered DTI within a region of interest (ROI) was used as an overall measure of the tensor alignment. There is a risk of circularity in defining a data fidelity term to assess how well the tensors are registered, as tensor agreement is formulated operationally using the registration cost function. We selected the eigenvalue—eigenvector metric as it is not identical to the cost function being used to drive the registration (other possible choices might include the inner products of the principal eigenvectors, after normalization). As noted earlier, our registration cost function is not based on the eigenvectors directly, but rather on the symmetrized KL-divergence of the diffusion PDFs (here simplified using Gaussian covariance functions to represent each PDF). As such it is a slightly more general model, given that the PDF approach is readily extensible to nonparametric or higher-order representations of the diffusion process, and to cases of fiber-crossing where fiber topography is not well represented by the eigenvectors of a first-order model.

Dyadic Coherence of Principal Eigenvectors

The variability within a group of tensors can be approximated by measuring the dispersion (or coherence) of individual principal eigenvectors from their mean vector. (Clearly, this is one of several measures of variance that can be defined for tensor-valued data.) This is achieved by analyzing the shape of the mean dyadic tensor of principal eigenvectors, given by [21], [54]

〈 e_{1} e_{1}^{T} 〉 = \frac{1}{N} \sum_{i = 1}^{N} e_{1}^{i} e_{1}^{i T}

(18)

where $e_{1}^{i}$ is the principal eigenvector of subject i. Then the dyadic coherence κ is defined as

κ = 1 - \sqrt{\frac{β_{2} + β_{3}}{2 β_{1}}}

(19)

where β₁, β₂, and β₃ are eigenvalues of $〈 e_{1} e_{1}^{T} 〉$ in descending order. κ = 0 is the minimum, indicating that the principal eigenvectors are randomly scattered, and the mean dyadic tensor $〈 e_{1} e_{1}^{T} 〉$ represents a sphere (β₁ = β₂ = β₃). When the principal eigenvectors are perfectly aligned, reaches a maximum value of 1.

Variance of FA

We compared the differences in the intersubject variance of FA values between affine—registered and fluidly registered DTI (by affine—registered, we mean data aligned with only a nine-parameter linear spatial transformation, which adjusts for differences in overall brain scale). Although this measure does not invoke an independent source of ground truth, its use is based on the premise that any misregistration on nonhomologous fiber tracts will result in increases in the FA voxel-level variance in a group. As such it makes sense to see if this FA variance can be reduced using diffeomorphic flows that preserve the overall FA values in the image (as only local rotations are ever applied to the tensors—they are not scaled or sheared). These comparisons help distinguish specific brain regions where cross-subject mismatches in fiber geometry are reduced by nonlinear registration. At each voxel, an F-statistic $f = σ_{1}^{2} / σ_{2}^{2}$ was obtained based on the intersubject variance of FA under two different conditions (fluid registered versus affine-registered only), and was converted to a Student’s t-statistic with (n−2) degrees-of-freedom, for paired comparison of data processed in the two different ways [55]

t = \frac{(f - 1) \sqrt{n - 2}}{\sqrt{4 f (1 - r^{2})}} .

(20)

Here, n denotes number of subjects; r is the Pearson correlation coefficient. Values of t < 0 indicate that the variance of the numerator (fluid-registered) is less than that of the denominator (affine–registered). We used color-coded maps to visualize t-values obtained from the last equation. Overall significance, correcting for the multiple spatial comparisons implicit in computing an image of statistics, was assessed by pFDR methods [56] (see Appendix E for details) for strong control over the likelihood of false rejections of the null hypotheses with multiple comparisons.

B. IC Error of the Deformation Fields

If h_ST(x)and h_TS(x) are the forward and backward displacement vector fields computed for registering two images S and T, then the average and maximum IC errors within a ROI Ω are defined as [52]

\begin{array}{l} Average IC error & = & \frac{1}{| Ω |} \int_{Ω} ‖ h_{S T} [h_{T S} (x)] - x ‖ d x \\ Maximum IC error & = & max_{x \in Ω} ‖ h_{S T} [h_{T S} (x)] - x ‖ . \end{array}

(21)

Equation (21) was discretized for implementation, and trilinear interpolation was used when h_ST(x) did not fall on a lattice point. This measure was used, following [52] to assess the internal consistency of the deformation fields, assuming that an ideal registration algorithm should set up correspondences that do not depend on the order of the images.

V. RESULTS

A. Inverse–Consistent Fluid Registration for DTI

Fig. 3 displays the forward (source to target) and backward (target to source) deformations of DTI based on the inverse—consistent fluid algorithm. The Jacobian determinant of the forward and backward deformation fields was confirmed to be greater than 0.6 everywhere in all subjects examined, confirming that the mappings were indeed diffeomorphic. The Jacobian tensor field, along with its determinant and rotation angle, is visualized on Fig. 4. Fiber tracts surrounding the lateral and third ventricles, e.g., the corpus callosum and cingulum, experienced the greatest magnitude of deformation. This is also shown by the tensor glyphs in Fig. 5, in which the diffusion tensors in the source image are displaced and reoriented to map the target. The average registered DTI (source to target) in Fig. 6 shows increased image sharpness (assessed visually) after fluid registration following affine registration, relative to the sharpness obtained after using affine registration only. Moreover, the variance of FA across the subjects is reduced after fluid registration (pFDR = 0.0061), especially in major WM fiber tracts with high diffusion anisotropy. These suggest that fluid registration further reduces the intersubject morphological variability. Fig. 7 depicts how registration performance depends on the fluid viscosity constants µ and λ by measuring the aOVL between the deformed source and the target DTI and the average IC error of the deformation fields. A high setting of µ and λ (e.g., µ = 0.9 and λ = 6.0) is favorable as it results in high aOVL and low IC error. Fig. 8 shows that our algorithm on average achieves subvoxel error in terms of IC when applied to noisy images and allowing large deformations (average IC error = 0.09±0.02 voxel; maximum IC error = 1.24±0.39 vowel). The IC error maps indicate that greater errors occur in the regions where greater magnitudes of deformation are required to map the target, such as the periventricular areas (see Fig. 3 and Fig. 4). Nevertheless, IC error is small (typically less than a voxel) in these areas too. At each voxel, the average IC error increases with the magnitude of the deformation vector [Fig. 8(b) and (c)]. This implies that in some cases there may be a trade-off between IC and registration accuracy. For instance, if the registration does nothing (i.e., the displacement field is zero everywhere) then the maximum and average IC errors will be zero too. Therefore, while low IC error is evidence of correct implementation of a registration strategy, it does not necessarily indicate high registration accuracy.

Fig. 4 — Visualization of the Jacobian tensor of the registration deformation field (for the forward direction of fluid registration) in one subject (upper row), and the average tensor across all subjects (lower row). Because the Jacobian tensor is not generally symmetric, we computed the corresponding strain tensor, defined as ∑ = J_T J = RS² R^T. The average tensor was obtained by averaging the strain tensors of all subjects in the log-Euclidean space, and exponentiating the result. (a) The determinant of the Jacobian tensor is det(Σ^1/2) = det(S). (b) The rotation angle was estimated by finding the smallest angle between the principal eigenvector of Σ and the x-, y-, and z- axes. If the deformation does not contain any local rotations, the principal eigenvector of Σ will be aligned with the x-, y-, and z- axes. The principal eigenvector of Σ is poorly defined when the Jacobian tensor is approximately isotropic, so we displayed the rotation angles only in regions where the Jacobian tensor was somewhat anisotropic, which we defined as when the fractional anisotropy of the Jacobian (strain) tensor > 0.3. (Note that the fractional anisotropy of the Jacobian tensor is different from the fractional anisotropy of the diffusion tensors). (c) To better visualize the Jacobian (strain) tensor, we modeled the tensor as an ellipsoid, using the visualization software “BrainSuite” (http://www.loni.ucla.edu/Software/) [62]. The RGB color coding is the same as in Fig. 3. Compared with the Jacobian determinant and rotational matrix maps from the single subject, the determinant of the *average* Jacobian (strain) tensor is closer to 1, the rotation angle is smaller, and the average tensor is more isotropic. Recent work has used this average strain map as a criterion to optimize a population average brain imaging template [63], and the accompanying average strain maps can be used, in general, as an index of how much a registration target deviates from the population.

Fig. 5 — Ellipsoidal tensor glyphs visualize fluid registration in tensor fields. The color of a tensor ellipsoid was determined by projections of the principal eigenvector onto the medial—lateral (coded in red), anterior—posterior (coded in green), and superior—inferior (coded in blue) axes of the anatomical reference frame, such that the color of the ellipsoid is independent of the view. We show tensor glyphs overlaid on the corresponding FA image, using the DTI visualization software, Camino, developed at University College, London [64] (left column), and also the ones that are fitted as colored ellipsoids (right column), as the former ones made using Camino present more information by overlaying tensor glyphs on the corresponding FA images, while the latter ones make it easier to see consistent patterns of fiber directions. The source and the target DT images are the same as in Fig. 3. Close-up views of the posterior corpus callosum for source, deformed source, and target DT images demonstrate that the tensor ellipsoids are displaced during the shape deformation.

Fig. 6 — (a) Average DTI obtained from registering all the subjects (the subject who served as the target was not included) to the target DTI displayed in Fig. 3, based on affine alignment (nine-parameter) only (top row), or affine followed by fluid registration (bottom row). Increased sharpness of the average DTI after fluid registration suggests that the fluid algorithm helps to eliminate local nonlinear morphological variability across subjects. (b) Color-coded t-maps obtained from comparing the variance of FA before and after fluid registration. Although a t-value greater than 4 already indicates a highly significant reduction in variance (this corresponds to light-blue areas, with P < 0.0005 in one-tailed paired comparisons; the degrees-of-freedom = number of subjects −2 = 31, and the subject serving as the target is not included), the color bar of the t-values is exaggerated to show that fluid registration reduces the variance of FA to a greater degree in the periventricular areas.

Fig. 7 — Average IC error and the aOVL at different settings of fluid viscosity constants µ and λ. The performance of the fluid algorithm is stable, in terms of high aOVL and low IC error, when µ is greater than 0.5.

Fig. 8 — (a) Color-coded maps are intended to display the spatial distribution of the IC error (upper panel: transverse view; lower panel: sagittal view). The greatest IC error is found in the splenium of the corpus callosum where a large warping of the underlying anatomy is required to map the target (see Fig. 3). This suggests, as might be naturally expected, that the IC error increases with the magnitude of the required deformation. This trend is more clearly demonstrated in the color-coded maps for the correlation coefficient r between the IC error and the magnitude of deformation in (b), and when the average IC error is plotted against the magnitude of deformation in (c)

B. Evaluation of Registration Performance

The OVL maps in Fig. 9 show that the major fiber structures are well matched in the cerebrum (e.g., the corpus callosum, cingulum, and internal capsule), as is the WM of the brainstem where the values of OVL are close to 1. The aOVL measure of correspondence in the DTI data increases after fluid registration based on minimization of the symmetrized KL-divergence (aOVL before/after fluid registration =0.448 ± 0.015/0.494 ± 0.012, P < 0.0001; paired-t tests). Moreover, the OVL maps demonstrate that the improvements in tensor alignment between linear and nonlinear registration occur mostly in the WM regions (defined as FA > 0.3, see [57]), with aOVL (before/after fluid registration)=0.567±0.032/0.672±0.021 (P < 0.0001 based on paired-t tests).

Fig. 10 shows that fluid registration increases the dyadic coherence of the principal eigenvectors. These findings are complementary to the OVL findings in Fig. 9.

We further compared the registration accuracy based on symmetrized KL-divergence versus other two cost functions, the inner product of diffusivity functions in (B3) with and without the linear (l=0) term [21], [58], denoted by IP_L(inner product with the linear term) and IP_D (inner product with the linear term deleted).IP_D is designed to compare only the anisotropic part of the diffusivity functions [21]. Under the assumption of Gaussian diffusion, IP_L and IP_D have a closed form in terms of the diffusion tensors [21]. See Appendix F for the derivation of the driving force for fluid registration. Fig. 11 shows that the aOVL measure of correspondence in DTI data after fluid registration is higher using symmetrized KL-divergence as the cost functional than using IP_D (symmetrized KL-divergence: 0.494±0.012,IP_D: 0.489 ± 0.012, P < 0.0001; paired-t tests). This difference is greater in the WM regions (defined as voxels where FA > 0.3), with aOVL (symmetrized KL-divergence/IP_D)=0.672±0.021/0.658±0.022 (P < 0.0001 based on paired-t tests). For the registration experiments using IP_L as the cost function, the gradient field was not robust enough to keep IP_L constantly increasing.

Fig. 11 — Comparisons of tensor overlap after fluid registration between using symmetrized KL-divergence and using IP_D as the cost function. (a) The bar graphs show that the aOVL is higher using symmetrized KL-divergence than IP_D (0.494 ± 0.012 versus 0.489 ± 0.012,P < 0.0001), especially in the WM regions, defined as regions where FA > 0.3 (0.672 ± 0.021 versus 0.658 ± 0.022,P < 0.0001). (b) This difference in the average OVL (averaged across all voxels with a given FA value) is plotted against the values of FA.

C. Application of the Symmetrized KL-Divergence to Higher-Order Diffusivity Models

Synthetic examples

We constructed a two-fiber diffusivity function using two orthogonal Gaussian tensors,T₀ and T₁, with typical eigenvalues for WM fibers [25], [26]. We set T₀ = diag(200,200,1700) × 10⁻⁶ (mm²/s); T₁ was obtained by rotating T₀ 90° around the y-axis. We also generated an isotropic gray-matter (GM) diffusivity function by setting T_iso = diag(700,700,700) × 10⁻⁶(mm²/s). Then the diffusivity function D is given by

D = - \frac{1}{b} \sum_{i = 0}^{n - 1} p_{i} exp (- b g^{T} T_{i} g)

(22)

where n is the number of fibers. p₀ = p₁ = 0.5 for two-fiber structures and p₀ = 1 for the isotropic diffusivity function. The gradient direction vector g consisted of 162 sampled points determined using an electrostatic approach [46]. Different b values (500, 1500, 3000) s/mm² were tested. The order of the SH series (l_m)was set to 8. We compared the symmetrized KL-divergence with the inner product of diffusivity functions, IP_L (with the linear term) and IP_D (without the linear term), on synthetic examples that were noise-free, or with Rician noise added to MR signals S such that the signal-to-noise ratio (SNR) was 35 or 10 [59]. Two identical two-fiber synthetic diffusivity functions served as the source and the target objects, with the source object rotating from 0 (complete overlap) to 90° . The three metrics were normalized for comparisons (see Fig. 12 for details).

Fig. 12 shows that at different noise levels, the symmetrized KL-divergence detects angular discrepancies in diffusivity functions more sensitively than IP_L and IP_D, in low b-settings. The symmetrized KL-divergence is also still comparable in performance with IP_D, at a high b-value. The performance of the symmetrized KL-divergence is therefore relatively stable under various diffusion weightings, and is applicable under ordinary MRI/DTI acquisition settings, though a high b-value can better detect higher-order angular structures in WM fibers at the cost of a decreased SNR [33].

HARDI Data

We further compared the three cost functions, which were summed over all voxels, in the 3-D HARDI data. The HARDI data was acquired from a healthy 22-year-old man imaged as part of a research study of twins on a 4 T Bruker Med-spec MRI scanner using an optimized diffusion tensor sequence [46]. Imaging parameters were: 21 axial slices (5 mm thick), FOV = 24 cm, TR/TE 6090/104.5 ms, 0.5-mm gap, with a 128 × 100 acquisition matrix and 30 images acquired at each location: 3 low (b = 0) and 27 high diffusion-weighted images in which the encoding gradient vectors were uniformly radially distributed in space (b = 1100s/mm²) using the electrostatic approach in [46]. The reconstruction matrix was 128 × 128, yielding a 1.875 × 1.875 mm² in-plane resolution. The total scan time was 3.09 min. We set l_m = 4 for the spherical harmonic analysis. The HARDI scan was compared with the same scan but rotated by angle φ(φ starts from 0° , which means complete overlap, up to ±20°), with the symmetrized KL-divergence and the inner products (IP_L and IP_D) computed at every 2°.

Fig. 13 shows that the symmetrized KL-divergence has a very sharp gradient near the optimal solution, and is sensitive enough to detect a 2° deviation in the images. The symmetrized KL-divergence is, therefore, a good candidate cost function for HARDI registration.

VI. DISCUSSION

In this paper, we proposed an inverse—consistent fluid registration algorithm, based on minimizing the symmetrized KL-divergence of the Gaussian PDFs between two DT-MR images, in which the diffusion tensor at each voxel was modeled as a spatial covariance matrix. Unlike the traditional tensor distance based on the Frobenius norm, this cost function based on the symmetrized KL-divergence (5) is affine invariant [15]. This invariance is helpful as it makes the computed correspondences between images robust to errors in estimating the linear (global alignment) term. Fiber structures between source and target DT images were well matched, especially in highly anisotropic areas, where the DTI signals convey the greatest information. Relative to approaches that aim to align the eigenvectors of the local diffusion process, we show that our approach is more generalizable, as the symmetrized KL-divergence is extendable to model higher-order water diffusivity profiles, by defining it as a PDF on the unit sphere and expanding the PDF using spherical harmonics. The symmetrized KL-divergence performs slightly but significantly better than IP_D (inner product of diffusivity functions with isotropic components deleted) as a cost functional for DTI registration in terms of tensor overlap. The OVL measure favors IP_D, but IP_D and OVL are somewhat correlated in gauging the alignment of eigenvalue—eigenvector pairs. We found that using IP_L, i.e., inner product with the linear term, did not provide a good gradient to maximize overlap between two DT images, which may indicate that introducing isotropic components of the diffusion tensors into a cost functional makes the cost functional less sensitive for accurate matching of DTI [21]. For HARDI, the symmetrized KL-divergence is more robust than standard inner product measures for detecting small rotational deviations between HARDI data at various diffusion weights and noise levels, making it an attractive measure for HARDI registration. In HARDI, we used the PPD procedure to adjust the direction of the diffusivity profiles, which is arguably more accurate than other methods (e.g., finite strain [21], [22]) given that it takes the originally estimated fiber directions into account. PCA determines one eigendirection, so it is appropriate both for diffusivity functions with a single global maximum, and for those with a dominant local maximum relative to other small local maxima. However, diffusivity functions may have multiple local maxima, e.g., in regions where fibers cross (e.g., the synthetic samples in Fig. 12), or when there is noise-free purely isotropic diffusion. In such cases, the principal direction estimated by PCA becomes arbitrary, and the simple PPD procedure may not be applicable in these regions. At this point, the finite strain method is advantageous as it extracts the rotational component from the Jacobian of the deformation fields, and the extracted rotational matrix can be used to reorient the diffusivity functions by resampling the SH series in (15). Nevertheless, the finite strain method discards some deformation components (e.g., shearing and stretching), and is accurate only when the image deformations are small, or when large deformations are decomposed into small ones (e.g., the piecewise affine algorithm in [21]).

Our algorithm is similar to that of Cao et al. [23], [24] for nonlinear warping of DTI data, using image deformations represented by flows that are guaranteed to be diffeomorphic (i.e., smooth mappings with smooth inverses). Cao et al.’s approaches have provided the key ingredient of adapting a scalar registration method to a tensor registration method. The cost function in our work is based on information theory, and we adopted an inverse—consistent fluid model to guarantee that our transformations are diffeomorphic. This work forms part of a general improvement in nonlinear registration methods that have recently been modified to allow 1) large diffeomorphic deformations (e.g., large deformation diffeomorphic metric mapping (LDDMM) [43] or geodesic shooting [44]), and 2) to achieve inverse–consistent mappings (i.e., maps where correspondences do not depend on the order of the two images being matched [39], [52]). The symmetrized KL-divergence, as a cost metric, can also be combined into the LDDMM scheme for DTI registration [24]. Nevertheless, our algorithm differs fundamentally from those by Cao et al. in that 1) the cost functional in our work is based on information theory rather than aligning principal fiber directions [23] or minimizing the Frobenius distance between two tensors [24]; 2) we showed that the symmetrized KL-divergence is extendable to any representation of the PDF that models higher-order water diffusion, whereas Cao et al.’s approaches assume that water diffusion in each voxel can be described by the dominant orientation or the three eigenvalue—eigenvector pairs of a rank-2 tensor; and 3) in Cao’s methods the tensor reorientation is explicitly optimized based on the PPD procedure, as eigenvectors of the tensors are part of the cost functional [21]. Because our approach can deal with generalized diffusion PDFs, it is likely to provide more reasonable assessments of signal matching in WM regions that may not have a single, well-defined eigenstructure. Ultimately, the best cost function should be the one that decreases error variance in the specific signal that is analyzed after the registration step, which differs in different applications of DTI (e.g., FA mapping versus statistics on fiber orientation parameters). In computing the symmetrized KL-divergence, it may help to incorporate tensor reorientation into the cost functional, e.g., by incorporating the dyadic expression for the diffusion tensor, $D = \sum_{i = 1}^{3} λ_{i} e_{i} e_{i}^{T}$ , into (6). This explicit reorientation, as used in [21] and [24], optimizes metrics comparing tensors as a whole, in which the gradients of the cost function incorporate the effect of tensor reorientation. Nevertheless, the symmetrized KL-divergence in (4) is a generalized metric and can be applied to higher-order diffusion imaging data (e.g., HARDI data), in which the closed form for explicit tensor reorientation may not exist. Moreover, in Fig. 2, we show that our implicit tensor reorientation scheme (where tensor reorientation is followed by each iteration of the fluid deformation) still ensures that the values of the KL-divergence for both the forward and backward deformations are smoothly decreasing, although we note that, compared with the explicit reorientation where the cost functional is minimized under the effect of reorientation, this implicit tensor reorientation scheme does not preclude the registration process from being trapped in local minima.

Finally, we should note some conceptual limitations of matching DTI. As with matching cortical gyri in standard structural images, there may not be a homology across subjects, between fibers or even major tracts at the gyrus by gyrus level, but some pragmatic normalization is necessary for multisubject mapping studies. Moreover, it might be problematic to match isotropic tensors, or those with very low anisotropy, e.g., in the gray matter and CSF, as good tract correspondence across subjects may not exist or may not be inferable from the available data [57]. In these cases, whether or not registration is successful depends on the goals and sensitivity of the analysis. Nevertheless, it is possible to quantify practical improvements in signal detection for clinically significant alterations in DTI parameters, without requiring a gold standard for matching fibers at the individual level.

Appendix

APPENDIX A

Derivation of (7)

If matrix Y is the function of matrix X, i.e., Y = Y(X), and A is a constant matrix, from linear algebra we have (d(tr(AY⁻¹))/dX) = −((Y^−TA^TY^−T)_c)^T(dY/dX), and (d(tr(YA))/dX = ((A^T)_c)^T(dY/dX). Using the above equations and the chain rule, we obtain

\begin{array}{l} {\frac{\partial [t r (T (x) S^{- 1} (x - (u + ε w)))]}{\partial ε} |}_{ε = 0} \\ = {\frac{d [t r (T S^{- 1} (z))]}{d z} |}_{z = x - (u + ε w)} ({\frac{\partial z (ε)}{\partial ε} |}_{ε = 0}) \\ = - {({(S^{- T} (x - u) T (x) S^{- T} (x - u))}_{c})}^{T} \\ \times \nabla S (x - u) (- w (x)) \\ = {({(S^{- 1} (x - >u) T (x) S^{- 1} (x - u))}_{c})}^{T} \\ \times \nabla S (x - u) (w (x) . \end{array}

(A1)

and

\begin{array}{l} {\frac{\partial [t r (S (x - (u + ε w)) T^{- 1} (x))]}{\partial ε} |}_{ε = 0} \\ = {\frac{d [t r (S (z) T^{- 1})]}{d z} |}_{z = x - (u + ε w)} ({\frac{\partial z (ε)}{\partial ε} |}_{ε = 0}) \\ = {({(T^{- T} (x))}_{c})}^{T} \nabla S (x - u) (- w (x)) \\ = - {({(T^{- 1} (x))}_{c})}^{T} \nabla S (x - u) (w (x) . \end{array}

(A2)

Here we use the property (S⁻¹)^T = (S^T)⁻¹ = S⁻¹ for a symmetric matrix S. Equation (7) is then derived by taking (A1) and (A2) together.

APPENDIX B

Computing sKL in terms of Spherical Harmonic Series

The diffusivity function D(θ, ϕ) can be expanded in terms of a spherical harmonic (SH) series [25], [26]

\begin{array}{l} D (θ, ϕ) & = & \sum_{t = 0}^{\infty} \sum_{m = - l}^{l} c_{l}^{m} Y_{l}^{m} (θ, ϕ), \\ Y_{l}^{m} (θ, ϕ) & = & \sqrt{\frac{(2 l + 1)}{4 π} \frac{(l - m)!}{(l + m)!}} P_{l}^{m} (\cos θ) e^{i m ϕ} . \end{array}

(B1)

Here, $P_{l}^{m} (cos θ)$ are the associated Legendre polynomials. D(θ, ϕ) is real and radially symmetric, so it is sufficient to adopt a real basis function set Y_lm while retaining the orthonormality of Y_lm [26]

\begin{array}{l} D (θ, ϕ) & = & \sum_{l = 0, even}^{\infty} \sum_{m = - l}^{l} c_{l m} Y_{l m} (θ, ϕ) \\ Y_{l m} & = & {\begin{array}{l} \frac{\sqrt{2}}{2} ({(- 1)}^{m} Y_{l}^{m} + Y_{l}^{- m}), \\ Y_{l}^{0}, \\ \frac{\sqrt{2 i}}{2} ({(- 1)}^{m + 1} Y_{l}^{m} + Y_{l}^{- m}), \end{array} & \begin{array}{l} if - l \leq m < 0 \\ if m = 0 \\ if 0 < m \leq l \end{array} \\ l = 0, 2, 4, \dots \end{array}

(B2)

As the are Y_lm orthonormal, the inner product of the real functions D₁ and D₂ can be expressed in terms of their SH coefficients ( $c_{l m}^{D_{1}}$ and $c_{l m}^{D_{2}}$ )

\int_{Ω} D_{1} (θ, ϕ) D {}_{2}{(θ, ϕ)} d Ω = \sum_{l = 0, even}^{\infty} \sum_{m = - l}^{l} c_{l m}^{D_{1}} c_{l m}^{D_{2}} .

(B3)

Moreover, it can be shown that (see Appendix C for derivations)

\int_{Ω} D (θ, ϕ) d Ω = 2 \sqrt{π} c_{00} .

(B4)

Given $D_{j} = \sum_{l = 0, even}^{\infty} \sum_{m = - l}^{l} c_{l m}^{j} Y_{l m}$ and $log (D_{j}) = \sum_{l = 0, even}^{\infty} \sum_{m = - l}^{l} d_{l m}^{j} Y_{l m}$ where j ∈ {p,q}, sKL in (13) can be expanded in terms of the SH series, as follows:

\begin{array}{l} sKL (p, q) = & \frac{1}{2} {\frac{1}{c_{00}^{p}} \sum_{l = 0, even}^{\infty} \sum_{m = - l}^{l} [c_{l m}^{p} d_{l m}^{p} - c_{l m}^{p} d_{l m}^{q}] \\ + \frac{1}{c_{00}^{q}} \sum_{l = 0, even}^{\infty} \sum_{m = - l}^{l} [c_{l m}^{q} d_{l m}^{q} - c_{l m}^{q} d_{l m}^{p}]} . \end{array}

(B5)

For numerical implementation, we use a truncated SH series with l ≤ l_m, where l_m is a positive even integer; the total length of the SH series is n_b=(l_m + 1)(l_m + 2)/2. We follow the least-squares method [25], [26] to solve (B2), yielding the SH coefficients

c = {(B^{T} B)}^{- 1} B^{T} d .

(B6)

Here, c = [c₀ c₁ … c_{n_b−1}]^T, and d = [D(θ₀,ϕ₀)D(θ₁,ϕ₁) … D(θ_{n_s−1}, ϕ_{n_s−1})]^T which represents the diffusivity function measured in n_s gradient directions, and B is the matrix of basis functions, with elements B_ij = Y_j(θ_i,ϕ_i. Here, we map (c_lm, Y_lm) in (B2) to (c_j, Y_j), using the relationship j=[(l²+l)/2] + m. Usually n_b ≪ n_s, so it is more cost-effective to compute sKL using the SH method (B5) rather than using (13) directly.

APPENDIX C

Derivation of (B4)

Given the orthonormality property of Y_lm

\int_{Ω} Y_{l}^{m} (θ, ϕ) {\bar{Y}}_{l^{'}}^{m^{'}} (θ, ϕ) d Ω = δ_{l l^{'}} δ m m^{'}

(C1)

where the overbar denotes complex conjugation and δ_ll′ is the Kronecker delta, we have

\begin{array}{l} \int_{Ω} Y_{l m} {\bar{Y}}_{00} d Ω \\ = {\begin{array}{l} \frac{\sqrt{2}}{2} \int_{Ω} ({(- 1)}^{m} Y_{l}^{m} + Y_{l}^{- m} \\ \times {\bar{Y}}_{0}^{0} d Ω = 0, & if - l \leq m < 0 \\ \int_{Ω} Y_{l}^{0} {\bar{Y}}_{0}^{0} d Ω = δ_{l 0}, & if m = 0 \\ \frac{\sqrt{2} i}{2} \int_{Ω} ({(- 1)}^{m + 1} Y_{l}^{m} + Y_{l}^{- m}) \\ \times {\bar{Y}}_{0}^{0} d Ω = 0, & if 0 < m \leq l . \end{array} \end{array}

(C2)

Here, we use the definition of Y_lm in (B2). Therefore

\begin{array}{l} \int_{Ω} D (θ, ϕ) d Ω & = & \int_{Ω} \sum_{l = 0, even}^{\infty} \sum_{m = - l}^{l} c_{l m} Y_{l m} (θ, ϕ) d Ω \\ = & \sum_{l = 0, even}^{\infty} \sum_{m = - l}^{l} c_{l m} \\ \times \int_{Ω} Y_{l m} (θ, ϕ) ({\bar{Y}}_{00} (θ, ϕ) \cdot 2 \sqrt{π}) d Ω \\ = & 2 \sqrt{π} \sum_{l = 0, even}^{\infty} c_{l 0} δ_{l 0} \\ = & 2 \sqrt{π} c_{00} . \end{array}

(C3)

Here, we note that $Y_{00} = Y_{0}^{0} = 1 / (2 \sqrt{π})$ [see (B1)].

APPENDIX D

Iterative Computations for Forward and Backward Displacement Fields

With the mapping initialized as the identical mapping at t = 0, the forward and the backward displacement fields at the time step n,h_n(x) and $h_{n}^{- 1} (x)$ , are given by

\begin{array}{l} h_{n} (x) & = & h_{n - 1} (x) + ε η_{1} (x) + ε η_{2} (x) \\ h_{n}^{- 1} (x) & = & h_{n - 1}^{- 1} (x) + ε ξ_{1} (x) + ε ξ_{2} (x) . \end{array}

(D1)

Here, ε is an infinitesimally small positive time step. η₁(x) and ξ₁(x) are the gradient directions that arise from the minimization of E₁ in (16) in the forward and backward mappings, respectively, and similarly η₂(x) and ξ₂(x) are defined in terms of E₂. $h_{0} (x) = h_{0}^{- 1} (x) = x$ are the initial mappings (and these are initially set to the identity, or zero displacement). E₁ and E₂were minimized by solving for the velocity fields v₁(x) and v₂(x) in the following Navier-Stokes equations [40], with the driving force F given in (8)

\begin{array}{l} \nabla E_{1, 2} & = & L v_{1, 2} + F (x, u_{1, 2}) \\ = & µ \nabla^{2} v_{1, 2} + (λ + μ) \nabla (\nabla \cdot v_{1, 2}) + F (x, u_{1, 2}) \\ = & 0 \end{array}

(D2)

corresponding to v₁(x) and v₂(x), respectively. This velocity field formulation, based on the kinematics of a compressible fluid with negligible inertia, was first developed by Christensen et al. [40]. By constructing the deformation as the integral of a time-varying velocity field, this approach is advantageous as it always guarantees diffeomorphic mappings (i.e., differentiable maps with differentiable inverses), even when large deformations are required; simpler formulations, such as those based on linear elasticity, thin-plate splines, or Laplacian penalties can lead to folding or tearing of the image when large deformations are necessary to match the images [60]. Equation (D2) was solved iteratively in time by convolving the force field with a filter kernel of size D × D × D. This filter is essentially a discretized approximation of the Green’s function of the linear operator L, derived using its eigenfunctions [61]. In the Eulerian reference frame, as used in [40], η₁(x) can be obtained from v₁(x) by

\begin{array}{l} η_{1} (x) & = & \frac{\partial h (x)}{d t} = - \frac{\partial u (x)}{\partial t} \\ = & - (v_{1} (x) - \sum_{i = 1}^{3} υ_{1, i} (x) \frac{\partial u (x)}{\partial x_{i}}) \\ = & - (\nabla h (x)) v_{1} (x) \end{array}

(D3)

where h = x − u is the displacement vector at time step n. Similarly, ξ₂(x)=−(∇h⁻¹(x))v₂(x). By concatenating h(x) and h⁻¹(x) using and (h⁻¹ o h)(x) = x and following the derivations in [39], we obtain

\begin{array}{l} η_{2} (x) & = & - (\nabla h (x)) ξ_{2} (h (x)) \\ = & - (\nabla h (x)) [- (\nabla h^{- 1} (h (x))) v_{2} (h (x))] \\ = & v_{2} (h (x)) \\ ξ_{1} (x) & = & - (\nabla h^{- 1} (x)) η_{1} (h^{- 1} (x)) \\ = & - (\nabla h^{- 1} (x)) [- (\nabla h (h^{- 1} (x))) v_{1} (h^{- 1} (x))] \\ = & v_{1} (h^{- 1} (x)) . \end{array}

(D4)

Here, we used the properties I = ∇ o h)(x) =[∇h⁻¹(h(x))]∇h(x) and I = ∇(h o h⁻¹)(x) = [∇h(h⁻¹(x))]∇h⁻¹(x).

APPENDIX E

The pFDR Method for Testing Overall Significance of Voxel Statistics

The pFDR method [56] estimates the probability that the null hypothesis is true, from the empirical distribution of observed p-values. Briefly, pFDR is the false discovery rate conditioning on the event that positive findings, rejecting the null hypothesis, have occurred, and is given by

pFDR (γ) = \frac{π_{0} Pr (P \leq γ | H = 0)}{Pr (P \leq γ)}

(E1)

where π₀ = Pr(H = 0) is the probability that the null hypothesis is true, and γ is the rejection threshold for the individual hypothesis, which was set to 0.01 in our experiments. We refer readers to [56] for the details of the estimation procedures to obtain pFDR for statistical maps. By convention, a statistical map with pFDR < 0.05, i.e., the false discovery rate is less than 5%, was considered to be significant.

APPENDIX F

Derivation of the Driving Force Based on IP_L and IP_D

Given two diffusivity functions, D_S and D_T, under the assumption of Gaussian diffusion, IP_L and IP_D of D_S and D_T can be extended in terms of their corresponding diffusion tensors and S, T, given by [21]

IP = \frac{8 π}{15} [t r (ST) + α \cdot t r (S) t r (T)] .

(F1)

IP = IP_L, if α = 1/2; IP = IP_D, if α = −1/3. Let h(x) = x − u(x) be the displacement field mapping image S to image T, the cost functional based on IP_L or IP_D is

\begin{array}{l} E (u, S, T) = & - \int_{Ω} \frac{8 π}{15} {t r (S (x - u) T (x)) \\ + α \cdot t r (S (x - u)) t r (T (x))} d x . \end{array}

(F2)

The negative sign before the cost functional integral is to transform a maximization problem into a minimization. Taking the gradient of E(u, S, T), as in (7) and (8) and (A1) and (A2), the driving force F(x, u) can be easily derived

F (x, u) = \frac{8 π}{15} {[\nabla S (x - u)]}^{T} {(T (x) + α \cdot t r (T (x)) \cdot I)}_{c}

(F3)

where I is the identity tensor.

Contributor Information

Ming-Chang Chiang, Laboratory of Neuro Imaging, Department of Neurology, UCLA School of Medicine, Los Angeles, CA 90095 USA..

Alex D. Leow, Laboratory of Neuro Imaging, Department of Neurology, UCLA School of Medicine, Los Angeles, CA 90095 USA.

Andrea D. Klunder, Laboratory of Neuro Imaging, Department of Neurology, UCLA School of Medicine, Los Angeles, CA 90095 USA.

Rebecca A. Dutton, Laboratory of Neuro Imaging, Department of Neurology, UCLA School of Medicine, Los Angeles, CA 90095 USA.

Marina Barysheva, Laboratory of Neuro Imaging, Department of Neurology, UCLA School of Medicine, Los Angeles, CA 90095 USA..

Stephen E. Rose, Centre for Magnetic Resonance, University of Queensland, 4072 Brisbane, Australia

Katie L. McMahon, Centre for Magnetic Resonance, University of Queensland, 4072 Brisbane, Australia

Greig I. de Zubicaray, Centre for Magnetic Resonance, University of Queensland, 4072 Brisbane, Australia

Arthur W. Toga, Laboratory of Neuro Imaging, Department of Neurology, UCLA School of Medicine, Los Angeles, CA 90095 USA.

REFERENCES

1.Mori S, Zhang J. Principles of diffusion tensor imaging and its applications to basic neuroscience research. Neuron. 2006;vol. 51:527–539. doi: 10.1016/j.neuron.2006.08.012. [DOI] [PubMed] [Google Scholar]
2.Catani M. Diffusion tensor magnetic resonance imaging tractography in cognitive disorders. Curr. Opin. Neurol. 2006;vol. 19:599–606. doi: 10.1097/01.wco.0000247610.44106.3f. [DOI] [PubMed] [Google Scholar]
3.Arsigny V, Fillard P, Pennec X, Ayache N. Fast and simple calculus on tensors in the log-Euclidean framework. Int. Conf. Med. Image Comput. Comput. Assist. Intervention (MICCAI); 2005. pp. 115–122. [DOI] [PubMed] [Google Scholar]
4.do Carmo MP. Riemannian Geometry. Boston, MA: Birkhauser; 1992. [Google Scholar]
5.Pennec X, Fillard P, Ayache N. A Riemannian framework for tensor computing. Int. J. Comput. Vis. 2006;vol. 66:41–66. [Google Scholar]
6.Fletcher PT, Joshi S. Comput. Vis. Math. Methods Med. Biomed, Image Anal.: ECCV 2004 Workshops CVAMIA MMBIA. Prague: Czech Republic; 2004. Principal geodesic analysis on symmetric spaces: Statistics of diffusion tensors; pp. 87–98. [Google Scholar]
7.Chung MK, Lazar M, Alexender AL, Lu Y, Davidson R. Probabilistic connectivity measure in diffusion tensor imaging via anisotropic kernel smoothing. Dept. Statist., Univ. Wisconsin; 2003. Tech. Rep. 1081. [Google Scholar]
8.Lenglet C, Deriche R, Faugeras O. Inferring white matter geometry from diffusion tensor MRI: Application to connectivity mapping. Eur. Conf. Comput. Vis. (ECCV); Prague. Czech Republic; 2004. pp. 127–140. [Google Scholar]
9.Chefd’hotel C, Tschumperlé D, Deriche R, Faugeras O. Constrained flows of matrix-valued functions: Application to diffusion tensor regularization. Int. Conf. Med. Image Comput. Comput. Assist. Intervention (MICCAI); 2002. pp. 251–265. [Google Scholar]
10.O’Donnell L, Haker S, Westin C-F. Med. Image Comput. Computer-Assisted Intervention (MICCAI) Tokyo, Japan: 2002. New approaches to estimation of white matter connectivity in diffusion tensor MRI: Elliptic PDEs and Geodesics in a tensor-warped space; pp. 459–466. [Google Scholar]
11.Fillard P, Arsigny V, Pennec X, Ayache N. Clinical DT-MRI estimation, smoothing and fiber tracking with log-Euclidean metrics. IEEE Int. Symp. Biomed. Imag.: From Nano to Macro (ISBI); Arlington, VA: 2006. [DOI] [PubMed] [Google Scholar]
12.Khurd P, Verma R, Davatzikos C. On characterizing and analyzing diffusion tensor images by learning their underlying manifold structure. 2006 Conf. Comput. Vis. Pattern Recognition Workshop; 2006. pp. 61–68. [Google Scholar]
13.Amari SI. Information geometry on hierarchy of probability distributions. Information Theory, IEEE Trans.; 2001. pp. 1701–1711. [Google Scholar]
14.Tsuda K, Akaho S, Asai K. The EM algorithm for kernel matrix completion with auxiliary data. J. Mach. Learn. Res. 2003;vol. 4:67–81. [Google Scholar]
15.Wang Z, Vemuri BC. DTI segmentation using an information theoretic tensor dissimilarity measure. IEEE Trans. Med. Imag.; 2005. Oct. pp. 1267–1277. [DOI] [PubMed] [Google Scholar]
16.Lenglet C, Rousson M, Deriche R. Segmentation of 3-D probability density fields by surface evolution: Application to diffusion MRI. the Int. Conf. Med. Image Comput. Comput. Assist Intervention (MICCAI); Saint-Malo, France. 2004. [Google Scholar]
17.Guimond A, Guttmann CRG, Warfield SK, Westin CF. Deformable registration of DT-MRI data based on transformation invariant tensor characteristics. Proc. IEEE Int. Symp. Biomed. Imag.; 2002. pp. 761–764. [Google Scholar]
18.Thirion J-P. Non-rigid matching using demons. IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognition (CVPR’96); 1996. pp. 245–251. [Google Scholar]
19.Ruiz-Alzola J, Westin CF, Warfield SK, Alberola C, Maier S, Kikinis R. Nonrigid registration of 3-D tensor medical data. Med. Image Anal. 2002;vol. 6:143–161. doi: 10.1016/s1361-8415(02)00055-5. [DOI] [PubMed] [Google Scholar]
20.Alexander DC, Gee JC. Elastic matching of diffusion tensor images. Comput. Vis. Image Understand. 2000;vol. 77:233–250. [Google Scholar]
21.Zhang H, Yushkevich PA, Alexander DC, Gee JC. Deformable registration of diffusion tensor MR images with explicit orientation optimization. Med. Image Anal. 2006;vol. 10:764–785. doi: 10.1016/j.media.2006.06.004. [DOI] [PubMed] [Google Scholar]
22.Alexander DC, Pierpaoli C, Basser PJ, Gee JC. Spatial transformations of diffusion tensor magnetic resonance. IEEE Trans. Med. Imag.; 2001. Nov. pp. 1131–1139. [DOI] [PubMed] [Google Scholar]
23.Cao Y, Miller MI, Winslow RL, Younes L. Large deformation diffeomorphic metric mapping of vector fields. IEEE Trans. Med. Imag.; 2005. Sep. pp. 1216–1230. [DOI] [PubMed] [Google Scholar]
24.Cao Y, Miller MI, Mori S, Winslow RL, Younes L. Diffeomorphic matching of diffusion tensor images. 2006 Conf. Comput. Vis. Pattern Recognition Workshop; 2006. pp. 67–74. [DOI] [PMC free article] [PubMed] [Google Scholar]
25.Alexander DC, Barker GJ, Arridge SR. Detection and modeling of non-Gaussian apparent diffusion coefficient profiles in human brain data. Magn. Reson. Med. 2002;vol. 48:331–340. doi: 10.1002/mrm.10209. [DOI] [PubMed] [Google Scholar]
26.Descoteaux M, Angelino E, Fitzgibbons S, Deriche R. Apparent diffusion coefficients from high angular resolution diffusion images: Estimation and applications INRIA. France: Sophia Antipolis; 2005. Res. Rep. 5681. [DOI] [PubMed] [Google Scholar]
27.Basser PJ, Mattiello J, LeBihan D. MR diffusion tensor spectroscopy and imaging. Biophys. J. 1994;vol. 66:259–267. doi: 10.1016/S0006-3495(94)80775-1. [DOI] [PMC free article] [PubMed] [Google Scholar]
28.Higham NJ. Computing a nearest symmetric positive semidefinite matrix. Linear Algebra Appl. 1988;vol. 103:103–118. [Google Scholar]
29.Parker GJ, Schnabel JA, Symms MR, Werring DJ, Barker GJ. Nonlinear smoothing for reduction of systematic and random errors in diffusion tensor imaging. J. Magn. Reson. Imag. 2000;vol. 11:702–710. doi: 10.1002/1522-2586(200006)11:6<702::aid-jmri18>3.0.co;2-a. [DOI] [PubMed] [Google Scholar]
30.Hermosillo G. Univ. Nice (INRIA-ROBOTVIS) France: Sophia Antipolis; 2002. Variational methods for multimodal image matching,” Ph.D. dissertation. [Google Scholar]
31.Tuch DS, Reese TG, Wiegell MR, Wedeen VJ. Diffusion MRI of complex neural architecture. Neuron. 2003;vol. 40:885–895. doi: 10.1016/s0896-6273(03)00758-x. [DOI] [PubMed] [Google Scholar]
32.Tuch DS. Q-ball imaging. Magn. Reson. Med. 2004;vol. 52:1358–1372. doi: 10.1002/mrm.20279. [DOI] [PubMed] [Google Scholar]
33.Tuch DS, Reese TG, Wiegell MR, Makris N, Belliveau JW, Wedeen VJ. High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity. Magn. Reson. Med. 2002;vol. 48:577–582. doi: 10.1002/mrm.10268. [DOI] [PubMed] [Google Scholar]
34.Ozarslan E, Vemuri BC, Mareci TH. Generalized scalar measures for diffusion MRI using trace, variance, and entropy. Magn. Reson. Med. 2005;vol. 53:866–876. doi: 10.1002/mrm.20411. [DOI] [PubMed] [Google Scholar]
35.Stejskal EO, Tanner JE. Spin diffusion measurements: Spin echoes in the presence of a time-dependent field gradient. J. Chem. Phys. 1965;vol. 42:288–292. [Google Scholar]
36.Xu D, Mori S, Shen D, van Zijl PC, Davatzikos C. Spatial normalization of diffusion tensor fields. Magn. Reson. Med. 2003;vol. 50:175–182. doi: 10.1002/mrm.10489. [DOI] [PubMed] [Google Scholar]
37.Alexander DC. Inf. Process. Med. Imag. Glenwood Springs, CO: 2005. Maximum entropy spherical deconvolution for diffusion MRI; pp. 76–87. [DOI] [PubMed] [Google Scholar]
38.Fukunaga K. Introduction to Statistical Pattern Recognition. 2nd ed. Boston, MA: Academic; 1990. [Google Scholar]
39.Leow A, Huang S, Geng A, Becker J, Davis S, Toga A, Thompson P. Inverse consistent mapping in 3-D deformable image registration: Its construction and statistical properties. the Inf. Proces. Med. Imag. (IPMI); Glenwood Springs, CO. 2005. [DOI] [PubMed] [Google Scholar]
40.Christensen GE, Rabbitt RD, Miller MI. Deformable templates using large deformation kinematics. IEEE Trans. Image Process.; 1996. Oct. pp. 1435–1447. [DOI] [PubMed] [Google Scholar]
41.Miller MI. Computational anatomy: Shape, growth, and atrophy comparison via diffeomorphisms. NeuroImage. 2004;vol. 23:S19–S33. doi: 10.1016/j.neuroimage.2004.07.021. [DOI] [PubMed] [Google Scholar]
42.Avants B, Gee JC. Geodesic estimation for large deformation anatomical shape averaging and interpolation. NeuroImage. 2004;vol. 23:S139–S150. doi: 10.1016/j.neuroimage.2004.07.010. [DOI] [PubMed] [Google Scholar]
43.Beg MF, Miller MI, Trouvé A, Younes L. Computing large deformation metric mappings via geodesic flows of diffeomorphisms. Int. J. Comput. Vis. 2005;vol. 61:139–157. [Google Scholar]
44.Vaillant M, Miller MI, Younes L, Trouve A. Statistics on diffeomorphisms via tangent space representations. NeuroImage. 2004;vol. 23:S161–S169. doi: 10.1016/j.neuroimage.2004.07.023. [DOI] [PubMed] [Google Scholar]
45.Rose SE, McMahon KL, Janke AL, O’Dowd B, de Zubicaray G, Strudwick MW, Chalk JB. Diffusion indices on magnetic resonance imaging and neuropsychological performance in amnestic mild cognitive impairment. J. Neurol. Neurosurg. Psychiatry. 2006;vol. 77:1122–1128. doi: 10.1136/jnnp.2005.074336. [DOI] [PMC free article] [PubMed] [Google Scholar]
46.Jones DK, Horsfield MA, Simmons A. Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magn. Reson. Med. 1999;vol. 42:515–525. [PubMed] [Google Scholar]
47.Studholme C, Cardenas V, Schuff N, Rosen H, Miller B, Weiner M. Detecting spatially consistent structural differences in Alzheimer’s and fronto temporal dementia using deformation morphometry. Int. Conf. Med. Image Comput. Comput. Assist. Intervention (MICCAI); 2001. pp. 41–48. [Google Scholar]
48.Joshi S, Davis B, Jomier M, Gerig G. Unbiased diffeomorphic atlas construction for computational anatomy. NeuroImage. 2004;vol. 23:S151–S160. doi: 10.1016/j.neuroimage.2004.07.068. [DOI] [PubMed] [Google Scholar]
49.Kochunov P, Lancaster JL, Thompson P, Woods R, Mazziotta J, Hardies J, Fox P. Regional spatial normalization: Toward an optimal target. J. Comput. Assist. Tomogr. 2001;vol. 25:805–816. doi: 10.1097/00004728-200109000-00023. [DOI] [PubMed] [Google Scholar]
50.Kochunov P, Lancaster J, Thompson P, Toga AW, Brewer P, Hardies J, Fox P. An optimized individual target brain in the Talairach coordinate system. NeuroImage. 2002;vol. 17:922–927. [PubMed] [Google Scholar]
51.Chiang M-C, Dutton RA, Hayashi KM, Lopez OL, Aizenstein HJ, Toga AW, Becker JT, Thompson PM. 3-D pattern of brain atrophy in HIV/AIDS visualized using tensor-based morphometry. NeuroImage. 2007;vol. 34:44–60. doi: 10.1016/j.neuroimage.2006.08.030. [DOI] [PMC free article] [PubMed] [Google Scholar]
52.Christensen GE, Johnson HJ. Invertibility and transitivity analysis for nonrigid image registration. J. Electron. Imag. 2003;vol. 12:106–117. [Google Scholar]
53.Thevenaz P, Ruttimann UE, Unser M. A pyramid approach to subpixel registration based on intensity. IEEE Trans. Image Process; 1998. pp. 27–41. [DOI] [PubMed] [Google Scholar]
54.Basser PJ, Pajevic S. Statistical artifacts in diffusion tensor MRI (DT-MRI) caused by background noise. Magn. Reson. Med. 2000;vol. 44:41–50. doi: 10.1002/1522-2594(200007)44:1<41::aid-mrm8>3.0.co;2-o. [DOI] [PubMed] [Google Scholar]
55.Sheskin D. Handbook of Parametric and Nonparametric Statistical Procedures. 3rd ed. Boca Raton, FL: Chapman Hall/CRC; 2004. [Google Scholar]
56.Storey JD. A direct approach to false discovery rates. J. R. Stat. Soc.: Series B (Stat. Methodol.) 2002;vol. 64:479–498. [Google Scholar]
57.Smith SM, Jenkinson M, Johansen-Berg H, Rueckert D, Nichols TE, Mackay CE, Watkins KE, Ciccarelli O, Cader MZ, Matthews PM, Behrens TE. Tract-based spatial statistics: Voxelwise analysis of multi-subject diffusion data. NeuroImage. 2006;vol. 31:1487–1505. doi: 10.1016/j.neuroimage.2006.02.024. [DOI] [PubMed] [Google Scholar]
58.Descoteaux M, Angelino E, Fitzgibbons S, Deriche R. A fast and robust ODF estimation algorithm in Q-ball imaging. the IEEE Int. Symp. Biomed. Imag.: From Nano to Macro (ISBI); Arlington, VA. 2006. [Google Scholar]
59.Sijbers J, den Dekker AJ, Van Audekerke J, Verhoye M, Van Dyck D. Estimation of the noise in magnitude MR images. Magn. Reson. Imag. 1998;vol. 16:87–90. doi: 10.1016/s0730-725x(97)00199-9. [DOI] [PubMed] [Google Scholar]
60.Christensen GE, Rabbitt RD, Miller MI, Joshi SC, Grenander U, Coogan TA, Essen DCV. Topological properties of smooth anatomic maps. the 14th Conf. Inf. Process. Med. Imag. Ile de Berder; France. 1995. [Google Scholar]
61.Bro-Nielsen M, Gramkow C. Fast fluid registration of medical images. Proc. 4th Int. Conf. Vis. Biomed. Comput.; 1996. pp. 267–276. [Google Scholar]
62.Shattuck DW, Leahy RM. BrainSuite: An automated cortical surface identification tool. Med. Image Anal. 2002;vol. 6:129–142. doi: 10.1016/s1361-8415(02)00054-3. [DOI] [PubMed] [Google Scholar]
63.Leporé N, Brun CA, Pennec X, Chou YY, Lopez OL, Aizenstein HJ, Becker JT, Toga AW, Thompson PM. Mean template for tensor-based morphometry using deformation tensors. Int. Conf. Med. Image Comput. Comput. Assist. Interv. (MICCAI); Brisbane, Australia: 2007. [DOI] [PubMed] [Google Scholar]
64.Cook PA, Bai Y, Nedjati-Gilani S, Seunarine KK, Hall MG, Parker GJ, Alexander DC. Camino: Open-source diffusion-MRI reconstruction and processing. 14th Sci. Meeting Int. Soc. Magn. Reson. Med.; Seattle, WA. 2006. p. 2759. [Google Scholar]

[R1] 1.Mori S, Zhang J. Principles of diffusion tensor imaging and its applications to basic neuroscience research. Neuron. 2006;vol. 51:527–539. doi: 10.1016/j.neuron.2006.08.012. [DOI] [PubMed] [Google Scholar]

[R2] 2.Catani M. Diffusion tensor magnetic resonance imaging tractography in cognitive disorders. Curr. Opin. Neurol. 2006;vol. 19:599–606. doi: 10.1097/01.wco.0000247610.44106.3f. [DOI] [PubMed] [Google Scholar]

[R3] 3.Arsigny V, Fillard P, Pennec X, Ayache N. Fast and simple calculus on tensors in the log-Euclidean framework. Int. Conf. Med. Image Comput. Comput. Assist. Intervention (MICCAI); 2005. pp. 115–122. [DOI] [PubMed] [Google Scholar]

[R4] 4.do Carmo MP. Riemannian Geometry. Boston, MA: Birkhauser; 1992. [Google Scholar]

[R5] 5.Pennec X, Fillard P, Ayache N. A Riemannian framework for tensor computing. Int. J. Comput. Vis. 2006;vol. 66:41–66. [Google Scholar]

[R6] 6.Fletcher PT, Joshi S. Comput. Vis. Math. Methods Med. Biomed, Image Anal.: ECCV 2004 Workshops CVAMIA MMBIA. Prague: Czech Republic; 2004. Principal geodesic analysis on symmetric spaces: Statistics of diffusion tensors; pp. 87–98. [Google Scholar]

[R7] 7.Chung MK, Lazar M, Alexender AL, Lu Y, Davidson R. Probabilistic connectivity measure in diffusion tensor imaging via anisotropic kernel smoothing. Dept. Statist., Univ. Wisconsin; 2003. Tech. Rep. 1081. [Google Scholar]

[R8] 8.Lenglet C, Deriche R, Faugeras O. Inferring white matter geometry from diffusion tensor MRI: Application to connectivity mapping. Eur. Conf. Comput. Vis. (ECCV); Prague. Czech Republic; 2004. pp. 127–140. [Google Scholar]

[R9] 9.Chefd’hotel C, Tschumperlé D, Deriche R, Faugeras O. Constrained flows of matrix-valued functions: Application to diffusion tensor regularization. Int. Conf. Med. Image Comput. Comput. Assist. Intervention (MICCAI); 2002. pp. 251–265. [Google Scholar]

[R10] 10.O’Donnell L, Haker S, Westin C-F. Med. Image Comput. Computer-Assisted Intervention (MICCAI) Tokyo, Japan: 2002. New approaches to estimation of white matter connectivity in diffusion tensor MRI: Elliptic PDEs and Geodesics in a tensor-warped space; pp. 459–466. [Google Scholar]

[R11] 11.Fillard P, Arsigny V, Pennec X, Ayache N. Clinical DT-MRI estimation, smoothing and fiber tracking with log-Euclidean metrics. IEEE Int. Symp. Biomed. Imag.: From Nano to Macro (ISBI); Arlington, VA: 2006. [DOI] [PubMed] [Google Scholar]

[R12] 12.Khurd P, Verma R, Davatzikos C. On characterizing and analyzing diffusion tensor images by learning their underlying manifold structure. 2006 Conf. Comput. Vis. Pattern Recognition Workshop; 2006. pp. 61–68. [Google Scholar]

[R13] 13.Amari SI. Information geometry on hierarchy of probability distributions. Information Theory, IEEE Trans.; 2001. pp. 1701–1711. [Google Scholar]

[R14] 14.Tsuda K, Akaho S, Asai K. The EM algorithm for kernel matrix completion with auxiliary data. J. Mach. Learn. Res. 2003;vol. 4:67–81. [Google Scholar]

[R15] 15.Wang Z, Vemuri BC. DTI segmentation using an information theoretic tensor dissimilarity measure. IEEE Trans. Med. Imag.; 2005. Oct. pp. 1267–1277. [DOI] [PubMed] [Google Scholar]

[R16] 16.Lenglet C, Rousson M, Deriche R. Segmentation of 3-D probability density fields by surface evolution: Application to diffusion MRI. the Int. Conf. Med. Image Comput. Comput. Assist Intervention (MICCAI); Saint-Malo, France. 2004. [Google Scholar]

[R17] 17.Guimond A, Guttmann CRG, Warfield SK, Westin CF. Deformable registration of DT-MRI data based on transformation invariant tensor characteristics. Proc. IEEE Int. Symp. Biomed. Imag.; 2002. pp. 761–764. [Google Scholar]

[R18] 18.Thirion J-P. Non-rigid matching using demons. IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognition (CVPR’96); 1996. pp. 245–251. [Google Scholar]

[R19] 19.Ruiz-Alzola J, Westin CF, Warfield SK, Alberola C, Maier S, Kikinis R. Nonrigid registration of 3-D tensor medical data. Med. Image Anal. 2002;vol. 6:143–161. doi: 10.1016/s1361-8415(02)00055-5. [DOI] [PubMed] [Google Scholar]

[R20] 20.Alexander DC, Gee JC. Elastic matching of diffusion tensor images. Comput. Vis. Image Understand. 2000;vol. 77:233–250. [Google Scholar]

[R21] 21.Zhang H, Yushkevich PA, Alexander DC, Gee JC. Deformable registration of diffusion tensor MR images with explicit orientation optimization. Med. Image Anal. 2006;vol. 10:764–785. doi: 10.1016/j.media.2006.06.004. [DOI] [PubMed] [Google Scholar]

[R22] 22.Alexander DC, Pierpaoli C, Basser PJ, Gee JC. Spatial transformations of diffusion tensor magnetic resonance. IEEE Trans. Med. Imag.; 2001. Nov. pp. 1131–1139. [DOI] [PubMed] [Google Scholar]

[R23] 23.Cao Y, Miller MI, Winslow RL, Younes L. Large deformation diffeomorphic metric mapping of vector fields. IEEE Trans. Med. Imag.; 2005. Sep. pp. 1216–1230. [DOI] [PubMed] [Google Scholar]

[R24] 24.Cao Y, Miller MI, Mori S, Winslow RL, Younes L. Diffeomorphic matching of diffusion tensor images. 2006 Conf. Comput. Vis. Pattern Recognition Workshop; 2006. pp. 67–74. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R25] 25.Alexander DC, Barker GJ, Arridge SR. Detection and modeling of non-Gaussian apparent diffusion coefficient profiles in human brain data. Magn. Reson. Med. 2002;vol. 48:331–340. doi: 10.1002/mrm.10209. [DOI] [PubMed] [Google Scholar]

[R26] 26.Descoteaux M, Angelino E, Fitzgibbons S, Deriche R. Apparent diffusion coefficients from high angular resolution diffusion images: Estimation and applications INRIA. France: Sophia Antipolis; 2005. Res. Rep. 5681. [DOI] [PubMed] [Google Scholar]

[R27] 27.Basser PJ, Mattiello J, LeBihan D. MR diffusion tensor spectroscopy and imaging. Biophys. J. 1994;vol. 66:259–267. doi: 10.1016/S0006-3495(94)80775-1. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R28] 28.Higham NJ. Computing a nearest symmetric positive semidefinite matrix. Linear Algebra Appl. 1988;vol. 103:103–118. [Google Scholar]

[R29] 29.Parker GJ, Schnabel JA, Symms MR, Werring DJ, Barker GJ. Nonlinear smoothing for reduction of systematic and random errors in diffusion tensor imaging. J. Magn. Reson. Imag. 2000;vol. 11:702–710. doi: 10.1002/1522-2586(200006)11:6<702::aid-jmri18>3.0.co;2-a. [DOI] [PubMed] [Google Scholar]

[R30] 30.Hermosillo G. Univ. Nice (INRIA-ROBOTVIS) France: Sophia Antipolis; 2002. Variational methods for multimodal image matching,” Ph.D. dissertation. [Google Scholar]

[R31] 31.Tuch DS, Reese TG, Wiegell MR, Wedeen VJ. Diffusion MRI of complex neural architecture. Neuron. 2003;vol. 40:885–895. doi: 10.1016/s0896-6273(03)00758-x. [DOI] [PubMed] [Google Scholar]

[R32] 32.Tuch DS. Q-ball imaging. Magn. Reson. Med. 2004;vol. 52:1358–1372. doi: 10.1002/mrm.20279. [DOI] [PubMed] [Google Scholar]

[R33] 33.Tuch DS, Reese TG, Wiegell MR, Makris N, Belliveau JW, Wedeen VJ. High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity. Magn. Reson. Med. 2002;vol. 48:577–582. doi: 10.1002/mrm.10268. [DOI] [PubMed] [Google Scholar]

[R34] 34.Ozarslan E, Vemuri BC, Mareci TH. Generalized scalar measures for diffusion MRI using trace, variance, and entropy. Magn. Reson. Med. 2005;vol. 53:866–876. doi: 10.1002/mrm.20411. [DOI] [PubMed] [Google Scholar]

[R35] 35.Stejskal EO, Tanner JE. Spin diffusion measurements: Spin echoes in the presence of a time-dependent field gradient. J. Chem. Phys. 1965;vol. 42:288–292. [Google Scholar]

[R36] 36.Xu D, Mori S, Shen D, van Zijl PC, Davatzikos C. Spatial normalization of diffusion tensor fields. Magn. Reson. Med. 2003;vol. 50:175–182. doi: 10.1002/mrm.10489. [DOI] [PubMed] [Google Scholar]

[R37] 37.Alexander DC. Inf. Process. Med. Imag. Glenwood Springs, CO: 2005. Maximum entropy spherical deconvolution for diffusion MRI; pp. 76–87. [DOI] [PubMed] [Google Scholar]

[R38] 38.Fukunaga K. Introduction to Statistical Pattern Recognition. 2nd ed. Boston, MA: Academic; 1990. [Google Scholar]

[R39] 39.Leow A, Huang S, Geng A, Becker J, Davis S, Toga A, Thompson P. Inverse consistent mapping in 3-D deformable image registration: Its construction and statistical properties. the Inf. Proces. Med. Imag. (IPMI); Glenwood Springs, CO. 2005. [DOI] [PubMed] [Google Scholar]

[R40] 40.Christensen GE, Rabbitt RD, Miller MI. Deformable templates using large deformation kinematics. IEEE Trans. Image Process.; 1996. Oct. pp. 1435–1447. [DOI] [PubMed] [Google Scholar]

[R41] 41.Miller MI. Computational anatomy: Shape, growth, and atrophy comparison via diffeomorphisms. NeuroImage. 2004;vol. 23:S19–S33. doi: 10.1016/j.neuroimage.2004.07.021. [DOI] [PubMed] [Google Scholar]

[R42] 42.Avants B, Gee JC. Geodesic estimation for large deformation anatomical shape averaging and interpolation. NeuroImage. 2004;vol. 23:S139–S150. doi: 10.1016/j.neuroimage.2004.07.010. [DOI] [PubMed] [Google Scholar]

[R43] 43.Beg MF, Miller MI, Trouvé A, Younes L. Computing large deformation metric mappings via geodesic flows of diffeomorphisms. Int. J. Comput. Vis. 2005;vol. 61:139–157. [Google Scholar]

[R44] 44.Vaillant M, Miller MI, Younes L, Trouve A. Statistics on diffeomorphisms via tangent space representations. NeuroImage. 2004;vol. 23:S161–S169. doi: 10.1016/j.neuroimage.2004.07.023. [DOI] [PubMed] [Google Scholar]

[R45] 45.Rose SE, McMahon KL, Janke AL, O’Dowd B, de Zubicaray G, Strudwick MW, Chalk JB. Diffusion indices on magnetic resonance imaging and neuropsychological performance in amnestic mild cognitive impairment. J. Neurol. Neurosurg. Psychiatry. 2006;vol. 77:1122–1128. doi: 10.1136/jnnp.2005.074336. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R46] 46.Jones DK, Horsfield MA, Simmons A. Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magn. Reson. Med. 1999;vol. 42:515–525. [PubMed] [Google Scholar]

[R47] 47.Studholme C, Cardenas V, Schuff N, Rosen H, Miller B, Weiner M. Detecting spatially consistent structural differences in Alzheimer’s and fronto temporal dementia using deformation morphometry. Int. Conf. Med. Image Comput. Comput. Assist. Intervention (MICCAI); 2001. pp. 41–48. [Google Scholar]

[R48] 48.Joshi S, Davis B, Jomier M, Gerig G. Unbiased diffeomorphic atlas construction for computational anatomy. NeuroImage. 2004;vol. 23:S151–S160. doi: 10.1016/j.neuroimage.2004.07.068. [DOI] [PubMed] [Google Scholar]

[R49] 49.Kochunov P, Lancaster JL, Thompson P, Woods R, Mazziotta J, Hardies J, Fox P. Regional spatial normalization: Toward an optimal target. J. Comput. Assist. Tomogr. 2001;vol. 25:805–816. doi: 10.1097/00004728-200109000-00023. [DOI] [PubMed] [Google Scholar]

[R50] 50.Kochunov P, Lancaster J, Thompson P, Toga AW, Brewer P, Hardies J, Fox P. An optimized individual target brain in the Talairach coordinate system. NeuroImage. 2002;vol. 17:922–927. [PubMed] [Google Scholar]

[R51] 51.Chiang M-C, Dutton RA, Hayashi KM, Lopez OL, Aizenstein HJ, Toga AW, Becker JT, Thompson PM. 3-D pattern of brain atrophy in HIV/AIDS visualized using tensor-based morphometry. NeuroImage. 2007;vol. 34:44–60. doi: 10.1016/j.neuroimage.2006.08.030. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R52] 52.Christensen GE, Johnson HJ. Invertibility and transitivity analysis for nonrigid image registration. J. Electron. Imag. 2003;vol. 12:106–117. [Google Scholar]

[R53] 53.Thevenaz P, Ruttimann UE, Unser M. A pyramid approach to subpixel registration based on intensity. IEEE Trans. Image Process; 1998. pp. 27–41. [DOI] [PubMed] [Google Scholar]

[R54] 54.Basser PJ, Pajevic S. Statistical artifacts in diffusion tensor MRI (DT-MRI) caused by background noise. Magn. Reson. Med. 2000;vol. 44:41–50. doi: 10.1002/1522-2594(200007)44:1<41::aid-mrm8>3.0.co;2-o. [DOI] [PubMed] [Google Scholar]

[R55] 55.Sheskin D. Handbook of Parametric and Nonparametric Statistical Procedures. 3rd ed. Boca Raton, FL: Chapman Hall/CRC; 2004. [Google Scholar]

[R56] 56.Storey JD. A direct approach to false discovery rates. J. R. Stat. Soc.: Series B (Stat. Methodol.) 2002;vol. 64:479–498. [Google Scholar]

[R57] 57.Smith SM, Jenkinson M, Johansen-Berg H, Rueckert D, Nichols TE, Mackay CE, Watkins KE, Ciccarelli O, Cader MZ, Matthews PM, Behrens TE. Tract-based spatial statistics: Voxelwise analysis of multi-subject diffusion data. NeuroImage. 2006;vol. 31:1487–1505. doi: 10.1016/j.neuroimage.2006.02.024. [DOI] [PubMed] [Google Scholar]

[R58] 58.Descoteaux M, Angelino E, Fitzgibbons S, Deriche R. A fast and robust ODF estimation algorithm in Q-ball imaging. the IEEE Int. Symp. Biomed. Imag.: From Nano to Macro (ISBI); Arlington, VA. 2006. [Google Scholar]

[R59] 59.Sijbers J, den Dekker AJ, Van Audekerke J, Verhoye M, Van Dyck D. Estimation of the noise in magnitude MR images. Magn. Reson. Imag. 1998;vol. 16:87–90. doi: 10.1016/s0730-725x(97)00199-9. [DOI] [PubMed] [Google Scholar]

[R60] 60.Christensen GE, Rabbitt RD, Miller MI, Joshi SC, Grenander U, Coogan TA, Essen DCV. Topological properties of smooth anatomic maps. the 14th Conf. Inf. Process. Med. Imag. Ile de Berder; France. 1995. [Google Scholar]

[R61] 61.Bro-Nielsen M, Gramkow C. Fast fluid registration of medical images. Proc. 4th Int. Conf. Vis. Biomed. Comput.; 1996. pp. 267–276. [Google Scholar]

[R62] 62.Shattuck DW, Leahy RM. BrainSuite: An automated cortical surface identification tool. Med. Image Anal. 2002;vol. 6:129–142. doi: 10.1016/s1361-8415(02)00054-3. [DOI] [PubMed] [Google Scholar]

[R63] 63.Leporé N, Brun CA, Pennec X, Chou YY, Lopez OL, Aizenstein HJ, Becker JT, Toga AW, Thompson PM. Mean template for tensor-based morphometry using deformation tensors. Int. Conf. Med. Image Comput. Comput. Assist. Interv. (MICCAI); Brisbane, Australia: 2007. [DOI] [PubMed] [Google Scholar]

[R64] 64.Cook PA, Bai Y, Nedjati-Gilani S, Seunarine KK, Hall MG, Parker GJ, Alexander DC. Camino: Open-source diffusion-MRI reconstruction and processing. 14th Sci. Meeting Int. Soc. Magn. Reson. Med.; Seattle, WA. 2006. p. 2759. [Google Scholar]

PERMALINK

Fluid Registration of Diffusion Tensor Images Using Information Theory

Ming-Chang Chiang

Alex D Leow

Andrea D Klunder

Rebecca A Dutton

Marina Barysheva

Stephen E Rose

Katie L McMahon

Greig I de Zubicaray

Arthur W Toga

Paul M Thompson

Abstract

I. INTRODUCTION

II. REGISTRATION SCHEME FOR DIFFUSION TENSORS

A. DTI Registration Using the Symmetrized KL-Divergence

B. Applicability of the Symmetrized KL-Divergence to Higher-Order Diffusion Models

Fig. 1.

C. Reorientation of Diffusion Profiles

D. Inverse—Consistent Fluid Regularization

III. EXPERIMENTAL METHODS

A. DTI Acquisition

B. Preprocessing of DTI

C. Target Selection

D. Fluid Registration of Diffusion Tensor Images

Fig. 2.

IV. EVALUATION OF THE REGISTRATION ALGORITHM

A. Evaluation of Registration Accuracy

Measurement of Tensor Overlap

Dyadic Coherence of Principal Eigenvectors

Variance of FA

B. IC Error of the Deformation Fields

V. RESULTS

A. Inverse–Consistent Fluid Registration for DTI

Fig. 3.

Fig. 4.

Fig. 5.

Fig. 6.

Fig. 7.

Fig. 8.

B. Evaluation of Registration Performance

Fig. 9.

Fig. 10.

Fig. 11.

C. Application of the Symmetrized KL-Divergence to Higher-Order Diffusivity Models

Synthetic examples

Fig. 12.

HARDI Data

Fig. 13.

VI. DISCUSSION

Appendix

APPENDIX A

Derivation of (7)

APPENDIX B

Computing sKL in terms of Spherical Harmonic Series

APPENDIX C

Derivation of (B4)

APPENDIX D

Iterative Computations for Forward and Backward Displacement Fields

APPENDIX E

The pFDR Method for Testing Overall Significance of Voxel Statistics

APPENDIX F

Derivation of the Driving Force Based on IPL and IPD

Contributor Information

REFERENCES

ACTIONS

PERMALINK

RESOURCES

Similar articles

Cited by other articles

Links to NCBI Databases

Derivation of the Driving Force Based on IP_L and IP_D