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

**|**HHS Author Manuscripts**|**PMC2865691

Formats

Article sections

- Abstract
- 1. Introduction
- 2. The space of Diffusion Tensors
- 3. Tensor Spline Approximation
- 4. Experimental results
- 5. Conclusions
- References

Authors

Related links

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

Published in final edited form as:

Proc IEEE Comput Soc Conf Comput Vis Pattern Recognit. 2006 June 17; 2006: 86.

doi: 10.1109/CVPRW.2006.179PMCID: PMC2865691

NIHMSID: NIHMS144130

University of Florida, Gainesville, FL 32611, USA

Angelos Barmpoutis: ude.lfu.esic@uopmraba; Baba C. Vemuri: ude.lfu.esic@irumev; John R. Forder: ude.lfu.ibm@redrofj

See other articles in PMC that cite the published article.

In this paper, we present a novel and robust spline approximation algorithm given a noisy symmetric positive definite (SPD) tensor field. Such tensor fields commonly arise in the field of Medical Imaging in the form of Diffusion Tensor (DT) MRI data sets. We develop a statistically robust algorithm for constructing a tensor product of B-splines – for approximating and interpolating these data – using the Riemannian metric of the manifold of SPD tensors. Our method involves a two step procedure wherein the first step uses Riemannian distances in order to evaluate a tensor spline by computing a weighted intrinsic average of diffusion tensors and the second step involves minimization of the Riemannian distance between the evaluated spline curve and the given data. These two steps are alternated to achieve the desired tensor spline approximation to the given tensor field. We present comparisons of our algorithm with four existing methods of tensor interpolation applied to DT-MRI data from fixed heart slices of a rabbit, and show significantly improved results in the presence of noise and outliers. We also present validation results for our algorithm using synthetically generated noisy tensor field data with outliers. This interpolation work has many applications e.g., in DT-MRI registration, in DT-MRI Atlas construction etc. This research was in part funded by the NIH ROI NS42075 and the Department of Radiology, University of Florida.

Data processing and analysis of matrix-valued image data is becoming quite common as imaging sensor technology advances allow for the collection of matrix-valued data sets. In Medical Imaging, in the last decade, it has become possible to collect magnetic resonance image (MRI) data that measures the apparent diffusivity of water in tissue in vivo. A 2-tensor has been commonly used to approximate the diffusivity profile at each lattice point of the image lattice [2]. This approximation yields a diffusion tensor (DT-MRI) data set which is a matrix-valued image. These tensors are elements of the space of the 3×3 positive-definite matrices denoted by *P*(3). Mathematically, these positive definite Diffusion Tensors belong to a Riemannian Symmetric space, where a Riemannian metric assigns an inner product to each point of this space. By using this metric, one can compute geodesic distances between the elements of this space (diffusion tensors) and compute various statistics on this space [4, 6, 9, 10].

Processing of DT-MRI data sets has scientific significance for health, information and other sciences. Most of the applications require preprocessing of the DT-MRI data which more often than not involves interpolation of the diffusion tensor fields. Directly performing interpolation of the coefficients of the Diffusion Tensor matrices, does not preserve most of the properties, such as the values of the determinant and eigenvalues etc., of the Diffusion Tensors. Thus, it motivates one to seek alternative methods to achieve interpolation.

In [4], a Riemannian metric was proposed for geodesic computation between two tensors. Having a tensor field, e.g. a volume image of DT-MRI, we can use the geodesic curves between spatially consecutive tensors in order to interpolate the dataset. However, this geodesic curve computation between tensors (in [4]) did not impose higher order smoothness constraints in achieving the interpolation. Thus, although there is continuity of the obtained interpolated dataset, higher order continuity is lacking.

