A Framework for Linear and Non-Linear Registration of Diffusion-Weighted MRIs Using Angular Interpolation

Registration of diffusion-weighted magnetic resonance images (DW-MRIs) is a key step for population studies, or construction of brain atlases, among other important tasks. Given the high dimensionality of the data, registration is usually performed by relying on scalar representative images, such as the fractional anisotropy (FA) and non-diffusion-weighted (b0) images, thereby ignoring much of the directional information conveyed by DW-MR datasets itself. Alternatively, model-based registration algorithms have been proposed to exploit information on the preferred fiber orientation(s) at each voxel. Models such as the diffusion tensor or orientation distribution function (ODF) have been used for this purpose. Tensor-based registration methods rely on a model that does not completely capture the information contained in DW-MRIs, and largely depends on the accurate estimation of tensors. ODF-based approaches are more recent and computationally challenging, but also better describe complex fiber configurations thereby potentially improving the accuracy of DW-MRI registration. A new algorithm based on angular interpolation of the diffusion-weighted volumes was proposed for affine registration, and does not rely on any specific local diffusion model. In this work, we first extensively compare the performance of registration algorithms based on (i) angular interpolation, (ii) non-diffusion-weighted scalar volume (b0), and (iii) diffusion tensor image (DTI). Moreover, we generalize the concept of angular interpolation (AI) to non-linear image registration, and implement it in the FMRIB Software Library (FSL). We demonstrate that AI registration of DW-MRIs is a powerful alternative to volume and tensor-based approaches. In particular, we show that AI improves the registration accuracy in many cases over existing state-of-the-art algorithms, while providing registered raw DW-MRI data, which can be used for any subsequent analysis.

1. Introduction

Diffusion-weighted magnetic resonance imaging (DW-MRI) is a non-invasive imaging technique that measures the diffusion of water molecules in biological tissues, such as the brain’s white matter, along several gradient directions. This information can be used to estimate the local orientation of fiber bundles, since water diffusion is favored along the fiber orientation, providing critical information for neuroscience and clinical studies (Jones et al., 1999; Basser and Jones, 2002; Filler, 2009).

Registration of DW-MRIs is a key step in population studies or brain atlases construction, among other important tasks. However, registration of such data is more challenging than for three-dimensional scalar images, not only because of its high dimensionality (hundreds of volumes) or noise and artifacts present in MR scans (Jones et al., 1999; Behrens et al., 2003; Le Bihan et al., 2006; Stobbe and Beaulieu, 2011), but because the signal, and any subsequently estimated local fiber orientation model, must remain consistent with the underlying fiber geometry after image transformations (Zhang et al., 2006; Barmpoutis et al., 2007; Cheng et al., 2009; Dhollander et al., 2010; Verma and Bloy, 2010; Yap et al., 2010; Du et al., 2011; Geng et al., 2011; Raffelt et al., 2011).

The simplest approach to DW-MRI registration consists in registering these datasets using a single transformation obtained from a representative volume, such as the Fractional Anisotropy (FA) (Pierpaoli and Basser, 1996), or the non-diffusion-weighted (b0) image. However, DW-MR registration using a single volume can be unsatisfactory, since it disregards the fiber orientation information provided by the diffusion-weighted volumes.

Another approach consists in registering instead the DTI estimated from the DW-MRI using the Stejskal and Tanner (1965) equation and assuming that the diffusion process follows a Gaussian distribution (Basser et al., 1994). Besides the spatial registration of DTI data, the images must also be reoriented so as to be consistent with the transformations made to the anatomy (Alexander et al., 2001; Zhang et al., 2006).

In order to overcome known limitations of the diffusion tensor model (Skare et al., 2000; Descoteaux et al., 2006; Zhang et al., 2006; Barmpoutis et al., 2007; Hess and Mukherjee, 2007; Koay et al., 2009), higher order models have been proposed (Barmpoutis et al., 2007; Cheng et al., 2009; Dhollander et al., 2010; Verma and Bloy, 2010; Yap et al., 2010; Du et al., 2011; Geng et al., 2011; Raffelt et al., 2011). Nevertheless, relying on such diffusion models might not completely capture the information contained in the raw data and could therefore affect the registration accuracy. In addition, these models are not generic and rely on specific assumptions.

The concept of “Angular Interpolation” (AI) was proposed by Tao and Miller (2006) to perform linear registration of the raw DW-MRIs, without imposing any specific diffusion model. The spatial registration with an angular interpolation correction of the image intensities attempts to match the registration of the underlying fiber structures. The evaluation of AI was entirely qualitative in Tao and Miller (2006). In this work, we extend the AI algorithm to non-linear registration and perform a wide range of tests that include human brain DW-MRIs undergoing known synthetic linear and non-linear transformations. We also evaluate our algorithm on pairs of human brain DW-MRIs to quantify the registration accuracy of AI versus b0-based registration (FLIRT¹ and FNIRT² in FSL³), and two well-known libraries that support registration of DTIs: DTI-TK⁴ (Zhang et al., 2006) and MedINRIA⁵. In addition, we provide an efficient C++ implementation of the algorithms in FSL, an open source and publicly available library that provides several analysis tools for FMRI, MRI, and diffusion MRI brain imaging data. Our implementation also handles the registration of DW-MRI datasets obtained with different gradient tables.

A recent work (Wang et al., 2011) compares the registration accuracy of several state of the art non-linear registration algorithms: FNIRT, MedINRIA, DTI-TK, and several demons algorithms using real brain images and an atlas. The conclusions of this work were that DTI-TK has the best registration performance. We will demonstrate that AI has a competitive performance, while providing registration of the raw DW-MRI data, not just the tensors.

Some other works have been recently published addressing the issue of gradient reorientation of DW-MRIs, in the context of image registration (Dhollander et al., 2010; Yap et al., 2010), and atlas building (Bouix et al., 2010). In particular, Yap et al. (2010) compute rotated gradient directions to estimate the rotated tensors from the spatially registered DW-MRIs. Dhollander et al. (2010) and Bouix et al. (2010) model the DW-MR datasets as a weighted superposition of spherical harmonic polynomials, updating the intensity of the registered DW-MRIs with the estimated weights and affine transformed basis. In fact, Dhollander et al. (2010) shows, using a synthetic example, that AI can produce undesirable rotation effects on fibers that have suffered pure shearing. This is a well-known problem in DTI registration (Alexander et al., 2001): pure shearing contains a complex rotational component not accounted for, when extracting the rigid rotational component of the affine matrix. In order to account for the complex rotation effect that shearing, stretching, or non-uniform scaling can have on the fiber orientation, Alexander et al. (2001) proposed the preservation of principal direction (PPD) algorithm. Dhollander et al. (2010) work can be seen as an extension of the PPD algorithm to their model of DW-MR datasets. However, as Zhang et al. (2006) and Alexander et al. (2001) indicate, there is little difference between using PPD or pure rigid rotation to register real DT-MR images and that is why it has been successfully used in the past (Zhang et al., 2006; Thomas Yeo et al., 2009), with a much lower computational cost.

Based on the previous considerations, we choose in this work to use the pure rigid rotational component of the affine matrices⁶. Nevertheless, AI can also be implemented using PPD to extract the rotation of general affine transformations to rotate independently the gradients in the image, making it similar to the work of Dhollander et al. (2010) and Bouix et al. (2010), with the advantage of not having to model the image as a superposition of polynomials.

The main contributions of this work are:

• The extension of AI to non-linear registration (Section 2.2).

• Improved computational performance of AI and quality of the associated FA estimate (see Section 2 in Appendix).

• Handling of DW-MR images with different gradient tables (Sections 2.1, 2.3).

• Comparative study of AI versus b0-based registration and two well-known tensor-based registration methods, MedINRIA and DTI-TK (Section 2.3).

• Public domain software to run AI within FSL⁷, including AI correction to b0-based and DTI-based spatial only registration. This code will be released upon publication as a plug-in for FSL⁸.

The concept of angular interpolation for linear registration is introduced in Section 2.1. In Section 2.2, we extend AI to non-linear registration by extending the FNIRT model to register diffusion-weighted MR images using AI. Some important details on the implementation of the registration algorithms in FSL are provided in the Appendix. The experiments performed, as well as the images used in this work are presented in Section 2.3. The results of these experiments are presented and discussed in Section 3. We present the conclusions of this work in Section 4. Finally, note that the Appendix includes more details about the implementation and complementary figures.

