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 2013 October 13.
Published in final edited form as:
Med Image Comput Comput Assist Interv. 2012; 15(0 2): 138–145.
PMCID: PMC3796182

Fast Diffusion Tensor Registration with Exact Reorientation and Regularization*


Diffusion tensor imaging is widely used in brain connectivity study. As more and more group studies recruit a large number of subjects, it is important to design registration methods that are not only theoretically rigorous, but also computationally efficient, for processing large data sets. However, the requirement of reorienting diffusion tensors complicates and slows down the registration, especially for those methods whose scalar-image versions have linear complexity, for example, the Demons algorithm. In this paper, we propose an extension of the Demons algorithm that incorporates exact reorientation and regularization into the calculation of deforming velocity, yet preserving its linear complexity. This method restores the computational efficiency of the Demons algorithm to diffusion images, but does not sacrifice registration goodness. In our experiments, the new algorithm achieved state-of-art performance at a ten-fold decrease of computational time.

1 Introduction

Water molecules in biology tissue tend to diffuse faster along, relative to across, obstacle structures. Diffusion tensor imaging (DTI) noninvasively measures this anisotropy (approximated with second-order tensors), providing information about structures such as neural fibers. As more and more large group studies, such as those sponsored by the Human Connnectome Project [1], require nonlinear normalization of diffusion tensor images, it is important to design deformation registration methods that are not only theoretically rigorous, but also computationally efficient, for processing large data sets.

Registration of diffusion tensor images is more complicated than scalar images, because displacing a voxel to a new place does not only change its own diffusion tensor value, but also reorients those of its adjacent voxels [2]. As a result, the deformation force is no longer independent for adjacent voxels, but has a sparse 3D grid structure.

Because reorientation significantly complicates the computation, methods based on rotation-invariant features have been developed. Many of these methods first extract scalar-valued or vector-valued features from tensors, for instance, the fractional anisotropy (FA), and then register diffusion images by aligning these features with multi-channel registration [3,4].

Some other methods directly work on tensors, with similarity metric defined on tensors, and indirect or direct involvement of reorientation into velocity field optimization. Alexander and Gee [5] reorientated tensors according to the new displacement field after each iteration, though not directly accounted the reorientation into optimization. Zhang et al. [6] applied piece-wise affine registration to image subblocks, and then fused these transformations together by smoothing. Since finite-strain (FS) reorientation [2] can be analytically incorporated into affine transformation, the method efficiently estimates the optimal local affine parameters, but it is not clear how the fusion step affects the total registration energy. Cao et al. [7] extended the Large Deformation Diffeomorphic Metric Mapping framework, to analytically embed the preservation-of-principal-direction (PPD) reorientation [2] into the deforming force. In 2009, Yeo et al. 8] embedded exact FS reorientation into the diffeomorphic Demons algorithm [9,10], showing that exact reorientation considerably improved registration accuracy.

Yeo et al.’s work [8] was based on the Demons algorithm, a fast algorithm with O(n) complexity for scalar images, where n is the number of voxels. However, because a large sparse linear system was solved at every iteration, the algorithm became considerably slower. This effect becomes more exacerbated for such an algorithm whose scalar-image version enjoys linear complexity. For example, for our DTI images of size 128×128×128, the method took about 7 hours and more than 20 Giga Bytes (GB) memory on a desktop with an Intel Xeon 2.80 GHz CPU. On the other hand, the regularization on the displacement field was separated as a Gaussian smoothing step after updating the displacement field. This may not be always consistent with the diffeomorphic framework [10].

This raises the question “can the Demons algorithm still enjoy its O(n) complexity when the deformation force is coupled between adjacent voxels and regularization is incorporated?” Though this question is raised from diffusion tensor images, it is generally applicable to other diffusion images, such as high-angular-resolution diffusion images (HARDI) and diffusion spectrum images (DSI). We are particularly interested in the Demons algorithm, because its original linear complexity in [9,10] is well suited for processing large group data sets.

In this paper, we extend the Demons algorithm to incorporate both exact reorientation and regularization into velocity calculation, but without directly solving a large non-separable linear system. This method restores the O(n) computational efficiency of the original Demons algorithm to diffusion images, but does not sacrifice registration goodness, in comparison with solving a large linear system at each iteration. In our experiments, it introduced a 10-fold reduction of the computation time, and achieved state-of-art registration performance.

