Sylvain Merlet

Emmanuel Caruyer

Aurobrata Ghosh

Rachid Deriche

A computational diffusion MRI and parametric dictionary learning framework for modeling the diffusion signal and its features

Keywords: Dictionary Learning, Sparse Reconstruction, Compressive Sensing, Compressed Sensing, MRI, Diffusion MRI, Orientation Distribution Function, Ensemble Average Propagator, White Matter

In this work, we first propose an original and efficient computational framework to model continuous diffusion MRI (dMRI) signals and analytically recover important diffusion features such as the Ensemble Average Propagator (EAP) and the Orientation Distribution Function (ODF). Then, we develop an efficient parametric dictionary learning algorithm and exploit the sparse property of a well-designed dictionary to recover the diffusion signal and its features with a reduced number of measurements. The properties and potentials of the technique are demonstrated using various simulations on synthetic data and on human brain data acquired from 7-T and 3-T scanners. It is shown that the technique can clearly recover the dMRI signal and its features with a much better accuracy compared to state-of-the-art approaches, even with a small and reduced number of measurements. In particular, we can accurately recover the ODF in regions of multiple fiber crossing, which could open new perspectives for some dMRI applications such as fiber tractography.

Introduction

Diffusion MRI (dMRI) assesses the integrity of brain anatomical connectivity and is very useful for examining and quantifying white matter (WM) microstructure and organization not available with other imaging modalities. dMRI determines the WM structure by exploiting the way the water molecules diffuse. The first diffusion images were obtained in the mid-1980s [START_REF] Le Bihan | Mr imaging of intravoxel incoherent motions: application to diffusion and perfusion in neurologic disorders[END_REF], and was based on the pioneering work of [START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF], who introduced the pulsed gradient spin-echo (PGSE) sequence. It allows the quantification of the water diffusion by estimating the displacement of water particles from the phase change that occurs during the acquisition process. More importantly, under the so called narrow pulse assumption, we can show that the normalized signal attenuation E(q) is written as the Fourier transform of the Ensemble Average Propagator (EAP) P(R)

E(q) = R∈R 3 P(R) exp(-2πiq • R)dR, (1)
where q and R are both 3D-vectors that respectively represent the effective gradient direction and the displacement direction. We can decompose them as q = qu and R = Rr, where u and r are 3D unit vectors. Diffusion Tensor Imaging (DTI) [START_REF] Basser | Mr diffusion tensor spectroscopy and imaging[END_REF][START_REF] Basser | Estimation of the effective self-diffusion tensor from the nmr spin echo[END_REF] method characterizes the diffusion by a Gaussian distribution, and is known to be a limited model. In particular DTI is not able to resolve crossing fibers. Resolving crossing fibers helps to disambiguate between several possible tracts in regions of crossing fibers and reconstruct more accurate anatomical connectivity through fiber tractography. Recently, more complex models appeared to overcome this limitation. Nevertheless, these techniques often require many acquisitions in particular when High Angular Resolution Diffusion Imaging (HARDI) [START_REF] Tuch | Q-ball imaging[END_REF][START_REF] Anderson | Measurements of fiber orientation distributions using high angular resolution diffusion imaging[END_REF][START_REF] Vega | Estimation of fiber orientation probability density functions in high angular resolution diffusion imaging[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF][START_REF] Jian | A novel tensor distribution model for the diffusion-weighted mr signal[END_REF][START_REF] Aganj | Reconstruction of the ODF in single and multiple shell q-ball imaging within constant solid angle[END_REF] or Diffusion Spectrum Imaging (DSI) [START_REF] Wedeen | Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging[END_REF] are used. HARDI techniques allow the estimation of the Orientation Distribution Function (ODF) [START_REF] Tuch | Q-ball imaging[END_REF][START_REF] Anderson | Measurements of fiber orientation distributions using high angular resolution diffusion imaging[END_REF][START_REF] Vega | Estimation of fiber orientation probability density functions in high angular resolution diffusion imaging[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF][START_REF] Aganj | Reconstruction of the ODF in single and multiple shell q-ball imaging within constant solid angle[END_REF], which gives the probability that a water molecule diffuses in a given direction. Several authors (Aganj et (

In [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], the authors propose to estimate the fiber orientation distribution called the fiber ODF (fODF). The fODF is able to resolve up to 30 degrees crossings consistently [START_REF] Tournier | Resolving crossing fibres using constrained spherical deconvolution: Validation using diffusion-weighted imaging phantom data[END_REF], which makes this model a promising ressource to estimate fiber orientation. In our work, we reconstruct the ODF as described in Eq. 2 and we do not aim to compare our approach with the family of method estimating the fODF. A review of method reconstructing the ODF and the fODF can be found in [START_REF]Diffusion MRI : from quantitative measurement to in-vivo neuroanatomy[END_REF]. Another HARDI technique has been proposed in [START_REF] Jian | A novel tensor distribution model for the diffusion-weighted mr signal[END_REF], where the authors characterize the diffusion signal by a Wishart distribution. Jian et al. shows improvements over the classical DTI technique and present an estimation scheme for the fiber orientation and EAP. Among the HARDI techniques, [START_REF] Zhang | Noddi: Practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF] introduces Neurite Orientation Dispersion and Density Imaging (NODDI), which allows the estimation of the microstructural complexity of dendrites and axons. Diffusion Spectrum Imaging (DSI) was developed in parallel to the HARDI techniques [START_REF] Wedeen | Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging[END_REF]. In DSI, the EAP P(R) is directly obtained by taking the inverse Fourier transform of the normalized signal E(q) measured in the q-space (see Eq. 1). However, the high resolution EAP obtained with DSI requires many measurements. HARDI and DSI are impractical for clinical use in MRI systems commonly found in hospital. Accelerated acquisitions, relying on a smaller number of sampling points, are thus very welcome to efficiently estimate the complex features of the diffusion process. Sparse reconstruction approaches were found to successfully reduce the number of acquisitions in dMRI [START_REF] Merlet | Compressed sensing for accelerated eap recovery in diffusion mri[END_REF][START_REF] Menzel | Accelerated diffusion spectrum imaging in the human brain using compressed sensing[END_REF][START_REF] Michailovich | Spatially regularized compressed sensing for high angular resolution diffusion imaging[END_REF][START_REF] Rathi | Sparse multi-shell diffusion imaging[END_REF][START_REF] Vega | Probabilistic odf estimation from reduced hardi data with sparse regularization[END_REF][START_REF] Bilgic | Accelerated diffusion spectrum imaging with compressed sensing using adaptive dictionaries[END_REF][START_REF] Gramfort | Sparse dsi: Learning dsi structure for denoising and fast imaging[END_REF][START_REF] Merlet | Parametric dictionary learning for modeling eap and odf in diffusion mri[END_REF][START_REF] Merlet | Continuous diffusion signal, eap and odf estimation via compressive sensing in diffusion mri[END_REF][START_REF] Ye | An over-complete dictionary based regularized reconstruction of a field of ensemble average propagators[END_REF]. These techniques are usually based on a l 1 minimization of the diffusion signal with respect to a sparse representation. Merlet and Deriche, Menzel et al. [START_REF] Merlet | Compressed sensing for accelerated eap recovery in diffusion mri[END_REF][START_REF] Menzel | Accelerated diffusion spectrum imaging in the human brain using compressed sensing[END_REF] combine the Compressive Sensing (CS) theory and DSI to accelerate the acquisition. Merlet and Deriche [START_REF] Merlet | Continuous diffusion signal, eap and odf estimation via compressive sensing in diffusion mri[END_REF] use orthonormal bases to sparsely describe the diffusion signal. In [START_REF] Michailovich | Spatially regularized compressed sensing for high angular resolution diffusion imaging[END_REF] and in [START_REF] Vega | Probabilistic odf estimation from reduced hardi data with sparse regularization[END_REF], the authors elegantly design dictionaries for sparse modeling in dMRI. They provide an overcomplete dictionary computed from a discretized version of predefined functions, i.e. the Spherical Ridgelets in [START_REF] Michailovich | Spatially regularized compressed sensing for high angular resolution diffusion imaging[END_REF] (see [START_REF] Rathi | Sparse multi-shell diffusion imaging[END_REF] for the multiple shells version) and the Spherical Wavelets in [START_REF] Vega | Probabilistic odf estimation from reduced hardi data with sparse regularization[END_REF]. Learning a dictionary provides an alternative way to design sparse dictionaries [START_REF] Bilgic | Accelerated diffusion spectrum imaging with compressed sensing using adaptive dictionaries[END_REF][START_REF] Gramfort | Sparse dsi: Learning dsi structure for denoising and fast imaging[END_REF][START_REF] Merlet | Parametric dictionary learning for modeling eap and odf in diffusion mri[END_REF][START_REF] Ye | An over-complete dictionary based regularized reconstruction of a field of ensemble average propagators[END_REF].

Some approaches have been recently proposed in order to design dictionaries that enable sparse representations (A good overview can be found in [START_REF] Aharon | K-svd: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF]). For instance, Bilgic et al., Gramfort et al. [START_REF] Bilgic | Accelerated diffusion spectrum imaging with compressed sensing using adaptive dictionaries[END_REF][START_REF] Gramfort | Sparse dsi: Learning dsi structure for denoising and fast imaging[END_REF] learn dictionaries from DSI like acquisitions and use it to either denoise full DSI data or to perform undersampled DSI acquisitions and reconstructions. In particular, Gramfort et al. nicely exploit the symmetry of the signal in order to assess free parameters of the dictionary learning problem. However, these two latter works lead to non-parametric dictionaries, which does not provide continuous representations of the diffusion signal nor allow the determination of analytical formulae for diffusion features. The strength of the parametric dictionary learning approach, as the one we propose in this article, lies in its ability to address these weaknesses. A work regarding parametric dictionary learning was published in [START_REF] Ye | An over-complete dictionary based regularized reconstruction of a field of ensemble average propagators[END_REF], in which the dictionary atoms are formed by a weighted combination of 3 rd order B-splines. It proved that the method is efficient on synthetic data simulated with 81 gradient directions. The work of [START_REF] Ye | An over-complete dictionary based regularized reconstruction of a field of ensemble average propagators[END_REF] appears promising in reconstructing the diffusion signals, and further enhancement could be done regarding the development of analytical formulae to estimate other diffusion features. This would make this work a good resource in the context of dictionary learning. More recently, we proposed in [START_REF] Merlet | Parametric dictionary learning for modeling eap and odf in diffusion mri[END_REF] to learn a dictionary where each atom is constrained to be a parametric function. In [START_REF] Merlet | Parametric dictionary learning for modeling eap and odf in diffusion mri[END_REF], this parametric function is a combination of a radial part and an angular part represented by the symmetric and real Spherical Harmonics (SH) [START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF]. The radial part is a polynomial weighted by an exponential. 50 measurements were sufficient to reconstruct very good quality diffusion signals, ODFs and EAPs. However, this approach essentially handles the learning of the radial part, i.e. the polynomial coefficients and a scale parameter in the exponential, whereas we observed (see [START_REF] Merlet | Parametric dictionary learning in diffusion mri[END_REF]) that the angular part could make the dictio-nary much sparser if we adequately combine several SH functions instead of only one.

