Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Med Image Comput Comput Assist Interv. Author manuscript; available in PMC 2010 April 29.
Published in final edited form as:
Med Image Comput Comput Assist Interv. 2009 October 1; 5761: 640–647.
doi:  10.1007/978-3-642-04268-3_79
PMCID: PMC2861508

Groupwise Registration and Atlas Construction of 4th-Order Tensor Fields Using the R+ Riemannian Metric*


Registration of Diffusion-Weighted MR Images (DW-MRI) can be achieved by registering the corresponding 2nd-order Diffusion Tensor Images (DTI). However, it has been shown that higher-order diffusion tensors (e.g. order-4) outperform the traditional DTI in approximating complex fiber structures such as fiber crossings. In this paper we present a novel method for unbiased group-wise non-rigid registration and atlas construction of 4th-order diffusion tensor fields. To the best of our knowledge there is no other existing method to achieve this task. First we define a metric on the space of positive-valued functions based on the Riemannian metric of real positive numbers (denoted by R+). Then, we use this metric in a novel functional minimization method for non-rigid 4th-order tensor field registration. We define a cost function that accounts for the 4th-order tensor re-orientation during the registration process and has analytic derivatives with respect to the transformation parameters. Finally, the tensor field atlas is computed as the minimizer of the variance defined using the Riemannian metric. We quantitatively compare the proposed method with other techniques that register scalar-valued or diffusion tensor (rank-2) representations of the DWMRI.

1 Introduction

Group-wise image registration is a challenging task in medical imaging which is related to the problem of computing an atlas, i.e. the image of the average subject from a set of co-registered subjects. There are two prevalent approaches for atlas construction. The first one is based on group-wise alignment of 3D shapes [1,2], while the second one is uses alignment of 3D image intensity maps.

In this paper we focus on the second category, and therefore we review only techniques that are based on intensity map registration. Joshi et al. [3] proposed a method for group-wise image registration and simultaneous atlas construction. In this method the atlas is formed by minimizing the distance between the displacement fields that warp the images and therefore it is not biased toward a specific subject data. The estimated atlas does not belong to the set of registered subjects unlike the method presented in [4], which perform pair-wise registration of all the subjects and select the least biased target as the atlas.

The aforementioned methods perform scalar-valued image registration. It has been shown, however, that registration of diffusion tensor-valued images (DTI) produces more accurate alignments of fibrous tissues [5]. In this approach the tensors should be re-oriented appropriately after the warping of the DTI images in order to preserve the micro-structural geometry in the subjects. One way to avoid the tensor re-orientation is to register rotation invariant quantities or other highly structured features extracted from DTI [6]. A DTI similarity measure that uses the full information in the tensors and performs their re-orientation using locally affine transformations was employed in [7]. Furthermore, two methods for diffeomorphic non-rigid DTI registration were proposed in [8] and [9] both of which use analytic derivatives of the reorientation term in the corresponding energy functions.

All the above techniques perform pair-wise DTI registration. Multi-subject registration for DTI atlas construction was proposed in [10] by extending the scalar-image framework in [3]. Another group-wise DTI registration technique which unfolds the manifold described by the Geodesic-Loxodromes metric on diffusion tensors and produces vector-valued images that are being warped in order to estimate the DTI atlas was recently proposed in [11].

Although the methods for DTI registration and atlas construction yield richer representations than the corresponding scalar-image based techniques, they fail in regions of fiber crossings and other complex tissue geometries since 2nd-order tensors cannot account for multiple peaks in the diffusivity function. This problem can be resolved by using 4thorder tensor fields and registering them using the recently proposed method in [12]. In their work, it was shown that the alignment of 4thorder tensor fields produces more accurate results compared to those obtained by DTI registration. This technique performs tensor comparison using Hellinger's distance, which is however defined between probabilities. Since diffusion tensors are not probabilities, Hellinger's distance is not a suitable measure, unless we perform tensor normalizations which are unnatural and we avoid in this paper. Furthermore, the method in [12] performs pair-wise tensor field registration and hence cannot be directly employed for group-wise registration or statistical atlas construction.

