Improved fiber tractography with Bayesian tensor regularization

www.elsevier.com/locate/ynimg NeuroImage 31 (2006) 1061 – 1074

Improved fiber tractography with Bayesian tensor regularization Yonggang Lu,a,* Akram Aldroubi,b John C. Gore,a Adam W. Anderson,a,c and Zhaohua Dinga a

Vanderbilt University Institute of Imaging Science, USA Department of Mathematics, Vanderbilt University, USA c Department of Biomedical Engineering, Vanderbilt University, USA b

Received 18 August 2005; revised 24 January 2006; accepted 25 January 2006 Available online 24 March 2006

Diffusion tensor tractography suffers from the effects of noise and partial volume averaging (PVA). For reliable reconstruction of fiber pathways, tracking algorithms that are robust to these artifacts are called for. To meet this need, the present study establishes a novel Bayesian regularization framework for fiber tracking that takes into account the effects of noise and PVA, thereby improving tracking accuracy and precision. With this framework, the propagation of a fiber path follows an optimal vector determined by Bayes decision rule; the probability functions involved are modeled on the basis of multivariate normal distributions of diffusion tensor elements, which allows the optimal solution with maximum a posteriori probability to be derived analytically. Parameters for the probability functions are estimated from the uncertainty of tensor elements and the variance among tensors within an oriented sampling volume weighted by fractional anisotropy. Experiments with Monte Carlo simulations, synthetic, and in vivo human diffusion tensor data demonstrate that this specialized scheme enhances the immunity of fiber tracking to noise and PVA, and hence enables fibers to be more faithfully reconstructed. D 2006 Elsevier Inc. All rights reserved. Keywords: Diffusion tensor imaging; Fiber tracking; Bayes decision rule; Regularization

Introduction Diffusion tensor imaging (DTI) based fiber tractography is a promising technique for characterizing fiber pathways in vivo (Le Bihan et al., 2001; Mori and van Zijl, 2002). To date, a number of fiber tracking algorithms have been proposed, which can be categorized into the following: (1) deterministic fiber tracking which provides a single connection between voxels, including

* Corresponding author. Wiscom System Co. Ltd., 100 Jiangjun Road, Jiangning District, Nanjing 211100, China. Fax: +86 25 52762929. E-mail address: [email protected] (Y. Lu). Available online on ScienceDirect (www.sciencedirect.com). 1053-8119/$ - see front matter D 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2006.01.043

simple fiber tracking that only follows the direction of the measured major eigenvector (Basser et al., 2000; Mori et al., 1999) and complex fiber tracking with built-in immunity to the effect of noise and/or partial volume averaging (PVA) (Go¨ssl et al., 2002; Lazar et al., 2003; Parker et al., 2002; Poupon et al., 2000; Staempfli et al., 2005; Tench et al., 2002; Zhukov and Barr, 2002); (2) probabilistic fiber tracking which provides maps of fiber connectivity (Hagmann et al., 2003; Koch et al., 2002; Parker et al., 2003). These have offered a repository of utilities for studying the connectivity of fibrous tissues such as white matter tracts in the brain. In spite of the variety of methods that come with different levels of sophistication, problems remain due to the ill-posed nature of fiber tracking, i.e., recovering linear structure from a noisy and voxel averaged tensor field. Therefore, the practice of fiber tractography still primarily relies on more basic methods such as the Euler (Basser et al., 2000) or FACT method (Mori et al., 1999), which are favored as a trade-off between reliability and computational complexity. These streamline methods simply propagate fiber paths by following the measured major eigenvector of tensors with the assumption that it represents the fiber local tangent direction. Hence, they suffer from a few limitations, the most notable of which is sensitivity to the effects of noise and PVA. Noise and PVA are inherent artifacts in DTI data; the former may significantly deflect the measured major eigenvector from its true direction, and the latter creates ambiguities in defining the direction of the underlying fascicle. Thus, the major eigenvector measured may be an incorrect estimate of the true local tangent direction of the underlying fascicle and is necessarily associated with a certain degree of uncertainty. Fiber tracking by simply following the direction of the measured major eigenvector is therefore error prone and frequently yields grossly aberrant pathways. The limitations of basic streamline fiber tracking have often led to questions of the fidelity and reproducibility of reconstructed fiber pathways. This study is thus motivated with the goal to develop a robust fiber tracking algorithm that takes into account the uncertainty of the diffusion tensor caused by image noise and PVA, aiming at improving the accuracy and precision

1062

Y. Lu et al. / NeuroImage 31 (2006) 1061 – 1074

of fiber tracking in the presence of these artifacts. This goal is pursued by tensor regularization using the traditional Bayes decision rule, but with novel formulations of the probability functions involved. Briefly, the uncertainty of the diffusion tensor measured and the variance among tensors within a sampling volume are incorporated into models of the conditional and a priori probabilities on the basis of multivariate normal distributions. The true local tangent direction of the underlying fascicle is defined to be an optimal solution to the Bayes decision problem that has a maximum a posteriori probability. The Bayesian tracking method proposed allows fiber pathways to be reconstructed more faithfully than the Euler method. In addition, the new formulation of probability functions naturally leads to the optimal solution with the maximum a posteriori probability derived analytically; this has important implications for the issues of computational complexity and convergence typically encountered in optimization procedures. The remainder of the paper is organized as follows. Section 2 describes the framework of the proposed algorithm in detail. Section 3 presents the results from Monte Carlo simulations and from experiments using synthetic and in vivo human DTI data. Section 4 summarizes the main contributions of this work and discusses some relevant technical issues, followed by concluding remarks.

The Bayes decision is the process of determining the value m that maximizes the a posteriori probability P(m j |d), i.e.:   1 0     p djmj I P mj A ¼ arg max @ m ¼ arg max P mj jd mj mj pðd Þ       ¼ arg max p djmj I P mj ; j ¼ 1; 2N ; n : mj

ð3Þ

Note in the above equation, since p(d) does not depend on m j , maximizing P(m j |d ) amounts to maximizing the product p(d|m j ) I P(m j ). In contrast to other methods such as nearest decision, Bayes decision is an optimal one in that it incorporates all available data and makes a decision with a minimum error rate (Duda and Hart, 1973). Bayesian fiber tracking

Methods The Bayesian fiber tracking method proposed is essentially a process of line propagation with an integration scheme similar to that used in the Euler fiber tracking method (Basser et al., 2000), i.e.: si þ 1 ¼ si þ a I ei ;

conditional probability density of a variable d if the state is m j . The a posteriori probability P(m j |d), i.e., the probability of being in the state m j given observed variable d, can be computed with the following equation:       p djmj I P mj : ð2Þ P mj jd ¼ pðd Þ

ð1Þ