2. Materials and Methods

2.1. Affine Registration with Angular Interpolation

For an n-voxel DW-MRI dataset taken with m different gradient directions, the MR attenuation due to water diffusion in organized tissues (such as the brain white matter) can be modeled as (Stejskal and Tanner, 1965; Basser and Pierpaoli, 1996; Jones et al., 1999; Hasan et al., 2001)

$S\left(x_i,g_j,\tau \right)=S_0{\int }_{{ℝ}^3}P\left(x_i+r|\tau \right)e^{-2\pi \text{i}{q}_j\cdot r}dr,$ (1)

where, $x_i\in {ℝ}^3,i=1,\dots ,n$ are the discrete spatial samples, $g_j\in {ℝ}^3,j=1,\dots ,m$, are the diffusion gradients, τ the molecular diffusion time, S₀ the anatomical image taken with no gradient applied, q _j = γδg _j/2π is the displacement reciprocal vector (with γ the gyromagnetic ratio of water protons and δ the duration of the diffusion gradients), r the displacement vector relative to x _i, (·) is the dot product between vectors, and P(x _i + r|τ) the ensemble average propagator (EAP) at x _i, whose shape depends on the underlying fiber bundle structures. Without any loss of generality, we can consider that ||g _j||₂ = 1, hence from (1), the DW-MRI signal can be seen as samples of the Fourier transform of the EAP taken along different gradient directions on the unit sphere. Hence, the intensity I(x _i, g _j) of the DW-MRI is a function of the position, gradient orientation, and the underlying fiber structure (given by the EAP at each voxel).

From now on, we assume that the volumes of the DW-MR datasets have been properly aligned to correct for possible mismatches (Section 2.3.1), due to eddy current distortions, geometric distortions, and head motion during the data acquisition. Hence, it is assumed here that the DW-MRI datasets measure the true underlying three-dimensional structures (through diffusion), and any linear or non-linear spatial transformation should be applied to every volume in the data set to reflect the corresponding spatial transformation of those structures.

Let I_ref(x _i, g _j) and I_test(x _i, g _j) be the reference and test DW-MRIs, respectively. Also, let A be a given transformation matrix that registers the test to the reference DW-MRI. The spatial registration of the test DW-MRI would be ${Ĩ}_{test}\left(A^{-1}{\mathbf{\text{x}}}_i,{\mathbf{\text{g}}}_j\right),$ requiring interpolation⁹ to compute the intensity at the spatial coordinates A⁻¹ x _i. Hence, all the m-volumes in the DW-MRI can be spatially registered using the same affine registration. However, performing only spatial registration ignores the fact that the underlying fibers have changed their orientation with respect to the gradient orientations g _j and hence, in general, the intensity of the reference image will not be the same as the intensity obtained after spatial registration of the test image.

The most widely used reorientation strategy in DTI consists in reorienting the tensors as RDR ^T (Alexander et al., 2001; Zhang et al., 2006), where D is the second order diffusion tensor and R is the orthonormal matrix corresponding to the rotational component of A. The rotational component of a matrix can be obtained using the Finite Strain (FS) method, based on the polar decomposition (Higham, 1986), A = RS, where S is the strain component of A. A more accurate estimation of the rotational component of A can be obtained with the preservation of principal direction (PPD) algorithm, proposed by Alexander et al. (2001). Here, and as in Zhang et al. (2006), we use the FS method given its good performance by comparison with PPD, when the images are real DW-MRIs (Alexander et al., 2001; Zhang et al., 2006). Also, FS provides an analytical decomposition of A that allows the use of closed form derivatives, necessary when using gradient descent minimization methods (as in FNIRT), and it is simpler and computationally more efficient than PPD. We must point out here that in the original paper on AI (Tao and Miller, 2006), the authors do not use the rotational component of A, but use instead Ag /||Ag||₂ to modify the gradient directions (see Section 2.1). However, as Alexander et al. (2001) indicate, the shape of regions in the image can change (due to spatial registration), but the underlying tissue microstructure can only change through pure rotation.

We will explain now what angular interpolation is. From the fact that ${\mathbf{\text{g}}}_j\cdot \left(\mathbf{\text{Rr}}\right)=\left({\mathbf{\text{R}}}^T{\mathbf{\text{g}}}_j\right)\cdot \mathbf{\text{r}}$ and (1), we can see that sampling the reoriented fibers using the gradients g _j is equivalent to sampling the unrotated fibers using the rotated gradients ${\mathbf{\text{R}}}^T{\mathbf{\text{g}}}_j$. Hence, the key idea in AI is that under an affine transformation, the rotation of the underlying fiber orientations can be obtained by rotating the applied gradient ${\mathbf{\text{g}}}_j^{\prime }={\mathbf{\text{R}}}^T{\mathbf{\text{g}}}_j$, so that in the new coordinate system provided by ${\mathbf{\text{g}}}_j^{\prime }$, the underlying fiber orientations are rotated by R, as desired. This seems trivial, but the advantage of using AI is that we do not need to compute tensors from the (noisy) images or assume any other diffusion model. Also, we do not need to know a priori the fiber orientations, or if there are fiber crossings on a given voxel. We simply rotate the known gradients as ${\mathbf{\text{R}}}^T{\mathbf{\text{g}}}_j$. and compute, by angular interpolation, the DW-MRI corresponding to those new orientations, based on the full original test DW-MRI dataset. Notice here that AI is done on the unit sphere (||g _j||₂ = ||R ^Tg _j||₂ = 1) and that is why it is called angular interpolation. Notice also that AI does not spatially affect the intensity of the volumes, it interpolates among the volumes corresponding to the known gradient directions g _j. Hence, AI provides a natural way to register DW-MR images. The usual spatial interpolation is complemented by angular interpolation to reflect the change in intensity due to fiber reorientation.

In practice, we do not have an oracle that provide us with the optimal affine transformation to register the test image. We need to estimate the affine matrix by minimizing a cost function that measures the distance between the reference and registered test images. Hence, we want to find the optimal affine transformation A that minimize the cost function,

$C\left(A\right)=\sum _{i=1}^n\sum _{j=1}^m\mathcal{𝒟}\left(I_{ref}\left(x_i,g_j\right),{Ĩ}_{test}\left(A^{-1}{x}_i,R^T{g}_j\right)\right),$ (2)

where, D is an appropriate distance metric, such as Euclidean distance, cross-correlation, or Mutual Information.

In general, angular interpolation of DW-MR images can be done as (Tao and Miller, 2006),

$\dot{I}\left(x_i,g_j\right)=\frac{{\sum }_{k=1}^p{\alpha }_{k}I\left(x_i,g_k\right)}{{\sum }_{k=1}^p{\alpha }_k},i=1,\dots ,n,j=1,\dots ,m,$ (3)

where α_k are the interpolation weights assigned to each of the p volumes in the test image, defined as non-linear functions of $co{s}^{-1}\left({\mathbf{\text{g}}}_j\cdot {\mathbf{\text{g}}}_k\right)$, the geodesic distance on the unit sphere between the desired target gradient, g _j, and the available gradients in the image, g _k (see Tao and Miller (2006) for details on the non-linear function used). Since g _j can be any gradient direction, we can use AI to match two DW-MR images with different gradient tables, by defining g _j as the gradients in the reference image and g _k the available gradients in the test image. From now on, we assume that the gradients in the test image have been matched to the gradients in the reference image (using AI) and hence, p = m.

The linear registration algorithm in FSL is FLIRT (Jenkinson et al., 2002), which provides single volume registration using Euclidean distance, Cross-correlation, and Mutual Information metrics. Our extension of FLIRT to register DW-MRIs is based on AI, where, besides simultaneous spatial registration of all the volumes, AI is also applied on all gradient directions, at each step of the minimization, as indicated in (3), where can be any of the distance metrics provided by FLIRT. We will explain in much more detail the extension to non-linear registration in the next section, since there, the implementation is more involved with the particular cost function and small displacement model used in FSL.

2.2. Non-Linear Registration with Angular Interpolation

In the previous section, we introduced the concept of angular interpolation to linearly register DW-MRIs. This idea can be extended to non-linear registration of DW-MRIs, if the non-linear transformation is locally affine. Let T(x) = x + u(x) be a non-linear transformation applied at x ∈ ℛ³, where u(x) is the so called deformation field, then (Irfanoglu et al., 2008),