In this work, we present a method, which exploits the sparse property of a well designed dictionary based on a computational dMRI framework, in order to recover the diffusion signal with a reduced number of measurements. This framework enables a continuous modeling of the diffusion signal and leads to analytical formulae to estimate important diffusion features, namely the ODF and the EAP. To improve our previous work in [START_REF] Merlet | Parametric dictionary learning for modeling eap and odf in diffusion mri[END_REF], we modify the parametric function, describing the atoms, to learn both the radial and the angular part, which provide a very sparse representation of diffusion signals and further reduce the number of measurements (15 measurements are found to be sufficient to start recovering the EAP and some derived diffusion features whereas 50 measurements are used in [START_REF] Merlet | Parametric dictionary learning for modeling eap and odf in diffusion mri[END_REF]). Furthermore, we extend the experimental part of [START_REF] Merlet | Parametric dictionary learning for modeling eap and odf in diffusion mri[END_REF] by learning and validating our approach on the synthetic data proposed in the HARDI contest at ISBI 2012 1 , and on real data acquired from both 3T and 7T scanners. A preliminary work [START_REF] Merlet | Parametric dictionary learning in diffusion mri[END_REF] regarding the learning of both the radial part and the angular part of the diffusion signal was published in the proceedings of the HARDI contest at ISBI 2012 and we obtained the best results in our category. Our approach presented in this paper indicates an increase in terms of reconstruction accuracy compared to the results presented in [START_REF] Merlet | Parametric dictionary learning in diffusion mri[END_REF].

The article is structured as follows : we start by introducing the dMRI framework together with the proposed dictionary, then we focus on the parametric dictionary learning algorithm and finally we conclude with an experimental part illustrating the added-value of our approach with promising results showing how our approach allows the accurate reconstruction of the diffusion signal and some of its features. This experimental part is completed by a comparison with state of the art approaches, and is performed on synthetic and real data from 3T and 7T scanners.

A computational framework for the recovery of the complete diffusion MRI process

In this section, we introduce a new dMRI framework for modeling the diffusion signal. From this continuous representation, we derive analytical formulae that enable the estimation of important diffusion features such as the Ensemble Average Propagator (EAP), the Orientation Distribution Function (ODF). We give full derivations for these formulae in the appendixes.

Continuous diffusion modeling with a constrained dictionary

We propose to design an overcomplete dictionary Ψ = {Ψ k } k=0,...,K , such that the diffusion signal E is expressed as a truncated linear combination of K 3D atoms Ψ k , i.e

E(qu) = K k=0 c k Ψ k (qu), (3)
with c k the transform coefficients, q the norm of the effective gradient and u an unitary vector.

In diffusion multiple shells imaging [START_REF] Assemlal | Efficient and robust computation of pdf features from diffusion mr signal[END_REF][START_REF] Ozarslan | Simple harmonic oscillator based reconstruction and estimation for threedimensional q-space mri[END_REF][START_REF] Descoteaux | Multiple q-shell diffusion propagator imaging[END_REF], previous works have proposed the modeling of the diffusion signal in bases where each atom Ψ k is a combination of a radial part and an angular part. These works primarily handle the question of the radial part. For instance, Descoteaux et al. [START_REF] Descoteaux | Multiple q-shell diffusion propagator imaging[END_REF] describes the radial part of the atom as a combination of two rational functions coming from the total solution of the Laplace equation, whereas Ozarslan et al., Assemlal et al. [START_REF] Ozarslan | Simple harmonic oscillator based reconstruction and estimation for threedimensional q-space mri[END_REF][START_REF] Assemlal | Efficient and robust computation of pdf features from diffusion mr signal[END_REF] use polynomials weighted by an exponential. Merlet et al. [START_REF] Merlet | Parametric dictionary learning for modeling eap and odf in diffusion mri[END_REF] increased the sparsity of the representation in [START_REF] Ozarslan | Simple harmonic oscillator based reconstruction and estimation for threedimensional q-space mri[END_REF] by designing a dictionary where the polynomial coefficients and the scale parameters in the exponentials are learned from a training data set [START_REF] Merlet | Parametric dictionary learning for modeling eap and odf in diffusion mri[END_REF]. Although the radial attenuation of the diffusion signal is more or less well fitted with these bases/dictionaries, the accurate and sparse estimation of the full diffusion signal is still limited by the angular part described by the real and symmetric Spherical Harmonic basis functions (SH). The SH have been proved useful in many settings but are not sparse enough in modeling the directional features of the diffusion process [START_REF] Michailovich | On approximation of orientation distributions by means of spherical ridgelets[END_REF]. However, our preliminary work published in [START_REF] Merlet | Parametric dictionary learning in diffusion mri[END_REF] shows that a well-chosen combination of SH could sparsely represent the angular part. Therefore, we propose to model this angular part with such a combination of SH. As for the radial part, we model it with a combination of exponential functions weighted by a monomial in order to ensure the continuity of the function at zero [START_REF] Caruyer | Diffusion MRI signal reconstruction with continuity constraint and optimal regularization[END_REF]. The complete description of each atom Ψ k of the dictionary is given by :

Ψ k (q) = Ψ k (qu) = 1 √ χ k I i=0 α ki exp -ν ki q 2 J j=0 β ki j q l( j) Y j (u) = 1 √ χ k I i=0 J j=0
γ ki j exp -ν ki q 2 q l( j) Y j (u)

= Ψ k (γ k , ν k , qu), (4)
with q the 3D effective gradient, u an unitary vector and q the norm of the effective gradient such that q = qu. I and J are, respectively, the radial order and the angular order of the dictionary. J also corresponds the total number of SH taken into account in the modeling not to be confused with the maximal SH order L. Indeed, J is directly related to the maximal SH order L as J = (L + 1)(L + 2)/2. Y j (u) is the SH of order l( j) = 0 for j = 1, l( j) = 2 for j ∈ {2, ..., 6}, l( j) = 4 for j ∈ {7, ..., 15} ... . γ k = {γ ki j } i=0...I, j=0...J and ν k = {ν ki } i=0...I are two vectors of parameters, which will be set during the learning process. The term q l( j) ensures the continuity of Ψ k at zero. χ k is a constant, which ensures the normalization of Ψ k for the ℓ 2 norm, i.e R 3 Ψ 2 k (q)dq = 1, and is expressed as

χ k = I i ′ =0 I i=0 J j=0 γ ki j γ ki ′ j 2(ν ki + ν ki ′ ) ℓ( j)+3/2 Γ ℓ( j) + 3 2 , (5)
with Γ the gamma function. We derive eq. 5 in Appendix A.

If we consider the simple case where all the coefficients γ ki j and ν k j are zero but γ k00 and ν k0 , we get

Ψ k (qu) = 1 √ χ k γ k00 2 √
π exp(-ν k0 q 2 ). This mono-exponential representation of the atom Ψ k does not depend on the angular direction u and is, for instance, appropriate to describe isotropic diffusion configuration.

Note that in [START_REF] Merlet | Parametric dictionary learning in diffusion mri[END_REF] each atom is described by a combination of SHORE basis functions [START_REF] Ozarslan | Simple harmonic oscillator based reconstruction and estimation for threedimensional q-space mri[END_REF] with a predefined scale parameter. The main advantage of the atom description in Eq. 4 lies in the possibility to learn the scale parameters ν ki , which provides a sparser signal estimation than in [START_REF] Merlet | Parametric dictionary learning in diffusion mri[END_REF]. We see in the experiment part that this new framework leads to an increase in terms of reconstruction accuracy compared to the results presented in [START_REF] Merlet | Parametric dictionary learning in diffusion mri[END_REF], which were already the best in their category.

Closed formulae for diffusion features

Using the dictionary, Ψ = {Ψ k } k=0,...,K , proposed in the previous section to reconstruct the diffusion signal (see Eq. 4), we derive important and analytical closed formulae for estimating the EAP and the ODF. We describe these formulae in the following.

The Ensemble Average Propagator

The EAP, denoted P(Rr), represents the full 3D displacement probability function of water molecules in every voxel and underlies the derivation of the ODF. It is the inverse Fourier transform of the normalized diffusion signal, denoted E(qu), P(Rr) = ∞ q=0 u∈S 2 E(qu) exp(+2πiqRu.r)duq 2 dq. [START_REF] Basser | Estimation of the effective self-diffusion tensor from the nmr spin echo[END_REF] q and R are, respectively, the norm of the effective gradient and the radius of the 3D location in every voxel, u and r are unit vectors. From eq. 3 and 6, we derive in Appendix B the following expression for the EAP :

P(Rr) = K k=0 c k √ χ k I i=0 J j=0 γ ki j (-1) l( j)/2 π ν ki l( j)+3/2 R l( j) exp -(πR) 2 ν ki Y j (r) (7)
If we consider the special case of the monoexponential representation of the atom Ψ k , we get

P(Rr) = K k=0 c k √ χ k γ k00 2 √ π π ν k0 3/2 exp -(πR) 2 ν k0
, i.e an isotropic propagator described by an mono-exponential decay similar in every direction.

Solid angle ODF

The ODF represents the full angular distribution of P(Rr). One relies on the ODF to perform fiber tractography [START_REF] Descoteaux | Deterministic and probabilistic tractography based on complex fibre orientation distributions[END_REF], then an accurate and fast computation of this diffusion feature is very appreciated. From Eq. 2, we derive in Appendix C the following closed form for the ODF :

Υ(r) = K k=0 c k √ χ k I i=0 J j=0 γ ki j (-1) l( j)/2 π ν ki l( j)+1 (l( j) + 1)!! 2 2π 2 v ki l( j)/2+1 Y j (r) (8)
Considering the precedent example, the ODF is represented by a scalar, i.e.

Υ(r) = K k=0 c k √ χ k γ k00
8π 3/2 , meaning that the ODF has the same value for every direction.

A parametric dictionary learning for sparse dMRI

Here, we introduce a parametric dictionary learning (PDL) method that enables a sparse representation of any diffusion signal from continuous and parametric functions. There are four advantages to consider a parametric approach for dictionary learning:

• A parametric dictionary is defined by a set of parameters (γ k and ν k in Sec. 2), which gives a continuous representation of each atom and, thus, enables a continuous modeling of the diffusion signal. This is suitable for data interpolation and extrapolation.

• Analytical formulae can be derived to estimate important diffusion features as the EAP and the ODF.

• PDL is acquisition independent, i.e. the sampling scheme used for learning the dictionary does not have to be the same as the sampling scheme used for reconstructing the signals.

• PDL enables one to reduce the dimensionality of the dictionary atoms.

These four advantages, together with the quality of the obtained results (see Sec. 4), makes our parametric dictionary approach very attractive compared to non parametric methods [START_REF] Bilgic | Accelerated diffusion spectrum imaging with compressed sensing using adaptive dictionaries[END_REF][START_REF] Gramfort | Sparse dsi: Learning dsi structure for denoising and fast imaging[END_REF].

Concerning the development of our algorithm, we started by considering the K-SVD algorithm [START_REF] Aharon | K-svd: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF] as a model for our own method. Although the K-SVD method appears powerful in designing sparse dictionaries, this technique as described in [START_REF] Aharon | K-svd: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF] designs only non-parametric dictionaries, which do not present the advantages described above. Hence, we developed our own algorithm, which overcomes the limitation of the K-SVD algorithm. Our algorithm alternates between 2 steps: a sparse coding step and a dictionary update step, where the vectors of parameters γ k and ν k (see Sec. 2) are estimated for every atom d k of the dictionary, using the non linear Levenberg-Marquardt (LM) algorithm. The section 3.1 presents our dictionary learning algorithm and the section 3.2 describes the method we use to reconstruct any diffusion signal using the dictionary previously learned.

Dictionary learning algorithm

Notation : Suppose the training data set consists in M observations {s i } M i=1 (i.e. M voxels). For each observation s i we have m s samples in the q-space, i.e. s 1..M ∈ R m s . We represent {s i } M i=1 in matrix form S ∈ R m s ×M where s i is the i th column. The algorithm searches for the dictionary D ∈ R m s ×K , that enables the sparsest representation for every column of S. The dictionary consists in K atoms {d k } K k=1 with d k ∈ R m s a column of D. We constrain d k to be an instance of the 3D function Ψ k (γ k , ν k , qu) in Eq. 4. Here, we do not try to directly estimate d k but the vectors of parameters γ k and ν k , that characterize the atom d k . For each observation s i , we define a coefficient vector c i ∈ R K , which forms the i th column of the coefficient matrix C ∈ R K×M .

Problem statement: Given a training data set S, we search for the dictionary D that gives the sparsest representation of this set (i.e. for each column s i of S). Mathematically, the problem is to find the dictionary D and the vectors c i in C by solving :

arg min c i ,D { S -DC 2 2 } subject to ∀ i c i 1 ≤ ǫ (9)
with ǫ a small real defining the degree of sparsity of the dictionary. The minimization of the first term in Eq. 9 enables the signal estimation DC to remain close to the training data set S and the constraint imposes the sparsity of each signal representation c i with respect to the dictionary D. The method to solve Eq. 9 is described in the following and a summary of the algorithm is given in Alg. 1. This algorithm iteratively alternates between sparse signal estimations (i.e. {c i } M i=1 ) and updates of the dictionary (i.e. D) so to better fit the training data set (i.e.

S).