In this paper we present a novel method for unbiased 4th-order tensor field atlas construction. Our method (significantly) generalizes the unbiased diffeomorphic scalar image atlas construction framework in [3] to the case of symmetric positive definite higher-order tensors. The atlas is computed simultaneously with the non-rigid deformation fields using a functional minimization procedure. We define a novel cost function using the Riemannian metric on positive valued functions which is a generalization of the Riemannian metric on R+. This metric appropriately handles the positive nature of the symmetric positive-definite high-order tensors and their re-orientation is performed analytically using the Gram-Schmidt orthogonalization process of the local Jacobian matrices. The method is validated using synthetic and real DW-MRI data from isolated human hippocampi.

The key contributions of this work are: To the best of our knowledge, this is the first report in literature for higher-order tensor field atlas construction. Our method outperforms the existing methods that register derived scalar images or 2nd-order tensor fields from DWMRI, both of which fail to accurately warp datasets with complex local tissue structures such as fiber crossings. Furthermore, we employ a novel metric based on the Riemannian geometry of positive-valued spherical functions and we show that it produces more accurate results compared to the standard Euclidean metric. Finally our cost function has analytic derivatives with respect to the unknown transformation parameters that lead to an efficient and easily scalable implementation of our framework.

2 Riemannian Metric for Positive-Valued Real Functions

Assume a, b [set membership] R+, i.e. are elements of the space of positive real numbers. The Logarithmic map at location a is given by Loga(x) = log(x/a) and corresponds to the local tangent vector toward x. Its inverse function is the Exponential map, which is given by Expa(t) = exp(t)a and projects the tangent t [set membership] R back to the space R+. The corresponding Riemannian distance between a and b [set membership] R+ is given by the length of the tangent


which satisfies scale invariance, i.e. dist(sa, sb) = dist(a, b) ∀, b, s [set membership] R+, additionally to the properties of distance measures.

The Riemannian metric in R+ can also be used to define distances between positive-valued functions fa(x) and fb(x) x [set membership] Ω as follows: dist2(fa, fb) = ∫Ω dist2(fa(x), fb(x))dx. In the particular case of parametric spherical functions d(g; D1) and d(g; D2), where g [set membership] S2 and D1 and D2 are the corresponding parameter vectors, the distance is given by


Note that the integral in Eq. 2 is over S2, i.e. the space of unit vectors g. This distance function is invariant with respect to 3D rotations and scale, i.e. dist(sR [composite function (small circle)] D1, sR [composite function (small circle)] D2) = dist(D1, D2) ∀s [set membership] R+ and R [set membership] SO3.

Similarly, the distance between ordered n-tuples whose elements are positive real numbers can be defined using the Riemannian metric in R+. In this case the distance between A = {a1, a2, … , an} and B = {b1, b2, … , bn} ai, bi [set membership] R+ is given by dist2(A,B)=i=1Ndist2(ai,bi). This can also be seen as a discrete approximation of Eq. 2 by taking ai = d(gi; D1) and bi = d(gi; D2), where gi is a predefined set of vectors in S2.

In the next section we will employ the above distance measure in order to achieve simultaneous group-wise registration and atlas construction of fields of spherical functions modeled using Cartesian tensor bases of order 4.

3 Groupwise Registration of 4th-Order Tensor Fields

Cartesian tensor bases of various orders have been used for approximating physical quantities computed from DW-MRI datasets. 4th -order tensors d(g; D) = Σi,j,k,l Di,j,k,l gigjgkgl have been employed to approximate the diffusivity function in generalized diffusion tensor images [13], and the kurtosis component of the diffusion in diffusion kurtosis images [14].

In the case of 4th-order generalized diffusion tensors, the diffusivity is a positive-valued function and can be computed using the parametrization in [15]. This produces fields of positive-valued spherical functions whose processing can be achieved using the Riemannian metric presented in Sec. 2.

