Document Sample

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 1 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, 4 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 brain. 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 3 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). n 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) S Where: s is the surface of the sphere, f is a delta function along orientation i, and m 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 ( ) −1 (4) 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 [17]. 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. 5 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, http://www.brain.org.au/software/) [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. 7 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 results. 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. 9 References 1. Tournier, J., Calamante, F., Gadian, D.G., Connelly, A.: Direct Estimation of the ﬁber 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 ﬁbre 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)

DOCUMENT INFO

Shared By:

Categories:

Tags:

Stats:

views: | 15 |

posted: | 9/30/2011 |

language: | English |

pages: | 10 |

OTHER DOCS BY wuyunqing

How are you planning on using Docstoc?
BUSINESS
PERSONAL

Feel free to Contact Us with any questions you might have.