Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2863327

Formats

Article sections

- Abstract
- 1 Introduction
- 2 Non-rigid Registration of Gaussian Mixture Fields
- 3 Experiments
- 4 Conclusion
- References

Authors

Related links

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_24PMCID: PMC2863327

NIHMSID: NIHMS140469

Guang Cheng: ude.lfu.esic@gnehcg; Baba C. Vemuri: ude.lfu.esic@irumev; Paul R. Carney: ude.lfu.sdep@rpenrac; Thomas H. Mareci: ude.lfu@iceramht

See other articles in PMC that cite the published article.

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.

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 [3–5] 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 [8–12] 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 *L*2 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.

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 *R*^{3} 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: R*^{3} → *R*^{3} 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

$$E(T)={\int}_{{R}^{3}}{\mathit{dist}}^{2}({I}_{T}(\mathbf{x}),J(T(\mathbf{x})))d\mathbf{x}$$

(1)

where *I _{T}*(

In this paper we use the *L*_{2} distance as a measure of the dissimilarity between zero mean GMs. Let
$f(\mathit{r})={\mathrm{\sum}}_{i=1}^{M}{\eta}_{i}G(\mathit{r};0,{\mathit{\sum}}_{i})$ and
$g(\mathit{r})={\mathrm{\sum}}_{j=1}^{N}{\rho}_{j}G(\mathit{r};0,{\mathit{\Gamma}}_{j})$ be two GM density functions, where **r** *R*^{3} is the displacement vector and *η _{i}, ρ_{j}* denote the mixture weights of the corresponding Gaussian components

$${\mathit{dist}}^{2}(f,g)={\int}_{{R}^{3}}{(f(\mathbf{r})-g(\mathbf{r}))}^{2}d\mathbf{r}={\mathit{\eta}}^{t}\mathit{A}\mathit{\eta}+{\mathit{\rho}}^{t}\mathit{B}\mathit{\rho}-{2\mathit{\eta}}^{t}\mathit{C}\mathit{\rho}$$

(2)

where ** η** = (

If a 3D covariance matrix **Σ** has eigen values *λ*_{1} ≥ *λ*_{2} = *λ*_{3} (cylindrical symmetry), its rank one decomposition can be written as **Σ** = (*λ*_{1} − *λ*_{2})**uu*** ^{t}*+

$$\mathit{det}(\mathbf{\sum}+\mathbf{\Gamma})=\alpha -\beta {({\mathbf{u}}^{\mathbf{t}}\mathbf{v})}^{2}$$

(3)

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 = *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(\mathbf{x})={\mathrm{\sum}}_{i=1}^{M}{\eta}_{i}(\mathbf{x})G(\mathit{r};0,({\lambda}_{1}^{i}-{\lambda}_{2}^{i}){\mathit{u}}_{i}{\mathit{u}}_{i}^{t}+{\lambda}_{2}^{i}\mathbf{I})$ and target GMF
$I(\mathbf{x})={\mathrm{\sum}}_{j=1}^{N}{\rho}_{j}(\mathbf{x})G(\mathit{r};0,({\xi}_{1}^{j}-{\xi}_{2}^{j}){\mathit{v}}_{i}{\mathit{v}}_{j}^{t}+{\xi}_{2}^{j}\mathbf{I})$. Here we can assume *λ ^{i}*.,

$$E(T)={\int}_{{R}^{3}}\mathit{\eta}{(\mathbf{T}(\mathbf{x}))}^{t}\mathit{A}\mathit{\eta}(\mathbf{T}(\mathbf{x}))+\mathit{\rho}{(\mathbf{x})}^{t}{\mathit{B}}_{\mathit{T}}\mathit{\rho}(\mathbf{x})-{2\mathit{\eta}}^{t}(\mathbf{T}(\mathbf{x})){\mathit{C}}_{\mathit{T}}\mathit{\rho}(\mathbf{x})d\mathbf{x}$$

(4)

where *C _{Tij}* =

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 {**x**_{1,} **x**_{2,} …**x*** _{n}*}, the transform is given by

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 *w _{kc}* in the matrix

$$\frac{E{wkc}_{}={\int}_{{R}^{3}}2\frac{{\mathbf{T}t}^{}{wkc}_{}{\mathit{\eta}}^{\mathit{t}}\mathit{A}\mathit{\eta}-2\frac{{\mathbf{T}t}^{}{wkc}_{}{\mathit{\eta}}^{\mathit{t}}{\mathit{C}}_{\mathit{T}}\mathit{\rho}+{\mathit{\rho}}^{t}\frac{{\mathit{B}\mathit{T}}_{}{wkc}_{}\mathit{\rho}-{2\mathit{\eta}}^{t}\frac{{\mathit{C}\mathit{T}}_{}{wkc}_{}\mathit{\rho}d\mathbf{x}}{}}{}}{}}{}}{}$$

where the derivative of the TPS function
${\scriptstyle \frac{\mathbf{T}{wkc}_{}}{}}$ is a 3-by-1 vector, with its *k*-th element being the *c*-th kernel function value at *x*, namely *H _{c}*(

$${\mathit{\rho}}^{t}\frac{{\mathit{B}\mathit{T}}_{}{wkc}_{}\mathit{\rho}=Tr(\frac{{\mathbf{F}t}^{}{wkc}_{}\sum _{{j}_{1},{j}_{2}}{\rho}_{{j}_{1}}{\rho}_{{j}_{2}}\frac{{B{j}_{1}{j}_{2}}_{}\mathbf{F}}{)}}{}}{}$$

where
${\scriptstyle \frac{{B{j}_{1}{j}_{2}}_{}\mathbf{F}}{=}}$, and
${\scriptstyle \frac{F{wkc}_{}}{=}}$, where *y** ^{kc}* is a 3-by-3 matrix with each element given by
${y}_{nm}^{kc}={\scriptstyle \frac{{Hc}_{}{xm}_{}}{\delta}}$ and

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.

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,

$$\mathit{Err}=\frac{{\sum}_{\mathbf{x}\mathbf{R}}{\sum}_{\mathbf{x}\mathbf{R}}}{}$$

(5)

where **R** is a user selected region to evaluate the errors; *I; J*^{′} denotes the target and deformed source image; *I*** _{x}**(

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.6*ms* and 4.8*ms*) with a *b* = 800*s/mm*^{2} for 21 of the 27 gradients and the remaining with a *b* = 100*s/mm*^{2}. The image size was 100 × 100 × 12, with each voxel of size, 0.3*mm ×* 0.3*mm ×* 0.9*mm*. 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.

Real data experiments:(a) and (b) are the source and targets *S*_{0} 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.

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]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |