1 Introduction
The study of the brain structure and how it functions is an old yet essential topic. As an invivo, noninvasive and noradiationincluded medical imaging technology, Diffusionweighted MRI (DMRI) has become widely used by neuroscientists to probe architecture of biological tissues. The detailed anatomy of white matter fiber tracts in the brain can be revealed by making use of diffusion displacement measures of water molecules in DMRI, since in white matter, water usually diffuses faster along the fiber bundles which makes the water diffusion anisotropic. Reconstructing white matter fiber tracts is of great value in understanding the brain structure and its functionality (Mori, 2007; Sporns, 2010)
. One of the earliest DMRI technology is the wellknown diffusion tensor imaging (DTI)
(Basser et al., 1994). Recently, high angular resolution diffusion imaging (HARDI) has become of great interest in accurate extraction of fiber tracts information, e.g., the orientation of intravoxel crossing fibers, through diffusion measurements along a large number of gradient directions (Tuch et al., 2002). The white matter fiber tracts can then be reconstructed by applying tracking algorithms (tractography) (Basser et al., 2000; Wong et al., 2016).There are many models developed for HARDI data, including angular reconstruction models such as fiber orientation distribution (FOD) (Tournier et al., 2004, 2007) and diffusion orientation distribution function (ODF) (Tuch, 2004; Descoteaux et al., 2007, 2009). In this paper, we focus on one particular model: the FOD, which is designed to describe the angular dependence of the axonal fiber bundles by making use of sharp geometric features (Jensen et al., 2016). White matter fiber reconstruction through tractography can greatly benefit from a good estimation of the FOD.
One commonly used method in FOD estimation is based on the spherical deconvolution framework described in Tournier et al. (2004, 2007), where the key to its success is a parsimonious representation of the FOD in a suitable basis. Early studies of the FOD estimation are mostly based on spherical harmonics (SH) representations of the FOD (Tournier et al., 2004, 2007). However, due to the global support of the SH basis functions, they do not represent the FOD efficiently when spiky features exist, which is usually the case in welldefined fiber tracts. Recently, Yan et al. (2018) has proposed a new FOD estimation method called SNlasso based on spherical needlets (SN) representations of the FOD. As smooth spherical functions localized in both space and frequency, spherical needlets form a tight spherical frame and thus provide parsimonious and stable representations for spherical functions with localized sharp peaks such as the FOD (Narcowich et al., 2006b, a). By introducing spherical needlets, the FOD estimation problem is reformulated as a linear regression problem where the needlet coefficients of the FOD are the variables to estimate. An
penalty is imposed on these needlet coefficients to ensure the sparse representation. Moreover, a nonnegative constraint on the estimated FOD is also imposed to guarantee the nonnegativeness of the probability density function FOD. In this paper, we focus on the
SNlasso method in FOD estimation, since it would lead to more accurate and stable reconstruction results of the FOD and preserve its sharp features.Most of the existing methods for fiber reconstruction, including the FOD estimation method SNlasso, are developed in each voxel and perform voxelwise fiber reconstruction independently. However, the spatial constraint, which is the structural nature of the HARDI data, is essentially ignored by these methods. Following Rao et al. (2016), the direction of single fiber bundles are expected to change smoothly along the dominant fiber direction from one voxel to its neighborhood voxels, until it reaches the boundaries between fiber tracts. Similarly, the directions of crossing fiber bundles are expected to vary smoothly along their own dominant fiber directions from one voxel to its surrounding voxels in the fiber crossing region.
Until recently, a few methods have been proposed for fiber reconstruction where spatial constraint is incorporated for potential improvement. Some of these methods are based on the wellknown propagationseparation method (Polzehl & Spokoiny, 2000, 2006), which is developed specifically for model functions with large homogeneous regions and sharp discontinuities. Although the original propagationseparation method is only applicable in denoising the usual images, several extensions of this method are developed in the area of fiber reconstruction. For example, smoothing and estimation on DTI data via an anisotropic structural adaptive smoothing procedure (Tabelow et al., 2008), smoothing raw HARDI images based on structural adaptive smoothing in both voxel space and diffusiongradient space (Becker et al., 2012, 2014), spatial and adaptive analysis of neuroimaging data using a multiscale adaptive regression model (MARM) which integrates propagationseparation method with voxelwise statistical modeling (Li et al., 2011), spatial and adaptive estimation of ODF using a penalized multiscale adaptive regression model (PMARM) which integrates MARM, propagationseparation method and Lasso (Rao et al., 2016). Other spatial smoothing methods using strategies different from propagationseparation method include spatial smoothing and estimation of DTI (Liu et al., 2013; Yu & Li, 2013), spatial smoothing of HARDI with Bayesian spatial regularization (Raj et al., 2011), spatial smoothing of estimated diffusion directions via clustering (Wong et al., 2016) and so on.
In this paper, we propose a novel Nearestneighbor Adaptive Regression Model (NARM) to conduct spatially smoothed and locally adaptive estimation of FODs across all voxels. NARM can be regarded as a natural integration of the extended propagationseparation method and the voxelwise FOD estimation method SNlasso. Specifically, for the FOD estimation in each voxel, instead of optimizing a single likelihood function based on the voxel itself (i.e., the linear regression model in SNlasso), we incorporate its neighborhood information by optimizing the sum of weighted likelihood functions based on the voxel and its neighborhood voxels. The weights are chosen appropriately to account for spatial proximity and potential heterogeneity due to different fiber configurations, measured by both the spatial distance between voxels and the similarity between FOD estimates.
Similar to the PMARM method proposed by Rao et al. (2016), we adaptively update these weights and also use a stopping rule to avoid oversmoothing. However, NARM differs from PMARM
in three major aspects. First, we use Hellinger distance in measuring similarity between FODs, where Hellinger distance is proposed specifically for probability distribution function and is more sensitive to peak value difference. In contrast, the
PMARM method uses the distance in FOD similarity measurement. Second, we use a simple yet effective stopping rule based on the minimum distance between the voxel and its nearest neighbors, whereas the PMARM method relies on a more complicated stopping rule with hyperparameters not easy to tune in practice. Third, we implement weightrescaling for voxels which are either extremely easy or extremely hard to be smoothed, to avoid potential oversmoothing or undersmoothing issues. On the contrary, the PMARM method does not cover this. In addition to the extension of the propagationseparation method, we also extend the SNlasso to the cases where diffusion measurements are sampled with multiple values. This is actually the format of DMRI dataset in the WUMinn Human Connectome Project database (Van Essen et al., 2013).We conduct extensive simulation studies based on twodimensional and threedimensional regions of synthetic data to demonstrate the effectiveness of our proposed NARM method, and to make comparisons with the voxelwise SNlasso method and the PMARM method. We also show the advantages of FOD estimation based on multiple values over FOD estimation based on single value. Moreover, we implement NARM to one real 3T DMRI dataset of healthy individual from the WUMinn Human Connectome Project database (Van Essen et al., 2013). The tractography results based on the estimated FODs from our proposed method tend to show more clear patterns of fiber tracts with less noise, compared to the estimated FODs from voxelwise method.
The rest of paper is organized as follows. Section 2 formally defines the problem and introduces the voxelwise FOD estimation method SNlasso. Section 3 proposes the Nearestneighbor Adaptive Regression Model (NARM) in detail. Section 4 and 5 present the experimental results on synthetic and real datasets respectively. Section 6 concludes our work and discusses future directions. Technical details can be found in an Appendix. Additional details are deferred to a Supplementary Material.
2 Voxelwise FOD estimation
2.1 Spherical convolution model
We estimate the FOD function for each voxel based on the observed diffusion weighted measurements corresponding to gradient directions. These diffusion weighted measurements can be regarded as the noise corrupted evaluations of the diffusion signal function at these gradient directions. Thus, it is crucial to understand the relationship between the FOD and the diffusion signal function . Following Tournier et al. (2004, 2007), the diffusion signal function can be represented by the spherical convolution of the FOD function : , a symmetric spherical probability density function, and the response function : , an azimuthal symmetric kernel function which represents the diffusion signal along a single fiber bundle and is assumed to be the same across voxels. Hereafter we assume the response function is known. In practice, can be estimated by applying tensor models onto voxels with a single dominant fiber bundle.
The formal definition of the spherical convolution model at voxel follows:
(1) 
Following Descoteaux et al. (2007); Yan et al. (2018), the real symmetric spherical harmonic (SH) functions form an orthonormal basis for the square integrable real symmetric functions on . By the orthonormality of the real symmetric SH basis, the diffusion signal function can be expressed as:
where , and are the spherical harmonic coefficients of the diffusion signal function, response function and FOD respectively.
Now we assume the real symmetric SH basis up to level can well approximate the diffusion signal function, response function and FOD. The observed diffusion weighted measurements can then be modeled as:
(2) 
where is an matrix ( is the number of basis functions) consisting of SH function evaluations at gradient directions; is an diagonal matrix with diagonal elements in blocks of length ; is an vector of FOD’s SH coefficients to be estimated; and is an error vector.
2.2 SNlasso method
In SNlasso, we represent the FOD by symmetric spherical needlets (SN) basis. The real symmetric SN functions form a tight frame for the square integrable real symmetric functions on , and they are localized in both scale and space. Thus, the FOD can be well approximated by a small set of SN functions since most of the FODs are either constant or with a few localized sharp peaks.
Now we set the maximum level of SN basis and being the corresponding number of SN basis functions. is selected such that the first SN functions can linearly represent the SH functions up to level . In this case, the SH coefficient vector in (2) can be further expressed as , where is an transition matrix and can be precomputed (Narcowich et al., 2006a; Yan et al., 2018), and is an vector of FOD’s SN coefficients. (2) can then be rewritten accordingly:
(3) 
Based on the sparsity nature of the FOD’s SN coefficients , SNlasso proposes to estimate by an penalized linear regression model:
(4) 
Here is the lasso regularizer with tuning parameter , which leads to sparse representation of . is the matrix of SH functions evaluated on a dense spherical grid, and guarantees that the estimated FOD on this dense grid is nonnegative. Notice both the design matrix and the constraint matrix are the same across all the voxels and can be computed in advance.
A computationally efficient algorithm based on ADMM (Alternating Directions Method of Multipliers) (Boyd et al., 2011) is developed for the optimization problem in (4), where the algorithm details are given in Appendix 7.1. A peak detection algorithm is also proposed to efficiently identify the major fiber directions according to the FOD estimation, where more details are given in Yan et al. (2018).
3 Spatiallysmoothed FOD estimation
The FOD estimation method SNlasso is a voxelwise method which is implemented for each voxel independently by using only the DMRI information from the voxel itself. However, the spatial smoothness nature of HARDI data across voxels is essentially ignored by this method. For example, for single fiber bundle in a voxel, its direction is expected to vary smoothly along the major fiber direction from the voxel to its neighborhood, except at the boundaries between fiber tracts. For crossing fiber bundles in a voxel, their directions are also expected to smoothly follow their own major fiber directions respectively from the voxel to its surrounding voxels. The spatial constraint in HARDI data is indeed very important and worthy exploring to further improve the reconstruction of fiber tracts.
To take the spatial constraint into consideration, we propose a new spatial smoothing method in FOD estimation, named Nearestneighbor Adaptive Regression Model (NARM). The fundamental idea is the implementation of local likelihood modeling (Fan & Gijbels, 1996), that is, instead of optimizing a single likelihood function at voxel based on the observed diffusion measurements only from voxel itself (i.e., Eq. (4) up to a constant), we incorporate signals from its neighborhood voxels as well and optimize the sum of weighted likelihood functions within a neighborhood centered at voxel . Here without loss of generality, we assume all the voxels are located in the Cartesian grid with unit distance between adjacent voxels. Let denote a (spherical) neighborhood with radius centered at voxel , and denote a weight function of pair such that and . The FOD’s SN coefficients corresponding to neighborhood can then be estimated from the penalized weighted linear regression model:
(5) 
where is the weighted average of observed diffusion measurements within the neighborhood . Note that when is set to be 0, the neighborhood contains only one voxel , and Eq. (5) actually coincides with Eq. (4) in the voxelwise method SNlasso. Following Yan et al. (2018), we select the largest value of tuning parameter in penalty such that the residual sum of squares will not be significantly improved if the value of is further decreased (See Appendix 7.2 for more details).
One popular approach in local likelihood modeling is the propagationseparation method (Polzehl & Spokoiny, 2000, 2006), which aims to figure out a maximal local neighborhood of each sample point where the local homogeneity assumption of parameter is justified. Inspired by the propagationseparation method, we adaptively determine the weight and estimate the FOD’s SN coefficients across all voxels. For each voxel , we construct a sequence of nested (spherical) neighborhoods with increasing radii . At , the estimated coefficients are exactly estimated voxelwisely by Eq. (4). At the th step, the weights are calculated by the method as detailed in the following sections where the information of the previously estimated coefficients is included, and the current estimated coefficients can be directly obtained afterwards. In practice, we set the series of as , , , where the default choice of is 1.15 and the default choice of is 10 in twodimensional region and 6 in threedimensional region. The intuition of different choices of is that we make the numbers of voxels within the largest neighborhood in twodimensional region and threedimensional region similar.
3.1 Weight selection
The selection of the weights , is crucial in finding a good estimation of FOD and hence a reasonable reconstruction of fiber tracts. Ideally, the weights should reflect the spatial distance between voxels, assuming that the closer the voxels are, the more useful signal information they can borrow from each other. Moreover, the weights are also expected to provide a similarity measure between the underlying true fiber directions in voxels. That is, the weight should be large if the fiber directions between two voxels are similar, so that the truly useful information can be provided by one voxel to the other; on the other hand, if the true fiber directions between two voxels are very different, the weight should be decreased to prevent incorporating spurious signals which leads to bad estimations. However, these true fiber directions are never known to us, and a natural idea is to use the estimated FODs to help determine these weights.
Let denote the similarity measure between the estimated FODs in voxels and at step . The formula for calculating the weights is written as
where and are two nonnegative kernel functions.
The weight can be regarded as the combination of two independent weights: the location weight and the similarity weight . In NARM, we choose which only considers the voxels within the neighborhood . Moreover, it gives more weight to the voxel closer to the voxel , indicating that the voxel would probably contain more useful information to help FOD estimation in the voxel . We also choose which takes small value if the estimated FODs in the voxel and differ significantly. is a positive scalar which gives less weight to on when it increases.
We set the similarity measure at step , , as the Hellinger distance between the previously estimated FOD (represented on a dense grid specified by ) in voxel , , and the previously estimated FOD in voxel , . Following the formula of Hellinger distance, can be expressed as
where the square root of vector is applied elementwisely. The Helllinger distance measures specifically on distribution functions such as FOD. Compared with a commonly used distance, the distance, the Hellinger distance includes the square root operation in its formula which emphasizes more on the peak values of distribution functions. As a result, a single fiber voxel will probably be more similar to its single fiber neighbors rather than its crossing fiber neighbors, so that the signals from crossing fiber neighbors may not affect the smoothing procedure of this single fiber voxel too much. On the other side, a crossing fiber voxel can be limitedly influenced by the signals from its single fiber neighbors as well. In this case, the transferral of neighborhood information is more targeted and efficient, and the boundary between the estimated single fiber region and crossing fiber region will be more distinct.
3.2 Weight rescaling
The proposed weight , however, does not handle the extreme cases very well. For example, if the FOD estimate at voxel is already wellsmoothed, when increases, both the location weight and the similarity weight for increase (the increment of is due to the increment of similarity within neighborhood), which somewhat overemphasizes the smoothness nature and results in potential oversmoothing. On the contrary, for badlysmoothed FOD estimate at voxel , we expect it would largely incorporate its neighborhood’s information to make the estimation consistent with its neighbors. However when increases, does not increase too much due to the dissimilarities between voxel and its neighbor , and is usually much lower than those weights on wellsmoothed voxels, which in some degree overlooks the smoothness nature and leads to potential undersmoothing.
Our proposed NARM method can well adjust these extreme smoothing circumstances. The main idea is that we rescale the similarity weight so that it is not very sensitive to those extreme cases. Specifically, at step , we calculate the minimum value among the Hellinger distances between the estimated FOD in voxel and the estimated FODs in its nearest neighboring voxels , which are the nearest four voxels (left, right, top, bottom) in twodimensional region or six voxels (left, right, top, bottom, front, end) in threedimensional region. We denote this minimum nearestneighbor distance as . Due to the smoothness nature of the HARDI data, it is natural to assume that each voxel has at least one similar nearest neighbor. In this case, a relatively large value of may indicate undersmoothness and a relatively small value may indicate wellsmoothness or oversmoothness. Thus, if we set to be the percentile of across all voxels , then for those percent voxels ( is usually set below ) with , we need to boost their smoothing procedure since they are considered to be undersmoothed, and we should slow down the smoothing procedure for those percent voxels with since it is very likely that they are already wellsmoothed or even oversmoothed.
The detailed method is presented in the following. Suppose for the voxel , , where is one of the nearest neighbor of , then is set to be if , or if , and all the other , , , are decreased/increased proportionally. In other words, the new similarity measure between the estimated FODs in voxels and at step is
and the new weight is
Here we can see that, for those undersmoothed voxels with , we have and thus an increased similarity weight, which could effectively accelerate the smoothing procedure. On the contrary, for those already wellsmoothed or even oversmoothed voxels, we have a decreased similarity weight leading to the deceleration of the smoothing procedure.
3.3 Stopping rule
A stopping rule is also essential in NARM to prevent potential oversmoothing. Our proposed stopping rule is based on the minimum nearestneighbor distance as well, and the basic idea is that we terminate smoothing procedure for voxel when stops decreasing as increases. Concretely, when we observe , we stop the smoothing procedure for voxel and set for . The intuition behind this stopping rule is straightforward: when stops decreasing, the extent of similarity between voxel and its neighbors cannot be improved anymore, which means voxel has already been wellsmoothed and the smoothing procedure should be terminated for voxel to prevent oversmoothing. In practice, the smoothing procedure for different voxels may stop at different steps, indicating the heterogeneity nature of the degree of smoothness across the voxel region.
The computation costs of weight calculation and stop checking are both much cheaper than the cost of FOD estimation by SNlasso. Hence, compared with voxelwise FOD estimation method when applied across all voxels, the NARM method costs at most times more, where is the total number of smoothing steps. In fact, the cost can be much less than times in practice (usually around times), since no FOD estimation is required if the smoothing procedure of the voxel is already terminated by the stopping rule.
3.4 Related works
Rao et al. (2016) proposed Penalized Multiscale Adaptive Regression Model (PMARM) for spatial and adaptive estimation of ODF. As a general smoothing framework, it can be easily migrated to FOD estimation. Similar to NARM, PMARM also replaces the single diffusion measurement within each voxel by the weighted average of diffusion measurements, adaptively updates the weights within nested neighborhoods with increasing radius, and uses a stopping rule to avoid oversmoothing. However, PMARM has three major drawbacks compared with our method. Here for convenience, we use the same notations in PMARM as those in our NARM method.
Firstly, the similarity measure in PMARM is defined as the distance between the previously estimated FODs at voxel and , whereas in our method, the Hellinger distance is used for calculating the similarity measure. As pointed out before, compared with the Hellinger distance, the distance performs less effectively in quantifying the similarity between estimated FODs. Specifically, the distance may struggle in distinguishing the single fiber voxel and crossing fiber voxel, resulting in blurred boundary between the estimated single fiber region and crossing fiber region, since spurious signals could be transferred more easily between these two regions. As a consequence, we may have difficulties in the following fiber orientation detection and fiber tracking.
The second drawback is that PMARM does not cover the rescaling process for the similarity weight . As mentioned before, some extreme cases may not be well treated by PMARM. For example, for those already wellsmoothed voxels, the weights seem to overemphasize the smoothness nature and may lead to oversmoothing. On the other hand, the weights of those undersmoothed voxels seem too low to incorporate enough neighborhood information and may lead to undersmoothing.
Finally, PMARM uses a different stopping rule. At step , let denote the distance between previously estimated FOD and currently estimated FOD . Rao et al. (2016) states that, if , where is a prespecified positive scalar, then some “bad” signals from neighborhood voxels may be incorporated and lead to a significant change in the estimated FOD, and we thus should stop the smoothing procedure and set and for the voxel . On the contrary, if , we continue with the smoothing procedure for the voxel . In practice, is set as , where is the upper percentile of the distribution, and is a positive scalar, which is chosen to be the median of where and are preselected voxels.
However, as pointed out by Rao et al. (2016), in practice it can be very difficult to determine a good value. Actually, it is hard to determine a global criterion based on a uniform measure such as which fits every voxels well. In some voxels where the original FOD estimates are quite reasonable, the FOD estimations can become oversmoothed very quickly, however the global stopping rule may not be able to stop the smoothing procedure in time. While in some other voxels where the FODs are not originally wellestimated, they could have a long way towards good estimates, whereas the global stopping rule may stop the smoothing procedure too early.
In fact, the selection of stopping rule is crucial in spatial smoothing method, and extensive simulation results show that the ultimate FOD estimates are quite sensitive to the hyperparameter selection in the stopping rule. However, in the expression of : , the distribution is actually not justifiable, and the tuning parameter selection (both the percentile of distribution and the constant value ) is somewhat adhoc in practice, where the way of preselecting voxels used in estimating value is not described in detail. Compared with the stopping rule in PMARM, the stopping rule in NARM based on minimum nearestneighbor distance is more reasonable and applicable, where no additional hyperparameters are needed to be tuned. The essential difference between these two methods is that instead of using a global criterion based on a uniform measure, we use a local criterion adjusted to each voxel, thus can well handle different smoothing cases for different voxels.
3.5 Further extension: multiple values
In this section, we briefly introduce a simple but powerful extension of the SNlasso method, which can be easily incorporated into our NARM method. Recall that the original SNlasso is built on diffusion weighted signals sampled with a single value. However, in many recent DMRI dataset such as the WUMinn Human Connectome Project (HCP) database (Van Essen et al., 2013), the diffusion weighted signals are sampled with multiple values. Specifically, a list of values is prespecified. For each value , we have its corresponding observed diffusion weighted measurements based on gradient directions. Using the original SNlasso, we model the measurements as:
(6) 
where is an matrix constructed on the gradient directions, is an diagonal matrix based on the value , is an matrix independent with , is an vector of FOD’s SN coefficients to be estimated, and is an error vector.
To combine the sampled diffusion signals from different values, a straightforward way is to combine the linear models (6) by stacking the response vectors, design matrices and error vectors. Let , be an response vector, be an design matrix and be an error vector. The combination of the linear models becomes:
(7) 
Now the FOD’s SN coefficients can be estimated by
(8) 
We can see that, compared with the original SNlasso using signals sampled from a single value, more samples can be included in the extended SNlasso, therefore the statistical efficiency in estimating is improved and its variation is decreased. Simulation results in the next section will show that, by stacking the responses and design matrices for different values, the FOD estimation result becomes not only more accurate than the result based on a single value, but also more robust to the selection of tuning parameter in penalty.
4 Synthetic experiments
In this section, we use extensive simulation studies based on twodimensional and threedimensional regions of interest (ROIs) with synthetic DMRI data to evaluate the performance of our proposed NARM method, and to compare it with other FOD estimation methods including the voxelwise SNlasso method and the PMARM method.
4.1 Simulation setup
In ROIs, we simulate the noiseless diffusion weighted signals according to the spherical convolution model (1). For voxels with at least one fiber bundles, the true FOD function is constructed by a weighted summation of dirac functions: , , , with being the volume fraction of the th fiber bundle and being the spherical coordinate of th fiber direction. In our simulation studies, we set for voxels with two fiber bundles. For voxels with zero fiber bundles (i.e., isotropic diffusion), the true FOD function is constant on the sphere: . The response function is built based on the single tensor model (Le Bihan & Warach, 1995; Basser & Jones, 2002; Mori, 2007): , where we set , value or , , and the ratio between and as 10 for voxels with at least one fiber bundles and 1 for voxels with zero fiber bundles throughout our simulation studies. We generate the noiseless diffusion weighted signals along gradient directions, where the gradient directions are based on the icosphere meshes and is set to be either 41 or 81. Finally, we add independent Rician noise (Gudbjartsson & Patz, 1995; Hahn et al., 2006)
to the noiseless diffusion weighted signals to generate the simulated diffusion weighted measurements. The signaltonoise ratio
, where is the Rician noise level, is set to be 20 in our simulation studies.In the simulation studies, we generate true fiber directions and their corresponding observed diffusion weighted measurements on a twodimensional or threedimensional voxel grid. These fiber directions are changing smoothly across voxels, and within each voxel, there is either one dominant fiber direction, or two crossing fiber directions, or no fiber directions at all (i.e., isotropic diffusion). We estimate the FOD at each voxel by using several FOD estimation methods. Specifically, in 2D ROI simulations, we conduct a comprehensive exploration of the FOD estimation methods, including the voxelwise SNlasso method (referred to as SNlasso(voxelwise)), the NARM method with different selections of in weight rescaling (referred to as NARM()), the NARM method with certain components missing (referred to as NARM(no stopping&rescaling), NARM(no stopping) and NARM(no rescaling)), the NARM method with the stopping rule replaced by the one in PMARM (referred to as NARM(stopping, )), and the PMARM method with different selections of in stopping rule (referred to as PMARM()). We also include the voxelwise SNlasso method applied on the noiseless diffusion signals (referred to as SNlasso(noiseless)) as a gold standard. In 3D ROI simulations, for convenience, we only focus on a small set of FOD estimation methods which are representative of the aforementioned methods. Fiber peaks can be further extracted from the estimated FODs by using the peak detection algorithm proposed in Yan et al. (2018).
We use a variety of numerical metrics to examine the performance of each method. They include: (1) The proportions of voxels where the number of detected fiber bundles are either correct, overestimated or underestimated, reported under “Co.”, “Ov.” and “Un.” respectively in the tables. (2) Among the “correct” voxels, the median angular error between the identified fiber direction(s) and true fiber direction(s), shown under “Err.” in the tables. (3) The mean and standard deviation (across all voxels in ROI) of the Hellinger distance between the estimated FOD and the true FOD (the true FOD function projected onto an SH basis with
), and between the estimated FOD and the estimated FOD from noiseless signals. They are shown under “Mean Hdist.” and “SD Hdist.” in the tables.4.2 2D ROI simulation I
In the first simulation study of ROI, we generate two fiber bundles on a 2D voxel grid, where onefiber voxels (72 such voxels) and twofiber voxels (28 such voxels) are included. Figure 1(a) shows the true FODs represented by SH basis on this 2D ROI. We can see clearly the two fiber bundles, where one is from the top left to the bottom right, and the other is from the top right to the bottom left, with the crossing fiber region in the middle of the ROI. We set the value and the number of gradient directions . In some real datasets, the values and the numbers of gradient directions may be larger than our settings here, leading to less challenging FOD estimation problems. However, the goal of this simulation study is to examine the performance of various FOD estimation methods and reveal the effectiveness of each component in the smoothing procedure, thus using the relatively difficult case can help our understanding better. We refer to this simulation as the 2D ROI simulation I.
We examine the performance of FOD estimation methods listed in Section 4.1. Specifically, we set in NARM(), in NARM(stopping, ), and in PMARM(). Extensive experiments show that the results are not very sensitive to the choice of when it varies between 0.1 and 0.2, whereas they are quite sensitive to the choice of . Here we preselect the values of which already leads to relatively wellperformed results. In practice, values are not easy to be selected and the selection approach is somewhat adhoc and not described explicitly by Rao et al. (2016). We set in for NARM and for PMARM (adjusted for scale difference).
Figure 1(b)(f) present the estimated FOD images on the ROI. In general, the FOD estimates with smoothing procedure (Figure 1(c)(f)) significantly outperform the voxelwise FOD estimates (Figure 1(b)) in terms of recovering the consistently orientated FODs across voxels especially in the crossing fiber region and removing the spurious peaks. The NARM estimates without stopping rule and weight rescaling (Figure 1(d)) seem heavily oversmoothed, where the boundary between single fiber region and crossing fiber region looks unclear. Compared with the PMARM estimates (Figure 1(e)(f)), the NARM estimates (Figure 1(c)) perform better in identifying the peak directions. For example, only the NARM estimates can successfully recover the crossing fiber structure in the voxel (6, 9) (assume the coordinates of the voxel in the bottom left is (0, 0)), which seems very challenging to recover due to its small crossing angle and its boundary location between single fiber region and crossing fiber region. Take single fiber voxels (1, 6) and (10, 6) as another example. the NARM method leads to much less noisy estimates than the PMARM method does, where the PMARM estimates seem to misidentify these two single fiber voxels as crossing fiber voxels.
Table 1 and S.1 (Table S.1 is in the Supplementary Material) present the quantitative summaries of the FOD estimation accuracy by using the success rate&angular error and the Hellinger distance respectively, where the detailed description can be found in Section 4.1. We can see that the NARM estimates are not very sensitive to the selection of (the proportion of extreme cases whose weights need to be rescaled), and performs best due to its lowest median angular error in the 2fiber cases. Both tables reveal that the NARM estimates substantially reduce the failure rates, median angular errors and mean Hellinger distances from the voxelwise SNlasso estimates, indicating that NARM is able to efficiently exploit the spatial smoothness and reduce noise for better peak detection and angle estimation. Compared with the PMARM estimates with three preselected values in stopping rule, NARM outperforms all of them due to its perfect success rates among both 1fiber and 2fiber voxels and smaller mean Hellinger distances, whereas PMARM seems to have blurred boundary between single fiber region and crossing fiber region due to its relatively high proportion of overestimated 1fiber voxels and underestimated 2fiber voxels. Finally, if the stopping rule is not implemented, the NARM estimates will have significantly lower successful rates and larger angular errors in the 1fiber cases due to oversmoothing. Even if we replace the original stopping rule in NARM by the stopping rule from PMARM (referred to as stopping) with preselected values, it is still less effective in distinguishing the single fiber region and crossing fiber region, indicating the superior performance of the stopping rule in NARM not only in higher estimation accuracy but also in more convenient estimation process (with no hyperparameters hard to tune in practice).
(a) True FODs  (b) SNlasso (voxelwise) 
(c) NARM ()  (d) NARM (no stopping&rescaling) 
(e) PMARM ()  (f) PMARM () 
1Fiber  2Fiber  
Estimator  Co./Ov./Un.  Err.  Co./Ov./Un.  Err. 
SNlasso(noiseless)  1.00/0.00/0.00  1.55  1.00/0.00/0.00  1.00 
SNlasso(voxelwise)  1.00/0.00/0.00  2.62  0.78/0.18/0.04  8.51 
NARM()  1.00/0.00/0.00  2.61  1.00/0.00/0.00  3.81 
NARM()  1.00/0.00/0.00  2.61  1.00/0.00/0.00  3.59 
NARM()  1.00/0.00/0.00  2.61  1.00/0.00/0.00  3.72 
NARM(no stopping&rescaling)  0.81/0.19/0.00  5.34  0.96/0.04/0.00  3.09 
NARM(no stopping)  0.76/0.24/0.00  4.62  0.96/0.04/0.00  3.04 
NARM(no rescaling)  1.00/0.00/0.00  2.61  1.00/0.00/0.00  3.81 
NARM(stopping, )  0.97/0.03/0.00  2.69  0.93/0.03/0.04  3.72 
NARM(stopping, )  0.94/0.06/0.00  2.72  0.96/0.00/0.04  3.34 
PMARM()  0.93/0.07/0.00  2.53  0.96/0.00/0.04  4.08 
PMARM()  0.92/0.08/0.00  2.61  0.96/0.00/0.04  3.44 
PMARM()  0.89/0.11/0.00  2.64  0.96/0.00/0.04  3.44 
4.3 2D ROI simulation II
In the second simulation study of ROI (2D ROI simulation II), two fiber bundles are generated on a 2D voxel grid, including 22 zerofiber (isotropic) voxels, 66 onefiber voxels and 12 twofiber voxels. The true FODs represented by SH basis are shown in Figure 2(a). We still set the value and the number of gradient directions .
In Figure 2 and Table 2, S.2 (Table S.2 is in the Supplementary Material), we focus on the FOD reconstruction results particularly on the zerofiber (isotropic) regions. As can be seen from Figure 2, the NARM estimates (Figure 2(c)) and the PMARM estimates with (Figure 2(e)) can identify more consistent zerofiber regions than the voxelwise SNlasso estimates (Figure 2(b)) do, where the isotropic structures in the voxels (6, 2) and (10, 6) are successfully recovered by these two estimates. Moreover, NARM outperforms PMARM due to its perfect success rate and much lower median angular error in recovering the crossing fiber region as well as its lower mean Hellinger distance. Note that the PMARM estimates with (Figure 2(f))) perform much worse than the PMARM estimates with (Figure 2(e)) in terms of successful rate in zerofiber region, indicating the difficulty in selection for PMARM in practice. Finally, removing stopping rule and weight rescaling from NARM will result in serious overfitting problem, where no zerofiber voxels are detected, as shown in Figure 2(d). Replacing the original stopping rule in NARM by the stopping rule from PMARM will worsen the success rate of FOD reconstruction in all three regions with zero fiber voxels, single fiber voxels and crossing fiber voxels.
(a) True FODs  (b) SNlasso (voxelwise) 
(c) NARM ()  (d) NARM (no stopping&rescaling) 
(e) PMARM ()  (f) PMARM () 
0Fiber  1Fiber  2Fiber  

