Jian Cheng
email: jian_cheng@med.unc.edu

Rachid Deriche
email: rachid.deriche@inria.fr

Tianzi Jiang
email: jiangtz@nlpr.ia.ac.cn

Dinggang Shen
email: dgshen@med.unc.edu

Pew-Thian Yap
email: ptyap@med.unc.edu

Non-Negative Spherical Deconvolution (NNSD) for Fiber Orientation Distribution Function Estimation

In diffusion Magnetic Resonance Imaging (dMRI), Spherical Deconvolution (SD) is a commonly used approach for estimating the fiber Orientation Distribution Function (fODF). As a Probability Density Function (PDF) that characterizes the distribution of fiber orientations, the fODF is expected to be non-negative and to integrate to unity on the continuous unit sphere S 2 . However, many existing approaches, despite using continuous representation such as Spherical Harmonics (SH), impose non-negativity only on discretized points of S 2 . Therefore, nonnegativity is not guaranteed on the whole S 2 . Existing approaches are also known to exhibit false positive fODF peaks, especially in regions with low anisotropy, causing an over-estimation of the number of fascicles that traverse each voxel. This paper proposes a novel approach, called Non-Negative SD (NNSD), to overcome the above limitations. NNSD offers the following advantages. First, NNSD is the first SH based method that guarantees non-negativity of the fODF throughout the unit sphere. Second, unlike approaches such as Maximum Entropy SD (MESD), Cartesian Tensor Fiber Orientation Distribution (CT-FOD), and discrete representation based SD (DR-SD) techniques, the SH representation allows closed form of spherical integration, efficient computation in a low dimensional space resided by the SH coefficients, and accurate peak detection on the continuous domain defined by the unit sphere. Third, NNSD is significantly less susceptible to producing

false positive peaks in regions with low anisotropy. Evaluations of NNSD in comparison with Constrained SD (CSD), MESD, and DR-SD (implemented using L1regularized least-squares with non-negative constraint), indicate that NNSD yields improved performance for both synthetic and real data. The performance gain is especially prominent for high resolution (1:25 mm) 3 data.

Introduction

Diffusion MRI (dMRI) non-invasively reveals the microstructure of white matter by capturing the diffusion patterns of water molecules. The most widely used dMRI approach, Diffusion Tensor Imaging (DTI), is however incapable of describing complex diffusion processes due to its reliance on the assumption of Gaussian diffusion [START_REF] Johansen-Berg | Diffusion MRI: From quantitative measurement to In vivo neuroanatomy[END_REF]. Hence, many methods have been developed in recent years to harness the power of High Angular Resolution Diffusion Imaging (HARDI) to characterize the non-Gaussian diffusion by evaluating quantities such as the Ensemble Average Propagator (EAP) [START_REF] Wedeen | Mapping Complex Tissue Architecture With Diffusion Spectrum Magnetic Resonance Imaging[END_REF][START_REF] Cheng | Nonnegative Definite EAP and ODF Estimation via a Unified Multi-Shell HARDI Reconstruction[END_REF], the diffusion Orientation Distribution Function (dODF) [START_REF] Tuch | Q-ball imaging[END_REF][START_REF] Descoteaux | Regularized, Fast and Robust Analytical Q-ball Imaging[END_REF][START_REF] Cheng | Nonnegative Definite EAP and ODF Estimation via a Unified Multi-Shell HARDI Reconstruction[END_REF], and the fiber Orientation Distribution Function (fODF) [START_REF] Tournier | direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution[END_REF][START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF][START_REF] Alexander | Maximum entropy spherical deconvolution for diffusion MRI[END_REF][START_REF] Jian | A Unified Computational Framework for Deconvolution to Reconstruct Multiple Fibers from Diffusion Weighted MRI[END_REF].

A popular method, called Spherical Deconvolution (SD), has been shown to be effective for estimating the fODF by assuming that the measured diffusion-weighted signal can be obtained via spherically convolving a latent fODF with a fiber response function, which can be estimated from voxels known to be traversed by a single fascicle [START_REF] Tournier | direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution[END_REF][START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF][START_REF] Jian | A Unified Computational Framework for Deconvolution to Reconstruct Multiple Fibers from Diffusion Weighted MRI[END_REF]. The fODF can hence be recovered via an inverse problem by deconvolving the signal with the fiber response function. SD methods can be classified into two categories, 1) continuous representation based SD (CR-SD), which is normally based on representation using Spherical Harmonics (SH) [START_REF] Tournier | direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution[END_REF][START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF][START_REF] Alexander | Maximum entropy spherical deconvolution for diffusion MRI[END_REF], and 2) discrete representation based SD (DR-SD), which is based on a discrete mixture of rotated versions of the fiber response function [START_REF] Jian | A Unified Computational Framework for Deconvolution to Reconstruct Multiple Fibers from Diffusion Weighted MRI[END_REF][START_REF] Dell'acqua | A model-based deconvolution approach to solve fiber crossing in diffusion-weighted MR imaging[END_REF][START_REF] Landman | Resolution of crossing fibers with constrained compressed sensing using diffusion tensor mri[END_REF].

Existing SD methods in both categories have some common limitations. First, they often result in false-positive fODF peaks [START_REF] Tournier | direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution[END_REF][START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF][START_REF] Alexander | Maximum entropy spherical deconvolution for diffusion MRI[END_REF][START_REF] Johansen-Berg | Diffusion MRI: From quantitative measurement to In vivo neuroanatomy[END_REF][START_REF] Jian | A Unified Computational Framework for Deconvolution to Reconstruct Multiple Fibers from Diffusion Weighted MRI[END_REF][START_REF] Landman | Resolution of crossing fibers with constrained compressed sensing using diffusion tensor mri[END_REF][START_REF] Weldeselassie | Symmetric positive semi-definite cartesian tensor fiber orientation distributions (CT-FOD)[END_REF], especially in low anisotropic regions, such as the gray matter and cerebrospinal fluid (CSF) regions. Second, they normally fall short in ensuring that the estimated fODF is a proper probability density function, because non-negativity and unit integral over the unit sphere are not always enforced. Most SD methods, including the popular Constrained SD (CSD) [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF], consider non-negativity only on discretized points of the unit sphere, but not the whole S 2 [START_REF] Jian | A Unified Computational Framework for Deconvolution to Reconstruct Multiple Fibers from Diffusion Weighted MRI[END_REF][START_REF] Dell'acqua | A model-based deconvolution approach to solve fiber crossing in diffusion-weighted MR imaging[END_REF][START_REF] Landman | Resolution of crossing fibers with constrained compressed sensing using diffusion tensor mri[END_REF]. Thus, negative values are still a problem at points where non-negativity is not explicitly enforced. To our knowledge, Maximum Entropy SD (MESD) [START_REF] Alexander | Maximum entropy spherical deconvolution for diffusion MRI[END_REF] and Cartesian Tensor Fiber Orientation Distribution (CT-FOD) [START_REF] Weldeselassie | Symmetric Positive-Definite Cartesian Tensor Orientation Distribution Functions (CT-ODF)[END_REF][START_REF] Weldeselassie | Symmetric positive semi-definite cartesian tensor fiber orientation distributions (CT-FOD)[END_REF] are the only existing methods that ensure non-negativity throughout S 2 . However, these methods are inefficient without the closed form of spherical integration needed to compute the fODF and are susceptible to error associated with the numerical spherical integration. Ad-hoc normalization is also required in these methods to arrive at the fODF with unit integral. Third, for estimation of the fODF with reasonable accuracy, DR-SD methods require a significant number of rotated versions of the fiber response function in directions that are distributed densely on the unit sphere, significantly increasing the dimensionality and time cost of the optimization problem.

This paper proposes a novel CR-SD method, called Non-Negative Spherical Deconvolution (NNSD), for fODF estimation. Non-negativity is achieved by representing the square root of the fODF as a linear combination of SH basis functions. To our knowledge, NNSD is the first SH based SD method to guarantee non-negativity on the whole S 2 , not only on discretized points on the unit sphere, as in CSD [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF]. Compared with MESD and CT-FOD, NNSD ensures unit integral and is more efficient due to the existence of a closed-form expression needed to arrive at the fODF. Moreover, NNSD significantly reduces the spurious peaks in regions with low anisotropy.

The rest of the paper is organized as follows. Section 2 provides an overview of three SD methods, i.e. CSD [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF], MESD [START_REF] Alexander | Maximum entropy spherical deconvolution for diffusion MRI[END_REF], and DR-SD (L1-regularized nonnegative least squares with non-negative constraint) [START_REF] Jian | A Unified Computational Framework for Deconvolution to Reconstruct Multiple Fibers from Diffusion Weighted MRI[END_REF][START_REF] Landman | Resolution of crossing fibers with constrained compressed sensing using diffusion tensor mri[END_REF]. Section 3 describes NNSD and the associated Riemannian gradient descent algorithm. In Section 4, NNSD is empirically evaluated in comparison with the three methods discussed in Section 2.

Background on SD Methods

If the fODF is represented as Φ(u) = ∑ L l=0 ∑ l m= l c lm Y m l (u), and the axisymmetric fiber response function along the z-axis as

H(uj(0; 0; 1)) = L ∑ l=0 h l Y 0 l (u); (1)
where Y m l (u), u P S 2 , is the l order and m degree real Spherical Harmonic basis function [START_REF] Johansen-Berg | Diffusion MRI: From quantitative measurement to In vivo neuroanatomy[END_REF][START_REF] Descoteaux | Regularized, Fast and Robust Analytical Q-ball Imaging[END_REF], then the diffusion signal can be represented as the convolution of the response function with the fODF [START_REF] Tournier | direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution[END_REF]:

E(u) = Z S 2 H(ujr)Φ(r)dr = L ∑ l=0 l ∑ m= l r 4π 2l + 1 c lm h l Y m l (u): (2)
Note that convolution here in conveniently performed via multiplication of the SH coefficients c lm and h l . The SH coefficient vector c = (c 00 ; : : : ; c LL ) T of the fODF can be estimated by minimizing kMc Ek 2 ;

