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

**|**HHS Author Manuscripts**|**PMC2920614

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Formulation of the Problem
- 3. First order variation
- 4. Numerical implementation
- 5. Experimental Results and Discussions
- 6. Conclusion
- References

Authors

Related links

Proc IEEE Comput Soc Conf Comput Vis Pattern Recognit. Author manuscript; available in PMC 2010 August 12.

Published in final edited form as:

Proc IEEE Comput Soc Conf Comput Vis Pattern Recognit. 2006 July 5; 2006: 67.

doi: 10.1109/CVPRW.2006.65PMCID: PMC2920614

NIHMSID: NIHMS189851

Yan Cao, Center for Imaging Science, Johns Hopkins University;

Yan Cao: ude.uhj.sic@nay; Michael I. Miller: ude.uhj.sic@mim; Susumu Mori: ude.uhj.irm@umusus; Raimond L. Winslow: ude.uhj.emb@wolsniwr; Laurent Younes: ude.uhj.sic@senuoy

See other articles in PMC that cite the published article.

This paper proposes a method to match diffusion tensor magnetic resonance images (DT-MRI) through the large deformation diffeomorphic metric mapping of tensor fields on the image volume, resulting in optimizing for geodesics on the space of diffeomorphisms connecting two diffusion tensor images. A coarse to fine multi-resolution and multi-kernel-width scheme is detailed, to reduce both ambiguities and computation load. This is illustrated by numerical experiments on DT-MRI brain and images.

Diffusion Tensor Magnetic Resonance Imaging (DT-MRI) is a non-invasive technique that measures the diffusion of water molecules in biological tissues [16, 14, 20, 9]. Diffusion refers to the random translational movement of molecules (also called Brownian motion). The rate and direction of the movement is closely related to the structural anisotropy of the medium, hence it can be used to infer the tissue micro structure. Diffusion is described mathematically by a symmetric second order tensor. Orientation of the principal eigenvector of the diffusion tensor is known to align with fiber tracts [16, 20]. DT-MRI is rapidly gaining importance in neuroscience and other fields but tools for exploring and analyzing DT-MRI data are still not widely available.

Image registration is inevitable whenever images acquired from different subjects, at different times, or from different scanners need to be combined or compared for analysis or visualization. Alexander and Gee [1] adapted the multi-resolution elastic matching algorithm using tensor similarity measures for matching DT-MRIs. Tensor reorientation is not included in the energy function, but tensors are reoriented in each iteration according to the estimated displacement field. Ruiz-Alzola [19] proposed a unified framework for registration of medical volumetric multi-valued data by matching local areas with high structure. Local structure is detected by measuring “cornerness”. Then the Kriging estimator was used to interpolate the displacement field in the remaining areas. Tensor reorientation is performed after the registration. Guimond et al. [11] proposed a deformable registration method using tensor characteristics that are invariant to rigid transformations to avoid the complexity of tensor reorientation. The invariants used in their experiments are the tensor eigenvalues. Rohde et al. [18] proposed an intensity based registration method capable of performing affine and nonlinear registration of multi-channel images. The multi-channel similarity measure considered both the similarity between same-index image channels and the similarity between image channels with different index. But tensor reorientation is performed after affine registration and after the nonlinear registration is complete. Zhang et al [21] proposed a piecewise local affine registration algorithm. Leemans et al [12] proposed an affine multi-channel matching method for DT-MRIs based on mutual information. The Diffusion Weighted Images was used as the channels.

This paper proposes a method to match diffusion tensor magnetic resonance images through an extension of the Large Deformation Diffeomorphic Metric Mapping (LD-DMM) framework [5] to this type of dataset, resulting in optimizing for geodesics on the space of diffeomorphisms connecting two tensor fields. In [7], LDDMM framework was extended to matching vector fields. In the paper, we extends it further to the full tensor fields. Establishing this methodology continues our efforts on establishing the construction and computation of metric spaces for studying anatomical structures. Our method includes tensor reorientation in the optimization, hence the measure of difference is calculated correctly during the optimization. Convergence of the algorithm is guaranteed.

In our approach, anatomy is modeled as a deformable template. The images are defined on an open bounded set and are viewed as the orbit under the diffeomorphic transformations **G** acting on a template. In detail, let the background space Ω be a bounded domain in **R*** ^{d}* and

Elements in **G** are defined through the solutions of the nonlinear Eulerian transport equation

$${\stackrel{.}{\varphi}}_{t}={v}_{t}({\varphi}_{t}),\phantom{\rule{0.38889em}{0ex}}{\varphi}_{0}(x)=x,\phantom{\rule{0.38889em}{0ex}}t\in [0,1].$$

(1)

For each *t*, *v _{t}* is a vector field on Ω which belongs to some Hilbert space

For any *M* **I** and *ϕ* **G**, *M* can be decomposed as
$M={\lambda}_{1}{e}_{1}{e}_{1}^{T}+{\lambda}_{2}{e}_{2}{e}_{2}^{T}+{\lambda}_{3}{e}_{3}{e}_{3}^{T}$, where *λ*_{1}, *λ*_{2} and *λ*_{3} are 3 eigenvalues, *λ*_{1} ≥ *λ*_{2} ≥ *λ*_{3}, and *e*_{1}, *e*_{2} and *e*_{3} are the corresponding eigenvectors of *M*. The action is defined as follows [2]:

$$\varphi \xb7M\doteq \left({\lambda}_{1}{\widehat{e}}_{1}{\widehat{e}}_{1}^{T}+{\lambda}_{2}{\widehat{e}}_{2}{\widehat{e}}_{2}^{T}+{\lambda}_{3}{\widehat{e}}_{3}{\widehat{e}}_{3}^{T}\right)\circ {\varphi}^{-1},$$

(2)

where

$${\widehat{e}}_{1}=\frac{D\varphi {e}_{1}}{\left|\right|D\varphi {e}_{1}\left|\right|},$$

(3)

