Spherical Deconvolution of Multichannel Diffusion MRI Data with Non-Gaussian Noise Models and Spatial Regularization

Spherical deconvolution (SD) methods are widely used to estimate the intra-voxel white-matter fiber orientations from diffusion MRI data. However, while some of these methods assume a zero-mean Gaussian distribution for the underlying noise, its real distribution is known to be non-Gaussian and to depend on many factors such as the number of coils and the methodology used to combine multichannel MRI signals. Indeed, the two prevailing methods for multichannel signal combination lead to noise patterns better described by Rician and noncentral Chi distributions. Here we develop a Robust and Unbiased Model-BAsed Spherical Deconvolution (RUMBA-SD) technique, intended to deal with realistic MRI noise, based on a Richardson-Lucy (RL) algorithm adapted to Rician and noncentral Chi likelihood models. To quantify the benefits of using proper noise models, RUMBA-SD was compared with dRL-SD, a well-established method based on the RL algorithm for Gaussian noise. Another aim of the study was to quantify the impact of including a total variation (TV) spatial regularization term in the estimation framework. To do this, we developed TV spatially-regularized versions of both RUMBA-SD and dRL-SD algorithms. The evaluation was performed by comparing various quality metrics on 132 three-dimensional synthetic phantoms involving different inter-fiber angles and volume fractions, which were contaminated with noise mimicking patterns generated by data processing in multichannel scanners. The results demonstrate that the inclusion of proper likelihood models leads to an increased ability to resolve fiber crossings with smaller inter-fiber angles and to better detect non-dominant fibers. The inclusion of TV regularization dramatically improved the resolution power of both techniques. The above findings were also verified in human brain data.


Introduction
After decades of developments in diffusion Magnetic Resonance Imaging (MRI), the successful implementation of a variety of advanced methods has shed light on the complex patterns of brain organization present at micro [1] and macroscopic scales [2][3][4]. Among these methods, Diffusion Tensor Imaging (DTI) [5] has become a classic in both clinical and research studies. DTI can deliver quantitative results, it may be easily implemented in any clinical MRI system and, thanks to its short acquisition time, it may be suitable for studying a wide range of brain diseases. Unfortunately, it is now well recognized that due to its simplistic assumptions, the DTI model does not adequately describe diffusion processes in areas of complex tissue organization, like in areas with kissing, branching or crossing fibers [6].
Spherical Deconvolution (SD) is a class of multi-compartment reconstruction technique that can be implemented using both parametric and nonparametric signal models [36][37][38][39][40][41][42][43][44][45][46][47][48][49]. SD methods have become very popular owing to their ability to resolve fiber crossings with small inter-fiber angles in datasets acquired within a clinically feasible scan time. This resolving power is driven by the fact that, as opposed to model-free techniques that estimate the diffusion Orientational Distribution Function (ODF), the output from SD is directly the fiber ODF itself.
Among the different SD algorithms, Constrained Spherical Deconvolution (CSD) [39,40] has been received with special interest due to its good performance and short computational time. In CSD, the average signal profile from white-matter regions of parallel fibers is first estimated, and afterwards, the fiber ODF is estimated by deconvolving the measured diffusion data in each voxel with this signal profile, which is also known as the single-fiber 'response function'.
More recently and as an alternative to CSD, a new SD method based on a damped Richardson-Lucy algorithm adapted to Gaussian noise (dRL-SD) has been proposed [37,42]. An extensive evaluation of both CSD and dRL-SD algorithms has revealed a superior ability to resolve low anisotropy crossing-fibers by CSD but a lower percentage of spurious fiber orientations and a lower over-all sensitivity to the selection of the response function by the dRL-SD approach [50]. This later feature is of great relevance since the assumption of a common response function for all brain tracts is a clear over-simplification of both methods, with the consequences of it minimized by the dRL-SD.
From an algorithmic perspective dRL-SD inherits the benefits of the standard RL deconvolution method applied with great success in diverse fields ranging from microscopy [51] to astronomy [52]. Remarkably, RL deconvolution is robust to the experimental noise and the obtained solution can be constrained to be non-negative without the need for including additional penalization functions in the estimation process. Moreover, from a modeling point of view, dRL-SD is implemented using an extended multi-compartment model that allows considering the partial volume effect in brain voxels with mixture of white matter (WM), gray matter (GM) and cerebrospinal fluid (CSF), a strategy that has been shown to be effective in reducing the occurrence of spurious fiber orientations [37].
However, in spite of the good properties of dRL-SD and other SD methods some methodological issues remain. These methods, to some extent, assume additivity and zero mean Gaussianity for the underlying noise and are potentially vulnerable to significant departures from such an assumption. Indeed, it is well known that the MRI noise is non-Gaussian and depends on many factors, including the number of coils in the scanner and the multichannel image combination method. Real experiments have shown that noise follows Rician [53] and noncentral Chi (nc-χ) distributions [54] evidencing the inappropriateness of the Gaussian model. This issue is especially relevant in diffusion MRI data where the high b-values required to enhance the angular contrast lead to extremely low signal-to-noise (SNR) ratios. A recent study [55] has shown that different multichannel image combination methods can changes the properties of the signal and can have an effect on fiber orientation estimation.
On the other hand, the standard reconstruction in SD, based on a voxel-by-voxel fiber ODF estimation, although reasonable it may not be optimal in a global sense as it does not take into account the underlying spatial continuity of the image. Recent research on the inclusion of spatial continuity into SD methods via regularization has yielded promising results [9,56,57]. Among these, spatially regularized SD methods based on Total Variation (TV) [9] are very appealing due to their outstanding ability to simultaneously smooth away noise in flat regions whilst preserving edges, and due to their robustness to high levels of noise [58].
This work has two main aims: (1) the study and quantification of the benefits of the adequate modelling of the noise distribution in the context of spherical deconvolution, and (2) the study and quantification of the effects of including a TV spatial regularization term in the proposed estimation framework.
To address the first objective we developed a new SD methodology which, following a more realistic view, deals with non-Gaussian noise models. Specifically, the estimation framework is based on a natural extension to the RL algorithm for Rician and nc-χ likelihood distributions. We had chosen the RL algorithm as a starting point for our work because this algorithm has proven to be highly efficient in diverse applications, and because the performance of the resulting method can be directly compared to the state-of-the-art dRL-SD method, which employs a nearly-equivalent SD estimation algorithm but based on a Gaussian noise model. The second aim was addressed by including TV regularization to the developed formulation. Moreover, for completeness we have extended also the dRL-SD method via the spatially-regularized proposed estimation.
To compare the relative performance between the SD methods based on Gaussian and non-Gaussian noise models, and their respective implementations including the TV regularization, the different algorithms were applied to several synthetic phantom datasets which had been contaminated with noise patterns mimicking the Rician and nc-χ noise distributions produced in multichannel scanners. To the best of our knowledge, this is the first evaluation of such methods in a scenario where Rician and nc-χ noise are explicitly created as a function of the number of coils, their spatial sensitivity maps, the correlation between coils, and the reconstruction methodology used to combine the multichannel signals. As a final analysis, the new method is also applied to real multichannel diffusion MRI data from a healthy subject.
Following this introduction there is a 'Theory' section providing an overview of the different topics relevant to the study and the derivation of the new SD reconstruction algorithm. Description about computer simulations, image acquisition strategies and metrics designed to evaluate the performance of the reconstructions is provided in the Materials and Methods section. Relevant findings are succinctly described in the Results section. Finally, main results, contributions and limitations of this work are addressed in the Discussion and Conclusions section.