Recently, a Log-Euclidean metric was proposed in [1] for computing with tensors. In their work, the elements from the space of positive definite tensors *P*(3), are mapped to their tangent space, *Sym*(3), using the matrix logarithmic mapping. The tangent space forms a vector space of dimension ^{6} for diffusion tensors in ^{3}. Therefore, one can use the Euclidean norm for computations in this tangent space and finally by using the inverse mapping, the data are mapped back to the space of positive definite matrices *P*(3). More recently, a similar approach for interpolation of tensor fields was proposed in [3]. In their work, the tensors are mapped to the tangent space of a reference tensor, using the Riemannian Logarithmic mapping. The elements of the tangent space are treated as vectors, which are interpolated using radial basis functions and finally are mapped back using the inverse mapping. Both frameworks are quite interesting and have advantages due to their high computational efficiency in comparison to the geodesic distance framework. However, a statistically robust version of these frameworks need to be developed and although the geodesics obtained by such procedures are “similar” to those obtained by using the Riemannian metric, quantitatively they are not the same. This is because by mapping the data to the tangent space, we lose the information about the curvature of the Riemannian manifold of diffusion tensors.

Interpolation of matrix-valued images can be done by filtering such datasets using various regularization methods. Regularization of matrix-valued images can be achieved using a PDE-based approach proposed in [12]. Although this approach can be used for denoising and interpolation of matrix-valued images, the lack of use of an appropriate metric that is defined on the space of positive definite tensors could lead to undesirable limitations in the solution. Another tensor field regularization method was proposed in [13] using Normalized Convolution and Markov Random Fields (MRF) in a Bayesian Framework. The SPD tensors are treated as vectors in 6D and their components are treated independently. This can lead to inaccuracies in the regularized solution.

In this paper, we present a novel and robust spline approximation algorithm given a noisy diffusion tensor field. We develop a robust algorithm for constructing a tensor product of B-splines using the Riemannian metric of the manifold of symmetric positive definite tensors. We present comparisons of our algorithm with four existing methods of interpolation for Diffusion Tensor MRI data from fixed heart slices of a rabbit, and show much improved results in the presence of noise and outliers. We also present validation results for our algorithm using synthetically generated noisy tensor field data with outliers.

The rest of the paper is organized into the following sections: In section 2, we present the geometry of the space of Diffusion Tensors and the related mathematical background. In section 3, we make a brief review of B-splines. This is followed by the presentation of a novel algorithm for computing splines on a given symmetric positive definite tensor field in section 3. Next, we present tensor splines using the Log-Euclidean metric as an improvement over recent work in [1]. Then, we present a robust tensor spline approximation algorithm. Finally, in section 4, comparisons of our algorithm with four existing methods of interpolation for DT MRI data are presented. Validation results using synthetically generated noisy tensor field data with outliers are also presented.

In this section we present the geometry of the space of Diffusion Tensors and the related mathematical background that is required for doing computations in this space. More detailed expositions on this material maybe found in [9, 6, 4]. The space of Diffusion Tensors can be formulated as a Riemannian Symmetric space [5], where a Riemannian metric assigns an inner product to each point of this space. By using this metric we can compute geodesic distances between diffusions tensors and calculate statistics on this space [4, 6, 9].

A symmetric positive definite 2-tensor is a *n*×*n* matrix that belongs to the space of all symmetric positive-definite matrices denoted by *P*(*n*). DT-MRI datasets are a field of such tensor usually in two or three-dimensions. *P*(*n*) space is a connected Riemannian manifold [5]. The Lie group *GL*^{+}(*n*) of all *n*×*n* real matrices **g** with |**g**| > 0 acts on *P*(*n*) and smoothly maps its elements back to *P*(*n*) via the group action *ϕ*(**g**, **p**) = **gpg ^{T}**.