First step (Sparse signal estimation): In the first step, the estimation of the column vector c i is performed separately for each signal s i , i.e for each column of S. Sparse estimation is achieved by solving the LASSO (Least Absolute Shrinkage and Selection Operator) problem [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. It consists in minimizing the following objective function min

c i s i -Dc i 2 2 + λ c i 1 . (10)
λ is a constant that controls the degree of sparsity of the coefficients estimated. Note that we relax the constraint in Eq. 9 by using a Lagrangian multiplier (λ). There exist numerous iterative algorithms for efficiently solving such kind of constrained problems. These include coordinate descent, least-angle regression (LARS) [START_REF] Efron | Least angle regression[END_REF], fast iterative thresholding shrinkage algorithm (FISTA) [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], etc. A number of these methods are available under the Python TM library Scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF]. We use a Python TM implementation [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF] of coordinate descent to solve Eq. 10.

Second step (Dictionary update): In the second step, we update the dictionary D. For this purpose, we compute an absolute averaged coefficient vector ĉ ∈ R n c , such that ĉ = 1/M i |c i | (|.| denotes the absolute value of each vector component), and find the atoms associated with the non zeros values of ĉ. It gives a rough idea of which atoms are used for modeling the signal and enables one to discard some unnecessary atoms and, thus, to enforce sparsity. Then, in this set of atoms, we update one atom at a time, while fixing all the others. This process is repeated for all the atoms associated with the non-zero coefficients of ĉ.

The in-update atom is denoted d k 0 . To update this atom, we begin by decomposing the error term in eq. 9 as in [START_REF] Aharon | K-svd: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF], i.e.

S -DC

2 ℓ2 = S - K k=1 d k c r k 2 ℓ2 =         S - k k 0 d k c r k         -d k 0 c r k 0 2 ℓ2 = E k 0 -d k 0 c r k 0 2 ℓ2 , (11)
where c r k is the k th row of C. The error matrix, denoted E k 0 , contains the error between each observation s i (the i th column of S) and its respective estimation with the dictionary where the k th 0 atom is removed. We could directly use the LM algorithm in order to fit the atom d k 0 to the error matrix E k 0 . However, because it takes into account all the observations s i , this dictionary update would not impose sparsity. Instead, we enforce the sparsity by constraining the atom d k 0 to fit only a subset of observations and not the entire data set. For this purpose, we define the group of observations that use the atom

d k 0 , i.e w k 0 = {i, 1 ≤ i ≤ M, c i (k 0 )
0}. In other words, they are the observations whose coefficients, associated with the atom d k 0 , are non zeros. Then, we compute the error matrix E w k 0 ∈ R m s ×card(w k 0 ) . It corresponds to the estimation error between the observation vector {s i } i∈w k 0 that forms the columns of S w k 0 ∈ R m s ×card(w k 0 ) and the signal estimated for the group of observation w k 0 (The k 0 th atom is still removed from the dictionary), i.e Sw

k 0 = k k 0 d k c r k (i), i ∈ w k 0 . Mathematically speak- ing, we have E w k 0 = S w k 0 -k k 0 d k c r k (i), i ∈ w k 0 .
Finally, we estimate the vector of parameters γ k and ν k by constraining d k 0 to fit the error matrix E w k 0 . This part is performed via the non linear Levenberg-Marquardt algorithm (LMA) . The atom update procedure is repeated for every atom d k 0 associated with the non zeros coefficients of ĉ.

The method is given in Algorithm 1 as a whole.

Convergence : The sparse coding step (Eq. 10) is well known to be convex and the coordinate descent algorithm allows one to converge to the unique solution specific to the current dictionary D. The dictionary update step, where γ k and ν k are estimated using the LMA may converges to a local minima, depending on the initial solutions. Then, Eq. 9 is convex for c and converge for D, which do not ensure a convergence to a global minimum. Nevertheless, in our experiments, a stationary point has been reached after few iterations and the resulting dictionaries were proved very good experimentally.

Initialization : The problem in Eq. 9 admits local minima, and the solution may vary depending on the initial parameters.

Algorithm 1 Semi-parametric dictionary learning 1. Initialize the dictionary by fixing its dimension K and the vectors of parameters γ k and ν k for k = 1...K as random.

2. Sparse estimation of the observations {s i } M i=1 . We use the coordinate descent algorithm to solve for c i associated to each observation :

min c i s i -Dc i 2 2 + λ c i 1 .
3. Updating the dictionary. Compute the absolute averaged coefficients vector ĉ = 1/M i |c i |. Repeat until all the atoms of the dictionary, with non zeros value in ĉ, have been scanned :

• Let note the current atom, the k th 0 .

• Define the group of observation that use this atom :

w k 0 = {i, 1 ≤ i ≤ M, c k 0 (i) 0}.
• Compute the error matrix

E w k 0 = S w k 0 -k k 0 d k c r k (i), i ∈ w k 0 . S w k 0 contains the observation s i , i ∈ w k 0 .
• Apply the Levenberg -Marquardt algorithm to estimate the vectors of parameters γ k 0 and ν k 0 , which constrain d k 0 to best fit E w k 0

• Update d k 0 according to γ k 0 and ν k 0 . 4. Go back to the step 2 unless the overall error does not vary anymore

We tried several ways to initialize the algorithm among which random initialization and initialization from signals selected at random in the training data set.. However, these kinds of initialization were not satisfactory since the corresponding solutions were too different between several attempts to build the dictionary. After many experiments, we finally selected each initial atom as a random combination of several training signals. This gave us the best satisfactory results, with the smallest sensitivity to the initialization.

Signal estimation via the learned dictionary

The purpose of section 3.1 was to design a parametric dictionary D. Now, we are able to recover any sparse representation c of diffusion signal s regarding the dictionary D by solving the LASSO problem :

min c s -Dc 2 2 + λ c 1 . (12)
We use the same algorithm as in the learning step to solve Eq. 10, i.e. the coordinate descent algorithm.

Experiments on synthetic data

We first train and validate our parametric dictionary on synthetic data. We assume the normalized diffusion signal E(q) is generated from the multi-tensor model for F fibers,

F f =1 p f exp(-4π 2 τq 2 u T T f u), (13)
where a fibre f is defined by a tensor matrix T f and weight p f , such that f p f = 1. q denotes the norm of the effective gradient and u is a unitary vector in Cartesian coordinate.

The analytical ground truth of the EAP for any radius R is then given by

P(Rr) = F f =1 p f 1 (4πτ) 3 |T f | exp -R 2 r T T f -1 r 4τ , (14)
with r a unitary vector in Cartesian coordinate.

We can also derive the ODF feature using the solid angle closed form expression [START_REF] Aganj | Reconstruction of the ODF in single and multiple shell q-ball imaging within constant solid angle[END_REF][START_REF] Vega | Estimation of fiber orientation probability density functions in high angular resolution diffusion imaging[END_REF],

ψ(r) = F f =1 p f 1 4π|T f | 1 2 (r T T f -1 r) 3 2 . ( 15
)
In the remainder of this section, we describe the steps required to correctly design our parametric dictionary, i.e. the choice of the dictionary radial and angular orders and the learning phase, then we validate the learned dictionary on synthetic data using the analytical formulae we have just described.

Which radial and angular order for the dictionary ?

We need to fix the angular and radial order, respectively denoted J and I in Eq. 4, for the dictionary generation.

We begin by defining the dictionary angular order J, which is related to the SH order L as J = (L + 1)(L + 2)/2. For this purpose, we generate synthetic ODFs (see Eq. 15) corresponding to two fibers crossing at different degrees : 0 • , 40 • , 60 • , 90 • . Then, we fit each ODF with the Spherical Harmonic basis of order L = 0, 2, 4, 6, 8, 10, 12 using a least squares fitting technique, and we compute the Normalized Mean Square Error (NMSE) between the original synthetic ODFs (ψ) and the estimated ODFs in terms of SH (ψ e ), i.e.

N MS E = ψ-ψ e 2 2 ψ 2 2
. The resulting NMSEs are shown in Fig. 1. Regarding this figure, we find that a SH order of L = 8 is sufficient to correctly estimate an ODF (In particular for fibers crossing at 40 • or more). This leads to set the dictionary angular order to J = 45.

The radial order is related to the number of fiber crossing in the voxel of interest. We assume a maximum of three crossing fibers in each voxel. This is the case in the region where the corticospinal, corpus callosum and superior longitudinal fasciculus fiber bundles are crossing. Voxel with more than three fibers are considered as part of the noisy background. Consequently, setting the radial order at three appears to be the best choice and was found experimentally satisfactory. We show in Fig. 2 some examples of radial attenuation in arbitray chosen directions simulated using a multi-Gaussian model mimicking three crossing fibers (plain line). We fit each of these radial attenuation with the radial part of the model in Eq. 4 and radial orders I = 1, 2, 3, using the Levenberg-Marquardt algorithm (The angular order J is set to zero). Fig. 2 indicates that we need a radial order I = 3 to accurately estimate radial attenuation characterized by three crossing fibers.

