# Fiber Direction Estimation, Smoothing and Tracking in Diffusion MRI

## Abstract

Diffusion magnetic resonance imaging is an imaging technology designed to probe anatomical architectures of biological samples in an in vivo and non-invasive manner through measuring water diffusion. The contribution of this paper is threefold. First it proposes a new method to identify and estimate multiple diffusion directions within a voxel through a new and identifiable parametrization of the widely used multi-tensor model. Unlike many existing methods, this method focuses on the estimation of diffusion directions rather than the diffusion tensors. Second, this paper proposes a novel direction smoothing method which greatly improves direction estimation in regions with crossing fibers. This smoothing method is shown to have excellent theoretical and empirical properties. Lastly, this paper develops a fiber tracking algorithm that can handle multiple directions within a voxel. The overall methodology is illustrated with simulated data and a data set collected for the study of Alzheimer’s disease by the Alzheimer’s Disease Neuroimaging Initiative (ADNI).

## 1Introduction

Diffusion magnetic resonance imaging (dMRI) is an in vivo and non-invasive medical imaging technology that uses water diffusion as a proxy to probe anatomical structures of biological samples. The most important application of dMRI is to reconstruct white matter fiber tracts in brain – large axonal bundles with similar destinations. In white matter, water diffusion appears to be anisotropic as water tends to diffuse faster along the fiber bundles. Therefore, white matter fiber structures can be deduced from the diffusion characteristics of water. Mapping white matter fiber tracts is of great importance in the study of structural organization of neuronal networks and the understanding of brain functionality [?]. Moreover, dMRI also has many clinical applications, including detecting brain abnormality in white matter due to axonal loss or deformation, which are thought to be related to many neuron degenerative diseases including Alzheimer’s disease, and also in surgical planning by resolving complex neuronal connections between white and gray matter [?].

dMRI techniques sensitize signal intensity with the amount of water diffusion by applying pulsed magnetic gradient fields on the sample. Specifically, water diffusion along the gradient field direction leads to signal loss and the amount of loss at a voxel equals to the summation (across locations within the voxel) of the sinusoid waves with shifted signal phases weighted by the proton density at their respective locations. In other words, signal loss (referred to as diffusion weighted signal) is the inverse Fourier transform of the diffusion probability density function of water molecules and thus can be used to recover water diffusion characteristics. The amount of signal loss is also influenced by various experimental parameters including the gradient field intensity (the stronger, the more loss), the duration of gradient fields (the longer, the more loss), etc. Their effects are aggregatively reflected by an experimental parameter called the “-value” which is often fixed throughout the experiment (though multiple -values are used in Q-space imaging). Since only water motion along the gradient field direction can be detected, multiple gradient directions need to be applied [?].

In its raw form, dMRI provides diffusion weighted signal measurements on a 3D spatial grid (of the sample) along a set of predetermined gradient directions [?]. For example, a typical data set from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) has diffusion measurements along gradient directions for each voxel on a 3D grid of the brain. The first step of dMRI analysis is to summarize these measurements into estimates of water diffusion at each voxel. A popular model for water diffusion is the so called single tensor model where the diffusion process is modeled as a 3D Gaussian process described by a positive definite matrix, referred to as a diffusion tensor; see [?] for an introduction to diffusion tensor imaging (DTI) techniques. Figure 1 depicts a tensor map on a 2D grid, where each diffusion tensor is represented by an ellipsoid, estimated from diffusion weighted measurements from an ADNI data set using a single tensor model. One then extracts the local diffusion direction as the principal eigenvector of the (estimated) diffusion tensor at each voxel and reconstructs the white matter fiber tracts by computer aided tracking algorithms via a process called tractography [?].

However, DTI cannot resolve multiple fiber populations with distinct orientations, i.e., crossing fibers, within a voxel since a tensor only has one principal direction. Consequently, in crossing fiber regions, estimated diffusion tensors may lead to low anisotropy estimation or oblate tensor estimation. Poor tensor estimation results in poor direction estimation which adversely affects fiber reconstruction; e.g., early termination of or biased fiber tracking.

