A Novel Richardson-Lucy Model with Dictionary Basis and Spatial Regularization for Isolating Isotropic Signals

Abstract

Diffusion-weighted magnetic resonance imaging is a non-invasive imaging method that has been increasingly used in neuroscience imaging over the last decade. Partial volume effects (PVEs) exist in sampling signal for many physical and actual reasons, which lead to inaccurate fiber imaging. We overcome the influence of PVEs by separating isotropic signal from diffusion-weighted signal, which can provide more accurate estimation of fiber orientations. In this work, we use a novel response function (RF) and the correspondent fiber orientation distribution function (fODF) to construct different signal models, in which case the fODF is represented using dictionary basis function. We then put forward a new index Piso, which is a part of fODF to quantify white and gray matter. The classic Richardson-Lucy (RL) model is usually used in the field of digital image processing to solve the problem of spherical deconvolution caused by highly ill-posed least-squares algorithm. In this case, we propose an innovative model integrating RL model with spatial regularization to settle the suggested double-models, which improve noise resistance and accuracy of imaging. Experimental results of simulated and real data show that the proposal method, which we call iRL, can robustly reconstruct a more accurate fODF and the quantitative index Piso performs better than fractional anisotropy and general fractional anisotropy.

About the Authors:

Tiantian Xu

Affiliation: Institute of Information Processing and Automation, Zhejiang University of Technology, Hangzhou, Zhejiang, China

Yuanjing Feng

* E-mail: [email protected]

Affiliation: Institute of Information Processing and Automation, Zhejiang University of Technology, Hangzhou, Zhejiang, China

Ye Wu

Affiliation: Institute of Information Processing and Automation, Zhejiang University of Technology, Hangzhou, Zhejiang, China

Qingrun Zeng

Affiliation: Institute of Information Processing and Automation, Zhejiang University of Technology, Hangzhou, Zhejiang, China

Jun Zhang

Affiliation: Institute of Information Processing and Automation, Zhejiang University of Technology, Hangzhou, Zhejiang, China

Jianzhong He

Affiliation: Institute of Information Processing and Automation, Zhejiang University of Technology, Hangzhou, Zhejiang, China

Qichuan Zhuge

Affiliation: Zhejiang Provincial Key Laboratory of Aging and Neurological Disorder Research, Wenzhou Medical University, Wenzhou, Zhejiang, China

Introduction

Magnetic resonance imaging (MRI) can offer important insights into brain disease [1]. Only diffusion-weighted MRI (DW-MRI) can provide a unique, non-invasive technique to study the microscopic structure of brain white matter (WM) in vivo [2–4]. DW-MRI provides valuable information about the fiber architecture of tissue by measuring the diffusion of water in three-dimensional space [5, 6]. An early form of this technique, i.e., diffusion tensor imaging (DTI) [7], is widely used in clinics and provides fiber orientations of WM based on principal eigenvector of that tensor [8] and many useful quantitative indexes, including mean diffusivity (MD), fractional anisotropy (FA) [9, 10], and so on. The major shortcoming of the representative DTI is its inability to characterize the orientations of crossing and branching neural tracts in brain, especially fiber tracts with intersected diffusion orientations or partial volume averaged within a voxel [11–13]. Many recent high angular resolution diffusion imaging (HARDI) techniques have been proposed to recover the complex WM geometry [14]. Most of these methods consider water-molecule diffusion as a function of direction, such as Q-ball imaging (QBI) [15], diffusion spectrum imaging (DSI) [16] and spherical deconvolution (SD) [17], which have all conquered the limitations of DTI. However, the data acquisition times for QBI and DSI are exorbitant [18] because of the high sampling numbers required to construct the full diffusion propagator. Given the linearity and sensitivity to multi-model diffusion [11], considerable interests have been generated with the model-free SD, which is based on convolution between fiber response function (RF) and fiber orientation distribution function (fODF). Although the SD shows both good angular resolution and short computational time, the defects emerge when facing PVEs and the imaging quality is degraded by spurious directions and negative orientations caused by the truncation of high-order harmonics and ill-posed solution, even in noise-free data [19].

Partial volume effects (PVEs) were put forward by Timo Roine et al. firstly [20]. It usually appears on the border of different tissues. The brain contains complex WM and non-WM tissues, such as gray matter (GM) and cerebrospinal fluid (CSF), which have different diffusion properties. Thus, the PVEs phenomenon is particularly obvious in human brain [12, 21, 22]. For PVEs, the SD method induces some changes on RF, but this does not solve the PVEs in essence. An informed constrained spherical deconvolution (iCSD) has been proposed to improve the estimation of fODF under non-WM PVEs by modifying RF to account for non-WM PVEs locally [23]. However, the iCSD method can’t correctly resolve fiber crossing angles of less than 60° under significant non-WM PVEs. Some authors have included an isotropic compartment in their signal models but these methods both require multiple b-value acquisitions and distinguish the signal of different tissues [24]. In other methods based on spherical deconvolution, the isotropic signal is dampened by using an iterative RL deconvolution algorithm [25]. Falvio et al. [19] infer that fODF can be represented by several discrete Dirac delta functions on unit sphere and propose a new spherical model based deconvolution approach to solve fiber crossing. They consider isotropic diffusion and anisotropic diffusion signal and combine both of two components. Dell’Acqua et al. suggest a new term, fiber orientation function (FOF) to represent the weights of anisotropic and isotropic diffusion [26]. However, the FOF, as a combination of anisotropic and isotropic diffusions, can’t really take them apart. Consequently, The use of FOF is difficult. Isotropic signal existing in GM or CSF misleads the algorithms to produce spurious peaks in FOF. In this framework, Dell’Acqua et al. further combine RL spherical deconvolution algorithm with an adaptive regularization technique to yield damped Richardson-Lucy (dRL) algorithm in spherical deconvolution, aiming to attenuate isotropic signal while reducing spurious and non-physical fiber orientations in regions affected by PVEs [27]. dRL has its limitations. Given the different degrees of attenuation in each voxel, small FOF portions are more likely to be preserved in a low isotropic volume fraction, which leads to spurious fiber orientations [26]. Notably, the method based on RL has settled the highly ill-conditioned problem of least squares algorithm. However, in the absence of constrains of solution, even small changes in the acquired signal (e.g., MR noise) can lead to nonphysical results [17, 28]. A number of regularization algorithms have thus been developed. Yap et al. [29] develop a spatially non-negative sparse representation framework and then present an algorithm for solving l0 sparse group representation problem and apply it to tissue signal separation problem [30]. While the computational cost and intractable computation will arise when the models are more sophisticated. To make full use of spatially constraints of brain fibers, many global tractography methods considered PVEs [4, 31]have been proposed in the last two years. But there are always many disadvantages, including computing space occupied, convergence property, sub-optimal solution and so on. There is a long way to realize global tractography perfectly.