Theory
This section contains a description of the forward/generative model used to relate the local diffusion process with the measured diffusion MRI data. It also provides a brief review of MRI noise models. Finally, the diffusion and MRI noise models are used to derive the new SD reconstruction algorithms.

Generative signal and fiber ODF model
The diffusion MRI signal measured for a given voxel can be expressed as the sum of the signals from each intra-voxel compartment. The term 'compartment' is defined as a homogeneous region in which the diffusion process possesses identical properties in magnitude and orientation throughout, and which is different to the diffusion processes occurring in other compartments. One example of this approach is the multi-tensor model that allows considering multiple WM parallel-fiber populations within the voxel. In this model the diffusion process taking place inside each compartment of parallel fibers is described by a second-order self-diffusion tensor [6].
In real brain data, in addition to the different WM compartments, voxels might also contain GM and CSF components. This issue was considered by [37], who extended the multi-tensor model by incorporating the possible contribution from these compartments. This is the generative multi-tissue signal model that will be used in the present study. In the absence of any source of noise, the resulting expression for the signal is: where M is the total number of WM parallel fiber bundles; f j denotes the volume fraction of the j th fiber-bundle compartment; f GM and f CSF are the volume fractions of the GM and CSF com- is the diffusion-sensitization factor (i.e., b-value) used in the acquisition scheme to measure the diffusion signal S i along the diffusion-sensitizing gradient unit vector v i , i = 1, . . ., N; D GM and D CSF are respectively the mean diffusivity coefficients in GM and CSF; S 0 is the signal amplitude in the absence of diffusionsensitization gradients (b i = 0); D j ¼ R T j AR j denotes the anisotropic diffusion tensor of the j th fiber-bundle, where R j is the rotation matrix that rotates a unit vector initially oriented along the x-axis toward the j th fiber orientation (θ j , ϕ j ) and A is a diagonal matrix containing information about the magnitude and anisotropy of the diffusion process inside that compartment: where λ 1 is the diffusivity along the j th fiber orientation, λ 2 and λ 3 are the diffusivities in the plane perpendicular to it. It is assumed that λ 1 >λ 2 %λ 3 . At each voxel, the measured diffusion signals S i for N different sampling parameters (i.e., v i and b i , i[1, . . ., N]) can be recast in matrix form as:  In the framework of model-based spherical deconvolution, H is created by specifying the diffusivities, which are chosen according to prior information, and by providing a dense discrete set of equidistant M-orientations O ¼ fðy j ; j Þ; j 2 ½1; . . .; Mg uniformly distributed on the unit sphere. Previous studies have used different sets of orientations, ranging from M = 129 [43] to M = 752 [42]. Then, the goal is to infer the volume fraction of all predefined oriented fibers, f, from the vector of measurements S and the 'dictionary' H of oriented basis signals. Under this reconstruction model, f can be interpreted to as the fiber ODF evaluated on the set O Matrix H is also known as the 'diffusion basis functions' [43], or the 'point spread function' [37][38][39] that blurs the fiber ODF to produce the observed measurements.
It should be noticed that solving the deconvolution problem given by Eq (3) is not simple because the resulting system of linear equations is ill-conditioned and ill-posed (i.e., there are more unknowns than measurements and some of the columns of H are highly correlated), which can lead to numerical instabilities and physically meaningless results (e.g., volume fractions with negative values). A common strategy to avoid such instabilities is to use robust algorithms that search for solutions compatible with the observed data but which also satisfy some additional constraints. Thus, in SD it is typical to estimate the fiber ODF by constraining it to be non-negative and symmetric around the origin (i.e., antipodal symmetry). As mentioned in the introduction, though, all these reconstruction algorithms may not be necessarily optimal when dealing with non-Gaussian noise, as it is the case for MRI noise.