Training phase

We train the dictionary on multi-Gaussian signals, used for the HARDI contest at ISBI 2012 2 . The contest was organized with the aim to provide a way for different groups to propose their own reconstruction algorithms and to fairly compare their methods against the others on a common set of ground-truth data.

Our training data set S train is composed of the first M = 5000 instances of diffusion signal contained in the file TestingIV, in which the multi-Gaussian synthetic signals are generated with parameters taken at random (number of fibers, fractional anisotropy related to a fiber and crossing angle between these fibers). The dictionary angular and radial orders are respectively set to J = 45 and I = 3 (see Sec. 4.1). We take m s = 1000 q-space samples for each instance of signal spread between b min = 0 and b max = 10000s/mm 2 .

One difficulty in dictionary learning is the choice of the regularization parameter λ in Eq. 9. In order to assess λ, we use a cross validation (CV) procedure. For this purpose, we consider another set of signals S val , called the validation data set, and composed of 1000 signals, which have not been used for training, and we repeat the following procedure for a range of λ,

1. Design a dictionary D λ using Algorithm 1 with regularization parameter λ and the training data set S train .

2. Using D λ , solve the LASSO problem for c (see Eq. 12) with S val as entry, and compute the validation error

ǫ val = S val -D λ c 2 2 / S val 2 2
We keep λ that minimizes ǫ val . This procedure is repeated after adding Rician noise, with SNR=10, 20 and 30, to the validation set. Rician noise is added in the following way We show in Fig. 3, ǫ val for λ in the range [1 • 10 -5 , 1 • 10 -4 ], in case of noisy and noiseless validation data. We also show the training error (blue curve in Fig. 3), i.e. ǫ train = S train -

D λ c 2 2 / S train 2
2
where c is the solution of the LASSO problem with S train as entry. We observe, in Fig. 3, four different λ (2.9 • 10 -5 , 4.4 • 10 -5 , 4.8 • 10 -5 , 5.8 • 10 -5 ) that minimize the validation error depending on the amount of noise we add to the validation data set. Because the noiseless case is not observed in practice, we discard the corresponding value of λ and consider the average value of (4.4 • 10 -5 , 4.8 • 10 -5 , 5.8 • 10 -5 ), which gives λ = 5 • 10 -5 . This value is used to generate a new dictionary. We obtain a dictionary containing 659 atoms. Fig. 4 shows the first 200 ODFs of these atoms. The atom ODFs are sorted in decreasing energy order from left to right and top to bottom. We observe various shape ranging from single fiber structures to more complex fiber configurations.

Validation

We validate the dictionary on the reconstruction of noisy multi-Gaussian signals, used for the HARDI contest at ISBI 2012 3 .

We consider 1000 signals, which have not been used for training the dictionary. Our preliminary and promising results on parametric dictionary learning (PDL) were published in the proceedings of this event [START_REF] Merlet | Parametric dictionary learning in diffusion mri[END_REF]. At this stage, we obtained the best results in our category. Here, we enrich these previous results with a comparison of three different sampling schemes (displayed in Fig. 5) :

• A single shell sampling scheme with 64 measurements uniformly spread on one shell at a b-value b = 3000 s • mm -2 (Fig. 5a).

• A multiple shells sampling scheme with 15 measurements spread on 2 shells at b-values b = 1500, 2500 s • mm -2 (Fig. 5b).

• A multiple shells sampling scheme with 64 measurements spread on 2 shells at b-values b = 1500, 2500 s • mm -2 (Fig. 5c).

To obtain the single shell (SS) sampling scheme, we use the algorithm given in [START_REF] Jones | Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging[END_REF][START_REF] Deriche | Optimal real-time q-ball imaging using regularized kalman filtering with incremental orientation sets[END_REF] to uniformly distribute points on a sphere. For the multiple shells (MS) sampling schemes, we use the algorithm given in [START_REF] Caruyer | Optimal Design of Multiple Q-shells experiments for Diffusion MRI[END_REF] by setting the parameters in such a way that the the number of points on each shell is proportional to q 1 . These particular parameters have been proved efficient in [START_REF] Caruyer | Optimal Design of Multiple Q-shells experiments for Diffusion MRI[END_REF][START_REF] Merlet | Impact of radial and angular sampling on multiple shells acquisition in diffusion mri[END_REF]. An important advantage of this algorithm is that the points from each shell have staggered directions and follow a near-optimal uniform distribution.

You can generate and download sampling schemes for multiple Q-shell diffusion MRI with this web application : http://www-sop.inria.fr/members/Emmanuel.Caruyer/q-spacesampling.php. We perform the experiments with two metrics used in the contest, i.e. the weighted difference in the number of fiber compartments (DNC) and the mean angular error (AE) at each voxel. For these two metrics we extract the maxima on the estimated ODFs and compare them to the ground truth maxima. Then, the DNC becomes the difference between the number of maxima extracted on the estimated ODFs M e and the true number of maxima M gt , weighted by the true number of maxima at each voxel, i.e., DNC = |M e -M gt | M gt . The AE is the mean angular error between the maxima extracted on the estimated ODFs and the respective maxima within the ground truth. We also enrich the results with a comparison of two other features, which have not been used in the contest validation, i.e. the diffusion signal and the Ensemble Average Propagator (EAP). To compare these features, we compute the Normalized Mean Square Error (NMSE) between the ground truth feature x and its estimation x e given by N MS E =

x-x e 2 2
x 2 2 . The DNC, AE and NMSE are, then, averaged on all the voxels. We add Rician noise to the normalized diffusion signal in the following way :

E noisy = (E + ǫ 1 ) 2 + ǫ 2 2 , where ǫ 1 , ǫ 2 ∼ N(0, σ) with σ = 1/S NR.
The three following sections present the results for the three proposed sampling schemes, i.e. the SS sampling scheme with 64 measurements, the MS sampling scheme with 15 measurements, and the MS sampling scheme with 64 measurements. We also compare with state of the art techniques such that QBI, using the solid angle ODF [START_REF] Aganj | Reconstruction of the ODF in single and multiple shell q-ball imaging within constant solid angle[END_REF], and the SHORE reconstruction using a sparse prior (ℓ 1 -SHORE) [START_REF] Ozarslan | Simple harmonic oscillator based reconstruction and estimation for threedimensional q-space mri[END_REF][START_REF] Cheng | Theoretical analysis and practical insights on eap estimation via a unified hardi framework[END_REF][START_REF] Merlet | Diffusion and multiple orientations from 1.5 mr systems with limited gradient tables[END_REF][START_REF] Merlet | Ensemble average propagator reconstruction via compressed sensing: Discrete or continuous bases ?[END_REF][START_REF] Merlet | Parametric dictionary learning for modeling eap and odf in diffusion mri[END_REF]. We show quantitative results in Tab. 1, 2, 3 and 4, and qualitative results in Fig. 6, 7 and 8. In each table, we write in blue letters the best score for a given SNR and metric.

Single shell sampling scheme with 64 measurements

In these experiments, we first use the SS sampling scheme, i.e. 64 measurements uniformly spread on a shell at a b-value b = 3000 s • mm -2 , and we compare our dictionary based ODF estimation (D-ODF) to the solid angle ODF via QBI (SA-ODF) developed in [START_REF] Aganj | Reconstruction of the ODF in single and multiple shell q-ball imaging within constant solid angle[END_REF]. For the SA-ODF, we set a SH order equal to the one used for the dictionary construction, i.e. L = 8. We adjust the Laplace-Beltrami regularization parameter [START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF][START_REF] Aganj | Reconstruction of the ODF in single and multiple shell q-ball imaging within constant solid angle[END_REF] using the generalized cross validation algorithm [START_REF] Craven | Smoothing noisy data with spline functions[END_REF]. For our PDL approach, because we deal with ℓ 1 norm and not ℓ 2 norm, we use a simple cross validation procedure [START_REF] Ward | Compressed sensing with cross validation[END_REF] to find the regularization parameter λ in Eq. 12.