where s i is a position vector at discrete step i, a the step size, and ( i a propagation vector. The propagation vector ( is the key parameter that determines the outcome of fiber tracking. Intuitively the major eigenvector of the diffusion tensor may be used as the propagation vector, but as mentioned, the major eigenvector may be an incorrect estimate of the true direction of the underlying fascicle due to inherent image noise and PVA, and therefore provides an unreliable basis for fiber path propagation. To reliably track fiber pathways with Eq. (1), the true direction of the underlying fascicle needs to be inferred from its estimate (Poupon et al., 2000; Tench et al., 2002). In this study, this problem is approached by using the traditional Bayes decision rule which gives an optimal solution that has a minimum rate of error. In the sections below, the Bayes decision rule is first briefly reviewed, followed by formulation of fiber tracking as a Bayes decision problem. Models of the probability functions necessary for the Bayes decision rule are presented, and the analytical solution that maximizes the a posteriori probability is provided. Finally, estimations of relevant model parameters are given. Bayes decision rule Let m j denote a possible state of nature ( j = 1, 2, . . ., n), P(m j ) is the a priori probability of state m j , and p(d|m j ) is the

Fiber tracking with Eq. (1) can be cast as a decision problem to which an optimal solution is determined by Bayes decision rule. In this context, the goal is to infer from an estimate (observed variable) the true fiber direction (true state of nature), which is presumably the solution with maximum a posteriori probability. To determine the true direction, the a priori and conditional probability functions need to be known or defined. These functions are not known a priori, and limited information can be gained from measured DTI data. In such cases, modeling is commonly used to approximate the true probability functions. In this study, the probability functions are modeled using the multivariate normal distribution. As will soon become clear, the convenient assumption of normality offers important practical advantages in the determination of an optimal solution. The Bayes decision rule is applied to the diffusion tensor instead of its major eigenvector. This avoids certain complications such as sorting bias, symmetrical ambiguity and error estimation of the major eigenvector. To begin with, let us reorganize the six independent elements of each diffusion tensor into a column vector (called tensor element vector henceforth): d = (D 11 D 22 D 33 D 12 D 13 D 23)t , and assume d to be multivariate normal: d ¨ N(m,S), where m represents the unknown true tensor element vector, and S is the covariance matrix. Variance of the measured tensor element vector reflects the uncertainty associated with noise and PVA. The conditional probability density of d given m j is given by    t   1 1 1 d  mj S exp  d  mj ; p djmj ¼ P 2 ð2pÞ3 j j1=2 ð4Þ where t denotes matrix transpose, S 1 is the inverse of S, and |S| is its determinant. The a priori probability P(m j ) is also assumed to be multivariate normal: m j ¨ N(m, F), where m is the mean tensor element vector and

Y. Lu et al. / NeuroImage 31 (2006) 1061 – 1074

1063

F is the covariance matrix. Similarly to p(d|m j ), P(m j ) has the following form:    t   1 1 1 m exp  ;  m F m  m : ð5Þ P mj ¼ j j 2 ð2pÞ3 jFj1=2

propagation vector should propagate the fiber path with a minimum rate of error and decreased uncertainty. Therefore, fiber pathways generated with this Bayesian tracking method should have simultaneously improved accuracy and precision.

Therefore, given a measured tensor element vector d, the a posteriori probability of being the true tensor element vector m j is:

Parameter estimation

pðmj jdÞ”pðdjmj ÞPðmj Þ ¼

1 ð2pÞ3 jRj1=2 

  1 exp  ðd  mj Þt R1 ðd  mj Þ 2

  1 t 1 ðm exp   mÞ F ðm  mÞ j 2 j ð2pÞ3 jFj1=2 1

 1  t   1 mj  S1 þ F1 S1 d þ F1m ” exp  2   1    S1 d þ F1 m  S1 þ F1  mj  S1 þ F1 : ð6Þ Consequently, P(m j |d) is still a normal distribution: P(m j |d) ¨ N(W,W) The expected value (mean) W and covariance matrix W are: 

1 1 1 W ¼ S1 þ F1 S d þ F m ; W1 ¼ S1 þ F1 :

Bayesian fiber tracking as described above involves probability functions for which model parameters need to be estimated. These parameters include the covariance matrix S of the tensor element vector for the conditional probability density, and the mean m and covariance matrix F of m j for the a priori probability. In this study, the conditional probability density is defined on the basis of the uncertainty of each measured tensor, and the a priori probability is defined on the basis of the variance among tensors within a sampling volume that is oriented with the direction of fiber propagation. Methods for estimating the model parameters are described below. Estimation of S The covariance matrix S of a tensor element vector d at each voxel can be estimated directly from measured diffusion-weighted signals, experiment parameters and signal noise profile. In a pulsed gradient spin echo experiment, the pixel signal intensity for a diffusion-weighted image (DWI) with normally distributed noise g ¨ N(0, r g2) is (Anderson, 2001; Basser et al., 1994): ! 3 X

 ˜ bij Dij þ g S b ¼ S0 exp  i;j ¼ 1

ð7Þ

¼ S0 expð  b11 D11  b22 D22  b33 D33  2b12 D12 It immediately follows that the solution with the maximum a posteriori probability can be found analytically without resorting to an optimization procedure, i.e.:

1 1  m ¼ W ¼ S1 þ F1 S d þ F1 m :

ð8Þ

This analytical solution has several important practical advantages. First, the need for searching for an optimal solution through a time-consuming optimization procedure can be eliminated, which gives considerable savings computationally. Second, convergence is not an issue for the analytical solution, although it may be for numerical optimization procedures. Third, the analytical solution is continuous, while optimization procedures typically give discretized, spatially limited solutions (Poupon et al., 2000; Tench et al., 2002). Eq. (7) also interestingly shows that W is a weighted sum of d and m, and |W| is always smaller than |S| and |F| (since W, S and F are all nonnegative definite). This indicates that the optimal solution, the one with a minimum rate of error, has a decreased uncertainty as well. The tensor matrix corresponding to the optimal solution m (l 1, l 2, l 3, l 4, l 5, l 6) is: 2 3 1 4 5 ð9Þ D 0 ¼ 4 4 2 6 5: 5 6 3 The major eigenvector of D 0 presumably represents the dominant fiber orientation within the voxel with maximum a posteriori probability. Fiber tracking with this eigenvector as the

 2b13 D13  2b23 D23 Þ þ g;

ð10Þ

where D ij are elements of a diffusion tensor matrix D, and b ij are elements of a diffusion weighting matrix b˜ (Mattiello et al., 1994). Taking the logarithm of both sides of Eq. (10), and defining f(b˜ ) = ln(S)  ln(S 0), s = g/ ( is the statistical expection of S(b˜ )), then we have:

 f b˜ ¼  b11 D11  b22 D22  b33 D33  2b12 D12  2b13 D13  2b23 D23 þ s:

ð11Þ

If we define B to be a design matrix of n independent experiments with each row corresponding to one experiment and having the form (b 11, b 22, b 33, b 12, b 13, b 23), then the relation between a n  1 signal vector f and a n  1 noise vector t is: f ¼ B I d þ t:

ð12Þ

A weighted least-square solution to Eq. (12) gives the estimated vector d and covariance matrix S respectively as (Anderson, 2001; Basser et al., 1994): d ¼ ðBt Ss1 BÞ1 I Bt S1 s If;

t 1 1 S ¼ B Ss B :

ð13Þ ð14Þ

Assuming noise at each experiment is independent, S s is a diagonal matrix with diagonal elements:

 ðSs Þnn ¼ r2g =2 ; ð15Þ where b˜ n is the b˜ matrix for the nth measurement.

1064

Y. Lu et al. / NeuroImage 31 (2006) 1061 – 1074

Estimation of m and F The a priori probability expresses the probability of a fiber orientation in the absence of experimental data at the particular point of interest. In this study, the a priori probability is estimated from the tensor element vectors measured in the voxels neighboring the point of interest. This is done by assuming the tensor element vectors in a sampling volume are random samples of a population that has a multivariate normal distribution. The sampling volume is designed to be a parallelepiped oriented with the direction of fiber propagation. This oriented sampling volume presumably contains tensors from the same population (i.e., from the same fiber bundle) and hence satisfies the distribution condition. Also it has a relatively large size which allows the model parameters to be more reliably estimated. To suppress the confounding effect of PVA, due to, for instance, fiber crossing or branching, the tensor element vectors in the sampling volume are weighted with fractional anisotropy (FA) of the tensor as the weighting factor. As FA is an indicator of fiber directional preference, tensors representing greater PVA tend to have smaller FA and thus should carry less weight in the calculation of model parameters. With FA weighting, the mean m and covariance matrix F (with element / ij ) are estimated as follows: m¼

1 L

~ FAl

L

~ FAl d l ;

ith element of m. L is the sample size which is 16 according to the design in this study. Fig. 1 shows two examples of oriented sampling volumes. A sampling volume at a point contains eight nearest voxels (first level neighborhood) and eight voxels in its second level neighborhood along the direction of the propagation vector. The sampling volume is determined using a ‘‘face-intersecting’’ approach. Specifically, for a given propagation vector, a line is drawn bidirectionally along the vector. The two points where this line intersects the plane formed by the voxels in the second level neighborhood are identified. The four closest voxels to each intersecting point are included in the sampling volume as the second level neighborhood. It should be noted that the sampling volume is chosen to be oriented with the direction of fiber propagation (see the section below) instead of the local major eigenvector, since the former is an optimally determined path direction, whereas the latter is a potentially erroneous estimate. Fiber tracking and visualization

ð16Þ

With probability functions defined and model parameters estimated, fiber tracking may be launched bidirectionally from seed points in a preselected region of interest. Implementation details at each step of line propagation are as follows:

ð17Þ

(1). Compute the tensor element vector d and its covariance matrix S at the current point according to Eqs. (13) and (14). (2). Estimate the mean m and covariance matrix F of the tensor element vector according to Eqs. (16) and (17). (3). Calculate the optimal tensor element vector m given by the Bayes decision rule according to Eq. (8). (4). Derive the major eigenvector, ei, of the tensor matrix corresponding to m. (5). Proceed to the next point with a small step size along direction ei according to Eq. (1).

l¼1

l¼1

/ij ¼

1 L

~ FAl

l¼1

L

  ~ ð FAl Idil  mi Þ FAl Idjl  mj ; ði; j ¼ 1; 2; N ; 6Þ; l¼1

where d l is the lth tensor element vector in the sampling volume, d il (i = 1, 2, . . ., 6) is the ith element of d l , and m i (i = 1, 2,. . .,6) is the

Fig. 1. Two examples of oriented sampling volumes. Each sampling volume contains eight voxels in the first level neighborhood (denoted as solid spheres) and eight voxels in the second level neighborhood (denoted as solid ellipsoids) around point O. Example one contains eight voxels denoted as solid spheres and eight voxels denoted as horizontal ellipsoids in the second level neighborhood. A and B are the points where a line along the direction of the propagation vector intersects the plane formed by voxels in the second level neighborhood in the forward and backward direction respectively. Example two contains eight voxels denoted as solid spheres and eight voxels denoted as vertical ellipsoids in the second level neighborhood. AVand BV are the points where a line along the direction of the propagation vector intersects the plane formed by voxels in the second level neighborhood in the forward and backward direction respectively.

The above procedures are iterated until certain termination criteria are reached (e.g., low FA, large curvature, or low image intensity). During the process of fiber path propagation, a simple linear interpolation is used to define the local values of the tensor element vector d, covariance matrix S, and FA. In each iteration, the sampling volume is oriented with the propagation vector from the previous iteration except for the seed point, for which the sampling volume is oriented with the major eigenvector. All the programs were written in Matlab 6.5 (The MathWorks Inc.) running on a Linux workstation. A visualization tool developed in-house based on IDL6.0 and windows XP was used for 3D fiber visualization (Li, 2001). FA maps were used as the 3D volume image to assist in visualizing the fibers. Tracking experiments To evaluate comprehensively the performance of the Bayesian tracking method proposed, four synthetic datasets with increasing complexity were designed. To mimic closely in vivo conditions, the synthetic tensors were constructed to have a trace of 2.1  105 cm2/s, comparable to that in normal brain parenchyma; diffusion weighting was simulated along six noncollinear directions with a b value of 1000 s/mm2. In addition to the synthetic

Y. Lu et al. / NeuroImage 31 (2006) 1061 – 1074

1065

levels of zero mean Gaussian noise were added to the original noise free synthetic data, generating noisy data with a signal-tonoise ratio (SNR) of 10, 20, 30, and 40 respectively (for the nondiffusion-weighted images). For each combination of eigenvalue contrast and SNR, a ‘‘fiber’’ was tracked 1000 times from the same seed point with a step size of one voxel, each time with a different realization of Gaussian noise. To measure tracking accuracy and precision, the mean displacement and standard deviation of ‘‘fiber tracts’’ in terms of Euclidean distance to the ‘‘true tract’’ were calculated point-wise along the tracking direction. The mean displacement is inversely proportional to tracking accuracy, and the standard deviation is inversely proportional to tracking precision. Test of tracking robustness to PVA with synthetic data For this purpose, a second synthetic dataset consisting of a region of PVA sandwiched by straight ‘‘fibers’’ was designed. Tensors for the straight ‘‘fibers’’ were cylindrically symmetric (k 1 > k 2 = k 3) but for the PVA region were spherically symmetric (k 1 = k 2 = k 3). The PVA region simulates the worst scenario of complete loss of directional preference by crossing of three orthogonal bundles or other likely configurations. The width of the PVA

Fig. 2. Variation of tract mean displacement (a) and standard deviation (b) with the number of tracking steps. The eigenvalue contrast is 0.8  105 cm2/s and SNR varies from 10 to 40 in increments of 10. Line styles denote different levels of SNR, and symbols denote different tracking methods.

data, a set of in vivo DTI data was acquired from a healthy human volunteer. Fiber tracking on all these datasets was performed using both Bayesian and Euler methods, and comparisons between them were made with respect to tracking accuracy and precision, robustness to noise and PVA, curvature dependency, and tract smoothness. To compare with other regularization methods aimed at achieving the same goal, fiber tracking on the synthetic datasets was also carried out using the well-known TEND method (tensor deflection method, Lazar et al., 2003), and quantitative comparisons of the above measures were made among the Euler, Bayesian, and TEND methods. Test of tracking accuracy and precision with synthetic data The first synthetic dataset designed contained parallel, vertically oriented ‘‘fiber tracts’’. Tensors for these tracts are identical with each being cylindrically symmetric (k 1 > k 2 = k 3) and having the same eigenvalue contrast (defined as Dk = k 1  k 3). Monte Carlo simulations were performed to test the effect of noise on fiber tracking accuracy and precision and the effect of eigenvalue contrast. In these simulations, the eigenvalue contrast was varied from 0.4  105 cm2/s to 1.4  105 cm2/s with an increment of 0.2  105 cm2/s. At each eigenvalue contrast, four different

Fig. 3. Variation of tract mean displacement (a) and standard deviation (b) at 100 tracking steps with eigenvalue contrast. Symbols and line styles are the same as in Fig. 2.

1066

Y. Lu et al. / NeuroImage 31 (2006) 1061 – 1074

region ranged from 1 to 4 voxels. Similar to the abovementioned, Monte Carlo simulations were performed with the eigenvalue contrast of straight ‘‘fibers’’ varied from 0.4  105 cm2/s to 1.4  105 cm2/s with an increment of 0.2  105 cm2/s, and zero mean Gaussian noise with SNR of 10, 20, 30, and 40 respectively added to the synthetic data. For each combination of eigenvalue contrast and SNR, a ‘‘fiber’’ was tracked with a step size of one voxel from the same start point 1000 times each with a different realization of Gaussian noise. To measure tracking robustness to PVA, a target point was specified at the other side of the PVA region, as shown in Fig. 4(a). The PVA robustness was defined to be the percent number of ‘‘fiber tracts’’ reaching within one voxel’s distance to the target point. A higher percentage represents a greater capability to penetrate the PVA region. Test of tracking curvature dependency with synthetic data A common limitation to regularized fiber tracking is the bias toward straight paths and hence a loss of sensitivity to follow subtle changes in the fiber direction. To examine the bias and sensitivity of the Bayesian fiber tracking proposed, a third synthetic dataset was designed to assess the dependency of tracking accuracy and precision on fiber curvature. The dataset

was comprised of several quadrantal, concentric ‘‘fibers’’ at different curvatures, with the radius ranging from 3 to 12 voxels (curvature = 1/radius), as shown in Fig. 5(a). Simulated tensors for the ‘‘fibers’’ were cylindrically symmetric (k 1 > k 2 = k 3) with an eigenvalue contrast of 0.8  105 cm2/s. Gaussian noise was added to the dataset to yield an SNR of 30. For each noise realization, fiber tracking was launched from six seed points (A, B, C, D, E, F in Fig. 5(a)) with a step size of 0.5 voxel. After 1000 times of noise realization and fiber tracking, the mean displacement of the endpoint of each putative ‘‘fiber’’ from its target point (AV, BV, C ,V DV, EV, FV respectively) and the standard deviation of the displacement were calculated. These measures were used to assess the dependency of tracking accuracy and precision on fiber curvature. Test of tracking error, noise sensitivity, and tract smoothness with synthetic data A fourth synthetic dataset was designed to contain six identical layers of straight and low planar curvature ‘‘fibers’’. A PVA region of two voxel width was synthesized around the middle portion of the straight and curved ‘‘fibers’’ respectively, as illustrated in Figs. 6(a,b) and 7(a,b). Simulated tensors for the straight ‘‘fibers’’ were cylindrically symmetric (k 1 > k 2 = k 3) with an eigenvalue contrast

Fig. 4. Variation of tracking robustness to PVA with eigenvalue contrast. Variation is shown for PVA region widths of one (b), two (c), three (d), and four (e) voxels. Locations of the start point, PVA region and target point are illustrated in panel a. Symbols and line styles are the same as in Fig. 2.

Y. Lu et al. / NeuroImage 31 (2006) 1061 – 1074

of 0.8  105 cm2/s, but for the PVA regions, the tensors were spherically symmetric (vide ante). For each PVA case, two noisy datasets at SNR of 20 and 30, respectively, were obtained by corrupting the original data with zero mean Gaussian noise. Fiber tracking on these noisy datasets was launched from four seed points as shown in Figs. 6(a,b) and 7(a,b) with a step size of 0.5 voxel. Four ‘‘true tracts’’ were obtained by tracking from the seed points with the Euler method on the noise and PVA free dataset. Tracking error, noise sensitivity, and tract smoothness, as defined in the Appendix A, were measured to quantitatively evaluate and compare the performance of the Euler, Bayesian, and TEND methods.

1067

steps. In this figure, the eigenvalue contrast is 0.8  105 cm2/s, and the SNR varies from 10 to 40 with an increment of 10. It can be seen that tract mean displacement and standard deviation from the Euler, Bayesian, and TEND methods increase with the

Fiber tracking with in vivo human DTI data To assess the performance of the Bayesian fiber tracking on in vivo data, a set of DWI data was acquired from a healthy human volunteer on a 3-T GE Signa MR scanner using a single shot, echoplanar pulsed gradient spin-echo imaging sequence. Diffusion weighting was performed along six noncollinear directions with b value of 1000 s/mm2. Timing parameters were TR of 9000 ms and TE of 88.5 ms. The imaging field of view was 240  240 mm2. Thirty contiguous, 4-mm-thick slices with a matrix size of 128  128 were acquired and subsequently interpolated into a matrix size of 256  256, yielding an in-plane pixel size of 0.94  0.94 mm2. A total of 37 repeated scans were obtained, and images showing motion and other artifacts were discarded. This high SNR (>100) DWI dataset was used for tensor calculation. Diffusion tensor elements were fitted using a weighted least square approach (Basser et al., 1994), from which FA maps were computed. T1weighted images were acquired as well during the same session for generation of brain mask images. Prior to the study, the subject gave informed consent for the study protocol that was approved by the local ethics committee. As it is not possible to measure quantitatively the accuracy of fiber tracking on in vivo data due to the lack of a ‘‘gold standard’’, performance evaluations and method comparisons were based on qualitative judgment of fiber tracts reconstructed and quantitative assessment of tracking noise sensitivity. To quantify noise sensitivity, a noisy DWI dataset was first generated by adding zero mean, 0.05 standard deviation Gaussian noise to the original DWI data. Fibers were then tracked from the noisy as well as the high SNR dataset using both the Bayesian and Euler methods. Finally, for each method, five fibers randomly selected from the high SNR reconstruction and five seed point matched fibers selected from the noisy reconstruction were used for quantification of noise sensitivity. This study focused on one major fiber bundle: the superior longitudinal fasciculus. This tract was chosen because it has wellknown anatomy, thus making it possible to judge the quality of reconstructed fibers and identify any erroneous connections. Furthermore, it passes through areas of high and low fractional anisotropy (possible area of fiber crossing), thus providing a relatively rigorous test of fiber tracking methods.

Results Tracking accuracy and precision Figs. 2(a – b) show respectively the variation of tract mean displacement and standard deviation with the number of tracking

Fig. 5. Variation of the mean displacement (b) of the tract endpoint and its standard deviation (c) with curvature. (a) Tracts at one realization of SNR = 30. A – F are seed points, and AV– FVare their corresponding target points. Red, blue and green curves are ‘‘fibers’’ from the Euler, Bayesian, and TEND methods respectively.

1068

Y. Lu et al. / NeuroImage 31 (2006) 1061 – 1074

number of tracking steps but decrease with increasing SNR. These trends attest to the cumulative nature of the impact of noise on fiber tracking, and that this impact is proportional to the noise level (Anderson, 2001). At all tracking steps and levels of SNR, the TEND method has comparable or even larger mean displacement and standard deviation (most notably at SNR = 10) than the Euler method, but the Bayesian method has much smaller values (less than half) of these measures from the other two methods. Fig. 3 shows the tract mean displacement and standard deviation at 100 tracking steps at several levels of eigenvalue contrast. In this figure, the SNR varies from 10 to 40, and the eigenvalue contrast varies from 0.4  105 cm2/s to 1.4  105 cm2/s with an increment of 0.2  105 cm2/s. Again, it shows that, at all eigenvalue contrasts and noise levels, the Bayesian method significantly outperforms the Euler and TEND method, with ¨50% overall improvement in both the mean displacement and standard deviation. In particular, at SNR of 30 or higher, the Bayesian method has ¨0.5 voxel mean displacement and ¨0.25 voxel standard deviation at all contrast levels; this demonstrates superior performance of Bayesian tracking in both accuracy and precision. It can also be observed from Fig. 3 that the tract mean displacement and standard deviation from all three methods tend to

decrease with increasing eigenvalue contrast. This is because deviations of the major eigenvector are proportional to tensor component errors divided by the eigenvalue contrast (Anderson, 2001). Tracking robustness to PVA Figs. 4b – e show the PVA robustness of the Euler, Bayesian and TEND methods at four different noise levels versus eigenvalue contrast at PVA width of one to four voxels respectively. It shows that PVA robustness of the Euler method is reasonably high (and increases with eigenvalue contrast) at PVA width of only one voxel and SNR 20, but consistently poor (from ¨20% to ¨7%) at PVA width of two or more voxels, regardless of SNR or eigenvalue contrast. The PVA robustness of the TEND method for all PVA widths is poor at low SNR (20) but is significantly improved with increased SNR and eigenvalue contrast. The Bayesian method, however, has a higher PVA robustness than the Euler and TEND method at all PVA widths and all noise levels (except at PVA width of one voxel and SNR of 10). In particular, when SNR is 30 or above, the Bayesian method has nearly 100% PVA robustness for all PVA widths and all eigenvalue contrasts. These results clearly demonstrate that the Bayesian fiber tracking proposed has greatly enhanced PVA robustness.

Fig. 6. Axial (a, b) and sagittal (c, d) views of ‘‘fibers’’ tracked from synthetic data at SNR of 30 (a, c) and 20 (b, d) with a PVA region around the middle portion of the straight ‘‘fibers’’. Red, blue, and green curves are ‘‘fibers’’ from the Euler, Bayesian, and TEND methods respectively, and black curves are the ‘‘true fibers’’. Four seed points are denoted with black dots in panels a, b. The black line segments denote the direction of the major eigenvector at each voxel. Note that the abscissa and ordinate are not drawn to the same scale in panels c and d.

Y. Lu et al. / NeuroImage 31 (2006) 1061 – 1074

1069

Fig. 7. Axial (a, b) and sagittal (c, d) views of ‘‘fibers’’ tracked from synthetic data at SNR of 30 (a, c) and 20 (b, d) with a PVA region around the middle portion of the curved ‘‘fibers’’. See captions in Fig. 6 for explanations.

Tracking curvature dependency Figs. 5b and c show the variation of mean displacement and standard deviation of the Euler, Bayesian and TEND methods with ‘‘fiber’’ curvature. Both measures from all the three methods increase with curvature, indicating all the methods track fibers with low curvature more faithfully than with high curvature. The Bayesian tracking has lower values of these measures than the other two methods, which suggests it is more sensitive to the change of fiber directions. Fig. 5a shows reconstructed tracts after one realization of Gaussian noise with SNR = 30. It can be seen that ‘‘fibers’’ from all three methods are biased outward (toward the low curvature side). It is not surprising to see this phenomenon with the Bayesian and TEND methods, since the regularized fiber

tracking tends to have a ‘‘straightening effect’’. The fact that the Euler method also has the straightening effect is primarily attributable to the relatively large step size used (see the next section for more discussions). Tracking error, noise sensitivity, and tract smoothness Figs. 6 and 7 show the axial (left columns) and sagittal (right columns) views of four ‘‘fiber tracts’’ tracked with the Euler (red lines), Bayesian (blue lines), and TEND (green lines) methods from the synthetic data with a PVA region around the middle portion of the straight and curved ‘‘fibers’’ respectively. The top rows have an SNR of 30, and the bottom rows have an SNR of 20. In both figures, all the ‘‘fibers’’ successfully penetrate the

Table 1 Comparison of tracking error, noise sensitivity, and tract smoothness among the Euler, Bayesian, and TEND methods with synthetic data Euler

Tracking error (voxels) Tracking noise sensitivity (voxels) Tract smoothness

Without PVAa With PVAb

Bayesian

TEND

SNR = 20

SNR = 30

SNR = 20

SNR = 30

SNR = 20

SNR = 30

0.21 3.11 (0.55) 5.30 (0.56)

0.11 2.66 (2.07)

0.15 0.35 (0.35) 0.45 (0.42)

0.07 0.30 (0.20)

0.22 0.32 (0.45) 0.62 (0.60)

0.11 0.61 (0.34)

0.9906 (0.9944)

0.9966 (0.9946)

0.9995 (0.9994)

0.9998 (0.9997)

0.9993 (0.9994)

0.9995 (0.9996)

Note. Numbers in the parentheses are the results from Fig. 7. a Calculations based on 20-point-long segments surrounding seed points prior to entering the PVA region. b Calculations based on entire ‘‘fiber’’ tracts.

1070

Y. Lu et al. / NeuroImage 31 (2006) 1061 – 1074

PVA regions with the Bayesian and TEND methods but not with the Euler method; by referring to the ‘‘true tracts’’ (black lines), it can be seen that ‘‘fibers’’ tracked with the Bayesian method are more faithful than those of the Euler and TEND methods. Table 1 summarizes measurements of tracking error, noise sensitivity, and tract smoothness from the three methods. As expected, the Bayesian method has smaller tracking errors and noise sensitivity, and greater tract smoothness than the Euler method. Compared to the TEND method, the Bayesian method generally has smaller tracking errors (except for the case of PVA in straight ‘‘fibers’’ at SNR = 20) and noise sensitivity, although the tract smoothness is comparable.

method to noise and PVA. Second, it can be appreciated with close inspection that, due to its noise immunity, the tracts from the Bayesian method appear to be smoother. Fig. 9 shows both views of the fiber tracts tracked with the same seed points as in Fig. 8 but from the dataset generated with added noise. Again, fibers from the Bayesian method appear more favorable. Compared with Fig. 8, noise seems to induce minor changes to fiber tracts from the Bayesian method but adds a few additional erroneous connections (pointed to by bright green arrows in panel c to fiber tracks from the Euler method. Quantitative analysis of tracking noise sensitivity from five pairs of randomly selected fibers shows that the Euler method is more sensitive to noise, as summarized in Table 2.

Fiber tracking with in vivo human data Discussion and conclusions Fig. 8 shows the axial (top row) and sagittal (bottom row) views of fiber tracts seeded in the superior longitudinal fasciculus tracked with the Euler (left column) and Bayesian (right column) methods. Although the major patterns are grossly similar from both methods, the results from the Bayesian method are more plausible. First, the Euler method produces a few erroneous connections (pointed to by green arrows in panel (c)) that are absent from the Bayesian method, which may be explained by the sensitivity of the Euler

To meet the pressing need for fiber tracking algorithms that are robust to noise and partial volume averaging (PVA), a novel Bayesian fiber tracking method was proposed and investigated. In this study, the framework of this algorithm was established, and its performance with respect to noise and PVA was comprehensively examined. Experiments with Monte Carlo simulations, synthetic and in vivo human diffusion tensor data demonstrate this new

Fig. 8. Axial (a, b) and sagittal (c, d) views of reconstructed fiber tracts seeded in the superior longitudinal fasciculus from the original high quality data (SNR > 100). Left column (a, c) shows fibers from the Euler method (red curves), and right column (b, d) shows fibers from the Bayesian method (blue curves). Green arrows point to erroneous pathways generated by the Euler method.

Y. Lu et al. / NeuroImage 31 (2006) 1061 – 1074

1071

Fig. 9. Axial (a, b) and sagittal (c, d) views of reconstructed fiber tracts seeded in the superior longitudinal fasciculus from the noisy dataset generated. Left column (a, c) shows fibers from the Euler method (red curves), and right column (b, d) shows fibers from the Bayesian method (blue curves). Bright green arrows point to additional erroneous pathways induced by noise.

algorithm has superior immunity to noise and PVA. Compared to the Euler as well as the TEND method, the Bayesian method proposed has significantly improved accuracy and precision. Noise and PVA are major hurdles in DTI based fiber tractography (Alexander et al., 2001; Basser and Pajevic, 2000). This is not only because both are inherent artifacts, due to random signal perturbations and finite imaging resolution, but also because they interact with each other. As alluded to earlier, noise deflects the direction of the tensor major eigenvector, and PVA creates directional ambiguities. With these artifacts, the major eigenvector always comes with a certain degree of uncertainty. Typically, the noise induced uncertainty increases with decreased diffusion anisotropy. Most notably, in regions where PVA has significantly

Table 2 Comparison of tracking noise sensitivity (expressed in voxels) between the Euler and Bayesian methods with in vivo human data

Fiber pair Fiber pair Fiber pair Fiber pair Fiber pair Average

1 2 3 4 5

Euler

Bayesian

1.08 1.05 0.20 0.25 0.10 0.53

0.38 0.76 0.12 0.06 0.24 0.31

reduced the diffusion anisotropy, the uncertainty caused by noise is particularly pronounced. The uncertainty induced by noise and exacerbated by PVA brings special complications to diffusion tensor tractography since errors may be accumulated throughout the process of line propagation. Therefore, the robustness of fiber tracking algorithms to these interwoven artifacts is essential to the reliability of reconstructed fiber tracts. Although there exist many noise reduction methods for regularizing DTI data (Parker et al., 2000; Pajevic et al., 2002; Wang et al., 2003; Coulon et al., 2004; Ding et al., 2005), and most fiber tracking algorithms have indeed considered noise inherently (see Introduction), to date only a limited few have built in the functionality of correcting for the effect of both noise and PVA, particularly with a unified framework. A representative example of such unified functionality is the tensor deflection model by Lazar et al. (2003), but the notion of tensor deflection remains empirical. To achieve a similar goal, Tuch (2002) introduced the traditional Bayes decision rule into fiber tracking by formulating the tracking problem as that of optimal path determination. In his pioneering work, the a priori probability was modeled with a fiber stiffness function and the conditional probability modeled with a fiber orientational distribution function measured from diffusion spectrum imaging (Tuch, 2002); an optimal fiber path was obtained by a numerical optimization procedure. The Bayesian fiber tracking proposed herein employs a similar concept of fiber path

1072

Y. Lu et al. / NeuroImage 31 (2006) 1061 – 1074

optimization but has novel formulations of probability functions that are tailored specifically for DTI. Compared to the work by Tuch, and others that use the Bayesian formulism for fiber path optimization (e.g., Poupon et al., 2000; Behrens et al., 2003; Hosey et al., 2005), this study possesses the following salient features. First, the probability models are predicated on multivariate normal distributions of diffusion tensors. In contrast to the simplifying univariate exponential models of major eigenvector or path direction, the multivariate model for the conditional probability may better approximate local complex fiber orientation distributions (Pajevic and Basser, 2003). Second, as proven previously, modeling of two probability functions based on normal distributions naturally leads to an exact analytical solution to the path optimization problem. Compared to numerical optimizations (Poupon et al., 2000; Tench et al., 2002), this gives considerable computational savings and eliminates the problem of convergence. Third, the tracking method proposed does not require any ad hoc parameters, which makes it operationally simple. The versatility and convenience of these features allow this method to be used as a standard routine utility. It should be pointed out that, according to the framework designed, the mechanism for noise immunity is quite obvious, but for PVA is subtler. The noise immunity comes from the fact that noise induced deflection of the major eigenvector can be corrected for by the a priori probability that draws upon parameters estimated more reliably from a volume of neighboring samples; this sampling volume is oriented with the direction of fiber propagation and is further weighted by FA to suppress PVA confounds. By the same token, the ambiguity due to PVA may also be resolved by the directional ‘‘inertia’’ provided by the oriented sampling volume (which favors straight fiber paths in regions of ambiguity). It should be noted, in worst case scenarios (e.g., fiber crossing), the PVA immunity is also partly provided by the wider spreading of the conditional probability density. Admittedly, this mechanism is neither very specific nor sensitive. The limitation stems from the inability of tensor model to resolve fiber crossing and also the insufficiency of a (unimodal) multivariate normal distribution to model such a situation. Remedies for this problem rest on multifiber modeling (Zeng et al., 2005; Hosey et al., 2005) or more fundamentally the use of high angular resolution diffusion imaging (Anderson, 2005; Tuch, 2004) which promises an ultimate solution (in this case, the proposed Bayesian tracking method can be similarly applied). Nonetheless, given the fact that tensor model remains the de facto standard in current clinical settings, and that multifiber modeling may incur other complications, the framework proposed represents a viable practical compromise. A number of questions might be raised regarding the estimation of the a priori probability. First, the use of oriented sampling volume is heuristic, and, as a result, the Bayesian decision is also heuristic. Second, the assumption that tensors within this sampling volume are normally distributed may be inadequate. Although it is valid to assume a normal distribution for tensors arising from single fiber populations (Pajevic and Basser, 2003), it may be not for tensors in regions with more complex fiber architectures. To examine this, normality tests on in vivo human data were carried out along different fiber tracts in the corpus callosum (CC) and superior longitudinal fasciculus (SLF), the common testbeds for fiber tracking algorithms. Specifically, the Lilliefors normal distribution goodness-of-fit testing (Conover, 1980) was used to test the normality of the tensors within the oriented, and FA

Table 3 Normal distribution rate of tensor elements in the corpus callosum (CC) and superior longitudinal fasciculus (SLF) (%)

D 11 (%)

CC a=2 a=5 a = 10

88 94 100

82 82 88

SLF a=2 a=5 a = 10

83 94 94

94 100 100

D 22 (%)

D 33 (%)

D 12 (%)

D 13 (%)

D 23 (%)

65 71 82

82 94 100

76 76 82

82 94 100

83 94 100

83 94 94

72 83 89

83 83 83

weighted sampling volume for each point along the fiber, and the rate of conformity to normal distribution to the total number of points tested was calculated. Normal distribution rates for the CC and SLF at different significance levels are summarized in Table 3. It shows that at significance level a = 5%, the vast majority of points (¨88%) have a normal tensor distribution in their sampling volume; the rate is higher for significance level a = 10%, reaching an overall rate of ¨93%. In spite of the imperfect normal distribution of tensors sampled heuristically, experiments on the synthetic data that may not satisfy the normality assumption demonstrated the desired functionality could still be achieved; this indicates a violation of normality may not imply a significant performance degradation of the tracking algorithm. Therefore, the normality assumption should be reasonable, particularly in light of the computational simplicity it offers. It has been previously mentioned that regularized fiber tracking tends to have a ‘‘straightening effect’’, i.e., a bias toward straight fibers. This has been indeed observed in Fig. 5(a) which shows the tendency of ‘‘fiber’’ tracts toward the low curvature side. Interestingly, the Euler method as implemented also has the straightening effect, but as pointed out, this is due to the large step size used. It is conceivable that a smaller step size (and higher imaging spatial resolution) will allow the Euler method to follow the true fiber path more faithfully (therefore, the step size here can be thought of as a regularizing factor). In a similar vein, if the Bayesian method proposed regularizes the tracking direction with a smaller sampling volume for the prior probability, a better tract sensitivity can also be gained, but in the meantime, the regularization power against the effect of noise and particularly PVA will be decreased. Hence, there comes a trade-off between the robustness to noise/PVA and tract sensitivity. As comprehensive exploration of this trade-off is beyond the scope of this study, it is reserved for future research. Another limitation of a more practical nature is the time efficiency of this method, since it involves the estimation of mean/ covariance matrices and several matrix inverse calculations. For example, on an Intel Pentium M 1.4-GHz computer, tracking a 64mm-long fiber with a step size of 1 mm, the Bayesian method takes 1.3 s, but the Euler method only 0.2 s. As fiber tracking is typically performed in an off-line fashion, and the speed of desktop computers is rapidly evolving, the issue of time efficiency should not be a major concern. In summary, in this study, a novel fiber tracking framework robust to the effects of noise and PVA was developed. The framework is based on a new formulation of Bayes decision rule and possesses a number of useful features from the new formulation. Notwithstanding the limitations mentioned above,

Y. Lu et al. / NeuroImage 31 (2006) 1061 – 1074

experiments using Monte Carlo simulations, with synthetic and in vivo diffusion tensor data have demonstrated the tracking method developed can significantly improve the accuracy and precision of fiber tracking. This offers the potential of using it as a reliable routine tool for DTI-based fiber tractography.

Acknowledgments The authors are very grateful to Rui Li (Department of Electrical Engineering and Computer Science, Vanderbilt University) for generously providing the 3D visualization tool. Thanks are also given to Drs. Jane Park and Jing Qi (Department of Molecular Physics and Biochemistry, Vanderbilt University) for reviewing the manuscript and offering valuable comments. This work was supported by NIH grants R01EB000461 funded to John C. Gore and R01EB02777 to Adam W. Anderson.

Appendix A. Definitions of tracking error, noise sensitivity and tract smoothness A.1. Tracking error Tracking error is defined to be the mean Euclidean distance between the corresponding segments of a putative fiber pathway and its true pathway (see (Ding et al., 2003) for the definition of corresponding segment). Let P = ( p 1, p 2, . . ., p N ) be a segment of a putative pathway, where p i is the coordinate of point i in Euclidean space, N is the total number of points, and P 0 = ( p 10, p 20, . . ., p n0) be the corresponding segment of the true fiber pathway, the tracking error is: e¼

N 1 X pi  p0i : N i¼ 1

ð18Þ

A large e indicates a high tracking error. A.2. Tracking noise sensitivity Tracking noise sensitivity is defined to be the mean Euclidean distance between the corresponding segments of two putative pathways tracked from the same seed point but on data at two different noise levels. Let P 1 = ( p 11, p 21, . . ., p N1) be a segment of the putative pathway at noise level 1, and P 2 = ( p 12, p 22, . . ., p N2) be the corresponding segment of the putative pathway at noise level 2, the tracking noise sensitivity is: N 1 X n¼ p1i  p2i : ð19Þ N i¼1 A large n indicates a high tracking noise sensitivity. A.3. Tract smoothness Tract smoothness is defined as: s¼

1 N 1 pi þ 1  pi  1  ; ~ N  2 i ¼ 2 pi þ 1  pi  þ pi  pi  1 

ð20Þ

where p i is the coordinate of point i at a putative pathway, and N is the total number of points. The range of s is between 0 and 1, with 1 representing the smoothest (a straight line).

1073

References Alexander, A.L., Hasan, K.M., Lazar, M., Tsuruda, J.S., Parker, D.L., 2001. Analysis of partial volume effects in diffusion-tensor MRI. Magn. Reson. Med. 45, 770 – 780. Anderson, A.W., 2001. Theoretical analysis of the effects of noise on diffusion tensor imaging. Magn. Reson. Med. 46, 1174 – 1188. Anderson, A.W., 2005. Measurement of fiber orientation distributions using high angular resolution diffusion imaging. Magn. Reson. Med. 54, 1194 – 1206. Basser, P.J., Pajevic, S., 2000. Statistical artifacts in diffusion tensor MRI (DT-MRI) caused by background noise. Magn. Reson. Med. 44, 41 – 50. Basser, P.J., Mattiello, J., Le Bihan, D., 1994. Estimation of the effective self-diffusion tensor from the NMR spin echo. J. Magn. Reson. 103, 247 – 254. Basser, P.J., Pajevic, S., Pierpaoli, C., Duda, J., Aldroubi, A., 2000. In vivo fiber tractography using DT-MRI data. Magn. Reson. Med. 44, 625 – 632. Behrens, T.E., Woolrich, M.W., Jenkinson, M., Johansen-Berg, H., Nunes, R.G., Clare, S., Matthews, P.M., Brady, J.M., Smith, S.M., 2003. Characterization and propagation of uncertainty in diffusion-weighted MR imaging. Magn. Reson. Med. 44, 1077 – 1088. Conover, W.J., 1980. Practical Nonparametric Statistics. Wiley, pp. 44 – 45. Coulon, O., Alexander, D.C., Arridge, S., 2004. Diffusion tensor magnetic resonance image regularization. Med. Image Anal. 8, 47 – 67. Ding, Z., Anderson, A.W., Gore, J.C., 2003. Classification and quantification of neuronal fiber pathways using diffusion tensor MRI. Magn. Reson. Med. 49, 716 – 772. Ding, Z., Anderson, A.W., Gore, J.C., 2005. Reduction of noise in diffusion tensor images using anisotropic smoothing. Magn. Reson. Med. 49, 485 – 490. Duda, R.O., Hart, P.E., 1973. Pattern Classification and Scene Analysis. John Wiley and Sons, New York, pp. 10 – 15. Go¨ssl, C., Fahrmeir, L., Pu¨tz, B., Auer, L.M., Auer, D.P., 2002. Fiber tracking from DTI using linear state space models: detectability of the pyramidal tract. NeuroImage 16, 378 – 388. Hagmann, P., Thiran, J.P., Jonasson, L., Vandergheynst, P., Clarke, S., Maeder, P., Meuli, R., 2003. DTI mapping of human brain connectivity: statistical fiber tracking and virtual dissection. NeuroImage 19, 545 – 554. Hosey, T.P., Harding, S.G., Green, H.A., Ansorge, R.E., Carpenter, T.A., Williams, G.B., 2005. Probabilistic determination and model based regularization of intra-voxel multiple fibre architecture in high angular resolution diffusion imaging data sets. Proceedings of 13th ISMRM, Miami, Florida, pp. 221. Koch, M.A., Norris, D.G., Hund-Georgiadis, M., 2002. An investigation of functional and anatomical connectivity using magnetic resonance imaging. NeuroImage 16, 241 – 250. 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., 2003. White matter tractography using diffusion tensor deflection. Hum. Brain Mapp. 18, 306 – 321. Le Bihan, D., Mangin, J.F., Poupon, C., Clark, C.A., Pappata, S., Molko, N., Chabriat, H., 2001. Diffusion tensor imaging: concepts and applications. J. Magn. Reson. Imaging 13, 534 – 546. Li, R., 2001. Automatic placement of regions of interest in medical images using image registration. Master thesis. Vanderbilt University. Mattiello, J., Basser, P.J., LeBihan, D., 1994. Analytical expressions for the b matrix in NMR diffusion imaging and spectroscopy. J. Magn. Reson. 108, 131 – 141. Mori, S., van Zijl, P.C.M., 2002. Fiber tracking: principles and strategies— A technical review. NMR Biomed. 15, 468 – 480. Mori, S., Crain, B.J., Chacko, V.P., van Zijl, P.C.M., 1999. Three dimensional tracking of axonal projections in the brain by magnetic resonance imaging. Ann. Neurol. 45, 265 – 269. Pajevic, S., Basser, P.J., 2003. Parametric and non-parametric statistical analysis of DT-MRI data. J. Magn. Reson. 161, 1 – 14.

1074

Y. Lu et al. / NeuroImage 31 (2006) 1061 – 1074

Pajevic, S., Aldroubi, A., Basser, P.J., 2002. A continuous tensor field approximation of discrete DT-MRI data for extracting microstructural and architectural features of tissue. J. Magn. Reson. 154, 85 – 100. Parker, G.J.M., Schnabel, J.A., Symms, M.R., Werring, D.J., Barker, G.J., 2000. Nonlinear smoothing for reduction of systematic and random errors in diffusion tensor imaging. J. Magn. Reson. Imaging 11, 702 – 710. Parker, G.J.M., Wheeler-Kingshott, C.A.M., Barker, G.J., 2002. Estimating distributed anatomical connectivity using fast marching methods and diffusion tensor imaging. IEEE Trans. Med. Imaging 21, 505 – 512. Parker, G.J.M., Haroon, H.A., Wheeler-Kingshott, C.A.M., 2003. A framework for a streamline-based probabilistic index of connectivity (PICo) using a structural interpretation of MRI diffusion measurements. J. Magn. Reson. Imaging 18, 242 – 254. Poupon, C., Clark, C.A., Frouin, V., Re’gis, J., Bloch, I., Le Bihan, D., Mangin, J.F., 2000. Regularization of diffusion-based direction maps for the tracking of brain white matter fascicles. NeuroImage 12, 184 – 195. Staempfli, P., Jaermann, T., Crelier, G.R., Meier, D., Valavanis, A., Kollias, S., Boesiger, P., 2005. Advanced fast marching tractography based on

SENSE-DTI: an attempt to resolve fiber crossing in artificial and invivo data. Proceedings of 13th ISMRM, Miami, Florida, pp. 392. Tench, C.R., Morgan, P.S., Blumhardt, L.D., Constantinescu, C., 2002. Improved white matter fiber tracking using stochastic labeling. Magn. Reson. Med. 48, 677 – 683. D.S., Tuch, 2002. Diffusion MRI of complex tissue structure Doctoral dissertation. Harvard University – Massachusetts Institute of Technology. Tuch, D.S., 2004. Q-ball imaging. Magn. Reson. Med. 52, 1358 – 1372. Wang, Z., Vemuri, B.C., Chen, Y., Mareci, T., 2003. A constrained variational principle for direct estimation and smoothing of the diffusion tensor field from DWI. In: Taylor, C.J., Noble, J.A. (Eds.), Information Processing in Medical Imaging. Springer-Verlag, Heidelberg, pp. 660 – 671. Zeng, Q., Guo, W., Chen, Y., He, A.G., Liu, Y., 2005. White matter fiber tracking based on multi-directional vector field. Proceedings of 13th ISMRM, Miami, Florida, pp. 218. Zhukov, L., Barr, A.H., 2002. Orientated tensor reconstruction: tracing neural pathways from diffusion tensor MRI. Proceedings IEEE Visualization.