(

where

M = [ p 4π=2l + 1h l Y m l (u i )]
is an N by (L + 1)(L + 2)=2 matrix, and E = (E 1 ; : : : ; E N ) T is a vector consisting of the signal measurements. The least squares solution is (M T M)1 M T E. However, the least squares formulation does not take into account the non-negativity of the fODF. Constrained SD (CSD) [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF] attempts to address this issue by iteratively suppressing the negative values by a discrete reconstruction on discrete points of the fODF, i.e., they solve

c (k+1) = arg min c kMc Ek 2 + λ 2 CSD kL (k) ck 2 ; (4)
where

L (k) i; j = P i; j if Pc (k) i < τ, L (k)
i; j = 0 if

Pc (k)
i > τ, and P is the SH matrix for reconstructing the fODF at discretized points on S 2 . τ is a threshold normally chosen as 0:1 [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF]. Note that CSD imposes non-negativity only on a set of discretized points, not on the whole unit sphere. Although CSD significantly reduces the negative values compared with the original version of the SD technique described in [START_REF] Tournier | direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution[END_REF], there is still a significant amount of negative values, even on points where non-negativity is imposed.

In MESD [START_REF] Alexander | Maximum entropy spherical deconvolution for diffusion MRI[END_REF], a CR-SD method, a representation of the fODF is derived based on the maximum entropy principle:

Φ(rjfλ i g K i=0 ) = exp(λ 0 + K ∑ i=1 λ i H(v i jr)); (5)
which is an exponential function of mixture of fiber response functions determined by a set of parameters fλ i g K i=0 . These parameters can be determined by applying the Levenberg-Marquardt method to solve the following nonlinear least-squares problem:

min fλ i g K i=0 N ∑ j=1 Z S 2 H(u j jr)Φ(rjfλ i g K i=0 )dr E j 2 : (6)
The exponential representation naturally ensures non-negativity on S 2 . However, unlike (2), the spherical integration in (6) needs to be approximated numerically. An ad-hoc normalization of fODF is also needed to achieve unit integral after determining fλ i g K i=0 . Based on the method of Lagrange multipliers [START_REF] Alexander | Maximum entropy spherical deconvolution for diffusion MRI[END_REF], the directions fv i g K i=1 should ideally be the N sampling points corresponding to the signal measurements. However, to compromise between accuracy and speed, in practice K < N evenly distributed points on S 2 are used 1 .

In DR-SD, the fODF was represented using discretized points fw i = Φ(u i )g [START_REF] Jian | A Unified Computational Framework for Deconvolution to Reconstruct Multiple Fibers from Diffusion Weighted MRI[END_REF][START_REF] Dell'acqua | A model-based deconvolution approach to solve fiber crossing in diffusion-weighted MR imaging[END_REF][START_REF] Landman | Resolution of crossing fibers with constrained compressed sensing using diffusion tensor mri[END_REF]. The SD problem in this case is normally formulated as min

w kAw Ek 2 + λ L1 jjwjj 1 ; s.t. w ! 0; (7)
where A is the matrix with each column containing a rotated version of the fiber response function. For simplicity, we call this method the L1-regularized Non-Negative Least Squares (L1NNLS). When λ L1 = 0, Eq. ( 7) becomes the Non-Negative Least Squares (NNLS) [START_REF] Jian | A Unified Computational Framework for Deconvolution to Reconstruct Multiple Fibers from Diffusion Weighted MRI[END_REF][START_REF] Dell'acqua | A model-based deconvolution approach to solve fiber crossing in diffusion-weighted MR imaging[END_REF]. Since the fODF is represented using dis-cretized points, its maxima are restricted to these points, not the continuous S 2 , eventually limiting the angular resolution for peak detection.

Non-Negative Spherical Deconvolution (NNSD)

The square root representation has been proposed for dODFs [START_REF] Cheng | A Riemannian Framework for Orientation Distribution Function Computing[END_REF] and EAPs [START_REF] Cheng | Diffeomorphism Invariant Riemannian Framework for Ensemble Average Propagator Computing[END_REF] and has been shown to be effective for non-negative estimation of these quantities [START_REF] Cheng | Nonnegative Definite EAP and ODF Estimation via a Unified Multi-Shell HARDI Reconstruction[END_REF]. In this work, we propose to employ this representation for non-negative estimation of the fODF. This is achieved by representing the square root of the fODF as a linear combination of SH basis functions: [START_REF] Descoteaux | Impact of Rician Adapted Non-local Means Filtering on HARDI[END_REF] where

Φ(ujc) = 2 L ∑ l=0 l ∑ m= l c lm Y m l (u) 3 2 = 2L ∑ α=0 α ∑ β = α 2 L ∑ l;m L ∑ l H ;m H c lm c l H m H Q mm H β ll H α 3 Y β α (u);
Q mm H β ll H α = R S 2 Y m l (u)Y m H l H (u)Y β α (u)
du is a constant tensor computed by the integration of three sets of SH basis functions [START_REF] Cheng | Nonnegative Definite EAP and ODF Estimation via a Unified Multi-Shell HARDI Reconstruction[END_REF]. Q mm H β ll H α can be calculated from the Wigner 3-j symbol. Note that the sum over α is up to 2L, because if α > 2L, then Q mm H β ll H α = 0 based on the property of Wigner 3-j symbol. In contrast to CSD, we use the SH basis function to represent the square root of the fODF, instead of the fODF itself in CSD. Note that since R S 2 Φ(u)du = 1 and due to the orthogonality of the SH basis functions, we have 2), the convolved diffusion signal is