The problem of group-wise registration of N tensor-fields and simultaneous atlas estimation can be formulated as a functional minimization problem. By using Eq. 2 the energy function to be minimized is given by


where Dμ is the 4th-order tensor coefficients of the estimated atlas, [var phi]n is the estimated deformation to be applied to the nth tensor field, and cost() is a cost function that adds smoothing constraints to the estimated deformations.

Note that the tensor coefficients are dependent on the local rotation of the coordinate system [12]. Hence, given a deformation [var phi]n the transformed spherical function field at location x can be computed as


where Rx is the rotation of deformation [var phi]n at location x, and the notation (Rxg)i represents the ith component of the rotated vector g.

The deformation can be parametrized as a time varying vector field such that [partial differential][var phi]n(x,t)/[partial differential]t = vn(x, t), t [set membership] [0, 1], where vn(x, t) is the velocity field at time t. In this formulation the estimated deformation is given by φn=φn(x,1)=01vn(x,t)dt. Furthermore, the cost() function in Eq. 3 can be defined as 01Lvn(x,t)2dt, where L is a differential operator on the velocity fields [3].

We will minimize the energy function (Eq. 3) by evolving the deformation fields [var phi]n using a greedy iterative scheme which approximates the solution to the above minimization problem, similar to the technique in [3]. For this purpose we will construct a field of forces by computing the first order variation of the first term in Eq. 3 with respect to the transformation parameters as follows


where the variation [nabla]trans is related with the local translation (i.e. variation of Dni,j,k,l(xφn) in Eq. 4) and [nabla]rot is related with the local rotation (i.e. variation of (Rxg)i(Rxg)j(Rxg)k(Rxg)l in Eq. 4). The computation of these terms is discussed in Sec. 3.1.

After the estimation of the fields of forces Fn, n = 1 … N we compute the update vector fields vn = ∫Ω K(x)Fn(x)dx, where K is a kernel applied to the field of forces. In our experiments we employed the kernel K(x) = η(x)G(x), where G is a Gaussian kernel centered at x and η is a smooth function that takes zero value at the boundaries and therefore imposes zero boundary conditions on the kernel K as was done in [8]. Note that the integration of K with Fn is a convolution that becomes multiplication in the frequency domain, hence it can be efficiently computed using the discrete Fourier transform [16]. Then, the deformation fields are updated as φnnew=φnold(x+εvn) using a small step ε.

Finally, the tensor coefficients of the atlas can be updated by also minimizing the first term in Eq. 3 with respect to the parameters of a positive definite 4th-order tensor using the parametrization in [15].

3.1 Implementation Details

In general, the integral over the sphere in Eq. 5 cannot be computed analytically when the Cartesian tensor parametrization is used for modeling the diffusivity function. On the other hand the Riemannian space of ordered n-tuples (see Sec. 2) leads to analytic calculations and therefore we used it in our implementation. We constructed an m-tuple space by using a set of unit vectors gm m = 1 … M uniformly distributed on the sphere. This set of vectors can be constructed by tessellating the icosahedron and then projecting the vectors on the unit hemisphere (we consider only a hemisphere due to antipodal symmetry of diffusivity functions). We use this set of vectors in order to evaluate the spherical functions In,m = log(d(gm; Dn [composite function (small circle)] [var phi]n)), m = 1 … M and n = 1 … N. This creates N vector valued images In, whose vectors contain the M elements of the m-tuples. Note that in this m-tuple space the integrals over the sphere in Eqs. 3 and 5 become summations over m.

The above discretization helps us also in reducing the time complexity of atlas computation, which can now be efficiently computed by


where dμ(gm) is also in the form of a vector valued image, whose vectors contain M elements. Note that log(d(gm; Dn[composite function (small circle)][var phi]n)) is an already computed image (In, m), and therefore there is no need to re-deform the images and re-compute the log maps. The corresponding driving forces in Eq. 5 are now computed as follows


