Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Med Image Comput Comput Assist Interv. Author manuscript; available in PMC 2010 May 4.
Published in final edited form as:
Med Image Comput Comput Assist Interv. 2009; 5761: 190–197.
doi:  10.1007/978-3-642-04268-3_24
PMCID: PMC2863327

Non-rigid Registration of High Angular Resolution Diffusion Images Represented by Gaussian Mixture Fields


In this paper, we present a novel algorithm for non-rigidly registering two high angular resolution diffusion weighted MRIs (HARDI), each represented by a Gaussian mixture field (GMF). We model the non-rigid warp by a thin-plate spline and formulate the registration problem as the minimization of the L2 distance between the two given GMFs. The key mathematical contributions of this work are, (i) a closed form expression for the derivatives of this objective function with respect to the parameters of the registration and (ii) a novel and simpler re-orientation scheme based on an extension to the “Preservation of Principle Directions” technique. We present results of our algorithm’s performance on several synthetic and real HARDI data sets.

1 Introduction

Diffusion-Weighted MR Imaging (DWMRI) is a unique non-invasive technique that makes the MR signal sensitive to the water molecule diffusion in the imaged tissue and infer its structure in vivo. When using DWMRI in longitudinal studies, multi-subject studies, and such tasks, non-rigid registration is essential to quantify the differences between the acquired data sets.

The diffusivity function at a lattice point in a DWMRI when approximated by a rank-2 symmetric positive definite tensor leads to what is popularly dubbed a diffusion tensor image (DTI)[1]. Several DTI registration methods have been reported in literature to date. Alexander et al. [2] developed a registration technique for DTI and observed that a re-orientation of the tensors was necessary as part of the registration of the DWMRI data sets. Some existing methods [35] used one or more rotation invariant characteristics of diffusion tensors to perform the registration and thereby avoiding the re-orientation task at each iteration of the registration algorithm. A variational model for the diffeomorphic DTI registration was proposed in [6] along with an analytic gradient of the Preservation of Principle Direction (PPD) reorientation strategy, which is very complex and cumbersome requiring substantial compute time. In [7] the analytic gradient of the Finite Strain (FS) strategy was developed for diffeomorphic DTI registration. None of these methods are directly applicable to the case of registering GMFs.

Rank-2 SPD tensors are a good approximation for the single fiber geometry. However, this model is known to fail at the locations with complex tissue geometry such as fiber crossings. Some higher order models [812] have been reported in literature to tackle this problem. Research on registration of the resulting representations from these methods has not been addressed to date with the exception of work in [13] wherein, a piece-wise affine approximation to the nonrigid registration of fourth order tensor fields was developed. Most recently, a fluid-flow based model of the non-rigid transformation was introduced in [14] for registering two given DTI data sets. Authors claimed that their work was easily extendable to registration of higher order representations derived from HARDI data but no such extension has been reported to date.

In this paper, we present a novel non-rigid registration method for HARDI datasets represented by GMFs generated by using the algorithm described in [9]. The non-rigid registration between the GMFs in this paper is represented by a thin plate spline (TPS) and the dissimilarity between the GMFs being registered is expressed by the L2 distance between Gaussian mixtures at corresponding points between the image lattices, and can be viewed as an extension of the work in [15]. The re-orientation operation is handled by a novel and sig-nificantly improved extension (developed here) of the PPD strategy in [2]. The key contributions of the work reported here are, (i) it is the first attempt at the non-rigid registration of HARDI data sets represented by GMFs. (ii) An objective function involving the L2 distance between Gaussian mixtures is used that leads to a closed form expression for the distance. (iii) A significant extension of the PPD strategy is to handle re-orientation of Gaussian mixtures involving a derivation of the objective function and its analytic gradient both of which are in a much simpler form than those derived for DTI re-orientation in [6]. Experiments on synthetic and real data along with comparisons demonstrate the performance of our algorithm.

The rest of the paper is organized as follows: In section 2, we present registration framework including a derivation of the analytic form of the L2 distance between two Gaussian mixtures and our re-orientation strategy. We then present the TPS formulas followed by the derivation of a closed form expression of the cost function derivatives w.r.t. the deformation parameters. In section 3, we present synthetic and real data experiments along with comparisons to other methods and conclude in section 4.

