|Home | About | Journals | Submit | Contact Us | Français|
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  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  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.  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.  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  proposed a piecewise local affine registration algorithm. Leemans et al  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  to this type of dataset, resulting in optimizing for geodesics on the space of diffeomorphisms connecting two tensor fields. In , 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 Rd and G a group of diffeomorphisms on Ω. Let the images be functions M: Ω → Rd×d that associate to each point x Ω the diffusion tensor, a symmetric second order tensor. G acts on the set I of all images.
Elements in G are defined through the solutions of the nonlinear Eulerian transport equation
For each t, vt is a vector field on Ω which belongs to some Hilbert space V. The meaning of this equation is that we have a path on the space of transformations ϕt, where t is the time. Under the action of this flow, a particle which starts at a point x(0), traces out a path x(t) = ϕt(x(0)), where (t) = vt(x(t)), vt is the speed at time t.
For any M I and ϕ G, M can be decomposed as , where λ1, λ2 and λ3 are 3 eigenvalues, λ1 ≥ λ2 ≥ λ3, and e1, e2 and e3 are the corresponding eigenvectors of M. The action is defined as follows :
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. e1 changes to Dϕ e1/||Dϕ e1||. The plane generated by the first and the second principal vectors e1 and e2 changes to the plane generated by Dϕ e1 and Dϕ e2. 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 . 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 denote the composition . Then is the solution at time t of the equation with initial condition ys = x. The function is called the flow associated to v starting at s. It is defined on [0, 1] × Ω. In physical terms, is the motion of a particle which starts at x at time s. For s, r, t [0, 1], we have . This equation states that the movement from s to t is the same as the movements from s to r, and then from r to t.
Given two elements template Mtemp and target Mtarg, we would like to find an optimal matching ϕ in the space of diffeomorphisms between Mtemp and Mtarg which minimizes a distance function dist(ϕ.Mtemp, Mtarg) on the tensor fields. We use the Frobenius or Hilbert-Schmidt matrix norm,
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 vt V via Equation (1).
Estimation of the optimal diffeomorphic transformation connecting two images of tensor fields Mtemp and Mtarg is done via the following variational problem:
where ||·||V is an appropriate functional norm on the velocity field vt(·).
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 M1 and M2 is given by the sum of square of the logarithm of the eigenvalues of . Even if a direct use of this distance is not feasible here, a first order approximation of it, namely
(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, vt, ht V and . For a function f(v), we define
Here hfv is called the Gâteaux differential of f at v with increment h and v f is called the Fréchet derivative of f .
where is continuous, positive definite and has compact support in Ω. The inner product on V is defined as <v, v′>V =<u, u′>L2. V is a RKHS with reproducing kernel K(x, y) = ∫Ω (x, z)(z, y)dz. To take advantage of the Fourier transform in the calculations related to the kernel, we define in our experiments, (x, y) = η(x)K1(x, y)η(y)|Ω, where K1 is a Gaussian function and η(x) is a function which smoothes K1 near the boundary and has zero boundary conditions.
where wt1, wt2 and wt3 are transformed eigenvectors of e01, e02 and e03 under the action of . Then
Here we denote Kx(y) = K(y, x) and K(α) the α-th column of K. See appendix for the details on calculations.
In our definition of the Hilbert space V, the element v is defined by means of L2 functions u. In our experiments, we find it is more efficient to work on u directly. We define
So, for hα(x) L2(Ω, Rd), we have
We extend the method used in  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 tn = nΔt, n = 0, 1, …, N and . Let denote the value of f at xj, in time step tn for kth iteration.
The variational optimization of the energy functional is performed in a steepest descent scheme
where (uE)t is computed using Equation (19). The derivatives are calculated via central differences with enforcement that Dxf (x) = 0 if the function f is a function of I0 and f(x) = 0. The integrations are calculated by Fourier transformation since they can be viewed as convolutions. We use the Gaussian kernel functions as K1 in our experiments. Other kernel functions can be used. The step size ε is decided by the golden section line search algorithm  in the direction of steepest descent. We use a second order semi-Lagrangian scheme  to integrate the velocity field v to generate the map ϕt0 and ϕt1.
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:
where (x, y) = η(x)K1(x, y) η(y)|Ω, K1 is a Gaussian function and η(x) is a function which smoothes K1 near the boundary and has zero boundary conditions.
Suppose V1 and V2 are two vector spaces with kernel size σ1 and σ2, σ1 > σ2 > 0. Then , and , where Z1 and Z2 are constants to make the integration of the kernel be 1.
Equation (27) provides the transition from u1 to u2 between the spaces V1 and V2.
At a given resolution and given kernel width:
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.
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.
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.
where wt1, wt2 and wt3 are transformed eigenvectors of e01, e02 and e03 under the action of . Let v and h be time-dependent vector fields, t [0, 1], such that for each t, vt, ht V and .
Equation (31) defines the function f. Let g = id + εf, then g−1 = id − εf + o(ε) and Dg = id + εDf.
where and are transformed eigenvectors of wt1, wt2 and wt3 under the action of g. It is easy to show that
Let v and h in χ(Ω), then for x Ω,
This lemma is proved in , but stated in a different way. Using this lemma, we can prove
Suppose K is the self-reproducing kernel of V, then K is a d × d matrix. We denote Kx(y) = K(y, x). Let K(α) be the α-th column of K. We have
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.