MRI noise models
In conventional MRI systems, the data are measured using a single quadrature detector (i.e., coil with two orthogonal elements) that gives two signals which, for convenience, are treated as the real and imaginary parts of a complex number. The magnitude of this complex number (i.e. the square root of the sum of their squares) is commonly used because it avoids different kinds of MRI artifacts [53]. Given that the noise in the real and imaginary components follows a Gaussian distribution, the magnitude signal S i will follow a Rician distribution [53] with a probability function given by where S i denotes the true magnitude signal intensity in the absence of noise, σ 2 is the variance of the Gaussian noise in the real and imaginary components, I 0 is the modified Bessel function of first kind of order zero and u is the Heaviside step function that is equal to 0 for negative arguments and to 1 for non-negative arguments.
Modern clinical scanners are usually equipped with a set of 4 to 32 multiple phased-array coils, the signals of which can be combined following different strategies that, in turn, will give rise to different statistical properties for the noise [54]. One frequent strategy uses the spatial matched filter (SMF) approach linearly combining the complex signals of each coil and producing voxelwise complex signals [59]. Since the noise in the resulting real and imaginary components remains Gaussian a Rician distribution is expected in the final combined magnitude image. An alternative to the SMF is to create the composite magnitude image as the root of the sum-of-squares (SoS) of the complex signals of each coil. Under this approach the combined image follows a nc-χ distribution [60] given by, where n is the number of coils and I n-1 is the modified Bessel function of first kind of order n-1. This expression is strictly valid when the different coils produce uncorrelated noise with equal variance, and when noise correlation cannot be neglected it provides a good approximation if effective n eff and s 2 eff values are considered [61], with n eff being a non-integer number lower than the real number of coils and s 2 eff is higher than the real noise variance in each coil. A related SoS image combination method that increases the validity of Eq (5) is the covariance-weighted SoS. This method is equivalent to pre-whitening (i.e., decorrelate) the measured signals before applying the standard SoS image combination. The covariance-weighted SoS approach requires the estimation of the noise covariance matrix of the system which, in practice, may be carried out by digitizing noise from the coils in the absence of excitations [62].
It is important to note that there are additional factors that can change the noise characteristics described above, including the use of accelerated techniques based on under-sampling approaches such as those used in parallel MRI (pMRI) and partial Fourier, certain reconstruction filters in k-space, and some of the preprocessing steps conducted after image reconstruction.
Empirical data suggest that some of these factors do not substantially change the type of distribution of the noise. On the one hand, [54] investigated the effects of the type of filter in kspace, the number of receiving channels and the use of pMRI reconstruction techniques, and found that noise distributions always followed Rician and nc-χ distributions with a reasonable accuracy-although their standard deviations and effective number of receiver channels were altered when fast pMRI and subsequent SoS reconstructions were used. On the other hand, [55] showed real diffusion MRI data noise to also follow Rician and nc-χ noise distributions after a preprocessing that included motion and eddy currents corrections. Unfortunately, the combined effect of all factors has not, to the best of our knowledge, being studied. In this regard, a complete evaluation should include the study of the effects of additional data manipulation processes routinely applied in many clinical research studies, such as B0-unwarping due to magnetic field inhomogeneity and partial Fourier reconstructions. Although the latter has been investigated in terms of signal-to-noise ratio, its influence on the shape of the noise distribution remains unknown.
However, while it is impossible to ensure that Rician and nc-χ distributions are the optimal noise models for all possible strategies used for sampling, reconstructing and preprocessing diffusion MRI data, such models are flexible enough to adapt to deviations from the initial theoretical assumptions. Their parameterization in terms of spatial-dependent effective parameters (i.e., n eff ðx; y; zÞ; s 2 eff ðx; y; zÞ, as in [61,63,64]) allows characterizing the spatially varying nature of the noise observed in accelerated MRI reconstructed data, as well as the spatial correlation introduced by reconstruction algorithms, whilst preserving the good theoretical properties of the models with standard parameters, i.e. the null probability to obtain negative signals and the ability to characterize the signal-dependent non-linear bias of the data.