In this study, we consider a new spherical deconvolution model (hereafter denoted as iRL), which can effectively isolate isotropic signal from each DW signal. A new quantitative index is put forward to distinguish WM and non-WM of human brain, and the quantitative results of that index are better than those FA and GFA. We also propose a novel method, based on RL to efficiently reconstruct the above fiber architecture and yield high-quality fODF results. The true fODFs are gathers of delta function pointing along fiber orientations, and zero in all other orientations [32]. Thus, a dictionary basis is introduced to represent the fODF, which effectively helps to separate isotropic signal and renders the coefficients of fODF sparse. Finally, we integrate total variation regularization and ℓ1 norm regularization on the above framework to smooth noise and suppress spurious fiber orientations at the same time. To compare the performances with existing methods, the experiments are conducted on simulated and real data using the proposed method in compared with several kinds of methods, which are introduced in detail in the following sections.

Materials and Methods

Spherical deconvolution

Spherical deconvolution based on a relatively simple model of signal generation has been recently developed to overcome the limitations of diffusion tensor model in resolving multiple fiber orientations and to improve tractography reconstruction. The motivation of this proposed method is to recover fODF directly from DW signal without prior assumption or estimation about the number of fibers representing the information about diffusion [31]. The DW signal S can be assumed as a superposition of anisotropic and isotropic signal, which can be regard as two different diffusion models for three reasons. Firstly, the sampling voxels have a relatively large volume. On the border of WM and non-WM, the signal of each WM is affected by isotropic signal from non-WM, such as GM and CSF, which is known as PVEs phenomenon. The second is that isotropic diffusion exists in WM. Given that isotropic diffusion is weaker than anisotropic diffusion in WM, the diffusion in WM is always considered as anisotropy. The third one is that the complex structures of fibers such as orthogonal fibers lead to increased isotropy. Generally, the signal contributed by isotropic tissue is usually not included in spherical deconvolution models [21]. However, to facilitate calculation, researchers often try not to differentiate between the two parts and instead only make some changes in RF. The best solution is to put the two parts of DW signal segregated. In this work, we try to separate the two different parts of DW signal which would produce better imaging results especially in the DW signal existed PVEs.

Let be unit sphere domain and SO(3) be rotation group in . The anisotropic diffusion signal is modeled by convolution between a kernel R ∈ L2(SO(3)) and a function , which respectively represent the signal response function (RF) and fODF ideally composed of N Dirac delta functions for n bundles of fibers [33]. We assume that the isotropic signal in each voxel is the same, thus the spherical deconvolution operator can be expressed as:(1)where g are diffusion gradient orientations containing I directions and , S(g) are diffusion attenuation signal along g, are isotropic signal, which are equal along each gradient orientation and overlooked in most medical imaging cases, the dot stands for standard (Euclidean) dot product in , v is unit sphere (v also represent the discretized directions of unit sphere in the following parts), R(g ⋅ v) is the RF describing DW signal intensity and f(v)dv is probability measure used to model fODF over [31]. The fODF contains all desired anisotropic information regarding both orientations of various fiber populations that may be presented and their respective volume fractions [34]. A common case is that we have N fibers in one voxel, where N is a limited natural number, and the corresponding fODF is the sum of N Dirac delta functions on sphere weighted by corresponding volume fractions. The form of fODF enables the separation of two diffusion models. Regarding the anisotropic signal, RF and fODF are defined as usual.

The novel fODF estimation with double models

Basser et al. [7] indicate that the signal in a pulsed gradient spin echo depends on diffusion sensitive coefficient b and diffusion tensor D, the relation is:(2)

This relation relies on assumptions that the compartments have equal relaxation rates and water density, and the exchange between volumes can be neglected on the time scale of measurement [35]. Where S(gi) denotes the diffusion attenuation signal along i-th diffusion gradient orientation gi. D is diffusion tensor, which describes the simplest model of diffusion in axon fiber bundles. The value of D is the extremum direction of diffusion, which can decide the degree of water diffusion. The RF [36–38] derived from the above signal relation has a certain inaccuracy. Improving the precision of RF is of great advantage in the subsequent RL iterative model. Thus, we use the original Eq (2) as our RF.

Our final goal is to construct the fODF which characterizes the relative likelihood of water diffusion along a given direction. Most of HARDI methods do not account for PVEs caused by non-WM tissues and orthogonal fibers. Signals contributed by GM or CSF both are actually isotropic compartments and are included in the existing model of spherical deconvolution. To accurately reconstruct brain connections from DW signal, we should properly model the different types of water diffusion signal [39]. In order to make calculate easy, we discretize the process of spherical deconvolution (the discretized directions are still expressed as v). The reconstruction of SD method is computed as linear combination of the diffusion measurements [11]. The fODF can be reasonably considered as two main terms, viz. anisotropic and isotropic parts. Thus, incorporating these contributions by using double models is possible. Based on algebraic theory, we can combine the parts of anisotropy and isotropy. Thus at each voxel, the special deconvolution can be expressed as:(3)where v are unit direction vectors which are acquired by averaging discretization of unit spherical surface along J directions and , and , D1 is diffusion tensor of fibers (FA = 1, MD = 0.0007mm2/s), whose value is to ensure the maximum anisotropy, Rani(g ⋅ v)(j) is the RF along j-th sample direction vj, which is a disc-shaped RF generated by the model presented in Eq (2) for a single fiber. There are J RFs oriented along each sampling direction. Thus, Rani(g ⋅ v) is an I × J matrix, Riso = S0 e−tr(bgT D2 g) is a column vector of length I containing the signal of isotropic compartment. D2 is isotropic tensor of DW signal (FA = 0, MD = 0.0007mm2/s). Thus, the final RF is an I × (J + 1) matrix. Naturally, fODF can be expressed as f(v) = [fani fiso], and it consists of two parts, the first J rows fani stand for the anisotropy. The last row fiso provides information related to isotropy. The fODF can be expressed more clearly as f = f(v) = [fani(v1), fani(v2),…fani(vJ), fiso].

To simplify the numerical solution, the fODF constructed by SD is originally formulated using spherical harmonics basis. Actually, the proposed method can be implemented using a number of well-characterized dictionary basis sets, which are flexible unimodal basis functions. This relationship can be expressed as:(4)where u are unit direction vectors along L (with L ≥ J) directions and , which are used to increase the accuracy of fiber directions, Φ(v, u) is a (J + 1) × (L + 1) matrix which will be illustrated in the next step. f(v) and c = c(v) denote (J + 1) × 1 and (L + 1) × 1 column vectors composed of estimated values of fODF and the coefficient of fODF, respectively. Notably that the diffusion measurements c also consist of two parts, the first L rows cani show the information about anisotropy; the last raw ciso represents the information related to isotropy. We can use this variable denoted as Piso to quantify the intensity of isotropy of each voxel. Piso can take place of the value of FA and GFA as well as conveys the message even better than them to some extent. Removing the isotropic part of each voxel inevitably increases the accuracy of fiber imaging. Once we have acquired the diffusion signal S(g) and , the unknown part fODF f can be computed using the iRL model.

Dictionary basis representation