In order to resolve intravoxel orientational heterogeneity, several approaches have been proposed. [?] propose a multi-tensor model which assumes a finite number of homogeneous fiber directions within a voxel. However, it has been shown that the parameters in the multi-tensor model are not identifiable [?]. Imaging techniques such as Q-ball and Q-space and the corresponding nonparametric methods have also been proposed [?]. However such methods often require high angular resolution diffusion imaging (HARDI) [?] where a large number of gradients is sampled. In light of these facts, the goal of this paper is to develop a new fiber direction estimation and tracking method that can handle crossing fibers without requiring any high resolution techniques. The proposed method, named DiST, short for **Di**ffusion Direction **S**moothing and **T**racking, is completely automated and improves existing methods in several aspects. Particularly, it is applicable either when there is a large number of gradient directions (as in the HARDI setting) or when only a relatively small number of gradient directions are available (as in most clinical settings).

The DiST method can be divided into three major steps.

**Step 1:** Estimate the tensor directions within each voxel under a multi-tensor model. A new parametrization is proposed which makes the tensor directions identifiable. An efficient and numerically stable computational procedure is developed to obtain the maximum likelihood (ML) estimate of the tensor directions. Here we highlight that, this method focuses on the estimation of the tensor directions rather than the actual tensors themselves.

**Step 2:** Using the voxel-wise tensor direction estimates from Step 1 as input, a new direction smoothing procedure is applied to further improve the diffusion direction estimates by borrowing information from neighboring voxels. A distinctive and unique feature of this smoothing procedure is that it handles crossing fibers through the clustering of directions into homogeneous groups. We note that, although various tensor smoothing methods have been proposed [?], little work has been done on direct diffusion direction smoothing. One notable exception is [?], which harnesses diffusion directions directly to construct a map of test statistics for detecting differences between diffusion direction maps from two groups of subjects,while the spatial smoothness of the test statistics is being considered. Also note that approaches to averaging unsigned directions in the real projective space are known in the directional statistics literature.

**Step 3:** Lastly, a fiber tracking algorithm is applied to reconstruct fiber tracts using the smoothed diffusion direction estimates obtained in Step 2. This tracking algorithm is designed to explicitly allow for multiple directions within a voxel.

We apply DiST to an ADNI data set measured on a healthy elderly person with 41-direction dMRI scan on a 3 Tesla GE Medical Systems MRI scanner. ADNI is a longitudinal study (since 2005) that collects serial MRI, cognitive assessments, and numerous additional measurements approximately twice per year from hundreds of elderly individuals spanning a range from cognitive health to clinically-diagnosed Alzheimer’s disease. We also examine DiST using simulated data sets which mimic the most commonly encountered experimental situations in terms of number of gradient directions and signal to noise ratio. DiST is shown to lead to superior results than those based on the single tensor model in the simulation study, as well as more biologically sensible results in the real data application.

The rest of the paper is organized as follows. Section 2 provides background material for some common tensor models. The proposed methods for tensor direction estimation, smoothing of estimated directions, and fiber tracking are presented in, respectively, Sections Section 3, Section 4 and Section 5. Section 6 summarizes simulation results. The application to an ADNI data set is presented in Section 7. Section 8 provides some concluding remarks, while additional simulation results and technical details are collected in an online Supplemental Material.

## 2Tensor models

Suppose dMRI measurements are made on voxels on a 3D grid representing a brain. For each voxel, we have measurements of diffusion weighted signals along a fixed set (i.e., the same for all voxels) of unit-norm gradient vectors . We write the set of measurements as , where is the 3D coordinate of the center of this voxel.

Assuming Gaussian diffusion, the noiseless signal intensity is given by [?]

where is the non-diffusion-weighted intensity, is an experimental constant referred to as the -value and is a covariance matrix referred to as the diffusion tensor. This model is called the single tensor model and suits for the case of at most one dominant diffusion direction within a voxel.

Although the single tensor model is the most widely used tensor model in practice, it is not suitable for crossing fiber regions. To deal with crossing fibers, this model has been extended to a multi-tensor model [?]:

where and for . Here represents the number of fiber populations and ’s denote weights of the corresponding fibers.

## 3Voxel-wise estimation of diffusion directions