Spherical deconvolution of diffusion MRI data
Eq (5) based on either conventional (i.e., n, σ 2 ) or effective (i.e., n eff , s 2 eff ) parameters provides a very general MRI noise model, which includes the Rician distribution (given in Eq (4)) as a special case with n = 1. Consequently, if we derive the spherical deconvolution reconstruction corresponding to Eq (6) any particular solution of interest will become available.
Specifically, if we assume the linear model given by Eqs (1)-(3) the likelihood model for the vector of measurements S under a nc-χ distribution is where S i and S i ¼ ðHf Þ i are the measured and expected signal intensities for ith sampling parameters, respectively. 3.1 Unbiased and positive recovery: the multiplicative Richardson-Lucy algorithm for nc-χ noise. The maximum likelihood (ML) estimate in Eq (6) is obtained by differentiating its negative log-likelihood J(f) = −log p(S|H,f,σ 2 , n) with respect to f and equating the derivative to zero, which after some algebraic manipulations becomes where '' stands for the Hadamard component-wise product, and the division operators are applied component-wise to the vector's elements. Eq (7) is nonlinear in f and its solution can be obtained through a modified version of the expectation maximization technique, originally developed by Richardson and Lucy for a Poisson noise model [65,66] and known as the RL algorithm. When we applied this technique to nc-χ and Rician distributed noise it naturally led to the following iterative estimation formula: in which the solution calculated at the k th iteration step (f k ) gradually improves (i.e. its likelihood increases after each step) until a final, stationary solution f kþ1 f k ¼ 1, is reached. As shown Appendix A in S1 File, this formula can also be related to the RL algorithm for Gaussian noise, employed in the undamped RL-SD technique [42].
Under the absence of any prior knowledge about f, the initial estimate (f°) can be fixed to a non-negative constant density distribution [42]. In that case, the algorithm transforms a perfectly smooth initial estimate into sharper estimates, with sharpness increasing with the number of iterations. So, roughly speaking, the number of iterations can be considered as a regularization parameter controlling the angular smoothness of the final estimate. Notably, if f°is non-negative, the successive estimates remain non-negative as well, and the algorithm always produces reconstructions with positive elements. Moreover, as in [37,42] the estimation does not involve any matrix inversion, thus avoiding related numerical instabilities.
In order to evaluate Eq (8) an estimates 2 of σ 2 is required. Although obtaining it from a region-of-interest (ROI) is feasible [67] its accuracy may be compromised by systematic experimental issues such as ghosting artifacts, signal suppression by the scanner outside the brain, zero padding and by filters applied in the k-space. Moreover, with the use of fast parallel MRI sequences, where each coil records signals with partial coverage in the k-space, properties of the noise become spatially heterogeneous (i.e. they change from voxel to voxel across the image). While some authors have proposed alternatives to overcome these limitations [68] here we have estimated the noise variance at each voxel from the same data used to infer the fiber ODF.
Specifically, by minimizing the negative log-likelihood with respect to σ 2 we have obtained an iterative scheme analogous to Eq (8): where α k is the estimate of σ 2 at the k th iteration (starting with an arbitrarily initial estimate α°) and 1 N is a N×1 vector of ones. The resulting algorithm based on Eqs (8) and (9) is termed RUMBA-SD, which is the abbreviation of 'Robust and Unbiased Model-Based Spherical Deconvolution'. The spatially-regularized extension to this algorithm is described in the following section.
3.2 Towards a robust recovery: Total variation regularization. When considering the TV model [58] the maximum a posteriori (MAP) solution at voxel (x, y, z) is obtained by minimizing the augmented functional: where the first term is the negative log-likelihood defined in previous sections and the second term is the TV energy, defined as the sum of the absolute values of the first-order spatial derivative (i.e., gradient "r") of the fiber ODF components over the entire brain image, TVðf Þ ¼ X j jr½f 3D j j, evaluated at voxel (x, y, z); [f 3D ] j is a 3D image created in a way that each voxel contains the element at position j of their corresponding estimate vector f, and α TV is a parameter controlling the level of spatial regularization. Importantly, and in contrast to the previous ML estimate, now the solution at a given voxel is not independent from the solutions in other voxels. The spatial dependence introduced by the TV functional promotes smooth solutions in homogeneous regions (discourages the solution from having oscillations), yet it does allow the solution to have sharp discontinuities [58]. This property is highly relevant for SD because, while it promotes continuity and smoothness along individual tracts, it prevents partial volume contamination from adjacent tracts. In this work, the MAP estimate from Eq (10) is obtained using an iterative scheme similar to that proposed in [51], where the estimate at each iteration is calculated by the multiplication of two terms: the standard ML estimate, and the regularization term derived from the TV functional with the TV regularization vector R k at voxel (x, y, z), and at the k th iteration, computed element-by-element as where (R k ) j is the element j of vector R k and div is the divergence operator. In practice, to correct for potential singularities at r f k where ε is a small positive constant. Moreover, any negative value in R k is replaced by its absolute value to preserve the non-negativity of the estimated fiber ODF. Notice that by setting α TV = 0 the estimator in Eq (11) becomes equal to the unregularized version in Eq (8).
In the current implementation the simultaneous estimation of all the parameters is carried out via an alternating iterative scheme summarized in Table 1. Briefly, it minimizes the functional in Eq (10) with respect to the fiber ODF while assuming that σ 2 is known and fixed, and then it updates the noise variance using the new fiber ODF estimate. While for SMF-based data all the equations are evaluated using n = 1, for SoS-based data n is fixed to the real number of coils, or to the effective value n eff if provided. (But see Appendix B in S1 File) The regularization parameter α TV is adaptively adjusted at each iteration following the discrepancy principle. Specifically, it is selected to match the estimated variance [69] using two alternative strategies: (i) assuming a constant mean parameter over the entire brain image, α TV = E{α k+1 } (see Table 1), potentially increasing the precision and robustness of the estimator; or (ii) assuming a spatial dependent parameter, α TV (x, y, z) = α k+1 (x, y, z), which may be more appropriate in situations where a differential variance across the image is expected, like in accelerated pMRI based-data.
It should be remembered that the accuracy of the reconstruction for the SoS case depends on the variant used to combine the images. In this work it is assumed that the data is combined using the covariance-weighted SoS method. However, even if the available data were combined using the conventional SoS approach (i.e., without taking into account the noise correlation matrix among coil elements), the method could still provide a reasonable approximation (for more details see Appendix B in S1 File).   3)). This step would make sense when the fiber response signal used to create the dictionary matches the real signal from the compartments, whereas it may be omitted when the latter cannot be guaranteed. Notice that the original The evaluation of the ratio of modified Bessel functions of first kind involved in the updates of Eqs (11) and (9) is best computed by considering the ratio as a new composite function, and not by means of the simple evaluation of the ratio of the individual functions. Specifically, this ratio is computed here in terms of Perron continued fraction [70]. All the details are provided in Appendix C in S1 File.
Following a similar estimation framework, the TV regularization was also included into the dRL-SD method. All the relevant equations are provided in Appendix D in S1 File.

Materials and Methods 1 Synthetic fiber bundles with different inter-fiber angles
In order to test the resolving power of the methods as a function of the underlying inter-fiber angle, various synthetic phantoms including two fiber bundles were generated. The inter-fiber angles were gradually modified from 1 to 90 degrees applying one degree increases, eventually yielding 90 different phantoms with 50 x 50 x 50 voxels. The volume fractions of the two fiber bundles were assumed to be equal (f 1 = f 2 = 0.5) in the fiber crossing region.
The intra-voxel diffusion MRI signal was generated via the multi-tensor model [6] using N = 70 sampling orientations with constant b = 3000 s/mm 2 plus one additional image with b = 0 (i.e., S 0 = 1 was assumed in all voxels). The diffusion tensor diffusivities of both fiber groups were assumed to be identical and equal to λ 1 = 1.7 10 −3 mm 2 /s and λ 2 = λ 3 = 0.3 10 −3 mm 2 /s respectively.

Synthetic fiber bundles with different volume fractions
To test the ability of the different methods to detect non-dominant fibers, various synthetic phantoms containing two fiber bundles were generated. In these phantoms the inter-fiber angle was fixed to 70 degrees (an angle presumably detectable by all the methods) and the volume fraction of the non-dominant fiber bundle was gradually changed from 0.1 to 0.5 in 0.01steps, generating 41 different phantoms with 50 x 50 x 50 voxels each. The intra-voxel diffusion MRI signal was created using the same generative multi-tensor model, b-value and sampling orientations as in the previous section.

