Non-Linear Spatial Normalization of High Angular Resolution

Document Sample
Non-Linear Spatial Normalization of High Angular Resolution Powered By Docstoc
					       Non-Linear Spatial Normalization of High Angular
        Resolution Diffusion Imaging Data using Fiber
                  Orientation Distributions

       David Raffelt1,2, J-Donald Tournier3,4, Jurgen Fripp1, Stuart Crozier2, Alan
                            Connelly3,4and Olivier Salvado1
    The Australian E-Health Research Centre, CSIRO, Brisbane, QLD, Australia, 2Department
   of Biomedical Engineering, University of Queensland, Brisbane, QLD, Australia, 3Brain
   Research Institute, Florey Neuroscience Institutes (Austin), Melbourne, VIC, Australia,
         Department of Medicine, University of Melbourne, Melbourne, VIC, Australia

       Abstract. Using high angular resolution diffusion weighted imaging data, fiber
       orientation distributions (FODs) can be computed to provide information on the
       structure and organization of white matter architecture. In order to investigate
       differences between subjects, we present a method for non-rigid registration of
       FODs that incorporates a novel technique for correcting the orientation of each
       FOD based on a local affine model of the non-rigid transformation. The
       registration optimizes directly on the FODs, and is independent of the FOD
       estimation method used (Spherical Deconvolution or Q-Ball Sharpening
       Deconvolution Transform). Dataset from ten healthy young adults were used
       for qualitative and quantitative assessment of the spatial normalization results.
       Unlike existing techniques that employ scalar or diffusion tensor information,
       our method exploits additional information on the orientation of crossing fibers
       to provide more accurate spatial correspondence.

1 Introduction

Diffusion-Weighted Imaging (DWI) provides a method to investigate the architecture
of oriented tissue such as white matter (WM) in vivo. Techniques such as Spherical
Deconvolution [1] and Q-Ball Sharpening Deconvolution Transform [2], provide a
means of estimating Fiber Orientation Distributions (FODs) using High Angular
Resolution Diffusion-weighted Imaging (HARDI) [3]. FODs can then be employed in
tractography to reveal the trajectories and connectivity of white matter fibers in the
   Spatial normalisation is an important step required to compare white matter
between individuals and across groups. Previous methods have focused on
coregistering DWI data to high resolution T1 weighted (T1W) images that are then
registered to a common template. However, T1W images contain no structural
information within the homogenous intensity of white matter and therefore the
registration of deep white matter is unguided. An additional problem is introduced by
the geometric distortions often present in DWI that may need to be corrected to allow
accurate coregistration.
   Other methods [4][5] have utilized invariant scalar information, such as Fractional
Anisotropy (FA), obtained from the diffusion tensor (DT) [6]. As different white
matter structures contain heterogeneous FA values this contrast is more informative
than T1W images. A number of other methods have been proposed whereby
registration is performed directly on the DT itself and therefore additional orientation
information is included in the registration [7][8]. These methods are more
complicated due to the required re-orientation of the diffusion tensor to ensure the
principal eigenvector remains aligned with the containing fiber bundle [9]. The
registration and normalization of the DT improves the detection of white matter
differences [10], and provides a means to perform group average fiber tracking [4].
   A well known limitation of the DT is its inability to adequately model crossing
fibers. At the current DWI resolution it has been suggested that up to one third of
white matter voxels contain more than a single fiber population [11]. Higher order
DTs can be used to better describe the complex diffusion propagator in regions with
crossing fibres [12]. Barmpoutis et al. [13] demonstrated that registering 4th order
DTs gave more accurate spatial normalisation in regions with two fibre populations
compared to the traditional rank-2 DT.
   With models such as higher order tensors and diffusion ODFs, it is often assumed
that the orientation of underlying fibre populations correspond to peaks in diffusion
profile. However, this has been shown to be true only for fibres crossing at 90
degrees, and at all other angles the orientation of the diffusion maxima deviates from
the true fibre orientation [14][15]. Recent work has demonstrated that the sharper
peaks of the FOD, that correspond to the true fibre directions, improve the spatial
localization of fibres generated by tractography [2].
   In this paper we describe a method to spatially normalize HARDI data using a non-
rigid registration method that optimizes directly on FODs and therefore, unlike
previous techniques, exploits the additional information on fibre orientations provided
in voxels with crossing fibres. In order to correct the FOD orientations we present a
novel method for transforming FODs based on the local affine transform.

2 Method

2.1 FOD Transformation

The transformation of scalar-valued images involves mapping the intensity at a point
x through some transform T to a new location x’.
                                  x' = T ( x)                                       (1)
Like DTs, transforming FODs involves correcting orientation information to ensure
corresponding fiber populations of neighbouring voxels are similarly aligned. This is
performed by modifying the FOD at point x in the original image using an affine
transform F. When T(x) is a rigid or affine transform F is independent of x. When
T(x) is a non-rigid transformation, a local affine model can be computed at each point

x by F = I + Ju where I is the identity matrix and Ju is the Jacobian matrix of the
displacement vector field u at point x [9].
   FODs are modeled by a series of NSH spherical harmonic (SH) coefficients up to a
maximum harmonic degree l (lmax). In order to correct the FOD orientation using F,
we approximate the FOD by the weighted sum of SH delta functions (with an
identical lmax) along a set of n > NSH vectors equally distributed over a hemisphere
(FOD symmetry is assumed).
                           fl ≈   ∑ w δ (θ , φ ) = Πw
                                   i =0
                                           i    l       i   i                                    (2)

where:    fl is a vector of SH coefficients representing the FOD with lmax of l; wi is the
weight of the delta function along vector I;    δ l (θ i , φi ) is a vector containing the SH
coefficients for a delta function along orientation i, up to order l; w is a vector
containing w1 to wn ; and Π = [ δl (θ1 , φ1 ) , δl (θ 2 ,φ2 ) , ... , δl (θ n , φn ) ]

Each SH coefficient     klm of the delta function vector δ l (θ i , φi ) is computed using the
harmonic expansion:

                              klm = ∫ f ( s ) ylm ( s )ds                                        (3)

Where:     s is the surface of the sphere, f is a delta function along orientation i, and
y l is   the real valued spherical harmonic function of degree l and order m.

Prior to transforming each FOD the following steps are performed:
   1. Form {v i } a set of n equally distributed 3D unit vectors
   2. Convert {v i } to spherical coordinates θ , φ
   3. Form matrix Π using the orientations from {v i }
   4. Calculate the pseudo inverse of Π :

                                  Π −1 = Π T Π Π T  (           )

To transform the FOD:
  5. Compute the delta function weights w to fit the FOD               fl :

                                          w = Π −1f l                                            (5)

    6. Transform each of the vectors           {v i } using the local affine transformation F:

                                           v ' = Fv i
                                           ˆ      ˆ                                              (6)
  7. Generate   Π' by performing steps 2 and 3 using {v i '}
  8. Calculate the transformed FOD:

                                  f l ' = Π' w                                      (6)

  9. Repeat steps 5 to 8 on all voxels in the moving image within a tissue mask.

Note that we make the assumption that the delta function weights are invariant to the
transformation in order to preserve the volume fraction of each fiber population.

2.2 Registration Method

The initial step in the registration involves performing a rigid and then affine
transformation. The affine transform is subsequently set as a bulk transform in a B-
spline based free form deformation (FFD) non-rigid registration [16]. Each step was
implemented in the Insight Segmentation and Registration Toolkit (ITK) framework
   Optimization was performed by regular step gradient descent and sought to
minimize the sum of squared differences of the FOD SH coefficients over all voxels.
This metric was chosen based on the generalization of Parseval's theorem to spherical
harmonics, which states that the sum of the squared SH coefficients is equal to the
power spectrum of the corresponding amplitudes. In other words, the chosen metric is
proportional to the squared amplitude difference between the FODs, averaged over all
orientations. It is therefore appropriate for this particular problem, and provides
considerable computational advantages as it avoids the need to explicitly compute the
FOD amplitudes.
   In order to eliminate gaps in the normalized image, the registration optimizes a
transformation that maps a point at each voxel in the template image to a point in the
moving image using tri-linear interpolation of the SH coefficients. Note that since the
mappings between the FOD amplitude and their respective SH coefficients is a linear
operation, performing the trilinear interpolation over any of these quantities will
produce equivalent results
   To decrease the computation time, registration was performed on FOD data with a
lmax of 4. This removes a large proportion of small unwanted noisy peaks, and reduces
the number of delta functions, n, required to sufficiently approximate the FOD
(during registration experiments we used n = 60). Once registration is completed, the
combined affine and B-spline transformation parameters are saved and applied to
FOD data with a larger lmax = 8 and n = 300, using cubic B-Spline interpolation.

3 Experiments and Results

3.1 Simulated Phantom

To demonstrate the transformation described in 2.1, a phantom with crossing fibers
was generated Fig.1(a). A non-rigid deformation was created that varies in the x
direction dependent on the y position. When the FODs were interpolated without
being transformed they were no longer aligned with the overall fiber direction as seen
in Fig.1(b). However when the FODs were transformed using the method above their
orientations matched those of the containing fiber bundle Fig.1(c). Note that even
though the vertical fiber has been re-orientated the horizontal fiber distributions
remain consistently aligned.

              a)                              b)                              c)
Fig. 1. The FOD orientations are corrected using a novel delta function approximation. a) A
generated phantom consisting of two fibres crossing at 90 degrees. b) Results after applying a
deformation in the y axis without correcting the orientations of each FOD. c) The phantom after
applying the same deformation as in b) with correction of the FOD orientation.

3.2 Self Registration

To assess whether the non-rigid registration is able to reach a global minimum
solution, a number of individual images were registered to themselves starting from a
perturbed deformation field. To ensure the initial deformation represented a realistic
transformation each subject was registered to a different subject and the resulting
deformation used to initialize the self registration.
   FODs were computed using MRtrix (J-D Tournier, Brain Research Institute,
Melbourne, Australia, [18] from HARDI data
acquired on a Siemens 3T scanner, with b=3000s/mm2, 60 DW directions and a
2.5mm isotropic voxel size. The self-registration was performed on 10 subjects, using
CP spacing of 10 mm and 40 iterations.
   Average errors were computed for those control points that lay within a tissue
mask for both the initial deformation and the resulting self-registration. The results
were 1.27±1.86 mm and 0.07±0.27 mm respectively. The errors for the top 1% of the
most severely deformed control points were 5.83±2.04 mm before and 0.64±0.25 mm
after self-registration.

3.3 Comparison with Existing Techniques

An experiment was performed in order to compare results of the FOD based
registration to those obtained using T1W and FA images. An identical ITK FFD was
implemented for scalar data. T1W, FOD and FA images were obtained from 9 of the
previously mentioned subjects and registered to a 10th representative template subject.
   We selected parameters for each method that were found to give optimal results.
Three CP spacing levels of 24, 16, and 8mm were selected for FOD and FA (with 60,
40, 20 and 300, 200, 100 iterations for each method respectively). Due to the higher
resolution of the T1W images (0.9mm isotropic) a finer CP spacing of 20, 10, and
5mm was selected with a multi resolution image pyramid with subsampling of 4, 2,
and 1 with 800, 400 and 200 iterations respectively. A full brain mask on the template
image was used for all methods. Both scalar registrations were performed using
normalized mutual information. T1W images were rigidly co-registered to b=0
images and all transforms (b0 template to T1W to T1W to b0 subject) were composed
into a single deformation.
   The resulting transformations from each method were applied to FOD (lmax = 8) and
FA volumes and then averaged. Average whole brain FA maps are shown in figure 2.
The FOD registration results (Fig. 2b) are slightly sharper than the FA driven (Fig.
2c) and T1W driven (Fig 2d) results, particularly in the superior right hemisphere. If
normalised FODs have similar orientations, the average FOD at any given voxel
should contain sharper and larger peaks. Figure 3 illustrates a close up region
indicated by the box in figure 2. Overall the fibre bundles are less blurred and contain
sharper FODs that are more similar to the template.
   To quantitatively assess the FOD sharpness in the average images from each
method, the average entropy of the FODs was computed within a masked region. The
mask was generated by thresholding the template FA map at 0.2. To compute the
entropy, each FOD was normalized and sampled in 1000 equally distributed
directions and binned in a histogram of size 128. The average entropy for the FOD,
FA and T1W registrations was 4.616, 4.642 and 4.675 respectively.

           a                       b                        c                       d
Fig. 2. The average FA map of 9 spatially normalized subjects. a) The template image (boxed
region is shown in fig 3). b) FOD driven registration average. c) FA driven registration average.
d) T1W driven registration average.