One important goal of dMRI studies is to estimate principal diffusion directions, referred to as diffusion directions hereafter, at each voxel. They may be interpreted as tangent directions along fiber bundles at the corresponding voxel. The estimated diffusion directions are then used as an input for tractography algorithms to reconstruct fiber tracts. This section explores the diffusion direction estimation within a single voxel. For notational simplicity, dependence on voxel index is temporarily dropped. Moreover, for ease of exposition, we assume that and are known

and delay the discussion of their estimation to Section 7.

Under the single tensor model, various methods for tensor estimation have been proposed including linear regression, nonlinear regression and ML estimation; e.g., see [?] for a comprehensive review. Then diffusion directions are derived as principal eigenvectors of (estimated) diffusion tensors. However, for the estimation of multi-tensor models, severe computational issues have been observed and additional prior information and additional assumptions are usually imposed to tackle these issues. For instance, [?] use shrinkage priors and [?] assume all tensors to be axially symmetric (i.e., the two minor eigenvalues are the same) and have the same set of eigenvalues. [?] show that the multi-tensor model is indeed non-identifiable in the sense that there exist multiple parameterizations that are observationally equivalent. These authors suggest to use multiple -values in data acquisition to make the model identifiable. However, due to practical limitations, most of the current dMRI studies are obtained under a fixed -value and so render their suggestion inapplicable. Below we show that the identifiability issue does not prevent one from estimating the diffusion directions and so neither strong assumptions nor special experimental settings are necessary if one is only interested in diffusion directions rather than the diffusion tensors themselves.

### 3.1Identifiability of multi-tensor model

Model (Equation 1) can be re-written as

where for such that , is positive definite and . When , one can easily derive the explicit conditions for to fulfill these criteria, and see that there are infinite sets of such ’s. However, note that shares the same set of eigenvectors with . Thus, one may still be able to estimate diffusion directions, which correspond to the major eigenvectors of the tensors. This motivates us to consider estimating diffusion directions directly instead of the tensors themselves.

Now we assume that ’s are axially symmetric; that is, the two minor eigenvalues of are equal. This is a common assumption for modeling dMRI data and it implies that diffusion is symmetric around the principal diffusion direction [?]. By not differentiating the two minor eigenvectors, we obtain a clear meaning of diffusion direction. In addition, this reduces the number of unknown parameters by one for each tensor in the multiple tensor model and thus facilitates estimation. In the following, we propose a new parametrization of the multi-tensor model which is identifiable and thus can be used for direction estimation.

Write as the space of the unit principal eigenvector; i.e., the 3D unit sphere with equivalence relation . Let , and be the difference between the larger and smaller eigenvalue, smaller eigenvalue and the standardized principal eigenvector of , respectively. Since , model (Equation 1) becomes

where . From the above, one can see that and are not simultaneously identifiable, so we cannot estimate the tensors. However, as stated in the following theorem, are identifiable and hence we can estimate the principal diffusion directions ’s.

The proof of this theorem can be found inSection S5.1 of the Supplemental Material. Note that,compared to the model in [?],model () allows for different eigenvalues and shapes of the tensors within a voxel,and thus is much more flexible.

### 3.2Parameter estimation using maximum likelihood (ML)

We first consider the case when is known and delay the selection of to Section 3.3. By assuming Gaussian additive noise on both real and imaginary parts of thecomplex signal, the observed signal intensity can be modeled as[?] where is the intensity of the noiseless signal, is a unit vector in representing the phase of the signal, is the noise random variable following and denotes the noise level. Note that both and may depend on . The observed signal intensity then follows a Rician distribution [?]: Moreover, we assume the noise ’s are independent across different voxels and gradient directions.

Under the Rician noise assumption, the log-likelihood of in model (Equation 2) is:

where is the zeroth order modified Bessel function of the first kind. The ML estimate is obtained through maximizing (Equation 3). Although the above new parametrization avoids the identifiability issue, the likelihood function usually has multiple local maxima, which makes the computation of ML estimate difficult and unstable.

The method that we used to overcome this issue can be briefly described as follows. We first develop an approximation of model (Equation 2) whose likelihood can be globally maximized via a grid search. We utilize the geometry of the problem so that the grid search can be done efficiently. Then we use the ML estimate of this approximated model as the initial value in a gradient method to obtain the ML estimate of model (Equation 2). This method provides very reliable estimates. To speed up the pace of this article, its full description is given in Section S1 of the Supplemental Material.