SD has been proven to produce a good imaging result. [17] proposed to express SD directly in spherical harmonics (SH) domain, so the operation can be reduced to a simple set of matrix multiplications. Simultaneously, the presence of SH basis in the process of SD has been proven to be of great importance. From a signal processing perspective, high-order SH basis is needed if we want to represent or reconstruct crossing fibers accurately with really small separated angles. However, the higher harmonic components are more sensitive to noise. Considering numerical difficulties, typically spherical harmonic up to the order of eight is used, which limits their capability in reliably resolving fiber crossing of small angles [40]. An inverse relationship exists between high frequency term and angular resolution. Thus, we cannot obtain the highest resolution and the best resistance to noise simultaneously.

On account of the above defects of SH basis, we use a new double-lobe basis function to build an over-completed dictionary basis. In this work, a set of over-completed orientation distribution basis {d(v, ul)|l = 1,…L} with discrete direction sets v ∈ RJ and positional direction sets u ∈ RL are introduced to represent fiber architecture in a voxel. The basis functions are uniformly distributed in unit sphere, thereby creating a predefined fODF field. A linear weighted combination of basis can be represented as ϕ = [d(v, u1),…, d(v, uL)]. By introducing an over-completed dictionary with cardinality L which is larger than unit sampling direction vectors J, we can construct a wide-ranged basis to map the fODF. Generally, fODF can be sparsely represented by the dictionary. Hence, most of the coefficients c are zero. The novel basis function quoted by [41] is proposed to establish the over-completed dictionary:(5)where ϑv, ul represent intersected angles between v and ul, and the other parameters, κ1, κ2 and τ are used to normalize the novel basis function. Detailed interpretation and specific parameters setting are described in [41]. Thus, we can obtain a novel dictionary basis which avoids high order’s truncation of SH function and guarantees the sparsity of coefficient at a certain extent(Fig 1). To make dictionary basis be applicable to above isotropic model, we have to make some deformations on the dictionary. The final dictionary basis can be represented as , where 1 represent J × 1 column vectors composed of 1.

[Figure omitted. See PDF.]

Fig 1. Weighting distributions.

Distributions of weighting estimated by solving the Least Norm method with Spherical harmonics (SH) and Dictionary basis (DB).

https://doi.org/10.1371/journal.pone.0168864.g001

A new Richardson-Lucy model

RL model is usually used in the field of astronomical imaging. This method has two advantages: the one is that it avoids the appearance of negative values in solutions because it satisfies non-negativity constraint of solution inherently; the other is that it well controls the instabilities in the process of solving and reduces the presence of noise artifacts in the solution for its robustness to noise [19]. Thus, the RL model has already been prevalent to settle the problem of fiber imaging in neurosciences field, as originally proposed by Daube-Whitherspoon and Muehllehner in [42].

Richardson-Lucy model with dictionary basis.

The RL model, also known as expectation maximization (EM) algorithm, follows a statistical Bayesian approach to deconvolution problem and implements an iterative estimation scheme for approximating the solutions of a maximum-likelihood problem in the case of different noises [19]. Therein, to establish a necessary foundation for the presentation and development of the proposed method, a brief overview of RL model is provided firstly. Like common approach of image restoration uses a probabilistic framework: given a sampling degraded image S, we can obtain the best image (actually is the fODF) when maximizing the probability of sampling image S. The probability obeys Bayes’ rule: . The magnitude signal of MR data is considered as Rician distribution [43], the likelihood is then:(6)where I0 is the modified Bessel functions of the first kind of zero order, denotes the true magnitude signal intensity in the absence of noise, S is measured value of each voxel with noise, and σ2 is the variance of noise. When the Rician distribution is acquired with large SNR(i.e.,), the process is better known as Gaussian approximation [19, 44, 45].

(7)

Considering the premise of hypothesis that the sampling signal submitted to spherical deconvolution, we optimize the Eq (7). The RL model finds f from the observation S, knowing response function by maximizing the likelihood distribution. The result can be derived by minimizing the function . We suppose that noise is independent from one voxel to another. When consider the whole brain region, the log-likelihood becomes a summation of the likelihood of all voxels. The multiplicative-type algorithm is equivalent to minimize J1(f) given by(8)where x is the voxel index. Given that the function J1(f) is convex in f, looking for a minimum is equal to searching for a zero value of gradient of J1(f). We set the derivative of J1(f) with respect to f to be zero and get . There are some mathematical operations which can be founded in [46]. Using the dictionary basis Φ to represent the fiber orientation f = Φc, we get:(9)

Richardson and Lucy suggested a multiplicative iterative method to solve Eq (9) (10)

Regularization with coefficient of fiber orientation.

For obvious reasons, the operation of spherical deconvolution is a NP-hard problem. To render the reconstruction perfectly and stably, we use regularization on the coefficient of fiber orientation, such as, total variation (TV) and sparse regularization. Putting a priori information on the coefficient of fiber orientation seems reasonable.

One of such information is spatial consistency. Despite many advantages of RL model, the fiber detail and noise interference are contradictory during the RL iteration process. This problem is generic for all maximum likelihood techniques because we usually want to attempt to fit the data as closely as possible. Thus, a trade-off exists between quality of image and the degree of noise interference when using RL method. In the intravoxel fODF field, voxels within a small neighbourhood usually consist of similar signals. Thus, the fODF derived from voxels ought to have a correlation in spatial structure. The advantages of using TV regularization are that it reserves the similarity of coefficient and avoids noise amplification by smoothing to certain extent. Here, we introduce TV constraint on the coefficient of fODF of the entire brain image to solve the above problem by adding energy function , defined as:(11)where λTV is the TV regularization parameter. Regularization is conducted in the entire field along each special gradient direction, which can be seen as I + 1 brain images. Although the images and cI+1 have different statistical properties, regularization processing in the neighbouring voxel is not prevented. The spatial dependence introduced by TV function promotes smooth solutions in homogeneous regions and prevents the solution from having oscillations. However the process of regulation will allow the solution to have sharp discontinuities [47, 48], we need to increase another constraints.

Sparse reconstruction method is broadly applied in the field of digital image processing. The sparsity constraint of the coefficient of fODF and the sparse recovery process lead us to estimate a sharp fODF from limited acquisitions. Notably, fiber orientation representation in the proposed basis is indeed sparse. The true distribution of fiber orientation can be considered sparse with the assumption that only a small number of elements of fODF are non-zero physically [49]. However, the introduction of TV regularization induces excessive smoothness between neighbouring voxels. To ensure each fiber sparse, the sparsity constraint is often added to fODF in spherical deconvolution problem. We make full use of ℓ1 norm to ensure the sparsity of coefficient in neighboring voxels. Here, we introduce the energy function of sparsity term , defined as(12)where λℓ1 is the sparse regularization parameter. The two regularization terms based on maximum likelihood estimation can get the derivatives of Jreg with respect to c, which can be expressed as and , respectively, where div and ∇ stand for divergence and differentiation, and x is voxel index indicating that regularization is conducted between voxels. The term |∇c| is replaced by its approximate value , where ε is a small positive constant [47]. The total energy function is known as(13)

