Approximating Anatomical Brain Connectivity with Diffusion Tensor ...
Approximating Anatomical Brain Connectivity with Diffusion Tensor MRI Using Kernel-Based Diffusion Simulations Jun Zhang1 , Ning Kang1 , and Stephen E. Rose2 1
Laboratory for High Performance Scientific Computing and Computer Simulation, Department of Computer Science, University of Kentucky, Lexington, KY 40506-0046, USA {jzhang, nkang2}@cs.uky.edu http://www.cs.uky.edu/∼jzhang 2 Centre for Magnetic Resonance, University of Queensland, Brisbane, QLD, 4072 Australia [email protected]
Abstract. We present a new technique for noninvasively tracing brain white matter fiber tracts using diffusion tensor magnetic resonance imaging (DT-MRI). This technique is based on performing diffusion simulations over a series of overlapping three dimensional diffusion kernels that cover only a small portion of the human brain volume and are geometrically centered upon selected starting voxels where a seed is placed. Synthetic and real DT-MRI data are employed to demonstrate the tracking scheme. It is shown that the synthetic tracts can be accurately replicated, while several major white matter fiber pathways in the human brain can be reproduced noninvasively as well. The primary advantages of the algorithm lie in the handling of fiber branching and crossing and its seamless adaptation to the platform established by new imaging techniques, such as high angular, q-space, or generalized diffusion tensor imaging.
1
Introduction
A number of fiber tracking algorithms have been developed since the appearance of diffusion tensor magnetic resonance imaging (DT-MRI). Typical fiber tracking schemes, including the streamline-based technique [1, 2, 10], reconstruct the white matter tracts by tracing down in a voxel-by-voxel manner, using an estimate of the local fiber orientation determined by the principal eigenvector in each voxel that is assumed to align with the mean fiber direction in that voxel. These techniques appear to give excellent results in many instances if the principal eigenvector field is smooth. However, it suffers from several significant limitations. The vector field is error prone in the sense that the noise of DT-MRI data will influence the direction of the principal eigenvector, yielding an accumulation of orientational errors and thus an erroneous fork of the trajectory. Another major restriction is that these techniques may also be affected by partial volume G.E. Christensen and M. Sonka (Eds.): IPMI 2005, LNCS 3565, pp. 64–75, 2005. c Springer-Verlag Berlin Heidelberg 2005 °
Approximating Anatomical Brain Connectivity
65
effects [19], leading to unstable tracking through the primary eigenvector field in regions of fiber crossing, branching, or merging. Under the diffusion tensor imaging platform, a variety of methods have been proposed aiming to palliate the difficulties with more information incorporated from the diffusion tensor data. The algorithm presented in [8] uses a deflection term obtained from the diffusion tensor to improve the image noise immunity. Other schemes [7, 15] use predefined knowledge to group together neighboring voxels based on a similarity measure. Taking into account the uncertainty of fiber direction, probabilistic and statistical approaches [3, 5, 12] have been developed to mitigate the effects of fiber crossing and diverging as well as the sensitivity to noise. The level set theory is also utilized to find fiber paths connecting different brain regions [13]. Another front evolution algorithm proposed in [17] utilizes the fiber orientation function to reconstruct fiber tracts. As the measured quantity in DT-MRI is for water diffusion, an intuitive way to gain insights from the diffusion tensor data is to treat the brain volume as a physical system and simulate a virtual water diffusion process over it, which is anisotropic and governed by the diffusion equation. The shape of the anisotropic diffusion, represented by diffusion fronts, can be used to estimate the directional arrangement of the underlying white matter fiber bundles. This reflection is based upon the principle that the faster the diffusion, the longer the distance will be traveled on average by water molecules within the same amount of diffusion time. The fiber tracts are thus expected to proceed along the direction where the diffusion is the greatest. The fiber tractography presented in this paper performs simulations of the diffusion process stemming from a series of diffusion starting voxels, within corresponding overlapped 3D diffusion kernels. The diffusion simulation initiated from a diffusion root node is utilized to construct a diffusion front in its associated kernel. The next set of diffusion root nodes, where a seed will be placed, are located on the diffusion front which is generated by the diffusion process initiated from a previous seeded root voxel. They are picked up according to the created distance map and the local orientation information involving these voxels and the diffusion root node. For the next round, each of the newly selected diffusion root voxels will be used to generate a front by starting a diffusion process in its own kernel. Given below is the detailed description of the diffusion-based fiber tractography theory and algorithm.
2 2.1
Methods The Anisotropic Diffusion Equation
The anisotropic diffusion process simulated in this work is governed by the equation ∂C = ∇ · (D∇C), (1) ∂t where D is the so-called diffusion coefficient, which is a second-order tensor in the presence of anisotropy, C the concentration, and t the independent time variable. The coefficient used in the anisotropic diffusion simulation is the diffusion tensor
66
J. Zhang, N. Kang, and S.E. Rose
calculated from the diffusion-weighted imaging data, which is represented by a three-by-three symmetric positive definite matrix, Dxx Dxy Dxz D = Dyx Dyy Dyz , Dzx Dzy Dzz where the subscripts xx, xy, xz, etc., denote the values of the individual coefficients in the matrix. 2.2
Extracting Front in Diffusion Kernel
The first step to reconstruct fiber pathways starting from a pre-chosen root node s involves simulating the diffusion process in its associated diffusion kernel, initiated from a seed (an initial concentration value) in this voxel. A diffusion kernel defined here is a cube with six rectangular sides, which covers only a small bulk of the whole 3D data volume and is geometrically centered upon the diffusion starting voxel. The virtual concentration seed of water spreads from the root node through the neighboring nodes, within a limited amount of time, forming a diffusion front which is the surface of a diffusion volume containing nodes with nonzero concentration values. The expansion of the diffusion volume originated from the root node is achieved by integrating Eq. (1) over a certain amount of time, subject to the following initial condition, ½ 1 at the root node, = (2) C 0 elsewhere in the diffusion kernel. t=0 The boundary of the diffusion kernel is assumed to be insulated, i.e., (D∇C) · n = 0,
(3)
where n is the direction normal to the boundary. This condition implies that the normal part of the gradient of the concentration on the boundary is zero, in other words, nothing escapes out of the domain. We have developed an unsteady state anisotropic diffusion solver framework, which is adapted to the cerebral circumstance and runs in both sequential and parallel computing environments. In the current paper, (1) was solved sequentially under the initial condition (2) and boundary condition (3) by resorting to the established computational framework. We used a diffusion kernel with dimensions 11 × 11 × 7 and a voxel size same as in the original data volume, which proved to be very efficient in time cost and showed no impairment on tracking performance. Fig. 1 shows a concentration distribution map of the anisotropic diffusion simulated in a diffusion kernel using the human brain tensor data. Once the time integration for solving (1) is done, a discrete approximation to the diffusion front can be calculated in terms of whether the concentration value is zero in a voxel. Thus all nodes in the diffusion kernel can be partitioned into two groups, one with zero concentration and the other with nonzero values. Since only one seed is diffused over the root node, the diffusion-swept volume, denoted
Approximating Anatomical Brain Connectivity
67
Fig. 1. The concentration distribution map of the anisotropic diffusion simulated in a diffusion kernel (bounded by the white rectangle). The profile is superimposed on a grey-scale axial map of the fractional anisotropy
as V (r), is comprised of voxels with nonzero concentration values, where r is the position of the root node. For each member of V (r), we consider its surrounding 26 closest neighboring nodes in a 3 × 3 × 3 cube. Let i, j, k index the relative coordinates of the 26 nearest neighbors of r with i, j, k ∈ {−1, 0, 1}. If F (r) is the set of voxels that form the diffusion front of r, then for any node p ≡ (px , py , pz ) ∈ V (r), we define p ∈ F (r) if ∃ (i, j, k), such that (px − i, py − j, pz − k) ∈ / V (r), which implies that if any of the 26 nodes is not in V (r), then p ∈ F (r). 2.3
The Criteria Set
In order to store and handle the front nodes dynamically produced in each diffusion kernel, we set up a queue Q, a first-in first-out data structure. Q is initialized to contain just the starting node s, i.e., Q = {s}, thereafter, Q always contains the set of diffusion front nodes which will be the subsequent diffusion root voxels. Once F (r) is computed for the root node r, we further apply the criteria in the set C (see below) to the nodes of F (r) and pick up those that meet the corresponding thresholds. We define I(r) to be the set of nodes selected from F (r) that satisfy the criteria in C, i.e., I(r) = { p ∈ F (r) | p meets all the criteria in C }. The qualified nodes in F (r) are inserted into I(r) in a non-ascending order of α (see the first criteria c1 below). I(r) is then appended to the tail of the queue Q. The set C bears a number of criteria, which determine the connection of fiber pathways. There are five criteria in C used to evaluate the information about distance S and orientation between the root r and its front nodes in F (r). Let C = i∈{1,··· ,5} {ci }. The first criterion, c1 , is the threshold for distance ratio measure α, which is defined as α = d/dmax , where d = kv(r)k is the Euclidean distance in R3 of two points, connected by the vector v(r) pointing from r to a node in F (r), while dmax is the maximum value among the d’s. We set c2 to be a threshold of an invariant anisotropy index, the fractional anisotropy (f a) [14]. The next criterion, c3 , is a curvature constraint introduced to secure the tracks yielded moving forward consistently and smoothly without erratically turning
68
J. Zhang, N. Kang, and S.E. Rose
back on themselves. A threshold is used to restrain the angle between v(π(r)) and the current direction of tracking, v(r). Here, π(r) is the predecessor voxel of r, i.e., r ∈ I(π(r)). v(π(r)) is an established vector pointing from π(r) toward r, which implies the presence of a trajectory passing in this direction. c4 is used to judge the coherence of fiber directions along the reconstructed trajectories passing through r. One threshold is set on three inner products, φ1 , φ2 , and φ3 , where φ1 = |ˆ v (r) · e1 (r)|, φ2 = |ˆ v (r) · e1 (f )|, and φ3 = |e1 (r) · e1 (f )|. Here, ˆ (r) = v(r)/kv(r)k; e1 (r) and e1 (f ) are principal eigenvectors (corresponding v to the largest eigenvalue of D) at the voxel r and f ∈ F (r), respectively. The last criterion, c5 , specifies the maximum number of voxels I(r) allowed to have if there are more voxels than expected satisfying all previous four criteria, which controls the overall computational time for simulating the diffusion process in diffusion kernels. 2.4
Recovering Fiber Pathways
When Q is not empty, the current head node of Q is removed off the queue and is considered to be a new root r′ where a seed is diffused. r′ is positioned at the geometrical center of the diffusion kernel, which is then initialized using the global-to-local mapping to retrieve necessary information from the original data volume for carrying out the new diffusion simulation. The diffusion front F (r′ ) is calculated in the same way as that of F (r). As in the derivation of I(r), the set I(r′ ) is determined as well by checking each member of F (r′ ) based upon the criteria in C, then it is added to the tail of Q. We continue in this way by repeatedly taking off the head node of Q and processing it as a new root to diffuse a seed over it, until the queue becomes empty. Aimed at recovering the fiber pathways after constructing the diffusion front in each kernel, each voxel p in the global data grid owns a memory of its predecessor voxel, π(p), where p ∈ I(π(p)). Since every voxel in the grid can be taken as a diffusion root no more than once, π(p) is the sole predecessor of p if there is one. Thus, back propagation from the voxels on diffusion fronts by following continuously the corresponding predecessor voxels leads to paths that merge to the starting voxel s. This merging corresponds to the procedure that can be viewed in the reverse direction as fiber tracts branch outwards from s. Finally, the pathways are smoothed out using B-spline least-square approximations. The diffusion equation-based fiber tractography procedure is outlined in Alg. 1. 2.5
Connectivity Index
Since the size of the set I(r) can be larger than one for any root r, there may exist a bunch of reconstructed fiber pathways that branch outwards from the voxel where tracking starts. We utilize a heuristic connectivity index, ξ, as a confidence measure to estimate the odds that any generated path well approximates a true anatomical connection. For a given putative pathway, the index is defined as Qn i ξ = i=1 ( αi +β 2 ), where α is the distance ratio, β = (φ1 + φ2 + φ3 )/3, and n the number of diffusion kernels used to produce the pathway. The definition of ξ may
Approximating Anatomical Brain Connectivity
69
Algorithm 1. Fiber tracking using kernel-based diffusion simulations 1: specify a starting node s and initialize Q such that Q = {s} 2: while Q is not empty do 3: remove the head node r off Q and take it as a root 4: initialize the 3D diffusion kernel that is geometrically centered on r 5: get V (r) by solving the diffusion equation (1) over the diffusion kernel with the initial and boundary conditions (2) and (3) imposed 6: compute F (r), then determine I(r) and append it to the tail of Q 7: end while 8: record π values for voxels during front construction 9: retrieve fiber pathways using back propagation
be construed as evaluating how faithful the computed pathways are to follow the fastest diffusion direction, yet adjusted by coherence with local fiber orientations.
3 3.1
Data Acquisition Synthetic Tensor Fields
In order to assess fidelity and robustness of the tracking algorithm, we generated synthetic DT-MRI data with a uniform voxel size of 1 mm3 , where the true path of a fiber tract is known. The tensor field was constituted upon an anisotropic and an isotropic tensors taken out of real DT-MRI data. The shape of the diffusion tensor in synthetic fibers was described by the anisotropic one such that λ1 : λ2 : λ3 was approximately 2.5 : 1 : 1, while the isotropic one was used to forge the background of the simulated tensor field. The vector field for fiber orientations was derived by sampling discretely the trajectories which were analytically defined. To make the simulated field more realistic, an approximation to Rician noise [4] was added in the diffusion-weighted images which were calculated from the Stejskal-Tanner equation using the gradient sequence in [18] and a b-value of 1000. The noisy realization led to a signal-to-noise ratio of 10. A compact analytic solution to the Stejskal-Tanner equation [18] was employed to yield the desired noisy synthetic diffusion tensor data. 3.2
Real Diffusion Tensor Data
Real diffusion tensor imaging data were acquired from a single healthy male subject. A 1.5T Siemens Sonata scanner was used to do the measurement using an optimized diffusion tensor imaging sequence described in [6]. The imaging parameters were 43 axial slices, FOV = 230 mm, TR = 6000 ms, TE = 106 ms, 2.5 mm slice thickness with 0.25 mm gap, acquisition matrix 128 × 128, and 60 images acquired at each location consisting of 16 with low diffusion weighting (b = 0) and 44 diffusion images in which the encoding gradient vectors were uniformly distributed in space (b = 1100 s/mm2 ) using the electrostatic approach described elsewhere. The reconstruction matrix was 256 × 256, resulting in an
70
J. Zhang, N. Kang, and S.E. Rose
in-plane resolution of 0.898 × 0.898 mm2 . The diffusion tensor was calculated according to the Stejskal-Tanner equation [16]. The resolution of the original calculated tensor data volume was 256 × 256 × 43 with a voxel size of 0.898 × 0.898 × 2.75 mm3 defined on a Cartesian mesh. It has been recomputed using trilinear interpolation, leading to a uniform voxel size of (0.898 mm)3 .
4
Results
Five single-turn helical fiber bundles were synthetically generated with radius being 25 mm, 20 mm, 15 mm, 10 mm, and 5 mm, respectively. For each helix, trajectories were traced from a single voxel at the lower end of the tract. In Fig. 2, the tracking result is presented, showing the simulated helical curves are closely reproduced. Fig. 3 delineates the tracing results on crossing fiber tracts synthetically constructed with two straight-line fiber bundles. It can be seen that the algorithm is capable of getting through the crossing area with planar tensors. Fig. 4 demonstrates the capability of the tracking algorithm on real human brain DT-MRI data, showing computed tracks launched from two starting voxels that are placed in the corticospinal tract approximately at the level of the left and right pons area, respectively. It is apparent that the calculated tracks emerging from the starting points branch into different cortical motor regions. In Fig. 5, the pathways reconstructed stem from three starting voxels located approximately in the midline of the corpus callosum, which interconnect the two hemispheres. Depicted in Figure 6 is another tracking example, which shows computed fiber pathways of the cingulum. We also generated pontocerebellar tracts with two different settings of the criteria set C, as shown in Figs. 7 and 8, respectively. One can observe that with appropriate selection of the thresholds, the medial lemniscus erroneously yielded in Fig. 7 can be eliminated when segmenting the pontocerebellar tracts in the pons which bear an entangled fiber crossing structure.
Fig. 2. Synthetic helical fiber tracts with varying radii (left, shown as diffusion tensor ellipsoid map) and the tracing results (right) yielded by the fiber tracking algorithm
Approximating Anatomical Brain Connectivity
71
Fig. 3. Synthetic crossing fiber tracts (left, shown as diffusion tensor ellipsoid map) and the tracing results (right) yielded by the fiber tracking algorithm
(b) (a)
(c)
(d)
(e)
Fig. 4. Fiber pathways of corticospinal tract computed from a starting voxel positioned approximately at the level of the left pons. The colors on tracks correspond to the connectivity index as shown in (e), the color bar legend. Fibers are incorporated into grey-scale fractional anisotropy (f a) maps for anatomical reference, where bright greyscale regions reflect high diffusion anisotropy. (a) Viewed from front, superimposed on a coronal f a map. (b) Viewed from left, overlaid on a midline sagittal f a map. (c) A 3D view, shown together with an axial f a map at the level of the internal capsule. (d) Viewed from top, overlaid on an axial f a map at the level of the motor cortex
72
J. Zhang, N. Kang, and S.E. Rose
(a)
(b)
Fig. 5. Fiber pathways generated from three starting voxels located in the midline of the body of the corpus callosum. The tracks are color-scaled as in Fig. 4. (a) Viewed from right, overlaid on a sagittal f a map at the midline. (b) Viewed from top-front, overlaid on an axial f a map
(a)
(b)
Fig. 6. Fiber pathways of cingulum calculated from two starting voxels which are slightly above the body of the corpus callosum. The tracks are color-coded as in Fig. 4. (a): Viewed from top, shown together with an axial f a map. (b): A lateral view from left, overlaid on a midline sagittal f a map
5
Discussion and Conclusion
We have conducted tracking experiments on synthetic as well as on real human brain diffusion tensor data, utilizing the tractography algorithm based on simulating the diffusion process in diffusion kernels. An accurate replication of the ideal track geometries has been presented in the tracing results on simulated tensor fields, while the estimated pathways on several major white matter tracts are faithful to the corresponding neuroanatomy performed with postmortem dissections and compatible to those obtained by using other reported tracking techniques. The demonstration shows that with the diffusion tensor imaging data, it is possible to employ the diffusion equation-based tracking technique to
Approximating Anatomical Brain Connectivity
(a)
73
(b)
Fig. 7. Fiber pathways of pontocerebellar tract. The tracks are color-coded as in Fig. 4. Here C = {0.8, 0.25, 0.7, 0.65, 4}. The pink pathways are corticospinal tracts as in Fig. 4. The tracks pointed to by the black arrow belong to medial lemniscus. (a): Viewed from top, shown together with an axial f a map. (b): A 3D projection with an axial and a coronal f a map superimposed
(a)
(b)
Fig. 8. Fiber pathways of pontocerebellar tract. The tracks are color-coded as in Fig. 4. Here C = {0.8, 0.2, 0.7, 0.7, 4}. The pink pathways are corticospinal tracts as in Fig. 4. (a): Viewed from top, shown together with an axial f a map. (b): A 3D projection with an axial and a coronal f a map superimposed
noninvasively follow the major white matter fiber tracts and construct maps of connectivity in the living human brain. The demonstrations have shown that the tracking algorithm has the capability of elucidating branched fiber pathways naturally from a single starting voxel, without using multiple interpolated starting points within the starting voxel or specifying regions of interest defined from anatomical landmarks. It relies on simulating the diffusion process to construct diffusion front in kernels, which is truly a physical phenomenon and the magnitude of the tensor contributes to fiber tracking, instead of fully relying on the orientation of the tensor in each voxel. So the primary advantage of the algorithm is its ability to accommodate branching fibers with the connectivity index assigned representing uncertainty.
74
J. Zhang, N. Kang, and S.E. Rose
Another desirable feature is its capability to behave correctly in crossing regions with reduced tensor information. It also reveals that the properties of the generated tracts are dependent on the threshold values in the criterion set C, which bears flexibility to improve the tracking reliability and robustness to noise. In fact, diffusion tensor imaging (DTI), as used by our diffusion equationbased tracking method, is unable to truly resolve the crossing of multiple axon directions within a single voxel. However, it has been suggested to get around the inadequacy by using newly developed imaging approaches, like high angular resolution diffusion imaging (HARDI), q-space imaging (QSI), or generalized diffusion tensor imaging (GDTI). An outstanding feature of fiber reconstruction using diffusion simulations is that it can be seamlessly adapted to the platform established by the new imaging techniques. Studies have shown that the generalized diffusion tensor model is able to not only accommodate HARDI and GDTI methods but QSI as well, due to the relationships among DTI, HARDI, GDTI, and QSI [9, 11]. This makes it possible for the diffusion simulation based tractography to become independent of the imaging techniques used, while the fiber tracking will require a more sophisticated diffusion simulation, which is governed by a generalized diffusion equation associated with generalized diffusion tensors, according to a generalization of Fick’s second law.
References 1. Basser, P.J., Pajevic, S., Pierpaoli, C., Duda, J., Aldroubi, A.: In vivo fiber tractography using DT-MRI data. Magn. Reson. Med. 44 (2000) 625-632 2. Conturo, T.E., Lori, N.F., Cull, T.S., Akbudak, E., Snyder, A.Z., Shimony, J.S., McKinstry, Burton, H., Raichle, M.E.: Tracking neuronal fiber pathways in the living human brain. Proc. Natl. Acad. Sci. USA 96 (1999) 10422-10427 3. G¨ ossl, C., Fahrmeir, L., P¨ utz, B., Auer, L.M., Auer, D.P.: Fiber tracking from DTI using linear state space models: detectability of the pyramidal tract. NeuroImage 16 (2002) 378-388 4. Gudbjartsson, H., Patz, S.: The Rician distribution of noisy MRI data. Magn. Reson. Med. 34 (1995) 910-914 5. Hagmann, P., Thiran, J.P., Jonasson, L., Vandergheynst, P., Clarke, S., Maeder, P., Meuli, R.: DTI mapping of human brain connectivity: statistical fibre tracking and virtual dissection. NeuroImage 19 (2003) 545-554 6. Jones, D.K., Horsefield, M.A., Simmons, A.: Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magn. Reson. Med. 42 (1999) 515-525 7. Jones, D.K., Simmons, A., Williams, S.C.R., Horsfield, M.A.: Non-invasive assessment of axonal fiber connectivity in the human brain via diffusion tensor MRI. Magn. Reson. Med. 42 (1999) 37-41 8. Lazar, M., Weinstein, D.M., Tsuruda, J.S., Hasan, K.M., Arfanakis, K., Meyerand, M.E., Badie, B., Rowley, H.A., Haughton, V., Field, A., Alexander, A.L.: White matter tractography using diffusion tensor deflection. Human Brain Mapping 18 (2003) 306-321 9. Liu, C., Bammer, R., Acar, B., Moseley, M.E.: Characterizing non-Gaussian diffusion by using generalized diffusion tensors. Magn. Reson. Med. 51 (2004) 924-937
Approximating Anatomical Brain Connectivity
75
10. Mori, S., Crain, B., Chacko, V.P., van Zijl, P.C.M.: Three dimensional tracking of axonal projections in the brain by magnetic resonance imaging. Ann. of Neurol. 45 (1999) 265-269 ¨ 11. 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 (2003) 955-965 12. Parker, G.J.M., Alexander, D.C.: Probabilistic Monte Carlo based mapping of cerebral connections utilising whole-brain crossing fibre information. LNCS 2737 (2003) 684-695 13. Parker, G.J.M., Wheeler-Kingshott, C.A.M., Barker, G.J.: Estimating distributed anatomical connectivity using fast marching methods and diffusion tensor imaging. IEEE Trans. Med. Imag. 21 (2002) 505-512 14. Pierpaoli, C., Basser, P.J.: Toward a quantitative assessment of diffusion anisotropy. Magn. Reson. Med. 36 (1996) 893-906 15. Poupon, C., Clark, C.A., Frouin, V., Regis, J., Bloch, I., Le Bihan, D., Mangin, J.: Regularization of diffusion-based direction maps for the tracking of brain white matter fascicles. NeuroImage 12 (2000) 184-95 16. Stejskal, E.O., Tanner, J.E.: Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient. J. Chem. Phys. 42 (1965) 288-292 17. Tournier, J.-D., Calamante, F., Gadian, D.G., Connelly, A.: Diffusion-weighted magnetic resonance imaging fibre tracking using a front evolution algorithm. NeuroImage 20 (2003) 276-288 18. Westin, C.-F., Maier, S.E., Mamata, H., Nabavi, A., Jolesz, F.A., Kikinis, R.: Processing and visualization for diffusion tensor MRI. Med. Imag. Anal. 6 (2002) 93-108 19. Wiegell, M.R., Larsson, H.B., Wedeen, V.J.: Fiber crossing in human brain depicted with diffusion tensor MR imaging. Radiology 217 (2000) 897-903
Laboratory for High Performance Scientific Computing and Computer Simulation, Department of Computer Science, University of Kentucky, Lexington, KY 40506-0046, USA {jzhang, nkang2}@cs.uky.edu http://www.cs.uky.edu/∼jzhang 2 Centre for Magnetic Resonance, University of Queensland, Brisbane, QLD, 4072 Australia [email protected]
Abstract. We present a new technique for noninvasively tracing brain white matter fiber tracts using diffusion tensor magnetic resonance imaging (DT-MRI). This technique is based on performing diffusion simulations over a series of overlapping three dimensional diffusion kernels that cover only a small portion of the human brain volume and are geometrically centered upon selected starting voxels where a seed is placed. Synthetic and real DT-MRI data are employed to demonstrate the tracking scheme. It is shown that the synthetic tracts can be accurately replicated, while several major white matter fiber pathways in the human brain can be reproduced noninvasively as well. The primary advantages of the algorithm lie in the handling of fiber branching and crossing and its seamless adaptation to the platform established by new imaging techniques, such as high angular, q-space, or generalized diffusion tensor imaging.
1
Introduction
A number of fiber tracking algorithms have been developed since the appearance of diffusion tensor magnetic resonance imaging (DT-MRI). Typical fiber tracking schemes, including the streamline-based technique [1, 2, 10], reconstruct the white matter tracts by tracing down in a voxel-by-voxel manner, using an estimate of the local fiber orientation determined by the principal eigenvector in each voxel that is assumed to align with the mean fiber direction in that voxel. These techniques appear to give excellent results in many instances if the principal eigenvector field is smooth. However, it suffers from several significant limitations. The vector field is error prone in the sense that the noise of DT-MRI data will influence the direction of the principal eigenvector, yielding an accumulation of orientational errors and thus an erroneous fork of the trajectory. Another major restriction is that these techniques may also be affected by partial volume G.E. Christensen and M. Sonka (Eds.): IPMI 2005, LNCS 3565, pp. 64–75, 2005. c Springer-Verlag Berlin Heidelberg 2005 °
Approximating Anatomical Brain Connectivity
65
effects [19], leading to unstable tracking through the primary eigenvector field in regions of fiber crossing, branching, or merging. Under the diffusion tensor imaging platform, a variety of methods have been proposed aiming to palliate the difficulties with more information incorporated from the diffusion tensor data. The algorithm presented in [8] uses a deflection term obtained from the diffusion tensor to improve the image noise immunity. Other schemes [7, 15] use predefined knowledge to group together neighboring voxels based on a similarity measure. Taking into account the uncertainty of fiber direction, probabilistic and statistical approaches [3, 5, 12] have been developed to mitigate the effects of fiber crossing and diverging as well as the sensitivity to noise. The level set theory is also utilized to find fiber paths connecting different brain regions [13]. Another front evolution algorithm proposed in [17] utilizes the fiber orientation function to reconstruct fiber tracts. As the measured quantity in DT-MRI is for water diffusion, an intuitive way to gain insights from the diffusion tensor data is to treat the brain volume as a physical system and simulate a virtual water diffusion process over it, which is anisotropic and governed by the diffusion equation. The shape of the anisotropic diffusion, represented by diffusion fronts, can be used to estimate the directional arrangement of the underlying white matter fiber bundles. This reflection is based upon the principle that the faster the diffusion, the longer the distance will be traveled on average by water molecules within the same amount of diffusion time. The fiber tracts are thus expected to proceed along the direction where the diffusion is the greatest. The fiber tractography presented in this paper performs simulations of the diffusion process stemming from a series of diffusion starting voxels, within corresponding overlapped 3D diffusion kernels. The diffusion simulation initiated from a diffusion root node is utilized to construct a diffusion front in its associated kernel. The next set of diffusion root nodes, where a seed will be placed, are located on the diffusion front which is generated by the diffusion process initiated from a previous seeded root voxel. They are picked up according to the created distance map and the local orientation information involving these voxels and the diffusion root node. For the next round, each of the newly selected diffusion root voxels will be used to generate a front by starting a diffusion process in its own kernel. Given below is the detailed description of the diffusion-based fiber tractography theory and algorithm.
2 2.1
Methods The Anisotropic Diffusion Equation
The anisotropic diffusion process simulated in this work is governed by the equation ∂C = ∇ · (D∇C), (1) ∂t where D is the so-called diffusion coefficient, which is a second-order tensor in the presence of anisotropy, C the concentration, and t the independent time variable. The coefficient used in the anisotropic diffusion simulation is the diffusion tensor
66
J. Zhang, N. Kang, and S.E. Rose
calculated from the diffusion-weighted imaging data, which is represented by a three-by-three symmetric positive definite matrix, Dxx Dxy Dxz D = Dyx Dyy Dyz , Dzx Dzy Dzz where the subscripts xx, xy, xz, etc., denote the values of the individual coefficients in the matrix. 2.2
Extracting Front in Diffusion Kernel
The first step to reconstruct fiber pathways starting from a pre-chosen root node s involves simulating the diffusion process in its associated diffusion kernel, initiated from a seed (an initial concentration value) in this voxel. A diffusion kernel defined here is a cube with six rectangular sides, which covers only a small bulk of the whole 3D data volume and is geometrically centered upon the diffusion starting voxel. The virtual concentration seed of water spreads from the root node through the neighboring nodes, within a limited amount of time, forming a diffusion front which is the surface of a diffusion volume containing nodes with nonzero concentration values. The expansion of the diffusion volume originated from the root node is achieved by integrating Eq. (1) over a certain amount of time, subject to the following initial condition, ½ 1 at the root node, = (2) C 0 elsewhere in the diffusion kernel. t=0 The boundary of the diffusion kernel is assumed to be insulated, i.e., (D∇C) · n = 0,
(3)
where n is the direction normal to the boundary. This condition implies that the normal part of the gradient of the concentration on the boundary is zero, in other words, nothing escapes out of the domain. We have developed an unsteady state anisotropic diffusion solver framework, which is adapted to the cerebral circumstance and runs in both sequential and parallel computing environments. In the current paper, (1) was solved sequentially under the initial condition (2) and boundary condition (3) by resorting to the established computational framework. We used a diffusion kernel with dimensions 11 × 11 × 7 and a voxel size same as in the original data volume, which proved to be very efficient in time cost and showed no impairment on tracking performance. Fig. 1 shows a concentration distribution map of the anisotropic diffusion simulated in a diffusion kernel using the human brain tensor data. Once the time integration for solving (1) is done, a discrete approximation to the diffusion front can be calculated in terms of whether the concentration value is zero in a voxel. Thus all nodes in the diffusion kernel can be partitioned into two groups, one with zero concentration and the other with nonzero values. Since only one seed is diffused over the root node, the diffusion-swept volume, denoted
Approximating Anatomical Brain Connectivity
67
Fig. 1. The concentration distribution map of the anisotropic diffusion simulated in a diffusion kernel (bounded by the white rectangle). The profile is superimposed on a grey-scale axial map of the fractional anisotropy
as V (r), is comprised of voxels with nonzero concentration values, where r is the position of the root node. For each member of V (r), we consider its surrounding 26 closest neighboring nodes in a 3 × 3 × 3 cube. Let i, j, k index the relative coordinates of the 26 nearest neighbors of r with i, j, k ∈ {−1, 0, 1}. If F (r) is the set of voxels that form the diffusion front of r, then for any node p ≡ (px , py , pz ) ∈ V (r), we define p ∈ F (r) if ∃ (i, j, k), such that (px − i, py − j, pz − k) ∈ / V (r), which implies that if any of the 26 nodes is not in V (r), then p ∈ F (r). 2.3
The Criteria Set
In order to store and handle the front nodes dynamically produced in each diffusion kernel, we set up a queue Q, a first-in first-out data structure. Q is initialized to contain just the starting node s, i.e., Q = {s}, thereafter, Q always contains the set of diffusion front nodes which will be the subsequent diffusion root voxels. Once F (r) is computed for the root node r, we further apply the criteria in the set C (see below) to the nodes of F (r) and pick up those that meet the corresponding thresholds. We define I(r) to be the set of nodes selected from F (r) that satisfy the criteria in C, i.e., I(r) = { p ∈ F (r) | p meets all the criteria in C }. The qualified nodes in F (r) are inserted into I(r) in a non-ascending order of α (see the first criteria c1 below). I(r) is then appended to the tail of the queue Q. The set C bears a number of criteria, which determine the connection of fiber pathways. There are five criteria in C used to evaluate the information about distance S and orientation between the root r and its front nodes in F (r). Let C = i∈{1,··· ,5} {ci }. The first criterion, c1 , is the threshold for distance ratio measure α, which is defined as α = d/dmax , where d = kv(r)k is the Euclidean distance in R3 of two points, connected by the vector v(r) pointing from r to a node in F (r), while dmax is the maximum value among the d’s. We set c2 to be a threshold of an invariant anisotropy index, the fractional anisotropy (f a) [14]. The next criterion, c3 , is a curvature constraint introduced to secure the tracks yielded moving forward consistently and smoothly without erratically turning
68
J. Zhang, N. Kang, and S.E. Rose
back on themselves. A threshold is used to restrain the angle between v(π(r)) and the current direction of tracking, v(r). Here, π(r) is the predecessor voxel of r, i.e., r ∈ I(π(r)). v(π(r)) is an established vector pointing from π(r) toward r, which implies the presence of a trajectory passing in this direction. c4 is used to judge the coherence of fiber directions along the reconstructed trajectories passing through r. One threshold is set on three inner products, φ1 , φ2 , and φ3 , where φ1 = |ˆ v (r) · e1 (r)|, φ2 = |ˆ v (r) · e1 (f )|, and φ3 = |e1 (r) · e1 (f )|. Here, ˆ (r) = v(r)/kv(r)k; e1 (r) and e1 (f ) are principal eigenvectors (corresponding v to the largest eigenvalue of D) at the voxel r and f ∈ F (r), respectively. The last criterion, c5 , specifies the maximum number of voxels I(r) allowed to have if there are more voxels than expected satisfying all previous four criteria, which controls the overall computational time for simulating the diffusion process in diffusion kernels. 2.4
Recovering Fiber Pathways
When Q is not empty, the current head node of Q is removed off the queue and is considered to be a new root r′ where a seed is diffused. r′ is positioned at the geometrical center of the diffusion kernel, which is then initialized using the global-to-local mapping to retrieve necessary information from the original data volume for carrying out the new diffusion simulation. The diffusion front F (r′ ) is calculated in the same way as that of F (r). As in the derivation of I(r), the set I(r′ ) is determined as well by checking each member of F (r′ ) based upon the criteria in C, then it is added to the tail of Q. We continue in this way by repeatedly taking off the head node of Q and processing it as a new root to diffuse a seed over it, until the queue becomes empty. Aimed at recovering the fiber pathways after constructing the diffusion front in each kernel, each voxel p in the global data grid owns a memory of its predecessor voxel, π(p), where p ∈ I(π(p)). Since every voxel in the grid can be taken as a diffusion root no more than once, π(p) is the sole predecessor of p if there is one. Thus, back propagation from the voxels on diffusion fronts by following continuously the corresponding predecessor voxels leads to paths that merge to the starting voxel s. This merging corresponds to the procedure that can be viewed in the reverse direction as fiber tracts branch outwards from s. Finally, the pathways are smoothed out using B-spline least-square approximations. The diffusion equation-based fiber tractography procedure is outlined in Alg. 1. 2.5
Connectivity Index
Since the size of the set I(r) can be larger than one for any root r, there may exist a bunch of reconstructed fiber pathways that branch outwards from the voxel where tracking starts. We utilize a heuristic connectivity index, ξ, as a confidence measure to estimate the odds that any generated path well approximates a true anatomical connection. For a given putative pathway, the index is defined as Qn i ξ = i=1 ( αi +β 2 ), where α is the distance ratio, β = (φ1 + φ2 + φ3 )/3, and n the number of diffusion kernels used to produce the pathway. The definition of ξ may
Approximating Anatomical Brain Connectivity
69
Algorithm 1. Fiber tracking using kernel-based diffusion simulations 1: specify a starting node s and initialize Q such that Q = {s} 2: while Q is not empty do 3: remove the head node r off Q and take it as a root 4: initialize the 3D diffusion kernel that is geometrically centered on r 5: get V (r) by solving the diffusion equation (1) over the diffusion kernel with the initial and boundary conditions (2) and (3) imposed 6: compute F (r), then determine I(r) and append it to the tail of Q 7: end while 8: record π values for voxels during front construction 9: retrieve fiber pathways using back propagation
be construed as evaluating how faithful the computed pathways are to follow the fastest diffusion direction, yet adjusted by coherence with local fiber orientations.
3 3.1
Data Acquisition Synthetic Tensor Fields
In order to assess fidelity and robustness of the tracking algorithm, we generated synthetic DT-MRI data with a uniform voxel size of 1 mm3 , where the true path of a fiber tract is known. The tensor field was constituted upon an anisotropic and an isotropic tensors taken out of real DT-MRI data. The shape of the diffusion tensor in synthetic fibers was described by the anisotropic one such that λ1 : λ2 : λ3 was approximately 2.5 : 1 : 1, while the isotropic one was used to forge the background of the simulated tensor field. The vector field for fiber orientations was derived by sampling discretely the trajectories which were analytically defined. To make the simulated field more realistic, an approximation to Rician noise [4] was added in the diffusion-weighted images which were calculated from the Stejskal-Tanner equation using the gradient sequence in [18] and a b-value of 1000. The noisy realization led to a signal-to-noise ratio of 10. A compact analytic solution to the Stejskal-Tanner equation [18] was employed to yield the desired noisy synthetic diffusion tensor data. 3.2
Real Diffusion Tensor Data
Real diffusion tensor imaging data were acquired from a single healthy male subject. A 1.5T Siemens Sonata scanner was used to do the measurement using an optimized diffusion tensor imaging sequence described in [6]. The imaging parameters were 43 axial slices, FOV = 230 mm, TR = 6000 ms, TE = 106 ms, 2.5 mm slice thickness with 0.25 mm gap, acquisition matrix 128 × 128, and 60 images acquired at each location consisting of 16 with low diffusion weighting (b = 0) and 44 diffusion images in which the encoding gradient vectors were uniformly distributed in space (b = 1100 s/mm2 ) using the electrostatic approach described elsewhere. The reconstruction matrix was 256 × 256, resulting in an
70
J. Zhang, N. Kang, and S.E. Rose
in-plane resolution of 0.898 × 0.898 mm2 . The diffusion tensor was calculated according to the Stejskal-Tanner equation [16]. The resolution of the original calculated tensor data volume was 256 × 256 × 43 with a voxel size of 0.898 × 0.898 × 2.75 mm3 defined on a Cartesian mesh. It has been recomputed using trilinear interpolation, leading to a uniform voxel size of (0.898 mm)3 .
4
Results
Five single-turn helical fiber bundles were synthetically generated with radius being 25 mm, 20 mm, 15 mm, 10 mm, and 5 mm, respectively. For each helix, trajectories were traced from a single voxel at the lower end of the tract. In Fig. 2, the tracking result is presented, showing the simulated helical curves are closely reproduced. Fig. 3 delineates the tracing results on crossing fiber tracts synthetically constructed with two straight-line fiber bundles. It can be seen that the algorithm is capable of getting through the crossing area with planar tensors. Fig. 4 demonstrates the capability of the tracking algorithm on real human brain DT-MRI data, showing computed tracks launched from two starting voxels that are placed in the corticospinal tract approximately at the level of the left and right pons area, respectively. It is apparent that the calculated tracks emerging from the starting points branch into different cortical motor regions. In Fig. 5, the pathways reconstructed stem from three starting voxels located approximately in the midline of the corpus callosum, which interconnect the two hemispheres. Depicted in Figure 6 is another tracking example, which shows computed fiber pathways of the cingulum. We also generated pontocerebellar tracts with two different settings of the criteria set C, as shown in Figs. 7 and 8, respectively. One can observe that with appropriate selection of the thresholds, the medial lemniscus erroneously yielded in Fig. 7 can be eliminated when segmenting the pontocerebellar tracts in the pons which bear an entangled fiber crossing structure.
Fig. 2. Synthetic helical fiber tracts with varying radii (left, shown as diffusion tensor ellipsoid map) and the tracing results (right) yielded by the fiber tracking algorithm
Approximating Anatomical Brain Connectivity
71
Fig. 3. Synthetic crossing fiber tracts (left, shown as diffusion tensor ellipsoid map) and the tracing results (right) yielded by the fiber tracking algorithm
(b) (a)
(c)
(d)
(e)
Fig. 4. Fiber pathways of corticospinal tract computed from a starting voxel positioned approximately at the level of the left pons. The colors on tracks correspond to the connectivity index as shown in (e), the color bar legend. Fibers are incorporated into grey-scale fractional anisotropy (f a) maps for anatomical reference, where bright greyscale regions reflect high diffusion anisotropy. (a) Viewed from front, superimposed on a coronal f a map. (b) Viewed from left, overlaid on a midline sagittal f a map. (c) A 3D view, shown together with an axial f a map at the level of the internal capsule. (d) Viewed from top, overlaid on an axial f a map at the level of the motor cortex
72
J. Zhang, N. Kang, and S.E. Rose
(a)
(b)
Fig. 5. Fiber pathways generated from three starting voxels located in the midline of the body of the corpus callosum. The tracks are color-scaled as in Fig. 4. (a) Viewed from right, overlaid on a sagittal f a map at the midline. (b) Viewed from top-front, overlaid on an axial f a map
(a)
(b)
Fig. 6. Fiber pathways of cingulum calculated from two starting voxels which are slightly above the body of the corpus callosum. The tracks are color-coded as in Fig. 4. (a): Viewed from top, shown together with an axial f a map. (b): A lateral view from left, overlaid on a midline sagittal f a map
5
Discussion and Conclusion
We have conducted tracking experiments on synthetic as well as on real human brain diffusion tensor data, utilizing the tractography algorithm based on simulating the diffusion process in diffusion kernels. An accurate replication of the ideal track geometries has been presented in the tracing results on simulated tensor fields, while the estimated pathways on several major white matter tracts are faithful to the corresponding neuroanatomy performed with postmortem dissections and compatible to those obtained by using other reported tracking techniques. The demonstration shows that with the diffusion tensor imaging data, it is possible to employ the diffusion equation-based tracking technique to
Approximating Anatomical Brain Connectivity
(a)
73
(b)
Fig. 7. Fiber pathways of pontocerebellar tract. The tracks are color-coded as in Fig. 4. Here C = {0.8, 0.25, 0.7, 0.65, 4}. The pink pathways are corticospinal tracts as in Fig. 4. The tracks pointed to by the black arrow belong to medial lemniscus. (a): Viewed from top, shown together with an axial f a map. (b): A 3D projection with an axial and a coronal f a map superimposed
(a)
(b)
Fig. 8. Fiber pathways of pontocerebellar tract. The tracks are color-coded as in Fig. 4. Here C = {0.8, 0.2, 0.7, 0.7, 4}. The pink pathways are corticospinal tracts as in Fig. 4. (a): Viewed from top, shown together with an axial f a map. (b): A 3D projection with an axial and a coronal f a map superimposed
noninvasively follow the major white matter fiber tracts and construct maps of connectivity in the living human brain. The demonstrations have shown that the tracking algorithm has the capability of elucidating branched fiber pathways naturally from a single starting voxel, without using multiple interpolated starting points within the starting voxel or specifying regions of interest defined from anatomical landmarks. It relies on simulating the diffusion process to construct diffusion front in kernels, which is truly a physical phenomenon and the magnitude of the tensor contributes to fiber tracking, instead of fully relying on the orientation of the tensor in each voxel. So the primary advantage of the algorithm is its ability to accommodate branching fibers with the connectivity index assigned representing uncertainty.
74
J. Zhang, N. Kang, and S.E. Rose
Another desirable feature is its capability to behave correctly in crossing regions with reduced tensor information. It also reveals that the properties of the generated tracts are dependent on the threshold values in the criterion set C, which bears flexibility to improve the tracking reliability and robustness to noise. In fact, diffusion tensor imaging (DTI), as used by our diffusion equationbased tracking method, is unable to truly resolve the crossing of multiple axon directions within a single voxel. However, it has been suggested to get around the inadequacy by using newly developed imaging approaches, like high angular resolution diffusion imaging (HARDI), q-space imaging (QSI), or generalized diffusion tensor imaging (GDTI). An outstanding feature of fiber reconstruction using diffusion simulations is that it can be seamlessly adapted to the platform established by the new imaging techniques. Studies have shown that the generalized diffusion tensor model is able to not only accommodate HARDI and GDTI methods but QSI as well, due to the relationships among DTI, HARDI, GDTI, and QSI [9, 11]. This makes it possible for the diffusion simulation based tractography to become independent of the imaging techniques used, while the fiber tracking will require a more sophisticated diffusion simulation, which is governed by a generalized diffusion equation associated with generalized diffusion tensors, according to a generalization of Fick’s second law.
References 1. Basser, P.J., Pajevic, S., Pierpaoli, C., Duda, J., Aldroubi, A.: In vivo fiber tractography using DT-MRI data. Magn. Reson. Med. 44 (2000) 625-632 2. Conturo, T.E., Lori, N.F., Cull, T.S., Akbudak, E., Snyder, A.Z., Shimony, J.S., McKinstry, Burton, H., Raichle, M.E.: Tracking neuronal fiber pathways in the living human brain. Proc. Natl. Acad. Sci. USA 96 (1999) 10422-10427 3. G¨ ossl, C., Fahrmeir, L., P¨ utz, B., Auer, L.M., Auer, D.P.: Fiber tracking from DTI using linear state space models: detectability of the pyramidal tract. NeuroImage 16 (2002) 378-388 4. Gudbjartsson, H., Patz, S.: The Rician distribution of noisy MRI data. Magn. Reson. Med. 34 (1995) 910-914 5. Hagmann, P., Thiran, J.P., Jonasson, L., Vandergheynst, P., Clarke, S., Maeder, P., Meuli, R.: DTI mapping of human brain connectivity: statistical fibre tracking and virtual dissection. NeuroImage 19 (2003) 545-554 6. Jones, D.K., Horsefield, M.A., Simmons, A.: Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magn. Reson. Med. 42 (1999) 515-525 7. Jones, D.K., Simmons, A., Williams, S.C.R., Horsfield, M.A.: Non-invasive assessment of axonal fiber connectivity in the human brain via diffusion tensor MRI. Magn. Reson. Med. 42 (1999) 37-41 8. Lazar, M., Weinstein, D.M., Tsuruda, J.S., Hasan, K.M., Arfanakis, K., Meyerand, M.E., Badie, B., Rowley, H.A., Haughton, V., Field, A., Alexander, A.L.: White matter tractography using diffusion tensor deflection. Human Brain Mapping 18 (2003) 306-321 9. Liu, C., Bammer, R., Acar, B., Moseley, M.E.: Characterizing non-Gaussian diffusion by using generalized diffusion tensors. Magn. Reson. Med. 51 (2004) 924-937
Approximating Anatomical Brain Connectivity
75
10. Mori, S., Crain, B., Chacko, V.P., van Zijl, P.C.M.: Three dimensional tracking of axonal projections in the brain by magnetic resonance imaging. Ann. of Neurol. 45 (1999) 265-269 ¨ 11. 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 (2003) 955-965 12. Parker, G.J.M., Alexander, D.C.: Probabilistic Monte Carlo based mapping of cerebral connections utilising whole-brain crossing fibre information. LNCS 2737 (2003) 684-695 13. Parker, G.J.M., Wheeler-Kingshott, C.A.M., Barker, G.J.: Estimating distributed anatomical connectivity using fast marching methods and diffusion tensor imaging. IEEE Trans. Med. Imag. 21 (2002) 505-512 14. Pierpaoli, C., Basser, P.J.: Toward a quantitative assessment of diffusion anisotropy. Magn. Reson. Med. 36 (1996) 893-906 15. Poupon, C., Clark, C.A., Frouin, V., Regis, J., Bloch, I., Le Bihan, D., Mangin, J.: Regularization of diffusion-based direction maps for the tracking of brain white matter fascicles. NeuroImage 12 (2000) 184-95 16. Stejskal, E.O., Tanner, J.E.: Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient. J. Chem. Phys. 42 (1965) 288-292 17. Tournier, J.-D., Calamante, F., Gadian, D.G., Connelly, A.: Diffusion-weighted magnetic resonance imaging fibre tracking using a front evolution algorithm. NeuroImage 20 (2003) 276-288 18. Westin, C.-F., Maier, S.E., Mamata, H., Nabavi, A., Jolesz, F.A., Kikinis, R.: Processing and visualization for diffusion tensor MRI. Med. Imag. Anal. 6 (2002) 93-108 19. Wiegell, M.R., Larsson, H.B., Wedeen, V.J.: Fiber crossing in human brain depicted with diffusion tensor MR imaging. Radiology 217 (2000) 897-903