$\begin{array}{llll}T\left(x\right)& =x+u\left(x\right)=A\left(x\right)x+t\left(x\right),& & \\ A\left(x\right)& =I+J_u\left(x\right),J_u\left(x\right)=\frac{du\left(x\right)}{dx},& \end{array}$ (4)

where, A(x) is the equivalent local affine transformation at the spatial coordinate x, t(x) a pure translation vector, I the identity matrix, and J_u the Jacobian of u. Hence, the deformation field can be decomposed into a local affine transformation plus a translation field.

The non-linear registration tool in FSL is FNIRT (Andersson et al., 2007). However, FNIRT can only be used to register scalar volumes¹⁰ and it works under the assumption of a small deformation field. Here, we will extend FNIRT to the non-linear registration of DW-MRIs using angular interpolation (AI). We extend here FNIRT, based on the derivation made by Andersson et al. (2007), for the case of scalar volumes. In this section, we will refer only to our FNIRT-based AI model and hence, the equations apply to DW-MRIs, not just single volumes. As in the case of affine registration, where the affine matrix applies to all the volumes in the image, the deformation field here is three-dimensional and also applies to the whole DW-MR image.

FNIRT models the deformation field as a linear superposition of three-dimensional cubic B-splines (Andersson et al., 2007),

$u\left(x_i\right)=\sum _k{w}_k{B}_k\left(x_i\right),k=1,\dots ,s,$ (5)

where w _k is the weight along each coordinate for the B-spline B_k, and s the number of splines, which depends on the spatial resolution of the field. Then, the non-linearly warped test image ${Ĩ}_{test}\left({\mathbf{\text{x}}}_i+\mathbf{\text{u}}\left({\mathbf{\text{x}}}_i\right),{\mathbf{\text{g}}}_j\right)={Ĩ}_{test}\left({\mathbf{\text{x}}}_i,{\mathbf{\text{g}}}_j,\mathbf{\text{w}}\right)$ is a function of the position, the gradient orientations, and w, the vector of all the B-spline coefficients defining the deformation field. Notice that all the diffusion images are in ℛ^{m × n}, where n correspond to the number of voxels in the image and m the number of diffusion gradients. Hence, the entire DWI series can be expressed in vector form as

$\begin{array}{c}{I}_{ref}={\left[I_{ref}^T\left(g_1\right)\dots I_{ref}^T\left(g_m\right)\right]}^T\in {ℝ}^{m\times n},\\ {\widetilde{\mathbf{\text{I}}}}_{test}\left(w\right)={\left[{\widetilde{\mathbf{\text{I}}}}_{test}^T\left(g_1,w\right)\dots {\widetilde{\mathbf{\text{I}}}}_{test}^T\left(g_m,w\right)\right]}^T\in {ℝ}^{m\times n},\end{array}$ (6)

where

$\begin{array}{llll}{I}_{ref}\left(g_i\right)& ={\left[I_{ref}\left(x_1,g_i\right)\dots I_{ref}\left(x_n,g_i\right)\right]}^T\in {ℝ}^n,& & \\ {\widetilde{\mathbf{\text{I}}}}_{test}\left(g_i,w\right)& ={\left[{Ĩ}_{test}\left(x_1,g_i,w\right)\dots {Ĩ}_{test}\left(x_n,g_i,w\right)\right]}^T\in {ℝ}^n.& & \end{array}$

The cost function to be minimized in FNIRT is the mean square error¹¹

$C\left(w\right)=\frac{1}{mn}{E}^{T}E,E={\widetilde{\mathbf{\text{I}}}}_{test}\left(w\right)-I_{ref},$ (7)

where E is the error between the reference and test images. We are interested in estimating the deformation field, defined by the coefficients w, such that it minimizes C(w).

We are going to derive next a solution to (7) that extends the equations derived in (Andersson et al., 2007) for single volume registration to registering diffusion data sets. FNIRT solves (7) using a Gauss-Newton minimization method,

$w^{l+1}=H{\left(w^l\right)}^{-1}\nabla C\left(w_l\right),\in {ℝ}^s$ (8)

where l is the iteration index and H is the Hessian, which can be approximated in terms of the Jacobian as $H\left(w^k\right)\approx \frac{2}{mn}J{\left(w^k\right)}^{T}J\left(w^k\right)$ (Fletcher, 1987), with the Jacobian given by

$\begin{array}{llll}J\left(w\right)& =\frac{dE}{dw}=\frac{d\left({\widetilde{\mathbf{\text{I}}}}_{test}\left(w\right)\right)}{dw}& & \\ & ={\left[\frac{d\left({\widetilde{\mathbf{\text{I}}}}_{test}^T\left(g_1,w\right)\right)}{dw}\dots \frac{d\left({\widetilde{\mathbf{\text{I}}}}_{test}^T\left(g_m,w\right)\right)}{dw}\right]}^T\in {ℝ}^{m\left(n\times s\right)}.& \end{array}$ (9)

Let $x^{\prime }=x+u\left(w\right)$,

$\frac{d\left({\widetilde{\mathbf{\text{I}}}}_{test}\left(g_j,w\right)\right)}{dw}={\left[\frac{\partial {Ĩ}_{test}\left({x^{\prime }}_i,g_j\right)}{\partial w_k}\right]}_{ik}=J^j\left(w\right)\in {ℝ}^{n\times s}.$ (10)

By the chain rule and (5),

$\begin{array}{llll}{J}^j\left(w\right)& ={\left[\frac{\partial {Ĩ}_{test}\left(x^{\prime },g_j\right)}{\partial x^{\prime }}{\left|\right.}_{{x^{\prime }}_i}\odot \frac{\partial x^{\prime }}{\partial w_k}{\left|\right.}_{x_i}\right]}_{ik}& & \\ & ={\left[\frac{\partial {Ĩ}_{test}\left(x,g_j\right)}{\partial x}{\left|\right.}_{{x^{\prime }}_i}{B}_k\left({x^{\prime }}_i\right)\right]}_{ik},& \end{array}$ (11)

where ⊙ is the Hadamard product and $\frac{d\left({\widetilde{\mathbf{\text{I}}}}_{test}\left(g_j,w\right)\right)}{dw}$ corresponds to the gradient of the moving image being warped by the deformation field, conveniently expressed as a vector derivative to highlight the derivative chain rule, and following the derivation in Andersson et al. (2007). This gradient is computed in FNIRT using finite differences on an image that is the non-linearly warped current test image (using the deformation field). Since the deformation field applies the same warp to all the diffusion images, we compute this gradient on each diffusion gradient, g_j, as indicated in (11).

From (9) and (10),

$H\left(w\right)=\frac{2}{mn}\sum _j{\left(J^j\left(w\right)\right)}^T{J}^j\left(w\right)\in {ℝ}^{s\times s},$ (12)

where J ^j(w) is given by (11). Now, the gradient of the cost function is given by (see (7) and (9))

$\nabla C\left(w\right)=\frac{dC\left(w\right)}{dw}=\frac{2}{mn}\frac{d{E}^T}{dw}E=\frac{2}{mn}{J}^T\left(w\right)\left({\widetilde{\mathbf{\text{I}}}}_{test}\left(w\right)-I_{ref}\right)\in {ℝ}^s.$ (13)

Until now, we have extended FNIRT to DW-MRIs and provided the equations to update w, from (8), (12), and (13), based on the deformation field model in FNIRT. However, this minimization algorithm still does not account for the intensity correction provided by angular interpolation. Hence, we need to extract the rotational component from the local affine matrix (see (4)),

$I+J_u\left(x_i\right)=R\left(x_i\right)S\left(x_i\right),i=1,\dots ,n,$ (14)

where R(x _i) and S(x _i) are respectively the rotational and strain components at each voxel. Hence, the transformed test image becomes

${\dot{\mathbf{\text{I}}}}_{test}\left(w\right)={\left[{\dot{\mathbf{\text{I}}}}_{test}^T\left(R^T\left(x_1\right)g,w\right)\dots {\dot{\mathbf{\text{I}}}}_{test}^T\left(R^T\left(x_n\right)g,w\right)\right]}^T,$ (15)

where, g = [g₁ … g _m] and ${\dot{I}}_{test}\left({\text{R}}^T\left({\text{x}}_i\right)\text{g},\text{\hspace{0.17em}w}\right)={\left[{\dot{I}}_{test}\left({\text{x}}_i,{\text{R}}^T\left({\text{x}}_i\right){\text{g}}_j\text{w}\right)\right]}_j$ is the spatially and angularly interpolated test DW-MRI at x _i. Now, the Jacobian is given by