2 Method

2.1 Diffeomorphic Demons Registration for Diffusion Tensors

Let F and M respectively denote the fixed and moving image, u and v respectively denote the transformation field and velocity field, * denote the compositive operator, •(i) denote the image value at voxel i, •i and •x denote the velocity, displacement or a certain value associated respectively with a voxel i or a point x. The diffeomorphic Demons algorithm for scalar images in [10] is as follows:

Repeat until convergence

  1. For each voxel i, at iteration k, given current transformation field uk:
    1. Define bi(v) = F(i) − M * [uk * exp(v)](i).
    2. Let bi = bi(v) |v=0 = F(i) − M * uk (i), and gi=[partial differential]bi(v)[partial differential]vi|v=0
    3. Let vi=gibigigi+1/δ12, where 1/δi2 is the velocity regularizer at voxel i.
  2. Update displacement field with the velocity field: uk+1 = uk * exp(v).

The algorithm minimizes the energy function

x2130(u)=[F(x)M[composite function (small circle)]u(x)]2dx+vxδx2dx

with the displacement-velocity relationship v=dudt=[partial differential]u[partial differential]t+([nabla]u)v (actually discretized as uk+1 = uk * exp(v)), and a discrete time interval τ = 1. Smoothing can be applied to the displacement field or the velocity field, to impose regularization. The minimization is by approximating the energy function as x2130=i(i[partial differential]bi(v)[partial differential]vjvj)|v=0bi2+iviδi2 with the Gauss-Newton method. For scalar images, [partial differential]bi(v)[partial differential]vj|v=0 equals zero if ij, so the optimization can be separated for each voxel as in Step (1) and has O(n) complexity at each iteration, where n is the number of voxels. For details, please refer to [9,10].

For diffusion tensor images, tensors must be reorientated during transformation to align their directions with the transformed space. Consequently, the displacement at a voxel not only changes its own warped tensor value, but also reorientates those of its adjacent voxels, that is, [partial differential]bi(v)[partial differential]vj|v=0 is not zero for adjacent voxels i and j. In this case, the energy function cannot be separated for each voxel. For details on the formulas of tensor reorientation, please refer to [2].

2.2 Regularization

Standard regularizers for displacement fields include the elastic, diffusion, and curvature regularizers [11] which all can be formatted as the squared norm of a certain derivative value of the displacement field. When discretized, these regularizers can be formulated as i=1nAregiu2 where Areg i is a differential operator at voxel i, represented as a matrix. In this paper, we used ∑l ∫ ‖Hlu2 dx as the regularizer, where Hl is the Hessian operator applied to the lth dimensional component of u. We call it affineness regularizer because the affine transformation is the only non-trivial kernel not punished by this regularizer.

2.3 Local Gauss-Newton Optimization

As shown in Sections 2.1 and 2.2, the energy function of the Demons algorithm, including both the dissimilarity and regularization terms, is the sum of squares of many linear functions. Therefore, the energy function can be written as x2130 = ‖Aimgvbimg2 + ‖Aregv − breg2 = ‖Av − b‖2. Due to tensor reorientation, A is not separable for adjacent voxels, and the optimization cannot be separated as in Step (1) of the Demons algorithm for scalar images. In [8], Yeo minimized the energy by solving this large sparse linear system at each iteration, excluding the regularization term. Though the system is sparse, the computational cost still goes up sharply as the image size increases. For images of size 128×128×128, we observed 20 GB memory usage and about 7 hours computation. This raises the question “can the Demons algorithm still have its O(n) complexity when the deformation force is coupled between adjacent voxels?”

In this paper, we propose a method that incorporates both tensor reorientation and displacement regularization into its optimization iterations, while restores the O(n) complexity of the Demons algorithm. The algorithm is as follows, where we rewrite the energy function as x2130 = ‖Avb2 = ‖(∑Aivi) − b2 and Ai is the columns of A related to the velocity at voxel i:

  1. For each voxel i:
    1. If lHlu2 < < stop threshold, then let di = 0, else let di=argminzAizb2=(A1Ai)1Aib
  2. Let [d1[vertical ellipsis]dn]=[(A1A1)[ddots, three dots, descending](AnAn)]1Ab
  3. Let τ=argmint(Adtb2)=bAdAd2=bAdAimgd2+Aregd2
  4. Let v = dτ.