where Lm,n=log(In,m(x)d(gm;Dμ(xφn))), [nabla] In, m is simply the spatial gradient of a scalar valued image and the second term in Eq. 7 correspond to the gradient related to the tensor re-orientation. In this term the rotation Ry at location ‘y’ can be efficiently computed by the Gram-Schidt algorithm as in [8]. Using this orthogonalization technique the components of the rotation matrix are expressed as functions of the displacement vectors in [var phi]n, hence we can easily compute analytic derivatives with the unknown transformation parameters denoted as [nabla]Rygm. The computed derivatives are non zero for those voxels ‘y’ which are in the neighborhood of our current voxel ‘x’. Furthermore, the gradient of the tensor with respect to the rotation is given by Rygmd(gm;Dn(yφn))=4i,j,kDni,j,k,l(yφn)(Ryg)i(Ryg)j(Ryg)k.

Finally, after the termination of the iterative minimization procedure, the 4th-order tensor coefficients can be computed by fitting the tensorial model to the estimated values dμ(gm) using the positive-definite parametrization in [15].

4 Experimental Results

In order to compare the Riemannian metric presented in Sec. 2 with a Euclidean metric in terms of fiber orientation accuracy of the atlas estimated by each metric, we performed the following experiment. We synthesized a 2-fiber crossing DW-MRI dataset (in a single voxel) using the realistic adaptive kernel model shown in Fig.3 of [17] (81 gradient directions and b = 1250s/mm2). We computed a 4th-order tensor (shown in Fig. 1 upper left) from the synthetic dataset using the algorithm in [15]. Then we generated 100 more datasets by applying small rotations to the simulated crossing and by adding outliers (few of them are shown in Fig. 1 left). The computed atlases (average tensors) are compared in the bar chart of Fig. 1. As expected, the Riemannian mean outperforms the Euclidean mean since the physical space of the data is that of positive-valued functions.

Fig. 1
Comparison of the 4th-order tensor atlases computed by various metrics: a) Euclidean mean, b) Riemannian mean (computed in the space presented in Sec. 2)

To motivate the use of 4th-order tensors in registering DW-MRI, we also simulated a fiber crossing dataset and synthesized a deformation field (Fig. 2). Then we computed the corresponding FA, DTI and 4th-order tensor fields and their deformed images as well. We registered the obtained datasets using the scalar image registration method in [3], its DTI modification [10], and the 4th-order tensor field algorithm in [12] respectively as well as our proposed method. After that, the displacement field produced by each algorithm was used to warp the deformed 4th-order tensor field and it was then compared to the ground truth field shown in Fig. 2(left) using Eq. 2. The results demonstrate that our method produced more accurate mappings and registered successfully the data.

Fig. 2
Comparison of registration methods using a synthetic fiber crossing dataset. The errors were measured by evaluating Eq. 2 on the whole field (blue).

Finally, we computed the 4th-order tensor field atlas from four hippocampal datasets. Each dataset consists of 21 diffusion-weighted images collected with a 415 mT/m diffusion gradient (Td = 17 ms, δ = 2.4ms, b = 1250 s/mm2). Figure 3 shows the original misalignment of the corresponding S0 images and the aligned images after applying our method. The 4th-order tensor field atlas is depicted at the bottom of this figure and contains all the known structures of the hippocampal anatomy. The variations in the dataset can be explored by observing the standard deviation field computed by the proposed Riemannian metric (shown in an ROI on the bottom left).

Fig. 3
Real datasets from hippocampus before and after alignment using our method. The constructed 4th-order tensor field atlas is shown at the bottom. The field of standard deviations can show the variations in the dataset.

5 Conclusions

In this paper we presented a novel groupwise registration and atlas construction algorithm for DWMRI data sets each of which is represented by a 4th order tensor field. To the best of our knowledge, there is no existing literature on this topic. The key contribution of this work is the definition of a novel metric for positive valued spherical functions which was then used in the unbiased group-wise registration and atlas construction. Experimental results on comparisons with scalar and DTI registration techniques are favourable to our method.


*This research was in part funded by the NIH grant EB007082 to BCV.