$J\left(w\right)={\left[\frac{\partial \left({\dot{\mathbf{\text{I}}}}_{test}^T\left(R^T\left(x_1\right)g,w\right)\right)}{\partial w}\dots \frac{\partial \left({\dot{\mathbf{\text{I}}}}_{test}^T\left(R^T\left(x_m\right)g,w\right)\right)}{\partial w}\right]}^T.$ (16)

As before, by the chain rule and (5),

$J^j(w)={\left[{\frac{\partial {\dot{I}}_{test}(x_i,R^T(x_i)g_j,w)}{\partial x}|}_{{x^{\prime }}_i}{B}_k({x^{\prime }}_i)\right]}_{jk}.$ (17)

The Hessian can be computed as indicated in (12), and the gradient of the cost function as in (13), but with J as defined in (16), (17), and ${\dot{\mathbf{\text{I}}}}_{test}\left(w\right)$ as defined in (15). It remains to define the local angular interpolation, given by,

${\dot{I}}_{test}\left(x_i,R^T\left(x_i\right)g_j,w\right)=\frac{{\sum }_{k=1}^m{\alpha }_k{Ĩ}_{test}\left(x_i,g_k,w\right)}{{\sum }_{k=1}^m{\alpha }_k},j=1,\dots ,m.$ (18)

As indicated in Equations (7), (12), and (13), the cost function, Hessian, and gradient are functions of the B-splines representing the displacement field. Hence, we can use the multi-scale Levenberg-Marquart (LM) minimization algorithm used in FNIRT (Andersson et al., 2007) to iteratively estimate the displacement field. At every iteration of the LM algorithm, the diffusion volumes are spatially registered using the local affine matrix (14) and their intensity interpolated using AI (15). This minimization method is commonly employed by other algorithms in FSL. The algorithm for non-linear registration using AI is detailed in Algorithm 1.

Algorithm 1. Non-linear Registration with Angular Interpolation

Require: Initial deformation field u0(x)
 (zero if not provided).
  while C(w) decreases do {OPTIMIZATION}
   for i = 1 → n do
    Compute the Jacobian and Hessian of
      the error E using (17) and (12).
    Compute the rotational component
      R(xi) of I + J u (xi) as indicated
      in (14).
    Compute the overall cost C(w) as
      indicated in (7).
    Compute the gradient of the cost
      function ∇C(w) as indicated
      in (13).
    w ← H(w)−1∇ C(w).
    for j = 1 → m do
     Compute Ĩtest(xi, gj, w) using spatial
      interpolation.
    end for
    Compute ${\dot{I}}_{test}\left({\mathbf{\text{X}}}_i,{\mathbf{\text{R}}}^T\left({\mathbf{\text{x}}}_i\right){\mathbf{\text{g}}}_j,\mathbf{\text{w}}\right)$ using
     angular interpolation (15).
   end for
  end while

2.3. Experimental Details

2.3.1. Diffusion-weighted MRI acquisitions

Five diffusion-weighted MRI datasets were obtained on a 7T scanner (Magnex Scientific, UK) driven by a Siemens console (Erlangen, Germany), and using a Siemens head-gradient insert capable of 80 mT/m using a single refocused 2D single shot spin echo EPI sequence, 1.5 mm isotropic voxels, 128 directions at b = 1500 s/mm² and 15 b0s (Lenglet et al., 2012). We also use two DW-MRI datasets obtained on a 3T scanner (Siemens MAGNETOM Trio) using whole-body gradients capable of 40 mT/m, a single refocused 2D single shot spin echo EPI sequence, 2 mm isotropic voxels, 128 directions at b = 1000 s/mm² and 16 additional b0s (Zhan et al., 2013). Each DWI dataset was corrected for eddy current distortions and head motion using a 12 degree of freedom linear registration with the initial b0 image. They were also corrected for geometric distortions, by unwrapping the images using a field map. The diffusion gradients were reoriented to take into account the spatial transformations used for eddy current distortions correction [see details in Lenglet et al. (2012), Zhan et al. (2013)].

The nomenclature used here refers to these images as Subject1 to Subject6, the first four subjects correspond to the 7T data, while the last two correspond to the 3T data. Subject3 was scanned twice, at two different dates. We refer here to Subject3 and Subject3-2 to differentiate them. The research protocols used for the acquisition of the 7T and 3T data were approved by the Institutional Review Board (IRB) of the University of Minnesota. All subjects provided informed written consent prior to participating in the research protocols. We performed inter-subject registrations among all five datasets (subjects 1-4) acquired at 7T and separated inter-subject registrations between the pair of datasets acquired at 3T (Subjects 5 and 6).

There are 20 possible registrations between Subjects 1–4. Taking into account that Subject3 and Subject3-2 correspond to the same subject, we choose seven registration pairs (see Table 1), which cover all the intra-subject and inter-subject registration cases, without reversing the roles of the reference and test images, nor registering to the same subject twice. We also performed two registrations between Subjects 5 and 6 as indicated in Table 1.

TABLE 1

Table 1. Tested registration pairs.

Linear registrations in FSL are performed with FLIRT and the new FLIRT-based AI using cross-correlation as the distance measure between images (see Section 2.1).

2.3.2. Synthetic transformations

In order to obtain unbiased synthetic transformations, we use the advanced normalization tools (ANTs)¹² to warp spatially a DW-MRI dataset Subject2 using a pre-defined random set of affine and non-linear transformations. ANTs uses one of the top five MRI non-linear registration tools identified in a recent evaluation study Klein et al. (2009). We generated 20 random affine matrices (see Kannala et al., 2005, for instance) of the form A = RSD, where R is a pure rotation matrix, S an upper triangular skew matrix, and D a diagonal matrix providing scale change (see details in Section 1.1 of the Appendix). We also generated 20 random non-linear deformation fields, based on the work of Noblet et al. (2006) (see details in Section 1.2 of the Appendix). These non-linear deformation fields ensures a topology preserving transformation (Noblet et al., 2006). As indicated in the introduction, registration of DW-MRIs requires also intensity corrections to reflect the change in fiber orientation. To avoid bias in favor of our proposed technique, we do not use here our AI method to perform the intensity correction [Equation (18)]. Instead, we use an eighth order spherical harmonics (SH) representation of the DW-MRI datasets. Using this SH basis, it is straightforward to interpolate the intensities at different gradient directions (see Dhollander et al. (2010); Bouix et al. (2010); Kamath et al. (2012), for instance). We must point out here that while the SH representation of DW-MRIs can smooth out certain features of the diffusion signal (at each particular b-value), it provides at least an unbiased correction of the image intensities that will be used for all the other algorithms tested. In addition, we only use ANTs to perform the spatial warping of each diffusion direction with a known pre-defined transformation. We do not use ANTs to try to register the warped images, since that would bias the results toward ANTs, as it was used to perform the warp in the first place.

Figure 1A shows the main direction of diffusion in the reference image, color coded with an RGB colormap (red for left-right, green for antero-posterior, and blue for superior-inferior), with brightness modulated by the value of the FA. Figure 1B shows the main direction of diffusion (same color code as before) for a linearly transformed test image. Figure 1C shows the main direction of diffusion for a non-linearly warped test image, and Figure 1D shows the corresponding non-linear deformation field, color coded by the magnitude of the displacement.

FIGURE 1

Figure 1. (A) Reference image. (B) Affinely warped reference image. (C) Non-linearly warped reference image using a small displacement field, (D) Non-linear small deformation field.

2.3.3. Intra- and inter-subject registration

In the previous section, we used known linear and non-linear transformations, since it allows us to consider the reference image as ground truth and therefore evaluate the registration accuracy. On the other hand, there is no consensus on how to evaluate the registration accuracy for any pair of images (Kannala et al., 2005; Noblet et al., 2006). In fact, the optimal registration of a pair of images depends largely on the application and registering two different subjects should not reduce or eliminate the structural differences between different brains. Nevertheless, and for completeness, we also tested the registration algorithms for pairs of different DW-MRIs, even though registration errors and anatomical variability cannot be differentiated from each other in this case.

Non-linear registration between any pair of DW-MR images requires to conduct first an affine registration. We will use here the best affine registration found, so that all the non-linear registration algorithms start with the same linearly registered image.

2.3.4. Evaluation of results