This algorithm first chooses a descent direction d with local Gauss-Newton optimizations at Step (1a), and then determines the step length τ with a higher-level Gauss-Newton optimization of the energy function given direction d. The algorithm can be understood as the following procedure. Each particle (or voxel) tries to take the shortest path to minimize the global energy as much as possible (Step (1a)), without knowing other particles’ movements. Because their movements change the energy function jointly, rather than independently, a global step length τ is used to coordinate their movement. This strategy is fully compatible with Thirion’s Demons algorithm [9] when A is separable.

At each optimization iteration, we do not choose as the descent direction the steepest direction A′b, for the following reason. The steepest direction A′b sometimes is dominated by large Aib at some voxels (as illustrated in Fig. 1), suppressing other voxels from displacing. Theoretically this phenomenon will disappear as Aib approaching zero. However, in practice, on a discrete grid, this cannot be perfectly achieved. We noticed this phenomenon when we initially tried to use the steepest descent optimization.

Fig. 1
Descent Directions. The steepest can be dominated by certain components.

This local Gauss-Newton strategy converges slower than the Gauss-Newton strategy, but it takes much less computation at each iteration. As shown in Section 2.4, its computational complexity is linear. In our experiments, it achieved equal registration goodness as the Gauss-Newton algorithm in [8], but significantly reduced the registration time and memory usage.

Though this local Gauss-Newton strategy is motivated to solve the problem of diffusion-tensor image registration, its application is not limited to diffusion-tensor images. Generally, it restores the O(n) complexity to the Demons algorithm for situations where the gradient matrix is not separable. It also allows the Demons algorithm to directly incorporate regularization without increasing computational complexity. Application to HARDI image registration, a computation intensive problem, is another interesting extension.

2.4 Computational Complexity

The computational complexity of the algorithm is O(n) where n is the number of voxels of interest. For a particular voxel, its tensor reorientation as well as its local energy function involves only the displacement of the voxels in its neighborhood, so A has a block sparse structure and the number and size of the non-zero blocks in Ai are constant numbers fully determined by image dimension and the regularizer’s differential operator. Therefore, Step (1a) has O(1) complexity, Step (3) has O(n) complexity and the total complexity is O(n).

3 Experiment

3.1 Data Acquisition and Preprocessing

Diffusion weighed images (DWI) of 120 pediatric subjects were collected with 30-direction isotropic DTI sequence (b = 1000s/mm2, voxel size= 2 × 2 × 2 mm3, dimension = 128 × 128 × 128). FSL brain extraction tool (BET) was applied to the B0 images to mask the brain region, with the “R” option turned. Due to the long computation time (about 7 hours) and high memory demand (about 20 GB) of the algorithm in [8], we randomly selected 12 pairs of the subjects for pair-wise registration.

3.2 Image Registration