Overall, in Fig. 6 the D-ODFs are sharper than the SA-ODFs. Furthermore, the SA-ODFs appear more sensitive to noise than the ODFs based on our PDL estimation. Indeed, we observe than the D-ODFs are very robust to noise, even at SNR=10 where they are still correctly aligned with the underlying structure shown by the ground truth (On the right of Fig. 6), whereas the maxima extracted from the SA-ODFs mostly give corrupted or completely false fiber orientation estimation even for simple configuration as single fibers.

The quantitative results, in Tab. 1, confirms our previous remarks. In particular, at SNR=10, the SA-ODFs are not able to provide proper diffusion directions. The DNC mean value is higher than 0.5 meaning that, in average, more than half of the maxima in each voxel are not detected. The results regarding our PDL approach give more accurate diffusion directions at every SNR.

Besides the good directional information given by the ODFs estimated with our PDL approach, we can also estimate the EAP and interpolate/extrapolate the diffusion signal on the entire q-space, whereas QBI only estimates the ODFs. More results regarding the estimation of the latter features are given in Tab. 4, and we will discuss these results in Sec. 4.3.3.

Multiple shells sampling scheme with 15 measurements

In these experiments, we compare a SHORE reconstruction using a sparse prior (ℓ 1 -SHORE) to our dictionary reconstruction. The SHORE basis has been introduced by [START_REF] Ozarslan | Simple harmonic oscillator based reconstruction and estimation for threedimensional q-space mri[END_REF] and was used in [START_REF] Merlet | Diffusion and multiple orientations from 1.5 mr systems with limited gradient tables[END_REF][START_REF] Merlet | Ensemble average propagator reconstruction via compressed sensing: Discrete or continuous bases ?[END_REF][START_REF] Merlet | Parametric dictionary learning for modeling eap and odf in diffusion mri[END_REF] in the context of sparse recovery. The ℓ 1 -SHORE method consists in solving the LASSO problem (see Eq. 12) using coordinate descent while replacing the dictionary D by the SHORE basis [START_REF] Ozarslan | Simple harmonic oscillator based reconstruction and estimation for threedimensional q-space mri[END_REF]. In order to provide a fair comparison a SH order of L = 8 is used for the generation of the SHORE basis. In both methods (ℓ 1 -SHORE and our PDL approach) we use cross validation [START_REF] Ward | Compressed sensing with cross validation[END_REF] to assess the regularization parameter λ in Eq. 12. [START_REF] Ozarslan | Simple harmonic oscillator based reconstruction and estimation for threedimensional q-space mri[END_REF] and [START_REF] Cheng | Theoretical analysis and practical insights on eap estimation via a unified hardi framework[END_REF] respectively provide closed formulae to estimate the ODF and the EAP when the diffusion signal is modeled in the SHORE basis.

We see in Tab. 2 that our PDL approach outperforms the ℓ 1 -SHORE reconstruction in terms of angular error, difference in the number of compartments (DNC), signal NMSE and EAP NMSE.

Fig. 7 shows the ODFs estimated via our dictionary approach (D-ODF) and the ODFs estimated via the ℓ 1 -SHORE method (SHORE-ODF). It qualitatively indicates an improvement of the angular information given by the D-ODFs over the SHORE-ODFs.

Regarding the PDL approach proposed at the ISBI contest, we obtained an angular error equal to 14.5 • at SNR=10, equal to 11 • at SNR=20, equal to 9.5 • at SNR=30. The results based on our new proposed framework, shown in Tab. 2, indicates an improvement on the accuracy of the maxima estimation compared to the results obtained in the contest. Note that the PDL approach proposed in the contest was already the best in its category.

In the following we show the impact of an increase of the number of samples while keeping the two shells at b-values b = 1500, 2500 s • mm -2 .

Multiple shells sampling scheme with 64 measurements

Now, we use the MS sampling scheme with 64 measurements and we compare once again the ℓ 1 -SHORE method to our PDL approach. Our dictionary approach still outperforms the SHORE reconstruction in terms of angular error, DNC, signal NMSE and EAP NMSE. Fig. 8 shows the ODFs estimated via our PDL approach (D-ODF) and the ODFs estimated via the ℓ 1 -SHORE method (SHORE-ODF). We observe that the D-ODFs give a very accurate estimation of the underlying fiber structure where the SHORE-ODFs fail to provide coherent fiber 1: Dictionary based ODF estimation (D-ODF) versus solid angle ODF via QBI (SA-ODF) using a single shell sampling scheme with 64 measurements. We added rician noise from SNR=10 to 30. Two metrics are shown : the angular error and the difference in the number of compartiments (DNC).

Figure 6: Dictionary based ODF estimation (D-ODF,top) versus solid angle ODF via QBI (SA-ODF,bottom). We added rician noise from SNR=10 to 30 (left to right). We also show the maxima extracted from the ODFs. On the right, we show the ground ODFs and their respective maxima. We added rician noise from SNR=10 to 30 (left to right). We also show the maxima extracted from the ODF. On the right, we show the ground ODFs and their respective maxima.

MS sampling scheme

Figure 8: Dictionary based ODF estimation (D-ODF, top) versus ℓ 1 -SHORE based ODF estimation (SHORE-ODF, bottom) using a multiple shells sampling scheme with 64 measurements. We added rician noise from SNR=10 to 30 (left to right). We also show the maxima extracted from the ODFs. On the right, we show the ground ODFs and their respective maxima.

direction estimation, especially at SNR=10. We confirm this by the quantitative results shown in Tab. 3. We also compare the results of our PDL approach using this MS sampling scheme , i.e. 2 shells at b-values b = 1500, 2500 s • mm -2 with 64 measurements, to our PDL approach using the SS sampling scheme studied in Sec. 4.3.1, i.e. one shell at b-value b = 3000 s • mm -2 with 64 measurements. All the results are given in Tab. 4. Regarding the directional results (angular error and DNC), we observe a slight advantage with the SS sampling scheme. Moreover, if we want to reconstruct the full diffusion signal and the EAP , we notice that a MS sampling scheme is more adequate. This is because the radial information of the diffusion process is better considered when using 2 shells instead of only one shell.

Discussion on experiments with synthetic data

We studied the choice of the free parameters in the learning process, i.e. the dictionary angular and radial orders (4.2), and λ in Eq. 9 (4.1). Our study led to a very good estimation of the diffusion direction (via the computation of the ODF), the diffusion signal and the EAP.

In particular, our PDL approach was shown to better estimate the diffusion directions than the solid angle ODF via QBI does and to compute the diffusion signal and the EAP in a more accurate way than a SHORE reconstruction (using a sparse prior) does.

Regarding the sampling scheme, we observe a slight advantage for the SS sampling scheme considering the directional features, but we found that the MS sampling at b-values b = 1500, 2500 s • mm -2 is more appropriate, because it enables one to deal with the radial information in addition to the angular information. In Sec. 4.3.3 the dictionary reconstruction combined with a multiple shells sampling scheme was proved to efficiently approximate the diffusion signal, the EAP and the ODF.

Regarding the minimal number of measurements required before a large decrease of the reconstruction accuracy, we also perform a last experiment on synthetic data, in which we study the impact of the number of samples N on the error metrics presented in this experiment part with synthetic data. For this purpose, we still consider the same set of signals (different from the training data set) and add Rician noise with S NR = 20. We use the MS sampling scheme with 2 shells at b-values b = 1500, 2500 s • mm -2 and vary the number of samples between N = 5 and N = 100. Fig. 9 shows the resulting values. In this figure, we also plot a vertical line, which represents the number of samples N where the metric errors show a large increase of their values (which means that the estimations are not correct anymore). Overall, we define this bound to N min ≈ 15. In Fig. 9, we also represent the variance of each metric error. We observe from the variance, that the estimation of the diffusion signal and the EAP are robust to noise. However, the estimation of the diffusion direction is more sensitive to noise when we reduce the number of samples.

In conclusion, these synthetic experiments show that our PDL approach can sparsely model multi-fiber compartments

DNC