We compare the results from linear and non-linear registration using b0-based, DTI-based, and DWI-based registration. The registration error in terms of the mean squared error (MSE) is given by

$MSE=\frac{1}{mn}\sum _{i=1,j=1}^{n,m}{\left|\right.\left|\right.\frac{I_{test}\left(x_i,g_j\right)-I_{ref}\left(x_i,g_j\right)}{I_{max}}\left|\right.\left|\right.}_2^2\times 100,$ (19)

where n is the total number of voxels, m the number of gradient directions in the reference image, and we are normalizing the error by the maximum intensity in the reference image I_max. Based on Wang et al. (2011), we also define the mean fiber orientation error as

$foe=\frac{1}{n}\sum _{i=1}^{n}co{s}^{-1}\left|⟨v_i^{test},v_i^{ref}⟩\right|,$ (20)

where ${\mathbf{\text{v}}}_i^{ref}$, ${\mathbf{\text{v}}}_i^{test}$ are respectively, the estimated fiber orientations of the reference and registered DW-MR images at the i-th voxel.

The fiber orientations of the linearly and non-linearly registered images correspond to the first eigenvector of the second order tensor model for the reference and registered test images. The second order tensors (DTI) required by DTI-TK were computed using DTIFIT from FSL. MedINRIA on the other hand, it uses its own algorithm to estimate the tensors from the diffusion volumes and diffusion gradients/b-value(s) information, so we use those within MedINRIA. DTI-TK and MedINRIA can register scalar and tensor volumes. Hence, the first eigenvector can be extracted from the registered tensors using DTI-TK interoperability tools, which also works for MedINRIA, since they are both NIFTI¹³ compliant. Given that FLIRT and FNIRT can only register scalar volumes, we generate registered DW-MRIs by applying the transformation found (using the b0 volume) by FLIRT and FNIRT to each one of the volumes in the DW-MR data set. Intensity correction was also applied to these spatially registered datasets using SH. From these registered DW-MRIs, we can compute the second order tensors using DTIFIT. Similarly, we can use the spatial transformation (affine matrix or deformation field) found by DTI-TK and MedINRIA to spatially register every volume in the DW-MR datasets, applying optionally SH intensity correction. This allows us to obtain registered DW-MR datasets for FLIRT, FNIRT, DTI-TK, and MedINRIA. The advantage of obtaining registered DW-MRIs for all methods considered here is that we can then compute, in an unbiased fashion, the registration MSE from the registered DW-MR datasets with and without intensity correction. It also allows us to compare the registration foe of these registered DW-MRIs with the foe of the registered DTIs and the proposed FLIRT-based and FNIRT-based AI registration algorithms.

The second order tensors (DTI) can be computed from the linearly and spatially only registered DW-MRIs (i.e., without SH correction), by providing the reoriented gradients ${{\mathbf{\text{g}}}^{\prime }}_j={\mathbf{\text{R}}}^T{\mathbf{\text{g}}}_j,j=1,\dots ,m$, to DTIFIT (see Section 2.1), where R is the rotational component of the affine matrix used to spatially register the DW-MRIs.

The registration foe was computed on forty¹⁴ regions of interest (ROIs) taken from the John Hopkins University (JHU) white matter atlas (Hua et al., 2008)¹⁵. The JHU white matter regions were registered (using the FA volumes) from the atlas to each one of the DW-MR reference datasets. The MSE was also computed on the the same brain regions.

3. Results

3.1. Linear Registration

Table 2 summarizes the results of all the linear registration experiments, as described in Section 2.3, in terms of the mean foe and MSE, with their corresponding standard deviation, within parenthesis. Since only two 3T datasets were used to illustrate the performance of the algorithms at lower spatial resolution, we do not report here the standard deviation of the registration error. FLIRT4D-AI indicates the proposed registration of (4D) DW-MR datasets using angular interpolation. The appended suffix – tensors indicates that the foes were estimated directly from the tensors registered using the corresponding registration method. The appended suffix – volumes indicates that the affine transformations found by the corresponding registration method were used to spatially register all the diffusion volumes of the test DW-MRIs. Finally, the appended suffix – SH indicates that SH intensity correction was performed on top of the spatial registration. There were 20 affine registrations performed, for all experiments performed using random synthetic affine warps (see Section 2.3.2) and 9 intra- and inter-subject affine registrations, seven of which are indicated in Table 1, and the remaining two correspond to registering back and forth the two 3T datasets (see Section 2.3.3).

TABLE 2

Table 2. Linear registration (best results are indicated in bold).

Conversion from the ITK¹⁶ matrix format used by MedINRIA to the standard 4 × 4 affine matrix was done using the c3d_affine_tool¹⁷. DTI-TK matrix format directly provides the 3 × 3 submatrix of the affine transformation plus the displacement. Hence, the rotational component of the affine matrix can be easily extracted in DTI-TK. The rotational component of the affine registration matrices found by all the algorithms was computed using the decompose_aff method in FSL. We used the baloo (Ourselin et al., 2000), which is the affine registration algorithm in MedINRIA 2.01¹⁸ and worked very well in all cases considered.

MedINRIA achieved the lowest foe as well as very low MSE using the affine transformation to spatially register each diffusion volume and then performing SH intensity correction. The MSE correlates for the most part with the foe, but not always. For instance, DTI-TK-volumes have the lowest MSE in Table 2, for the synthetic affinely warped datasets, but it also has the largest foe. This could indicate that minimizing the MSE does not necessarily reduces the foe. Also, The SH intensity correction helps reduce the foe in some, but not in all cases. This can be understood from the fact that if the spatial registration was not accurate, then the SH correction is also going to be inaccurate and could even increase the error, as can be appreciated in the FLIRT-volumes-SH results compared to FLIRT-volumes. Notice also that the foe is low in MedINRIA-volumes despite not having intensity correction using SH. This is due to the fact that we are providing the reoriented diffusion gradients, nonetheless, the foes and MSE are even lower using the SH intensity correction, which confirms the importance of adjusting the intensity of spatially registered diffusion volumes.

A paired t-test reveals that there is significant statistical difference in favor of SH intensity correction, which reduced the foe by almost 12% (p = 0.0015), in the synthetic transformations. FLIRT4D-AI has a much lower foe and MSE than FLIRT-volumes as well as a lower variance. However, no significant statistical difference was found in this case, given the large variance of FLIRT-volumes results. Overall, the results indicate that performing intensity correction is important, but the accuracy of the registration algorithm is obviously critical, as can be seen from the very good results of MedINRIA. The fact that MedINRIA did better than the other methods, including FLIRT4D-AI, could be attributed to the fact that its affine registration algorithm is superior to the other registration methods considered here.

Figure 2 shows a checkerboard comparison between the main diffusion direction (RGB color coded) of the reference DW-MRI (Subject2) and one of the synthetic affine transformations. Figure 3 shows the tensors (represented as RGB color coded glyphs) of the reference (Subject2) and the registered images indicated in Figure 2. The best performance of MedINRIA can be appreciated in these figures, where the checkerboard effect is lower than in the other methods, followed by DTI-TK, FLIRT4D-AI, and finally FLIRT-volumes. Figure 4 shows a checkerboard comparison between the main diffusion direction in Subject3-2 and Subject3 registered to Subject3-2. Figure 5 shows the tensors (represented as RGB color coded glyphs) of the reference (Subject3-2) and Subject 3 registered to Subject3-2. Since this is the same individual, the registration errors are much lower than between different subjects. In terms of the mean foe FLIRT4D-AI has the lowest error (10.1), followed by FLIRT4D (10.6), DTI-TK (11.0), and MedINRIA (11.7).

FIGURE 2

Figure 2. Checkerboard comparison of the RGB color coded main diffusion direction between Subject2 and a registered random affine transformation of Subject2 using (A) FLIRT-volumes-SH. (B) FLIRT4D-AI. (C) MedINRIA-volumes-SH, (D) DTI-TK-volumes-SH.

FIGURE 3

Figure 3. Comparison of the tensors (represented as RGB color coded glyphs) of Subject2 and the registered tensors of a random affine transformation of Subject2. (A) Subject2 FA indicating the region selected to show the tensors. (B) Subject2 tensors. Registered tensors using (C) FLIRT-volumes-SH, (D) FLIRT4D-AI, (E) MedINRIA-volumes-SH, and (F) DTI-TK-volumes-SH.

FIGURE 4

