Non-linear integration of DTI-based fiber tracts into standard 3D MR data

(1)

Non-linear Integration of DTI-based Fiber Tracts into Standard 3D MR Data

Dorit Merhof^1,2, Peter Hastreiter^1,2, Grzegorz Soza^1,2, Marc Stamminger², Christopher Nimsky¹

1Neurocenter, Dept. of Neurosurgery Schwabachanlage 6, 91054 Erlangen, Germany

2Computer Graphics Group, University of Erlangen-Nuremberg Am Weichselgarten 9, 91058 Erlangen, Germany

Abstract

Diffusion tensor imaging (DTI) provides information about the location of white matter tracts within the human brain. This information is essential for preoperative neurosurgical planning to achieve maximal tumor resection while avoiding postoperative neurological deficits. Due to the anatomical distortion of echo planar imaging, DT images - and as a result the fiber tracts computed from them - are distorted. In this paper, we present a novel approach to account for those distortions. All voxels containing fibers within the distorted DT dataset were marked. Subsequently, a non-linear registration with standard 3D MR data was performed. The marked voxels were re-extracted from the registered DT dataset and displayed within the 3D MR dataset.

The strategy introduced in this paper is an essential prerequisite for the integration of fiber tract data into 3D MR datasets. The fused data is of high value for neuronavigation and thereby a benefit for neurosurgery.

1 Introduction

Diffusion tensor imaging (DTI) is a non-invasive in vivo imaging modality and is used for the computation of fiber-tracts within tissues such as the human brain. This kind of imaging is based on the phenomenon of water diffusion known as Brown- ian motion. Neural fibers consist of long cylindri- cal cells that are filled with fluid and are bounded by a cell membrane. The diffusion behavior within neural cells is anisotropic, as the less water perme- able cell membranes restrict diffusion. Therefore, both, tissue structure and architecture at the micro-

scopic level can be reflected by DTI. Those structures can be revealed by computing tensors from a series of directional diffusion images. Techniques commonly referred to as tractography are then applied upon the vector field derived by determining the dominant eigenvector of each tensor. The tracking procedure between a source and target region results in curvilinear paths which most likely represent the course of fibers within the brain. Fiber tracking thus reveals the individual location of im- portant fiber tract systems known from anatomical studies.

A more detailed introduction to the basics of DTI is available in [1, 2, 3], a number of different possibilities to visualize DTI data is presented in [4, 5, 6]. The tractography algorithm used in this work is based on the strategies published by Basser and Mori [7, 8, 9, 10].

For neurosurgery these techniques are of particu- lar importance as fiber tracts are visually hardly dis- tinguishable from surrounding tissue. The course of fibers changes considerably between different indi- viduals and is especially in brain tumor cases un- predictable. DTI and fiber tracking are therefore essential tools for neurosurgery [11].

However, DTI does not show anatomical structures in detail due to limitations regarding the reso- lution of the data. In addition, DTI suffers from the fact that distortions occur at the proximity of the skull base and near air-filled spaces such as in the case of the brainstem and the frontal lobe (Fig. 1).

Those artifacts can be attributed to the small bandwidth of the echo planar imaging sequence in phase encoding direction and the inhomogeneity of the magnetic field due to a varying susceptibility within the measured tissue.

3D anatomical data, such as MPRAGE (Magne- First publ. in: Vision, Modeling and Visualization 2004. Proceedings / Bernd Girod,

Marcus Magnor, Hans-Peter Seidel (eds.). Heidelberg: AKA, 2004, pp. 371-378

Konstanzer Online-Publikations-System (KOPS) URN: http://nbn-resolving.de/urn:nbn:de:bsz:352-opus-92182

URL: http://kops.ub.uni-konstanz.de/volltexte/2009/9218/

(2)

Figure 1: Slice of distorted DT image (B0 image)(left)and corresponding MPRAGE slice(right)

tization Prepared Rapid Acquisition Gradient Echo) images on the other hand are able to represent even tiny anatomical structures in detail. The distortion of MPRAGE images is negligible so that they can be used for neuronavigation, i.e. in a stereotactic system.

Integration of fiber tracts into a MPRAGE dataset would allow visualizing the spatial relation of a tumor and adjacent fiber tracts for neurosurgical planning, thus allowing to avoid postoperative neurological deficits. There are several strategies to solve this problem: By measuring field maps describ- ing image distortions applying other MR sequence techniques or performing a non-linear registration of DTI and MPRAGE data so that distortions can be taken into account.

In this paper, we present an approach to realize this task via non-linear registration. First the fibers are computed from the DT images. After that all voxels containing fibers within the diffusion tensor dataset are marked. In section 2 and 3 these steps are outlined in detail and some information on image acquisition is given. Finally, a piecewise linear registration between diffusion tensor and MPRAGE dataset is performed which is presented in section 4.

The results derived from our approach are analyzed in section 5.

2 Image Data

All images were acquired at the Department of Neu- rosurgery, Erlangen, Germany using a Siemens MR Magnetom Sonata Maestro Class 1.5 Tesla scan-

ner equipped with a gradient system with a field strength of up to 40 mT/m (effective 69 mT/m) and a slew rate of up to 200 T/m/s (effective 346 T/m/s).

DT image data: For the computation of tensors and fibers a series of directional diffusion images is needed. The diffusion images differ in the direction of the gradient that is applied. For each slice a total amount of seven images was created: One image was measured without any gradient direction which is equivalent to common MR data. In the following, this image will be referred to as B0 image. The remaining six images were acquired using different gradient directions that correspond to the diffusion directions measured at that time. Those images will be referred to as B1, ... , B6 image. Each of the B0 - B6 images shows the distortion artifacts men- tioned above. As regards anatomical information, the B0 image represents anatomical structures best though in poor detail. The other images primarily show diffusion properties.

DT imaging parameters:TR = 9200, TE = 86 ms, bhigh = 1000 s/mm², blow = 0 s/mm², field of view 240 mm, voxel size 1.875×1.875×1.9 mm³, 1502 Hz/Px bandwidth, acquisition matrix 128×128, 60 slices with no intersection gap were measured (measurement time for 5 averages: 5 min and 31 sec), the diffusion-encoding gradients for the six diffusion weighted images were directed along the following axes: (±1,1,0), (±1,0,1) and (1,±1,0).

(3)

Figure 2: The left pyramidal tract in a healthy volunteer(left), respective voxels containing the pyramidal tract(middle), 3D visualization of marked voxels within the B0 volume(right).

MPRAGE image data: The advantages of con- ventional 3D imaging are a high signal to noise ra- tio, thin coherent slices and the possibility of multi- planar reconstruction.

MPRAGE imaging parameters: TR = 2020 ms, TE = 4.38 ms, field of view 250 mm, voxel size 0.488281×0.488281×1.0 mm³, acquisition matrix 256×256 interpolated (512×512), 160 slices, the time required for acquisition of the MPRAGE data was 8 min 39 sec.

3 Image Processing

In a first step, the computation of fibers from the DT data was performed using the DTI task card version 1.66 (Magnetic Resonance Center, Mas- sachusetts General Hospital, Boston) on a Siemens scanner using the MR software MRease N4 VA21B under syngo VB10I. A knowledge-based multiple- ROI (region of interest) approach was applied, initi- ating the tracking algorithm from at least one user- defined seed region. Tracking was performed in both retrograde and orthograde direction according to the direction of the principal eigenvector within the ROI. A vector step length of 0.5 mm, an an- gular threshold of 35^o and a fractional anisotropy threshold of 0.3 were used in this study. The computed tract coordinates were saved to a file in raw data format. This fiber information is later used for marking voxels within the B0 volume that contain fibers. Fig. 2 (left) shows fiber bundles computed from the diffusion datasets.

For further processing, the original DICOM data was transformed to the.hdrformat using in house software and FSL (FMRIB Software Library) [12].

The next step comprises the segmentation of the brain which in this case was automatically accom- plished using BET (Brain Extraction Tool) [13].

The fiber information derived from the DTI task cardis then used to determine all voxels containing fibers within the B0 dataset. Those voxels are marked by setting the respective grey value to a pre- defined, unambiguous value. The marked voxels derived from the fiber data are shown in Fig. 2 in the middle. On the right, the complete volume containing marked voxels is presented. The marked B0 dataset together with the MPRAGE imaging data are later used for registration.

4 Registration

In a first step, the B0 dataset was rigidly registered to the MPRAGE dataset. The result of the rigid registration was taken as an initial estimation for a nonrigid transformation which was then performed using Free-Form Deformation (FFD). To accomplish the deformation hardware-based three-dimensional B´ezier functions were considered since B´ezier functions provide a mechanism for their modification and are characterized by intuitive behavior on their change [14, 15, 16].

In our approach the source B0 dataset was associated with the object space OS. The function P : P S 7→ OS leads from the parameter spaceP S defined as[0,1]³ into this object space

(4)

Figure 3: Distorted B0 image before(left)and after(middle)registration – MPRAGE image(right)

and served as parameterization of the object space.

The second dataset, i.e. the MPRAGE dataset, re- mained fix.