1. De Craene M, du Bois d'Aische A, Macq B, Warfield SK. Multi-subject registration for unbiased statistical atlas construction. In: Barillot C, Haynor DR, Hellier P, editors. MICCAI 2004 LNCS. Vol. 3216. Springer; Heidelberg: 2004. pp. 655–662.
2. Yeo BTT, Sabuncu MR, Desikan R, Fischl B, Golland P. Effects of registration regularization and atlas sharpness on segmentation accuracy. Medical Image Analysis. 2008;12:603–615. [PMC free article] [PubMed]
3. Joshi S, Davis B, Jomier A, Gerig G. Unbiased diffeomorphic atlas construction for computational anatomy. NeuroImage. 2004;23:151–160. [PubMed]
4. Park H, Bland PH, Hero AO, III, Meyer CR. Least biased target selection in probabilistic atlas construction. In: Duncan JS, Gerig G, editors. MICCAI 2005 LNCS. Vol. 3750. Springer; Heidelberg: 2005. pp. 419–426. [PubMed]
5. Alexander DC, Pierpaoli C, Basser PJ, Gee JC. Spatial transformations of diffusion tensor magnetic resonance images. TMI. 2001;20(11):1131–1139. [PubMed]
6. Ruiz-Azola J, Westin CF, Warfield SK, Alberola C, Maier S, Kikinis R. Non rigid registration of 3d tensor medical data. Med Im Anal. 2002;6:143–161. [PubMed]
7. Zhang H, Yushkevich P, Gee J. Registration of DTI. CVPR. 2004;1:842–847.
8. Cao Y, Miller M, Mori S, Winslow R, Younes L. Diffeomorphic matching of diffusion tensor images. Computer Vision and Pattern Recognition Workshop. 2006 June;67 2006. [PMC free article] [PubMed]
9. Yeo B, Vercauteren T, Fillard P, Pennec X, Gotland P, Ayache N, Clatz O. DTI registration with exact finite-strain differential. ISBI. 2008:700–703.
10. Zhang H, Yushkevich PA, Rueckert D, Gee JC. Unbiased white matter atlas construction using diffusion tensor images. In: Ayache N, Ourselin S, Maeder A, editors. MICCAI 2007, Part II LNCS. Vol. 4792. Springer; Heidelberg: 2007. pp. 211–218. [PubMed]
11. Irfanoglu MO, Machiraju R, Sammet S, Pierpaoli C, Knopp MV. Automatic deformable diffusion tensor registration for fiber population analysis. In: Metaxas D, Axel L, Fichtinger G, Sz'ekely G, editors. MICCAI 2008, Part II LNCS. Vol. 5242. Springer; Heidelberg: 2008. pp. 1014–1022. [PubMed]
12. Barmpoutis A, Vemuri BC, Forder JR. Registration of high angular resolution diffusion mri images using 4th order tensors. In: Ayache N, Ourselin S, Maeder A, editors. MICCAI2007, Part I LNCS. Vol. 4791. Springer; Heidelberg: 2007. pp. 908–915. [PMC free article] [PubMed]
13. Özarslan E, Mareci TH. Generalized diffusion tensor imaging and analytical relationships between dti and hardi. MRM. 2003;50(5):955–965. [PubMed]
14. Jensen JH, Helpern JA, Ramani A, Lu H, Kaczynski K. DKI: the quantification of non-gaussian water diffusion by means of MRI. MRM. 2005;53(6):1432–1440. [PubMed]
15. Barmpoutis A, Hwang MS, Howland D, Forder JR, Vemuri BC. Regularized positive-definite fourth order tensor field estimation from DW-MRI. NeuroImage. 2009;45(1 supl 1):153–162. [PMC free article] [PubMed]
16. Joshi S, Grenander U, Miller MI. On the geometry and shape of brain sub-manifolds. IJPRAI. 1997;11(8):1317–1343.
17. Barmpoutis A, Jian B, Vemuri BC. IPMI 2009 LNCS. Vol. 5636. Springer; Heidelberg: 2009. Adaptive kernels for multi-fiber reconstruction; pp. 338–349. [PMC free article] [PubMed]