Figure 4. Checkerboard comparison of the RGB color coded main diffusion direction between Subject3-2 and Subject3 registered to Subject3-2 using (A) FLIRT-volumes-SH. (B) FLIRT4D-AI. (C) MedINRIA-volumes-SH, (D) DTI-TK-volumes-SH.

FIGURE 5

Figure 5. Comparison of the tensors (represented as RGB color coded glyphs) of Subject3-2 and the tensors Subject3 registered to Subject3-2. (A) Subject3-2 FA indicating the region selected to show the tensors. (B) Subject3-2 tensors. Registered tensors of Subject3 using (C) FLIRT-volumes-SH, (D) FLIRT4D-AI, (E) MedINRIA-volumes-SH, and (F) DTI-TK-volumes-SH.

3.2. Non-Linear Registration

Table 3 summarizes the results of all the non-linear registrations indicated in Section 2.3, in terms of the mean foe and MSE, with their corresponding standard deviation, except for the two 3T datasets, where only the mean is reported.

TABLE 3

Table 3. Non-linear registration (best results are indicated in bold).

The notation is the same as in Table 2, except that now FLIRT is replaced by FNIRT. Also, the non-linear registration used in MedINRIA corresponds to the diffeomorphic demons algorithm. In this case the appended suffix −volumes indicates that the corresponding registration algorithm used the estimated deformation field to non-linearly register each one of the diffusion volumes. The appended suffix −SH indicates that besides non-linear spatial registration, the intensity of the diffusion volumes is corrected using SH q-space interpolation and the rotational component of the equivalent local affine transformation [see Equation (14)]. The local affine transformation depends on the Jacobian of the deformation field [Equation (14)], which was computed using the deffield2jacobian method in FSL. The estimated deformation fields in MedINRIA and DTI-TK were converted to FSL format using the NIFTI Matlab toolbox¹⁹.

Table 3 indicates that the lowest foe for the synthetic non-linearly warped DW-MRIs correspond to FNIRT4D-AI, while the best foe for the intra and inter-subject registration was achieved by DTI-TK, which coincides with previous studies reporting that DTI-TK performed very well when compared to other tensor-based non-linear registration algorithms (Wang et al., 2011). MedINRIA (demons registration algorithm) did not perform as well as FNIRT or DTI-TK. Our results also indicate that SH correction did not improve any of the results obtained using non-linear spatial registration. This might be due to the rather coarse internal representation of deformation fields in FSL that is based on splines (as used in FNIRT). A paired t-test indicates that FNIRT4D-AI has a 6% lower error than FNIRT in the synthetic transformations (p = 0.001). The fact that DTI-TK performed better than FNIRT4D-AI for intra- and inter-subject registration possibly indicates that differences between the underlying algorithms are bigger than the gain obtained using AI or SH. This also corroborates our previous results on linear registration, where MedINRIA performed optimally. This is also inevitable for such a study, where various algorithms, with various strengths, and weaknesses are compared.

Figure 6 shows a checkerboard comparison between the main diffusion direction (RGB color coded) of the reference DW-MRI (Subject2) and one of the registered synthetic random non-linear warped DW-MRIs. Figure 7 shows the tensors (represented as RGB color coded glyphs) of the reference (Subject2) and the registered images indicated in Figure 6. The best performance of FNIRT4D-AI and FNIRT-volumes is clear over the other methods. Figure 8 shows a checkerboard comparison between the main diffusion direction in Subject3 and Subject2 registered to this subject. Figure 9 shows the tensors (represented as RGB color coded glyphs) of the reference (Subject3) and Subject2 registered to Subject3. As expected, there are many differences between the reference and registered image, since they correspond to different individuals. However, it is not difficult to see that the registered images using DTI-TK followed by FNIRT4D-AI have lower differences with respect to the reference, compared to FNIRT-volumes and MedINRIA.

FIGURE 6

Figure 6. Checkerboard comparison of the RGB color coded main diffusion direction between Subject2 and a registered random non-linear transformation of Subject2 using (A) FNIRT-volumes-SH. (B) FNIRT4D-AI. (C) MedINRIA-volumes-SH, (D) DTI-TK-volumes-SH.

FIGURE 7

Figure 7. Comparison of the tensors (represented as RGB color coded glyphs) of Subject2 and the tensors of a registered random non-linear transformation of Subject2. (A) Subject2 FA indicating the region selected to show the tensors. (B) Subject2 tensors. Registered tensors using (C) FLIRT-volumes-SH, (D) FLIRT4D-AI, (E) MedINRIA-volumes-SH, and (F) DTI-TK-volumes-SH.

FIGURE 8

Figure 8. Checkerboard comparison of the RGB color coded main diffusion direction between Subject3 and Subject2 registered to Subject3 using (A) FNIRT-volumes-SH. (B) FNIRT4D-AI. (C) MedINRIA-volumes-SH, (D) DTI-TK-volumes-SH.

FIGURE 9

Figure 9. Comparison of the tensors (represented as RGB color coded glyphs) of Subject3 and the tensors of Subject2 registered to Subject3. (A) Subject3 FA indicating the region selected to show the tensors. (B) Subject3 tensors. Registered tensors of Subject2 using (C) FLIRT-volumes-SH, (D) FLIRT4D-AI, (E) MedINRIA-volumes-SH, and (F) DTI-TK-volumes-SH.

4. Discussion

Previous work using higher order models (Barmpoutis et al., 2007; Cheng et al., 2009; Dhollander et al., 2010; Verma and Bloy, 2010; Yap et al., 2010; Du et al., 2011; Geng et al., 2011) compared their proposed model with DTI or single volume registration, using exactly the same underlying registration algorithm, where the difference is due only to the model used. We addressed here a slightly more challenging task using different registration algorithms (with the exception of FLIRT and FNIRT in FSL), in order to provide a useful and more practical comparison. The intrinsic differences among the different registration algorithms analyzed here can be larger than the differences found using AI or SH intensity correction. Consequently, AI should also be implemented within other registration algorithms, in order to make it competitive with other well-known and tested state of the art registration algorithms (based on second order tensors or higher order models such as the ODF), a representative set of which has been used in this study. Nonetheless, we have clearly demonstrated that AI improves the registration accuracy in many cases over existing state-of-the-art algorithms. We also provided a framework and software to perform FNIRT-based AI registration or just intensity correction, which will enable future comparative studies between various registration algorithms using different diffusion models.

Despite the strong differences between the different algorithms used here, the results indicate that registration of DW-MR datasets should include intensity correction within the minimization algorithm, as proposed in our extension of FNIRT. Intensity correction as a post-processing step after spatial correction does not necessarily improve results, and in fact can increase the foe. In addition, minimizing the MSE does not necessarily minimize the foe. We proposed here to use angular interpolation of q-space rather than SH interpolation, since it is simpler and uses the information in the image to perform intensity correction in the DW-MRIs, while SH uses a smooth basis that might lead to loss of information. Hence, AI should provide a better intensity correction than SH. Since AI is currently based on the intensity of the images, future work should address the possibility of extending the idea of AI to use a different feature.

Finally, no significant difference was found between the registration errors using the 3T dataset and the 7T datasets. This indicates that the algorithms studied can be applied to clinical data, where 3T datasets are more common than 7T data.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported in part by NIH grants R01 EB008432, R01 EB008645, R01 EB007813, S10 RR026783, P30 NS057091, P30 NS076408, P41 RR008079 (NCRR), P41 EB015894 (NIBIB), the Human Connectome Project (U54 MH091657) from the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research and the W.M. Keck Foundation.

Footnotes