At each point **p** of *P*(*n*) the Riemannian metric is given by the following inner product **X**, **Y**** _{p}** =trace(

In figure 1, **p**_{1} is a point in *P*(*n*) and **X** is a tangent vector at **p**_{1}. There is a unique geodesic *γ*(*t*) starting at **p**_{1} for *t* = 0 and having *γ*′(0) = **X**. The point on the geodesic *γ*(*t*) at *t* = 1 can be computed by using the Riemannian exponential map *γ*(1) = *Exp*_{p1} (**X**). The Riemannian exponential map *Exp*** _{p}**(

Tangent space of the manifold *M* of diffusion tensors at point **p**_{1}. The tangent vector **X** points to the direction of geodesic *γ*(*t*) between the points **p**_{1} and **p**_{2}.

$${\mathit{Exp}}_{\mathbf{\text{p}}}(\mathbf{\text{X}})=(\mathbf{\text{gv}})\mathit{\text{exp}}(\sum ){(\mathbf{\text{gv}})}^{\mathbf{\text{T}}}$$

(1)

where **g** *GL*^{+}(*n*) and **p=gg*** ^{T}*. Let

$$\mathit{\text{exp}}(\mathbf{\text{X}})={\sum}_{\mathbf{\text{i}}=\mathbf{0}}^{\infty}\frac{\mathbf{1}}{\mathbf{\text{i}}\mathbf{!}}{\mathbf{X}}^{\mathbf{\text{i}}}.$$

(2)

The inverse map is the Riemannian Logarithmic map *Log*_{p1}(**p _{2}**), which maps

$${\mathit{\text{Log}}}_{{\mathbf{\text{p}}}_{1}}({\mathbf{\text{p}}}_{2})=(\mathbf{\text{gu}})\mathit{\text{log}}(\mathbf{\Lambda}){(\mathbf{\text{gu}})}^{\mathbf{\text{T}}}$$

(3)

where, **g** *GL*^{+}(*n*) and **p**_{1} = **gg ^{T}**. Let

For any real number *t* the point *γ*(*t*) in the geodesic between **p**_{1} and **p**_{2} can be computed as

$$\gamma (t)=(\mathbf{\text{gu}})\mathit{\text{Diag}}({\mathbf{\Lambda}}^{t}){(\mathbf{\text{gu}})}^{T}$$

(4)

where **g**, **u** and **Λ** are the same as in the equation of Riemannian Logarithmic map, and *Diag*(**X**) are the diagonal elements of symmetric matrix **X**. Equation 4 can be easily derived from *γ*(*t*) = *Exp*** _{p}** (

Using the Riemannian metric one can also compute statistics of Diffusion Tensors [4, 6, 9]. The average tensor of a set of Diffusion Tensors, can now be computed as that tensor which minimizes the sum of squared Riemannian distances between itself and the given set of tensors.

In the following section the Riemannian exponential and logarithmic maps and the equation of a geodesic between two diffusion tensors will be used in order to define and compute tensor splines.

In this section we present a novel and robust spline approximation algorithm given a noisy symmetric positive definite tensor field. Our algorithm involves the use of the Riemannian distance between SPD tensors in order to evaluate a tensor spline by computing a weighted intrinsic average of tensors. By intrinsic average, we mean, an average that minimizes the Riemannian distance between the unknown average and the members of the population whose average is being sought. This module is then used in a robust tensor product B-spline fitting method involving the minimization of the Riemannian distance between the spline function and the data.

This section has three subsections. First we make a brief review of B-splines. Next, we present a novel algorithm for computing splines on a given symmetric positive definite tensor field. After that we present tensor splines using the Log-Euclidean metric (described earlier). Finally a robust tensor spline approximation technique is presented.

The equation for *k*-1 degree B-spline with *n*+1 control points (*c*_{0}, *c*_{1}, …, *c _{n}*) and

$$S(t)=\sum _{i=0}^{n}{N}_{i,k}(t){c}_{i}$$

(5)

where *t*_{0} ≤ *t* ≤ *t _{n}*

$${N}_{i,1}=\{\begin{array}{ll}1& \text{if}\phantom{\rule{0.2em}{0ex}}{t}_{i}\le t<{t}_{i+1}\\ 0& \text{otherwise}\end{array}$$

(6)

and

$${N}_{i,k}(t)={N}_{i,k-1}(t)\frac{t-{t}_{i}}{{t}_{i+k-1}-{t}_{i}}+{N}_{i+1,k-1}(t)\frac{{t}_{i+k}-t}{{t}_{i+k}-{t}_{i+1}}$$

(7)

*N _{i,k}*(

One useful property of the functions *N _{i,k}*(

In the symmetric positive definite tensor space *P*(*n*), a tensor *γ*(*t*), where *t* *R* is a scalar, that lies on the geodesic between two tensors **p**_{1} and **p**_{2} is given by equation (4). By definition, a geodesic curve is the shortest path between two points. Having a dataset of tensors, e.g., a volume of DT-MRI, we can use the geodesic curves between spatially consecutive tensors in order to interpolate the dataset. However, a geodesic curve does not contain information about the neighborhood of the two interpolating points. Thus, although there is continuity of the interpolated dataset, there is lack of higher degree of continuity. It is more natural to have a higher degree of continuity in the interpolant within smoothly varying regions. Recent work in [7], on continuous tensor field approximation achieves smoothness, however, a Riemannian framework is not employed for tensor calculations. In this section we define tensor splines which are curves in the input space *P*(*n*) × * ^{m}*, (

A 3rd degree tensor spline *S*(*t*), that passes through 5 given tensors **p**_{i} of a 1-D tensor field. The given points **p**_{i} and the points of the tensor spline *S*(*t*) are SPD matrices, elements of the Riemannian manifold *M*. However the tensor spline is in *P*(*n*) × **...**

Let us assume that we have a set of *N* diffusion tensors (**p**_{0}, **p**_{1}, …, **p**_{N}_{−1}) on a one dimensional grid, and we want to interpolate between them. Linear (1st degree) interpolation on the tensor space can be achieved by computing points on the geodesics connecting two consecutive diffusion tensors. Higher degree interpolation can be done by using a set of control points and a knots vector. A *k*-1 degree tensor spline that fits to our data requires *N*-1+*k*-1 control points (**c**_{0}, **c**_{1}, …, **c**_{N}_{−1+}_{k}_{−2}) which are also tensors and *N*+2(*k*-1) monotonically increasing knots (*t*_{−}_{k}_{+1}, *t*_{−}_{k}_{+2}, …, *t _{N}*

$$S(t)={\text{\u2211}^{\sim}}_{i=0}^{n}{w}_{i}{\mathbf{\text{c}}}_{i}$$

(8)

The intrinsic weighted average of tensors is defined using Riemannian distance instead of Euclidean, and it is the minimizer of the function

$$\underset{\mu \in P(n)}{min}\rho (\mu )=\underset{\mu \in P(n)}{min}\frac{1}{2}\sum _{\mathbf{\text{i}}=0}^{n}{w}_{i}\mathit{\text{dist}}{(\mu ,{\mathbf{\text{c}}}_{i})}^{2}$$

(9)

where *dist*(., .) is the Riemannian geodesic distance. The weighted average can be computed using a gradient descent algorithm which is an extension of the algorithm described by Pennec [8] for computing the mean of tensors. The gradient of *ρ*(*μ*) is given by

$${\nabla}_{\mu}\rho =-{\sum}_{i=0}^{n}{w}_{i}{\mathit{\text{Log}}}_{\mu}({\mathbf{\text{c}}}_{i})$$

(10)

Thus the intrinsic weighted average of a set of diffusion tensors can be computed by the following procedure:

input : c_{1}, …, c_{N}P(n) |

w_{1}, …, w weights_{N} |

output: μ P(n), the weighted mean |

μ_{0} ← I ; |

i ← 0 ; |

while ‖ X ‖ > _{i}e do |

X ← −_{i};_{μi} ρ |

μ_{i+1} ← Exp (_{μi}X) ;_{i} |

end |

In order to fit a tensor spline to the diffusion tensor data, we have to approximate the control tensors of such a spline. A tensor spline that fits to our data, minimizes the Riemannian distance of the given tensors from the tensor spline curve.

$$E=\frac{1}{2N}\sum _{i=0}^{N-1}\mathit{\text{dist}}{(S({t}_{i}),{\mathbf{\text{p}}}_{i})}^{2}$$

(11)

In this energy expression, the Riemannian metric should be used for the distance calculation, since the tensor space, where the data and control points belong, is a Riemannian manifold. We need to find a set of control points (**c**_{0}, **c**_{1}, …, **c**_{N}_{−1+}_{k}_{−2}) that form the spline *S*(*t*) which minimizes the energy *E*. The gradient of E with respect to **c*** _{j}* is then given by,

$${\nabla}_{{\mathbf{\text{c}}}_{j}}E=\frac{1}{2N}{\sum}_{i=0}^{N-1}{\nabla}_{S({t}_{i})}\mathit{\text{dist}}{(S({t}_{i}),{\mathbf{\text{p}}}_{i})}^{2}{\nabla}_{{\mathbf{\text{c}}}_{j}}S({t}_{i})$$

(12)

The gradient of the square distance between *S*(*t _{i}*) and

$${\nabla}_{S({t}_{i})}\mathit{\text{dist}}{(S({t}_{i}),{\mathbf{\text{p}}}_{i})}^{2}=-2{\mathit{\text{Log}}}_{S({t}_{i})}({\mathbf{\text{p}}}_{i})$$

(13)

where *Log _{S(ti)}* (

$${\nabla}_{S({t}_{i})}\mathit{\text{dist}}{(S({t}_{i}),{\mathbf{\text{p}}}_{i})}^{2}\approx -2{\mathbf{\Lambda}}_{{\mathbf{\text{c}}}_{j}}({\mathbf{\text{p}}}_{i},S({t}_{i}))$$

(14)

Furthermore, the gradient of *S*(*t _{i}*) with respect to

$${\nabla}_{{\mathbf{\text{c}}}_{j}}S({t}_{i})={N}_{j,k}\phantom{\rule{0.1em}{0ex}}({t}_{i})$$

(15)

Using equations 15 and 14 in 12 we obtain

$${\nabla}_{{\mathbf{\text{c}}}_{j}}E=-\frac{1}{N}{\sum}_{i=0}^{N-1}{\mathbf{\Lambda}}_{{\mathbf{\text{c}}}_{j}}({\mathbf{\text{p}}}_{i},S({t}_{i})){N}_{j,k}\phantom{\rule{0.1em}{0ex}}({t}_{i})$$

(16)

Starting with an initial guess of the control tensors, we can update them by using the gradient descent technique. The new values **c**′* _{j}* of control tensors will be

$$\mathbf{\text{c}}{\prime}_{j}={\mathit{\text{Exp}}}_{{\mathbf{\text{c}}}_{j}}\phantom{\rule{0.2em}{0ex}}\left(\frac{1}{N}{\sum}_{i=0}^{N-1}{\mathbf{\Lambda}}_{{\mathbf{\text{c}}}_{j}}({\mathbf{\text{p}}}_{i},S({t}_{i})){N}_{j,k}\phantom{\rule{0.1em}{0ex}}({t}_{i})\right)$$

(17)

where *Exp*. (.) is the Riemannian exponential map (eq. 1). The initial guess of the control tensors can be either the given data or the average tensor of the given tensors. The gradient descent algorithm is summarized below:

The error introduced by the approximation of equation 14 can be large, if the tensor spline approximation *S*(*t _{i}*) is far from the target

input : N tensors (p_{0}, …, p_{N−1}), |

N+2(k-1) monotonically increasing |

knots (t_{−k+1}, …, t_{N−1+k−1}) |

k, and a small value e |

output: N-1+k-1 control tensors |

(c_{0}, …, c_{N−1+k−2}) |

‖ X_{0} ‖ ← e+1; |

while
$\text{\u2211}_{j}\Vert {\mathbf{\text{X}}}_{j}\Vert >e$ do |

for j=0 to N-1+k-2 do |

X_{j} ← zero matrix ; |

for i=0 to N-1 do |

X_{j} ← X_{j} + |

Λ_{cj}(p_{i},S(t))_{i}N(_{j,k}t) ;_{i} |

end |

c′ ← _{j}Exp_{cj} X_{j} ; |

end |

c ← c′ ; |

end |

The knot sequence can be parameterized by different ways that preserve the monotonically increasing property of the knots series. In the experiments we used the series of integer numbers 1,2,3… as the knot sequence.

Tensor Splines can be fit to higher dimensional tensor fields, by simply extending the presented algorithms to the new dimensionality. For example consider the case of a 2D *N* × *M* tensor field. In this case we are doing interpolation on *P*(*n*) × ^{2}. A *k* − 1 degree tensor spline that fits on our data requires (*N* − 1 + *k* − l)×(*M* − 1 + *k* − 1) control tensors and (*N* + 2(*k* − l))×(*M* + 2(*k* − 1)) monotonically increasing in both the dimensions, the knots (*t*_{−}_{k}_{+1,−}_{k}_{+1}, …, *t _{N}*

Recently, in [1], Arsigny et. al., proposed a new Log-Euclidean metric for tensor calculations. By using this metric, the diffusion tensors are mapped, with the matrix logarithmic map to the space of the symmetric matrices *Sym*(*n*). Therefore, we can make Euclidean calculations in this space and finally by using the matrix exponential mapping, the data are mapped back to the space of positive definite matrices *P*(*n*). Although the tensors obtained by this procedure are “similar” to those obtained by using the Riemannian metric, quantitatively they are not the same. This is because by mapping the data to the tangent space, we lose the information about the curvature of the Riemannian manifold of diffusion tensors.

Using the Log-Euclidean metric we can also fit a spline to the logarithmically mapped data into the symmetric space, and after that we can map the interpolated tensors back to the tensor space using the exponential mapping. In the section 4, we provide a quantitative comparison between tensor splines and Log-Euclidean Splines.

The presence of outliers is common in DT-MRI data due to noise in the original data obtained from the DT-MRI sensors [11]. A robust interpolation algorithm should detect these outlier tensors and reject them from further consideration in any processing algorithms applied to the dataset.

A robust function can be used in the energy function, in order to weight the given data **p*** _{i}* appropriately. We need to use a robust function that assigns weights in the interval [0, 1], where weights which are almost zero imply rejection of the corresponding data point. Furthermore, high weights should be assigned to the data points whose distance from the corresponding points on the unknown spline curve is small and on the other hand lower weights should be assigned to the data points whose distance from the corresponding points on the unknown spline is larger. Let us consider the following function
$\phi (x)=-{\mathit{\text{se}}}^{-\frac{x}{s}}$. The derivative of this function is
$\psi (x)={e}^{-\frac{x}{s}}$ and it has the aforementioned properties. By using the above function

$$E=\frac{1}{2N}{\sum}_{i=0}^{N-1}\phi (\mathit{\text{dist}}{(S({t}_{i}),{p}_{i})}^{2})$$

(18)

The gradient of the energy with respect to the control tensors now becomes

$$\begin{array}{l}{\nabla}_{{\mathbf{\text{c}}}_{j}}E=-\frac{1}{N}{\sum}_{i=0}^{N-1}\psi (\mathit{\text{dist}}{(S({t}_{i}),{\mathbf{\text{p}}}_{i})}^{2})\times {\mathbf{\Lambda}}_{{\mathbf{\text{c}}}_{j}}({\mathbf{\text{p}}}_{i},S({t}_{i})){N}_{j,k}\phantom{\rule{0.1em}{0ex}}({t}_{i})\end{array}$$

In the above equation the quantity *ψ*(*dist*(*S*(*t _{i}*),

In this section, we present several experiments with noisy synthetic as well as real diffusion tensor data. We also present comparisons with four other existing methods to demonstrate the performance of our proposed algorithm.

We synthesized a tensor field on a 2D lattice of size 33×33 (see fig. 4(a)). For the generation of this field, we chose four diffusion tensors which were then positioned at the corners of the given rectangular lattice. These diffusion tensors were arbitrarily selected from a real DT-MRI data set of a fixed rabbit heart. The *tensors at the rest of the lattice points were obtained via geodesic interpolation* using the methods described earlier. The left plate of figure 4(a) depicts the Fractional Anisotropy (FA) map (see [2] for definition of FA) and the right plate depicts the dominant eigenvector field. The *x, y, z* components of the dominant eigenvectors of the tensors were assigned Red-Green-Blue (RGB) colors respectively for visualization purposes.

We then subsampled the tensor field by a factor of four and figure 3(a) is a depiction of the same. The Gaussian ellipsoids that represents the diffusion tensors of the synthetic tensor field are shown in the right plate of this figure, and the FA map is shown in the left plate. Then, we generated two noisy tensor fields by adding noise to the subsampled tensor field. The first noisy tensor field (fig. 3(b)) contains 10% outliers which were also arbitrarily selected from a set of outlier diffusion tensors that were contaminating the previously mentioned data set.

For the generation of the second noisy tensor field (fig. 3(c)) a 3×3 symmetric matrix **X*** _{i,j}* was randomly generated from a Gaussian distribution, for every lattice point of the subsampled tensor field

Finally we interpolated the noisy tensor fields by a factor of 4, using the four existing methods a) Riemannian geodesic interpolation [4], b) Log-Euclidean linear interpolation [1], c) Tangent space interpolation [3], d) PDE-based anisotropic non-linear diffusion/interpolation [12], e) our robust tensor spline method with degree 3 (cubic) and also using our own modification of the Log-Euclidean metric, namely, f) a Log-Euclidean spline. The results are presented in figure 4. We use two methods to measure the distance of the estimated tensor fields from the ground truth tensor field of figure 4(a); a) the Riemannian metric and b) the length of the error vector defined as the difference between the ground truth and the estimated dominant eigenvector of the diffusion tensors respectively. These errors are computed at each voxel and the mean and standard deviation of these errors are reported in table 1 for the two noisy data sets. As evident and expected, the error is much smaller for our algorithm in comparison to the others. Moreover, since our method is robust to the presence of outliers in the data, its performance is least affected among all the methods as the corrupting mechanism of the original tensor field was changed. This conclusively demonstrates superior interpolation performance of our algorithm over other existing methods.

In figure 5, the first row depicts three images. (a) Shows a slice of a T2-weighted MR image obtained using a zero gradient field. A white box has been used to identify the region of interest (ROI) containing a reasonably noisy section in the DT-MRI slice. The rest of the plates in this figure depict the processing limited to this ROI. Figure (b) depicts the FA map of the data. Higher brightness indicates higher FA values in this depiction. Note the random changes in brightness values in this image caused due to the presence of noise and outliers in the tensor field. Figure (c) depicts the FA computed after applying tensor spline approximation. Note the relatively smooth appearance in brightness corresponding to the FA values, which is biologically more consistent in the shown ROI. The part (d) of the same figure depicts an ellipsoidal visualization (in color) of the original tensor field in the ROI. Note the presence of the outliers in the bottom right part of this image. Image (e) of the same figure shows the result after applying our tensor spline approximation. Notice that all of the outliers have been rejected and the field has been interpolated smoothly.

a)A DT-MRI slice of a rabbit heart. The next figures come from the marked region of this image. b)The Fractional Anisotropy (FA) map of the original data. c)The FA map after the fitting of a cubic tensor spline in the data. d) and e) show the ellipsoid **...**

