|Home | About | Journals | Submit | Contact Us | Français|
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 . 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 , 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 ) 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  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 . 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 . 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  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 . 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 , 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 . 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) = gpgT.
At each point p of P(n) the Riemannian metric is given by the following inner product X, Yp =trace(g−1Xp−1Yg−T), where g GL+(n) and p = ggT and X,Y Sym(n) are two tangent vectors at the point p P(n). At any point p P(n) the tangent space is defined by the space of the n×n symmetric matrices.
In figure 1, p1 is a point in P(n) and X is a tangent vector at p1. There is a unique geodesic γ(t) starting at p1 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) = Expp1 (X). The Riemannian exponential map Expp(X) at a point p P(n) maps a tangent vector X at p to γ(1), where γ(t) is the geodesic starting at γ(0) = p having γ′(0) = X. Expp(X) can be computed using
where g GL+(n) and p=ggT. Let Y=g−1Xg−T, then v and Σ are the matrix of eigenvectors and the diagonal matrix of eigenvalues of Y respectively. Therefore in equation (1) exp(Σ) is the matrix exponential map which is given by,
The inverse map is the Riemannian Logarithmic map Logp1(p2), which maps p2 P(n) to a tangent vector X at p1, such that if γ(0) = p1 and γ(1) = p2 then X=γ′(0). Logp1 (p2) can be computed using,
where, g GL+(n) and p1 = ggT. Let y = g−1p2g−T, then u and Λ are the matrix of eigenvectors and the diagonal matrix of eigenvalues of y respectively. The distance between the two tensors p1 and p2 equals to ‖ Logp1 (p2)‖p1 where ‖.‖p1 is the Riemannian norm at the point p1.
For any real number t the point γ(t) in the geodesic between p1 and p2 can be computed as
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) = Expp (X, t) = (gv)exp(Σt)(gv)T  and 3.
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 (c0, c1, …, cn) and n+k+1 knots (t−k+1, t−k+2, …, tn+1), is
where t0 ≤ t ≤ tn+1−(k−1). Each control point is associated with a basis function Ni,k where
Ni,k(t) functions are polynomials of degree k-1. Cubic basis functions Ni,4 can be used for a 3rd degree B-spline. Knots must be series of monotonic increasing numbers.
One useful property of the functions Ni,k(t) is that Ni,k(t) ≥ 0, for all i and . Considering the above properties, functions Ni,k(t) behave as blending functions and eq. 5 is a weighted average of the control points ci.
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 p1 and p2 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 , 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, (P(3) × 2 in the case of a 2D diffusion tensor field) constructed using the geometry of the space of symmetric positive definite tensors. Note that doing spline interpolation on P(n) is not meaningful in itself, in this context, as there is no ordering of points (matrices) on P(n), which reflects the ordering imposed by the lattice in 2D/3D. Hence, interpolation on P(3) × 2 or P(3) × 3 is the right thing to do because the tensors in DTI are defined on a lattice in 2D or 3D. Another issue to keep in mind is that we are defining tensor-splines by doing weighted intrinsic averages on P(n) and choosing the weight functions to be B-splines. As an illustration of interpolation on a 1D grid of tensors, figure (2) depicts the idea of using weight functions (B-splines here) to perform weighted average of tensors using the Riemannian metric. This weighted averaging leads to the desired degree spline interpolant of the diffusion tensor data. Applications of tensor splines are numerous e.g., higher order interpolation between diffusion tensors, smoothing and filtering diffusion tensor data sets, and compression of such data sets.
Let us assume that we have a set of N diffusion tensors (p0, p1, …, pN−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 (c0, c1, …, cN−1+k−2) which are also tensors and N+2(k-1) monotonically increasing knots (t−k+1, t−k+2, …, tN−1+k−1). A tensor S(t), where t [tj, tj+1), which is a point on a tensor spline, can be now computed by generalizing the equation (5) to the space of tensors. We can compute the value S(t) of the k-1 degree B-spline of tensors for a particular t value, by calculating a weighted intrinsic average, , of the control tensors ci, where the weights equal to the basis functions wi = Ni,k(t), discussed earlier.
The intrinsic weighted average of tensors is defined using Riemannian distance instead of Euclidean, and it is the minimizer of the function
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  for computing the mean of tensors. The gradient of ρ(μ) is given by
Thus the intrinsic weighted average of a set of diffusion tensors can be computed by the following procedure:
|input : c1, …, cN P(n)|
|w1, …, wN weights|
|output: μ P(n), the weighted mean|
|μ0 ← I ;|
|i ← 0 ;|
|while ‖ Xi ‖ > e do|
|Xi ← −μi ρ;|
|μi+1 ← Expμi (Xi) ;|
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.
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 (c0, c1, …, cN−1+k−2) that form the spline S(t) which minimizes the energy E. The gradient of E with respect to cj is then given by,
The gradient of the square distance between S(ti) and pi with respect to S(ti) equals,
where LogS(ti) (pi) is the Riemannian logarithmic map (eq. 1). However, LogS(ti) (pi) is a tangent vector at S(ti). Since, the gradient of the energy (eq. 12) is with respect to cj, we need to express the gradient of eq. 13 by using tangent vectors at point cj. Taking this into consideration, equation 13 can be approximated by the formula Λcj (pi, S(ti)) = Logcj (pi) − Logcj (S(ti)), so we obtain
Furthermore, the gradient of S(ti) with respect to cj in the equation 12 is
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
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(ti) is far from the target pi. While S(ti) moves closer to pi during the spline fitting procedure, the error introduced by the approximation of equation 14 moves to zero. By setting a small number e, the outer loop of algorithm 2 will be iterated enough times in order the error of equation 14 to be as small as we need. Hence the control tensors cj, which are obtained as the output of algorithm 2, are estimated without the assumption that there is no curvature in the manifold of SPD matrices. On the contrary, this assumption was used in  and .
|input : N tensors (p0, …, pN−1),|
|N+2(k-1) monotonically increasing|
|knots (t−k+1, …, tN−1+k−1)|
|k, and a small value e|
|output: N-1+k-1 control tensors|
|(c0, …, cN−1+k−2)|
|‖ X0 ‖ ← e+1;|
|for j=0 to N-1+k-2 do|
|Xj ← zero matrix ;|
|for i=0 to N-1 do|
|Xj ← Xj +|
|c′j ← Expcj Xj ;|
|c ← c′ ;|
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, …, tN−1+k−1,M−1+k−1). Note that in this case the knots are vectors with 2 elements, one for each parametric dimension. Finally the new basis functions are formed by the tensor product of 1-dimensional basis functions Ni,j,k([t1t2]) = Ni,k(t1)Nj,k(t2).
Recently, in , 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 . 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 pi 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 . The derivative of this function is and it has the aforementioned properties. By using the above function ϕ, the energy function that we are going to minimize can be written as
The gradient of the energy with respect to the control tensors now becomes
In the above equation the quantity ψ(dist(S(ti),pi)2), weights the given data points pi, leading to a spline approximation that is robust to outliers. The distance function dist(., .), as it was previously mentioned, measures the Riemannian distance between the tensors.
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  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 Xi,j was randomly generated from a Gaussian distribution, for every lattice point of the subsampled tensor field pi,j. These matrices were used as vectors in the tangent space of SPD matrices. Then the noisy tensor field of fig. 3(c) was obtained by applying the Riemannian Exponential map Exppi,j (Xi,j).
Finally we interpolated the noisy tensor fields by a factor of 4, using the four existing methods a) Riemannian geodesic interpolation , b) Log-Euclidean linear interpolation , c) Tangent space interpolation , d) PDE-based anisotropic non-linear diffusion/interpolation , 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.
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 , 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].