2 Non-rigid Registration of Gaussian Mixture Fields

There are several ways to represent the HARDI data and the choice depends on the end goal. We choose to employ the diffusion propagator pdf P(r, t) – where r is the displacement vector from a given location at time t – as our representation and compute this field from the given HARDI data. The reason for this choice is the generality and accuracy of the continuous mixture model used here & described in [8, 9].

We now present our registration framework for GMFs obtained using the method presented in [9]. Here a GMF is a map from R3 to the space of smooth functions C. Given the target GMF I and the source GMF J, registration can be described as a process of estimating the transform T: R3R3 from the coordinates of I to the coordinates of J, which best aligns the two GMFs. This can be done by minimizing the dissimilarity between the target GMF and the transformed source GMF with respect to T. The energy function for this optimization process can be written as


where IT(x) and J(T(x)) for a fixed position x [set membership] R3 are two zero mean Gaussian mixtures(GMs), dist(ο, ο) denotes the dissimilarity measure between the two GMs and IT denotes the GMF I after re-orientation. The details of equation 1 are described in the remaining of this section.

L2 Distance between Zero Mean Gaussian Mixtures

In this paper we use the L2 distance as a measure of the dissimilarity between zero mean GMs. Let f(r)=i=1MηiG(r;0,i) and g(r)=j=1NρjG(r;0,Γj) be two GM density functions, where r [set membership] R3 is the displacement vector and ηi, ρj denote the mixture weights of the corresponding Gaussian components G(r; 0, Σi) and G(r; 0, Γj) with covariance matrices Σi and Γj respectively. The L2 distance between f and g can be written as a quadratic function of the mixture weights


where η = (η1, …, ηM)t and ρ = (ρ1, …, ρN)t are vectors representing the mixture weights, and AM×M, BN×N and CM×N are the matrices generated by the Gaussian components, with their elements to be Aiii2 = ((2π)3 det (Σi1 + Σi2))−1/2, Bjij2 = ((2π)3 det (Γj1 + Γj2))−1/2 and Cij = ((2π)3 det (Σi + Γj))−1/2

If a 3D covariance matrix Σ has eigen values λ1λ2 = λ3 (cylindrical symmetry), its rank one decomposition can be written as Σ = (λ1λ2)uut+λ2I, where u is the principal eigen vector of Σ. Given another cylindrically symmetric covariance matrix Γ = (ξ1ξ2)vvt + ξ2I, we can derive the following equation after applying some matrix algebra,


where α = (ξ1 + λ2) (λ1 + ξ2) (λ2 + ξ2), β = (λ1λ2) (ξ1ξ2) (λ2 + ξ2). It is reasonable to assume that all the Gaussian components in our GMFs are cylindrically symmetric, as the fibers have approximately cylindrical geometry. Thus, equation 3 can be used to compute the matrices A, B and C in equation 2. Using equation 3, reorientation can be easily achieved and the registration function along with its gradient can be simplified.


As a GMF is not rotation invariant, we need to perform a re-orientation during the registration process. In this context, DTI reorientation methods would be useful to consider and generalize. There are three DTI reorientation strategies reported in literature, namely, (i) Preservation of Principle Directions(PPD) [2], (ii) Finite Strain(FS) [2] and (iii) Re-transformation [13].

In the PPD, the principle direction v of the 2rd order tensor D is transformed to v = Fυ/|Fυ| where F is the Jacobian of the transform T. In this paper, we extend the PPD strategy to the mixture of Gaussians model in a direct way by applying it to the covariance matrices of each Gaussian component. With the cylindrical geometry assumption on the nerve fibers, let the source GMF J(x)=i=1Mηi(x)G(r;0,(λ1iλ2i)uiuit+λ2iI) and target GMF I(x)=j=1Nρj(x)G(r;0,(ξ1jξ2j)vivjt+ξ2jI). Here we can assume λi., ξj., ui vj to be constant across the lattice for simplicity. The re-orientation is applied to the target GMF I by applying the Jacobian F to each of its covariance matrices. Thus, IT(x)=j=1Nxρj(x)G(r;0,(ξ1jξ2j)vjFFtvjt/||Fvi||2+ξ2jI). Using Equation 3, the energy function can be written as