1. ^FMRIB’s linear image registration tool (http://fsl.fmrib.ox.ac.uk/fsl/flirt/).
2. ^FMRIB’s non-linear image registration tool (http://www.fmrib.ox.ac.uk/fsl/fnirt/index.html).
3. ^http://www.fmrib.ox.ac.uk/fsl.
4. ^http://www.nitrc.org/projects/dtitk.
5. ^http://www-sop.inria.fr/asclepios/software/MedINRIA.
6. ^The most recent version of DTI-TK uses the PPD reorientation algorithm.
7. ^Currently available at www.nitrc.org/projects/dwiregistration/.
8. ^https://www.fmrib.ox.ac.uk/phpwiki/index.php/FSL%20Plugins.
9. ^Usually nearest neighbor, trilinear, spline, or windowed sinc interpolation.
10. ^FNIRT accepts 4D images, but works only with the first volume it encounters and disregards the remaining volumes.
11. ^FNIRT minimization algorithm is a gradient descent approach, hence, its extension to other non analytically derivable cost functions is not straightforward.
12. ^http://www.picsl.upenn.edu/ANTS/
13. ^http://nifti.nimh.nih.gov/nifti-1.
14. ^We discarded 8 regions (see supplementary material) from the 48 white matter labels in the JHU atlas, located in the lower brain areas and that were not always fully visible in our datasets.
15. ^Available within the FSL library.
16. ^http://www.itk.org/
17. ^http://www.itksnap.org/pmwiki/pmwiki.php?n=Convert3D.Documentation
18. ^http://med.inria.fr/
19. ^http://www.mathworks.com/matlabcentral/fileexchange/8797.

References

Alexander, D., Pierpaoli, C. J. G., Basser, P. J., and Gee, J. C. (2001). Spatial transformations of diffusion tensor magnetic resonance images. IEEE Trans. Med. Imaging 20, 1131–1139.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Andersson, J. L. R., Jenkinson, M., and Smith, S. (2007). Non-linear registration aka spatial normalisation. Technical Report FMRIB Technical Report TR07JA2, Oxford: FMRIB Centre.

Arsigny, V., Fillard, P., Pennec, X., and Ayache, N. (2005). Fast and simple calculus on tensors in the log-Euclidean framework. Med. Image Comput. Comput. Assist. Interv. 3749, 115–122.

Barmpoutis, A., Vemuri, B. C., and Forder, J. R. (2007). Registration of high angular resolution diffusion MRI images using 4th order tensors. Med. Image Comput. Comput. Assist. Interv. 10, 908–915.

Pubmed Abstract | Pubmed Full Text

Basser, P. J., and Jones, D. K. (2002). Diffusion-tensor MRI: theory, experimental design and data analysis – a technical review. NMR. Biomed. 15, 456–467.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Basser, P. J., Mattielo, J., and Le Bihan, D. (1994). MR diffusion tensor spectroscopy and imaging. J. Biophys. 66, 259–267.

CrossRef Full Text

Basser, P. J., and Pierpaoli, C. (1996). Microstructural and physiological features of tissues elucidated by quantitative diffusion tensor MRI. J. Magn. Reson. 111, 209–219.

CrossRef Full Text

Behrens, T. E. J., Woolrich, M., Jenkinson, M., Johansen-Berg, H., Nunes, R., Clare, S., et al. (2003). Characterization and propagation of uncertainty in diffusion-weighted MR imaging. Magn. Reson. Med. 50, 1077–1088.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bouix, S., Rathi, Y., and Sabuncu, M. (2010). Building an average population HARDI atlas, in MICCAI Workshop on Computational Diffusion MRI, Beijing, 84–91.

Chao, T.-C., Choua, M.-C., Yang, P., Chunga, H.-W., and Wub, M.-T. (2009). Effects of interpolation methods in spatial normalization of diffusion tensor imaging data on group comparison of fractional anisotropy. Magn. Reson. Imaging 27, 681–690.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cheng, G., Vemuri, B. C., Carney, P. R., and Mareci, T. H. (2009). Non-rigid registration of high angular resolution diffusion images represented by Gaussian mixture fields. Med. Image Comput. Comput. Assist. Interv. 5761, 190–197.

Pubmed Abstract | Pubmed Full Text

Descoteaux, M., Angelino, E., Fitzgibbons, S., and Deriche, R. (2006). “Apparent diffusion profile estimation from high angular resolution diffusion images,” in Proceedings of the SPIE Image Processing: Medical Imaging, Vol. 6144 (San Diego: SPIE), D1–D12.

Dhollander, T., Hecke, W. V., Maes, F., Sunaert, S., and Suetens, P. (2010). Spatial transformations of high angular resolution diffusion imaging data in q-space. MICCAI Workshop on Computational Diffusion MRI, Beijin, 73–83.

Du, J., Goh, A., and Qiu, A. (2011). Diffeomorphic metric mapping of high angular resolution diffusion imaging based on Riemannian structure of orientation distribution functions. IEEE Trans. Med. Imaging 31, 1021–1033.

Pubmed Abstract | Pubmed Full Text

Filler, A. (2009). Magnetic resonance neurography and diffusion tensor imaging: origins, history, and clinical impact of the first 50,000 cases with an assessment of efficacy and utility in a prospective 5000-patient study group. Neurosurgery 65, A29–A43.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fletcher, R. (1987). Practical Methods of Optimization, 2nd Edn. Hoboken: John Wiley & Sons.

Geng, X., Ross, T. J., Gu, H., Shin, W., Zhan, W., Chao, Y.-P., et al. (2011). Diffeomorphic image registration of diffusion MRI using spherical harmonics. IEEE Trans. Med. Imaging 30, 747–758.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hasan, K. M., Parker, D. L., and Alexander, A. L. (2001). Comparison of optimization procedures for diffusion-tensor encoding directions. Magn. Reson. Imaging 13, 769–780.

CrossRef Full Text

Hess, C. P., and Mukherjee, P. (2007). Visualizing white matter pathways in the living human brain: diffusion tensor imaging and beyond. Neuroimaging Clin. N. Am. 17, 407–426.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Higham, N. (1986). Computing the polar decomposition-with applications. SIAM J. Sci. Stat. Comput. 7, 1160–1174.

CrossRef Full Text

Hua, K., Zhang, J., Wakana, H. J., Li, X., Reich, D. S., Calabresi, P. A., et al. (2008). Tract probability maps in stereotaxic spaces: analyses of white matter anatomy and tract-specific quantification. Neuroimage 39, 336–347.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Irfanoglu, M. O., Machiraju, R., Sammet, S., Pierpaoli, C., and Knopp, M. (2008). Automatic deformable diffusion tensor registration for fiber population analysis. Med. Image Comput. Comput. Assist. Interv. 2, 1014–1022.

Jenkinson, M., Bannister, P., Brady, M., and Smidth, S. (2002). Improved optimization for the robust and accurate linear registration and motion correction of brain images. Neuroimage 17, 825–841.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jones, D. K., Horsfield, M. A., and Simmons, A. (1999). Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magn. Reson. Med. 42, 515–525.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kamath, A., Aganj, I., Xu, J., Yacoub, E., Ugurbil, K., Sapiro, G., et al. (2012). Generalized constant solid angle odf and optimal acquisition protocol for fiber orientation mapping, in MICCAI Workshop on Computational Diffusion MRI. Nice.

Kannala, J., Rahtu, E., Heikkil, J., and Salo, M. (2005). Fast and simple calculus on tensors in the log-Euclidean framework. Lect. Notes Comput. Sci. 3540, 224–234.

CrossRef Full Text

Klein, A., Andersson, J. L. R., Ardekani, B., Ashburner, J., Avants, B., and Chiang, M. (2009). Evaluation of 14 nonlinear deformation algorithms applied to human brain MRI registration. Neuroimage 46, 786–802.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Koay, C. G., Ozarslan, E., and Pierpaoli, C. (2009). Probabilistic identification and estimation of noise (PIESNO): a self-consistent approach and its applications in MRI. J. Magn. Reson. 199, 94–103.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Le Bihan, D., Poupon, C., Amadon, A., and Lethimonnier, F. (2006). Artifacts and pitfalss in diffusion MRI. J. Magn. Reson. Imaging 24, 478–488.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lenglet, C., Abosch, A., Yacoub, E., De Martino, F., Sapiro, G., and Harel, N. (2012). Comprehensive in vivo mapping of the human basal ganglia and thalamic connectome in individuals using 7T MRI. PLoS ONE 7:e29153. doi:10.1371/journal.pone.0029153

CrossRef Full Text

Noblet, V., Heinrich, C., Heitz, F., and Armspach, J.-P. (2006). Retrospective evaluation of a topology preserving non-rigid registration method. Med. Image Anal. 10, 366–384.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ourselin, S., Roche, A., Prima, S., and Ayache, N. (2000). “Block matching: a general framework to improve robustness of rigid registration of medical images,” in MICCAI Workshop on Computational Diffusion MRI, Vol. 1935, Pittsburgh 557–566.

Pierpaoli, C., and Basser, P. J. (1996). Toward a quantitative assessment of diffusion anisotropy. Magn. Reson. Med. 36, 893–906.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Raffelt, D., Tournier, J.-D., Fripp, J., Crozier, S., Connelly, A., and Salvado, O. (2011). Symmetric diffeomorphic registration of fibre orientation distributions. Neuroimage 56, 1171–1180.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Skare, S., Li, T.-Q., Nordell, B., and Ingvar, M. (2000). Noise considerations in the determination of diffusion tensor anisotropy. Magn. Reson. Imaging 18, 659–669.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Stejskal, E., and Tanner, J. (1965). Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient. J. Chem. Phys. 42, 288–292.

CrossRef Full Text

Stobbe, R. W., and Beaulieu, C. (2011). Assessment of averaging spatially correlated noise for 3-d radial imaging. IEEE Trans. Med. Imaging 30, 1381–1390.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tao, X., and Miller, J. V. (2006). A method for registering diffusion weighted magnetic resonance images. Med. Image Comput. Comput. Assist. Interv. 9, 594–602.

Pubmed Abstract | Pubmed Full Text

Thomas Yeo, B., Vercauteren, T., Fillard, P., Peyrat, J.-M., Pennec, X., Golland, P., et al. (2009). DT-REFinD: diffusion tensor registration with exact finite-strain differential. IEEE Trans. Med. Imaging 28, 1914–1928.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Verma, R., and Bloy, L. (2010). “Demons registration of high angular resolution diffusion images,” in Biomedical Imaging: From Nano to Macro (Philadelphia: IEEE), 1013–1016.

Wang, Y., Gupta, A., Liu, Z., Zhang, H., Escolar, M. L., Gilmore, J. H., et al. (2011). DTI registration in atlas based fiber analysis of infantile Krabbe disease. Neuroimage 55, 1577–1586.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Yap, P.-T., Chen, Y., An, H., Gilmore, J. H., Lin, W., and Shen, D. (2010). “Non-parametric deformable registration of high angular resolution diffusion data using diffusion profile statistics,” in ISMRM 18 annual meeting, Stockholm.

Zhan, L., Mueller, B. A., Jahanshad, N., Jin, Y., Lenglet, C., Yacoub, E., et al. (2013). Magnetic resonance field strength effects on diffusion measures and brain connectivity networks. Brain Connect. 3, 72–86.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Zhang, H., Yushkevich, P. A., Alexander, D., and Gee, J. C. (2006). Deformable registration of diffusion tensor MR images with explicit orientation optimization. Med. Image Anal. 10, 764–785.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Appendix

Generation of Synthetic Deformations

The following matrices R, S, and D provide random three-dimensional rotations, skewness, and change of scale, respectively

$\begin{array}{llll}R& =\left[\begin{array}{cccc}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill cos\alpha \hfill & \hfill -sin\alpha \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill sin\alpha \hfill & \hfill cos\alpha \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill \hfill \end{array}\right]\left[\begin{array}{cccc}\hfill cos\beta \hfill & \hfill 0\hfill & \hfill sin\beta \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -sin\beta \hfill & \hfill 0\hfill & \hfill cos\beta \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill \hfill \end{array}\right]& & \\ & \times \left[\begin{array}{cccc}\hfill cos\gamma \hfill & \hfill -sin\gamma \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill sin\gamma \hfill & \hfill cos\gamma \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill \hfill \end{array}\right],& & \\ S& =\left[\begin{array}{cccc}\hfill 1\hfill & \hfill a\hfill & \hfill b\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill c\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill \hfill \end{array}\right],D=\left[\begin{array}{cccc}\hfill \delta \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \delta \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \delta \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill \hfill \end{array}\right],& & \end{array}$

where α, β, γ are uniform random rotation angles along the x, y, and z axis respectively in the [−15°, 15°] range; a, b, c are uniform random numbers in the −0.125, 0.125] range; and δ is a scale uniform random number in the [0.875, 1.125] range.

Random nonlinear warps

The following equation defines the displacement field u(x) along the x, y, and z axis,

$u\left(x\right)=\left\{\begin{array}{c}{a}_{x}sin\left(\frac{w\pi x}{n_x}\right)\hfill \\ a_{y}sin\left(\frac{w\pi y}{n_y}\right)\hfill \\ a_{z}sin\left(\frac{w\pi z}{n_z}\right)\hfill \\ \hfill \end{array}\right.s.t.\left(\frac{\left|a_x\right|}{n_x}+\frac{\left|a_y\right|}{n_y}+\frac{\left|a_y\right|}{n_y}\right)\le \frac{1}{w\pi },$

where n_x, n_y, n_z are the number of voxels along the x, y, z axis, respectively, and w, a_x, a_y, a_z are uniform random variables satisfying the indicated conditions.

Implementation Details

Effect of Interpolation on FA

The effect of spatial interpolation on DT-MRIs has been recently studied (Chao et al., 2009), indicating that the measured Fractional Anisotropy (FA) can be significantly reduced if linear interpolation is used. Basically, linear interpolation acts as a low-pass filter with equal weights that averages the tensors, leading to tensor swelling (Arsigny et al., 2005) and a reduction of the estimated FA values. We have not found an equivalent study for the effect of spatial interpolation of DW-MRIs on the estimated FAs, but we were able to see the same effect on the estimated FA values from the registered DW-MRIs when tri-linear interpolation is used, which is the default interpolation in FSL. In addition, we also found that using the recommended values in Tao and Miller (2006) for the standard deviation in (3) and (19), AI can also reduce the calculated FA values. Notice that AI can also be seen as a Gaussian filter, with zero mean and standard deviation σ defining the weights of the interpolating window.

We address this problem by using a Blackman windowed sinc interpolation filter (provided by FSL) for the spatial interpolation and reducing the standard deviation of the Gaussian coefficients for AI, in such a way that only the volumes closest to the specified gradient are accounted for. Only the closest gradients in the image should indeed be used to interpolate along a given gradient direction, otherwise, volumes with different gradients (and very different intensity responses) could affect the intensities computed for a new gradient direction.

Figure A1A shows the estimated FA for the registered test image using trilinear interpolation and the smoothing parameter σ computed from the average of the geodesic distances between neighboring gradients, as indicated in Tao and Miller (2006), Figure A1B shows the estimated FA for the registered test image using windowed sinc interpolation and a much lower σ (0.005) value, and Figure A1C compares the histogram of FA values for the reference and registered test images using the cited interpolation methods. As can be seen from these images, using sinc interpolation and a lower σ the FA values gets closer to those from the original image.

FIGURE A1

Figure A1. (A) FA of the registered image using trilinear interpolation and large σ. (B) FA registered using sinc interpolation and small σ, (C) Histograms of the FAs.

Registration

We extended FSL to register DW-MR images, trying always to minimize the impact on the library, by using as much as possible the algorithms and classes already implemented there. Given that working with the whole DW-MR image with AI increases the computation cost, especially for nonlinear registration, we parallelized the linear and nonlinear registration algorithms using OpenMP¹ and modify the AI algorithm to reduce computational costs. As can be seen from (3) and (18), we need to combine linearly hundreds of volumes, which is a costly operation. We reduced significantly the running time by fixing a low σ and using a kd-tree to search for the 16 closest gradients to a given target gradient orientation. The value of σ and the number of closest neighbors used were found experimentally as those values that produces a negligible difference with the original algorithm that uses all the volumes in the DW-MR image. The kd-tree is constructed only once, using all the gradients in the image, and it is used later in all the computations involving AI, which provides a reduction in the time complexity of the algorithm from O(mn) to O(n log(m)), where m is the number of volumes in the image and n the number of voxels.

We investigated also the effect of the number of volumes, for nonlinear registration, and found that using just m = 65 volumes were enough to register the whole DW-MRIs obtaining results very close to those using all the volumes in the original images (128 volumes), which provides a further reduction in the running time of nonlinear registration. Currently, linear registration can take 15 min using 4 processors and all the volumes in the image, while nonlinear registration can take 30–45 min using 8 processors and 65 volumes.

In MedINRIA, we used the b0 volume to obtain a transformation matrix and then perform linear and nonlinear registration of the estimated DT-MRIs. DTI-TK provides registration algorithms to register linearly and nonlinearly (Zhang et al., 2006) DT-MRIs as well as tools to convert from the tensor estimates from FSL-dtifit to the NIFTI² format that DTI-TK adheres to.

TABLE A1

Table A1. John Hopkins University (JHU) white matter selected labels.

Footnotes

1. ^http://openmp.org/wp.
2. ^http://nifti.nimh.nih.gov/