Finally figure 6 shows another real data example. In this figure the proposed Tensor Spline approximation algorithm is compared with the recently proposed Log-Euclidean geodesic interpolation [1], discussed earlier. The vector field of the dominant eigenvector of the original data is depicted in 6(b). Note that the data contain noise and outliers. Figures 6(c) and 6(d) show the interpolation results of Log-Euclidean geodesic and Tensor Spline algorithm respectively. Since the Log-Euclidean method is not robust the interpolation result still contains a large amount of the outliers seen in the original data. Notice that in the Tensor Spline approximation results these outliers have been rejected. Figure 6(a) depicts the FA map of the original and the interpolated data. In the case of Tensor Spline approximation (6(a) bottom) the FA map of the approximated data, has also been interpolated smoothly.

In this paper, we presented a novel and robust spline approximation algorithm that we dub, tensor-splines, given a noisy symmetric positive definite tensor field. We presented comparisons of our algorithm with four existing methods of interpolation for Diffusion Tensor MRI data from fixed heart slices of a rabbit, and presented significantly improved results in the presence of noise and outliers. We also presented validation results for our algorithm using synthetically generated noisy tensor field data with outliers. Our future work will be focused on applying this algorithm as a module in applications such as registration of DT-MRI data sets, tractography etc. [15, 14].