Synthetic "HARDI reconstruction challenge 2013" phantom
The reconstruction algorithms were also tested on the synthetic diffusion MRI phantom developed for the "HARDI Reconstruction Challenge 2013" Workshop, within the IEEE International Symposium on Biomedical Imaging (ISBI 2013). This phantom comprises a set of 27 fiber bundles with fibers of varying radii and geometry which connect different areas of a 3D image with 50 x 50 x 50 voxels. It contains a wide range of configurations including branching, crossing and kissing fibers, together with the presence of isotropic compartments mimicking the CSF contamination effects occurring near ventricles in real brain images.
The intra-voxel diffusion MRI signal was generated using N = 64 sampling points on a sphere in q-space with constant b = 3000 s/mm 2 , plus one additional image with b = 0. For pure GM and CSF voxels, signals were generated using two mono-exponential models: exp (−D GM b) and exp(−D CSF b) with D GM = 0.2 10 −3 mm 2 /s and D CSF = 1.7 10 −3 mm 2 /s. In voxels belonging to single-fiber WM bundles, the signal measured along the q-space unit direction q ¼ q=jqj was generated by a mixture of signals from intra-and extra-axonal compartments: f int s int ðq; v; t; L; RÞ þ f ext s ext ðq; v; b; l 1 ; l 2 Þ; where v denotes the local fiber orientation. Signal from the intra-axonal compartment S int was created following the theoretical model of a restricted diffusion process inside a cylinder of length L = 5 mm and radius R = 5 μm at the diffusion time τ = 20.8 s [19,71]. The extra-axonal signal S ext was generated using a diffusion tensor model with cylindrical symmetry (i.e., λ 1 = 1.7 10 −3 mm 2 /s, λ 2 = λ 3 = 0.2 10 −3 mm 2 /s). Mixture fractions were fixed to f int = 0.6 and f ext = 0.4. The noiseless dataset can be freely downloaded from the Webpage of the site: http://hardi.epfl.ch/static/events/2013_ISBI/.

Multichannel noise generation
The synthetic diffusion images from the above phantoms were contaminated with noise mimicking the SoS and SMF strategies used in scanners in order to combine multiple-coils signals. To that aim, the noisy complex-valued image measured from the kth coil was assumed to be equal to where C k is the relative sensitivity map [72] of kth coil, e R k $ Nð0; SÞ and e I k $ Nð0; SÞ are two different Gaussian noise realizations simulating the noise in the real and imaginary components with zero-mean value and covariance matrix S. For simplicity S was assumed to be given by where σ 2 is the noise variance of each coil and ρ indicates the correlation coefficient between any two coils. For the SoS reconstruction, magnitude images were generated as: where |S k | stands for the magnitude of S k . Notice that Eq (15) is the conventional SoS image combination and not the covariance-weighted variant. We have followed this approach in order to simulate the effect of any remaining residual correlation ρ present in real systems (we have assumed a ρ = 0.05) [63].
In the SMF reconstruction, magnitude images were generated as: with simulated sensitivity maps depicted in Fig A in S2 File satisfying the relationship X n k¼1 C 2 k ¼ 1; which holds in practice when the relative sensitivity maps are calculated as C k = |S k |/S SOS [72]. These sensitivity maps have been previously used in [73].
It should be noted that different scanner vendors can implement different SMF and SoS variants. In this work we have used the variants given in [55] for datasets acquired without undersampling in the k-space, i.e., R = 1, where R is the acceleration factor of the acquisition defined as the ratio of the total k-space phase-encoding lines over the number of k-space lines actually acquired. Notice that in the absence of noise S SOS = S SMF . Besides, for the particular case of a single coil with uniform sensitivity, i.e., n = 1 and C = 1, Eq (15) and Eq (16) become identical.
The 132 3D phantoms resulting from the procedures described in the three previous sections were contaminated with noise using a range of clinical signal-to-noise ratios (SNR) of 10, 15, 20 and 30, where SNR = S 0 /σ. In order to generate signals under equivalent conditions, for each value of σ the same noise realizations fe R k g and fe I k g were used to generate the final images S SOS and S SMF . All datasets were created simulating a scanner with 8 coils (see Figs A and B in S2 File).

Evaluation metrics
The performance of the algorithms was quantified by comparing the obtained reconstructions against the ground-truth via three main criteria: (i) the angular error in the orientation of fiber populations, (ii) the proper estimation of the number of fiber populations present in every voxel and (iii) the volume fraction error.
For the analyses, local peaks from the reconstructed fiber ODFs were identified as those vertices in the grid with higher values than their adjacent neighbors, considering only cases where magnitudes exceeded at least one tenth of the amplitude of the highest peak (i.e., 0.1Áf max ) [50]. From all identified peaks, the highest four where finally retained.
Next, we adopted some of the evaluation metrics widely used in the literature [9]. Specifically, we used the angular error, defined as the average minimum angle between the extracted peaks and the true fiber directions [74]: where M true is the true number of fiber populations, e m is the unitary vector along with the mth detected fiber peak and v k is the unitary vector along the kth true fiber direction. The volume fraction error of the estimated fiber compartments was assessed by means of the average absolute error between the estimated and the actual peak amplitudes: where f m is the normalized height of the mth detected fiber peak and f k is the volume fraction of the kth true fiber. As usual, the angular and volume fraction errors between each pair of fibers were measured by comparing the true fiber with the closest estimated fiber. Finally, the success rate (SR) was employed to quantify the estimation of the number of fiber compartments. The SR is defined as the proportion of voxels in which the algorithms estimate the right number of fiber compartments. To discriminate the different factors leading to an erroneous estimation, the mean number of over-estimated n + and under-estimated n − fiber populations were computed over the whole image [9].

Settings for the evaluation algorithms
Both RUMBA-SD and dRL-SD methods were implemented using in-house developed Matlab code, applying the same dictionaries H created from the signal generative model given in Eqs (1)-(3). These used M = 724 fiber orientations distributed on the unit sphere, with a mean angular separation between adjacent neighbor vertices of 8.36 degrees, and a standard deviation of 1.18 degrees.
To assess the effect of using dictionaries with optimal and non-optimal diffusivities, two different dictionaries were created and applied to the datasets described in subsections 1 (i.e., fiber bundles with different inter-fiber angles) and 2 (i.e., fiber bundles with different volume fractions). The first dictionary was generated by using the same diffusivities employed in the synthetic data, whereas the second dictionary was created from diffusivities estimated in regions of parallel fibers (outside the fiber crossing area) by means of a standard diffusion tensor fitting on the noisy data (i.e., dtifit tool in FSL package).
Similarly, two dictionaries were created to test the reconstructions on the data described in subsection 3 (i.e., "HARDI Reconstruction Challenge 2013" phantom data). In this case the model diffusivities and the 'true' diffusivities were deliberately set to different values in order to consider the possibility of model misspecification. The first dictionary was created with tensor diffusivities equal to λ 1 = 1.4 10 −3 mm 2 /s and λ 2 = λ 3 = 0.4 10 −3 mm 2 /s. Two isotropic compartments with diffusivities equal to 0.2 10 −3 mm 2 /s and 1.4 10 −3 mm 2 /s were also included. In the second dictionary the diffusivities were assumed to be equal to λ 1 = 1.6 10 −3 mm 2 /s and λ 2 = λ 3 = 0.3 10 −3 mm 2 /s, and the isotropic diffusivities were equal to 0.2 10 −3 mm 2 /s and 1.6 10 −3 mm 2 /s respectively.
The starting condition f°in all cases was set as a non-negative iso-probable spherical function [37]. The accuracy and convergence of both methods as a function of the number of iterations was investigated by repeating the calculations using 200 and 400 iterations, which is within the optimal range as suggested in [37] and [50]. The extended algorithms with TV regularization were also tested using 600 and 1000 iterations. The geometric damping and threshold parameters for dRL-SD were set to v = 8 and η = 0.06 respectively [37], see Appendix D in S1 File. For SoS-based data, n was fixed to the real number of coils in RUMBA-SD.
To differentiate the standard RUMBA-SD and dRL-SD algorithms from their regularized versions, we have added the term '+TV' to their names, i.e., RUMBA-SD+TV and dRL-SD +TV.

Real brain data
Diffusion MRI data were acquired from a healthy subject on a 3T Siemens scanner (Erlangen) located at the University of Oxford (UK). The subject provided informed written consent before participating in the study, which was approved by the Institutional Review Board of the University of Oxford. Whole brain diffusion images were acquired with a 32-channel head coil along 256 different gradient directions on the sphere in q-space with constant b = 2500 s/mm 2 . Additionally, 36 b = 0 volumes were acquired with in-plane resolution = 2.0 x 2.0 mm 2 and slice thickness = 2 mm. The acquisition was carried out without undersampling in the k-space (i.e., R = 1). Raw multichannel signals were combined using either the standard GRAPPA approach or the GRAPPA approach with the adaptive combination of the SMF available in the scanner, giving SoS and SMF-based datasets respectively. Then, the two resulting datasets were separately corrected for eddy current distortions and head motion as implemented in FSL [75].
A subset of 64 directions with nearly uniform coverage on the sphere was selected from the full set of 256 gradients directions, and measurements for this subset were used to 'create' an under-sampled version of the data, which also included 3 b = 0 volumes. The resulting HARDI sequence based on 64 directions is similar to those widely employed in clinical studies, thus results from this dataset are useful to evaluate the impact of the new technique on standard clinical data.