### 3.3Selecting the number of tensor components

Common model selection methods can be applied to select the number of components . Results from extensive numerical experiments suggest that the Bayesian information criterion (BIC) [?] is good choice; see Section S2 of the Supplemental Material.

Under model (), eachtensor corresponds to four free scalar parameters since is characterized by two free scalar parameters. The BIC for a model with tensors is

where is the number of gradient directions and is the ML estimate of under tensors. Then is chosen as where is a pre-specified upper bound for the number of components. Based on our experience, is a reasonable choice.

In practice, there are voxels with no major diffusion directions. This corresponds to the case where there is only one isotropic tensor. In the case of isotropic tensor, (Equation 2) reduces to Thus there is only one parameter . We write the corresponding likelihood function as and denote the ML estimate of by , which can be obtained by a generic gradient method. The corresponding BIC criterion is

where 0 represents no diffusion direction. Combined with the previous BIC formulation (Equation 4), one has a comprehensive model selection rule, which handles voxels with from zero to up to (here 4) fiber populations.

In practice, we follow the convention and use fractional anisotropy (FA) [?], where , and are the eigenvalues of the corresponding tensor in the single tensor model,to conduct an initial screening to speed up the whole procedure. The FA value lies between zero and one and the larger it is, the more anisotropic the water diffusion is at the corresponding voxel. Thus we first remove voxels with very small FA values and then apply the BIC approach over those suspected anisotropic voxels. Note that such removal is mainly for reducing computational cost as a typical dMRI data set consists of hundreds of thousands of voxels. From our experiences, this has little effect on the final tracking results. We also note that the proposed framework including selection of can be applied without such removal if enough computational resources are available.

We summarize our voxel-wise estimation procedure in Algorithm S2 in the Supplemental Material. A simulation study is conducted and the corresponding results are presented in Section S2 of the Supplemental Material. These numerical results suggest that our voxel-wise estimation procedure provides extremely stable and reliable results under various settings.

## 4Spatial smoothing of diffusion directions

Although model (Equation 2) provides a better modeling than the single tensor model for crossing fiber regions, it also leads to an increase in the number of parameters and thus the variability of the estimates. To further improve estimation, we consider borrowing information from neighboring voxels and develop a novel smoothing technique for diffusion directions.

In many brain regions, it is reasonable to model the fiber tracts as smooth curves at the resolution of voxels in dMRI ( 2mm). Therefore, we shall assume that the tangent directions of fiber bundles change smoothly. This leads to the spatial smoothness of diffusion directions that belong to the same fiber bundle.

### 4.1Smoothing along a single fiber

This subsection considers the simpler situation where there is only one homogeneous population of diffusion directions; i.e., there is only one single fiber bundle without crossing. Write as the total number of estimated diffusion directions from all voxels and as the set of all estimated diffusion directions. Also write as the corresponding voxel location associated with . Note that some ’s share the same value, as some voxels contain multiple estimated directions. Following the idea of kernel smoothing on Euclidean space [?], the smoothing estimate at voxel is defined as a weighted Karcher mean of the neighboring direction vectors:

where ’s are spatial weights and the metric is defined as

i.e., is the acute angle between and . The weights ’s place more emphasis on spatially closer observations. Here with as a 3D kernel function satisfying , and is a bandwidth matrix. In our numerical work, we choose as the standard Gaussian density, and set , where is chosen using the cross-validation (CV) approach described in Section S3 of the Supplemental Material. We adopt the leave-one-out CV idea to develop an ordinary CV score and two robust CV scores. Their practical performances are reported inSection S6 of the Supplemental Material.

### 4.2Smoothing over multiple fibers

When there are crossing fibers in a voxel , the above smoothing procedure will not work well. To address this issue, we first cluster the neighboring estimated directions of into groups that correspond to different fiber populations. Then we apply the above smoothing procedure to each individual cluster. This subsection describes this procedure in details.