Estimator  Co./Ov.  Co./Ov./Un.  Err.  Co./Ov./Un.  Err. 
SNlasso(noiseless)  1.00/0.00  1.00/0.00/0.00  1.92  1.00/0.00/0.00  1.00 
SNlasso(voxelwise)  0.86/0.14  1.00/0.00/0.00  2.72  0.92/0.00/0.08  8.88 
NARM()  0.95/0.05  1.00/0.00/0.00  2.93  1.00/0.00/0.00  4.76 
NARM()  0.95/0.05  1.00/0.00/0.00  2.86  1.00/0.00/0.00  4.76 
NARM()  0.95/0.05  0.98/0.02/0.00  2.86  1.00/0.00/0.00  4.76 
NARM(no stopping&rescaling)  0.00/1.00  0.79/0.21/0.00  4.00  1.00/0.00/0.00  4.37 
NARM(no stopping)  0.68/0.32  0.79/0.21/0.00  4.18  1.00/0.00/0.00  3.70 
NARM(no rescaling)  0.95/0.05  1.00/0.00/0.00  2.93  1.00/0.00/0.00  4.76 
NARM(stopping, )  0.73/0.27  0.97/0.03/0.00  2.48  0.92/0.00/0.08  5.69 
NARM(stopping, )  0.73/0.27  0.95/0.05/0.00  2.57  0.92/0.00/0.08  4.28 
NARM(stopping, )  0.73/0.27  0.94/0.06/0.00  2.57  0.92/0.00/0.08  3.70 
PMARM()  0.95/0.05  1.00/0.00/0.00  2.57  0.92/0.00/0.08  6.87 
PMARM()  0.41/0.59  0.97/0.03/0.00  2.55  0.92/0.00/0.08  6.41 
4.4 3D ROI simulation
In the third simulation study (3D ROI simulation), we generate two fiber bundles on a 3D voxel grid, where 52 zerofiber voxels, 163 onefiber voxels and 285 twofiber voxels are included. Figure S.1(a) shows the true FODs represented by SH basis on the 5 slices of the 3D ROI. Figure S.1(b) presents the fiber tracts based on the true FODs within the entire 3D ROI. Two fiber bundles can be easily recognized from these plots, where one follows top left to bottom right, into to outofpage direction, and the other follows top right to bottom left, into to outofpage direction. We consider value at and , and the number of gradient directions . Based on the simulation results in previous sections, we select a subset of representative FOD estimation methods for further examination in this simulation study, including the voxelwise SNlasso method, the NARM method with , the NARM method without stopping and rescaling, and the PMARM method with and . In NARM, we set in at and at . The intuition is that we need to slightly slow down the smoothing process as value increases, since stronger diffusion signals from larger value can themselves lead to better FOD estimates, and no smoothing process as intensive as that for weak diffusion signals is needed.
In the main text, we only report part of the simulation results at . For other results at and the results at , see Sections S.1.3 and S.1.4 of the Supplementary Material, respectively. We can see from Figure 3, S.2 and Table 3, S.3 (Table S.3 is in the Supplementary Material) that the NARM method leads to the most reasonable reconstruction of FODs: It has nearly perfect success rates among 0fiber, 1fiber and 2fiber voxels, the lowest median angular errors, and the smallest mean Hellinger distance. Among the other methods, the voxelwise SNlasso method performs poorly in accurate detection of peak direction in crossing fiber region, the PMARM method does not deal with many challenging cases very well and results in less faithful recovery of FODs in terms of success rate and angular error, and the NARM method without stopping and rescaling tends to heavily oversmooth the zero fiber and single fiber regions and results in noisy estimates of FODs with a number of spurious peaks.
We use the tractography results in Figure S.3 to further illustrate and justify our findings. To reconstruct the fiber trajectories, we use the peak detection algorithm (Yan et al., 2018) to extract peaks from the estimated FODs and then apply the deterministic tracking algorithm (Wong et al., 2016) to the extracted peaks. In Figure S.3, the fiber tracts are colored by a RGB color model with red for lefttoright, green for uptodown and blue for intotooutofpage. We also report the total number of fiber tracts and the median of fiber lengths in Figure S.3. As can be seen from Figure S.3, the tracking result based on the NARM method (Figure S.3(c)) seems to be the closest to the underlying truth (Figure S.3(a)), where it shows very clear fiber crossing region with almost no redundant fiber tracts. This implies that the NARM method is able to efficiently exploit the spatial smoothness nature and reduce noise for more accurate fiber orientation detection and angle estimation. On the contrary, the tracking result based on the voxelwise SNlasso method (Figure S.3(b)) looks very messy especially in the crossing fiber region, and the small value of its median fiber length 4.11 indicates the incoherence of fiber directions identified from its FOD estimates. Moreover, both the fiber tracking based on the PMARM method with and (Figure S.3(e)(f)) produce more noisy reconstruction of fiber tracts. Finally, Figure S.3(d) shows that many spurious fiber tracts will be generated due to oversmoothing if the NARM method estimates FODs without stopping and rescaling.
For value , NARM again has visually the best FOD estimates and fiber tracking result, and the highest success rates in all 0fiber, 1fiber and 2fiber regions. See Section S.1.4 of the Supplementary Material for more details. One more thing we notice is that the FOD estimation results from value at are better than those from . This indicates that stronger diffusion signals would bring more useful information and lead to more satisfactory FOD reconstruction.
(a) True FODs  (b) SNlasso (voxelwise) 
(c) NARM()  (d) NARM(no stopping&rescaling) 
(e) PMARM(c=0.05)  (f) PMARM(c=0.10) 
0Fiber  1Fiber  2Fiber  

Estimator  Co./Ov.  Co./Ov./Un.  Err.  Co./Ov./Un.  Err. 
SNlasso(noiseless)  1.00/0.00  1.00/0.00/0.00  1.61  1.00/0.00/0.00  1.77 
SNlasso(voxelwise)  0.94/0.06  1.00/0.00/0.00  2.47  0.84/0.15/0.01  6.75 
NARM()  0.96/0.04  0.98/0.02/0.00  2.35  1.00/0.00/0.00  3.31 
NARM(no stopping&rescaling)  0.13/0.87  0.75/0.25/0.00  3.27  1.00/0.00/0.00  2.74 
PMARM()  0.90/0.10  0.95/0.05/0.00  2.63  0.94/0.06/0.00  3.93 
PMARM()  0.40/0.60  0.90/0.10/0.00  2.81  0.97/0.03/0.00  3.42 
We further explore the effectiveness of NARM extension with multiple values. We first create gradient directions equally spaced on the sphere, and we generate diffusion weighted measurements at on 41 out of 81 gradient directions (equally spaced), and diffusion weighted measurements at on the rest 40 gradient directions. We compare the performance of voxelwise SNlasso method on three datasets: the measurements at only, the measurements at only, and the measurements at both and . We also set the parameter in the RSSbased selection (see Appendix 7.2 for more details) to be either or to investigate the effects of selection to the FOD estimation result, where larger will lead to larger and thus sparser estimate. The performance of the NARM method applied on the diffusion weighted measurements at both and is also included to mimic our treatments in real data applications.
Several important findings can be concluded from Figure S.7, S.8, S.9 and Table S.6, S.7 in the Supplementary Material. First, the SNlasso estimates on diffusion measurements at is robust to the selection, with almost identical estimation results at and , whereas the SNlasso estimates at performs poorly at with no identified crossing fibers at all, and the SNlasso estimates at tends to oversmooth the zero fiber region at . Second, even considering only the best case for each dataset, the SNlasso estimates at still significantly outperforms the estimates at either or in terms of success rate, median angular error and mean Hellinger distance. This indicates the effectiveness of multiple value stacking in accurate fiber orientation estimation. In fact, sampling at a low value only would probably cause the loss of direction information in the FOD estimates and make the fiber tracking hard, while sampling at a high value only would probably bring in too much noise into FOD estimation. Third, the spatial smoothing procedure in NARM can further improve the recovery of FODs with much smaller median angular error in crossing fiber region. Its resulting fiber tracking plot is almost visually identical to the fiber tracking plot based on true FODs.
4.5 Summary
We have the following findings from 2D ROI simulation I & II and 3D ROI simulation. Firstly, the FOD estimation methods with spatial smoothing procedure outperform the voxelwise SNlasso method, especially in crossing fiber regions. Secondly, compared with PMARM (with best chosen stopping parameter based on truth), the NARM method is able to produce more accurately orientated FOD estimates, and leads to more faithful recovery of fiber tracts. Thirdly, stopping rule is crucial for adaptive smoothing methods such as NARM and PMARM. Without stopping, the estimation can be overly smoothed, while stopping too early can lead to high variability. Our proposed stopping rule (including weight rescaling) for NARM is effective and robust to the hyperparameter . In contrast, the stopping rule in PMARM is highly sensitive to the hyperparameter , making it hard to use in practice. Lastly, it is beneficial to combine observed diffusion measurements at multiple values, which would improve not only the accuracy of fiber direction estimation, but also the robustness of performance to the selection.
5 Real DMRI data experiments: HCP
In this section, we implement our proposed method into a real DMRI dataset obtained from the WUMinn Human Connectome Project (HCP) database (https://www.humanconnectome.org/). The primary goal of HCP has been to map macroscopic human brain circuits and their relationship to behavior in a large population of healthy adults. The HCP dataset includes highresolution 3T MR scans from young healthy adult twins and nontwin siblings (ages 2235) using four imaging modalities: structural images (T1w and T2w), restingstate fMRI (rfMRI), taskfMRI (tfMRI), and high angular resolution diffusion imaging (dMRI). See Van Essen et al. (2013) for more details of its scanner design, parameter setting and implementation. In the following, we analyze DMRI dataset from one participant, a female in the age 3135, to illustrate our methodology.
We use eddycurrentcorrected diffusion images with diffusion weighted signals measured on a 3D grid (voxel sidelength: , , ) along 270 distinct gradient directions, including 90 gradient directions under , 90 gradient directions under and 90 gradient directions under . In addition, there are also 18 images and we estimate the image intensity and the Rician noise level based on them. The median signaltonoise ratio across voxels is 20.7. We estimate the response function by using the method proposed in Yan et al. (2018)
, where the estimated leading eigenvalue is
and the ratio between the leading and the minor eigenvalue is 6.1.5.1 Roi I
We focus on a subregion (ROI I) of the DMRI dataset, with slices 5468, slices 118132 and slices 61:75. This ROI contains the intersections of corpus callosum (CC), cingulum bundle (CB) and superior occipitofrontal fasciculus (SOF) according to the FSL Jülich Histological Atlas (Eickhoff et al., 2005, 2006, 2007). The left panel of Figure S.10(a) shows the fiber orientation colormap on the entire brain at
slice 71, where the ROI is indicated by the white box. The fiber orientation colormap depicts the direction of principle eigenvector under single tensor model, where red stands for lefttoright direction, green stands for uptodown direction and blue stands for intotooutofpage direction. The fiber orientation colormap is further modulated by the FA values such that it becomes darker as the FA value gets smaller. The rest plots in Figure
S.10(a) are the zoomedin version of the fiber orientation colormap, FA map and MD map on the ROI at slice 71. In the FA/MD maps, the color also becomes darker as the FA/MD value gets smaller. The details of FA/MD calculation are in Section S.2.1 of the Supplementary Material.We apply the voxelwise SNlasso method to three distinct datasets: the diffusion measurements along 90 gradient directions at , the diffusion measurements along 90 gradient directions at , and the diffusion measurements along 270 gradient directions at , and (i.e., the entire set of measurements). We also apply the NARM method onto the entire set of measurements to demonstrate the effectiveness of our proposed spatial smoothing method. We set in . Here we focus on the ROI at slice 71 since this subregion contains several crossing fiber tracts as well as CSF as shown below.
Figure S.10(b) shows the images of FOD estimates on slice 71. We can see that all the methods produce satisfactory FOD estimates in the single fiber region, with clear peak directions consistent with the ones suggested by the fiber orientation colormap. However, the method performance differs quite a lot in the crossing fiber region. We take the bottom right part of this region as an example (indicated by the white box), where we have low FA values in the middle as well as low MD values. This is indeed a strong indication of crossing fibers rather than CSF since CSF is believed to have large MD values due to its fast water diffusion. From the zoomedin subregion plots in Figure S.11, it can be seen that in the crossing fiber region, the SNlasso at leads to less pointy estimates with less direction information contained, whereas the SNlasso at tends to produce very noisy estimates without clear pattern. The SNlasso estimates based on all three values outperform the previous two estimates in more coherent fiber crossings in this subregion. And by applying the NARM method, the FOD estimates are further improved with more clearly orientated fiber directions and less spurious fiber peaks. We also look into the top right corner of the region in Figure S.10, which is likely to be CSF due to small FA values and large MD values. The best FOD estimates to capture this isotropic structure is from the NARM method, where the largest number of isotropic fibers are identified in this subregion. In contrast, the FOD estimates without spatial smoothing procedure contain less isotropic fibers, and among them the SNlasso estimates at performs the worst, and the possible reason is that it may be hard to get rid of the noisy information from the diffusion measurements at only.
Figure 4 and S.12 present the fiber tracking results at two different perspectives by applying the deterministic tracking algorithm (Wong et al., 2016) to the extracted peaks from FOD estimates on the entire ROI. These tracking results reveal the fiber tracts in the atlas, where the red fiber bundle corresponds to CC, the bigger green bundle on the left corresponds to SOF and the smaller green bundle on the top right corresponds to CB. Among the tracking results based on four different estimates, the result based on the NARM estimates (Figure 4(d) and S.12(d)) shows the clearest crossings between the red fiber tracts and the green fiber tracts. Moreover, both the red and green fiber tracts in the result based on the NARM estimates are more tightly clustered, smoothly curved and parallel to each other than those in other results. Compared with the tracking results based on the SNlasso estimates at or at only, the tracking result based on the SNlasso estimates at all three values shows more coherent and less noisy fiber tracts, with a larger value of median fiber length at 5.84. In contrast, the SNlasso estimates at leads to difficulties in fiber tracking due to its lack of direction information, thus the resulting fiber lengths have the lowest median value 4.30. Meanwhile, the SNlasso estimates at involves too much noise, and we can see clearly from Figure S.12(b) that the tracking result contains spurious fiber tracts in the top right part of the ROI which is believed to be CSF according to the FA and MD maps. In conclusion, when comparing with the SNlasso based on single value measurements, the SNlasso based on multiple value measurements can produce more coherent and clear tracking result, and the NARM method in which spatial smoothing procedure is implemented can further improve the smoothness and parallelism of the tracking result to make it more biologically meaningful and interpretable.
(a) SNlasso(voxelwise)()  (b) SNlasso(voxelwise)() 
(fiber number: 3556, fiber length: 4.30)  (fiber number: 5445, fiber length: 5.32) 
(c) SNlasso(voxelwise)()  (d) NARM() 
(fiber number: 4596, fiber length: 5.84)  (fiber number: 5007, fiber length: 6.80) 
To further demonstrate the effectiveness of our proposed spatial smoothing method NARM, we also apply the voxelwise SNlasso method and the NARM method to DMRI dataset in HCP retest data from the same participant. We use the FLIRT module in FSL to align the ROI I in test data to retest data, and it appears that the locations of ROI I in test data and retest data are almost identical, where we have slices 5468, slices 118132 and slices 6175 in test data, and slices 5468, slices 119133 and slices 6074 in retest data.
We also use a variety of numerical metrics to examine the consistency of FOD estimates in test and retest data by using the voxelwise SNlasso method and the NARM method respectively. Since test and retest data are obtained with independent errors, we believe that if the estimates in test and retest data are similar, they should be similar to the underlying truth. In other words, the FOD estimates should be more accurate if the estimates in test and retest data are more consistent. The numerical metrics we have used include: (1) The proportions of voxels where the numbers of detected fiber bundles match between test and retest data (number of voxels with a certain detected fiber number in both test and retest data, divided by number of voxels with this detected fiber number in either test or retest data), reported under “Matched” in the tables. (2) Among the “matched” voxels, the median angular difference between the identified fiber direction(s) in test and retest data, shown under “Diff.” in the tables. (3) The mean and standard deviation (across all voxels in ROI) of the Hellinger distance between the estimated FODs in test and retest data. They are shown under “Mean Hdist.” and “SD Hdist.” in the tables. (4) Relative difference of numbers of connected fiber tracks across subregions in ROI between test and retest data (for each pair of subregions, relative difference of numbers of connected fiber tracks across the corresponding two subregions in the pair between test and retest data is calculated, and the mean of relative differences is reported). They are shown under “Mean Track Diff.” in the tables.
0Fiber  1Fiber  2Fiber  

Estimator  Co./Ov.  Co./Ov./Un.  Err.  Co./Ov./Un.  Err. 
SNlasso(noiseless)  1.00/0.00  1.00/0.00/0.00  1.61  1.00/0.00/0.00  1.77 
SNlasso(voxelwise)  0.94/0.06  1.00/0.00/0.00  2.47  0.84/0.15/0.01  6.75 
NARM()  0.96/0.04  0.98/0.02/0.00  2.35  1.00/0.00/0.00  3.31 
NARM(no stopping&rescaling)  0.13/0.87  0.75/0.25/0.00  3.27  1.00/0.00/0.00  2.74 
PMARM()  0.90/0.10  0.95/0.05/0.00  2.63  0.94/0.06/0.00  3.93 
PMARM()  0.40/0.60  0.90/0.10/0.00  2.81  0.97/0.03/0.00  3.42 
Table 4 reports the quantitative summaries of the consistency of FOD estimates in test and retest data. We can see that the estimates by using NARM are more similar between test and retest data than the estimates by using voxelwise SNlasso, according to higher matching rate and lower median angular difference in zerofiber, onefiber and twofiber regions, as well as lower mean Hellinger distance. Moreover, the fiber tracks obtained from the deterministic tracking algorithm shows more consistent pattern between test and retest data in NARM than voxelwise SNlasso, where the mean relative difference of numbers of connected fiber tracks across subregions is smaller in NARM. As previously pointed out, the more consistent the FOD estimates in test and retest data are, the more close these estimates are to the underlying truth, and thus the more accurate the FOD estimation method is. Therefore, the NARM method has outperformed the voxelwise SNlasso method in this real data experiment.
5.2 Roi Ii
The second subregion (ROI II) we are going to explore in the same HCP dataset is a region with slices 5675, slices 86105 and slices 6685. According to the FSL Jülich Histological Atlas, this region contains the intersections of four tracts including corpus callosum (CC), cingulum bundle (CB), superior occipitofrontal fasciculus (SOF) and corticospinal tract (CT). The ROI at slice 82 is shown by the white box in the left panel of Figure S.13(a) and its corresponding zoomedin fiber orientation colormap, FA map and MD map are also shown in Figure S.13(a).
We present the FOD estimates on slice 82 in Figure S.13(b). According to the FA and MD maps in Figure S.13(a), the middle part of the region (highlighted by the white box) is very likely to contain crossing fibers due to both small FA and MD values. This is corroborated by the NARM estimates shown in the zoomedin subregion plots in Figure S.14, where they clearly capture the crossing between a fiber in the lefttoright direction and a fiber in the uptodown direction, as indicated by the colormap as well. The other voxelwise SNlasso estimates again contain many spurious peaks and do not show clear directionality. The right part of the region in Figure S.13 is probably CSF which shows large MD values and small FA values. Among all the four FOD estimates, the NARM estimates depict the isotropic region the best where FODs in most part of this subregion are estimated as isotropic. On the contrary, the SNlasso estimates at perform the worst in identifying isotropic voxels.
The reconstructed fiber tracks are presented in Figure 5 and S.15 at two different perspectives. All the four fiber tracks in the atlas are clearly shown in the tracking results, where the blue tracts correspond to CT, the red tracts correspond to CC, the green tracts in the top correspond to CB and the green tracts in the bottom correspond to SOF. Similar to our previous findings, the NARM estimates also produce the smoothest and most coherent trajectories, whereas the SNlasso estimates at struggles in finding appropriate fiber directions and results in the smallest median fiber length, and the SNlasso estimates at generate too many noisy tracts in the top right corner of the ROI.
(a) SNlasso(voxelwise)()  (b) SNlasso(voxelwise)() 
(fiber number: 7933, fiber length: 4.21)  (fiber number: 12249, fiber length: 6.13) 
(c) SNlasso(voxelwise)()  (d) NARM() 
(fiber number: 9374, fiber length: 6.43)  (fiber number: 9075, fiber length: 7.65) 
6 Conclusions
In this paper, we have proposed a Nearestneighbor Adaptive Regression Model (NARM) to adaptively estimate the FODs across all voxels using DMRI data. NARM is able to incorporate spatial constraint in DMRI data for fiber reconstruction by adaptively borrowing information from neighboring voxels. NARM has improved the existing PMARM method (Rao et al., 2016) in three major aspects: the Hellinger distance in FOD similarity measurement which makes FODs with different peak orientations more distinguishable, the weight rescaling process which well handles the extreme cases in smoothing, and the stopping rule based on the minimum nearestneighbor distance which is easy to implement in practice and is effective in preventing oversmoothing. Moreover, NARM is further extended to the cases where diffusion measurements are sampled with multiple values.
Experimental results on synthetic data demonstrate that NARM results in better FOD reconstruction than the voxelwise SNlasso method and the existing PMARM method, where it leads to more accurate and less noisy fiber orientation detection, and more coherent and smoother fiber tracking result. In experiments on real 3T DMRI datasets, the NARM estimates lead to clearer crossing fiber regions and more biologically sensible and interpretable tractography results than the voxelwise SNlasso method. Moreover, the SNlasso estimates based on diffusion measurements at multiple values perform better than those at single value in terms of the coherence and noise level of fiber tracts. Therefore, the NARM estimates based on diffusion measurements at multiple values would provide more useful information on brain connectivity and lead to better understanding of brain structure and functionality.
One direction for future research is to further improve the similarity measure between FODs in weight calculation. The Hellinger distance we are currently using is able to better identify the difference between peak values of FODs than the distance. However, the Hellinger distance is not very sensitive to the peak orientation shift of FODs since it is not very aware of the peak location. Some locationaware distances between density functions can be explored in the future, for example, the earthmover distance. However, the computational efficiency is a big issue since the computation cost of the earthmover distance is much expensive than that of the Hellinger distance.
7 Appendix
7.1 ADMM algorithm in SNlasso
We first simplify the notation in the optimization problem (4). Let
Then (4) can be rewritten as
To solve this optimization problem using ADMM algorithm, we notice that it can be further written as
where , , when (elementwise) and otherwise.
The scaled augmented Lagrangian is
where and are dual variables.
The ADMM algorithm for SNlasso is presented as follows. We first initialize , , , . We also need to specify , which in practice is recommended to be (Wahlberg et al., 2012). For step until convergence:

minimization.
Comments
There are no comments yet.