Fig. 3. A close up of a coronal view (as indicated by the box in figure 2a), of the average of 9
subjects registered to a single template subject. FODs are overlaid on a FA map. a) The
template image. b) FOD driven registration average. c) FA driven registration average. d) T1W
driven registration average.

4 Discussion

A new spatial normalization method for HARDI data has been described that
optimizes directly the FODs. The strength of this approach lies in its ability to utilize
fibre orientation information provided by FODs on crossing fibers to better align deep
white matter structures. We have demonstrated that this method provides better
spatial normalisation than FA and T1W images using both qualitative and quantitative
measures. The decreased measure of entropy over the template white matter mask
within the averaged images illustrates that overall the FODs are sharper. A possible
limitation is that the rigid coregistration of the T1W and b=0 images involved no
distortion correction; however these geometric distortions are limited to specific
regions of the brain and are unlikely to significantly affect the overall experimental
   Ultimately this method will be compared with state of the art tensor based
registration to highlight any advantages of using intra-voxel fiber orientation
information. However, the high b-value data that is used to estimate the FODs has
low SNR and increased Rician bias that is not suitable for tensor estimation. While
lower b-values are optimal for DT estimation, they cannot be used for accurately
estimating FODs. Therefore any direct comparison of these methods should be
approached with care.
   To our knowledge, the only registration algorithm to date that exploits crossing
fibre information available from higher order models was developed by Barmpoutis et
al. [13]. In this work, the Hellinger distance was selected as the similarity metric to
compare 4th order DTs. However, this metric is both scale and rotationally invariant
and therefore important orientation information is discarded.
   We have developed a novel method for transforming FODs using a delta function
approximation and incorporated it into an ITK based FFD registration framework to
enable fair comparisons to be made with alternative approaches. FOD transformation
requires that the deformation field is well-defined throughout the entire domain. At
the current CP spacing folding has not been a problem, however in future work we
intend to enforce this requirement using a Jacobian determinant constraint.
   An alternative method for transforming FODs was recently presented by Hong et
al [19]. Both methods assume that the volume fraction of each fiber population should
be preserved and theoretically similar results should be obtained. However, the
advantage of this method is the lack of a computationally expensive pseudoinverse
required to estimate the transformed FOD SH coefficients (step 5 in [19]). Future
work will focus on a comparison of these two methods.

5 Conclusion

We have proposed a new method for non-rigid registration and spatial normalization
of FODs. Our results demonstrate the advantages of exploiting intra-voxel fiber
orientation information compared to T1W and FA based scalar methods. This work
has direct implications for techniques such as atlas generation, tensor based
morphometry analysis, and group average fiber tracking. Improvements to these
methods will ultimately benefit group based clinical and neuroscientific studies.


1.    Tournier, J., Calamante, F., Gadian, D.G., Connelly, A.: Direct Estimation of the
      fiber Orientation Density Function from Diffusion-Weighted MRI Data Using
      Spherical Deconvolution. Neuroimage 23(3) 1176-1185. (2004)
2.    Descoteaux, M., Deriche, R., Knosche, T.R., Anwander, A.: Deterministic and
      Probabilistic Tractography Based on Complex Fibre Orientation Distributions.
      IEEE Trans Med Imaging 28(2) 269-86. (2009)
3.    Tuch, D.S., Reese, T.G., Wiegell, M.R., Makris, N., Belliveau, J.W., Wedeen,
      V.J.: High Angular Resolution Diffusion Imaging Reveals Intravoxel White
      Matter Fiber Heterogeneity. Magn. Reson. Med. 48(4) 577-82. (2002)
4.    Jones, D.K., Griffin, L.D., Alexander, D.C., Catani, M., Horsfield, M.A.,
      Howard, R., Williams, S.C.R.: Spatial Normalization and Averaging of Diffusion
      Tensor MRI Data Sets. NeuroImage 17(2) 592-617. (2002)