where CTij = k(ιij)and BTj1j2 = k(τjij2), with scalar ιij=(uitFvj)2/||Fvj||2,τjij2=(vj1tFtFvj2)2/(||Fvj1||2||Fvj2||2), and k(t) = ((2π)3(αβt))(−1/2), α and β are scalar constants defined using the eigen values as in equation 3.

2.1 The TPS-based Non-rigid Registration

We are now ready to present the estimation of the non-rigid registration T – between the GMFs I ad J – represented by a TPS [15]. Let the set of control grid points for computing the registration be {x1, x2,xn}, the transform is given byT(x) = WH(x) + Ax + t, where Ax + t is the affine part, and WH(x) is the non-rigid part. H(x) = (H1(x), H2(x).., Hn(x))t is a group of non-linear kernel functions with Hi(x) = ri(x)2 log(ri(x)) in the case of 2D and Hi(x) = ri(x)3 in the case of 3D, where ri(x) = ||x–xi||. This non-rigid transform is regularized by minimizing the bending energy of the TPS given by trace(WKWt). Here K is the dissimilarity matrix with elements Kij = ri(xj). To perform re-orientation, the Jacobian matrix F at every image grid point also needs to be computed, and is given by F = WJ(H(x)) + A, where J(H(x)) = ([nabla]H1t, [nabla]H2t,[nabla]Hnt)t|x is the Jacobian of H(x). The non-rigid registration here involves minimization of the energy function w.r.t. the parameters A, t and W. The affine parameters A and t can be solved for via an affine registration technique prior to solving for the non-rigid component W.

2.2 Analytic Derivative of the Objective Function

The derivative of the objective function is straight forward to compute by applying the chain rule. The partial derivative of the energy function (without the bending energy) w.r.t. the elements wkc in the matrix W (recall that W is weight matrix of the non-rigid part of the TPS) can be written as:

[partial differential]E[partial differential]wkc=R32[partial differential]Tt[partial differential]wkc[nabla]ηtAη2[partial differential]Tt[partial differential]wkc[nabla]ηtCTρ+ρt[partial differential]BT[partial differential]wkcρ2ηt[partial differential]CT[partial differential]wkcρdx

where the derivative of the TPS function [partial differential]T[partial differential]wkc is a 3-by-1 vector, with its k-th element being the c-th kernel function value at x, namely Hc(x), and two other elements being zeros; the derivative of CT and BT can be computed as:

ρt[partial differential]BT[partial differential]wkcρ=Tr([partial differential]Ft[partial differential]wkcj1,j2ρj1ρj2[partial differential]Bj1j2[partial differential]F);ηt[partial differential]CT[partial differential]wkcρ=Tr([partial differential]Ft[partial differential]wkci,jηiρj[partial differential]Cij[partial differential]F)

where [partial differential]Bj1j2[partial differential]F=2kτj1j2(F(vjivj2t+vj2vj1t)/(vj1tFtFvj2)Fvj1vj1t/||Fvj2||2Fvj2vj2t/||Fvj2||2),[partial differential]Cij[partial differential]F=2kιij(uivjt/uitFvjFvjvjt/||Fvj||2), and [partial differential]F[partial differential]wkc=ykc, where ykc is a 3-by-3 matrix with each element given by ynmkc=[partial differential]Hc[partial differential]xmδ(nk) and δ· being the discrete delta function.

Since we have the gradient of the energy function in the analytic form, any of the gradient based optimization strategies in the literature can be applied to optimize the registration cost function. In this paper we use the efficient limited-memory Broyden-Fletcher-Goldfarb quasi-Newton technique [16]. Note that, as claimed in [9], the weight vectors in the Gaussian mixtures are sparse. Thus, we can just compute the elements of the matrices corresponding to nonzero weights.

3 Experiments

We now report the experimental results obtained on synthetic and real HARDI datasets. In the synthetic data experiment, we registered two images with crossing fiber bundles. Firstly, a 2D synthetic image with two crossing fiber bundles was manually generated. Then, 20 randomly deformed images were generated from the crossing bundles by applying a b-spline based non-rigid deformation and a PPD based re-orientation scheme described here. The method described in [17] was used to generate the simulated MR signals from the fiber bundles. Rician noise was added to simulate data at 4 different noise levels with SNR = 50, 20, 10 and 5. The method in [9] was used to generate the GMF from the MR signals with 46 Gaussian components at each voxel. After the data generation, we registered each of the randomly deformed source images to the fixed target image separately. To evaluate the registration, the resulting deformation obtained from the registration was applied to the noise free source image, and then the relative dissimilarity between the deformed source and target image was computed as the relative error in registration given by,

Err=x[set membership]RS2(Ix(r)Jx(r))2drx[set membership]RS2(Ix(r))2dr

where R is a user selected region to evaluate the errors; I; J denotes the target and deformed source image; Ix(r) and Jx(r) are the displacement probability profiles (represented in spherical harmonic coefficients) at voxel x. Note that this dissimilarity measure is different from the one used in our registration algorithm. In this paper, two different regions are used for error computation: (1) The whole fiber bundle region and (2) The region containing just the fiber crossings. Also, a GA (generalized Anisotropy) based registration algorithm (using a sum of squared differences cost) and the DTI based registration algorithm in [18] were tested using the same data set for comparison purposes. The data sets and the results are displayed in Figure 1. Figure (f) and (g) show that our method has much lower mean and standard deviation of registration errors for both the chosen ROIs for first three noise levels, and has a comparable error to other two methods at SNR5 (which is very high amount of noise).

Fig. 1
Experimental results on synthetic dataset. (a) is the source, (b) is the target and (c) is the deformed source image. Figure (d) and (e) are the mean and standard deviation of the error for the 20 registrations from all the three methods at different ...

In the real data experiment, we used rat brain DWMRI. For data acquisition, in each scan 27 diffusion gradients were used (with Δ and δ set to be 17.6ms and 4.8ms) with a b = 800s/mm2 for 21 of the 27 gradients and the remaining with a b = 100s/mm2. The image size was 100 × 100 × 12, with each voxel of size, 0.3mm × 0.3mm × 0.9mm. In the first experiment, we show results from two distinct rats. For the second experiment, we depict results from a single rat, taken at different time points after implanting electrodes used to stimulate injury that results in the development of epilepsy. Consecutive scans were used in this registration experiment to capture structural changes via the non-rigid registration.

For the first real data experiment, the data and the registration result are shown in a checker-board view in Figure 2. (In the checker-board view alternate blocks come from the source and target images respectively). From the figure its evident that our registration algorithm successfully aligned the two given HARDI data sets. In sub-figure (e) and (f), there is a fiber bundle from the corpus callosum in the ROI, which contains most of the probability profile in green and brown (the stream tubes are a result of the tracking results using a simple vector field integration based tracker). In (e) (before registration) this fiber bundle is shown disrupted and depicted as two distinct disconnected tracts in the checker-board view, and in (f) (after registration) it is shown as one coherent connected tract indicating the accuracy of the registration visually.

Fig. 2
Real data experiments:(a) and (b) are the source and targets S0 images, (c) checker-board view of source and target images, (d) checker-board view of the deformed source and target images. (e) and (f) Displacement probabilities in the ROI shown in (c) ...

For quantitative evaluation, we compared the dissimilarity of the target and source image pair before and after registration using equation 5. Then, we applied our algorithm to consecutive DWMR time scans of the same rat, and plotted the dissimilarity for the image pairs before and after registration. Plates (g) and (h) in Figure 2 are the results for two different rats one with 9 consecutive scans and the other with 6. From the figure, it is apparent that the dissimilarity between them has decreased significantly after registration. This indicates the effectiveness of our registration method on real data.

4 Conclusion