First we define neighboring voxels for . We begin with computing the spatial weights defined in Section 4.1. We then remove those voxels with weights smaller than a threshold. By filtering out these voxels, we obtain tighter and better separated clusters of directions. Moreover, such voxels have little effects on smoothing due to their small weights. The artificial data set displayed in Figure ? provides an illustrative example. Each black dot in the left panel represents an estimated direction (from the center of the sphere). In the middle panel, the size of each dot is proportional to its spatial weight in equation (Equation 5). Lastly, the right panel shows all dots with spatial weights larger than a threshold. Notice that such a trimming operation leads to two obvious clusters of directions, which makes the subsequent task of clustering the directions much easier.

Next we need a clustering strategy to choose the number of clusters adaptively. With the distance metric (Equation 6), one can define dissimilarity matrix for a set of directions and make use of a generic clustering algorithm. Our choice is the Partition Around Medoids (PAM) [?] due to its simplicity. Also, we apply the average silhouette [?] to choose the number of clusters; see Algorithm S3 of the Supplemental Material. The silhouette of a datum measures the strength of its membership to its cluster, as compared to the neighboring cluster. Here, the neighboring cluster is the one, apart from cluster of datum , that has the smallest average dissimilarity with datum . The corresponding silhouette is defined as , where and represent the average dissimilarities of datum with all other data in the same cluster and that with the neighboring cluster respectively. The average silhouette of all data gives a measure of how good the clustering is. Thus we select the number of clusters via maximizing the average silhouette. The detailed smoothing procedure is given in Algorithm S4.

### 4.3Theoretical results

This subsection derives asymptotic properties of the proposed direction smoothing estimator. Note that, since the space of direction vectors has a non-Euclidean geometry and so the theoretical framework is different from that of classical smoothing estimators. Without loss of generality, suppose we observe at spatial locations respectively. Let be the 3D unit sphere. Then is the quotient space of with equivalence relation for any . This space is also identified with the so-called real projective space .

The theoretical results below were derived under the more convenient random design where ’s are independently and identically sampled from a distribution with density .

The below theorem (Theorem ?) remains valid even under a fixed, regular design setting, with the number of grid points increasing to infinity. In this case, in the statement of the asymptotic formulae and their proofs, the density function is replaced with a constant-valued function, representing a regular grid, with corresponding changes wherever derivatives of appear.

Given a spatial location , our target is to estimate , namely the diffusion direction at , which is defined as the minimizer of over , where . Here is a random unit vector representing a random diffusion direction and the expectation is taken over conditional on , where represents the location of where is observed. For simplicity, we assume and write it as thereafter. Thus, our estimator (Equation 5) at can be written as

where is the number of diffusion direction vectors and . Here, with slight notation abuse, represents a one dimensional kernel function throughout the theoretical developments.

We first describe a working coordinate system. For each , one can endow a tangent space with the metric tensor defined as . Note that the tangent space is identified with . The geodesics are great circles and the geodesic distance is , for any . The corresponding exponential map at , , is given by

while the corresponding logarithm map at , , is given by

One can use the exponential map and the logarithm map to define a coordinate system for the in the following way. Given , we define the logarithmic coordinate as

where and forms an orthonormal basis for . Write . In addition, we define

which aligns with , and for . Note that for any , we have where . Here first aligns a direction with the true diffusion direction and then represents it by its logarithmic coordinate.

We now present the asymptotic results. Now, write for , and . We have . Also, let be the -th order derivative of with respect to for . Let and . Also, denote .

Under the assumptions 1-10 laid out in Section S5.2 of the Supplemental Material which are all standard technical conditions (except for Assumption 1 which is to ensure the representation of the geodesic distance as a function of the working coordinate system), we have the following theorem.

The proof of the Theorem ? can be found in Section S5.2 of the Supplemental Material.

## 5Fiber tracking

For dMRI, fiber tractography can be classified as deterministic and probabilistic methods. Deterministic methods [?] track fiber bundles by utilizing the principal eigenvectors of tensors, while probabilistic methods [?] use the probability density of diffusion orientations. Most deterministic methods assume one single diffusion tensor in each voxel, and hence are unable to handle voxels with crossing fibers. In view of this, this section develops a deterministic tracking algorithm that allows for multiple or no principal diffusion directions in a voxel.