$${\widehat{e}}_{2}=\frac{D\varphi {e}_{2}-\langle {\widehat{e}}_{1},D\varphi {e}_{2}\rangle {\widehat{e}}_{1}}{\sqrt{{\left|\right|D\varphi {e}_{2}\left|\right|}^{2}-{\langle {\widehat{e}}_{1},D\varphi {e}_{2}\rangle}^{2}}},$$

(4)

$${\widehat{e}}_{3}={\widehat{e}}_{1}\times {\widehat{e}}_{2}.$$

(5)

This is the Gram-Schmidt ortho-normalization.

Here in Equation (2) to (5) everything depends on *x* and *Dϕ* denotes the Jacobian matrix of the function *ϕ*. In summary, the eigenvector deformation strategy assumes that the eigenvalues of the diffusion tensor do not change, the principal eigenvector changes according to the Jacobian, but preserve the length, i.e. *e*_{1} changes to *Dϕ e*_{1}*/*||*Dϕ e*_{1}||. The plane generated by the first and the second principal vectors *e*_{1} and *e*_{2} changes to the plane generated by *Dϕ e*_{1} and *Dϕ e*_{2}. Although the exact mechanism for diffusion anisotropy is not well understood, it is clear that this anisotropy directly reflects the presence of spatially oriented structures in the tissue [6]. Assume the tissue microstructure doesn’t change during the deformation, then the shape of the diffusion tensor should not change, hence we assume the eigenvalues of the diffusion tensor do not change. The eigenvectors reflect the orientation of the tissue, so they will be affected by the deformation.

It can be shown that the action is well-defined, i.e., the choice of the eigenvectors of *M* doesn’t affect the value of *ϕ.M*. It can also be shown that the action defined above is a group action, i.e., *ϕ*.(*ψ.M*) = (*ϕ*○*ψ*).*M*, using the chain rule of differentiation and the fact that a Jacobian acts linearly.

Let the notation
${\varphi}_{st}^{v}:\mathrm{\Omega}\to \mathrm{\Omega}$ denote the composition
${\varphi}_{st}^{v}={\varphi}_{t}\circ {\varphi}_{s}^{-1}$. Then
${\varphi}_{st}^{v}(x)$ is the solution at time *t* of the equation
${\scriptstyle \frac{dy}{dt}}={v}_{t}(y)$ with initial condition *y _{s}* =

Given two elements template *M ^{temp}* and target

$$\text{dist}(\varphi \xb7{M}^{\mathit{temp}},{M}^{\mathit{targ}})={\left|\right|\varphi \xb7{M}^{\mathit{temp}}-{M}^{\mathit{targ}}\left|\right|}_{F}$$

(6)

$$=\text{trace}\left[{(\varphi \xb7{M}^{\mathit{temp}}-{M}^{\mathit{targ}})}^{\ast}(\varphi \xb7{M}^{\mathit{temp}}-{M}^{\mathit{targ}})\right],$$

(7)

where the notation *A*^{*} denotes the conjugate transpose of *A*. The mapping *ϕ* is generated as the end-point *ϕ* = *ϕ*_{1} = *ϕ*_{01} of the flow of smooth time-dependent vector fields *v _{t}*

Estimation of the optimal diffeomorphic transformation connecting two images of tensor fields *M ^{temp}* and

$$\widehat{v}={\text{argmin}}_{v}E(v),$$

(8)

$$E(v)\doteq {\int}_{0}^{1}{\left|\right|{v}_{t}\left|\right|}_{V}^{2}dt+c{\int}_{\mathrm{\Omega}}\text{dist}{(\varphi \xb7{M}^{\mathit{temp}},{M}^{\mathit{targ}})}^{2}dx,$$

(9)

where ||·||* _{V}* is an appropriate functional norm on the velocity field

The Frobenius or Hilbert-Schmidt norm we have used to compare the tensor can be related to a noise model in which matrix coefficients are subject to additive independent Gaussian noise. From a metric point of view several authors [8, 15, 4] have recently introduced a specific affine invariant distance designed for the comparison of symmetric tensors. The squared distance between two tensors *M*_{1} and *M*_{2} is given by the sum of square of the logarithm of the eigenvalues of
${M}_{1}^{-1/2}{M}_{2}{M}_{1}^{-1/2}$. Even if a direct use of this distance is not feasible here, a first order approximation of it, namely

$$\text{dist}{({M}_{1},{M}_{2})}^{2}=\text{trace}[{({M}_{1}^{-1/2}{M}_{2}{M}_{1}^{-1/2}-I)}^{2}],$$

(10)

(where *I* is the identity matrix) can be used with slight modification of the algorithm (will be introduced later) in our framework.

To solve the variational problem (8), we need to compute the first order variation of the energy function with respect to the vector field *v*.

Let *v* and *h* be time-dependent vector fields, *t* [0, 1], such that for each *t*, *v _{t}*,

$${\partial}_{h}f(v)\doteq {\frac{d}{d\epsilon}f(v+\epsilon h)|}_{\epsilon =0}\doteq {\int}_{0}^{1}<{({\nabla}_{v}f)}_{t},{h}_{t}{>}_{V}dt.$$

(11)

Here * _{h}f^{v}* is called the Gâteaux differential of

We construct the Hilbert space *V* using the theory of reproducing kernel Hilbert spaces (RKHS) [3]. The Hilbert space *V* is defined as follows [7]:

$$V\doteq \left\{v\in {L}^{2}(\mathrm{\Omega},{\mathbf{R}}^{d}):\exists u\in {L}^{2}(\mathrm{\Omega},{\mathbf{R}}^{d}),v(\xb7)={\int}_{\mathrm{\Omega}}\stackrel{\sim}{K}(\xb7,y)u(y)dy\right\},$$

(12)

where is continuous, positive definite and has compact support in Ω. The inner product on *V* is defined as <*v*, *v*′>* _{V}* =<

Let

$${M}^{\mathit{temp}}={M}_{0}={\lambda}_{1}{e}_{01}{e}_{01}^{T}+{\lambda}_{2}{e}_{02}{e}_{02}^{T}+{\lambda}_{3}{e}_{03}{e}_{03}^{T},$$

(13)

$${M}_{t}={\varphi}_{0t}^{v}\xb7{M}_{0}={\lambda}_{1}\circ {\varphi}_{t0}^{v}{w}_{t1}{w}_{t1}^{T}+{\lambda}_{2}\circ {\varphi}_{t0}^{v}{w}_{t2}{w}_{t2}^{T}+{\lambda}_{3}\circ {\varphi}_{t0}^{v}{w}_{t3}{w}_{t3}^{T},$$

(14)

where *w _{t}*

$${({\nabla}_{v}E)}_{t}=2{v}_{t}+2c{\int}_{\mathrm{\Omega}}\left|\right|D{\varphi}_{t1}\left|\right|\sum _{\tau \zeta}\left({({M}_{t})}^{\tau \zeta}-{({M}_{1}^{\mathit{targ}})}^{\tau \zeta}\circ {\varphi}_{t1}\right)\{{K}^{(\beta )}{({D}^{2}{\varphi}_{t1}^{v})}_{\beta \gamma}^{\alpha}{({(D{\varphi}_{t1})}^{-1})}_{\eta}^{\gamma}[({\lambda}_{1}-{\lambda}_{2})\circ {\varphi}_{t0}{({w}_{t2})}^{\alpha}{({w}_{t1})}^{\eta}({({w}_{t2})}^{\tau}{({w}_{t1})}^{\zeta}+{({w}_{t1})}^{\tau}{({w}_{t2})}^{\zeta})+({\lambda}_{1}-{\lambda}_{3})\circ {\varphi}_{t0}{({w}_{t3})}^{\alpha}{({w}_{t1})}^{\eta}({({w}_{t3})}^{\tau}{({w}_{t1})}^{\zeta}+{({w}_{t1})}^{\tau}{({w}_{t3})}^{\zeta})+({\lambda}_{2}-{\lambda}_{3})\circ {\varphi}_{t0}{({w}_{t3})}^{\alpha}{({w}_{t2})}^{\eta}({({w}_{t3})}^{\tau}{({w}_{t2})}^{\zeta}+{({w}_{t2})}^{\tau}{({w}_{t3})}^{\zeta})]-{K}^{(\beta )}{({DM}_{t})}_{\beta}^{\tau \zeta}+{({DK}^{\beta})}_{\gamma}{(D{\varphi}_{t1}^{v})}_{\beta}^{\alpha}{({(D{\varphi}_{t1})}^{-1})}_{\eta}^{\gamma}[({\lambda}_{1}-{\lambda}_{2})\circ {\varphi}_{t0}{({w}_{t2})}^{\alpha}{({w}_{t1})}^{\eta}({({w}_{t2})}^{\tau}{({w}_{t1})}^{\zeta}+{({w}_{t1})}^{\tau}{({w}_{t2})}^{\zeta})+({\lambda}_{1}-{\lambda}_{3})\circ {\varphi}_{t0}{({w}_{t3})}^{\alpha}{({w}_{t1})}^{\eta}({({w}_{t3})}^{\tau}{({w}_{t1})}^{\zeta}+{({w}_{t1})}^{\tau}{({w}_{t3})}^{\zeta})+({\lambda}_{2}-{\lambda}_{3})\circ {\varphi}_{t0}{({w}_{t3})}^{\alpha}{({w}_{t2})}^{\eta}({({w}_{t3})}^{\tau}{({w}_{t2})}^{\zeta}+{({w}_{t2})}^{\tau}{({w}_{t3})}^{\zeta})]\}dx.$$

(15)

Here we denote *K _{x}*(

In our definition of the Hilbert space *V*, the element *v* is defined by means of *L*^{2} functions *u*. In our experiments, we find it is more efficient to work on *u* directly. We define

$${\partial}_{h}f(u)\doteq {\frac{d}{d\epsilon}f(u+\epsilon h)|}_{\epsilon =0}\doteq {\int}_{0}^{1}<{({\nabla}_{u}f)}_{t},{h}_{t}{>}_{{L}^{2}}dt.$$

(16)

So, for *h ^{α}*(

$${h}^{\alpha}(x)={\langle {\stackrel{\sim}{K}}_{x}^{(\alpha )}(y),h(y)\rangle}_{{L}^{2}},$$

(17)

$${({D}_{x}h)}_{\beta}^{\alpha}={\langle {({D}_{x}{\stackrel{\sim}{K}}_{x}^{(\alpha )})}_{\beta}(y),h(y)\rangle}_{{L}^{2}}.$$

(18)

$${({\nabla}_{u}E)}_{t}=2{u}_{t}+2c{\int}_{\mathrm{\Omega}}\left|\right|D{\varphi}_{t1}\left|\right|\sum _{\tau \zeta}\left({({M}_{t})}^{\tau \zeta}-{({M}_{1}^{\mathit{targ}})}^{\tau \zeta}\circ {\varphi}_{t1}\right)\{{\stackrel{\sim}{K}}^{(\beta )}{({D}^{2}{\varphi}_{t1}^{u})}_{\beta \gamma}^{\alpha}{({(D{\varphi}_{t1})}^{-1})}_{\eta}^{\gamma}[({\lambda}_{1}-{\lambda}_{2})\circ {\varphi}_{t0}{({w}_{t2})}^{\alpha}{({w}_{t1})}^{\eta}({({w}_{t2})}^{\tau}{({w}_{t1})}^{\zeta}+{({w}_{t1})}^{\tau}{({w}_{t2})}^{\zeta})+({\lambda}_{1}-{\lambda}_{3})\circ {\varphi}_{t0}{({w}_{t3})}^{\alpha}{({w}_{t1})}^{\eta}({({w}_{t3})}^{\tau}{({w}_{t1})}^{\zeta}+{({w}_{t1})}^{\tau}{({w}_{t3})}^{\zeta})+({\lambda}_{2}-{\lambda}_{3})\circ {\varphi}_{t0}{({w}_{t3})}^{\alpha}{({w}_{t2})}^{\eta}({({w}_{t3})}^{\tau}{({w}_{t2})}^{\zeta}+{({w}_{t2})}^{\tau}{({w}_{t3})}^{\zeta})]-{\stackrel{\sim}{K}}^{(\beta )}{({DM}_{t})}_{\beta}^{\tau \zeta}+{(D{\stackrel{\sim}{K}}^{\beta})}_{\gamma}{(D{\varphi}_{t1}^{u})}_{\beta}^{\alpha}{({(D{\varphi}_{t1})}^{-1})}_{\eta}^{\gamma}[({\lambda}_{1}-{\lambda}_{2})\circ {\varphi}_{t0}{({w}_{t2})}^{\alpha}{({w}_{t1})}^{\eta}({({w}_{t2})}^{\tau}{({w}_{t1})}^{\zeta}+{({w}_{t1})}^{\tau}{({w}_{t2})}^{\zeta})+({\lambda}_{1}-{\lambda}_{3})\circ {\varphi}_{t0}{({w}_{t3})}^{\alpha}{({w}_{t1})}^{\eta}({({w}_{t3})}^{\tau}{({w}_{t1})}^{\zeta}+{({w}_{t1})}^{\tau}{({w}_{t3})}^{\zeta})+({\lambda}_{2}-{\lambda}_{3})\circ {\varphi}_{t0}{({w}_{t3})}^{\alpha}{({w}_{t2})}^{\eta}({({w}_{t3})}^{\tau}{({w}_{t2})}^{\zeta}+{({w}_{t2})}^{\tau}{({w}_{t3})}^{\zeta})]\}dx.$$

(19)

We extend the method used in [7] for the numerical implementation of the estimation of the optimal vector field *v*. Let Ω = [0, 1]^{2}/[0, 1]^{3} be the 2D/3D background space. Ω is discretized into a uniform square/cube mesh. The time space is discretized as *t ^{n}* =

The variational optimization of the energy functional is performed in a steepest descent scheme

$${u}^{k+1}={u}^{k}-\epsilon {({\nabla}_{{u}^{k}}E)}_{t},$$

(20)

where (* _{u}E*)

A hierarchical multi-resolution matching strategy is used to reduce ambiguity issues and computation load. It is employed from coarse to fine, and results achieved on one resolution are considered as approximations for the next finer level. We generate the image pyramid by reducing the resolution from one level to the next by a factor of 2 using linear interpolation.

The width of the kernel significantly affects the deformation. When the width of the kernel is large, the deformation is smoother but the constraints on the deformation will also be stronger, hence the template may not have enough freedom to deform to the target, but it is also not easily trapped into the local minimum. When the width of the kernel is small, the template deforms more freely, hence the local matching will be better, but it may also be easily trapped into the local minimum. Hence, at the early stage, we should use bigger width to get a rough matching, then at later stage, use smaller width to refine the matching.

If we change the width of the kernel, the vector space *V* also changes. We need to find the relationship between two vector spaces with different kernel width. This is a difficult problem for general kernels, but with Gaussian kernels, we have a simple solution due to their special property. Recall that the vector space *V* is defined as follows:

$$V\doteq \left\{v\in {L}^{2}(\mathrm{\Omega},{\mathbf{R}}^{d}):\exists u\in {L}^{2}(\mathrm{\Omega},{\mathbf{R}}^{d}),v(\xb7)={\int}_{\mathrm{\Omega}}\stackrel{\sim}{K}(\xb7,y)u(y)dy\right\}.$$

(21)

where (*x*, *y*) = *η*(*x*)*K*_{1}(*x*, *y*) *η*(*y*)|_{Ω}, *K*_{1} is a Gaussian function and *η*(*x*) is a function which smoothes *K*_{1} near the boundary and has zero boundary conditions.

Suppose *V*_{1} and *V*_{2} are two vector spaces with kernel size *σ*_{1} and *σ*_{2}, *σ*_{1} > *σ*_{2} > 0. Then
${K}_{1}(x,y)={\scriptstyle \frac{1}{{Z}_{1}}}exp(-{\scriptstyle \frac{{\left|\right|x-y\left|\right|}^{2}}{2{\sigma}_{1}^{2}}})$, and
${K}_{2}(x,y)={\scriptstyle \frac{1}{{Z}_{2}}}exp(-{\scriptstyle \frac{{\left|\right|x-y\left|\right|}^{2}}{2{\sigma}_{2}^{2}}})$, where *Z*_{1} and *Z*_{2} are constants to make the integration of the kernel be 1.

$$v(x)={\int}_{\mathrm{\Omega}}\stackrel{\sim}{K}(x,y){u}_{1}(y)dy$$

(22)

$$={\int}_{\mathrm{\Omega}}\eta (x)\frac{1}{{Z}_{1}}exp\left[-\frac{{\left|\right|x-y\left|\right|}^{2}}{2{\sigma}_{1}^{2}}\right]\eta (y){u}_{1}(y)dy$$

(23)

$$={\int}_{\mathrm{\Omega}}\eta (x)\frac{1}{{Z}_{2}}exp\left[-\frac{{\left|\right|x-y\left|\right|}^{2}}{2{\sigma}_{2}^{2}}\right]{\int}_{\mathrm{\Omega}}\frac{1}{{Z}_{d}}exp\left[-\frac{{\left|\right|y-{y}_{1}\left|\right|}^{2}}{2({\sigma}_{1}^{2}-{\sigma}_{2}^{2})}\right]\eta ({y}_{1}){u}_{1}({y}_{1})d{y}_{1}dy$$

(24)

$$\doteq {\int}_{\mathrm{\Omega}}\eta (x)\frac{1}{{Z}_{2}}exp\left[-\frac{{\left|\right|x-y\left|\right|}^{2}}{2{\sigma}_{2}^{2}}\right]{u}_{2}(y)dy.$$

(25)

with

$${Z}_{d}=\frac{{Z}_{1}}{{Z}_{2}}{\left(\frac{{\sigma}_{2}}{{\sigma}_{1}}\sqrt{2\pi ({\sigma}_{1}^{2}-{\sigma}_{2}^{2})}\right)}^{d}$$

(26)

$${u}_{2}(x)={\int}_{\mathrm{\Omega}}\frac{1}{{Z}_{d}}exp\left[-\frac{{\left|\right|x-y\left|\right|}^{2}}{2({\sigma}_{1}^{2}-{\sigma}_{2}^{2})}\right]\eta (y){u}_{1}(y)dy.$$

(27)

Equation (27) provides the transition from *u*_{1} to *u*_{2} between the spaces *V*_{1} and *V*_{2}.

At a given resolution and given kernel width:

- 1) Initialize the algorithm with
*k*= 0, ${u}_{{t}_{j}}^{k}=0$,*j*= 0, 1, …,*N*or any initial guess. - For each iteration
*k*, - 4) Calculate the new velocity
*u*^{k}^{+1}=*u*−^{k}*ε*_{uk}*E*. Find the optimal step size*ε*using golden section line search algorithm [17] in the direction of steepest descent. If*ε*<*TOL*, stop the iteration, otherwise, return to step 2).

The multi-resolution matching algorithm is as follows: (1.) Rigid matching. (2.) Generate multi-resolution image pyramid by reducing the resolution from one level to the next by a factor of 2. (3.) Pick an initial kernel width, from the coarsest level to the finest level, do the following: i). Use the fixed resolution matching algorithm to find the optimal vector field *u* for that level. ii). Use linear interpolation to convert the optimal vector field *u* from one level to the next finer level. (4.) Reduce the kernel width, convert the optimal vector field *u* to the new space using Equation (27). (5.) Use the fixed resolution matching algorithm to find the optimal vector field *u* for the new kernel at the finest resolution. (6.) Calculate the corresponding *v* and *ϕ*_{01} and *ϕ*_{10}.

The algorithm is implemented in C++ for 3D tensor field images. In our experiments, the time interval of the flow is discretized into 20 steps, with step size Δ*t* = 0.05. We use *c* = 5. We first use 0.075 and then 0.025 as the standard deviation of the Gaussian kernel with the grid normalized to [0, 1]^{3}.

We have done experiments on various normal human brain diffusion tensor images. Figure 1 shows the tensor matching result for one pair of human diffusion tensor brain images in detail. We can see from the figures that the alignment of the diffusion tensor at each voxel improved a lot after the matching.

3D tensor matching of two normal human brains. Left column shows the tensor distribution of slice 30 before matching; Right column shows the tensor distribution of slice 30 after matching. Template and target are superimposed with blue color the template **...**

Figure 2 shows the comparison of several matching methods on one pair of normal human brain diffusion tensor images. The shape matching method refers to matching the geometrical shape of the two brains. The vector matching method refers to matching the principal eigenvector of the diffusion tensor at each voxel. All the methods are implemented under the LDDMM framework. It is clear that tensor matching is better than vector matching and vector matching is better than shape matching. This result is as expected since the tensor matching method uses the most information from the data and the shape matching method the least information among the three methods. Figure 3 shows the comparison of the deformations result from LD-DMM vector matching and tensor matching schemes. For the vector matching scheme, the deformation in areas with low fractional anisotropy (FA) values is very large since the principal directions of the tensor in those areas are not clear and kind of random. Forcing matching them causes unnecessary deformations. For the tensor matching scheme, the deformation in areas with low FA values is not forced. From Figure 2 and Figure 3, we can see that tensormatching scheme improves the matching quality in areas with low FA values. In the areas with high FA values, the two schemes produce similar results.

Comparison of the deformed template and the target for different LDDMM matching schemes. Top left panel shows the histogram of the tensor difference (Frobenius norm) at each voxel. Top right panel shows the histogram of the FA difference at each voxel. **...**

Comparison of the deformations result from LDDMM vector matching and tensor matching schemes. Top panel shows the FA weighted color-coded orientation map of the target. Second row shows the deformations from vector matching, third row shows the deformation **...**

In our tensor matching method, we assume the eigenvalues of the diffusion tensor doesn’t change. This is based on the assumption that the tissue micro-structure doesn’t change during the deformation and different subjects have similar tissue types. Left column of Figure 4 shows the histogram of the 3 eigenvalues of different normal brain diffusion tensor images. From the figure we can see that the distributions for different subjects are very similar. They are peaked at the similar values and have similar ranges. This shows our assumption is reasonable. However, due to changes of conditions, sometimes the eigenvalues did change for different scans. Right column of Figure 4 shows the histogram of the 3 eigenvalues of different normal canine heart diffusion tensor images. The distribution are quite different. This figure shows sometimes it is necessary to normalize the eigenvalues of the tensors before matching. We will investigate that problem in the future.

In conclusion, we have presented in this paper the large deformation diffeomorphic metric mapping of tensor fields. The optimal mapping is the endpoint of a geodesic path on the manifold of diffeomorphisms connecting two tensor fields. Finding the optimal mapping and the geodesic path is formulated as a variational problem over a vector field. A gradient descent based multi-resolution multi-kernel-width algorithm is implemented. We did the experiments on various 3D normal brain diffusion tensor magnetic resonance images. This method gives good results when the assumption that the eigenvalues of the diffusion tensor do not change during transformation holds. We expect this method be a useful tool for analysis diffusion tensor magnetic resonance images and other images which have similar properties.

This work was supported in part by NIH grant P41-RR15241, NIH grant P50 MH71616, NIH 1U24RR0211382-01, NIH grant 1-R01-HL70894 and the Falk Medical Trust.

To solve the variational problem (8), we need to compute the first order variation of the energy function with respect to the vector field *v*.

Let

$${M}^{\mathit{temp}}={M}_{0}={\lambda}_{1}{e}_{01}{e}_{01}^{T}+{\lambda}_{2}{e}_{02}{e}_{02}^{T}+{\lambda}_{3}{e}_{03}{e}_{03}^{T},$$

(28)

$${M}_{t}={\varphi}_{0t}^{v}\xb7{M}_{0}={\lambda}_{1}\circ {\varphi}_{t0}^{v}{w}_{t1}{w}_{t1}^{T}+{\lambda}_{2}\circ {\varphi}_{t0}^{v}{w}_{t2}{w}_{t2}^{T}+{\lambda}_{3}\circ {\varphi}_{t0}^{v}{w}_{t3}{w}_{t3}^{T},$$

(29)

where *w _{t}*

$${M}_{t}^{\epsilon}={\varphi}_{0t}^{v+\epsilon h}\xb7{M}_{0}=({\varphi}_{0t}^{v+\epsilon h}\circ {\varphi}_{t0}^{v})\xb7({\varphi}_{0t}^{v}\xb7{M}_{0})$$

(30)

$$=({\varphi}_{0t}^{v+\epsilon h}\circ {\varphi}_{t0}^{v})\xb7{M}_{t}\doteq (\text{id}+\epsilon f)\xb7{M}_{t}+o(\epsilon ).$$

(31)

Equation (31) defines the function *f*. Let *g* = id + *εf*, then *g*^{−1} = id − *εf* + *o*(*ε*) and *Dg* = id + *εDf*.

$${M}_{t}^{\epsilon}=({\lambda}_{1}\circ {\varphi}_{t0}^{v}{w}_{t1}^{\epsilon}{({w}_{t1}^{\epsilon})}^{T}+{\lambda}_{2}\circ {\varphi}_{t0}^{v}{w}_{t2}^{\epsilon}{({w}_{t2}^{\epsilon})}^{T}+{\lambda}_{3}\circ {\varphi}_{t0}^{v}{w}_{t3}^{\epsilon}{({w}_{t3}^{\epsilon})}^{T})\circ {g}^{-1}+o(\epsilon ),$$

(32)

where
${w}_{t1}^{\epsilon},{w}_{t2}^{\epsilon}$ and
${w}_{t3}^{\epsilon}$ are transformed eigenvectors of *w _{t}*

$$Dg{w}_{ti}={w}_{ti}+\epsilon Df{w}_{ti}+o(\epsilon ),$$

(33)

$$\frac{1}{\left|\right|Dg{w}_{ti}\left|\right|}=1-\epsilon <{w}_{ti},Df{w}_{ti}>+o(\epsilon )\phantom{\rule{0.38889em}{0ex}}i=1,2,3.$$

(34)

then

$${w}_{t1}^{\epsilon}={w}_{t1}+\epsilon <{w}_{t2},Df{w}_{t1}>{w}_{t2}+\epsilon <{w}_{t3},Df{w}_{t1}>{w}_{t3}+o(\epsilon ),$$

(35)

$${w}_{t2}^{\epsilon}={w}_{t2}-\epsilon <{w}_{t2},Df{w}_{t1}>{w}_{t1}+\epsilon <{w}_{t3},Df{w}_{t2}>{w}_{t3}+o(\epsilon ),$$

(36)

$${w}_{t3}^{\epsilon}={w}_{t3}-\epsilon <{w}_{t3},Df{w}_{t1}>{w}_{t1}-\epsilon <{w}_{t3},Df{w}_{t2}>{w}_{t2}+o(\epsilon ).$$

(37)

Let v and h in χ(Ω), then for x Ω,

$${({\partial}_{h}{\varphi}_{st}^{v}(x))}^{\alpha}={\int}_{s}^{t}{({D}_{{\varphi}_{su}^{v}(x)}{\varphi}_{ut}^{v})}_{\beta}^{\alpha}{({h}_{u})}^{\beta}\circ {\varphi}_{su}^{v}(x)du.$$

This lemma is proved in [5], but stated in a different way. Using this lemma, we can prove

$$f={\int}_{0}^{t}{({D}_{{\varphi}_{tu}^{v}(x)}{\varphi}_{ut}^{v})}_{\beta}^{\alpha}{({h}_{u})}^{\beta}\circ {\varphi}_{tu}^{v}(x)du$$

(38)

$$Df={\int}_{0}^{t}{({D}^{2}{\varphi}_{ut}^{v})}_{\beta \gamma}^{\alpha}\circ {\varphi}_{tu}^{v}{(D{\varphi}_{tu})}_{\eta}^{\gamma}{({h}_{u})}^{\beta}\circ {\varphi}_{tu}^{v}(x)+{(D{\varphi}_{ut}^{v})}_{\beta}^{\alpha}\circ {\varphi}_{tu}^{v}{(D{h}_{u})}_{\gamma}^{\beta}\circ {\varphi}_{tu}^{v}{(D{\varphi}_{tu})}_{\eta}^{\gamma}du$$

(39)

Let

$$\begin{array}{c}{E}_{1}(v)={\int}_{0}^{1}{\left|\right|{v}_{t}\left|\right|}_{V}^{2}dt,\\ {E}_{2}(v)={\int}_{\mathrm{\Omega}}{\left|\right|{M}_{1}-{M}_{1}^{\mathit{targ}}\left|\right|}^{2}dx\end{array}$$

(40)

We have

$${\partial}_{h}{E}_{1}=2{\int}_{0}^{1}<{v}_{t},{h}_{t}{>}_{V}dt,$$

(41)

$${\partial}_{h}{E}_{2}=2{\int}_{\mathrm{\Omega}}\sum _{\mathit{all}\phantom{\rule{0.16667em}{0ex}}\mathit{entries}}({M}_{1}-{M}_{1}^{\mathit{targ}}){\partial}_{h}{M}_{1}dx.$$

(42)

$${\partial}_{h}{M}_{t}={\frac{d}{d\epsilon}{M}_{t}^{\epsilon}|}_{\epsilon =0}$$

(43)

$$=({\lambda}_{1}-{\lambda}_{2})\circ {\varphi}_{t0}^{v}<{w}_{t2},Df{w}_{t1}>({w}_{t2}{w}_{t1}^{T}+{w}_{t1}{w}_{t2}^{T})+({\lambda}_{1}-{\lambda}_{3})\circ {\varphi}_{t0}^{v}<{w}_{t3},Df{w}_{t1}>({w}_{t3}{w}_{t1}^{T}+{w}_{t1}{w}_{t3}^{T})+({\lambda}_{2}-{\lambda}_{3})\circ {\varphi}_{t0}^{v}<{w}_{t3},Df{w}_{t2}>({w}_{t3}{w}_{t2}^{T}+{w}_{t2}{w}_{t3}^{T})-{DM}_{t}f.$$

(44)

Suppose *K* is the self-reproducing kernel of *V*, then *K* is a *d* × *d* matrix. We denote *K _{x}*(

$${h}^{\alpha}(x)={\langle {K}_{x}^{(\alpha )}(y),h(y)\rangle}_{V},$$

(45)

$${({D}_{x}h)}_{\beta}^{\alpha}={\langle {({D}_{x}{K}_{x}^{(\alpha )})}_{\beta}(y),h(y)\rangle}_{V}.$$

(46)

then

$${\partial}_{h}{E}_{2}=2{\int}_{0}^{1}\langle {\int}_{\mathrm{\Omega}}\left|\right|D{\varphi}_{u1}\left|\right|\sum _{\tau \zeta}\left({({M}_{u})}^{\tau \zeta}-{({M}_{1}^{\mathit{targ}})}^{\tau \zeta}\circ {\varphi}_{u1}\right)\{{K}^{(\beta )}{({D}^{2}{\varphi}_{u1}^{v})}_{\beta \gamma}^{\alpha}{({(D{\varphi}_{u1})}^{-1})}_{\eta}^{\gamma}[({\lambda}_{1}-{\lambda}_{2})\circ {\varphi}_{u0}{({w}_{u2})}^{\alpha}{({w}_{u1})}^{\eta}({({w}_{u2})}^{\tau}{({w}_{u1})}^{\zeta}+{({w}_{u1})}^{\tau}{({w}_{u2})}^{\zeta})+({\lambda}_{1}-{\lambda}_{3})\circ {\varphi}_{u0}{({w}_{u3})}^{\alpha}{({w}_{u1})}^{\eta}({({w}_{u3})}^{\tau}{({w}_{u1})}^{\zeta}+{({w}_{u1})}^{\tau}{({w}_{u3})}^{\zeta})+({\lambda}_{2}-{\lambda}_{3})\circ {\varphi}_{u0}{({w}_{u3})}^{\alpha}{({w}_{u2})}^{\eta}({({w}_{u3})}^{\tau}{({w}_{u2})}^{\zeta}+{({w}_{u2})}^{\tau}{({w}_{u3})}^{\zeta})]-{K}^{(\beta )}{({DM}_{u})}_{\beta}^{\tau \zeta}+{({DK}^{\beta})}_{\gamma}{(D{\varphi}_{u1}^{v})}_{\beta}^{\alpha}{({(D{\varphi}_{u1})}^{-1})}_{\eta}^{\gamma}[({\lambda}_{1}-{\lambda}_{2})\circ {\varphi}_{u0}{({w}_{u2})}^{\alpha}{({w}_{u1})}^{\eta}({({w}_{u2})}^{\tau}{({w}_{u1})}^{\zeta}+{({w}_{u1})}^{\tau}{({w}_{u2})}^{\zeta})+({\lambda}_{1}-{\lambda}_{3})\circ {\varphi}_{u0}{({w}_{u3})}^{\alpha}{({w}_{u1})}^{\eta}({({w}_{u3})}^{\tau}{({w}_{u1})}^{\zeta}+{({w}_{u1})}^{\tau}{({w}_{u3})}^{\zeta}){+({\lambda}_{2}-{\lambda}_{3})\circ {\varphi}_{u0}{({w}_{u3})}^{\alpha}{({w}_{u2})}^{\eta}({({w}_{u3})}^{\tau}{({w}_{u2})}^{\zeta}+{({w}_{u2})}^{\tau}{({w}_{u3})}^{\zeta})]\}dx,{h}_{u}\rangle}_{V}du.$$

(47)

So

$${({\nabla}_{v}E)}_{t}=2{v}_{t}+2c{\int}_{\mathrm{\Omega}}\left|\right|D{\varphi}_{t1}\left|\right|\sum _{\tau \zeta}\left({({M}_{t})}^{\tau \zeta}-{({M}_{1}^{\mathit{targ}})}^{\tau \zeta}\circ {\varphi}_{t1}\right)\{{K}^{(\beta )}{({D}^{2}{\varphi}_{t1}^{v})}_{\beta \gamma}^{\alpha}{({(D{\varphi}_{t1})}^{-1})}_{\eta}^{\gamma}[({\lambda}_{1}-{\lambda}_{2})\circ {\varphi}_{t0}{({w}_{t2})}^{\alpha}{({w}_{t1})}^{\eta}({({w}_{t2})}^{\tau}{({w}_{t1})}^{\zeta}+{({w}_{t1})}^{\tau}{({w}_{t2})}^{\zeta})+({\lambda}_{1}-{\lambda}_{3})\circ {\varphi}_{t0}{({w}_{t3})}^{\alpha}{({w}_{t1})}^{\eta}({({w}_{t3})}^{\tau}{({w}_{t1})}^{\zeta}+{({w}_{t1})}^{\tau}{({w}_{t3})}^{\zeta})+({\lambda}_{2}-{\lambda}_{3})\circ {\varphi}_{t0}{({w}_{t3})}^{\alpha}{({w}_{t2})}^{\eta}({({w}_{t3})}^{\tau}{({w}_{t2})}^{\zeta}+{({w}_{t2})}^{\tau}{({w}_{t3})}^{\zeta})]-{K}^{(\beta )}{({DM}_{t})}_{\beta}^{\tau \zeta}+{({DK}^{\beta})}_{\gamma}{(D{\varphi}_{t1}^{v})}_{\beta}^{\alpha}{({(D{\varphi}_{t1})}^{-1})}_{\eta}^{\gamma}[({\lambda}_{1}-{\lambda}_{2})\circ {\varphi}_{t0}{({w}_{t2})}^{\alpha}{({w}_{t1})}^{\eta}({({w}_{t2})}^{\tau}{({w}_{t1})}^{\zeta}+{({w}_{t1})}^{\tau}{({w}_{t2})}^{\zeta})+({\lambda}_{1}-{\lambda}_{3})\circ {\varphi}_{t0}{({w}_{t3})}^{\alpha}{({w}_{t1})}^{\eta}({({w}_{t3})}^{\tau}{({w}_{t1})}^{\zeta}+{({w}_{t1})}^{\tau}{({w}_{t3})}^{\zeta})+({\lambda}_{2}-{\lambda}_{3})\circ {\varphi}_{t0}{({w}_{t3})}^{\alpha}{({w}_{t2})}^{\eta}({({w}_{t3})}^{\tau}{({w}_{t2})}^{\zeta}+{({w}_{t2})}^{\tau}{({w}_{t3})}^{\zeta})]\}dx.$$

(48)

Yan Cao, Center for Imaging Science, Johns Hopkins University.

Michael I. Miller, Center for Imaging Science, Johns Hopkins University.

Susumu Mori, Department of Radiology and Kennedy Krieger Institute, Johns Hopkins University, School of Medicine.

Raimond L. Winslow, Center for Cardiovascular Bioinformatics & Modeling and the Whitaker Biomedical Engineering Institute, Johns Hopkins University.

Laurent Younes, Center for Imaging Science, Johns Hopkins University.

1. Alexander DC, Gee JC. Elastic matching of diffusion tensor images. Computer Vision and Image Understanding. 2000;77:233–250.

2. Alexander DC, Pierpaoli C, Basser P, Gee JC. Spatial transformations of diffusion tensor magnetic resonance images. IEEE TMI. 2001;20(11):1131–1139. [PubMed]

3. Aronszajn N. Theory of reproducing kernels. Trans Amer Math Soc. 1950 May;68(3):337–404.

4. Batchelor PG, Moakher M, Atkinson D, Calamante F, Connelly A. A rigorous framework for diffusion tensor calculus. Magnetic Resonance in Medicine. 2005;53:221–225. [PubMed]

5. Beg MF, Miller MI, Trouvé A, Younes L. Computing metrics via geodesics on flows of diffeomorphisms. IJCV. 2005;61(2):139–157.

6. Bihan DL, van Zijl P. From the diffusion coefficient to the diffusion tensor. NMR in Biomedicine. 2002;15:431–434. [PubMed]

7. Cao Y, Miller MI, Winslow RL, Younes L. Large deformation diffeomorphic metric mapping of vector fields. IEEE TMI. 2005;24(9):1216–1230. [PubMed]

8. Corouge I, Fletcher PT, Joshi S, Gilmore JH, Gerig G. Fiber tract-oriented statistics for quantitative diffusion tensor mri analysis. MICCAI. 2005 [PubMed]

9. Dou J, Tseng WYI, Reese TG, Wedeen VJ. Combined diffusion and strain mri reveals structure and function of human myocardial laminar sheets in vivo. Magnetic Resonance in Medicine. 2003;50:107–113. [PubMed]

10. Durran DR. Numerical Methods for Wave Equations in Geophysical Fluid Dynamics. Springer; 1999.

11. Guimond A, Guttmann CRG, Warfield SK, Westin C-F. Deformable registration of dt-mri data based on transformation invariant tensor characteristics. Proceedings of 2002 IEEE International Symposium on Biomedical Imaging; July 2002.

12. Leemans A, Sijbers J, Backer SD, Vandervliet E, Parizel PM. Affine coregistration of diffusion tensor magnetic resonance images using mutual information. Advanced Concepts for Intelligent Vision Systems: 7th International Conference, ACIVS 2005, LNCS 3708; September 2005.pp. 523–530.

13. Luenberger DG. Optimization by Vector Space Methods. John Wiley & Sons, Inc; 1969.

14. Mori S, Barker PB. Diffusion magnetic resonance imaging: Its principle and applications. The Anatomical Record (New Anat) 1999;257:102–109. [PubMed]

15. Pennec X, Fillard P, Ayache N. A riemannian framework for tensor computing. IJCV. 2005;66:41–66.

16. Pierpaoli C, Jezzard P, Basser PJ, Barnett A, Chiro GD. Diffusion tensor mr imaging of the human brain. Radiology. 1996;201(3):637–648. [PubMed]

17. Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical Recipes in C, The Art of Scientific Computing. 2. Cambridge University Press; 1992.

18. Rohde GK, Pajevic S, Pierpaoli C. Multi-channel registration of diffusion tensor images using directional information. Proceedings of 2004 IEEE International Symposium on Biomedical Imaging: From Nano to Macro; April 2004.

19. Ruiz-Alzola J, Westin CF, Warfield SK, Alberola C, Maier S, Kikinis R. Nonrigid registration of 3d tensor medical data. Medical Image Analysis. 2002;6:143–161. [PubMed]

20. Scollan DF, Holmes A, Winslow RL, Forder J. Histological validation of myocardial microstructure obtained from diffusion tensor magnetic resonance imaging. American Journal of Physiology (Heart and Circulatory Physiology) 1998;275:2308–2318. [PubMed]

21. Zhang H, Yushkevich PA, Gee JC. Towards diffusion profile image registration. Proceedings of the IEEE International Symposium on Biomedical Imaging (ISBI); 2004.

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. |