Signal NMSE EAP NMSE MS-SC SS-SC MS-SC SS-SC MS-SC SS-SC MS-SC SS-SC SNR 30 5.6233 4.9398 0.2187 0.2068 0.0035 0.0207 0.0066 0.0298 SNR 20 6.3080 5.6386 0.2309 0.2102 0.0054 0.0239 0.0104 0.0321 SNR 10 8.3224 8.2530 0.2511 0.2500 0.0144 0.0433 0.0222 0.0526 Table 4: Comparison between the multiple shells sampling scheme (MS-SC) with 64 measurements and the single shell sampling scheme (SS-SC) with 64 measurements, on reconstruction using the PDL approach. Four metrics are shown: the angular error, the difference in the number of compartments (DNC), the signal NMSE, the EAP NMSE. The simulations are performed with Rician noise at SNR=10, 20 and 30.

signals with the assumption of mono-exponential signal decay with b-value. In addition, it overcomes the preliminary PDL approach presented at the HARDI contest at ISBI 2012, already ranked first in its category.

Experiments on real data

In this section, we propose to validate our parametric dictionary learning (PDL) method on real data from human brains. For this purpose, we acquired three distinct sets of data :

• A first set of measurements coming from a 7T scanner, used both to learn the dictionary and to validate it.

• A second set of measurements coming from a 3T scanner, used for the learning process.

• A third set of measurements coming from a 3T scanner, used to validate the dictionary learned on the previous 3T scanner data.

Learning and reconstruction on a 7T scanner data

Training data was acquired on a 7T whole-body MR scanner (MAGNETOM, Siemens Healthcare) equipped with Siemens-AC072 whole body gradient coils, and an 24 channel phased array coil (Nova Medical).

12 axial slices were acquired with a 2D single shot DW-STEAM-EPI (TR/TE/∆/δ = 3000/58/120/15 ms) sequence giving 2 mm isotropic resolution. The echo time (TE), diffusion time (∆), and gradient duration (δ) were optimized to provide maximum SNR for a maximum b-factor. 8 different b-values b=500/1000/2000/3000/4000/5000/6000/7000 s/mm 2 are considered with 70 orientations at each b-value, and an imaging matrix of 96x96x12. The measurement locations are distributed using the algorithm given in [START_REF] Caruyer | Optimal Design of Multiple Q-shells experiments for Diffusion MRI[END_REF] by setting the parameters in such way that there are a constant number of measurements per b-value. We use 11 slices as training data set, in which only the voxels corresponding to a Fractional Anisotropy (FA) superior to 0.20 are taken into account. The FA is computed from the diffusion tensor estimated with the whole set of measurements.

We validate the reconstruction, based on the previously learned dictionary, on the twelfth slice. We compare the ODFs estimated via our PDL approach (D-ODF) to the ODFs estimated with the ℓ 1 -SHORE method (SHORE-ODF) on three different sampling schemes. The three sampling schemes consider measurements at b-values b=1000/2000 s/mm 2 , and we change the number of measurements as N=15, 30 and 60. To distribute the samples between the two b-values, we use the algorithm given in [START_REF] Caruyer | Optimal Design of Multiple Q-shells experiments for Diffusion MRI[END_REF] by setting the parameters in such a way that the the number of samples on each shell is proportional to q 1 [START_REF] Caruyer | Optimal Design of Multiple Q-shells experiments for Diffusion MRI[END_REF][START_REF] Merlet | Impact of radial and angular sampling on multiple shells acquisition in diffusion mri[END_REF].

We choose a region of interest and show the estimated ODFs (see Fig. 10) along with the extracted maxima. This region contains several crossing configurations with different degree of crossing, and thus is appropriate for ODF validation. In Fig. 10, the middle corresponds to the D-ODFs and the bottom to the SHORE-ODFs. From the left to the right, we see the results for N=15, 30 and 60 samples. We consider the ground truth as the estimated signal using the ℓ 1 -SHORE method when all the measurements are taken into account, i.e. 70x8 measurements. The corresponding ground truth ODF are shown at the top of Fig. 10.

Overall, the SHORE-ODFs lead to more false maxima than the D-ODFs. For instance, at N = 60 in the region A (in red), the SHORE-ODFs are not able to correctly resolve the crossing fiber configuration (erroneous number of detected maxima). This phenomenon is emphasized when the number of measurements decreases, whereas our dictionary estimation still provides a coherent map of ODFs.

Learning and reconstruction on a 3T scanner data

We also train our dictionary on data from a 3T Verio (MAGNETOM, Siemens Healthcare) scanner equipped with a 32-channel head coil. The data were acquired at a spatial resolution of 2 mm 3 isotropic, for 6 different b-values b=500/1000/2000/3000/5000/7000 s/mm 2 , 70 orientations at each b-value, and an imaging matrix of 128x128x60. The dictionary is learned from all the measurements on the axial slices 25 to 35, for the voxels with FA ≥ 0.20. The FA is computed from the diffusion tensor estimated with the whole set of measurements.

For the reconstruction, we use data from a 3T scanner (Philips Achieva) equipped with a 8-channel SENSE coil. The data were acquired at a spatial resolution of 2 mm 3 isotropic, for 6 different b-values b=500/1000/2000/4000/6000/8000 s/mm 2 , 70 orientations at each b-value, and an imaging matrix of 128x128x60. We consider three sampling schemes with N=15, 30 and 60 samples and distribute them proportionally to q 1 on two b-values b=1000/2000 s/mm 2 [START_REF] Caruyer | Optimal Design of Multiple Q-shells experiments for Diffusion MRI[END_REF][START_REF] Merlet | Impact of radial and angular sampling on multiple shells acquisition in diffusion mri[END_REF].

In Fig. 11, the top corresponds to the D-ODFs and the bottom to the SHORE-ODFs, on a selected region of interest. From the left to the right, we see the results for N=15, 30 and 60 samples. We also consider the ground truth as the estimated signal using the ℓ 1 -SHORE method when all the measurements are taken into account, i.e. 70x6 measurements. The corresponding ground truth ODFs are shown at the top of Fig. 11.

We obtain sharper ODFs estimation with our PDL approach on the selected region of interest. Again, we see that our dictionary is still able to model the fiber crossing configuration with only 15 samples, and provide less noisy ODFs than the ODFs estimated with SHORE.

Conclusion

We have proposed an original and efficient computational framework to model continuous diffusion MRI (dMRI) signal and to recover its important features such as the ODF and the EAP with a reduced number of measurements. The idea, we implemented, has been to use a parametric dictionary learning algorithm and to exploit the sparse property of a well designed dictionary to recover the diffusion signal and its features. Numerous experimental results have been carried out for validation on synthetic and human brain data acquired from 7-T and 3-T scanners. We have shown that we can clearly recover the diffusion signals and its features with a much better accuracy compared to state-of-the-art approaches and can accurately recover ODF in regions of multiple fiber crossing, even with a small and reduced number of measurements. This opens new perspectives for some dMRI applications, including for example tractography, where the improved characterization of the fiber orientations is likely to greatly and quickly help tracking through regions with and/or without crossing fibers

Appendix A. Derivation of the normalization constant

We define a dictionary of M functions f m , from R 3 to R. We want to normalize the atoms of the dictionary, for the classical ℓ 2 norm:

|| f m || 2 2 = R 3 f m (q) 2 dq = 1. (A.1)
The functions f m are constructed from elementary functions g (m) i, j (in what follows, we drop the index m for the sake of clarity):

f m = I i=1 J j=1 g i, j (q), (A.2)
where g i, j (q • u) = γ i, j e -ν i q 2 q ℓ( j) Y j (u).

(A.

3)

The normalization in Eq. A.1 rewrites as

1 = R 3         i, j g i, j (q)         2 dq = I i=1 J j=1 I i ′ =1 J j ′ =1
γ i, j γ i ′ , j ′ ∞ 0 e -ν i q 2 e -ν i ′ q 2 q ℓ( j) q ℓ( j ′ ) q 2 dq S 2

Y j (u)Y j ′ (u)d 2 u (A.4)

Provided that the spherical harmonics functions form an orthonormal basis, we have γ i, j γ i ′ , j ∞ 0 e -(ν i +ν i ′ )q 2 q 2ℓ( j)+2 dq. (A.6)

Let's use the substitution x = (u i +u i ′ )q 2 in the above integral, we have

1 = I i=1 J j=1 I i ′ =1 γ i, j γ i ′ , j 2(ν i + ν i ′ ) ℓ( j)+3/2
∞ 0 e -x x ℓ( j)+1/2 dx (A.7) where J ℓ+1/2 is the standard bessel function of the first kind and order ℓ + 1/2. Then we get