In this paper we presented a novel non-rigid registration algorithm for HARDI data represented by GMFs. The non-rigid transformation is modeled using a TPS function. The key contributions of this work are, (i) it is the first report on non-rigid registration of GMFs representing HARDI data sets. (ii) A non-trivial derivation leading to a closed form expression for the gradient of the registration was presented. (iii) A novel extension of the PPD strategy of re-orientation and a novel derivation of the reorientation that is simpler than the previously reported closed form for PPD-based re-orientation of DTI. Both real and synthetic data experiments along with comparisons were presented depicting the effectiveness of the presented registration algorithm.


This research was funded by the NIH grants EB007082 to BCV, & EB004752 to PC & TM.


1. Basser PJ, Mattiello J, LeBihan D. Estimation of the effective self-diffusion tensor from the NMR spin echo. J Magn Reson B. 1994;103:247–254. [PubMed]
2. Alexander DC, Pierpaoli C, Basser PJ, Gee JC. Spatial transformations of diffusion tensor magnetic resonance images. IEEE TMI. 2001;20:1131–1139. [PubMed]
3. Jones DK, Griffin LD, Alexander DC, Catani M, Horsfield MA, Howard R, Williams SCR. Spatial normalization and averaging of diffusion tensor MRI data sets. NeuroImage. 2002;17:592–617. [PubMed]
4. Park HJ, Kubicki M, Shenton ME, Guimond A, McCarley R, Maier SE, Kikinis R, Jolesz FA, Westin CF. Spatial normalization of diffusion tensor MRI using multiple channels. NeuroImage. 2003;20:1995–2009. [PMC free article] [PubMed]
5. Ziyan U, Sabuncu MR, O’Donnell LJ, Westin C. Nonlinear registration of diffusion MR images based on fiber bundles. MICCAI. 2007:351–358. [PubMed]
6. Cao Y, Miller M, Mori SRL, Winslow LY. Diffeomorphic matching of diffusion tensor images. Proc CVPRW’06. 2006:67–67. [PMC free article] [PubMed]
7. Yeo B, Vercauteren T, Fillard P, Pennec X, Golland P, Ayache N, Clatz O. DTI registration with exact finite-strain differential. ISBI. 2008:700–703. [PubMed]
8. Jian B, Vemuri BC, Ozarslan E, Carney PR, Mareci TH. A novel tensor distribution model for the DWMR signal. NeuroImage. 2007;37:164–176. [PMC free article] [PubMed]
9. Jian B, Vemuri BC. A unified computational framework for deconvolution to reconstruct multiple fibers from DWMRI. IEEE TMI. 2007;26:1464–1471. [PMC free article] [PubMed]
10. Barmpoutis A, Jian B, Vemuri BC, Shepherd TM. Symmetric positive 4th order tensors and their estimation from DWMRI. IPMI. 2007:308–319. [PMC free article] [PubMed]
11. Özarslan E, Mareci TH. Generalized diffusion tensor imaging and analytical relationships between diffusion tensor imaging and high angular resolution diffusion imaging. Magn Reson Med. 2003;50:955–965. [PubMed]
12. Descoteaux M, Angelino E, Fitzgibbons S, Deriche R. Apparent diffusion coefficients from high angular resolution diffusion imaging: Estimation and applications. Magn Reson Med. 2006;56:395–410. [PubMed]
13. Barmpoutis A, Vemuri BC, Forder JR. Registration of high angular resolution diffusion MRI images using 4 th order tensors. MICCAI. 2007:908–915. [PMC free article] [PubMed]
14. Chiang MC, Leow AD, Klunder AD, Dutton RA, Barysheva M, Rose SE, McMahon KL, de Zubicaray GI, Toga AW, Thompson PM. Fluid registration of DT images using information theory. IEEE TMI. 2008;27:442–456. [PMC free article] [PubMed]
15. Jian B, Vemuri BC. A robust algorithm for point set registration using mixture of gaussian. ICCV. 2005:1246–1251. [PMC free article] [PubMed]
16. Nocedal J, Wright S. Numerical Optimization. Springer; 1999.
17. Soderman O, Jonsson B. Restricted diffusion in cylindrical geometry. J Magn Reson A. 1995;117:94–97.
18. Zhang H, Yushkevich PA, Gee JC. Deformable registration of diffusion tensor MR images with explicit orientation optimization. MICCAI. 2005:172–179. [PubMed]