Gaussian versus non-gaussian noise models
The angular and volume fraction errors from the dRL-SD and RUMBA-SD reconstructions in the 90 synthetic phantoms with different inter-fiber angles are depicted in Fig 1, as well as in Fig C in S2 File. Fig 1 shows results using a dictionary created with the same diffusivities applied to generate the data (i.e., λ 1 = 1.7 10 −3 mm 2 /s and λ 2 = λ 3 = 0. 3  A set of patterns can be drawn from these results. First, RUMBA-SD was able to resolve fiber crossings with smaller inter-fiber angles (around 5 degrees and 10 degrees for datasets corrupted with Rician and nc-χ noise respectively). Second, RUMBA-SD produced volume fraction estimates with a higher precision (lower variance), even in phantoms where the fiber configuration was well-resolved by both methods. Third, although dRL-SD produced a relatively lower proportion of spurious fibers (n + ), RUMBA-SD produced a lower proportion of undetected fibers (n − ) leading to a higher success rate (SR). Finally, the performance of both methods was inferior when the dictionary was created using diffusivities estimated from a 'standard' DTI fitting in WM regions of parallel fibers. In line with previous findings on dRL-SD [50], optimal results were obtained from the sharper fiber response model.
All these points also hold for results obtained with other SNRs and with a differing number of algorithm iterations. Specifically, when a higher number of iterations was employed (i.e., 400), a lower proportion of n − and a higher proportion of n + was obtained with both methods.    The performance of RUMBA-SD was better in terms of the estimated volume fractions and the success rate, especially in the Rician noise case. In order to verify if the lower bias in the volume fractions estimated by RUMBA-SD could be explained only by its higher success rate, the calculation was repeated by considering only those voxels in each phantom where the two fiber populations were identified. After correcting for this factor, a noticeable advantage was still observed for RUMBA-SD (see left lower panel of Fig 2). When comparing these results with those from Fig 1 (and Fig D in S2 File), it is clear that TV regularization provides multiple benefits in both algorithms, including a superior ability to detect fiber crossings with smaller inter-fiber angles and a higher success rate. The later is due to a lower proportion of undetected fibers (n − ) and of spurious fibers (n + ). This pattern is evident in Fig 4, which depicts the peaks extracted from the fiber ODFs estimated from the SMF-based phantom with inter-fiber angle equal to 45 degrees and SNR = 15. Peaks are plotted as thin cylinders.

Original versus TV-regularized algorithms
As before, the above patterns were also observed in the analyses obtained with datasets based on other SNRs and with different number of iterations. In all cases RUMBA-SD+TV detected fiber crossings at lower inter-fiber angles. Fig 5 shows an example of this in the phantom with inter-fiber angle equal to 33 degrees corrupted with Rician noise and SNR = 15. 3 "HARDI reconstruction challenge 2013" phantom data  On the one hand, RUMBA-SD was able to resolve some fiber configurations that were not detected by dRL-SD, especially in voxels involving small inter-fiber angles or fiber crossings with a non-dominant tract. On the other hand, the spatially-regularized algorithms have substantially improved the performance of the original methods. Moreover, RUMBA-SD+TV was the method providing better reconstructions. These findings are in line with results from previous sections and remained valid also when different dictionaries, number of iterations and noise levels were employed.
Additional complementary results about the performance of RUMBA-SD in relation to several other reconstruction methods can be found in the website of the 'HARDI Reconstruction Challenge 2013': http://hardi.epfl.ch/static/events/2013_ISBI/workshop.html#results. An earlier version of RUMBA-SD took part in that Challenge, ranking number one in the 'HARDIlike' category (team name: 'Capablanca'). In the discussion section we provide additional information.

Real brain data
Fiber ODFs were estimated separately for each different SMF-and SoS-based dataset, including the original measured data with the full set of 256 gradient directions and its reduced form including a subset of 64 directions. In all cases, both RUMBA-SD and dRL-SD methods were implemented using the same dictionary, which was created assuming a sharp fiber response model with diffusivities equal to λ 1 = 1.7 10 −3 mm 2 /s and λ 2 = λ 3 = 0.3 10 −3 mm 2 /s, and two isotropic terms with diffusivities equal to 0.7 10 −3 mm 2 /s and 2.5 10 −3 mm 2 /s. Fig 7 shows the fiber ODFs estimated from the reduced SMF-and SoS-based datasets (i.e., data containing the reduced set of 64 gradient directions) in a coronal ROI on the right brain hemisphere. These results correspond to the reconstructions employing 200 iterations. Visual inspection of Fig 7 reveals that RUMBA-SD has produced sharper fiber ODFs than dRL-SD in both data. It has detected more clearly the fiber crossings. Interestingly, the fiber ODF profiles estimated by dRL-SD from the SoS-based data are smoother than those estimated from the SMF-based data. This behavior is less perceptible in the case of RUMBA-SD, suggesting that it could be more robust to different multichannel combination methods. Fig 8 depicts the fiber ODFs estimated with RUMBA-SD and RUMBA-SD+TV in a ROI of both the full and reduced SMF-based datasets. This region contains complex fiber geometries, including the mixture of the anterior limb of internal capsule (alic), the external capsule (ec) and part of the superior longitudinal fasciculus (slf) on the left brain hemisphere. Although in both cases multiple fibers were detected in the area of intersection, RUMBA-SD+TV has provided multi-directional fiber ODFs with a higher number of lobes, which may represent fiber crossings as well as intra-voxel fiber dispersion. The similarity of the reconstructions from the full and reduced datasets suggests that the method is robust with respect to the number of measurements, with the regularized version being the most robust.
In a subsequent analysis we examined the statistical properties of inter-fiber angles as estimated by all methods. Only white matter voxels where both methods detected one or two fibers are included in each plot, with inter-fiber angles in single fiber voxels assumed to be zero. Points on the main diagonal line characterize voxels where both methods gave identical inter-fiber angle estimates, whereas points above and below the main diagonal correspond to voxels where the two methods detected two fibers with different inter-fiber angle. The high density of points forming two secondary lines near the main diagonal indicates nearly similar reconstructions by both methods, with the angular differences being similar to the angular resolution of the reconstruction grid (i.e., about 8 degrees). The higher number of points above the main diagonal in both panels, especially for inter-fiber angles lower than 50 degrees, suggests that RUMBA-SD and RUMBA-SD+TV provide higher inter-fiber angles than dRL-SD and dRL-SD+TV, respectively. Finally, points located on the X and Y axes are voxels where one method detected two fibers while the other detected one. The very low density of points on the X axis of panel A for inter-fiber angles lower than 50 degrees (see the blue bracket) indicates that in nearly all voxels where dRL-SD detected two fibers, RUMBA-SD was also able to detect two fibers. In contrast, the high density of points on the Y axis in the same range of inter-fiber angles indicates that in many voxels RUMBA-SD detected two fibers whereas dRL-SD detected a single one. A similar but attenuated effect can be noticed in panel B, suggesting that TV regularization contributes to reduce the differences between both methods.

Discussion and Conclusions
In this study we propose a new model-based spherical deconvolution method, RUMBA-SD. In contrast to previous methods, usually based on Gaussian noise with zero mean, RUMBA-SD considers Rician and noncentral Chi noise models, which are more adequate for characterizing the non-linear bias introduced in the diffusion images measured in current 1.5T and 3T multichannel MRI scanners. Although recent progress has been made in new SD methods adapted to corrupted Rician data, e.g., see [41,76] to the best of our knowledge, our study provides the first SD extension to noncentral Chi noise. Furthermore, RUMBA-SD offers a very general estimation framework applicable to different datasets, with its flexibility emanating from two features: i) the explicit dependence between the likelihood model and the number of coils in the scanner and ii) the specific methodology employed to combine multichannel signals. In addition, the voxel-wise estimation of the noise variance adequately deals with potential deviations in the noise distribution, due, for instance, to accelerated MRI techniques or to preprocessing effects. We hope that the proposed technique will help extend SD methods to a wide range of datasets taken from different scanners and using different protocols.
This study adds to previous diffusion MRI studies trying to overcome the signal-dependent bias introduced by Rician and noncentral Chi noise. Apart from the robust DTI estimation methods in [77] and the earlier DTI study conducted by [78], the noise filtering techniques recently described in [79,80] are specially relevant. These techniques can be applied in the preprocessing steps prior to HARDI data estimation.
The performance of RUMBA-SD has been evaluated exhaustively against the state-of-theart dRL-SD technique. For that, we have used 132 different 3D synthetic phantoms, including 90 phantoms simulating fiber crossings with different inter-fiber angle, 41 phantoms simulating fiber crossings with different volume fraction, and the complex phantom designed for the "HARDI Reconstruction Challenge 2013" Workshop organized within the IEEE International Symposium on Biomedical Imaging. The comparison of these two methods has allowed us to weight the impact on the results of the Rician and noncentral Chi likelihood models included in RUMBA-SD, in relation to the Gaussian model assumed in dRL-SD. Since both approaches were implemented using the same dictionary of basis signals and similar reconstruction methods based on Richardson-Lucy algorithms adapted to Gaussian, Rician and noncentral Chi noise models, results should be considered comparable. Only voxels in white matter where both methods detected one or two fibers are included. The inter-fiber angle in voxels with a single fiber was assumed to be equal to zero. Points on the main diagonal line are those voxels where the inter-fiber angle estimated from both methods was identical, whereas points above and below the main diagonal correspond to voxels where the two methods detected two fibers but with different inter-fiber angle. Points located on the X and Y axes are voxels where one method detected two fibers whereas the other detected a single fiber. Taken together, findings from all synthetic datasets demonstrate the benefits of an adequate modelling of the noise distribution in the context of spherical deconvolution, and of the inclusion of TV regularization. Interestingly, RUMBA-SD resolved fiber crossings with smaller inter-fiber angles and smaller non-dominant fibers. Likewise, RUMBA-SD produced volume fraction estimates with higher accuracies and precision and produced, as well, a lower proportion of undetected fibers, resulting in a higher success rate (see Figs 1 and 2 in S2 File). On the other hand, the TV spatially-regularized versions of both dRL-SD and RUMBA-SD have substantially improved the performance of the original methods in all studied metrics (see Figs 3-6 in S2 File).
As previously mentioned, an earlier version of RUMBA-SD took part in the HARDI Reconstruction Challenge 2013, ranking number one in the 'HARDI-like' category. Notably, this position was shared with a reconstruction based on the CSD method included in Dipy software [81] (http://nipy.org/dipy/), which had applied a Rician denoising algorithm [82] to the raw diffusion MRI data prior to the actual CSD reconstruction. The superior performance of these two approaches strengthens the importance of taking into account the non-Gaussian nature of the noise. Moreover, it opens new questions on the optimal strategy to be followed. Should we divide the deconvolution process into to two disjoint steps (first denoising and then estimation) or it is more adequate to follow the unified approach proposed here?. The main advantage of the former is that it may benefit from including state-of-art denoising algorithms like the adaptative non-local means denoising method proposed in [82]. Conversely, the main advantage of the unified approach is that it provides a precise model to distinguish real signals from noise throughout all the 4D diffusion MRI data. Many of the advanced denoising algorithms that are currently applied in isolation were developed to filter volumetric (3D) data. Since each 3D image is processed individually, their mutual dependence in terms of orientation is ignored. In contrast, the unified approach described here provides a more general estimation framework that may be extended to include advanced similarity measures like those employed in [82], merging the benefits of both strategies. A new manuscript on the 'Challenge' phantom (currently under preparation) will provide additional information on the performance of RUM-BA-SD in relation to several reconstruction methods and in terms of connectivity metrics derived from fiber tracking analyses.
When applied to human brain data, RUMBA-SD has also achieved the best results, with its reconstructions showing the highest ability to detect fiber crossings (see . And although any conclusion derived from real data is hampered by the unknown anatomy at the voxel level, all previous results on synthetic data seem to support the validity of RUMBA-SD for real data.
Our findings can also be contrasted with those reported in [55]. In that study, the authors show that the SoS approach produces a signal-dependent bias that reduces the signal dynamic range and may subsequently lead to decreased precision and accuracy in fiber orientation estimates. Our study, however, suggests that the noncentral Chi noise in SoS-based data is not a major concern for the SD methods considered. Thus, heavier squashing of fiber ODFs when SoS reconstruction is used [55] is not that prominent with the SD if compared to diffusion ODF estimation methods [11]. This result may have different explanations for each technique. The robustness of dRL-SD may be explained, in part, by its lower over-all sensitivity to selection of the response function [50], which make it robust to the use of dictionaries estimated from either biased or unbiased signals. This behaviour may be additionally boosted by the inclusion of the damping factor in the RL algorithm. In contrast, the robustness of RUMBA-SD can be explained by the use of proper likelihood models that explicitly consider the bias as a function of the noise corrupting the data.
To finish, some limitations and future extensions of the study should be acknowledged. First, we have not evaluated the proposed method in synthetic data simulating partial Fourier k-space acquisitions and with parallel imaging using various acceleration factors (i.e., R>1), yet it may be interesting to consider it. Second, the RUMBA-SD estimation framework is based on a discrete approximation of the fiber ODF, which may be potentially extended to continuous functions on the sphere, like the spherical harmonics and the wavelets. Third, the TV regularization implemented in this study is based on a channel-by-channel first order scheme. New studies may be designed to compare different regularization techniques such as higher order TV, vectorial TV and the fiber continuity approach introduced in [56], to mention only a few examples. Fourth, different strategies for creating the signal dictionary could be explored, like using mixtures of intra-compartment models to capture different diffusion profiles, or applying more appropriate models to fit multi-shell data [31,83]. Fifth, the recursive calibration of the single-fiber response function proposed by [84] may be another possible add-on. Finally, it is worth mentioning that the inversion algorithm behind RUMBA-SD is not limited to fiber ODF reconstructions, but it can be also applied to solve other linear mixture models from diffusion MRI data. It was recently showed that some microstructure imaging methods such as ActiveAx and NODDI can be reformulated as convenient linear systems, however, the deconvolution methods proposed for them assume Gaussian noise and are performed on a voxel-by-voxel basis [85]. Here the iterative scheme proposed in RUMBA-SD could be used to address both limitations, potentially leading to improved reconstructions in microstructure imaging.