= I i=1 J j=1 I i ′ =1 γ i, j γ i ′ , j 2(ν i + ν i ′ ) ℓ( j)+3/2 Γ ℓ( j) +
P(Rr) = ∞ q=0 u∈S 2         q 2 K k=0 c k 1 √ χ k I i=0 J j=0
γ ki j exp -ν ki q 2 q ℓ( j) Y j (u)

                 2π √ qR ∞ j ′ =0
(±i) ℓ( j ′ ) j l (2πqR)Y j ′ (u)Y j ′ (r)

         dudq, (B.3)
We shorten this expression using the orthonormal property of the spherical harmonic basis, i.e. u∈S 2 Y j (u)Y j ′ (u)du = δ j j ′ . Then formula (B. exp -ν ki q 2 q ℓ( j)+3/2 J ℓ( j)+1/2 (2πqR)dq

I ki j (R)
, (B.4)

We use the formula from [START_REF] Ryzhik | Table of integrals, series and products[END_REF], i.e.

∞ q=0 x ν+1 exp(-αx 2 )J ν (βx) = β ν (2α) ν+1 exp( -β 2 4α )

I ki j (R) = ∞ q=0
exp -ν ki q 2 q ℓ( j)+3/2 J ℓ( j)+1/2 (2πqR)dq (B.5)

= (2πR) ℓ+1/2 (2ν ki ) ℓ+3/2 exp -(2πR) 2 4ν ki (B.6)
Finally, we get a closed form for the propagator :

P(R • r) = K k=0 c k √ χ k I i=0 J j=0
γ ki j (-1) ℓ( j)/2 π ν ki

ℓ( j)+3/2 R ℓ( j) exp -(πR) 2 ν ki Y j (r) (B.7)

Appendix C. Derivation of the Orientation Distribution Function

The ODF is given by Υ(r) = We use the formula from [START_REF] Ryzhik | Table of integrals, series and products[END_REF], i.e.

  al., Tristan-Vega et al., Wedeen et al.) express the ODF Υ(r) as the integration of the EAP over a solid angle, i.e. Υ(r) = ∞ 0 P(R.r)R 2 dR.

Figure 1 :

 1 Figure 1: NMSE between the ground truth ODFs and the estimated ODFs in terms of SH. The arrow indicates the SH order corresponding to a NMSE considered as sufficiently close to zero.

Figure 2 :

 2 Figure2: Examples of radial attenuation in arbitray chosen directions simulated using a multi-Gaussian model mimicking three crossing fibers (plain line). We fit each of these radial attenuation with the radial part of the model in Eq. 4 and radial orders I = 1, 2, 3, using the Levenberg-Marquardt algorithm (The angular order J is set to zero).

Figure 3 :

 3 Figure 3: Training and validation error (ǫ train and ǫ val ) computed for λ in the range [1 • 10 -5 , 1 • 10 -4 ]. The blue curve represents the training error and the other curves represent the validation error in a noiseless case (purple curve), after adding Rician noise to the validation data set with SNR=30 (green curve), SNR=20 (red curve), SNR=10 (black curve). The arrows indicate the minimum of each curve corresponding to the validation error.

Figure 4 :

 4 Figure 4: First 200 ODFs of the dictionary atoms. The atoms are sorted in decreasing energy order from left to right and top to bottom.

Figure 5 :

 5 Figure 5: Sampling schemes used for validation. (a) A single shell sampling scheme with 64 measurements uniformly spread on one shell at a b-value b = 3000 s • mm -2 . (b) A multiple shells sampling scheme with 15 measurements spread on 2 shells at b-values b = 1500, 2500 s • mm -2 . (c) A multiple shells sampling scheme with 64 measurements spread on 2 shells at b-values b = 1500, 2500 s • mm -2 .

Figure 7 :

 7 Figure 7: Dictionary based ODF estimation (D-ODF, top) versus ℓ 1 -SHORE based ODF estimation (SHORE-ODF, bottom).We added rician noise from SNR=10 to 30 (left to right). We also show the maxima extracted from the ODF. On the right, we show the ground ODFs and their respective maxima.

Figure 9 :

 9 Figure 9: Evolution of the angular error, the DNC, the signal NMSE and the EAP NMSE in function of the number of samples N. We added Rician noise with SNR=20.

Figure 11 :

 11 Figure 11: ODFs estimated from a 3T scanner data via our PDL approach (D-ODF, middle) and via the SHORE technique (SHORE-ODF, bottom). N=15, 30 and 60 samples are considered (left to right). We show the ground truth ODF at the top.

S 2 Y

 2 j (u)Y j ′ (u)d 2 u = δ j, j ′ .

u∈S 2 E

 2 (qu) exp(+2πiqRu.r)duq 2 dq. (B.1) We use the spherical plane wave expansion exp(±2πiqRu.r) = 2π √ qR ∞ j=0 (±i) ℓ( j) J ℓ+1/2 (2πqR)Y j (u)Y j (r), (B.2)

γ

  ki j (-1) ℓ( j)/2 Y j (r) ∞ q=0

∞ 0 P

 0 (R • r)R 2 dR. (C.1)We insert[START_REF] Basser | Mr diffusion tensor spectroscopy and imaging[END_REF] in (C.1) and gather all the R-dependant terms within the integral to get

γ

  ki , ℓ( j)) = (ℓ( j) + 1)ki j (-1) ℓ( j)/2 π ν ki ℓ( j)+1 (ℓ( j) + 1)

Table

	Single shell sampling	Angular error	DNC
	scheme with 64 measurements D-ODF SA-ODF D-ODF SA-ODF
	SNR 30	4.9398	6.6509	0.2068	0.2329
	SNR 20	5.6386	7.1436	0.2102	0.2524
	SNR 10	8.2530	12.419	0.2500	0.5993

15 measurements Dictionary

  ℓ 1 -SHORE Dictionary ℓ 1 -SHORE Dictionary ℓ 1 -SHORE Dictionary ℓ 1 -SHORE

		Angular error	DNC		Signal NMSE	EAP NMSE
	with SNR 30	8.6066	14.670	0.2472	0.4010	0.0104	0.0433	0.0176	0.1040
	SNR 20	9.7626	16.313	0.2540	0.4463	0.0170	0.0578	0.0251	0.1122
	SNR 10	13.344	22.354	0.2734	0.4836	0.0343	0.1027	0.0422	0.1350

Table 2 :

 2 Dictionary based reconstruction versus ℓ 1 -SHORE based reconstruction using a multiple shells sampling scheme with 15 measurements. We added rician noise from SNR=10 to 30. Four metrics are shown : the angular error, the difference in the number of compartiments (DNC), the signal NMSE and the EAP NMSE.

	MS sampling scheme	Angular error	DNC		Signal NMSE	EAP NMSE
	with 64 measurements Dictionary ℓ 1 -SHORE Dictionary ℓ 1 -SHORE Dictionary ℓ 1 -SHORE Dictionary ℓ 1 -SHORE
	SNR 30	5.6233	8.8950	0.2187	0.3106	0.0035	0.0368	0.0066	0.0732
	SNR 20	6.3080	9.6641	0.2309	0.3401	0.0054	0.0386	0.0104	0.0746
	SNR 10	8.3224	13.126	0.2511	0.3995	0.0144	0.0752	0.0222	0.0825

Table 3 :

 3 Dictionary based reconstruction versus ℓ 1 -SHORE based reconstruction. We added rician noise from SNR=10 to 30. Four metrics are shown : the angular error, the difference in the number of compartiments (DNC), the signal NMSE and the EAP NMSE.

http://hardi.epfl.ch/

http://hardi.epfl.ch/

http://hardi.epfl.ch/

Acknowledgments

We would like to express our thanks and gratitude to A. Anwander and D. Bibeck for the 7T whole-body MR scanner data, E. Bannier who is part of the Neurinfo Imaging Platform (IN-RIA, INSERM, FEDER, Région Bretagne et Rennes metropole funding) for the 3T Magnetom Verio scanner data and G. Gilbert for the 3T scanner (Philips Achieva) data. This work has been partially supported by the NucleiPark research project (ANR Program Maladies Neurologique et maladies Psychiatriques) and the France Parkinson Association. Part of this work has been done within the framework of the France-Quebec Cooperation and thanks to the support of "Service de Cooperation et dAction Culturelle du Consulat General de France a Quebec "