The deformation functionD:P S 7→ T Slead- ing from the parameter space into the texture space T S is represented by a B´ezier function. Thus, the shape of this deformation function is uniquely defined by the corresponding lattice of control points bi,j,k(i = 0, . . . , l, j = 0, . . . , m, k = 0, . . . , n) placed within texture space. This deformation function is then expressed by means of a trivariate tensor product whereas the basis func- tionsBi^l, Bj^m, Bⁿk are Bernstein polynomials of or- derl, mandn, respectively:

D(s, t, u) = Xl

i=0

Xm

j=0

Xn

k=0

Bi^l(s)Bj^m(t)Bkⁿ(u)bi,j,k . (1) In order to accelerate the computation of the FFD, the benefits of graphics hardware were exploited.

All trilinear interpolations needed for computing the similarity measure were performed in the texture processing unit of the graphics hardware. For this purpose, the image volume within the object space was recomputed and loaded into the 3D texture memory of the graphics card. In order to per- form texture mapping, the texture spaceT S being [0,1]³is associated with the texture memory. The advantage of this approach is that less computation time is needed [14] since there is no need to recalcu- late the whole 3D dataset voxel by voxel in software in order to obtain the new intensities.

Movements of the control pointsbi,j,kin the lattice are followed by immediate changes of the de-

formation functionDand thus of the deformed B0 volume. The main idea of the non-rigid registration is to manipulate free control points in the lattice in such a way that the similarity measure between the B0 volume which is deformed with FFD and the reference MPRAGE volume takes its maximum. Thus, the registration task is equivalent to a multidimensional optimization problem. The quality measure which was used to assess the quality of the alignment is based on normalized mutual information (NMI).

5 Results and Discussion

The processing steps described in section 3 were ex- ecuted for the left pyramidal tract of five healthy volunteers. All computations were performed on a PC equipped with an Intel Pentium 4 (2.4 GHz) and a NVIDIA Quadro4 graphics card with 128 MB graphics memory. The average computation time was 1.5 h. Thereby, 30 min were required for data conversion and brain extraction whereas the major part of the overall computation time was consumed by the registration.

The average number of marked voxels within the B0 dataset was 1490 and the total number of voxels amounted to 983040 (128×128×60 voxels). This means that only a percentage of about 0.15% voxels within the whole volume was marked. For this reason it can be assumed that the influence of the marked voxels on the registration process is negligible.

The registration of the marked B0 dataset with the MPRAGE dataset was used to compensate

(5)

Figure 4: Marked voxels within a B0 image before(left)and after(middle)registration and analysis (see arrow) of the maximal displacement(right)

for the distortion artifacts (up to 28 mm in the frontal lobes) of the B0 dataset. In Fig. 3 the benefits of our registration approach can be observed: The figure shows corresponding sagit- tal slices of the B0, the registered B0 and the MPRAGE dataset. The registered B0 dataset in the middle has compensated for the distortions that were visible before registration. Visual inspection of the registration result (using a semitransparent overlay) by the senior author having more than 14 years experience with MRI revealed a satisfy- ing correspondence in the gyri of the central region, the internal capsule and the brainstem.

For evaluation purposes corresponding axial slices within the B0 and the registered B0 dataset were analyzed. For every second fiber-containing slice we determined the maximum displacement of the fibers, i.e. the marked voxels. We measured the maximum distance between two corresponding boundary points. In Fig. 4 this approach is illus- trated.

Fig. 5 presents the results of our measurements with respect to the registered B0 volume: In the axial slices 0-24 and 124-160 the pyramidal tract was not visualized. So, the analysis started in slice 25 and ended in slice 123. In between, the displacement of the fibers, i.e. marked voxels, was measured for every second axial slice. The dia- gram on the left shows our measurements denoted in mm. The results are consistent with our expecta- tions showing a slight distortion at the brain surface (in the central region) and distinct distortions in the brainstem due to nearby air-filled structures (clivus

and sphenoid sinus). Accordingly, the distortions in the region of the internal capsule are only minor. In general, we observed only slight left-right displace- ments whereas the main displacement occurred in anterior-posterior direction, i.e. in phase-encoding direction.

Fig. 6 illustrates the distortion by blending an unregistered and a registered B0 dataset. Finally, the marked voxels were re-extracted from the registered B0 dataset and displayed within the MPRAGE dataset (Fig. 7).

An alternative approach would be to register B0 and MPRAGE volumes and apply the resulting transformation to the B1-B6 datasets. Computing the fiber tracts directly from the distortion corrected B0-B6 datasets provides similar results. However, the presented strategy clearly shows the volume of the fibers which is of importance for surgery.

6 Conclusion and Future Work