1. Arsigny V, Fillard P, Pennec X, Ayache N. Fast and Simple Calculus on Tensors in the Log-Euclidean Framework. Proceedings of MICCAI, LNCS. 2005:259–267. [PubMed]

2. Basser P, 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]

3. Fillard P, Arsigny V, Pennec X, Thompson P, Ayache N. Extrapolation of sparse tensor fields: Application to the modeling of brain variability. IPMI. 2005:27–38. [PubMed]

4. Fletcher P, Joshi S. Principal geodesic analysis on symmetric spaces: Statistics of diffusion tensors. Proc of CVAMIA. 2004:87–98.

5. Helgason S. Differential geometry, Lie groups, and symmetric spaces. American Mathematical Society; 2001.

6. Lenglet C, Rousson M, Deriche R, Faugeras O. Statistics on Multivariate Normal Distributions: A Geometric Approach and its Application to Diffusion Tensor MRI. INRIA Research Report 5242. 2004

7. Pajevic S, Aldroubi A, Basser P. Visualization and Image Processing of Tensor Fields. Springer; 2005. Continuous Tensor Field Approximation of Diffusion Tensor MRI Data.

8. Pennec X. Probabilities and statistics on Riemannian manifolds: basic tools for geometric measurements. IEEE Workshop on Nonlinear Signal and Image Processing. 1999