For the diffeomorphic Demons algorithm proposed by Yeo et al. in [8], its open source implementation, the Tensor Toolkit (TTK) ( was used, with the parameters recommended in [8]: three level multi-scaling and each with ten iterations. For the Local Gauss-Newton Demons algorithm, we applied the affineness regularization to the displacement field. We adjusted the regularization weight to match the displacement smoothness of TTK’s outputs. In Section 3.4, we evaluated the algorithms with three types of harmonic energy, include those that our algorithm did not directly optimize.

3.3 Evaluation Metric

Tensor rooted-mean-square (RMS) difference, FA correlation coefficients (Corr. Coef.), and MD Corr. Coef., between the fixed image and the transformed moving image, were used to indicate registration similarity. Three types of harmonic energy, affineness, curvature and diffusion, were used to indicate displacement smoothness. We used different harmonic energy to give a comprehensive evaluation, avoiding the bias introduced by algorithms’ preference. All these metrics were estimated within the brain masks.

3.4 Results

For each of the registration cases, the two algorithms achieved similar matching between the fixed images and the warped images, with similar regularization energy, while the Local Gauss-Newton algorithm took significantly less computation time and memory usage, as shown in Table 1 and Figures 2 and and3.3. Regarding FA Corr. Coef, both algorithms performed almost equally well; regarding MD Corr. Coef and tensor RMS, the Local Gauss-Newton algorithm performed better, with smoother displacement fields.

Fig. 2
Statistics of Each Registration Case
Fig. 3
Registration Results Example
Table 1
Summary Statistics of the Experiment

This experiment with twelve cases suggests that the Local Gauss-Newton method can achieve similar goodness at similar harmonic energy as the method in [8] does. However, the Local Gauss-Newton algorithm significantly reduced the computation time from about 7 hours to about 45 minutes, and decreased memory demand from 20 GB to 2 GB. The smoother displacement fields of the Local Gauss-Newton algorithm was possibly achieved by directly incorporating the affineness regularization into its optimization of the velocity fields.

4 Conclusion

As more and more brain-connectivity studies recruit a large number of subjects, it is important to design registration methods that are not only theoretically rigorous, but also computationally efficient, for processing large data sets. However, the requirement of reorienting diffusion tensors complicates the calculation of deforming force, and significantly slows computation.

In this paper, we propose a method for the Demons algorithm to exactly incorporate both reorientation and regularization into deforming velocity calculation, but without directly solving a large non-separable linear system. This method restores the O(n) computational efficiency of the original Demons algorithm to diffusion images, without sacrificing registration goodness, in comparison with directly solving a large linear system at each iteration. In our experiments, the new algorithm reduced the computation time and memory demand for ten times.


*This work is supported by grants K01EB013633, R01MH080892 9P41EB015922, P41 RR013642, R01MH71940, U54RR021813, U24RR025736 from NIH.


1. The Human Connectome Project.
2. Alexander DC, Pierpaoli C, Basser PJ, Gee JC. Spatial transformations of diffusion tensor magnetic resonance images. IEEE Transactions on Medical Imaging. 2001;20(11):1131–1139. [PubMed]
3. Park HJ, Kubicki M, Shenton ME, Guimond A, McCarley RW, Maier SE, Kikinis R, Jolesz FA, Westin CF. Spatial normalization of diffusion tensor MRI using multiple channels. Neuroimage. 2003;20(4):1995–2009. [PMC free article] [PubMed]
4. Yang J, Shen D, Davatzikos C, Verma R. Metaxas D, Axel L, Fichtinger G, Székely G. MICCAI 2008, Part II. LNCS, vol. 5242. Springer, Heidelberg; 2008. Diffusion Tensor Image Registration Using Tensor Geometry and Orientation Features; pp. 905–913.
5. Alexander DC, Gee JC. Elastic matching of diffusion tensor images. Computer Vision and Image Understanding. 2000;77(2):233–250.
6. Zhang H, Yushkevich PA, Alexander DC, Gee JC. Deformable registration of diffusion tensor MR images with explicit orientation optimization. Med. Image Anal. 10(5):764–785. [PubMed]
7. Cao Y, Miller M, Mori S, Winslow R, Younes L. Diffeomorphic matching of diffusion tensor images; Conference on Computer Vision and Pattern Recognition Workshop, CVPRW 2006; 2006. Jun, p. 7. [PMC free article] [PubMed]
8. Yeo T, Vercauteren T, Fillard P, Peyrat JM, Pennec X, Golland P, Ayache N, Clatz O. DT-REFinD: Diffusion tensor registration with exact finite-strain differential. IEEE Transactions on Medical Imaging. 2009;28(12):1914–1928. [PMC free article] [PubMed]
9. Thirion JP. Image matching as a diffusion process: an analogy with Maxwell’s demons. Med. Image Anal. 1998;2(3):243–260. [PubMed]
10. Vercauteren T, Pennec X, Perchant A, Ayache N. Diffeomorphic demons: Efficient non-parametric image registration. Neuroimage. 2009;45(1 suppl.):S61–S72. [PubMed]
11. Fischer B, Modersitzki J. A unified approach to fast image registration and a new curvature based registration technique. Linear Algebra and its Applications. 2004;380:107–124.