In this paper we presented a novel approach to account for the inherent distortions of DTI data and the fiber tracts computed from them. We observed slight distortions at the cortical surface and distinct distortions of up to 9 mm in the brainstem. The techniques developed in this context can be applied for DTI datasets in patients with space occupying lesions. This enables spatially correct preoperative neurosurgical planning and thus allows for the integration of fiber tracts in datasets used for neuronavigation.

(6)

0 2 4 6 8 10

020406080100120140160

distortion in mm

slice number

Figure 5: Distortion of left pyramidal tract in corresponding axial slices in a B0 and the registered B0 dataset with respect to the registered B0 volume. For each second axial slice the movement of the fibers was measured showing that most distinct distortions occurred near the brainstem

Future work will focus on the acceleration of the preprocessing steps, especially the registration approach, to achieve intra-operative applicability.

A clinical study is currently in progress to verify the value of the presented approach in pathological cases.

7 Acknowledgments

This work was supported by the Deutsche Forschungsgemeinschaft in the context of SFB 603, Project C9 and the Graduate Research Center “3D Image Analysis and Synthesis”. We are especially grateful to A. G. Sorensen (Department of Radi- ology/Nuclear Magnetic Resonance Center, Mas- sachusetts General Hospital, Boston, USA) for pro- viding theDTI task card softwarewe used for the computation of the fiber tracts.

References

[1] Westin, C., Maier, S.E., Mamata, H., Nabavi, A., Jolesz, F., Kikinis, R.: Processing and visualization for diffusion tensor MRI. Med Im- age Anal6(2002) 93–108

[2] Basser, P., Pierpaoli, C.: Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor MRI. J Magn Reson Imaging111(1996) 209–219

[3] Jones, D.K., Horsfield, M.A., Simmons, A.:

Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magn Reson Med42 (1999) 515–

525

[4] Kindlmann, G.: Superquadric tensor glyphs.

In: Proc. IEEE Symposium on Visualization.

(2004) (accepted)

[5] Zhukov, L., Barr, A.: Oriented Tensor Re- construction: Tracing Neural Pathways from Diffusion Tensor MRI. In: Proc. IEEE Visu- alization. (2002)

[6] Zhang, S., Kindlmann, G., Laidlaw, D.H.:

Diffusion tensor mri visualization. In: Visu- alization Handbook. (2004) inpress.

[7] Basser, P., Pajevic, S., Pierpaoli, C., Duda, J., Aldroubi, A.: In Vivo Fiber Tractography Using DT-MRI Data. Magn Reson Med 44 (2000) 625–632

[8] Mori, S., Crain, B., Chacko, V., van Zijl, P.:

Three-dimensional tracking of axonal projec- tions in the brain by magnetic resonance imaging. Ann Neurol45(1999) 265–269

(7)

[9] Mori, S., van Zijl, P.: Fiber tracking: princi- ples and strategies – a technical review. NMR Biomed15(2002) 468–480

[10] Stieltjes, B., Kaufmann, W., van Zijl, P., Fred- ericksen, K., Pearlson, G., Solaiyappan, M., Mori, S.: Diffusion tensor imaging and axonal tracking in the human brainstem. NeuroImage 14(2001) 723–735

[11] Nimsky, C., Ganslandt, O., Hastreiter, P., Wang, R., Benner, T., Sorensen, A., Fahlbusch, R.: Intraoperative diffusion tensor imaging: shifting of white matter tracts during neurosurgical procedures - initial experience.

Radiology [in press] (2004) in press.

[12] Image Analysis Group, Functional Magnetic Resonance Imaging of the Brain (FMRIB), Oxford, UK: FMRIB Software Library Release 3.1, University of Oxford. (2003) http://www.fmrib.ox.ac.uk/fsl/.

[13] Smith, S.M.: Fast robust automated brain extraction. Human Brain Mapping 17 (2002) 143–155 http://www.fmrib.ox.ac.uk/fsl/bet/.

[14] Soza, G., Bauer, M., Hastreiter, P., Nimsky, C., Greiner, G.: Non-rigid Registration with Use of Hardware-Based 3D B´ezier Functions.

In: Proc. MICCAI. Lect. Notes in Comp. Sc., Springer (2002) II,549–556

[15] Hastreiter, P., Rezk-Salama, C., Soza, G., Greiner, G., Fahlbusch, R., Ganslandt, O., Nimsky, C.: Strategies for Brain Shift Eval- uation. Med Image Anal (2004) accepted for publication.

[16] Rueckert, D., Sonoda, L.I., Hayes, C., Hill, D., Leach, M., Hawkes, D.: Nonrigid Reg- istration Using Free-Form Deformations: Ap- plication to Breast MR Images. IEEE Trans Med Img18(1999) 712–721