5.    Guimond, A., Guttmann, C., Warfield, S., Westin, C.: Deformable Registration
      of DT-MRI Data Based on Transformation Invariant Tensor Characteristics.
      Proceedings of the 1st Annual Meeting of the ISBI, Washington DC, USA 761-
      764. (2002)
6.    Basser, P.J., Mattiello, J., LeBihan, D.: MR Diffusion Tensor Spectroscopy and
      Imaging. Biophys J 66(1) 259-67. (1994)
7.    Park, H., Kubicki, M., Shenton, M.E., Guimond, A., McCarley, R.W., Maier,
      S.E., Kikinis, R., Jolesz, F.A., Westin, C.: Spatial Normalization of Diffusion
      Tensor MRI Using Multiple Channels. NeuroImage 20(4) 1995-2009. (2003)
8.    Zhang, H., Yushkevich, P.A., Alexander, D.C., Gee, J.C.: Deformable
      Registration of Diffusion Tensor MR Images with Explicit Orientation
      Optimization. Med. Image. Anal. 10(5) 764-85. (2006)
9.    Alexander, D.C., Pierpaoli, C., Basser, P.J., Gee, J.C.: Spatial Transformations of
      Diffusion Tensor Magnetic Resonance Images. IEEE Trans. Med. Imaging
      20(11) 1131-9. (2001)
10.   Zhang, H., Avants, B.B., Yushkevich, P.A., Woo, J.H., Wang, S., McCluskey,
      L.F., Elman, L.B., Melhem, E.R., Gee, J.C.: High-Dimensional Spatial
      Normalization of Diffusion Tensor Images Improves the Detection of White
      Matter Differences: An Example Study Using Amyotrophic Lateral Sclerosis.
      IEEE Trans. Med. Imaging 26(11) 1585-97. (2007)
11.   Behrens, T.E.J., Berg, H.J., Jbabdi, S., Rushworth, M.F.S., Woolrich, M.W.:
      Probabilistic Diffusion Tractography with Multiple Fibre Orientations: What Can
      We Gain? Neuroimage 34(1) 144-55. (2007)
12.   Ozarslan, E., Mareci, T.H.: Generalized Diffusion Tensor Imaging and
      Analytical Relationships Between Diffusion Tensor Imaging and High Angular
      Resolution Diffusion Imaging. Magn Reson Med 50(5) 955-965. (2003)
13.   Barmpoutis, A., Vemuri, B.C., Forder, J.R.: Registration of High Angular
      Resolution Diffusion MRI Images Using 4th Order Tensors. Med Image Comput
      Comput Assist Interv Int Conf Med Image Comput Comput Assist Interv 10(Pt
      1) 908-915. (2007)
14. Zhan, W., Yang, Y.: How Accurately Can the Diffusion Profiles Indicate
    Multiple Fiber Orientations? A Study on General Fiber Crossings in Diffusion
    MRI. J. Magn. Reson 183(2) 193-202. (2006)
15. Tournier, J., Yeh, C., Calamante, F., Cho, K., Connelly, A., Lin, C.: Resolving
    Crossing Fibres Using Constrained Spherical Deconvolution: Validation Using
    Diffusion-Weighted Imaging Phantom Data. Neuroimage 42(2) 617-25. (2008)
16. Rueckert, D., Sonoda, L., Hayes, C., Hill, D., Leach, M., Hawkes, D.: Nonrigid
    Registration Using Free-Form Deformations: Application to Breast MR Images.
    IEEE Trans. Med. Imaging 18(8) 712-721. (1999)
17. Ibanez, L., Schroeder, W., Ng, L., Cates, J.: The ITK Software Guide. (2005)
18. Tournier, J., Calamante, F., Connelly, A.: Robust Determination of the fibre
    Orientation Distribution in Diffusion MRI: Non-Negativity Constrained Super-
    Resolved Spherical Deconvolution. Neuroimage 35(4) 1459-1472. (2007)
19. Hong, X., Arlinghaus, L.R., Anderson, A.W.: Spatial Normalization of the Fiber
    Orientation Distribution Based on High Angular Resolution Diffusion Imaging
    Data. Magn Reson Med (2009)

Shared By: