4TH ORDER DIFFUSION TENSOR INTERPOLATION WITH ...
4TH ORDER DIFFUSION TENSOR INTERPOLATION WITH DIVERGENCE AND CURL ´ CONSTRAINED BEZIER PATCHES Inas Yassine, Tim McGraw Department of Computer Science and Electrical Engineering, West Virginia University ABSTRACT We present a tensor field interpolation method based on tensor-valued B´ezier patches. The control points of the patch are determined by imposing physical constraints on the interpolated field by constraining the divergence and curl of the tensor field. The method generalizes to Cartesian tensors of all orders. Solving for the control points requires the solution of a sparse linear system. Results are presented for order 2 and 4 tensors and comparisons are made with linear interpolation and a previously proposed subdivision scheme. Index Terms— Biomedical image processing, Magnetic resonance imaging, Diffusion processes, Interpolation, Spline functions. 1. INTRODUCTION The Gaussian diffusion process assumed by Basser et al. [1] in the development of diffusion tensor MRI (DT-MRI) is characterized by a symmetric, positive-definite 2nd order tensor. Recently, higher order tensors (order > 2) have been used to ¨ describe diffusion by Ozarslan et al. [2] and Liu et al. [3]. Diffusion within biological structures can be non-Gaussian, for example, when bifurcating or intersecting white-matter fiber bundles occur in the brain. In particular, the 4th order tensor model has recently been studied by Barmpoutis et al. [4], Moakher [5] and others for representing not only diffusivity, but also the orientation distribution function. Tensor interpolation is valuable in tractography and registration algorithms, or any application which may require computation of tensor values between voxels in the image lattice. Many alternatives to componentwise linear interpolation of tensors have been proposed. Geodesic [6, 7, 8], logEuclidean [9], tensor spline [10], and geodesic-loxodrome [11] approaches formulate interpolation in terms of intrinsic distances on some manifold. Some methods have the desirable property of monotonically interpolating some scalar measure, such as determinant [6, 7, 8, 9] or other orthogonal invariants [11]. In this work we propose an interpolation scheme based on minimizing the divergence and curl of the continuous tensor field which interpolates a given set of tensors. “Div-curl splines” [12, 13] use these constraints to impose smoothness on motion vector fields. A harmonic vector
field is both divergence-free and curl-free, and is known to be very smooth. Divergence constraints are commonly used in simulations [14] of incompressible fluid flows. Since these differential operators extend naturally to various spatial dimensions and tensor orders we can develop interpolation and approximation methods in a very general way. An approach to tensor field interpolation based on divergence and curl constraints was presented by Yassine and McGraw [15]. The method was implemented within a subdivision framework, meaning it was not possible to generate interpolated values at arbitrary domain locations. In contrast, the method we present does not suffer from this shortfall. 2. METHODS We begin by presenting the definitions of the divergence and curl operators. The definitions are commonly expressed for the special case, order n = 1 for vector fields. We will present the definitions as generalized for the case n = 4, although the definitions can be extended to higher orders. Tensor divergence. In general, the divergence of an order n tensor field is an order (n − 1) tensor field. For n = 4 the divergence is given in Einstein notation as div D = ∂i Dijkl
(1)
This notation indicates that for all possible values of index i, the tensor components are differentiated with respect to that index and summed over. Note that when the field consists of totally symmetric tensors the divergence tensor is also totally symmetric. Since the diffusion tensors we will be working with are totally symmetric, we can use this constraint to reduce the computation time of our calculations. Tensor curl. The curl of an order n tensor field is an order (n + d − 3) tensor field in d dimensions. For the order 4 case it is defined as curl D = εijk (∂j Dklmn )
(2)
where εijk is the Levi-Civita symbol (permutation tensor) +1 (i, j, k) is an even permutation of indices −1 (i, j, k) is an odd permutation of indices εijk = 0 otherwise. (3)
The curl of an order 4 tensor with minor (or total) symmetry is right minor symmetric (invariant to permutations of the rightmost pair or indices.) 2.1. B´ezier Curves and Surfaces Our spline will take the form of a tensor-valued B´ezier curve. The Bernstein polynomials [16] which form the basis of the B´ezier curve of degree n are given by n i Bin (t) = t (1 − t)n−i (4) i where the binomial coefficients are given by n! n if 0 ≤ i ≤ n i!(n−i)! = i 0 otherwise
(5)
A tensor-valued 3d volume can be defined in terms of the basis functions 4 as D(u, v, w) =
n X n X n X
D(i,j,k) Bin (u)Bjn (v)Bkn (w) (6)
been reshaped into 15(n+1)3 ×1 column vectors c, d and Dx is a sparse 15(n + 1)3 × 15(n + 1)3 matrix. In an analogous fashion, we can define matrices Dy and Dz which can be used to compute the derivative of the spline in the y and z directions. The rows of these matrices can then be used to create matrices which implement the curl and divergence operators. 2.2. Minimization We compute the control points for the tensor-valued patch by solving a system of equations in the least squares sense. The system contains 4 types of equations: 1. Interpolation. For interpolation within a cell, we impose a hard constraint that the 8 tensors at the corners of the cell are known. 2. Boundary conditions. The derivative across the boundary can be controlled to match derivatives with neighboring voxels, or clamp the derivative at the boundary of the dataset.
i=0 j=0 k=0
where D(i,j,k) are the control tensors. The Bernstein form of the B´ezier curve permits the derivatives of a patch to be computed by simply computing the differences of control points. For example, the derivative in the u-direction of the patch 6 is given by ∂u D(u, v, w) = n−1 n X n XX
(7)
(D(i+1,j,k) − D(i,j,k) )Bin−1 (u)Bjn (v)Bkn (w).
i=0 j=0 k=0
3. Divergence minimization. The divergence of the spline can be expressed in terms of the control vertices. Minimization of the divergence can be seen as a physical constraint which favors conservation of mass. 4. Curl minimization. The curl of the spline can also be expressed in terms of the control vertices. These constraints impose additional smoothness on the resulting tensor field.
Continuity across the boundary between adjacent patches C, D can be obtained by imposing the constraint C (n) = B (0) . Smoothness across the boundary between adjacent patches C, D can be obtained by imposing the constraint
The boundary, divergence and curl constraints have associated weights (σbc , σdiv , σcurl ) which permit the user to control the influence of each type of constraint independently. By reshaping d into a column vector the equations can be rearranged in the form
C (n) − C (n−1) = α(D(1) − D(0) )
0 = Ad0 + Bd
1
(8)
with α > 0 for C continuity and α = 1 for G continuity. Note that the derivative operation reduces the degree of the patch by one in the direction being differentiated. We can add and subtract patches by simply adding and subtracting control points as long as the patches being operated on have the same degree in each direction. In computing the divergence and curl of patches we will utilize the degree elevation operation on patches to achieve this condition. A degree n+1 curve with control points E which is equivalent to the given degree n curve with control points D can be obtained by E (i) =
(10)
1
i i D(i−1) + (1 − )D(i) . n+1 n+1
(9)
The differentiation and degree elevation operators for the order 4 tensor valued spline can be implemented as the matrixvector multiplication c = Dx d where the control points have
by moving all of the coefficients of the known interpolated data into A. Both matrices A and B are sparse. The equation is then solved for d in the least squares sense using sparse Cholesky factorization. 3. RESULTS Results for 2nd order tensor interpolation of synthetic data are shown in Figure (1). The polynomial degree, n, was varied from 1 to 7. For n = 1 the B´ezier patch we compute reduces to linear interpolation. Results are also compared with the subdivision approach in [15]. The background color indicates fractional anisotropy (FA). Note the behavior in the bottom row of voxels in Figures (1-7) as the degree of the interpolating curve increases the better that anisotropy is preserved
Fig. 1. Degree n=1 (left) and n=2 (right).
Fig. 2. Degree n=3 (left) and n=4 (right).
Fig. 5. Fourth order tensor interpolation. Degree n=7 (left) and subdivision (right). The background image is generalized anisotropy. Fig. 3. Degree n=5 (left) and n=6 (right). Tesla General Electric Medical Systems Horizon LX imaging system with a diffusion weighted spin echo pulse sequence. Imaging parameters were : effective TR = 9000 ms, TE = 78 ms, NEX = 1. Diffusion-weighted images were acquired with s 25 different gradient directions with b = 1000 mm 2 and a single image was acquired with b ≈ 0. The image field of view was 24×24cm and the acquisition matrix was 256×256×30.
Fig. 4. Degree n=7 (left) and subdivision results (right). 4. CONCLUSIONS during interpolation and the behavior of the subdivision solution is better approximated. In particular, the linear interpolation result has difficulty interpolating between the two tensors with perpendicular dominant orientations. The intermediate result is isotropic. Results for 4th order tensor interpolation of synthetic data are shown in Figure (5). Due to space limitations we only present the results for n = 7 and subdivision. The background image represents generalized anisotropy (GA), computed as described in [2]. The data in Figure (6) were acquired at the Center for Advanced Imaging at WVU on a 3.0
We have presented a scheme for tensor field interpolation which can be extended to tensors of arbitrary order. The method is computationally efficient - It requires only solution of a sparse linear system, and the matrix can be precomputed since it is independent of the data. The interpolation results show that this method is superior to linear interpolation of 4th order tensors in terms of anisotropy preservation. For higher degrees, such as n = 7 our results approximate those of the subdivision approach. Future work will involve the use of these splines for tensor field regularization, tractography and model-based tensor field segmentation.
[5] M. Moakher, “Fourth-order cartesian tensors: old and new facts, notions and applications,” The Quarterly Journal of Mechanics and Applied Mathematics, vol. 61, no. 2, pp. 181, 2008. [6] P.T. Fletcher and S. Joshi, “Riemannian geometry for the statistical analysis of diffusion tensor data,” Signal Processing, vol. 87, no. 2, pp. 250–262, 2007. [7] X. Pennec, P. Fillard, and N. Ayache, “A Riemannian Framework for Tensor Computing,” International Journal of Computer Vision, vol. 66, no. 1, pp. 41–66, 2006. [8] P.G. Batchelor, M. Moakher, D. Atkinson, F. Calamante, and A. Connelly, “A rigorous framework for diffusion tensor calculus.,” Magn Reson Med, vol. 53, no. 1, pp. 221–5, 2005. [9] V. Arsigny, P. Fillard, X. Pennec, and N. Ayache, “Fast and simple calculus on tensors in the Log-Euclidean framework,” Proceedings of MICCAI, pp. 259–267, 2005. Fig. 6. Interpolation of real data taken from sample within the blue box (top). Detail of the sample (bottom left), interpolated field (bottom right) 5. ACKNOWLEDGEMENTS We wish to thank Marc Haut (WVU Health Science Center) for providing the brain data. 6. REFERENCES [1] P.J. Basser, J. Mattiello, R. Turner, and D. Le Bihan, “Diffusion tensor echo-planar imaging of human brain,” in Proceedings of the SMRM, 1993, vol. 56, pp. 584– 562. ¨ [2] E. Ozarslan, B.C. Vemuri, and T.H. Mareci, “Generalized scalar measures for diffusion MRI using trace, variance, and entropy,” Magnetic Resonance in Medicine, vol. 53, no. 4, pp. 866–876, 2005. [3] C. Liu, R. Bammer, B. Acar, and M. E. Moseley, “Characterizing non-Gaussian diffusion by using generalized diffusion tensors,” Magnetic Resonance in Medicine, vol. 51, pp. 924–937, 2004. [4] A. Barmpoutis, B. C. Vemuri, and J. R. Forder, “Fast displacement probability profile approximation from HARDI using 4th-order tensors,” IEEE International Symposium on Biomedical Imaging, pp. 911–914, 14 17 May 2008.
[10] A. Barmpoutis, B.C. Vemuri, T.M. Shepherd, and J.R. Forder, “Tensor Splines for Interpolation and Approximation of DT-MRI With Applications to Segmentation of Isolated Rat Hippocampi,” Medical Imaging, IEEE Transactions on, vol. 26, no. 11, pp. 1537–1546, 2007. [11] G. Kindlmann, R. Estepar, M. Niethammer, S. Haker, and C.-F. Westin, “Geodesic-loxodromes for diffusion tensor interpolation and difference measurement,” in 10th Intl. Conf. on Medical Image Computing and Computer-Assisted Intervention, Brisbane, Australia, October 2007, pp. 1–9. [12] D. Suter, “Vector splines in computer vision.,” in Submitted : Australasian Workshop on Thin-plates, Sydney, 1994. [13] S. Gupta and J. Prince, “Stochastic models for div-curl optical flow methods,” IEEE Signal Processing Letters, vol. 3, no. 2, pp. 32–34, 1996. [14] C.-H. Kim J.-M. Hong, J.-C. Yoon, “Divergenceconstrained moving least squares for fluid simulation,” in Computer Animation and Virtual Worlds, September 2008, pp. 469–477. [15] I. Yassine and T. McGraw, “A subdivision approach to tensor field interpolation,” in Workshop On Computational Diffusion MRI, 2008, pp. 117–124. [16] G. Farin, Curves and surfaces for computer aided geometric design: a practical guide, Academic Press Professional, Inc. San Diego, CA, USA, 1993.
field is both divergence-free and curl-free, and is known to be very smooth. Divergence constraints are commonly used in simulations [14] of incompressible fluid flows. Since these differential operators extend naturally to various spatial dimensions and tensor orders we can develop interpolation and approximation methods in a very general way. An approach to tensor field interpolation based on divergence and curl constraints was presented by Yassine and McGraw [15]. The method was implemented within a subdivision framework, meaning it was not possible to generate interpolated values at arbitrary domain locations. In contrast, the method we present does not suffer from this shortfall. 2. METHODS We begin by presenting the definitions of the divergence and curl operators. The definitions are commonly expressed for the special case, order n = 1 for vector fields. We will present the definitions as generalized for the case n = 4, although the definitions can be extended to higher orders. Tensor divergence. In general, the divergence of an order n tensor field is an order (n − 1) tensor field. For n = 4 the divergence is given in Einstein notation as div D = ∂i Dijkl
(1)
This notation indicates that for all possible values of index i, the tensor components are differentiated with respect to that index and summed over. Note that when the field consists of totally symmetric tensors the divergence tensor is also totally symmetric. Since the diffusion tensors we will be working with are totally symmetric, we can use this constraint to reduce the computation time of our calculations. Tensor curl. The curl of an order n tensor field is an order (n + d − 3) tensor field in d dimensions. For the order 4 case it is defined as curl D = εijk (∂j Dklmn )
(2)
where εijk is the Levi-Civita symbol (permutation tensor) +1 (i, j, k) is an even permutation of indices −1 (i, j, k) is an odd permutation of indices εijk = 0 otherwise. (3)
The curl of an order 4 tensor with minor (or total) symmetry is right minor symmetric (invariant to permutations of the rightmost pair or indices.) 2.1. B´ezier Curves and Surfaces Our spline will take the form of a tensor-valued B´ezier curve. The Bernstein polynomials [16] which form the basis of the B´ezier curve of degree n are given by n i Bin (t) = t (1 − t)n−i (4) i where the binomial coefficients are given by n! n if 0 ≤ i ≤ n i!(n−i)! = i 0 otherwise
(5)
A tensor-valued 3d volume can be defined in terms of the basis functions 4 as D(u, v, w) =
n X n X n X
D(i,j,k) Bin (u)Bjn (v)Bkn (w) (6)
been reshaped into 15(n+1)3 ×1 column vectors c, d and Dx is a sparse 15(n + 1)3 × 15(n + 1)3 matrix. In an analogous fashion, we can define matrices Dy and Dz which can be used to compute the derivative of the spline in the y and z directions. The rows of these matrices can then be used to create matrices which implement the curl and divergence operators. 2.2. Minimization We compute the control points for the tensor-valued patch by solving a system of equations in the least squares sense. The system contains 4 types of equations: 1. Interpolation. For interpolation within a cell, we impose a hard constraint that the 8 tensors at the corners of the cell are known. 2. Boundary conditions. The derivative across the boundary can be controlled to match derivatives with neighboring voxels, or clamp the derivative at the boundary of the dataset.
i=0 j=0 k=0
where D(i,j,k) are the control tensors. The Bernstein form of the B´ezier curve permits the derivatives of a patch to be computed by simply computing the differences of control points. For example, the derivative in the u-direction of the patch 6 is given by ∂u D(u, v, w) = n−1 n X n XX
(7)
(D(i+1,j,k) − D(i,j,k) )Bin−1 (u)Bjn (v)Bkn (w).
i=0 j=0 k=0
3. Divergence minimization. The divergence of the spline can be expressed in terms of the control vertices. Minimization of the divergence can be seen as a physical constraint which favors conservation of mass. 4. Curl minimization. The curl of the spline can also be expressed in terms of the control vertices. These constraints impose additional smoothness on the resulting tensor field.
Continuity across the boundary between adjacent patches C, D can be obtained by imposing the constraint C (n) = B (0) . Smoothness across the boundary between adjacent patches C, D can be obtained by imposing the constraint
The boundary, divergence and curl constraints have associated weights (σbc , σdiv , σcurl ) which permit the user to control the influence of each type of constraint independently. By reshaping d into a column vector the equations can be rearranged in the form
C (n) − C (n−1) = α(D(1) − D(0) )
0 = Ad0 + Bd
1
(8)
with α > 0 for C continuity and α = 1 for G continuity. Note that the derivative operation reduces the degree of the patch by one in the direction being differentiated. We can add and subtract patches by simply adding and subtracting control points as long as the patches being operated on have the same degree in each direction. In computing the divergence and curl of patches we will utilize the degree elevation operation on patches to achieve this condition. A degree n+1 curve with control points E which is equivalent to the given degree n curve with control points D can be obtained by E (i) =
(10)
1
i i D(i−1) + (1 − )D(i) . n+1 n+1
(9)
The differentiation and degree elevation operators for the order 4 tensor valued spline can be implemented as the matrixvector multiplication c = Dx d where the control points have
by moving all of the coefficients of the known interpolated data into A. Both matrices A and B are sparse. The equation is then solved for d in the least squares sense using sparse Cholesky factorization. 3. RESULTS Results for 2nd order tensor interpolation of synthetic data are shown in Figure (1). The polynomial degree, n, was varied from 1 to 7. For n = 1 the B´ezier patch we compute reduces to linear interpolation. Results are also compared with the subdivision approach in [15]. The background color indicates fractional anisotropy (FA). Note the behavior in the bottom row of voxels in Figures (1-7) as the degree of the interpolating curve increases the better that anisotropy is preserved
Fig. 1. Degree n=1 (left) and n=2 (right).
Fig. 2. Degree n=3 (left) and n=4 (right).
Fig. 5. Fourth order tensor interpolation. Degree n=7 (left) and subdivision (right). The background image is generalized anisotropy. Fig. 3. Degree n=5 (left) and n=6 (right). Tesla General Electric Medical Systems Horizon LX imaging system with a diffusion weighted spin echo pulse sequence. Imaging parameters were : effective TR = 9000 ms, TE = 78 ms, NEX = 1. Diffusion-weighted images were acquired with s 25 different gradient directions with b = 1000 mm 2 and a single image was acquired with b ≈ 0. The image field of view was 24×24cm and the acquisition matrix was 256×256×30.
Fig. 4. Degree n=7 (left) and subdivision results (right). 4. CONCLUSIONS during interpolation and the behavior of the subdivision solution is better approximated. In particular, the linear interpolation result has difficulty interpolating between the two tensors with perpendicular dominant orientations. The intermediate result is isotropic. Results for 4th order tensor interpolation of synthetic data are shown in Figure (5). Due to space limitations we only present the results for n = 7 and subdivision. The background image represents generalized anisotropy (GA), computed as described in [2]. The data in Figure (6) were acquired at the Center for Advanced Imaging at WVU on a 3.0
We have presented a scheme for tensor field interpolation which can be extended to tensors of arbitrary order. The method is computationally efficient - It requires only solution of a sparse linear system, and the matrix can be precomputed since it is independent of the data. The interpolation results show that this method is superior to linear interpolation of 4th order tensors in terms of anisotropy preservation. For higher degrees, such as n = 7 our results approximate those of the subdivision approach. Future work will involve the use of these splines for tensor field regularization, tractography and model-based tensor field segmentation.
[5] M. Moakher, “Fourth-order cartesian tensors: old and new facts, notions and applications,” The Quarterly Journal of Mechanics and Applied Mathematics, vol. 61, no. 2, pp. 181, 2008. [6] P.T. Fletcher and S. Joshi, “Riemannian geometry for the statistical analysis of diffusion tensor data,” Signal Processing, vol. 87, no. 2, pp. 250–262, 2007. [7] X. Pennec, P. Fillard, and N. Ayache, “A Riemannian Framework for Tensor Computing,” International Journal of Computer Vision, vol. 66, no. 1, pp. 41–66, 2006. [8] P.G. Batchelor, M. Moakher, D. Atkinson, F. Calamante, and A. Connelly, “A rigorous framework for diffusion tensor calculus.,” Magn Reson Med, vol. 53, no. 1, pp. 221–5, 2005. [9] V. Arsigny, P. Fillard, X. Pennec, and N. Ayache, “Fast and simple calculus on tensors in the Log-Euclidean framework,” Proceedings of MICCAI, pp. 259–267, 2005. Fig. 6. Interpolation of real data taken from sample within the blue box (top). Detail of the sample (bottom left), interpolated field (bottom right) 5. ACKNOWLEDGEMENTS We wish to thank Marc Haut (WVU Health Science Center) for providing the brain data. 6. REFERENCES [1] P.J. Basser, J. Mattiello, R. Turner, and D. Le Bihan, “Diffusion tensor echo-planar imaging of human brain,” in Proceedings of the SMRM, 1993, vol. 56, pp. 584– 562. ¨ [2] E. Ozarslan, B.C. Vemuri, and T.H. Mareci, “Generalized scalar measures for diffusion MRI using trace, variance, and entropy,” Magnetic Resonance in Medicine, vol. 53, no. 4, pp. 866–876, 2005. [3] C. Liu, R. Bammer, B. Acar, and M. E. Moseley, “Characterizing non-Gaussian diffusion by using generalized diffusion tensors,” Magnetic Resonance in Medicine, vol. 51, pp. 924–937, 2004. [4] A. Barmpoutis, B. C. Vemuri, and J. R. Forder, “Fast displacement probability profile approximation from HARDI using 4th-order tensors,” IEEE International Symposium on Biomedical Imaging, pp. 911–914, 14 17 May 2008.
[10] A. Barmpoutis, B.C. Vemuri, T.M. Shepherd, and J.R. Forder, “Tensor Splines for Interpolation and Approximation of DT-MRI With Applications to Segmentation of Isolated Rat Hippocampi,” Medical Imaging, IEEE Transactions on, vol. 26, no. 11, pp. 1537–1546, 2007. [11] G. Kindlmann, R. Estepar, M. Niethammer, S. Haker, and C.-F. Westin, “Geodesic-loxodromes for diffusion tensor interpolation and difference measurement,” in 10th Intl. Conf. on Medical Image Computing and Computer-Assisted Intervention, Brisbane, Australia, October 2007, pp. 1–9. [12] D. Suter, “Vector splines in computer vision.,” in Submitted : Australasian Workshop on Thin-plates, Sydney, 1994. [13] S. Gupta and J. Prince, “Stochastic models for div-curl optical flow methods,” IEEE Signal Processing Letters, vol. 3, no. 2, pp. 32–34, 1996. [14] C.-H. Kim J.-M. Hong, J.-C. Yoon, “Divergenceconstrained moving least squares for fluid simulation,” in Computer Animation and Virtual Worlds, September 2008, pp. 469–477. [15] I. Yassine and T. McGraw, “A subdivision approach to tensor field interpolation,” in Workshop On Computational Diffusion MRI, 2008, pp. 117–124. [16] G. Farin, Curves and surfaces for computer aided geometric design: a practical guide, Academic Press Professional, Inc. San Diego, CA, USA, 1993.