We minimize Eq (13) using multiplicative gradient-based algorithm (or equivalently using EM algorithm for penalized criterion of Eq (13) and obtain the final result defined as(14)where c(k) is the estimated coefficient of fiber orientation, which is a ((L+1) × 1) dimension column vector at iteration k at voxel x, and and TV(k) are the ℓ1 and TV regularization vector at iteration k. The element at different gradient positions i of ℓ1 regularized vector is computed as:(15)

The element at different gradient position i of TV regularized vector is computed as:(16)

Numerically, we notice that the regularization parameter should be neither too small nor too large. In the simulated experiments, we will discuss the selection of regularization parameters.

Experiments and Results

Experimental data

Simulated data.

Datastes are generated assuming axially symmetric diffusion tensor profiles for each fiber population (MD = 0.7 × 10−3mm2/s) using a typical 81 directions sampling scheme [50]. To study the effect of each parameter separately in simulations, only one parameter at a time is varied. Details of these simulated datasets are provided in the following sections.

* Simulated data1: To guarantee the impartiality of comparative methods, We build the following simulated dataset which is the same with the data in [23], so that we can contrast iCSD method directly. Two crossing fibers are constructed, assuming the angle of crossing fiber is 70°, with varying PVEs values ranging from 0.1 to 1 (with a step of 0.1) and with different b-values of 1000 and 3000. The other dataset also reconstructs two crossing fibers, with varying crossing angles of fiber ranging from 40° to 90° (with a step of 10°) and with 50% isotropic signal. Complex Gaussian noise is added to obtain noisy signals with SNR = 20.

* Simulated data2: We create the synthetic data with two crossing fibers and different parameters which determine the imaging quality. Each simulated dataset is composed by 11 times 11 voxels whose fraction of isotropy is varied from 0.1 to 1 with intervals of 0.1 along x-axis, and SNR is changed from 10 to 30 with intervals of 2 along y-axis. The dataset is used to prove the validity of iRL to solve the PVEs under the condition of PVEs and noise changed. The representative angles are 40° and 90° between fiber populations in configurations, and the diffusion weighting b = 3000s/mm2.

IEEE international Symposium on Biomedical Imaging (ISBI) challenge phantom data.

The second simulated dataset coming from the ISBI 2013 Reconstruction Challenge is acquired from an open-source software library (http://hardi.ep.ch/static/events/2013-ISBI/), which creates realistic phantoms in structural and diffusion MRI. The synthetic datasets consist of 27 simulated ground truths, including branching, kissing, and crossing structures with angles between 30° and 90°. The dataset contains 64 gradient directions with b = 3000s/mm2 at SNR = 10, SNR = 20 and SNR = 30. The fODF mapping is color-coded by the standard DTI colour scheme (red: left-right; green: front-back; and blue: up-down).

In vivo human brain data.

Evaluation is performed using real human data which is published on Dipy (http://nipy.org/dipy/). The whole brain is covered with contiguous 2mm slices with an in-plane resolution of 2 × 2mm2. For preprocessing of diffusion data, we use MRIcron and SPM8 toolbox. First, the DICOM images sets (.dcm) are split into NIfTI (.nii), gradient sequence (.bvecs), and sensitive coefficient (.bvals) datasets using MRIcron software, where the NIfTI dataset contains scanned sequence corresponding to the gradient sequence. DW images are acquired along 150 uniformly distributed directions using b = 2000s/mm2 and a single b = 0s/mm2 (the size of the whole brain is 81 × 106 × 76).

Comparison metrics for phantom data

The performances of algorithms are quantified by comparing the obtained reconstructions with ground-truth. We adopt some of evaluation metrics widely used in the literatures [51–53].

* Average angular error (AAE): We compute the deviation between estimated fiber orientation and ground truth [54]:(17)where εx is the “ground truth” and is estimated fiber orientation, Ω is the local region used to compute angular error. we obtain one or more significant peaks of fODF (the number of peaks defined as nP) in each voxel x ∈ Ω, sum angular error of all peaks and finally get the average angular error. These operations are repeated about 100 times.

* Average probability of false direction (APFD): APFD is used to evaluate the probability of false directions compared with real fiber number inside a voxel x. The ratio of false positive (r+) and ratio of false negative (r−) are defined as(18)

In a region Ω, and denote the over-estimated and under-estimated number of fibers inside a voxel compared to ground truth.

* Fractional anisotropy (FA): The FA characterizes the degree of “out-of-roundness” of diffusion ellipsoid. It measures the fraction of total magnitude of diffusion tensor that is anisotropic(19)where λ1, λ2, λ3 are the eigenvalues provided by diffusion tensor, which is one of the most important rotationally invariant quantitative scalar parameters. is the arithmetic mean of the three eigenvalues.

* Generalized fractional anisotropy (GFA): Scalar measures on the fODF are useful in defining tissue contrast, performing statistical analyses, or summarizing the geometric properties of fODF. We define the scalar measures GFA as(20)where n is the number of fODF, std is the standard deviation, rms is the root-mean-square, and is the mean of the ODF. The GFA metric is automatically normalized to [0, 1].

* Generalized relative anisotropy (GRA): The GRA scalar represents a measurement of deviation from the isotropic state of the fODF of each voxel:(21)

It’s worth noting that peaks in clusters that are less than half of the crossing angle (with an upper limit of 35 degrees) from the true orientations are considered correct peaks.

Implementation details

All experiments of the proposed method are conducted on Inter(R)@2.4 GHz (48 G RAM, 64 bit). For measured signal, the obtained mask image is down-sampled to the dimensions of dMRI. Mask analysis is conducted on DSI Studio 1 (http://www.dsi-studio.labsolver.org). For the dictionary basis, the dimension of coefficients and the basis vectors are the same, representing the related percentage of each dictionary basis. For the positional direction sets u of dictionary basis, a tessellation scheme is distributed evenly on 321 points on a hemisphere and is generated by the subdivision of the face of an icosahedron. By avoiding repeated sampling, the discrete direction sets v are made to be identical with u. To perfectly reconstruct the fODF, the reconstructed dictionary basis is designed using a symmetric sphere with 10 242 vertices from Dipy (http://dipy.org/), which is an array of 10 242 fODF values corresponding to the vertices of sphere. To ensure the applicability of in vivo data, the two RFs in vivo data are acquired according to typical value of diffusion tensor signals in the corpus callosum and cortex respectively [17, 23]. We choose 50 voxels with the highest FA and use the average of signals whose principal eigenvectors are aligned along z-axis to acquire the anisotropic RF. Identically, we choose 50 voxels with the lowest FA and use the average of signals to acquire the isotropic RF.

We compare the proposed method iRL with the other state-of-the-art methods on simulated phantom and real data. The alternative approaches include Recursive calibration constrained spherical deconvolution (RC-CSD) [55], Sparse Fascicle Model (SFM) [56], damped Richardson Lucy (dRL) [26], information constrained spherical deconvolution (iCSD) [23] and Multi-shell multi-tissue constrained spherical deconvolution(MSMT-CSD) [24]. RC-CSD is an improvement of SD, which provides an accurately calibrated RF. SFM treats each MRI voxel as two types of compartments, non-oriented tissues and oriented fascicles considering the PVEs, which is implemented using Dipy (http://nipy.org/dipy/index.html) publish library [57]. The dRL is aiming at reducing isotropic background effects in spherical deconvolution, which is implemented using a software package provided in (http://neuroimagen.es/webs/hardi-tools/). The iCSD improves the estimation of fODF by modifying the RF to account for non-WM PVEs locally. MSMT-CSD uses CSD approach to estimate a multi-tissue ODF and implements in MRtrix (http://www.mrtrix.org/) [58]. It’s worth mentioning that MSMT-CSD can reconstruct brain fibers using single shell data, but the function of separating different tissues can not work well. The number of iterations of each method is set to 200 times. The related parameters used in compared methods are set to their optimal values according to the reference documents. For dRL algorithm, η acts as a threshold parameter and controls the damped amplitude of FOF, which is set to η = 0.08.

Results

Optimal regularization parameter.

The new deconvolution algorithm with TV and ℓ1 regularization has shown good imaging result with the elaborately chosen regularization parameters. The choice of good parameters value plays a crucial role in imaging results when using iRL. Thus, the first step of our experiment is to study if and how c estimation is influenced by setting different regularization parameters and by choosing different numbers of algorithm iterations during the process of our algorithm. To obtain the best regularization parameters and the number of iterations, we use different parameters to image the ISBI data with SNRs of 10, 20 and 30, and identify the quantitative index to evaluate image quality. To select regularization parameters, we use the AAE to be the quantitative index (Fig 2).

[Figure omitted. See PDF.]

Fig 2. Simulated results of parameter selection.

Average angular error using different regularization parameters in ISBI data with different SNRs.

https://doi.org/10.1371/journal.pone.0168864.g002

We performe 100 repetitions with simulated data. We find that the ℓ1 regularization parameter affects the angular resolution of imaging fiber and the TV regularization plays a vital role in resisting noise. We need only to increase the value of the TV or ℓ1 regularization parameters to improve the quality of imaging when the signals have low SNR or small angle, respectively. From Fig 2, the best regularization parameters can be set to λℓ1 = 0.01 and λTV = 0.5. The RL algorithm has certain superiority in resisting noise, but when the SNR is low, as shown in Fig 2, the imaging results are unsatisfactory and have a relatively large angular error.

The RL algorithm is known to have the property of ‘semi-convergence’ [59], i.e., the solution initially converges to the true value and then diverges as iterations proceed [19]. We choose 200 as the maximum iteration numbers to prevent noise amplification and generation of artifacts.

Simulated data in the presence of isotropic diffusion.

We use different simulated datasets to verify the effectiveness of iRL. Comparative tests are conducted by four kinds of methods. This experiment is used to verify the ability of imaging the signal with different volume fractions of isotropic signal (Fig 3). The other simulated datasets are generated in the same way, excepting that the diffusion weighting b is changed (Fig 4). We perform 100 repetitions with the simulated datasets that are generated randomly.

[Figure omitted. See PDF.]

Fig 3. Comparison of simulated results.

AAE, False peaks and Correct peaks for different proportions of anisotropic signal (diffusion weighting 3000s/mm2, angle 70°, and SNR 20).

https://doi.org/10.1371/journal.pone.0168864.g003

[Figure omitted. See PDF.]

Fig 4. Comparison of simulated results.

AAE, False peaks and Correct peaks for different proportions of anisotropic signal (diffusion weighting 1000s/mm2, angle 70°, and SNR 20).

https://doi.org/10.1371/journal.pone.0168864.g004

Compared with the low b value dataset, the high b value dataset shows fODF with a partial increase in angular resolution. However, no change exists in angular resolution for the signal with low proportion of anisotropy. Figs 3 and 4 show that when the signal has high proportion of anisotropy, the imagings of five kinds of methods are all accurate. The iCSD and iRL have a relative better angular resolution and less numbers of false peaks. When high isotropy exists in the simulated signal, iRL is advantageous over the other four kinds of methods in the aspect of angular resolution. Regardless of signal composition, iRL has the best and smallest angular resolution.

We perform simulated experiments to investigate the simulated datasets with different fiber crossing angles. We utilize five methods to image the above simulated signals respectively. This experiment is used to verify the ability of imaging the signal with different crossing angles (the results are shown in Fig 5). We also perform 100 repetitions with the simulated datasets which are generated randomly.

[Figure omitted. See PDF.]

Fig 5. Comparison of simulated results.

AAE, false peaks and correct peaks for different crossing angles(with 50% isotropic signal, diffusion weighting 3000s/mm2, and SNR 20).

https://doi.org/10.1371/journal.pone.0168864.g005

The five methods are all becoming more effective as the crossing angles increasing. In our method, the quantitative indexes of AAE and false peaks is lower for all angles and the precision is improved remarkably for angles larger than 50° (Fig 5). It’s worth mentioning that the 40° crossing angle could be identified with 50% PVEs using iRL.

We also perform simulated experiments to investigate simulated datasets with different PVEs and SNRs and utilize five methods to image the above simulated signals respectively (the results are shown in Fig 6). To verify the effectiveness of our method in aspect of the new isotropic quantitative index, we conduct the signal of simulated data2 and the imaging result is mapped to quantitative indexes, FA and GFA included.

[Figure omitted. See PDF.]

Fig 6. Comparison of the simulated results.

FA and GFA with RC-CSD(a), SFM (b), dRL (c), iCSD (d) and iRL (e). Figures a1-f1 are the imagings for 40° cross-angle. Figures a2-f2 are the imagings for 90° cross-angle. Figures f1-f2 are the new index Piso.

https://doi.org/10.1371/journal.pone.0168864.g006

In the case of anisotropy and SNR increased, the upper-left corner of each figure has the poorest simulated signal, and the lower-right corner of each figure has the best simulated signal. In Fig 6, we compare FA, GFA, and our new quantitative index Piso (because the quantification of Piso is the extent of isotropy, which is contrary to FA and GFA. For comparison, we make Piso a reverse imaging) using different methods (Figures a1-e1) and different crossing angles (Figures a1-a2). For 40° cross fibers, no significant difference is observed. For 90° cross fibers, the quantitative index FA has an obvious deficiency in which the degree of anisotropy is lower than the normal levels. However, the quantitative indexes GFA and Piso have a correct indication. Considering both experiments, Piso has better implementation in low anisotropy.

fODF estimation for ISBI data.

We compare several different methods using the authoritative ISBI simulated experiment data. Fig 7 compares the reconstructed fODF. We observe that the fODF estimations of each voxel are relatively independent and prone to noise. The fiber orientations reconstructed by standard RC-CSD, SFM, dRL, and MSMT-CSD methods always lack important information on fiber crossing.

[Figure omitted. See PDF.]

Fig 7. Visualization of the fODF reconstructed from ISBI dataset with HARDI data.

Depicted fODF profiles correspond to estimations from the RC-CSD (a), SFM (b), dRL (c), MSMT-CSD (d) and our method iRL (e).

https://doi.org/10.1371/journal.pone.0168864.g007

In the marked regions in Fig 7, the crossing angles are very small. The iRL can separate this part of crossing, but the results are imperfect. In the crossing fiber case, performances are assessed according to two criteria: (1) the effect of miscalibration on angular resolution, and (2) the over-estimated and under-estimated number of fibers. Fig 8 shows that iRL produces fewer angular errors. About the overestimation of false peaks, iRL has a better result when compared with RC-CSD and SFM. There is a better result about underestimation of false peaks when compared with dRL and MSMT-CSD. It’s mentioning that the iRL has fewer total numbers of false peaks than the other five methods, regardless of SNR.

[Figure omitted. See PDF.]

Fig 8. Quantification of the reconstruction accuracy.

The results of RC-CSD, SFM, dRL, MSMT-CSD, and iRL in terms of AAE, n+, and n− using ISBI data.

https://doi.org/10.1371/journal.pone.0168864.g008

fODF estimation for human data.

Evaluation is performed using real human data acquired on public datasets (http://nipy.org/dipy/). We select two representative areas, one of the areas contains multiple functional areas of the brain, such as the cortex and CSF (i.e, containing possible isotropic compartment).

Fig 9 compares the intravoxel fiber architecture estimated by five different methods on the human datasets. In the posterior thalamic radiation (refer to Human Brain in ICBM-152 Space) region (marked with a yellow box in Fig 9), the situation of fiber crossing is complex, containing single fiber and multiple fiber crossings. The iRL has a good imaging of multiple fiber crossing trends. The other methods always lack of some fiber directions. The same results can be seen in Fig 10. In addition, in the posterior thalamic radiation region, the isotropic signal is stronger, and the compared results are more obvious. In particular, the fibers (red ellipses) in the superior temporal gyrus WM (STG-WM) and the middle temporal gyrus WM (MTG-WM) regions are well represented by iRL.

[Figure omitted. See PDF.]

Fig 9. Visualization of fODFs reconstructed from real data.

Depicted fODF profiles correspond to the estimations from RC-CSD (a), SFM (b), dRL (c), MSMT-CSD (d) and iRL (e). The background images are fractional anisotropy images computed from each reconstruction.

https://doi.org/10.1371/journal.pone.0168864.g009

[Figure omitted. See PDF.]

Fig 10. Visualization of fODFs reconstructed from real data.

Depicted fODF profiles correspond to the estimations from RC-CSD (a), SFM (b), dRL (c), MSMT-CSD (d) and iRL (e). The background images are the fractional anisotropy images computed from each reconstruction.

https://doi.org/10.1371/journal.pone.0168864.g010

The quantitative indexes of GM and WM are carried out in above areas. We use Piso to quantify the difference between WM and GM in brain regions by using different indexes, including FA, GFA, and GRA. The three indexes are well-known and used in various occasions to describe the strength of anisotropic diffusion.

The degree of diffusion anisotropy is severely underestimated using the indexes calculated by diffusion coefficients acquired in fiber orientations. Some researchers present that water diffusivity in the directions parallel to the fiber is almost 10 times higher than the average diffusivity in directions perpendicular to them [9]. The marked area where the fibers have vertical distribution. The anisotropy is actually very strong, whereas the figure of FA (Fig 11b1) shows a strong isotropy. The figures of Piso (Fig 11a), GFA (Fig 11b2), and GRA (Fig 11b3) show similar results on anisotropy.

[Figure omitted. See PDF.]

Fig 11. Display of diffusion degree using four methods.

(a): Piso quantifies the intensity of isotropic signal, and (b1-b3) quantify the intensity of anisotropic signal.

https://doi.org/10.1371/journal.pone.0168864.g011

Conclusions and Discussion

We focus on PVEs in the reconstruction of fiber configuration, which rarely elicit interest of researchers. PVEs are some of the greatest obstacles in improving the accuracy of fiber imaging. We usually utilize the anisotropic signal to reconstruct fiber orientation, which is affected by the isotropic signal. Only by removing the isotropic signal from DW signal, can we obtain the best imaging results, as we have done in this paper. The contribution of our approach is that we initially propose a method based on the local maximum likelihood estimation to isolate the isotropic from DW signal in entire regions included in both non-WM and WM by rebuilding RF and fODF used to estimate the coefficients of fODF to account for tissues composition. At the same time, the separated parts can be used to quantify the degree of isotropic signal in each individual voxel. Secondly, the application of dictionary basis and RL model successfully solves the ill-posed problem and ringing effect. Finally, the spatial regularization of FOD is approximated by combining TV and ℓ1 norms that stabilize the deconvolution problem and promote sparsity in the solution. We also compare the performances of proposed method with several state-of-the-art algorithms on synthetic data and human brain datasets. Results show significant improvement over contrastive methods in its ability to reduce false positive fiber orientations and preserve angular resolution on both simulated and in vivo datasets.

Some of non-WM PVEs are due to the reduced SNR of WM compartment, which cannot be recovered, and the rest of effects are due to mostly isotropic diffusion from non-WM tissue [23]. In this paper, we extend PVEs’ influence, including the isotropic diffusion in WM and the increase in isotropy caused by complex fiber directions. By isolating the isotropic signal, the imaging results significantly improve, especially on the AAE, throughout the whole brain. From the Fig 8, we control the AAE within 8° using open ISBI data.

Simulated results show that with the reduction of isotropic signal, the AAE significantly increases. As regard 50° of crossing fiber, although the proportion of isotropic signal is as low as 0.1, the AAE remains within 30°. This is a complicated process because the imaging result is affected by many parameters, such as b value, regularization parameters, iterations, and so on. For different datasets, we should adjust the corresponding parameters to obtain the best imaging results. Notably, a lower b value leads to poorer imaging. We can find another defect, i.e., the decrease in fiber quantity is more outstanding than the overestimation of fiber in the simulated data. This problem is inherent in the method related to RL, which will be our concern in a future study.

Real experimental results indicate that iRL efficiently improves the ability of resolving crossing fibers in regions with high PVEs, whereas in high anisotropy regions, iRL and others produce roughly identical results. In the region of the internal capsule and the corpus callosum, the tracts have relatively larger amplitude, which is particularly useful in connectomics. Given the abandonment of least squares and spherical harmonic function, the spurious fODF peaks (consistent with well-known ringing artefacts) have a prominent reduction on Figs 9 and 10. The comparisons of the tract density image between iRL and others show increased tract density in the main WM tracts and decreased tract density in non-WM region, which are useful for fiber tracking.

Some open areas of researches exist in iRL. Firstly, for the two different diffusion models, different choices exist for regularization parameters. Considering the different diffusion regions, the strength of regularization should be discrepant. Secondly, a calibrated RF must be used to further reduce spurious peaks. Fortunately, the methods based on RL have a low over-all sensitivity to miscalibration. Thirdly, this method has potential to considerably reduce gradient directions, indicating a clinically feasible acquisition time. Thus, the application of this method is significant in clinical studies in the future. Finally, the assumed unimodal Gaussian diffusion model does not apply to MRI measurements, which are completely proven to be Rician distribution model [60]. These existing problems will be studied in our future work.

Acknowledgments

I thank my parents for allowing me to realize my own potential. All supports they have provided me over these years are the greatest gift I have ever received. Moreover, I need to thank my tutor for providing me good academic atmosphere and answering all my questions.

Author Contributions

1. Conceptualization: TTX YW YJF.

2. Data curation: TTX YW JZ.

3. Formal analysis: TTX QRZ YW.

4. Funding acquisition: YJF QCZG.

5. Investigation: YJF QCZG.

6. Methodology: TTX YW JZ.

7. Project administration: TTX QRZ.

8. Resources: TTX JZ YW.

9. Software: TTX YW QRZ.

10. Supervision: YJF QCZG.

11. Validation: TTX QRZ JZH.

12. Visualization: TTX YW JZH.

13. Writing – original draft: TTX YW JZ.

14. Writing – review & editing: TTX JZH YJF.

Citation: Xu T, Feng Y, Wu Y, Zeng Q, Zhang J, He J, et al. (2017) A Novel Richardson-Lucy Model with Dictionary Basis and Spatial Regularization for Isolating Isotropic Signals. PLoS ONE 12(1): e0168864. https://doi.org/10.1371/journal.pone.0168864

References

1. Ganguly D, Chakraborty S, Kim Th. A cognitive study on medical imaging. International Journal of Bio-Science and Bio-Technology. 2010;2(3):1–18.

2. Stejskal EO, Tanner JE. Spin Diffusion Measurements: Spin Echoes in the Presence of a Time-Dependent Field Gradient. Journal of Chemical Physics. 1965;42(1):288–292.

3. Dell’Acqua F, Catani M. Structural human brain networks: hot topics in diffusion tractography. Current Opinion in Neurology. 2012;25(25):375–83. pmid:22766720

4. Christiaens D, Reisert M, Dhollander T, Sunaert S, Suetens P, Maes F. Global tractography of multi-shell diffusion-weighted imaging data using a multi-tissue model. Neuroimage. 2015;123:89101.

5. Aganj I, Lenglet C, Sapiro G, Yacoub E, Ugurbil K, Harel N. Reconstruction of the orientation distribution function in single-and multiple-shell q-ball imaging within constant solid angle. International Journal of Bio-Science and Bio-Technology. Magnetic Resonance in Medicine. 2010;64(2):554–566.

6. Van Essen DC, Smith SM, Barch DM, Behrens TE, Yacoub E, Ugurbil K, et al. The WU-Minn human connectome project: an overview. Neuroimage. 2013;80:62–79. pmid:23684880

7. Basser PJ, Mattiello J, Lebihan D. Estimation of the effective self-diffusion tensor from the NMR spin echo. Journal of Magnetic Resonance. 1994;103(3):247–54. pmid:8019776

8. Basser PJ, Pajevic S, Pierpaoli C, Duda J, Aldroubi A. In vivo fiber tractography using DTMRI data. Magnetic Resonance in Medicine. 2000;44(4):625–32. pmid:11025519

9. Pierpaoli C, Basser PJ. Toward a quantitative assessment of diffusion anisotropy. Magnetic Resonance in Medicine Official Journal of the Society of Magnetic Resonance in Medicine. 1996;36(6):893–906. pmid:8946355

10. Basser PJ, Pierpaoli C. Microstructural and Physiological Features of Tissues Elucidated by Quantitative-Diffusion-Tensor MRI. Journal of Magnetic Resonance. 2011;213(2):560–70. pmid:22152371

11. Hess CP, Mukherjee P, Han ET, Xu D, Vigneron DB. Q-ball reconstruction of multimodal fiber orientations using the spherical harmonic basis. Magnetic Resonance in Medicine. 2006;56(1):104–117. pmid:16755539

12. Wiegell MR, Larsson HB, Wedeen VJ. Fiber crossing in human brain depicted with diffusion tensor MR imaging. Radiology. 2000;217(3):897–903. pmid:11110960

13. Frank LR. Anisotropy in high angular resolution diffusion-weighted MRI. Magnetic Resonance in Medicine Official Journal of the Society of Magnetic Resonance in Medicine. 2001;45(6):935–939.

14. Descoteaux M, Deriche R, Le Bihan D, Mangin JF, Poupon C. Multiple q-shell diffusion propagator imaging. Medical image analysis. 2011;15(4):603–621. pmid:20685153

15. Tuch DS. Q-ball imaging. Magnetic Resonance in Medicine. 2004;52(6):1358–1372. pmid:15562495

16. Wedeen VJ, Hagmann P, Tseng WYI, Reese TG, Weisskoff RM. Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging. Magnetic Resonance in Medicine. 2005;54(6):1377–1386. pmid:16247738

17. Tournier JD, Calamante F, Gadian DG, Connelly A. Direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution. Neuroimage. 2004;23(3):1176–1185. pmid:15528117

18. Li-Wei K, Jyh-Horng C, Van Jay W, Wen-Yih Isaac T. Optimization of diffusion spectrum imaging and q-ball imaging on clinical MRI system. Neuroimage. 2008;41(1):7–18. pmid:18387822

19. Flavio D, Giovanna R, Paola S, Rafael Alonso C, Giuseppe S, Ferruccio F. A model-based deconvolution approach to solve fiber crossing in diffusion-weighted MR imaging. IEEE transactions on bio-medical engineering. 2007;54(3):462–472. pmid:17355058

20. Alexander AL, Hasan KM, Lazar M, Tsuruda JS, Parker DL. Analysis of partial volume effects in diffusion-tensor MRI. Magnetic Resonance in Medicine. 2001;45. pmid:11323803

21. Dell’Acqua F, Scifo P, Rizzo G, Scotti G, Fazio F. A modified damped richardson-lucy algorithm to improve the estimation of fiber orientations in spherical deconvolution. In: Proc. Intl. Soc. Mag. Reson. Med. vol. 16; 2008.

22. Claudia MB, Sarah H, Jones DK, Alexander L, Aggleton JP, O’Sullivan MJ. Temporal association tracts and the breakdown of episodic memory in mild cognitive impairment. Neurology. 2012;79(23):2233–40. pmid:23175726

23. Roine T, Jeurissen B, Perrone D, Aelterman J, Philips W, Leemans A, et al. Informed constrained spherical deconvolution (iCSD). Medical Image Analysis. 2015;24(1):269–281. pmid:25660002

24. Jeurissen B, Tournier JD, Dhollander T, Connelly A, Sijbers J. Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data. Neuroimage. 2014;103:411–426. pmid:25109526

25. Roine T, Jeurissen B, Perrone D, Aelterman J, Leemans A, Philips W, et al. Isotropic non-white matter partial volume effects in constrained spherical deconvolution. Frontiers in Neuroinformatics. 2014;8:28. pmid:24734018

26. Flavio D, Paola S, Giovanna R, Marco C, Andrew S, Giuseppe S, et al. A modified damped Richardson-Lucy algorithm to reduce isotropic background effects in spherical deconvolution. Neuroimage. 2010;49(2):1446–1458. pmid:19781650

27. Canales-Rodríguez EJ, Daducci A, Sotiropoulos SN, Caruyer E, Aja-Fernández S, Radua J, et al. Spherical deconvolution of multichannel diffusion MRI data with non-Gaussian noise models and total variation spatial regularization; 2014.

28. Bertero M, Boccacci P, Robberto M. Inversion method for the restoration of chopped and nodded images. Proceedings of SPIE—The International Society for Optical Engineering. 1998;406(2):765–772.

29. Yap PT, Shen D. Spatial Transformation of DWI Data Using Non-Negative Sparse Representation. IEEE Transactions on Medical Imaging. 2012;31(11):2035–49. pmid:22711770

30. Yap PT, Zhang Y, Shen D. Multi-Tissue Decomposition of Diffusion MRI Signals via L0 Sparse-Group Estimation. IEEE Transactions on Image Processing A Publication of the IEEE Signal Processing Society. 2016; p. 1-1.

31. Zhou Q. Spatially Regularized Reconstruction of Fibre Orientation Distributions in the Presence of Isotropic Diffusion. April 2014;.

32. Jones DK, Knsche TR, Turner R. White matter integrity, fiber count, and other fallacies: The do’s and don’ts of diffusion MRI. Neuroimage. 2013;73:239–254. pmid:22846632

33. Assemlal HE, Tschumperl D, Brun L, Siddiqi K. Recent advances in diffusion MRI modeling: Angular and radial reconstruction. Medical Image Analysis. 2011;15(4):369–96. pmid:21397549

34. Tournier J, Calamante F, Connelly A, et al. Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution. NeuroImage. 2007;35(4):1459–1472. pmid:17379540

35. Anderson AW. Measurement of fiber orientation distributions using high angular resolution diffusion imaging. Magnetic Resonance in Medicine Official Journal of the Society of Magnetic Resonance in Medicine. 2005;54(5):1194–1206. pmid:16161109

36. Alexander DC. Maximum entropy spherical deconvolution for diffusion MRI. In: Information Processing in Medical Imaging. Springer; 2005. p. 76–87.

37. Clark C, Le-Bihan D. Water diffusion compartmentation and anisotropy at high b values in the human brain. Magnetic Resonance in Medicine. 2000;44(6):852–859. pmid:11108621

38. Geert DG, Marleen V, Colline P, Alexander L, Marcel E, Darras VM, et al. Structural changes between seasons in the songbird auditory forebrain. Frontiers in Behavioral Neuroscience. 2013;7(7):1–13. pmid:19864568

39. Tran G, Shi Y. Fiber orientation and compartment parameter estimation from multi-shell diffusion imaging. Medical Imaging, IEEE Transactions on. 2015;34(11):2320–2332. pmid:25966471

40. Tran G, Shi Y. Adaptively constrained convex optimization for accurate fiber orientation estimation with high order spherical harmonics. In: Medical Image Computing and Computer-Assisted Intervention-MICCAI 2013. Springer; 2013. p.485–492.

41. Wu Y, Feng Y, Li F, Westin CF. Global consistency spatial model for fiber orientation distribution estimation. In: Biomedical Imaging (ISBI), 2015 IEEE 12th International Symposium on; 2015.

42. Daube-Witherspoon ME, Muehllehner G. An Iterative Image Space Reconstruction Algorthm Suitable for Volume ECT. Medical Imaging IEEE Transactions on. 1986;5(2):61–66. pmid:18243988

43. Basu S, Fletcher T, Whitaker R. Rician noise removal in diffusion tensor MRI. Springer Berlin Heidelberg; 2006.

44. Behrens TE, Berg HJ, Jbabdi S, Rushworth MF, Woolrich MW. Probabilistic diffusion tractography with multiple fibre orientations: What can we gain? Neuroimage. 2007;34(1):144–155. pmid:17070705

45. Gudbjartsson H, Patz S. The Rician distribution of noisy MRI data. Magnetic Resonance in Medicine. 1996;34(2):910–4.

46. Raj A, Hess C, Mukherjee P. Spatial HARDI: Improved visualization of complex white matter architecture with Bayesian spatial regularization. Neuroimage. 2011;54(1):396–409. pmid:20670684

47. Camprubio L. Spherical Deconvolution of Multichannel Diffusion MRI Data with Non-Gaussian Noise Models and Spatial Regularization. Plos One. 2014;10(10):320. pmid:26470024

48. Jun Zhang YFYWYLJHSZ Tiantian Xu. A Self-adaptive Local Feature Extraction Based Magnetic Resonance Imaging. In: Control and Decision Conference (CCDC), 2016 28th Chinese. IEEE; 2016.

49. Feng Y, Wu Y, Rathi Y, Westin CF. Sparse deconvolution of higher order tensor for fiber orientation distribution estimation. Artificial Intelligence in Medicine. 2015;65(3):229–238. pmid:26428956

50. Leemans A, Sijbers J, Verhoye M, Linden AVD, Dyck DV. Mathematical framework for simulating diffusion tensor MR neural fiber bundles. Magnetic Resonance in Medicine Official Journal of the Society of Magnetic Resonance in Medicine. 2005;53(4):944–53. pmid:15799061

51. Daducci A, Canales-Rodriguez EJ, Descoteaux M, Garyfallidis E. Quantitative Comparison of Reconstruction Methods for Intra-Voxel Fiber Recovery From Diffusion MRI. IEEE Transactions on Medical Imaging. 2014;33(2):384–99. pmid:24132007

52. Seunarine KK, Alexander DC. Chapter 4Multiple Fibers: Beyond the Diffusion Tensor. Diffusion Mri. 2009; p. 55–72.

53. Ct MA, Girard G, Bor A, Garyfallidis E, Houde JC, Descoteaux M. Tractometer: Towards validation of tractography pipelines. Medical Image Analysis. 2013;17(7):844–857. pmid:23706753

54. Canales-Rodriguez E, Melie-Garcia Y, Liturria Medina, Martinez-Montes E, Aleman-Gomez Y, Lin C. Inferring multiple maxima in intravoxel white matter fiber distribution. Magnetic Resonance in Medicine. 2008;60(3):616–30. pmid:18727080

55. Tax CMW, Jeurissen B, Vos SB, Viergever MA, Leemans A. Recursive calibration of the fiber response function for spherical deconvolution of diffusion MRI data. Neuroimage. 2014;86(2):67–80. pmid:23927905

56. Rokem A, Yeatman JD, Pestilli F, Kay KN, Mezer A, Van d WS, et al. Evaluating the Accuracy of Diffusion MRI Models in White Matter. Plos One. 2015;10(4). pmid:25879933

57. Garyfallidis E, Brett M, Amirbekian B, Rokem A, Van d WS, Descoteaux M, et al. Dipy, a library for the analysis of diffusion MRI data. Frontiers in Neuroinformatics. 2014;8(8):8. pmid:24600385

58. Tournier JD, Calamante F, Connelly A. MRtrix: Diffusion tractography in crossing fiber regions. International Journal of Imaging Systems & Technology. 2012;22(1):53–66.

59. Bertero M, Boccacci P. Image restoration methods for the Large Binocular Telescope (LBT). Astronomy and Astrophysics Supplement. 2000;147(2):323–333.

60. Wiestdaessl N, Prima S, Coup P, Morrissey SP, Barillot C. Rician Noise Removal by Non-Local Means Filtering for Low Signal-to-Noise Ratio MRI: Applications to DT-MRI. In: International Conference on Medical Image Computing and Computer-assisted Intervention; 2008. p. 171-9.