9. Pennec X, Fillard P, Ayache N. A Riemannian framework for tensor computing. International Journal of Computer Vision. 2005:65.

10. San-Jose R, Estepar S, Haker S, Westin CF. Riemannian Mean Curvature Flow. ISVC, volume 3804 of LNCS. 2005 December;:613–620.

11. Wang Z, Vemuri B, Chen Y, Mareci T. A Constrained Variational Principle for Direct Estimation and Smoothing of the Diffusion Tensor Field from DWI. IPMI. 2003:660–671. [PubMed]

12. Weickert J, Welk M. Tensor field interpolation with PDEs. Preprint 142, Saarland University. 2005

13. Westin CF, Martin-Fernandez M, Alberola-Lopez C, Ruiz-Alzola J, Knutsson H. Tensor Field Regularization Using Normalized Convolution and markov random fields in a bayesian framework. In: Weickert J, Hagen H, editors. Visualization and Image Processing of Tensor Fields. Springer; 2005. In Press.

14. Xu D, Mori S, Solaiyappan M, Zijl PC, Davatzikos C. A framework for callosal fiber distribution analysis. Neuroimage. 2002 November;17(3):1131–1143. [PubMed]

15. Zhang H, Yushkevich PA, Gee JC. Registration of Diffusion Tensor Images. CVPR. 2004;1:842–847.