R S 2 Φ(u)du = kck 2 = 1. Based on Eq. (
E(u) = 2L ∑ α=0 α ∑ β = α L ∑ l;m L ∑ l H ;m H r 4π 2α + 1 c lm c l H m H Q mm H β ll H α h α Y β α (u) = c T K(u)c; (9)
where for any fixed u, K(u) is a (L + 1) ¢(L+2)=2 square matrix with the elements

K mm H ll H (u) = 2L ∑ α=0 α ∑ β = α r 4π 2α + 1 Q mm H β ll H α h α Y β α (u): (10)
Based on Eq. ( 9), we propose to estimate c from the measured fE i g N i=1 by minimizing

J(c) = 1 2 N ∑ i=1 c T K(u)c E i ¡ 2 + 1 2 c T Λ c; s.t. kck = 1; ( 11
)
where Λ is a diagonal matrix with elements Λ lm = λ NNSD l 2 (l + 1) 2 for the Laplace-Beltrami regularization [START_REF] Descoteaux | Regularized, Fast and Robust Analytical Q-ball Imaging[END_REF] . The Riemannian gradient ∇J(c) is the projection of Euclidean gradient ∂ J(c) ∂ c onto the tangent space of c, i.e.

∇J(c) = ∂ J(c) ∂ c c T ∂ J(c) ∂ c c; ∂ J(c) ∂ c = N ∑ i=1 2 c T K(u i )c E i ¡ K(u i )c + Λ c: (12) Then the Riemannian gradient descent of J(c) is c (k+1) = Exp c (k) dt ∇J(c) k∇J(c)k ;
Exp c (v) = c cos kvk+ v kvk sin kvk; [START_REF] Tournier | direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution[END_REF] where c (k) is the estimated c in the k-th step, dt is the step size chosen from the inexact line search in (0; dt 0 ] , and Exp c (v) is the exponential map on the sphere [START_REF] Cheng | A Riemannian Framework for Orientation Distribution Function Computing[END_REF].

The whole process is summarized in Algorithm 1. We use the isotropic fODF with c (0) = (1; 0; : : : ; 0) T for initialization, and experimentally choose dt 0 = 0:1. We found that NNSD is robust to noise and λ NNSD = 0 works in most cases. To obtain a local minimum of J(c), the standard stopping condition can be set as

J(c (k 1) ) J(c (k) ) J(c (k 1) ) < δ : (14)
A small value for δ should be chosen such that NNSD converges to a local minimum of Eq. [START_REF] Landman | Resolution of crossing fibers with constrained compressed sensing using diffusion tensor mri[END_REF]. Experimentally δ = 10 2 works well in most cases. If δ is much smaller than 10 2 , we found that although the cost function J(c) becomes a little lower and the fODFs in regions with high anisotropy become a little shaper, the fODFs in regions with low anisotropy however become fuzzier with more spurious peaks. Thus we propose an adaptive stopping condition such that if GFA(c (k) ) < T and J(c

(k 1) ) J(c (k) ) J(c (k 1) ) < δ 0 , or if GFA(c (k) ) ! T and J(c (k 1) ) J(c (k) )
J(c (k 1) ) < 0:01δ 0 , then the gradient descent stops, where

δ 0 = 10 2 , GFA(c) = r 1 c 2 00 kck 2 = q 1 c 2 00
is the Generalized Fractional Anisotropy (GFA) [START_REF] Tuch | Q-ball imaging[END_REF], and T P [0; 1] is a threshold depending on the noise level. For simplicity, we call NNSD with this adaptive stopping condition NNSD-ASC. When T = 0, NNSD-ASC becomes NNSD with δ = 0:01δ 0 = 10 4 , and when T = 1, NNSD-ASC becomes NNSD with δ = δ 0 = 10 2 . Both NNSD and NNSD-ASC converge fast, normally in a dozen of steps. Based on our C++ implementation, when L = 6, fODFs in 1000 voxels can be estimated within 8 seconds using an ordinary laptop, which is much faster than MESD in Camino [START_REF] Alexander | Maximum entropy spherical deconvolution for diffusion MRI[END_REF]. Note that in each iteration step, the fODF Φ(ujc (k) ) is naturally non-negative in the whole S 2 due to the square operator in Eq. ( 8), and it has unit integral since kc (k) k = 1.

Experiments

We compared the proposed NNSD and NNSD-ASC with CSD, MESD, and L1NNLS using synthetic and real data. Output: SH coefficient c and the fODF Φ(ujc). begin initialization: c (0) = (1; 0; :::; 0) T , k = 0 ; // the isotropic fODF repeat calculate v = ∇J(c (k) ) in Eq. ( 12) ; if kvk < ε then break; choose step size dt P (0; dt 0 ] via line search;

c (k+1) = Exp c (k) dt v kvk , k 2 k + 1;
until Under a certain stopping condition; Φ(ujc) = Eq. ( 8) end Parameters: Note that in all the following experiments, MRtrix2 with default parameters was used for CSD [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF], and 321 evenly distributed orientations were used to generate P for CSD as suggested [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF][START_REF] Jian | A Unified Computational Framework for Deconvolution to Reconstruct Multiple Fibers from Diffusion Weighted MRI[END_REF]. Accelerated MESD with K = 16 was used as suggested in Camino. We experimentally set λ L1 = 1 for L1NNLS. Note that ideally the optimal λ L1 in Eq. ( 7) should be different in each voxel depending on E. We normalize the columns of A and E to have unit norm so that we can use the same λ L1 in all cases. For fair comparison, we intentionally set λ NNSD = 0 to switch off Laplace-Beltrami regularization, and set T = 0:5 in NNSD-ASC. For CSD, MESD, and NNSD, NNSD-ASC, the local maxima of the estimated fODF larger than the mean of the minimal and maximal values are detected by gradient ascent on the continuous unit sphere S 2 [START_REF] Tournier | direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution[END_REF][START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF]. For L1NNLS, the peaks are detected from the 321 orientations.

Generation of Synthetic Data: Synthetic signals were generated using a mixture of tensor model E(u) = ∑ K k=1 p k exp bu T D k u ¡ with b = 1500 s=mm 2 and 60 evenly distributed directions fu i g. For simulating crossings, two tensors with equal weights p 1 = p 2 = 0:5 and common eigenvalues (λ 1 ; λ 2 ; λ 2 ) but different orientations were used. Rician noise with signal-to-noise ratio (SNR) 1=σ was added to the signal, where σ is the standard deviation of the complex Gaussian noise. Anisotropy and Non-Negativity: 1000 realizations of Rician noise corrupted signal were generated with SNR = 15 and 30, respectively from a ground truth isotropic signal profile with eigenvalues (0:7; 0:7; 0:7) ¢ 10 3 mm 2 =s and an anisotropic signal profile with eigenvalues (1:7; 0:2; 0:2) ¢10 3 mm 2 =s. The Gaussian fiber response function with eigenvalues (1:7; 0:2; 0:2) ¢10 3 mm 2 =s was used in all methods. From the estimated fODFs, we recorded the GFA values [START_REF] Tuch | Q-ball imaging[END_REF] and the proportion of the negative fODF values. We ignored negative values close to zero and only take into consideration negative values with absolute values larger than 1% of the maximal fODF value. This evaluation was performed on 5121 points generated via tessellation of the unit sphere. The top left sub-figure in Fig. 1 although CSD considers non-negativity, it still suffers from negative values. NNSD-ASC and MESD guarantee non-negativity throughout S 2 . Note that the fODFs estimated by L1NNLS are non-negative on the discrete 321 points where non-negativity is imposed; however, since fODF value for other points are unknown, we cannot compute the proportion of negative values for L1NNLS. The top right sub-figure shows that CSD, MESD and L1NNLS obtain very high GFA even for isotropic diffusion, indicating an over-estimation of the number of fascicles. Only NNSD-ASC yields large contrast between isotropic and anisotropic signal profiles. Note that due to the sparsity constraint used in L1NNLS, the GFA values of fODFs estimated by L1NNLS for both isotropic and anisotropic signals are always close to 1. Simulation of Crossing Fibers: Corrupted signals were generated from two tensor model with SNR=10, eigenvalues [1:7; 0:2; 0:2] ¢10 3 mm 2 =s, and with different crossing angles in [30 ; 90 ]. For all methods, we estimated the fODFs with the same Gaussian fiber response and detected the local peaks. We recorded the success ratio, which is the proportion of trials with the correct number (i.e. two) of detected peaks, then calculated the Mean Difference of Angles (MDA) only in successful trials [START_REF] Descoteaux | Regularized, Fast and Robust Analytical Q-ball Imaging[END_REF]. The second row in Fig. 1 shows the success ratio and the MDA values for different crossing angles. MESD and CSD with L = 10 yield higher success ratios for small crossing angles, but lower success ratios for larger crossing angles. However, we are more concerned with large crossing angles, since this is more likely to happen in real data. With L = 6, NNSD-ASC generally yields higher success ratios than CSD and the lowest MDA values among all methods. L1NNLS generally yields lower success ratios, probably because one has to tune λ L1 for different fiber configurations to achieve optimal results, while we had just simply set λ L1 = 1. L1NNLS yields larger MDA values than others, probably due to the limited angular resolution. The large MDA values given by MESD are probably due to a larger K is needed for more accurate peak detection; this however requires a much longer computation time.

Real Data: Evaluation was also performed using real human data with b = 2000 s=mm 2 , 120 gradient directions, TR/TE=12400=116 ms, 2 mm isotropic voxel dimensions. We set L = 6 for CSD and NNSD, NNSD-ASC. MRtrix was used for CSD and for the estimation of the fiber response function from the voxels with FA > 0:7 [START_REF] Tournier | direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution[END_REF][START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF]. The results are showed in Fig. 2. NNSD with δ = 10 4 and NNSD-ASC with T = 0:5 yield similar results, which means the outcome is insensitive to T for this data. We calculated the GFA maps [START_REF] Tuch | Q-ball imaging[END_REF] from the estimated fODFs, and set them as the background color in every sub-figure. Consistent with the results in Fig. 1, the results in Fig. 2 for CSD, MESD and L1NNLS show a significant amount of false positive peaks, which is especially evident in the regions with low anisotropy, as indicated by the low background GFA values. Because of the sparsity consideration in L1NNLS, normally only less than 5 out of the 321 samples of fODFs have non-zero values, thus the fODFs by L1NNLS have very high GFA values close to 1, based on the GFA definition. NNSD and NNSD-ASC dramatically reduced the spurious peaks of the estimated fODFs; see the yellow circles in Fig. 2. Note that in isotropic regions, the fODFs estimated by NNSD and NNDS-ASC are very close to being isotropic. Although all fODFs estimated by NNSD have unit integral, the values of isotropic fODFs are around 1 4π and the peaks of anisotropic fODFs are normally larger than 1. Hence, when visualizing fODFs at the same scale, the isotropic fODFs estimated by NNSD/NNSD-ASC look like dots which have much small sizes compared with anisotropic fODFs. We would like to point out that most existing papers on SD overlaid the estimated fODF glyphs on a FA or GFA map calculated from other models like DTI, not from the estimated fODF itself [START_REF] Tournier | direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution[END_REF][START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF][START_REF] Jian | A Unified Computational Framework for Deconvolution to Reconstruct Multiple Fibers from Diffusion Weighted MRI[END_REF][START_REF] Landman | Resolution of crossing fibers with constrained compressed sensing using diffusion tensor mri[END_REF][START_REF] Weldeselassie | Symmetric positive semi-definite cartesian tensor fiber orientation distributions (CT-FOD)[END_REF]. The proposed NNSD and NNSD-ASC obtain a very clean fODF field with similar peaks as detected by CSD in highly anisotropic regions and more isotropic fODFs in the isotropic regions.

High-Resolution Data

We also evaluated the SD methods using high resolution (1:25 mm) 3 data obtained from the Human Connectome Project (HCP) 3 . The data is challenging because the SNR is low due to the small voxel dimensions. It has three shells, 90 directions per shell, b = 1000=2000=3000; s=mm 2 . We performed NNSD (δ = 10 4 ), NNSD-ASC (δ = 10 2 ), CSD, MESD, L1NNLS on a single shell with b = 2000 s=mm 2 ; see Fig. 3. Although NNSD and NNSD-ASC yield results similar to the previous data (see Fig. 2), NNSD-ASC demonstrated its robustness to the spurious peaks in this high-resolution data. In the isotropic areas, only NNSD-ASC obtained isotropic fODFs, while other methods obtained many false positive peaks. See the yellow circles. In the anisotropic areas, NNSD and NNSD-ASC obtained the sharpest fODFs. fODFs by MESD are also very sharp in the anisotropic areas; however there are some small spurious lopes, as shown in the glyphs, similarly to the results shown in Fig. 2. fODFs estimated by CSD for this high-resolution data are not as sharp as those estimated by NNSD/NNSD-ASC and MESD. This can be observed in the areas marked by the yellow squares. L1NNLS as before obtained the most anisotropic GFA map, indicating spurious peaks in isotropic areas.

Discussion and Conclusion

In all the above experiments, for fair comparisons between methods, we did not preprocess the data via denoising to enhance the SNR. We also did not consider spatia regularization for local smoothness. In practical applications, good denoising methods such as non-local means [START_REF] Descoteaux | Impact of Rician Adapted Non-local Means Filtering on HARDI[END_REF] and good spatial regularization can improve the final results. Although in this paper, we have only considered diffusion signals sampled from a single shell (e.g. one b value), the proposed NNSD can be used in multiple shell data by considering a signal response function defined in the 3D space.

The SH coefficients fh l (b)g for the different b-values can then be used in Eq. (1).

We participated ISBI 2013 HARDI reconstruction challenge4 using the proposed NNSD with non-local mean spatial regularization [START_REF] Cheng | Non-Local Non-Negative Spherical Deconvolution for Single and Multiple Shell Diffusion MRI[END_REF]. This challenge considered 3 categories of sampling schemes (DTI, HARDI, and multiple shell DSI-like schemes) and 3 SNRs (SNR=10,20,30), amounting to a total of 9 datasets. Compared with other methods in this challenge, NNSD was ranked as the best technique in terms of local fiber orientation accuracy in all 9 datasets 5 .

In this paper, we proposed a novel SD method, called Non-Negative SD (NNSD). NNSD is based on the square root representation of the fODF and can be solved by a Riemannian gradient descent algorithm with an adaptive stopping condition. Comparison with existing SD methods, including CSD, MESD, and L1NNLS, demonstrated the advantages of NNSD: 1) NNSD ensures that the estimated fODF is nonnegative with unit integral throughout S 2 ; 2) NNSD significantly reduces the false positive peaks and yields high contrast between isotropic and anisotropic regions; and 3) Due to the SH representation, NNSD is efficient and allows accurate peak detection on the continuous unit sphere. The experiments showed that the proposed NNSD-ASC works well in both isotropic and anisotropic areas, especially for the high-resolution data. The clean fODF field obtained by NNSD will be helpful for important applications such as fiber tractography.

Algorithm 1 :

 1 Non-Negative Spherical Deconvolution (NNSD)Input: Diffusion signal attenuation measurements fE i g N i=1 .

= 6 Fig. 1

 61 Fig.1The first row shows the proportion of negative values and the GFA values of the fODFs estimated from the isotropic signal profiles (denoted as ISO) and anisotropic signal profiles (denoted as ANI). The error bars indicate the standard deviations. The second row shows the success ratios and the MDA values of the different methods for crossings generated with various angles.

1 Fig. 2

 12 Fig. 2 First row: the whole coronal view of the fODF field using NNSD-ASC with T = 0:5, and close-up views of results generated using NNSD-ASC and NNSD. Second row: close-up views for results generated using CSD, MESD, and L1NNLS. The colors of glyphs indicate directions. The backgrounds are the GFA maps.

1 Fig. 3

 13 Fig. 3 First row: coronal view of the fODF field estimated using NNSD-ASC with T = 0:5, and close-up views of results generated using NNSD-ASC and NNSD. Second row: close-up views of results generated using CSD, MESD and L1NNLS. The RGB color of glyphs indicates the directions. The backgrounds consist of the GFA maps.

Camino: http://cmic.cs.ucl.ac.uk/camino/

MRtrix: http://www.brain.org.au/software/mrtrix/

http://www.humanconnectome.org/data/

http://hardi.epfl.ch/static/events/2013 ISBI/

http://hardi.epfl.ch/static/events/2013 ISBI/ static/talk Max.pdf


Jian Cheng

Rachid Deriche

Tianzi Jiang

Dinggang Shen

Pew-Thian Yap

Non-Local Non-Negative Spherical Deconvolution for Single and Multiple Shell Diffusion MRI

I. Introduction

In diffusion MRI (dMRI), Spherical Deconvolution (SD) is a category of methods which estimate the fiber Orientation Distribution Function (fODF). Existing SD methods, including the widely used Constrained SD [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained superresolved spherical deconvolution[END_REF], normally have two common limitations: 1) the non-negativity constraint of the fODFs is not satisfied in the continuous sphere; 2) many spurious peaks are detected, especially in the regions with low anisotropy; In [START_REF] Cheng | Non-Negative Spherical Deconvolution for Fiber Orientation Distribution Estimation[END_REF], we proposed a novel SD method, called Non-Negative SD (NNSD), to avoid these two limitations. NNSD guarantees the non-negativity constraint of fODFs in the continuous sphere S 2 , and it is robust to the false positive peaks. In this abstract, we propose Non-Local NNSD (NLNNSD) which considers non-local spatial information and Rician noise in NNSD, and apply it to the testing data in ISBI contest.

II. Method

We represent the square root of fODF Φ(u) as a linear combination of real Spherical Harmonic (SH) basis Y m l (u) with even order, i.e.

Φ(u) = L l=0 l m=-l c lm Y m l (u) 2 = 2L α=0 α β=-α L l,m L l ,m c lm c l m Q mm β ll α Y β α (u), where u ∈ S 2 , Q mm β ll α = S 2 Y m l (u)Y m l (u)Y β α (
u)du is the integral constant of three SHs which can be calculated from the Wigner 3-j symbol. Then based on the closed form of spherical convolution using SH basis, for a given axisymmetric fiber response function along z-axis H(qu|(0, 0, 1)) = L l=0 h l (q)Y 0 l (u), the convolved diffusion signal is

E(qu) = 2L α=0 α β=-α L l,m L l ,m 4π 2α + 1 c lm c l m Q mm β ll α h α (q)Y β α (u) = c T K(qu)c (1)
where for any fixed vector q = qu, K(u) is a square matrix with the elements

K mm ll (qu) = 2L α=0 α β=-α 4π 2α+1 Q mm β ll α h α (q)Y β α (u).
Then NNSD [START_REF] Cheng | Non-Negative Spherical Deconvolution for Fiber Orientation Distribution Estimation[END_REF] is to estimate c by minimizing

J(c) = 1 2 N i=1 c T K(qu)c -E i 2 + 1 2 c T Λc, s.t. c = 1 (2)
where Λ is a diagonal matrix with elements Λ lm = λ NNS D l 2 (l + 1) 2 for the Laplace-Beltrami regularization. The constraint c = 1 is because of S 2 Φ(u)du = 1. In this abstract, we propose Non-local NNSD (NLNNSD) which considers the non-local spatial information and Rician noise. Non-local mean has been used in image denoise [START_REF] Buades | A non-local algorithm for image denoising[END_REF], [START_REF] Descoteaux | Impact of Rician Adapted Non-local Means Filtering on HARDI[END_REF] and regularization [START_REF] Peyré | Non-local regularization of inverse problems[END_REF]. The cost function in NLNNSD is

J({c x }) = 1 2 V x=1 N i=1 (c x ) T K(qu)c x -NLM(E x i ) 2 + 1 2 (c x ) T Λc x + 1 2 λ NLM c x -NLM(c x ) 2 (3)
where c x and E x i are the coefficient vector and diffusion signal at voxel x, V is the number of voxels, NLM(c x ) = arg min c y∈V w y d(c, c y ) 2 is the non-local Riemannian mean of c x [6],

NLM(E x i ) = y∈V p y (E y i ) 2 -2σ 2 is the non-local mean of E x i
considering Rician noise with standard deviation of σ. w y is the nonlocal weights determined by the distance of coefficient vectors, i.e.

w y = 1 Zy exp(- j∈Nx ,k∈Ny Ga c j -c k 2 2h 2
), where c j and c k are the coefficient vectors respectively in the neighborhood N x of x and the neighborhood N y of y, G a is the Gaussian weighting with standard deviation of a, and Z y is the normalization factor. p y is the non-local weight determined by the distance of {E x i } with another set of {a, h}. To minimize Eq. ( 3) with the constraint c = 1, we first set λ NLM = 0, and perform a Riemannian gradient descent on the sphere c = 1 [START_REF] Cheng | A Riemannian Framework for Orientation Distribution Function Computing[END_REF] to minimize J(c x ) individually for each voxel x.

(c x ) (k+1) = Exp (c x )(k) -dt ∇J(c x ) ∇J(c x ) , Exp c (v) = c cos v + v v sin v
The isotropic fODF with c = (1, 0, . . . , 0) T is chosen as the initialization. Then the non-local Riemannian mean is performed to calculate NLM(c x ) at each voxel. Then the Riemannian gradient descent is performed again with λ NLM and the estimated non-local mean NLM(c x ) in Eq. ( 3). Note that this procedure can be iteratively performed to update NLM(c x ) and c x , however in practice we found the result with just one iteration is enough.

III. ISBI HARDI Reconstruction Challenge

In the ISBI reconstruction challenge, the testing data was generated based on Numerical Fibre Generation toolbox [START_REF] Close | A software tool to generate simulated white matter structures for the assessment of fibre-tracking algorithms[END_REF]. we test the proposed NLNNSD in the data with three kind of sampling schemes: 1) single shell DTI scheme with 32 directions, b = 1200s/mm 2 ; 2) single shell HARDI scheme with 64 directions, b = 3000s/mm 2 ; 3) multiple shell DSI-like scheme with 514 directions, b ∈ (0, 4000]s/mm 2 . For all night datasets (three schemes with three SNR 10, 20, 30), we fixed L = 8, λ NNS D = 0, λ NLM = 1, and used the tensor fiber response function with FA of 0.8, mean diffusivity of 0.8. In the non-local mean of c x and E x i , we uses a 11 × 11 × 11 search window, and a 3 × 3 × 3 patch to define the weights, where the parameters a and h were tuned respectively for {c x } and {E x i } to obtain visually good results for each dataset.


Jian Cheng

Dinggang Shen

Pew-Thian Yap

Non-Negative Spherical Deconvolution for Fiber Orientation Distribution Estimation

Introduction:

In diffusion MRI, Spherical Deconvolution (SD) was proposed to estimate the fiber Orientation Distribution Function (fODF) based on spherical deconvolution using a single-fiber response function [1,2]. The peaks or the shape of fODFs can be used to infer local fiber directions. Constrained Spherical Deconvolution (CSD) [1], which takes into consideration the non-negative of the fODF, is the most widely used method among SD variants. Although CSD is capable of accurately determining fiber directions, it is susceptible to false positive peaks especially in the regions with low anisotropy. This is a common drawback of all existing SD-based methods. Moreover, in practice the fODF estimated using CSD still has significant negative values. We propose a method called Non-Negative Spherical Deconvolution (NNSD) to solve the above two problems. Based on a Riemannian framework of ODFs [3] and Square Root Parameterized Estimation for non-negative definite Ensemble Average Propagator [4], NNSD is formulated such that the non-negativity of the fODF is guaranteed with largely reduced false positive peaks.

Methods:

We represent both the fODF and the signal from a single fiber along the z-axis (i.e., the kernel) using Spherical Harmonic (SH) basis with the maximal order , i.e.

∑ c Y , and ∑ h Y , , respectively, where is a unit vector and Y • is the SH basis with order l and degree m. Based on the theoretical results in [1,2], the MR signal attenuation is reconstructed using the spherical convolution between the fODF and the

J ∑ E λ s.t. ‖ ‖ 1
where is the Laplace-Beltrami regularization diagonal matrix with element l l 1 . Since fODF is a probability density function and SH is an orthonormal basis, we have the constraint ‖ ‖ 1 [3,4]. Riemannian gradient descent on the sphere is performed to minimize the above cost function [4].

1,0, … ,0 , indicating isotropic fODF, is used for initialization. The fODF is recovered by squaring the square root of the fODF naturally resulting in non-negativity of the fODF.

Experiments and Results:

The proposed NNSD method was validated using both synthetic data and real data. The synthetic DWI data was generated based on a mixture of tensors with eigenvalues 1.7,03,0.3 10 mm /s, with a b-value of 1500s/mm , and with 81 directions non-collinear gradientsensitized directions. This data was perturbed with Rician noise with a low SNR=10.The SNR is defined as the ratio between the standard deviation of the complex Gaussian noise and the baseline signal without diffusion weighting. NNSD were compared with the state-of-the-art CSD using the same fiber kernel estimated from the voxels with FA larger than 0.8. The CSD and the kernel estimation implemented by their authors in MRtrix (http://www.brain.org.au/software/mrtrix/) were used. For fair comparison, we chose for both methods a maximum order of 6 for the SH basis. We intentionally set λ 0 so that NNSD can be compared with CSD with no smoothness regularization. The results are shown in the left two subfigures of Fig. 1. The intensities of the background image represent the Generalized Fractional Anisotropy (GFA) [5] values calculated from the fODFs. The fODF glyphs were directionally colored. It can be observed that in both single-fiber and crossing regions, NNSD and CSD both can correctly detect the local fiber orientations. The fODFs computed by NNSD however do not contain any negative values, as opposed to CSD. In more isotropic regions, the fODFs given by NNSD have significantly lower GFA than those by CSD, signifying the fact that NNSD removes a significant amount of false positive peaks. Evaluation was also performed using real human data with a b-value of 2000 s/mm , 120 gradient directions, 2mm isotropic voxel dimensions, TR/TE=124,000ms/116ms. We used the same parameters as those in the synthetic data for CSD and NNSD. The results are shown in the right two subfigures in Fig. 1.CSD results in a significant amount of false positive peaks, which is especially evident in the regions with low anisotropy, as indicated by the high background GFA values. NNSD obtains a very clean fODF field with similar peaks as detected by CSD in the high anisotropic regions and more isotropic fODFs in the isotropic regions, such as CSF and grey matter.

Conclusion:

We have proposed a novel Spherical Deconvolution (SD) method called Non-Negative SD (NNSD) to overcome the two main limitations of existing fODF estimation methods. First, NNSD guarantees that the estimated fODFs are non-negative across the whole continuous spherical domain. Second, NNSD reduces significantly the false positive peaks that are prevalent in methods like CSD [1]. The clean fODF field given by NNSD will help reducespurious fiber tracts in tractography algorithms.

References: [1]. Tournier, J., et al., Neuroimage, 34(4), 1459-1472, 2007. [2]. Anderson, A., et al., Magnetic Resonance in Medicine, 54(5), 1194-1206, 2005. [3] Cheng, J., et al., MICCAI 2009, LNCS 5761, 911-918. [4] Cheng, J. et al., MICCAI 2012, LNCS, 6892, 98-106. [5] Tuch, D. S., Magnetic Resonance in Medicine, 52, 1358-1372.

Fig. 1 .

 1 Fig. 1. CSD and NNSD results for synthetic data (two left panels) and real data (two right panels). Min-max normalization was not performed.