The proposed algorithm can be seen as a generalization of the popular Fiber Assignment by Continuous Tracking (FACT) [?] algorithm. A brief description of FACT is as follows. Tracking starts at the center of a voxel (Voxel 1 in Figure ? left panel) and continues in the direction of the estimated diffusion direction. When it enters the next voxel (Voxel 2 in Figure ? left panel), the track changes its direction to align with the new diffusion direction and so on. This tracking rule may produce many short and fragmented fiber tracts due to either a wrongfully identified isotropic voxel or spurious directions which go nowhere. In addition, it cannot determine which direction to follow in case there are multiple directions in a voxel, which happens in crossing fiber regions.

To address these issues, we modify the above procedure in the following manner. Given a current diffusion direction (we refer to the corresponding voxel as the current voxel), the voxel that it points to (we refer to this voxel as the destination voxel) may have (i) at least one direction; (ii) no direction (i.e., isotropic). In case (i), we will first identify the direction with the smallest angular difference with the current direction. If its separation angle is smaller than a pre-specified threshold (e.g., ), we enter the destination voxel and tracking will go on along this direction. See Figure ? (Middle). On the other hand, if the separation angle is greater than the threshold, or case (ii) happens, we deem that the destination voxel does not have a viable direction. In this case, tracking will go along the current direction if it finds a viable direction within a pre-specified number of voxels. The number of voxels that are allowed to be skipped is set to be in our numerical illustrations. See Figure ? (Right). On the other hand, the tracking stops at the current voxel if no viable directions within a pre-specified number of voxels can be found. The detailed tracking algorithm is described in Algorithm S5 in the Supplemental Material.

As for the choice of starting voxels, also known as seeds, there are two common strategies. One can choose seeds based on tracts of interest and start the tracking from a region of interest (ROI). This approach is based on knowledge on ROI and may not give a full picture of the tracts of interest if there are diverging branches. The other approach is the brute-force approach, where tracking starts from every voxel. It usually leads to a more comprehensive picture of tracts at a higher computational cost. The proposed algorithm can be coupled with either strategy.

Combining the voxel-wise estimation method in Section 3 and the direction smoothing procedure in Section 4 gives the proposed DiST method.

## 6Simulation study

Extensive simulation experiments have been conducted to evaluate the practical performances of DiST. They are reported inSection S6 of the Supplemental Material. Overall, the DiST method provided highly promising results.

## 7Real data application

In this section, we apply the proposed methodology to a real dMRI data set, which was obtained from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database (www.loni.ucla.edu/ADNI). The primary goal of ADNI has been to test whether serial MRI, positron emission tomography (PET), other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of mild cognitive impairment (MCI) and onset of Alzheimer’s disease (AD). In the following, we use an eddy-current-corrected ADNI data set of a normal subject for illustration of our technique.

This data set contains 41 distinct gradient directions with -value set as . In addition, there are 5 images (corresponding to ), forming in total 46 measurements for each of the voxels. To implement our technique, we require estimates of ’s and . We first estimate and for each voxel by ML estimation based on the 5 images. Then we fix as the median of estimated ’s for voxel-wise estimation of the diffusion directions. Since the original voxels contain volume outside the brain, we only take median over a human-chosen set of voxels. The estimated is 56.9.

In this analysis, we focus on a subset of voxels (), which contains the intersection of corpus callosum (CC) and corona radiata (CR). This region is known to contain significant fiber crossing [?]. See Figure (Left) for a fiber orientation color map of one of the five -planes. Within the whole focused region, ’s have mean 1860.1 and standard deviation 522.7.

We then apply voxel-wise estimation to individual voxels followed by the DiST-mcv procedure. Distributions of the estimated number of diffusion directions are summarized in Table ?. For comparison purposes, we also fit the single tensor model with the commonly used regression estimator [?].

The tracking results are produced by applying the proposed tracking algorithm to the estimated diffusion directions from DiST-mcv, which represents the DiST procedure with chosen by the median cross-validation score (see Sections S3 and S6 of the Supplemental Material), and those from the single tensor model estimation. For visualization purposes, we present the longest 900 tracts in Figure . From anatomy, the CC has a mediolateral direction while the CR has a superoinferior orientation. They are clearly shown in both tracking results. In these figures, reconstructed fiber tracts are colored by a RGB color model with red for left-right, green for anteroposterior, and blue for superior-inferior. Thus, one can easily locate the CC and the CR as the red fiber bundle and the blue fiber bundle respectively. Tracking result based on DiST-mcv shows clear crossing between mediolateral fiber and the superoinferior fiber (in the figure, the crossing of red and blue fiber tracts). From neuroanatomic atlases and previous studies, [?] conclude that there are several fiber populations with crossing structure in this conjunction region of CC and CR, which matches with the tracking based on DiST-mcv. However, the single tensor model estimation can only reconstruct one major diffusion direction in each voxel and thus the corresponding tracking result does not show crossing structure. Instead, the CC (red fiber bundle) is blocked by the CR (blue fiber bundle) and this leads to either termination of the CC fiber tracts or significant merging of the CC and the CR fiber tracts instead of the known crossing structure. To give further illustration, Figure 2 shows the locations of the CC, the CR and the region of crossing fibers (Cross). One can see that estimated directions based on DiST-mcv reproduces the crossing fiber structures between the CC and the CR, while the result based on single tensor model tends to connect the CC and the CR fibers.

Moreover, the green fiber on top of the CC represents the cingulum bundle. Both fiber tracking based on DiST and single fiber model produce clear and sensible reconstruction of cingulum bundle. All these features match with neuroanatomic atlases and provide a good demonstration of our proposed method.

As shown by Figures Figure 2 and ?, when comparing with the results obtained by the single tensor model,DiST produces more biologically sensible and interpretable tracking results. This provides more reliable information on brain connectivity and in turn could lead to better understanding of neuro-degenerative diseases such as Alzheimer’s disease and autism as well as better detection of brain abnormality, such as deformation and neuron loss in white matter regions.

## 8Discussion

Using tensor estimation to resolve crossing fiber can be problematic, due to the inability of estimating multiple diffusion directions by the single tensor model and the non-identifiability issue in multi-tensor model. In this paper, we take a different route by focusing on the estimation of diffusion directions rather than the diffusion tensors. We develop the corresponding direction smoothing procedure and fiber tracking strategy, together called DiST, along this route. Our technique gives promising empirical results in both simulation study (see Section S6 of the Supplemental Material) and real data analysis.

The procedure we presented works well even with moderate number of gradient directions (a few tens), as long as the number of distinct crossing fibers within a voxel is not larger than three. With HARDI data, which can have up to a couple of hundreds gradient directions, rather than modeling the direction distribution within a tensor framework, we can estimate the fiber orientation distribution nonparametrically [?].

ApplyingDiST to multiple images from ADNI (either from the same subject over time or from multiple subjects) and then relating the tracking results with clinical outcomes such as cognitive measures would provide valuable information about the role of white matter connectivity in initiation and progression of Alzheimer’s disease and dementia. Although this is an important direction of research, it is beyond the scope of this paper which focuses on developing a statistical procedure to denoise dMRI data and to provide better tracking results. We plan to explore more sophisticated applications of the proposed procedure in our future research.

## Acknowledgement

Data collection and sharing for this project was funded by the Alzheimer’s Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant U01 AG024904). ADNI is funded by the National Institute on Aging, the National Institute of Biomedical Imaging and Bioengineering, and through generous contributions from the following: Abbott, AstraZeneca AB, Bayer Schering Pharma AG, Bristol-Myers Squibb, Eisai Global Clinical Development, Elan Corporation, Genentech, GE Healthcare, GlaxoSmithKline, Innogenetics, Johnson and Johnson, Eli Lilly and Co., Medpace, Inc., Merck and Co., Inc., Novartis AG, Pfizer Inc, F. Hoffman-La Roche, Schering-Plough, Synarc, Inc., as well as non-profit partners the Alzheimer’s Association and Alzheimer’s Drug Discovery Foundation, with participation from the U.S. Food and Drug Administration. Private sector contributions to ADNI are facilitated by the Foundation for the National Institutes of Health (www.fnih.org). The grantee organization is the Northern California Institute for Research and Education, and the study is coordinated by the Alzheimer’s Disease Cooperative Study at the University of California, San Diego. ADNI data are disseminated by the Laboratory for Neuro Imaging at the University of California, Los Angeles. This research was also supported by NIH grants P30 AG010129, K01 AG030514, and the Dana Foundation. The authors would like to thank Professor Owen Carmichael for making available the data and his valuable comments.