Haz-Edine Assemlal

David Tschumperlé

Luc Brun

Kaleem Siddiqi

Recent Advances in Diffusion MRI Modeling: Angular and Radial Reconstruction

Keywords: Diffusion MRI Reconstruction, Local Modeling, Angular Sampling, Radial Sampling, Brain Tissue Features

Recent advances in diffusion magnetic resonance image (dMRI) modeling have led to the development of several state of the art methods for reconstructing the diffusion signal. These methods allow for distinct features to be computed, which in turn reflect properties of fibrous tissue in the brain and in other organs. A practical consideration is that to choose among these approaches requires very specialized knowledge. In order to bridge the gap between theory and practice in dMRI reconstruction and analysis we present a detailed review of the dMRI modeling literature. We place an emphasis on the mathematical and algorithmic underpinnings of the subject, categorizing existing methods according to how they treat the angular and radial sampling of the diffusion signal. We describe the features that can be computed with each method and discuss its advantages and limitations. We also provide a detailed bibliography to guide the reader.

Introduction 1.Context

Diffusion magnetic resonance imaging (dMRI) allows one to examine the microscopic diffusion of water molecules in biological tissue in-vivo. Water molecules are in constant thermal motion, with a locally random component, but this motion is constrained by the presence of surrounding structures including nerves, cells and surrounding tissue. Measurements of this diffusion, therefore, reveal micro-structural properties of the underlying tissue. In practice, this imaging modality requires the collection of successive images with magnetic field gradients applied in different directions. A reconstruction step is then used to estimate the 3D diffusion probability density function (PDF) from the acquired images.

Since the development of the first dMRI acquisition sequence in the mid nineteen-sixties, many applications of this modality have emerged. These can be classified into two main categories. The first category aims at diagnosing certain brain abnormalities that alter the dynamics of the diffusion of water in the brain, e.g., as in the case of a stroke. Such changes can be detected very easily in diffusion images, but remain invisible in other "static" imaging modalities, e.g., anatomical MRI and CT. The dMRI has been increasingly exploited by neurologists for the diagnosis of a wide variety of brain pathologies including: tumors (cerebral lymphoma, epidermoid and cholesteatoma cysts), infections (pyogenic brain abscess, encephalitis herpes); degenerative diseases (Creutzfeld-Jakob Disease); inflammatory conditions (multiple sclerosis), and trauma (shock, fracture) [START_REF] Moritani | Diffusion-Weighted MR Imaging of the Brain[END_REF]Ron and Robbins, 2003).

Email address: assemlal@cim.mcgill.ca (Haz-Edine Assemlal) The second class of applications of dMRI focus on the study of the neuroanatomy of the human brain and more specifically on the understanding of its microstructure. In a diffusion image, each voxel has a signal which results from the motion of a large number of water molecules, revealing features of a portion of tissue at an atomic scale. Early in the dMRI literature it became apparent that this summary of local diffusion depends on attributes of the vector magnetic field gradient, and how it affects the profile of the diffusion signal in both the angular and radial directions. However, the hardware requirements being demanding for MRI scanners of that time, it was not until the early nineties that this kind of imaging could be initiated.

The founding method in this area is Diffusion Tensor Imaging (DTI), where a second-order tensor D is used to model the PDF within a voxel, a model which is adequate when there is a single coherent fibres population present. This tensor may be visualized as an ellipsoid with its 3 axes given by the eigenvectors of D, scaled by their corresponding eigenvalues. The eigenvector with the largest eigenvalue reflects the orientation of the fibre population. The extraction of further information from DTI can then take several forms, such as scalar factional anisotropy indices which express the degree to which diffusion is restricted to particular directions at a voxel.

Advances in MRI now enable one to acquire dMRI images with a larger number of magnetic field gradient diffusion encoding directions, typically 64-100, compared to the 6 directions that are the minimum necessary for DTI. With such angular sampling, the tensor formulation can be replaced with mathematical models of higher dimension, leading to High Angular Resolution Diffusion Imaging (HARDI) reconstructions. These models allow for better detection and representation of complex sub-voxel fibre geometries. From HARDI it is in fact possible to obtain a more accurate fibre orientation distribution function (ODF) within a voxel without a prior assumption on the number of fibre populations present, via a spherical deconvolution method [START_REF] Tournier | Direct estimation of the fiber orientation density function from diffusion-weighted mri data using spherical deconvolution[END_REF]. There is also a growing interest in the development of tractography algorithms to group such local estimates to reconstruct complete fibre tract systems [START_REF] Mori | Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging[END_REF]. Typical examples of such algorithms are described in the recent proceedings of the 2009 MICCAI Diffusion Modeling and Fibre Cup Workshop. With such algorithms in hand, it becomes possible to not only recover connectivity patterns between distinct anatomical regions in the human brain, but also to define continuous measures of the degree of connectivity between two such regions.

Whereas the reconstruction of a more precise angular diffusion profile has been the subject of much work in the dMRI community, the study of the radial profile of the diffusion signal has also been carried out using nuclear magnetic resonance (NMR) spectroscopy in place of existing MRI scanners. Spectrometers can be used to acquire radial diffusion profiles with great accuracy, but at the expense of being limited to a small, spatially unresolved, sample. The use of this method reveals a phenomenon of "diffusion-diffraction" caused by the interference between water molecules and the walls of the microscopic structures within biological tissue. From the observation of this phenomenon, it is theoretically possible to extract features in vivo of brain structure at a microscopic scale, including average cell size, axon diameter, probability of diffusion permeability of the walls, etc.

This article presents an in depth review of the diffusion MRI reconstruction and modeling literature, including low angular resolution methods, high angular resolution methods, methods to sample the radial component of the diffusion signal and methods to combine both angular and radial sampling. Our intent is to provide a comprehensive treatment while emphasizing the mathematical and algorithmic underpinnings of each method. We begin by introducing the necessary mathematical background. We first define the notation used throughout this paper in Section 1.2. We briefly present the scanner sequence used in dMRI to acquire the signal in Sections 1.3 and 1.4, which is widely described using two formalisms: the Stejskal-Tanner equation (Section 1.5) and the q-space formalism (Section 1.6). We then review the typical shape of the diffusion signal in the literature in Section 1.7.

Notation

Let Ω x ⊂ R 3 be the position space (x-space), a vector space of dimension 3, whose natural orthonormal basis is the family of vectors {e x = (1, 0, 0), e y = (0, 1, 0), e z = (0, 0, 1)}. Let x ∈ Ω x be a vector so that the coordinates are expressed as x = (x, y, z) T . Similarly, we define the diffusion space Ω q ⊂ R 3 (q-space) so that the natural orthonormal basis is now {u x , u y , u z } and a vector q ∈ Ω q . We model the diffusion MR image as a continuous function E:

E : Ω x × Ω q → R (x, q) → E(x, q) (1)
The figure 1 illustrate the equation (1), so that each couple (x, q), respectively describing the spatial and the diffusion coordinates, Right: the diffusion q-space Ω q for one voxel of the x-space (depicted here as a yellow rectangle).

is associated with a diffusion value E(x, q). For the sake of simplify, we use the equivalent notation E(q) = E(x, q) to refer to the diffusion signal inside any voxel x ∈ Ω x of the image. In the next sections, we detail the acquisition of the position space Ω x and the diffusion space Ω q .

Imaging the k-space

The MRI scanner acquires the three dimensional anatomical image I on a slice-by-slice basis. Let Ω k ⊂ R 2 be the kspace [START_REF] Ljunggren | A simple graphical representation of fourier-based imaging methods[END_REF], a vector space of dimension 2, so that the acquired image I is modeled as the function:

I : Ω k × R × Ω q → R (k, z, q) → I(k, z, q) (2)
where k = (k x , k y ) ∈ Ω k is a vector of the k-space proportional to the areas of the imaging gradients and z is the third coordinate of the vector x ∈ Ω x . The gradients are applied along each slice z of the image I; this creates a gradient in the spin phases within this slice and enables one to locate each spin of the image I at coordinates (k, z). We are interested in the ensemble magnetization vectors image S which characterizes the brain tissue at each voxel x. The image S is given by the two-dimensional Fourier transform F 2D of the image I:

S (x, q) = k∈Ω k I(k, z, q) exp -i2πk • x xy dk (3)
where x xy = (x, y) is a two-dimensional vector so that x = (x xy , z). Therefore, if the image is acquired in the k-space sampling with sufficient samples of k in each slice z, we can recover the image S as illustrated in figure 2. For more details on the encoding of the image by MRI magnetic field gradients, readers can refer to [START_REF] Pipe | Diffusion MRI: from quantitative measurement to in vivo neuroanatomy[END_REF]. In addition to this anatomical MRI sequence presented in this section, the imaging of the diffusion requires an additional specific acquisition sequence.

Diffusion Gradient Sequence

The classical acquisition sequence used in dMRI is the gradient pulsed spin echo (PGSE), introduced by Stejskal and Tanner ( 1965). This sequence can be explained as follows (as illustrated in Figure 3):

• the first diffusion gradient pulse g(0) "labels" the spins of water molecules according to their initial positions at t = 0;

• after the 180 • radio-frequency signal (which causes an adjustment in spin phases at echo time T E), a second diffusion gradient pulse g(∆) is applied which labels the spins after a diffusion time t = ∆;

• the scanner coils receive the diffusion signal at echo time t = T E and two situations are possible at this point: either the water molecules did not move, so that the spin labels cancel out each other; or some molecules moved during the diffusion time lapse, which leads to a signal loss proportional to the displacement of water molecules.

Figure 3: Pulsed Gradient Spin Echo (PGSE) sequence proposed by [START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF]. The echo formation, resulting in the diffusion signal, is acquired at the echo time T E after the 90 • and 180 • radio-frequency pulses.

The diffusion gradient g sequence is expressed as a function of time [START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF]:

g(t) = H(t 1 ) -H(t 1 -δ) + H(t 2 ) -H(t 2 -δ) u, (4)
where t 2 = t 1 +∆, the symbol H being the Heaviside step function and u ∈ S 2 representing the gradient direction.

Unbounded Diffusion

The diffusion signal received by the scanner coils was first described by [START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF]. Indeed, the incorpo-ration of g from the equation ( 4) into the Bloch-Torrey equations [START_REF] Torrey | Bloch equations with diffusion terms[END_REF] leads to the signal attenuation equation as given by [START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF]:

E(b) = S (b) S (0) = exp -γ 2 δ 2 ∆ - δ 3 g T Dg , (5)
(6)

This equation relates the normalized signal decay E at echo time t = T E with the duration, time separation and strength of the magnetic field pulse gradients (δ, ∆ and g respectively), γ the gyromagnetic ratio, and the apparent diffusion coefficient (ADC) D. For simplicity of notation it is convenient to introduce the b factor which groups together the main parameters of the diffusion sequence (Le Bihan, 1991):

b = γ 2 δ 2 (∆ -δ/3)||g|| 2 (7)
Therefore the Stejskal-Tanner equation ( 5) is commonly written as:

E(q) = exp (-bD) . ( 8
)
where the signal notation was simplified to drop its dependance to the diffusion time ∆, i.e. E(q) = E(q, ∆). The result of the equation ( 8) is valid for diffusion in an unrestricted medium and is not restricted to the limit that ∆ δ as it will be the case in the next section [START_REF] Tanner | Restricted self-diffusion of protons in colloidal systems by the pulsed-gradient, spin-echo method[END_REF]).

Restricted Diffusion

The Stejskal-Tanner equation (8) links the observed diffusion signal to the underlying diffusion coefficient, under the assumption that the diffusion is purely Gaussian. However, this Gaussian hypothesis is often violated when the diffusion is hindered, as is the case in the brain due to the presence of white matter fibers [START_REF] Assaf | In vivo and in vitro bi-exponential diffusion of nacetyl aspartate (naa) in rat brain: a potential structural probe?[END_REF][START_REF] Niendorf | Biexponential diffusion attenuation in various states of brain tissue: implications for diffusion-weighted imaging[END_REF].

The propagator formalism enable one to characterize the diffusion without a prior Gaussian assumption [START_REF] Kärger | The propagator representation of molecular transport in microporous crystallites[END_REF]. Within the narrow pulse approximation (NPA) δ → 0, the diffusion gradient does not vary with time anymore g(t)dt = g(0)-g(∆). This greatly simplifies the relationship between the spin phase difference and the position of the molecules during the gradient pulses [START_REF] Callaghan | Principles of Nuclear Magnetic Resonance Microscopy[END_REF][START_REF] Cory | Measurement of translational displacement probabilities by nmr: an indicator of compartmentation[END_REF][START_REF] Tanner | Restricted self-diffusion of protons in colloidal systems by the pulsed-gradient, spin-echo method[END_REF]):

E(q) = p 0 ∈Ω x ρ(p 0 ) p t ∈Ω x P(p 0 |p ∆ )e i2πq•(p 0 -p ∆ ) dp t dp 0 (9)
where ρ(p 0 ) is the spin density at initial time t = 0, which is assumed to be constant in the voxel and is zero elsewhere. The propagator P(p 0 |p ∆ ) gives the probability that a spin at its initial position p 0 will have moved to position p ∆ after a time interval ∆. We introduce the wave-vector q which is defined as

q = γ 2π t 0 g(t)dt = γ 2π g δ (10)
where q ∈ Ω q . This simplification occurs in the NPA regime, and as a consequence q is not a function of time anymore and only depends on the encoding time δ. An excellent review of the relationship between the Stejskal-Tanner equation and the q-space formalism appears in [START_REF] Basser | Relationships between diffusion tensor and q-space mri[END_REF].

Remark. An alternative convention q = γδg is sometimes found in the diffusion MRI literature, which comes from the convention of the Fourier transform without the exponential term 2π. In this paper, we use the convention described in equation (10).

In this formalism, the b factor defined in the equation ( 7) can also be expressed as a function of the q wave vector:

b = (2π) 2 (∆ -δ/3)q 2 (11)
where q stands for the norm of the wave-vector, i.e. q = ||q||. Let p be the net displacement vector p = p ∆p 0 , then the diffusion Ensemble Average Propagator (EAP) P(p) of a voxel is defined as

P(p) = p 0 ∈Ω x ρ(p 0 )P(p 0 |p 0 + p)dp 0 (12)
Therefore, the diffusion signal E(q) at the diffusion time ∆ is linked to the EAP by the following relationship [START_REF] Callaghan | Principles of Nuclear Magnetic Resonance Microscopy[END_REF][START_REF] Cory | Measurement of translational displacement probabilities by nmr: an indicator of compartmentation[END_REF][START_REF] Stejskal | Use of spin echoes in a pulsed magnetic-field gradient to study anisotropic, restricted diffusion and flow[END_REF]:

E(q) = p∈Ω x P(p) exp i2πq T p dp = F -1 3D [P](q) (13)
Equation ( 13) is the inverse three dimensional Fourier transform F -1 3D of the average propagator, with respect to the displacement vector p ∈ Ω p of water molecules inside a voxel at position p. Similarly, the three dimensional Fourier transform F 3D of E relates to P:

P(p) = q∈Ω x E(q) exp -i2πq T p dq = F 3D [E](p) (14)
In other words, the diffusion signal E is acquired in the q-space which is the Fourier domain. The q-space formalism (or gspace [START_REF] Stejskal | Use of spin echoes in a pulsed magnetic-field gradient to study anisotropic, restricted diffusion and flow[END_REF]) was introduced by [START_REF] Callaghan | Principles of Nuclear Magnetic Resonance Microscopy[END_REF][START_REF] Cory | Measurement of translational displacement probabilities by nmr: an indicator of compartmentation[END_REF][START_REF] Stejskal | Use of spin echoes in a pulsed magnetic-field gradient to study anisotropic, restricted diffusion and flow[END_REF]. The q-space formalism is able to describe in a common framework numerous estimation methods. In this case, the diffusion signal is related to the displacement of the molecules between t = 0 and t = ∆ (c.f . figure 4a). However, this formalism is based on several assumptions about the diffusion signal acquisition [START_REF] Callaghan | Principles of Nuclear Magnetic Resonance Microscopy[END_REF][START_REF] Cory | Measurement of translational displacement probabilities by nmr: an indicator of compartmentation[END_REF]):

1. It is assumed that the gradient time δ ≈ 0 is negligible so that there is no displacement of water molecules during that period. In practice, the gradient duration δ is not always negligible and the measured displacement is related to the mean position of the molecules between time intervals t = [0, δ] and t = [∆, δ] (c.f . figure 4b); 2. It is assumed that that the gradient time δ is very short compared to the diffusion time ∆, so that the q wave-vector can be assumed to not be a function of time (c.f . equation 10). (a) Ideal case (δ ≈ 0), (b) Real case (δ 0). Adapted from [START_REF] Hagmann | From diffusion mri to brain connectomics[END_REF].

In the present paper we assume that these assumptions are valid, since they are necessary for the definition of q-space. Several studies [START_REF] Bar-Shir | The effect of the diffusion time and pulse gradient duration ratio on the diffraction pattern and the structural information estimated from q-space diffusion MR: experiments and simulations[END_REF][START_REF] Blees | The effect of finite duration of gradient pulses on the pulsedfield-gradient NMR method for studying restricted diffusion[END_REF][START_REF] Coy | Pulsed gradient spin echo nuclear magnetic resonance for molecules diffusing between partially reflecting rectangular barriers[END_REF][START_REF] Mair | The narrow pulse approximation and long length scale determination in xenon gas diffusion NMR studies of model porous media[END_REF][START_REF] Mitra | Effects of finite gradient-pulse widths in pulsedfield-gradient diffusion measurements[END_REF] have shown that when these assumptions do not hold, the Fourier relationship (equation ( 13) and equation ( 14)) is still valid, but that the interpretation of the signal should be somewhat different. Instead of the diffusion propagator formalism, a center of mass propagator formalism applies, due to a homothetic transformation of the signal features.

Remark. The fundamental explicit relationship between the acquired image I and the diffusion propagator P is expressed by combining equations (3) and (14):

P(x, p) = F 3D E(x, q) (x, p) = F 3D F 2D [I](k, z, q) F 2D [I](k, z, 0) (x, p) (15)
where E(x, q) = S (x, q)/S (x, 0) which normalizes the signal so that variations of E along q can be attributed solely to diffusion.

Modeling the Shape of the Diffusion Signal

Figure 5 provides a qualitative sense of the expected diffusion profile in several typical situations in brain white matter. When no fibers are present the diffusion is typically equal in all directions, with a Gaussian distribution radially (Figure 5a). When a single fibers bundle is present there is maximum diffusion along its direction (Figure 5b). When two fibers bundles cross there is preferential diffusion in the direction of each (Figure 5c). Finally, Figure 5d illustrates a case where there is equal diffusion all directions, but that there is a variation in speed radially which indicates the apparent radius of the fibers bundle.

It is presently impractical to acquire a dense radial and angular sampling of the diffusion space due to the significant acquisition time this would imply. As a result, several advanced methods exist in the dMRI literature to sample and process the signal. The analysis of the diffusion signal is closely related to the sampling of the q-space as illustrated in Figure 6. In the remainder of this article we review the state of the art methods, organized into three groups: low and high angular sampling methods (Sections 2 and 3), radial sampling methods (section 4), and methods which combine radial and angular sampling (Section 5). We also include an Appendix which reviews the mathematical concepts and tools required to understand and apply these methods. The description of each method is organized according to the following themes:

1. Local diffusion modeling: we express the mathematical models used for the interpretation of the diffusion signal in the q-space, while clarifying the underlying assumptions; 2. Model estimation from the data: we describe the model estimation methods from the data acquired from the MRI scanner; 3. Processing and extraction of diffusion features: we present common post-processing techniques on the estimated model, to extract features of tissue microstructure; 4. Advantages and limitations of the method: We enumerate the main advantages (denoted with the + symbol) and limitations (denoted with asymbol) of the method.

Low Angular Resolution Diffusion Imaging (DTI)

Diffusion anisotropy, as captured by diffusion nuclear magnetic resonance (NMR), was pointed out in early investigations involving controlled environments [START_REF] Stejskal | Use of spin echoes in a pulsed magnetic-field gradient to study anisotropic, restricted diffusion and flow[END_REF][START_REF] Tanner | Transient diffusion in a system partitioned by permeable barriers. application to nmr measurements with a pulsed field gradient[END_REF][START_REF] Tanner | Restricted self-diffusion of protons in colloidal systems by the pulsed-gradient, spin-echo method[END_REF]). This was followed by studies of the diffusion process within brain tissue using dMRI [START_REF] Moseley | Diffusion-weighted MR imaging of anisotropic water diffusion in cat central nervous system[END_REF]. The acquired diffusion image depends on the orientation u of the diffusion wave-vector q. Hence, the use of a tensor, a rotationally invariant object, is convenient to characterize the anisotropy of the Apparent Diffusion Coefficient (ADC) of brain tissue, as suggested in several studies [START_REF] Casimir | On onsager's principle of microscopic reversibility[END_REF][START_REF] De Groot | Non equilibrium thermodynamics[END_REF]Onsager, 1931a,b;[START_REF] Stejskal | Use of spin echoes in a pulsed magnetic-field gradient to study anisotropic, restricted diffusion and flow[END_REF]. In this section, we describe the major method used for the characterization of this anisotropy of diffusion, in the case of low angular sampling of the q-space.

Local diffusion modeling Whereas the scalar ADC measure is modeled with a zeroth-order tensor, DTI introduces the use of a second-order tensor D, allowing a more accurate angular characterization of the diffusion process in the brain (Basser and LeBihan, 1992;[START_REF] Filler | Images, apparatus, algorithms and methods[END_REF][START_REF] Moseley | Diffusion-weighted MR imaging of anisotropic water diffusion in cat central nervous system[END_REF][START_REF] Stejskal | Use of spin echoes in a pulsed magnetic-field gradient to study anisotropic, restricted diffusion and flow[END_REF]. The mathematical framework which explicitly relates the diffusion tensor to the NMR signal was demonstrated by (Basser and LeBihan, 1992;[START_REF] Basser | Estimation of the effective selfdiffusion tensor from the nmr spin echo[END_REF][START_REF] Stejskal | Use of spin echoes in a pulsed magnetic-field gradient to study anisotropic, restricted diffusion and flow[END_REF]:

E(q) = exp -4π 2 τq T Dq (16)
In this formalism, the average local diffusion process is described by a second-order tensor D, whose coordinates in the q-space basis {u x , u y , u z } are given by a 3 × 3 symmetric and positivedefinite matrix:

          D xx D xy D xz D xy D yy D yz D xz D yz D zz           (17)
In a environment such as water, the diffusion process D is assumed to be symmetric (i.e. D = D T ) according to the principles of thermodynamics [START_REF] Basser | Estimation of the effective selfdiffusion tensor from the nmr spin echo[END_REF][START_REF] Casimir | On onsager's principle of microscopic reversibility[END_REF][START_REF] De Groot | Non equilibrium thermodynamics[END_REF]Onsager, 1931a,b;[START_REF] Stejskal | Use of spin echoes in a pulsed magnetic-field gradient to study anisotropic, restricted diffusion and flow[END_REF].

Model estimation from the data The apparent diffusion tensor profile D is expressed as a function of the wave-vector diffusion q defined in the q-space, so that the logarithm of the equation ( 16) is:

D(q) = q T Dq = - ln (E(q)) 4π 2 τ (18)
Since the diffusion tensor D is symmetric, it is entirely defined by six components which can be grouped into a vector D [START_REF] Basser | Estimation of the effective selfdiffusion tensor from the nmr spin echo[END_REF]:

D = (D xx , D xy , D xz , D yy , D yz , D zz ) T . ( 19
)
The construction of the sampling matrix of the q-space requires at least n = 6 acquisitions q i , i ∈ [1, n] and one additional acquisition at q = 0 for normalization [START_REF] Basser | Estimation of the effective selfdiffusion tensor from the nmr spin echo[END_REF][START_REF] Stejskal | Use of spin echoes in a pulsed magnetic-field gradient to study anisotropic, restricted diffusion and flow[END_REF][START_REF] Tuch | Diffusion mri of complex tissue structure[END_REF]:

B = 4π 2 τ                  q x 1 q x 1 2q x 1 q y 1 2q x 1 q z 1 q y 1 q y 1 2q y 1 q z 1 q z 1 q z 1 q x 2 q x 2 2q x 2 q y 2 2q x 2 q z 2 q y 2 q y 2 2q y 2 q z 2 q z 2 q z 2 . . . . . . . . . . . . . . . . . . q x n q x n 2q x n q y n 2q x n q z n q y n q y n 2q y n q z n q z n q z n                  (20)
where the wave-vector is decomposed as q i = q x i u x + q y i u y + q z i u z . The sampling matrix defined in the equation ( 20) is traditionally named the B-matrix, in reference to its multiple b-factor entries b i j = 4π 2 τq i q j (c.f . equation ( 7)).

The logarithm of the data samples E i , i ∈ [1, n] are grouped in a vector Y: Y = (-ln(E 1 ), -. . . ,ln(E n )) T .

(21) Finally, equation ( 16) which links the model to the data, is expressed in the matrix form as:

Y = BD. ( 22
)
In the case where there are exactly 6 acquisitions in different orientation of the q-space, the components of the diffusion tensor can be computed by the relationship D = B -1 Y. However, such a process is very sensitive to the quality of the data and to perturbations due to acquisition noise. In practice, scanners are now able to acquire many mores images (typically up to n = 60 directions). From these n images, the tensor is estimated to be the one which minimizes a notion of error to the set of acquired data (Johansen-Berg and [START_REF] Johansen-Berg | Diffusion MRI: from quantitative measurement to in vivo neuroanatomy[END_REF]. There are several methods for the estimation and regularization of second order tensor fields, including: weighted least squares [START_REF] Basser | Estimation of the effective selfdiffusion tensor from the nmr spin echo[END_REF], variational methods for the estimation of the image volume with positivity and regularity constraints [START_REF] Chefd'hotel | Constrained flows of matrix-valued functions: Application to diffusion tensor regularization[END_REF][START_REF] Chefd'hotel | Regularizing flows for constrained matrix-valued images[END_REF][START_REF] Neji | A Convex Semi-definite Positive Framework for DTI Estimation and Regularization[END_REF]Tschumperlé and Deriche, 2003a,b), estimation in a Riemannian space [START_REF] Arsigny | Processing data in Lie groups: An algebraic approach. application to non-linear registration and diffusion tensor MRI[END_REF][START_REF] Fillard | Clinical dt-mri estimation, smoothing, and fiber tracking with log-euclidean metrics[END_REF][START_REF] Lenglet | Geometric and variational methods for diffusion tensor mri processing[END_REF] and the use of sparse representations [START_REF] Bao | Denoising of human cardiac diffusion tensor magnetic resonance images using sparse representation combined with segmentation[END_REF][START_REF] Luo | Denoising by averaging reconstructed images: Application to magnetic resonance images[END_REF].

Processing and extraction of diffusion features DTI data is often visualized using a field of ellipsoids, where each voxel is represented by an isosurface of the diffusion tensor [START_REF] Basser | Inferring microstructural features and the physiological state of tissues from diffusion-weighted images[END_REF]:

D(q) = q T Dq = constant. ( 23
)
In such a visualization, the eigenvectors of D give the principal axes of the ellipsoid, with their lengths scaled by the corresponding eigenvalues. Figure 7 illustrates DTI data of an adult human brain data. In regions where there is a dominant local fiber direction the ellipsoids appear elongated to reflect it. The formalism of the second order tensor contains rich information at a voxel scale. Nonetheless, this information is hard to grasp qualitatively in a volumetric image consisting of thousands of voxels. As a result, a wide variety of measures have been proposed to process and simplify the diffusion information reflected in the second order tensor. In the following paragraphs, we classify these measures into two groups: scalar features and vectorial features. Scalar features: There are various scalars features based on the second order tensor which summarize the diffusion information at each voxel, enabling the depiction of qualitative aspects of diffusion on a larger scale [START_REF] Papadakis | A study of rotationally invariant and symmetric indices of diffusion anisotropy[END_REF][START_REF] Pierpaoli | Toward a quantitative assessment of diffusion anisotropy[END_REF]. In contrast to the DWI and ADC methods, these scalar images are not dependent on the orientation of the patient with respect to the scanner, and they allow for straightforward inter-subject statistics to be computed. An example of such a scalar feature is the trace of the diffusion tensor, which is defined as the sum of the eigenvalues λ i , computed in the basis of tensor eigenvectors (see Appendix 7.2):

Trace = 3 i=1 λ i (24)
It is interesting to note that the trace/3 is the mean diffusion in a voxel. In a healthy patient this feature tends to lead to an image with relatively homogeneous intensity [START_REF] Pierpaoli | Toward a quantitative assessment of diffusion anisotropy[END_REF]. This is useful in clinical applications since it is likely that any diffusion anomaly (such as due to an acute brain stroke) will lead to a local hyper-or hypo-intensity [START_REF] Lythgoe | Effects of diffusion anisotropy on lesion delineation in a rat model of cerebral ischemia[END_REF].

A second popular feature is the fractional anisotropy (FA) which measures the normalized angular variance of diffusion within a voxel. It is defined as:

FA = 3 2 3 i=1 (λ i -λ ) 2 3 i=1 λ 2 i , ( 25
)
where λ is the mean of the eigenvalues λ i . At a voxel scale, the FA measures the alignment and the coherence of the underlying microstructure. Thus, voxels with compact fibers bundles have a high FA value (e.g. in the corpus callosum), whereas voxels with non aligned bundles tend to have a low FA value (e.g., as in the case of crossing fibers). The FA is likely the most used scaler anisotropy index in the dMRI literature because it is relatively robust to noise in comparison to other features [START_REF] Pierpaoli | Toward a quantitative assessment of diffusion anisotropy[END_REF].

Several features related to the geometric shape of the second order tensor have also been proposed in the literature [START_REF] Basser | New histological and physiological stains derived from diffusion-tensor MR images[END_REF][START_REF] Westin | Geometrical diffusion measures for MRI from tensor basis analysis[END_REF]. These are typically based on the ratio between the eigenvalues, sorted in a decreasing order, and can help distinguish between flat, elongated, and other profiles.

Vectorial features: It is often useful to extract vectorial features based on the DTI. A colour-coded orientation map of nerve fibers in the brain [START_REF] Pajevic | Color schemes to represent the orientation of anisotropic tissues from diffusion tensor data: application to white matter fiber tract mapping in the human brain[END_REF] is shown in Figure 8a). The direction of nerve fibers is also shown using random texture smoothing [START_REF] Tschumperlé | Vector-valued image regularization with PDEs: A common framework for different applications[END_REF][START_REF] Weickert | Visualization and Image Processing of Tensor Fields[END_REF]) (Figure 8c) and using three dimensional curves [START_REF] Conturo | Tracking neuronal fiber pathways in the living human brain[END_REF]Jones et al., 1999b;[START_REF] Mori | Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging[END_REF])(Figure 8b) The process of tracking to reconstruct fiber bundles in DTI is generally based on a continuous tensor field which requires an interpolation of tensors [START_REF] Arsigny | Processing data in Lie groups: An algebraic approach. application to non-linear registration and diffusion tensor MRI[END_REF][START_REF] Fillard | Clinical dt-mri estimation, smoothing, and fiber tracking with log-euclidean metrics[END_REF][START_REF] Pennec | A Riemannian framework for tensor computing[END_REF]. (b) The curvature encodes the direction of each fiber and the brightness the smoothness strength. (c) Fibertracking: A visualization of nerve fiber bundles as three dimensional curves overlaid on the diffusion ellipsoids. These images are adapted from [START_REF] Tschumperlé | Pde's based regularization of multivalued images and applications[END_REF].

Advantages and limitations of DTI With improvements in the quality and speed of MRI scanners, the DTI model is somewhat limited. Indeed, it is not uncommon to have an angular acquisition resolution of the q-space of up to 128 directions (in comparison to the minimum 7 directions required for a DTI reconstruction). When such angular precision is available, the model over-simplifies the diffusion process. Although the model is valid when there is a single dominant diffusion direction, such as in the corpus callosum (illustrated in Figure 7), is it not accurate enough to characterize more complex diffusion profiles in the brain. These limits, those of modeling the diffusion coefficient by a second order tensor, can be analyzed through a Taylor expansion of the diffusion signal E at q = 0. The diffusion signal E is usually considered to be a symmetric function, so that the odd terms are equal to zero [START_REF] Wedeen | Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging[END_REF]:

E(q) = exp 0 - (2π) 2 2 q T pp T -p p T q + O ||q|| 4 . ( 26
)
Here • is the expected value and ( pp Tp p T is the covariance matrix of the diffusion propagator. Furthermore, the DTI model discards the terms with order higher than two, under the assumptions that: (i) the diffusion vector q is small (q T pp T q 1 -typically this is true in brain tissues for b < 1000 s/mm 2 ) and (ii) the local diffusion is Gaussian (moments greater than two are zero). Under these assumptions equation ( 26) becomes:

E(q) ≈ exp - (2π) 2 2 q T pp T q = exp -4π 2 τq T Dq (27)
where the coefficient pp T = 2Dτ is given by the Einstein-Smoluchowski relationship [START_REF] Einstein | On the movement of small particles suspended in stationary liquids required by the molecular-kinetic theory of heat[END_REF], yielding equation ( 16) of DTI.

Since the ADC is estimated with the hypothesis that the local diffusion follows a zero mean Gaussian function, the DTI method is expressed in the diffusion propagator P formalism as:

P(p) = 1 (4πτ) 3 |D| exp p T D -1 p 4τ , (28)
where P(p) is the probability that a water molecule inside a voxel travels a distance p during the diffusion time τ. Under these assumptions, the second-order tensor captures the diffusion covariance, which corresponds to the second order moment of the diffusion propagator P.

The Gaussian hypothesis of the DTI model does not hold in the general case, e.g., when crossing fibers are present (c.f . figure 9b).) In such situations the moments of order greater than two are not equal to zero. As a result, the representation of a crossing or any complex sub-voxel geometry is excessively simplified by DTI [START_REF] Basser | In vivo fiber tractography using dt-mri data[END_REF]. Furthermore, fibertracking results based solely on DTI can be unreliable in the presence of crossing fibers or other complex geometries [START_REF] Campbell | Flow-based fiber tracking with diffusion tensor and q-ball data: validation and comparison to principal diffusion direction techniques[END_REF]. We now summarize the advantages and limitations of DTI.

+ The time during which the patient must lie motionless is short, since the reconstruction requires a minimum of only seven acquisitions in q-space.

+ The second order tensor appears to adequately model brain white matter regions where there is a single dominant direction of diffusion. The measured angle between the principal eigenvector and the ground-truth fiber orientation is approximately 13 • [START_REF] Lin | Validation of diffusion tensor magnetic resonance axonal fiber imaging with registered manganese-enhanced optic tracts[END_REF].

+ The DTI model has became a standard method in the dMRI community and there is a vast literature on topics such as noise robustness [START_REF] Basu | Rician noise removal in diffusion tensor mri[END_REF], data regularization [START_REF] Mangin | Distortion correction and robust tensor estimation for MR diffusion imaging[END_REF][START_REF] Tschumperlé | Vector-valued image regularization with PDEs: A common framework for different applications[END_REF], sampling (Jones et al., 1999a;[START_REF] Poupon | Détection des faisceaux de fibres de la substance blanche pour l'étude de la connectivité anatomique cérébrale[END_REF] and statistical studies [START_REF] Chung | Comparison of bootstrap approaches for estimation of uncertainties of DTI parameters[END_REF].

-The modeling of the angular profile of the diffusion signal by a Gaussian function is often violated in the presence of complex geometries, such as fiber crossings, Yconfigurations, bottlenecks, etc. [START_REF] Savadjiev | Labeling of ambiguous subvoxel fibre bundle configurations in high angular resolution diffusion MRI[END_REF].

-Similarly, the modeling of the radial profile of the diffusion signal by a Gaussian function prevents one from retrieving more complete radial information, thus excluding features such as mean cell size, axon diameter, etc. [START_REF] Assaf | In vivo and in vitro bi-exponential diffusion of nacetyl aspartate (naa) in rat brain: a potential structural probe?[END_REF][START_REF] Callaghan | Diffractionlike effects in nmr diffusion studies of fluids in porous solids[END_REF][START_REF] Cory | Measurement of translational displacement probabilities by nmr: an indicator of compartmentation[END_REF][START_REF] Kuchel | Nmr diffusion-diffraction of water revealing alignment of erythrocytes in a magnetic field and their dimensions and membrane transport characteristics[END_REF][START_REF] Regan | NMR studies of diffusion-coherence phenomena in red cell suspensions: current status[END_REF].

High Angular Resolution Diffusion Imaging (HARDI)

High Angular Resolution Diffusion Imaging (HARDI) was proposed by [START_REF] Tuch | High angular resolution diffusion imaging of the human brain[END_REF] to enable a more precise angular characterization of the diffusion signal, while keeping the acquisition time compatible with clinical constraints. HARDI reduces the diffusion signal sampling to a single sphere of the q-space, and has initiated a significant interest in the scientific community as indicated by the numerous methods to tackle the recovery the geometry of crossing fibers, as illustrated in Figure 9. In the following we describe the principal methods in the literature, which fall into two classes: parametric and non-parametric. The first class represents the diffusion signal as a sum of functions, each of which models a single fiber population. The second represents the diffusion signal as a mathematical series. A detailed quantitative comparison of several of these methods is available in (Jian and Vemuri, 2007;[START_REF] Ramirez-Manzanares | A comparison of methods for recovering intra-voxel white matter fiber architecture from clinical diffusion imaging scans[END_REF]. (c) Generalization of the ADC modeling with HARDI sampling. The maxima of the angular profile do not match the underlying fiber bundle directions given in (a). (d) Angular feature estimation of the diffusion propagator with HARDI sampling. Adapted from [START_REF] Descoteaux | High angular resolution diffusion mri: from local estimation to segmentation and tractography[END_REF].

Mixture models (parametric)

Local diffusion modeling A variety of models in the literature assume that the diffusion signal can be decomposed as a weighted sum of generic diffusion models h i :

E(q) = n i=1 f i h i (q) with n i f i = 1, ( 29
)
where f i stands for the weight of the i-th bundle of nerve fibers and n is the total number of bundles. The mixture models of the literature can be expressed in the formalism of the equation ( 29).

The multi-Gaussian (or multi-DTI) model proposed in [START_REF] Tuch | Diffusion mri of complex tissue structure[END_REF] generalizes DTI with the assumption that the diffusion signal is the sum of signals from several bundles of fibers, each modeled as a second-order tensor:

h i (q) = exp(-4π 2 τq T D i q). ( 30
)
The parametric "Ball and Stick" [START_REF] Behrens | Characterization and propagation of uncertainty in diffusion-weighted MR imaging[END_REF][START_REF] Hosey | Inference of multiple fiber orientations in high angular resolution diffusion imaging[END_REF] model describes the signal as a bi-Gaussian function (n = 2), where the two components h iso and h aniso correspond respectively to an isotropic ADC D and an anisotropic ADC D:

h iso (q) = exp -4π 2 q 2 τD (31) h aniso (q) = exp -4π 2 q 2 τu T Du (32) with D = RΛR T = R           1 0 0 0 0 0 0 0 0           R T (33)
with the rotation matrix R being that between the basis of the q-space and the basis defined by the eigenvectors of D.

The CHARMED model [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF] proposes a parametrization with several components, one for the hindered diffusion modeled by a Gaussian function, and the others for the restricted diffusion modeled by diffusion inside a cylinder [START_REF] Neuman | Spin echo of spins diffusing in a bounded medium[END_REF]):

E(q) = f hind h hind (q) + n i=2 f restr i h restr i (q) h hind (q) = exp -4π 2 τ((q ⊥ ) 2 λ ⊥ + (q ) 2 λ ) h restr i (q) = exp -4π 2 K i with K i = τ(q ) 2 D i + (q ⊥ ) 2 7 96 2 - 99 112 R 2 D ⊥ i τ R 4 D ⊥ i τ . ( 34
)
The symbols f hind and f restr are the population fractions of the hindered and restricted terms, respectively; q and q ⊥ are the components of the q vector parallel and perpendicular to the fibers, respectively; λ and λ ⊥ are the eigenvalues of the diffusion tensor parallel and perpendicular to the axons (for a single coherent fiber bundle), respectively; D parallel and D ⊥ are the parallel and perpendicular diffusion coefficients within the cylinder; R is the cylinder radius; and τ is half of the echo time [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF].

Model estimation from the data The estimation of these models is generally non linear and is obtained by an iterative computation in which various priors (e.g. positivity, tensor shape) on n, D i and f i are introduced to increase the numerical stability of the estimation [START_REF] Alexander | Maximum entropy spherical deconvolution for diffusion MRI[END_REF][START_REF] Behrens | Characterization and propagation of uncertainty in diffusion-weighted MR imaging[END_REF][START_REF] Chen | Estimation, smoothing, and characterization of apparent diffusion coefficient profiles from high angular resolution dwi[END_REF][START_REF] Hosey | Inference of multiple fiber orientations in high angular resolution diffusion imaging[END_REF][START_REF] Tuch | High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity[END_REF]. The noise and the number of samples influence the number of parameters and the estimation is carried out by a non-linear iterative numerical process such as gradient descent [START_REF] Tuch | High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity[END_REF], a Levenberg-Marquardt scheme [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF][START_REF] Maier | Biexponential diffusion tensor analysis of human brain diffusion data[END_REF], a Gauss-Newton scheme [START_REF] Peled | Geometrically constrained two-tensor model for crossing tracts in DWI[END_REF] or an unscented Kalman filter [START_REF] Malcolm | Filtered multi-tensor tractography[END_REF]. The choice of the number of components n is either arbitrarily set [START_REF] Tuch | High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity[END_REF], or is chosen by a statistical criterion [START_REF] Behrens | Probabilistic diffusion tractography with multiple fibre orientations: what can we gain?[END_REF][START_REF] Parker | Probabilistic Monte Carlo based mapping of cerebral connections utilising whole-brain crossing fibre information[END_REF][START_REF] Tuch | High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity[END_REF].

Processing and extraction of diffusion features The diffusion features which characterize the brain micro-architecture are directly obtained from the estimated model parameters. Thus, no additional postprocessing is required.

Advantages and limitations

+ The modeling of non-Gaussian angular profiles allows crossing fibers to be identified, which is a problematic case for DTI [START_REF] Alexander | Detection and modeling of nongaussian apparent diffusion coefficient profiles in human brain data[END_REF][START_REF] Frank | Characterization of anisotropy in high angular resolution diffusion-weighted mri[END_REF][START_REF] Özarslan | Generalized scalar measures for diffusion MRI using trace, variance, and entropy[END_REF].

+ The extraction of diffusion features does not require any post-processing once the reconstructions have been obtained.

+ The number of samples in the q-space is modest, e.g., typically 64 directions on a sphere suffice [START_REF] Alexander | Maximum entropy spherical deconvolution for diffusion MRI[END_REF].

-The estimation of the number of compartments requires additional processing [START_REF] Behrens | Probabilistic diffusion tractography with multiple fibre orientations: what can we gain?[END_REF]. Often the choice is somewhat arbitrary and under or over estimation of this parameter can be problematic.

-The estimation of the diffusion signal is based on an empirical hypothesis of the nature of diffusion within a healthy brain, in a context where the origin of the diffusion is still debated [START_REF] Cohen | High b-value q-space analyzed diffusion-weighted mrs and mri in neuronal tissues -a technical review[END_REF][START_REF] Niendorf | Biexponential diffusion attenuation in various states of brain tissue: implications for diffusion-weighted imaging[END_REF]. Consequently, there are no guaranties that these models remain valid for the diffusion in other tissues, e.g., in diseased brains, the heart, etc.

-The stability and the accuracy of the non-linear estimation step depends on the initial parameters, which vary with data [START_REF] Aubert | Mathematical problems in image processing: partial differential equations and the calculus of variations[END_REF]. The estimation time is notably slower than the linear estimation of DTI, since it requires a iterative minimization process, e.g., this takes approximately 2 hours for a slice of 64×64 pixels according to [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF]).

-As is the case with DTI (c.f . Section 2), the modeling of the radial profile by a Gaussian function is inadequate since it excludes a wide range of features, including mean cell size, axon diameters, etc.

Spherical Deconvolution (parametric)

Local diffusion modeling The dMRI estimation by the Spherical Deconvolution (SD) was proposed by [START_REF] Anderson | Sub-voxel of fiber orientation using high angular resolution diffusion tensor imaging[END_REF]Jian and Vemuri, 2007;[START_REF] Tournier | Direct estimation of the fiber orientation density function from diffusion-weighted mri data using spherical deconvolution[END_REF]. Let S 2 be the single sphere domain and SO(3) the rotation group in R 3 . The diffusion signal E is modeled by the convolution of a kernel h ∈ L 2 (S 2 ) and a function f ∈ L 2 (SO(3)) which respectively represent the signal response for a single bundle of nerve fibers and the fiber orientation density function (fODF), ideally composed of n Dirac delta functions for n bundles of fibers. The spherical deconvolution operator is expressed as [START_REF] Healy | Spherical deconvolution[END_REF]:

E(q) = ( f * h(q))(u) = R∈SO(3) f (R)h(q, R T u)dR, (35)
where u ∈ S 2 and the symbol u = R T u stands for the rotation of the vector u by the matrix R. Jian et al. (2007) proposed a kernel based on the Wishart distribution, which is a multidimensional generalization of the χ 2 distribution [START_REF] Wishart | The generalised product moment distribution in samples from a normal multivariate population[END_REF]:

h(q, u ) = 1 + 4π 2 q 2 τu T Du -p (36)
with D ∈ P n and P n the set of the 3 × 3 positive-definite matrices, i.e., u T Du > 0. The parameter p ∈ R should satisfy p ≥ n with n the degree of freedom and this generalizes other approaches in the literature (Jian et al., 2007). Indeed, when p → ∞, the kernel h given by the equation ( 36) is Gaussian as in several methods [START_REF] Anderson | Measurement of fiber orientation distributions using high angular resolution diffusion imaging[END_REF][START_REF] Anderson | Sub-voxel of fiber orientation using high angular resolution diffusion tensor imaging[END_REF]Descoteaux et al., 2009a;[START_REF] Seunarine | Linear persistent angular structure mri and non-linear spherical deconvolution for diffusion mri[END_REF][START_REF] Tournier | Direct estimation of the fiber orientation density function from diffusion-weighted mri data using spherical deconvolution[END_REF]. When p = 2, the kernel h follows a Debye-Porod distribution proposed for the use of dMRI in [START_REF] Sen | Debye-Porod law of diffraction for diffusion in porous media[END_REF].

The Persistent Angular Structure (PAS-MRI) method, proposed by [START_REF] Jansons | Persistent angular structure: new insights from diffusion magnetic resonance imaging data[END_REF], is formally expressed as the function defined on the sphere of radius k ∈ R for which the inverse Fourier transform best fits the signal. This method can also be expressed as a deconvolution of the signal E with a kernel h expressed as [START_REF] Alexander | Maximum entropy spherical deconvolution for diffusion MRI[END_REF][START_REF] Seunarine | Linear persistent angular structure mri and non-linear spherical deconvolution for diffusion mri[END_REF]:

h(q, u ) = k -2 cos kqu T u . ( 37
)
Some other studies propose the use of a kernel h which is dependent on the data, and which is computed by statistical estimation on the whole diffusion image [START_REF] Kaden | Parametric spherical deconvolution: Inferring anatomical connectivity using diffusion MR imaging[END_REF][START_REF] Tournier | Direct estimation of the fiber orientation density function from diffusion-weighted mri data using spherical deconvolution[END_REF].

Model estimation from the data Let f be the orientation density function of fibers. [START_REF] Tournier | Direct estimation of the fiber orientation density function from diffusion-weighted mri data using spherical deconvolution[END_REF] propose to express the spherical convolution expressed in equation ( 35) directly in the spherical harmonics domain, as a simple multiplication of the coefficients h lm and f lm [START_REF] Healy | Spherical deconvolution[END_REF]:

E = ( f * h) = ∞ l=0 ∞ m=-l f lm h lm . (38)
The matrix formulation of equation ( 38) is interesting as the f lm can be computed by a simple matrix inversion. The deconvolution operation is nonetheless unstable in the presence of noisy data, notably when h is an anisotropic kernel. [START_REF] Alexander | Maximum entropy spherical deconvolution for diffusion MRI[END_REF]; [START_REF] Tournier | Resolving crossing fibres using constrained spherical deconvolution: validation using diffusion-weighted imaging phantom data[END_REF] suggest to regularize the function f with a Tikhonov filter. This frequency regularization can favourably be expressed in the spherical harmonic basis as a simple linear relationship which is easy to implement [START_REF] Descoteaux | Apparent diffusion coefficients from high angular resolution diffusion imaging: Estimation and applications[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] Tournier | Direct estimation of the fiber orientation density function from diffusion-weighted mri data using spherical deconvolution[END_REF].

Processing and extraction of diffusion features

The visualization of the distribution of fiber bundles f as a spherical function is useful to indicate the orientations of the sub-voxel fiber bundles [START_REF] Jansons | Persistent angular structure: new insights from diffusion magnetic resonance imaging data[END_REF][START_REF] Seunarine | Linear persistent angular structure mri and non-linear spherical deconvolution for diffusion mri[END_REF]. This is often used as a preprocessing step in the context of using spherical deconvolution as an input to fiber tracking algorithms, as it reduces the angular estimation error (Descoteaux et al., 2009a;[START_REF] Savadjiev | Labeling of ambiguous subvoxel fibre bundle configurations in high angular resolution diffusion MRI[END_REF].

Advantages and limitations

+ When using a spherical harmonics basis the estimation is linear and is hence very fast.

+ In contrast to the mixture models, for which the number of compartiments has to be determined a priori (c.f . Section 3.1), this parameter can be obtained by simple thresholding of the reconstructed f distribution.

+ The non-Gaussian modeling of the angular profile enables the characterization of regions with crossing fibers.

+ The use of spherical harmonics for HARDI sampling generally requires approximately 60 samples on a sphere in the literature, which leads to a modest acquisition time [START_REF] Alexander | Maximum entropy spherical deconvolution for diffusion MRI[END_REF][START_REF] Tournier | Direct estimation of the fiber orientation density function from diffusion-weighted mri data using spherical deconvolution[END_REF].

-The choice of the deconvolution kernel h is empirical and is based on a priori assumptions about the data.

-As with DTI (c.f . Section 2), the Gaussian modeling of the radial profile is inadequate and this excludes several features including apparent cell size, axon diameter, etc.

High Angular ADC (non-parametric)

Local diffusion modeling The first experiments to characterize multiple fiber bundle configurations were based on the generalization of the modeling of the apparent diffusion coefficient (ADC) from a low to a high angular resolution, without the Gaussian assumption imposed by DTI [START_REF] Alexander | Detection and modeling of nongaussian apparent diffusion coefficient profiles in human brain data[END_REF][START_REF] Frank | Characterization of anisotropy in high angular resolution diffusion-weighted mri[END_REF][START_REF] Tuch | High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity[END_REF]. Recall that the normalized diffusion signal E is modeled by the Stejskal-Tanner equation [START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF]:

E(q) = exp -4π 2 τq 2 D(u) . (39)
In contrast to DTI, HARDI sampling allows high angular modeling of the ADC with functions of order higher than the second order diffusion tensor used in DTI. Thus, [START_REF] Descoteaux | Apparent diffusion coefficients from high angular resolution diffusion imaging: Estimation and applications[END_REF]; [START_REF] Özarslan | Generalized diffusion tensor imaging and analytical relationships between diffusion tensor imaging and high angular resolution diffusion imaging[END_REF] propose to model the ADC by a Higher Order Tensor (HOT):

D(u) = J l j=1 D j µ j l p=1 u j(p) , (40)
where the symbol J l = (l + 1)(l + 2)/2 corresponds to the number of terms of a l-order tensor, u j(p) is the p-th value of the j-th tensor component D j and µ j is the multiplicity index [START_REF] Özarslan | Generalized diffusion tensor imaging and analytical relationships between diffusion tensor imaging and high angular resolution diffusion imaging[END_REF].

In parallel, several studies suggested the use of spherical harmonics (SH) for angular ADC estimation [START_REF] Alexander | Detection and modeling of nongaussian apparent diffusion coefficient profiles in human brain data[END_REF][START_REF] Chen | Estimation, smoothing, and characterization of apparent diffusion coefficient profiles from high angular resolution dwi[END_REF][START_REF] Frank | Characterization of anisotropy in high angular resolution diffusion-weighted mri[END_REF][START_REF] Zhan | Circular spectrum mapping for intravoxel fiber structures based on high angular resolution apparent diffusion coefficients[END_REF]:

D(u) = J l j=0 D j y j (u), ( 41
)
where J l = (l + 1)(l + 2)/2 corresponds to the number of terms of an order l harmonic expansion (c.f . equation ( 132)), the symbols y j stands for the real and symmetric spherical harmonics, and D j are the weighting coefficients.

Remark. There is an mathematical equivalence between spherical expansions which use the same order of spherical harmonics (equation ( 40)) and cartesian tensors (equation ( 41)) [START_REF] Descoteaux | Apparent diffusion coefficients from high angular resolution diffusion imaging: Estimation and applications[END_REF][START_REF] Johnston | Cartesian tensor scalar product and spherical harmonic expansions in boltzmann's equation[END_REF][START_REF] Özarslan | Generalized diffusion tensor imaging and analytical relationships between diffusion tensor imaging and high angular resolution diffusion imaging[END_REF].

Model estimation from the data As with the DTI method detailed in equation ( 16), the high angular ADC coefficients D j result from the dot product between the log-normalized diffusion signal E and the higher order tensors:

D j = u∈S 2 - ln (E) 4π 2 τq 2 µ j l p=1 u j(p) du. (42)
Hence, the ADC is expressed in the spherical harmonics basis as:

D j = u∈S 2 - ln (E) 4π 2 τq 2 y j (u)du. ( 43
)
In practice, the estimation of the parameters D j of the HOT (equation ( 40)) and SH (equation ( 41)) models is solved by a linear least squares method because because this is simple and fast to compute [START_REF] Descoteaux | Apparent diffusion coefficients from high angular resolution diffusion imaging: Estimation and applications[END_REF], which can be seen as follows. Let us group the parameters into a vector D:

D T = D j 1≤ j≤J l . ( 44
)
The resulting system of linear equations can be written in the matrix form Y = MD, with E being the vector of samples given by equation ( 21) and M being the corresponding basis matrix, that is for the J-th order tensors

M = µ j l p=1 u i 1(p) 1≤i≤n q ,1≤ j≤J l , (45)
with n q being the number of samples. For a spherical harmonics basis, the basis matrix is written as:

M = E j (u i ) 1≤i≤n q ,1≤ j≤J l . (46)
Finally, the vector D of model parameters is given by the pseudoinverse of the matrix M:

D = (M T M) -1 M T E. ( 47
)
Other studies have proposed the use of estimation techniques that are more robust to noise, which can have a significant effect in MR images [START_REF] Gudbjartsson | The rician distribution of noisy mri data[END_REF], notably the use of regularized least squares [START_REF] Descoteaux | Apparent diffusion coefficients from high angular resolution diffusion imaging: Estimation and applications[END_REF] and of a variational estimation framework [START_REF] Chen | Estimation, smoothing, and characterization of apparent diffusion coefficient profiles from high angular resolution dwi[END_REF].

Whereas equation ( 47) enables a fast estimation, it does not ensure positivity of the vector D. Hence, in [START_REF] Barmpoutis | Symmetric positive 4 h order tensors & their estimation from diffusion weighted mri[END_REF], the estimation is expressed as a ternary quartic (a homogeneous polynomial of degree 4 with three variables), such that the Hilbert theorem [START_REF] Hilbert | Ueber die Darstellung definiter Formen als Summe von formenquadraten[END_REF] ensures positivity. Alternatively, Ghosh et al. (2008a) proposed the use of a Riemannian metric on the tensor space in order ensure positivity of the components of D.

Processing and extraction of diffusion features

Scalar features: DTI has become a standard method for clinical use. The indices (i.e. scalar features) based on higher order models therefore attempt to overcome the modeling limitations of DTI, particularly the detection of "crossing fibers" in brain tissue. In the following we enumerate three indices extracted from higher order models that are directly computed using the spherical harmonics or higher order tensor coefficients.

In the case of a spherical harmonics expansion of the diffusion signal, [START_REF] Frank | Characterization of anisotropy in high angular resolution diffusion-weighted mri[END_REF] demonstrates that the D coefficients can discriminate between the cases of isotropic diffusion (index F 0 ), anisotropic diffusion with a single intravoxel direction (index F 2 ) and anisotropic diffusion with several intravoxel directions (index F multi ):

F 0 = |D 0 | F 2 = j:l=2 |D j | F multi = j:l≥4 |D j | (48)
Some studies have proposed other scalar features based on a ratio between these three indices (F 0 , F 2 , F multi ) to obtain more accurate results [START_REF] Chen | Estimation, smoothing, and characterization of apparent diffusion coefficient profiles from high angular resolution dwi[END_REF][START_REF] Chen | Apparent diffusion coefficient approximation and diffusion anisotropy characterization in DWI[END_REF][START_REF] Frank | Characterization of anisotropy in high angular resolution diffusion-weighted mri[END_REF].

As mentioned in Section 2, fractional anisotropy (FA) is one of the most widely used scalar features of DTI, which motivates a generalization of this index to the higher order models. The generalized fractional anisotropy (GFA) is defined as the normalized variance of a spherical function, here the ADC, and is expressed for a discretization of the sphere in a set of n points by {u i } 1≤i≤n [START_REF] Tuch | Q-ball imaging[END_REF]:

GFA(D) = std(D) rms(D) = n n -1 n i=1 (D(k i ) -D ) 2 n i=1 D(k i ) 2 . ( 49
)
The GFA index can be directly computed using the spherical harmonics basis:

GFA(D) = 1 - a 2 0 J j=0 a 2 j . ( 50
)
Like the FA, the GFA is defined on the interval [0, 1]: a value of zero indicates a perfectly isotropic diffusion, while a value of one indicates a totally anisotropic diffusion. Özarslan et al. (2005) proposed an alternative generalization of the FA, for a signal estimation based on higher order tensors. The obtained index, namely the generalized anisotropy (GA), is defined by the generalized trace Tr gen and variance Var gen operators:

Tr gen (D) = 3 4π u∈S 2 D(u)du Var gen (D) = 1 3 Tr gen D 2 Tr gen (D) - 1 3 . ( 51
)
The GA is normalized in the interval [0, 1) using an ad hoc data-dependent mapping (Özarslan et al., 2005):

GA = 1 -1 + (250Var gen ) f (Var gen ) -1 f (Var gen ) = 1 + 1 + 5000Var gen -1 . (52)
Further details on this index are described in (Özarslan et al., 2005). Furthermore, a numerical comparison of selected scalar features of the high angular ADC is presented in [START_REF] Descoteaux | High angular resolution diffusion mri: from local estimation to segmentation and tractography[END_REF].

Vectorial features:

The modeling of the ADC profile by higher order models leads to a more accurate characterization of intravoxel anisotropy than that provided by DTI. However, as illustrated in Figure 9(a,b), the maxima of the estimated ADC do not coincide with the directions of the underlying fiber bundles (von dem [START_REF] Von Dem Hagen | Orientational diffusion reflects fiber structure within a voxel[END_REF]. Therefore the high angular resolution ADC can not be used directly to extract meaningful vectorial features, and a fortiori for fiber tracking in brain white matter.

Advantages and limitations

+ The non Gaussian modeling of the angular profile enables the characterization of brain regions with crossing fibers, which is a problematic case for the DTI [START_REF] Alexander | Detection and modeling of nongaussian apparent diffusion coefficient profiles in human brain data[END_REF][START_REF] Frank | Characterization of anisotropy in high angular resolution diffusion-weighted mri[END_REF][START_REF] Özarslan | Generalized scalar measures for diffusion MRI using trace, variance, and entropy[END_REF].

+ The use of parametric models for HARDI sampling generally requires approximately 60 samples on a sphere in the literature, which leads to a modest acquisition time [START_REF] Alexander | Maximum entropy spherical deconvolution for diffusion MRI[END_REF][START_REF] Tournier | Direct estimation of the fiber orientation density function from diffusion-weighted mri data using spherical deconvolution[END_REF]. However this number of samples is variable and depends on the model order in use.

The studies in the literature generally recommend an order equal to four or six, which is sufficient for the representation of crossing fibers [START_REF] Alexander | Detection and modeling of nongaussian apparent diffusion coefficient profiles in human brain data[END_REF][START_REF] Frank | Characterization of anisotropy in high angular resolution diffusion-weighted mri[END_REF].

-The fact that maxima of ADC profile are not aligned with the underlying fiber directions prevents the direct extraction of accurate vectorial features, and a fortiori the use of fiber tracking algorithms in brain white matter [START_REF] Von Dem Hagen | Orientational diffusion reflects fiber structure within a voxel[END_REF].

-As is the case with DTI (c.f . section 2), the Gaussian modeling of the radial profile is inadequate and this excludes a wide variety of features including apparent cell size, axon diameter, etc.

Q-Ball Imaging (non-parametric)

Local diffusion modeling [START_REF] Tuch | Q-ball imaging[END_REF] introduced the technique of Q-Ball Imaging (QBI), which is primarily aimed at extracting the diffusion orientation density function (ODF). The ODF is a spherical feature of the ensemble average propagator, which we shall explain in greater detail after describing the signal modeling process.

In practice, a few samples from a HARDI sampling can be used to infer a continuous approximation process using a spherical basis to reconstruct the signal. In its original version [START_REF] Tuch | Q-ball imaging[END_REF], the authors interpolate the diffusion signal with a spherical Radial Basis Function (sRBF) with a Gaussian kernel [START_REF] Fasshauer | Scattered data fitting on the sphere[END_REF]. Later work [START_REF] Anderson | Measurement of fiber orientation distributions using high angular resolution diffusion imaging[END_REF][START_REF] Descoteaux | Regularized, fast and robust analytical q-ball imaging[END_REF][START_REF] Hess | Q-ball reconstruction of multimodal fiber orientations using the spherical harmonic basis[END_REF] suggest the use of the spherical harmonics basis, which lowers the number of samples needed. In the latter basis, the normalized diffusion signal E, sampled by a HARDI acquisition on a sphere of radius q in the q-space, is expressed as:

E(q) = L l=0 l m=-l a lm y m l (u)δ(q -q ). ( 53
)
Several recent studies suggest the use of other spherical functions including: spherical ridgelets [START_REF] Michailovich | On approximation of orientation distributions by means of spherical ridgelets[END_REF], spherical wavelets [START_REF] Kezele | Multiresolution decomposition of HARDI and ODF profiles using spherical wavelets[END_REF][START_REF] Khachaturian | Boosting the sampling efficiency of q-ball imaging using multiple wavevector fusion[END_REF], or Watson, de von Mises and de la Vallé Poussin density functions [START_REF] Rathi | Directional functions for orientation distribution estimation[END_REF]. The diffusion signal expansion in these latter bases requires fewer samples, but at the expense of a stronger prior assumption on the signal (concentration parameter which narrows the support of the functions).

Model estimation from the data The coefficients a lm of the signal E in the real spherical harmonic basis y m l of order l are the results of the following projection:

a lm = E, y m l = u∈S 2 E(u)y m l (u)du. ( 54
)
In practice, the (l + 1)(l + 2)/2 system of equations ( 54) is overdetermined, and a least squares estimation method is used. [START_REF] Poupon | Real-time MR diffusion tensor and Q-ball imaging using Kalman filtering[END_REF] tackles the issue of "real-time" QBI estimation, during the repetition time (TR) of the diffusion MRI acquisition sequence. At each new acquisition, the authors recursively estimate the incomplete set of images by a Kalman filter, for which [START_REF] Deriche | Optimal real-time q-ball imaging using regularized kalman filtering with incremental orientation sets[END_REF] proposes an optimal resolution for Q-Ball Imaging.

Following the work of [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF], which uses compressed sensing to reduce the number of required samples in the k-space and the scanning time a fortiori, several studies recently propose the use of compressed sensing in the q-space [START_REF] Landman | Resolution of crossing fibers with constrained compressed sensing using traditional diffusion tensor mri[END_REF][START_REF] Lee | Compressed sensing based diffusion spectrum imaging[END_REF][START_REF] Menzel | Accelerated diffusion spectrum imaging in the humain brain using compressed sensing[END_REF][START_REF] Merlet | Compressed Sensing for Accelerated EAP Recovery in Diffusion MRI[END_REF]Michailovich and Rathi, 2010;[START_REF] Michailovich | On approximation of orientation distributions by means of spherical ridgelets[END_REF]Michailovich et al., , 2010;;[START_REF] Yin | Bregman iterative algorithm for l-1 minimization with application to compressed sensing[END_REF]. These estimation methods assume that the diffusion signal is sparse in the basis of modeling functions, i.e., that most coefficients are equal to zero. Thus, the ideal estimation would minimize the number of non-zero coefficients (also known as the L 0 norm). However, solving this problem is not computationally feasible since it is non-deterministic polynomial-time hard (NP-hard). Estimation methods generally approximate this problem by minimizing the L 1 norm instead of the L 0 norm, e.g., [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. Information Theory[END_REF][START_REF] Donoho | De-noising by soft-thresholding[END_REF].

Processing and extraction of diffusion features

Scalar Features: These are the same as the high angular ADC features which have been detailed in Section 3.3.

Spherical Features: Equation ( 55) links the diffusion orientation density function (ODF) that QBI aims to compute to the ensemble average propagator (EAP) formalism. The computation of the ODF overcomes the limits of ADC by correctly recovering the direction of crossing fibers, and thus is very often used as a preprocessing step for fiber tracking [START_REF] Anwander | Connectivity-based parcellation of Broca's area[END_REF][START_REF] Lenglet | Mathematical methods for diffusion mri processing[END_REF][START_REF] Savadjiev | Labeling of ambiguous subvoxel fibre bundle configurations in high angular resolution diffusion MRI[END_REF].

The exact diffusion ODF is defined as the radial projection of the diffusion EAP P on the unit sphere. Let k ∈ S 2 be a point on this sphere. The ODF is then expressed in the following general form [START_REF] Canales-Rodríguez | Diffusion orientation transform revisited[END_REF]:

ODF κ (k) = 1 Z κ ∞ p=0 P(pk)p κ d p, ( 55
)
where κ is the order of the radial projection and Z n is a normalization constant such that 1/Z κ k∈S 2 ODF κ (k) = 1. The classical ODF was introduced in the QBI method by [START_REF] Tuch | Q-ball imaging[END_REF] as the zero-order radial projection, and will be denoted as ODF 0 . Similarly, the second order radial projection will be referred to as ODF 2 . It has been used in other methods including the Diffusion Spectrum Imaging method (DSI, section 5.1), the Diffusion Orientation Transform (DOT, section 3.5) and the Spherical Polar Fourier expansion (SPF, section 5.8).

Despite the fact that all the above ODF definitions capture angular information of the diffusion propagator, each one leads to slightly different results. Higher values of order κ will favour a longer displacement length.

The classical ODF 0 can be reformulated in the q-space as [START_REF] Tuch | Q-ball imaging[END_REF]:

ODF 0 (k) = 1 Z 0 p∈R 3 P(p) δ 1 -r T k dp (56) = q∈R 3 E(q) δ u T k dq (57)
It is important to remark that the computation of the exact ODF 0 requires the value of the signal E in the whole q-space.

Since HARDI limits the samples to lie on a sphere of radius q , the QBI method cannot compute the exact ODF 0 but rather computes an approximate ODF 0 . More precisely, the QBI method assumes that the diffusion signal is equal to zero outside the HARDI sampling sphere in the q-space, and retrieves an approximate ODF using the Fourier transform, which is equivalent in this particular setting to the Funk-Radon transform FRT q [START_REF] Tuch | Q-ball imaging[END_REF]:

FRT q (k) = u∈S 2 E(q u) δ(u T k) du = 2πq p∈R 3 P(p)J 0 (2πq p) δ 1 -r T k dp, (58)
with k ∈ S 2 , and the symbol J 0 representing the zeroth order Bessel function of the first kind [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]:

J 0 (z) = ∞ n=0 (-1 n ) z 2n 4 n (n!) 2 . ( 59
)
Thus, the ODF given from the QBI equation ( 58) is an approximation of the exact ODF defined in the equation ( 55), i.e. ODF ≈ FRT q . Some studies propose the use of the spherical harmonics basis for the approximation of the diffusion signal [START_REF] Anderson | Measurement of fiber orientation distributions using high angular resolution diffusion imaging[END_REF][START_REF] Descoteaux | Regularized, fast and robust analytical q-ball imaging[END_REF][START_REF] Hess | Q-ball reconstruction of multimodal fiber orientations using the spherical harmonic basis[END_REF], which interestingly simplifies the computation of the Funk-Radon transform as stated by the Funk-Hecke theorem [START_REF] Andrews | Special Functions[END_REF]:

FRT q (k) = L l=0 l m=-l 2πP l (0)a lm y m l (k), (60)
The Funk-Radon formulation of equation ( 60) is interesting since the computation involves only a matrix multiplication between the coefficient vectors (a lm ) and (2πP l (0)). Note that similar analytical Funk-Radon relationships to the equation (60) were proposed for the spherical ridgelets functions [START_REF] Michailovich | On approximation of orientation distributions by means of spherical ridgelets[END_REF] and the Watson, de von Mises and de la Vallé Poussin density functions (for approximate ODFs) [START_REF] Rathi | Directional functions for orientation distribution estimation[END_REF].

Extraction of fiber bundle directions: As in the context of fiber tracking in DTI, fiber bundles are usually reconstructed by following the direction of the principle eigenvector of the second-order tensor. When using the QBI method or other models capable of representing complex geometries, such as fiber crossings (illustrated in Figure 9), the extraction of fiber bundle directions requires post processing.

The numerical lookup of the approximated ODF maxima, discretized on a grid of n points {k i } i≤n distributed on the sphere for equation (60), requires significant computing time and the precision is dependent on the chosen number of points n. More recently, [START_REF] Bloy | On computing the underlying fiber directions from the diffusion orientation distribution function[END_REF]; Ghosh et al. (2008b) have proposed an analytical space reduction for the lookup of the maxima, from the higher order tensor coefficients. These maxima are expressed as fixed points of homogeneous polynomials. Ghosh et al. (2008b) have suggested that these methods can be also applied to spherical harmonics coefficients.

Once the ODF maxima are extracted for each voxel of the image, fiber tracking is carried out on the maxima field using three dimensional curves [START_REF] Savadjiev | 3D curve inference for diffusion MRI regularization and fibre tractography[END_REF]. Some studies have proposed a validation of the fiber tracking based on the QBI method, using a comparison with fiber tracking based on DTI [START_REF] Campbell | Flow-based fiber tracking with diffusion tensor and q-ball data: validation and comparison to principal diffusion direction techniques[END_REF][START_REF] Perrin | Validation of q-ball imaging with a diffusion fibrecrossing phantom on a clinical scanner[END_REF].

Advantages and limitations

The original QBI method [START_REF] Tuch | Q-ball imaging[END_REF] assumes that P(p) ≈ P(p)J 0 (2πq p). This approximation of the diffusion propagator leads to the corruption of the neighbourhood of direction k by the Bessel function J 0 , which narrows in extent as the value of q grows [START_REF] Tuch | Q-ball imaging[END_REF]. As a consequence, the obtained ODF 0 using the QBI method generally has a low angular contrast.

To solve this problem, the methods in the literature traditionally use an ad hoc postprocessing step, such as min-max normalization [START_REF] Tuch | Q-ball imaging[END_REF] or a deconvolution with the Laplace-Beltrami operator (defined in section 7.3.2) (Descoteaux et al., 2009a). Aganj et al. (2009a) has suggested the use of a secondorder ODF 2 in place of the zero-order ODF 0 which suppresses the need for such an ad hoc normalization (more details in section 3.5). The advantages and limitations of QBI are summarized below.

+ The estimation of the angular diffusion profile is carried out without any biologic prior assumptions on the signal.

+ The linear estimation algorithm of the ODF in the spherical harmonics basis is very fast -As in the case of high angular ADC, described in Section 3.3, QBI requires approximately 60 samples on the sampling sphere.

-The method allows the extraction of only a single spherical feature.

-The estimation of the ODF 0 as proposed in the original QBI method [START_REF] Tuch | Q-ball imaging[END_REF]) is an approximation of the exact ODF. The corruption arises from the inappropriate Dirac modeling of the radial decay of the signal.

Diffusion Orientation Transform (DOT) (non-parametric)

Local diffusion modeling Özarslan et al. (2006) proposed the Diffusion Orientation Transform (DOT) method which computes the isoradius of the diffusion propagator P(p k), with radius p ∈ R + and k ∈ S 2 a point on the unit sphere. This results in a spherical feature of the diffusion, however it is important to note that the isoradius is different from the ODF feature proposed by the QBI method (detailed in the previous section). Nonetheless, several groups recently proposed to estimate an approximation of the ODF from the DOT modeling (Aganj et al., 2009a;[START_REF] Canales-Rodríguez | Diffusion orientation transform revisited[END_REF][START_REF] Vega | Estimation of fiber orientation probability density functions in high angular resolution diffusion imaging[END_REF]. We shall explain these features in greater details after describing the signal modeling process.

Assuming that the radial diffusion follows a Gaussian distribution, the normalized diffusion signal E is modeled by the Stejskal-Tanner equation (Basser et al., 1992;[START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF]):

E(q) = exp -4π 2 τq 2 D(u) . (61)
However, Özarslan et al. (2006) propose to estimate the Fourier-Bessel transform I l of the diffusion signal, which is the radial part of the Fourier transform decomposed as a plane wave expansion into the spherical harmonic basis and spherical Bessel of order l:

I l (u) = 4π ∞ q=0 E(q) j l (2πqp )q 2 dq avec j l (z) = (-1) l z l d zdz l sin z z . ( 62
)
Here the symbols j l stand for the spherical Bessel functions [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]. The use of the Gaussian prior on the radial diffusion profile, as expressed in the equation ( 61), enables one to solve the equation ( 62) (Özarslan et al., 2006) analytically:

I l (u) = R l 0 Γ l+3 2 1F1 l+3 2 ; l + 3 2 ; R 2 0 ln(E(u))(πq) -2 2 l+3 π 3/2 -ln(E(u))4 -1 (πq) -2 (l+3)/2 Γ(l + 3/2) , (63)
with 1F1 the confluent hypergeometric function of the first kind [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]:

1F1(a; b; x) = ∞ k=0 (a) k x k (b) k k! with (a) k = Γ(a + k) Γ(a) , (64)
where (a) k stands for the Pochhammer symbol [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF].

Model estimation from the data The a lm coefficient resulting from the expansion of equation ( 63) in the spherical harmonics basis y m l is expressed as:

a lm = I l , y m l = u∈S 2 I l (u)y m l (u)du. ( 65
)
The system of (l + 1)(l + 2)/2 equations given by equation ( 65) can be solved by a linear least square method. Özarslan et al. (2006) suggest to speed up the estimation of equation ( 63) by using a table of precomputed values of the confluent hypergeometric 1F1 function defined in equation ( 64).

Processing and extraction of diffusion features

Scalar Features: These are the same as the high angular ADC features (detailed in Section 3.3).

Spherical features:

As mentioned in section 3.3, the ADC formalism does not allow for the recovery of fiber bundle orientations.

On the other hand, the classical DOT method, as proposed by Özarslan et al. (2006), uses the diffusion EAP P in order to extract the isoradius, defined as the surface P(p k) with k ∈ S 2 at a constant arbitrary distance p ∈ R. This feature can be expressed directly using the signal coefficients of equation ( 65) in the spherical harmonics basis y m l (Özarslan et al., 2006):

P(p k) = L l=0 l m=-l a lm (-i) l y m l (k)dp with k ∈ S 2 . ( 66
)
The estimation of equation ( 66 2009) have proposed to extend the DOT method to approximate the ODF. As mentioned for the QBI method (Section 3.4), the exact resolution of the {ODF κ } κ={0,2} requires a complete sampling of the q-space, which is not available for HARDI (spherical sampling). These methods use a single Gaussian prior on the radial diffusion profile E. We detail the more general case of a multi-Gaussian prior in Section 5.3.

Advantages and limitations

+ The estimation of the angular diffusion profile is carried out without any prior biological assumptions about the signal.

+ The linear estimation algorithm of the isoradius in the spherical harmonics basis is fast.

-As with the case of the high angular ADC, described in Section 3.3, QBI requires approximately 60 samples on the sampling sphere.

-The precision of this method is limited by the mono-Gaussian assumption of the decay of radial diffusion, an assumption that does not always hold [START_REF] Assaf | In vivo and in vitro bi-exponential diffusion of nacetyl aspartate (naa) in rat brain: a potential structural probe?[END_REF][START_REF] Callaghan | Diffractionlike effects in nmr diffusion studies of fluids in porous solids[END_REF][START_REF] Niendorf | Biexponential diffusion attenuation in various states of brain tissue: implications for diffusion-weighted imaging[END_REF][START_REF] Tanner | Restricted self-diffusion of protons in colloidal systems by the pulsed-gradient, spin-echo method[END_REF]).

Radial Reconstruction in Diffusion MRI

Introduction

Radial diffusion analysis by Nuclear Magnetic Resonance (NMR), or q-space imaging [START_REF] Callaghan | Principles of Nuclear Magnetic Resonance Microscopy[END_REF][START_REF] Cory | Measurement of translational displacement probabilities by nmr: an indicator of compartmentation[END_REF], favours a dense q-space sampling at the expense of the k-space sampling (see section 1.3 for details on k-space).

The dMRI data can be acquired over a large volume, such as an entire human brain. Diffusion NMR data, on the other hand, is acquired for a much smaller volume, such as a specimen enclosed in a test tube. The resulting image in diffusion NMR is composed of a single voxel. Figure 10: A plot of the radial diffusion in an environment composed of red blood cells. It is important to note that the radial decay is not Gaussian, but a "diffusion-diffraction" pattern occurs [START_REF] Callaghan | Diffractionlike effects in nmr diffusion studies of fluids in porous solids[END_REF][START_REF] Cory | Measurement of translational displacement probabilities by nmr: an indicator of compartmentation[END_REF][START_REF] Kuchel | Nmr diffusion-diffraction of water revealing alignment of erythrocytes in a magnetic field and their dimensions and membrane transport characteristics[END_REF][START_REF] Regan | NMR studies of diffusion-coherence phenomena in red cell suspensions: current status[END_REF] where the position of the minima are related to cell size and extracellular average size. (a) A micrograph of red blood cells (adapted from Noguchi, Rodjers and Schechter, NIDDK, under free licence). (b)

Signal intensity E as a function of q (adapted from [START_REF] Regan | NMR studies of diffusion-coherence phenomena in red cell suspensions: current status[END_REF].)

Applications of diffusion NMR include the study of the microscopic geometry of mineralogical environments (sedimentary rock, ground, cement) and biological environments (proteins, lungs, the interstitial space of the skin). In theory, since the average diffusion length during typical measurement duration is close to a cell size, the restricted diffusion in these confined environments can accurately characterize properties of the underlying micro-architecture such as its porosity and tortuosity. Traditionally, this characterization is analyzed with a two dimensional plot arising from sampling along a line in q-space [START_REF] Callaghan | Principles of Nuclear Magnetic Resonance Microscopy[END_REF][START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF]. The mathematical methods to analyze this plot are divided into three groups:

1. The methods of the first group take an analytical approach and propose a exact expression of the diffusion signal in q-space (i.e., the Bloch-Torrey equation with the q-space formalism). Nonetheless, these methods only characterize diffusion in simple models: infinite parallel planes, cylinders and spheres [START_REF] Brownstein | Importance of classical diffusion in nmr studies of water in biological cells[END_REF][START_REF] Robertson | Spin-echo decay of spins diffusing in a bounded region[END_REF][START_REF] Tanner | Transient diffusion in a system partitioned by permeable barriers. application to nmr measurements with a pulsed field gradient[END_REF][START_REF] Tanner | Restricted self-diffusion of protons in colloidal systems by the pulsed-gradient, spin-echo method[END_REF], or bi-compartment models (having a bi-Gaussian distribution) [START_REF] Clark | Water diffusion and anisotropy at high b values in the human brain[END_REF][START_REF] Le Bihan | from the diffusion coefficient to the diffusion tensor[END_REF][START_REF] Niendorf | Biexponential diffusion attenuation in various states of brain tissue: implications for diffusion-weighted imaging[END_REF]. 2. The methods of the second group extend to more complex diffusion situations, which do not have a known analytical solutions. The equations describing the behaviour of the diffusion signal are solved using numerical schemes, namely either: (i) resolution of the Bloch-Torrey equation with finite differences schemes [START_REF] Putz | Edge enhancement by diffusion in microscopic magnetic resonance imaging[END_REF][START_REF] Wayne | Nuclear-magnetic-resonance study of self-diffusion in a bounded medium[END_REF][START_REF] Zielinski | Relaxation of nuclear magnetization in a nonuniform magnetic field gradient and in a restricted geometry[END_REF] or (ii) by simulating a Brownian motion with the Monte-Carlo method [START_REF] Coy | Pulsed gradient spin echo nuclear magnetic resonance for molecules diffusing between partially reflecting rectangular barriers[END_REF][START_REF] Grebenkov | Restricted diffusion in a model acinar labyrinth by NMR: Theoretical and numerical results[END_REF][START_REF] Hyslop | Effects of restricted diffusion on microscopic NMR imaging[END_REF][START_REF] Kuchel | Nmr diffusion-diffraction of water revealing alignment of erythrocytes in a magnetic field and their dimensions and membrane transport characteristics[END_REF][START_REF] Mendelson | Percolation model of nuclear magnetic relaxation in porous media[END_REF][START_REF] Mitra | Effects of finite gradient-pulse widths in pulsedfield-gradient diffusion measurements[END_REF]. 3. Unlike the methods of the two previous groups which approximate simulations of the diffusion signal with biologi-cal models, the methods of the last group propose a functional analysis of the in vivo diffusion signal, e.g., by using a discrete Fourier transform [START_REF] Callaghan | Diffractionlike effects in nmr diffusion studies of fluids in porous solids[END_REF][START_REF] Cory | Measurement of translational displacement probabilities by nmr: an indicator of compartmentation[END_REF], the Lévy distribution and Kohlrausch-Williams-Watts function [START_REF] Köpf | Anomalous diffusion of water in biological tissues[END_REF][START_REF] Köpf | NMR studies of anomalous diffusion in biological tissues: experimental observation of Lévy stable processes[END_REF], or cumulant expansion [START_REF] Callaghan | Principles of Nuclear Magnetic Resonance Microscopy[END_REF][START_REF] Mitra | Effects of microgeometry and surface relaxation on nmr pulsed-field-gradient experiments: Simple pore geometries[END_REF][START_REF] Stepisnik | Analysis of NMR self-diffusion measurements by a density matrix calculation[END_REF].

For a more detailed overview of these methods, we refer the reader to [START_REF] Grebenkov | NMR survey of reflected Brownian motion[END_REF]. Diffusion NMR studies have revealed that the radial profile of diffusion in brain structures is non-Gaussian [START_REF] Assaf | In vivo and in vitro bi-exponential diffusion of nacetyl aspartate (naa) in rat brain: a potential structural probe?[END_REF][START_REF] Callaghan | Diffractionlike effects in nmr diffusion studies of fluids in porous solids[END_REF][START_REF] Cory | Measurement of translational displacement probabilities by nmr: an indicator of compartmentation[END_REF][START_REF] Niendorf | Biexponential diffusion attenuation in various states of brain tissue: implications for diffusion-weighted imaging[END_REF][START_REF] Tanner | Transient diffusion in a system partitioned by permeable barriers. application to nmr measurements with a pulsed field gradient[END_REF][START_REF] Tanner | Restricted self-diffusion of protons in colloidal systems by the pulsed-gradient, spin-echo method[END_REF]. A limitations of these studies, though, is that they typically assume the diffusion signal to be isotropic in its angular profile. Numerous studies have shown that this hypothesis is clearly not valid in the brain (c.f . Section 2 and 3). In the following we describe two radial reconstruction methods in the literature, the diffusion kurtosis imaging [START_REF] Jensen | Diffusional kurtosis imaging: the quantification of non-gaussian water diffusion by means of magnetic resonance imaging[END_REF] and the simple harmonic oscillator (Özarslan et al., 2008). et al. (2005) proposed the Diffusion Kurtosis Imaging (DKI) method, which approximates the diffusion radial profile by a series of cumulants ln (E(q)) = -D app (2π) 2 τq 2 + 1 6 D 2 app K app (2π) 4 τ 2 q 4 + O(q 5 ) (67)

Diffusion Kurtosis Imaging

Local diffusion modeling

Jensen

where the symbols D app and K app refer to the Apparent Diffusion Coefficient (ADC) and the Apparent Kurtosis Coefficient (AKC), respectively. The diffusion time is denoted τ = ∆ -δ/3. The expansion of equation ( 67) is a special case of the Generalized DTI (GDTI) method restricted to a fourth-order isotropic expansion (88) (see Section 5.4).

Model estimation from the data

Jensen et al. ( 2005) fit the diffusion signal on a voxel-by-voxel basis to the model described in the following formula, using the non-linear least square Levenberg-Marquardt method:

E(q) =        η 2 E(0) 2 + -D app (2π) 2 τq 2 + 1 6 D 2 app K app (2π) 4 τ 2 q 4 2        1/2 , ( 68
)
where η is the background noise, which is estimated a priori as the mean signal intensity in air.

Processing and extraction of diffusion features

Scalar Features: In [START_REF] Jensen | Diffusional kurtosis imaging: the quantification of non-gaussian water diffusion by means of magnetic resonance imaging[END_REF], the authors compute an image of the apparent diffusion coefficient D app and kurtosis K app . The latter quantifies the degree to which the water diffusion is non-Gaussian. It is hypothesized that the kurtosis K app can be altered by the space between the barriers and their porosity (also refereed to as "water exchange").

Advantages and limitations

+ The apparent diffusion kurtosis coefficient quantifies the degree to which the water diffusion is non-Gaussian, and thus reveals brain microstructure information hidden to DTI;

+ The number of samples in the q-space required for the estimation of the apparent diffusion kurtosis coefficient is small (minimum of 6 samples) [START_REF] Jensen | Diffusional kurtosis imaging: the quantification of non-gaussian water diffusion by means of magnetic resonance imaging[END_REF];

-The Diffusion Kurtosis Imaging (DKI) method does not consider the angular profile of the diffusion signal, and essentially assumes that it is isotropic in direction. This assumption is clearly violated in the human brain [START_REF] Moseley | Diffusion-weighted MR imaging of anisotropic water diffusion in cat central nervous system[END_REF] and may lead to erroneous interpretation of the results.

Simple Harmonic Oscillator

Local diffusion modeling

Özarslan et al. (2008) propose an approximation of the diffusion radial profile in the basis of Φ n functions, which is made up of Hermite polynomials H n :

E(q) = N-1 n=0 a n Φ n (q, α), ( 69
) with Φ n (q, α) = i n (2 n n!) -1/2 exp -2π 2 q 2 α 2 H n (2παq). ( 70
)
where α ∈ R stands for the scale factor of the basis Φ n functions.

Model estimation from the data

The coefficients a n of the diffusion signal E in the basis of functions Φ n are obtained with a dot product:

a n = E, Φ n = ∞ q=0 E(q)Φ n (q, α)dq. ( 71
)
In practice, Özarslan et al. (2008) suggest that the coefficients a n should be estimated using least squares techniques.

Processing and extraction of diffusion features

Scalar Features: The m-order moments of the ensemble average propagator (EAP) P can be expressed as (Özarslan et al., 2008):

p m = α m N-1 k=0,2,... (k + m -1)!! k! N-k-1 l=0,2,... (-1) l/2 2 k-l (k + l)! (l/2)! a k+l ( 72
)
where k is even is m is even, otherwise k is odd. Besides, the functions Φ n are eigenfunctions of the Fourier transform, which allows the radial EAP to be directly expressed as a function of the signal coefficients a n :

P(p) = N-1 n=0 (-i) n √ 2πα a n Φ n p, (2πα) -1 . ( 73
)

Advantages and limitations

+ The radial diffusion profile is estimated without prior assumptions.

+ The diffusion EAP and its moments can be estimated using linear methods, which are very fast.

+ The Gauss-Hermite functions can approximate complex signals with relatively few samples of the q-space (33 according to (Özarslan et al., 2008)).

-This method does not consider the angular profile of the diffusion signal, and essentially assumes that it is isotropic in direction. This assumption is clearly violated in the human brain [START_REF] Moseley | Diffusion-weighted MR imaging of anisotropic water diffusion in cat central nervous system[END_REF].

Combining Angular and Radial Reconstruction

Diffusion Spectrum Imaging (DSI)

Local diffusion modeling Diffusion Spectrum Imaging (DSI) is based on the direct sampling of the diffusion signal E, leading to the EAP P using the q-space formalism [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF][START_REF] Hagmann | Diffusion spectrum imaging tractography in complex cerebral white matter: an investigation of the centrum semiovale[END_REF][START_REF] Kuo | Optimization of diffusion spectrum imaging and q-ball imaging on clinical mri system[END_REF][START_REF] Lin | Validation of diffusion spectrum magnetic resonance imaging with manganese-enhanced rat optic tracts and ex vivo phantoms[END_REF][START_REF] Tuch | Diffusion mri of complex tissue structure[END_REF][START_REF] Wedeen | Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging[END_REF][START_REF] Wedeen | Mapping fiber orientation spectra in cerebral white matter with fourier transform diffusion mri[END_REF][START_REF] Wu | Computation of Diffusion Function Measures in q-Space Using Magnetic Resonance Hybrid Diffusion Imaging[END_REF][START_REF] Wu | Hybrid diffusion imaging[END_REF]. Since the diffusion signal sampling is carried out using practically the entire q-space in three dimensions, DSI does not require local diffusion modeling. Nevertheless, in practice, the diffusion signal E samples are multiplied by an observation window to ensure a smooth decay of E at high values of the diffusion wave-vector norm q = ||q||. In [START_REF] Wedeen | Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging[END_REF], the authors use the Hann window, which is expressed as [START_REF] Blackman | The measurement of power spectra: from the point of view of communications engineering[END_REF]:

f (q) = 1 2 1 + cos πq 2a , ( 74
)
where the symbol a ∈ R stands for the full width at half maximum of the window (FWHM) f .

Model estimation from the data As mentioned above, DSI does not require modeling of the local diffusion so there is no model estimation step.

Processing and extraction of diffusion features

The diffusion propagator P is related to the diffusion signal E by the three dimensional Fourier transform F 3D [START_REF] Stejskal | Use of spin echoes in a pulsed magnetic-field gradient to study anisotropic, restricted diffusion and flow[END_REF]:

P(p) = q∈Ω q E(q) exp -i2πq T p dq. ( 75
)
In practice, [START_REF] Wedeen | Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging[END_REF] use the discrete threedimensional Fourier transform. Figure 11 illustrates the DSI method for a slice of size 64 × 64 × 1 taken from a human brain, acquired in 25 minutes with the parameters shown. The majority of the features defined in previous sections can be computed numerically using the DSI method, e.g., anisotropy indices, fiber bundle directions [START_REF] Wedeen | Diffusion spectrum magnetic resonance imaging (DSI) tractography of crossing fibers[END_REF], apparent axon diameter [START_REF] Weng | Maturation-dependent microstructure length scale in the corpus callosum of fixed rat brains by magnetic resonance diffusion-diffraction[END_REF], etc. Notably, using DSI one can compute the exact second-order ODF 2 as opposed to the approximate ODF 2 of the DOT method (section 3.5) or the zero-order ODF 0 of QBI (c.f . section 3.4) .

Advantages and limitations

Given that the origin of the diffusion process is still debated, DSI offers the advantage that the diffusion propagator P is estimated without a model or prior assumptions about the angular or radial diffusion profiles.

Although DSI is one of the most accurate methods, two major issues make it unfit for clinical use:

1. DSI requires a very long acquisition time for the entire qspace to be sampled with sufficient density [START_REF] Wedeen | Diffusion spectrum magnetic resonance imaging (DSI) tractography of crossing fibers[END_REF]. 2. DSI requires very high magnetic gradients to explore the q-space, e.g., b = 17000s/mm 2 according to [START_REF] Wedeen | Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging[END_REF].

In the literature, these two limitations restrict the use of DSI to animal ex vivo studies for which the acquisition time and the strength of the gradients are not limited. Alternatively, it is possible to limit the total acquisition time by reducing the kspace sampling, which is equivalent to scaling down the spatial resolution of the image. Typically, [START_REF] Wedeen | Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging[END_REF] suggests the use of a voxel size equal to 3.6 × 3.6 × 3.6 mm 3 , versus 2 × 2 × 2 mm 3 for current HARDI acquisitions [START_REF] Tournier | Resolving crossing fibres using constrained spherical deconvolution: validation using diffusion-weighted imaging phantom data[END_REF], leading to a total acquisition time of 25 min for an image size equal to 64 × 64 × 1. We summarize the advantages and limitations of DSI below.

+ The radial and angular diffusion profiles are estimated without prior assumptions.

+ All the features of the local diffusion propagator can be computed.

-The hardware requirements to image the brain image lead to a very long acquisition time (typically several hours), which is incompatible with clinical applications.

-The use of the FFT algorithm imposes a rectangular sampling of the q-space, which is typically not suited to MRI scanners.

The remaining methods we describe utilize sparse or hybrid sampling of the q-space [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF][START_REF] Wu | Hybrid diffusion imaging[END_REF], in contrast to the dense sampling of the DSI method. In the literature, this sampling is traditionally carried out on several spheres of different radii, discretized in q-space. This sparse sampling enables both an angular and radial analysis of the diffusion signal, as opposed to the methods of the previous sections (Sections 2, 3 and 4). Whereas we refer to "sparse sampling" in a broad sense, to our knowledge there is no mention of sampling other than multi-spherical sampling in the literature. We refer the reader to (Assemlal et al., 2009b) for a multi-sphere sampling model which unifies most of the existing sampling methods in the literature.

Hybrid Diffusion Imaging

Local diffusion modeling Hybrid diffusion imaging [START_REF] Wu | Computation of Diffusion Function Measures in q-Space Using Magnetic Resonance Hybrid Diffusion Imaging[END_REF] uses the raw diffusion data acquired in several spheres of the q-space, without any modeling or interpolation of the local diffusion signal.

Model estimation from the data Since this method does not use any local diffusion model, the diffusion propagator is obtained with the Fast Fourier Transform (FFT). Since the FFT requires the data to lie on a rectangular grid [START_REF] Wu | Computation of Diffusion Function Measures in q-Space Using Magnetic Resonance Hybrid Diffusion Imaging[END_REF]; [START_REF] Wu | Hybrid diffusion imaging[END_REF] propose to adjust the diffusion signal to a 9 × 9 × 9 cartesian grid, using bilinear interpolation. On the other hand, [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF] interpolate the signal with a bi-Gaussian parametric model (CHARMED), as detailed in Section 3.1. [START_REF] Wu | Computation of Diffusion Function Measures in q-Space Using Magnetic Resonance Hybrid Diffusion Imaging[END_REF] have also introduced an alternative approach, which does not require any adjustment of the diffusion signal to a cartesian grid. The authors propose the computation of three features of the data samples, two of which are scalar and one of which is spherical.

Processing and extraction of diffusion features

Scalar Features: The return to zero probability P(0) is defined as [START_REF] Wu | Computation of Diffusion Function Measures in q-Space Using Magnetic Resonance Hybrid Diffusion Imaging[END_REF]:

P(0) = p∈R 3 P(p)δ(p)dp = q∈R 3 E(q)dq. ( 76
)
The root mean square (RMS) is defined as [START_REF] Wu | Computation of Diffusion Function Measures in q-Space Using Magnetic Resonance Hybrid Diffusion Imaging[END_REF]:

RMS = p T p = p∈R 3 P(p)||p|| 2 dp = ∞ p=0 P s (p)p 2 dp = Var (P m ) , (77)
with P m (p) the mean probability of P on the sphere of radius p. According to [START_REF] Wu | Computation of Diffusion Function Measures in q-Space Using Magnetic Resonance Hybrid Diffusion Imaging[END_REF], the formulation of equation ( 77) in the q-space depends on the geometric mean (GM) of the signal E on the sphere of radius q: (c) DSI acquisition sequence parameters [START_REF] Wedeen | Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging[END_REF].

RMS = Var (F 1D [GM(E)]) GM(E)(q) = exp        1 n n i=1 ln (E[qu i ])        (78)
Figure 11: The DSI method for imaging in the coronal plane of a human brain. Note the non-Gaussian profile of the diffusion propagator P.

Adapted from [START_REF] Wedeen | Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging[END_REF].

where F 1D is the unidimensional Fourier transform and the family {qu i } 1≤i≤n constitutes the set of n sampling points of the sphere of radius q.

Spherical Features: The exact ODF 2 is defined by equation (56).

Tridimensional Features: [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF]; [START_REF] Wu | Hybrid diffusion imaging[END_REF] also compute the diffusion propagator P with the Fourier transform of E, first preprocessing the data with an interpolation step, as described above.

Advantages and limitations

The Fourier transform applied to the raw signal is quite sensitive to numerical inaccuracies since it is not based on a continuous approximation of the diffusion signal and the set of available samples is small. + The sparse sampling proposed in the literature results in similar acquisition times to HARDI (typically less than 30 minutes) and similar numbers of samples in the q-space (approximately 130) [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF][START_REF] Wu | Computation of Diffusion Function Measures in q-Space Using Magnetic Resonance Hybrid Diffusion Imaging[END_REF].

+ [START_REF] Wu | Computation of Diffusion Function Measures in q-Space Using Magnetic Resonance Hybrid Diffusion Imaging[END_REF] have proposed the computation of three features of the diffusion propagator, which characterize properties of both the angular and the radial profiles of the signal.

-The numerical accuracy of the computations suffers from thead hoc signal interpolation used to retrieve the diffusion propagator [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF][START_REF] Wu | Hybrid diffusion imaging[END_REF].

Generalized ODF

Local diffusion modeling As mentioned in our discussion of the QBI and DOT methods in Sections 3.4 and 3.5, the estimation of the exact ODF κ of any order κ requires a full radial sampling. This is not the case for HARDI sampling and therefore only an approximate ODF κ can be computed. Aganj et al. (2009b) model the diffusion signal E with N Gaussian functions α k weighted by the coefficients λ k :

E(q) = N k=1 λ k (u)α k (u) q 2 with α k (u) q 2 = exp q 2 ln α k (u) , (79)
where N ≤ (M -1)/2, with M standing for the number of sampling spheres in the q-space. This multi-Gaussian approximation (b) Parameters of the acquisition sequence of the HYDI method, proposed by [START_REF] Wu | Computation of Diffusion Function Measures in q-Space Using Magnetic Resonance Hybrid Diffusion Imaging[END_REF][START_REF] Wu | Hybrid diffusion imaging[END_REF].

Figure 12: Parameters of two acquisition sequences, resulting in a sparse sampling of the q-space. (a) CHARMED [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF]; (b) HYDI [START_REF] Wu | Computation of Diffusion Function Measures in q-Space Using Magnetic Resonance Hybrid Diffusion Imaging[END_REF][START_REF] Wu | Hybrid diffusion imaging[END_REF].

of the signal E allows one to analytically express the spherical Laplacian ∇ 2 b of the restriction of E to the plane orthogonal to the vector k:

q⊥k 1 q ∇ 2 b E(q)dq = - 1 2 u⊥k ∇ 2 b I(u)du, (80)
where the function I(u) is defined as a function of coefficients λ k and α k (Aganj et al., 2009b):

I(u) = N k=1 λ k (u) ln(-ln α k (u)). ( 81
)
Model estimation from the data Aganj et al. (2009b) propose the numerical estimation of the parameters λ k and α k with the confidence region algorithm [START_REF] Branch | A subspace, interior, and conjugate gradient method for large-scale bound-constrained minimization problems[END_REF] in the general case, i.e. for any number N of Gaussian functions. They also propose an analytical formulation for the special case of a mono or bi-Gaussian decay (N = {1, 2}).

Once the parameters λ k and α k are estimated, the function I defined in equation ( 81) is approximated in the basis of the spherical harmonics functions y m l :

a lm = I, y m l = u∈S 2 I(u)y m l (u)du (82)
In practice, (Aganj et al., 2009b) estimates equation ( 82) with the linear least squares method. Furthermore, the numerical computation of the l function as defined in the equation ( 81) can be instable, requiring the introduction of some regularization (Aganj et al., 2009b).

Processing and extraction of diffusion features

Exact ODF: Aganj et al. (2009b) propose a second-order ODF 2 estimation based on an multi-Gaussian interpolation:

ODF 2 (k) = ∞ p=0 P(pk)p 2 d p (83) = - 1 8π 2 q⊥k ∇ 2 E(q)dq = - 1 8π 2 q⊥k 1 q ∂ 2 ∂q 2 (qE) + 1 q 2 ∇ 2 E(q)dq. (84)
The radial Laplacian is simplified under the assumptions that the diffusion signal is equal to zero when q → ∞, and that the signal derivative has a supremum. Combining equation ( 80) and ( 83) leads to the following relationship:

ODF 2 (k) = 1 4π + 1 16π 2 FRT ∇ 2 b I(u) , (85)
where the FRT acronym stands for the Funk-Radon transform operator defined in equation ( 58). Since the spherical harmonics are the eigenfunctions of the Laplace-Beltrami operator ∇ 2 b (c.f . section 7.3.2), the combination of equation ( 82) and (85) leads to the simplification of the ODF expressions (Aganj et al., 2009b):

ODF 2 (k) = 1 4π + 1 8π L l l m=-l P l (0)a lm .y m l (k) (86)
As with the case of the QBI method based on HARDI sampling (one sphere in the q-space), the relationship described by equation ( 86) is interesting because it linearly relates the generalized ODF to the coefficients a lm of the function I.

Fusion ODF: [START_REF] Khachaturian | Boosting the sampling efficiency of q-ball imaging using multiple wavevector fusion[END_REF] propose an alternative for the computation of the ODF 0 for the special case of two sampling spheres in the q-space. They estimate the ODF independently in each sphere in the spherical Gabor wavelets domain [START_REF] Freeden | Constructive approximation on the sphere[END_REF], which results, respectively, in the set of coefficients a i and b i . They then merge these coefficients using the non-linear maxmod operator [START_REF] Burt | The Laplacian pyramid as a compact code[END_REF][START_REF] Li | Multisensor image fusion using the wavelet transform[END_REF]:

maxmod(a i , b i ) = a i if |a i | ≥ |b i | b i otherwise (87)
Advantages and limitations Aganj et al. (2009b) propose a linear and fast computation of the second-order ODF 2 in the spherical harmonics basis for a sparse sampling of the q-space.

+
+ The second-order ODF 2 does not require any postprocessing to sharpen the ODF angular profile.

+ In contrast to the approximate ODF, for which the shape depends on the radius of the HARDI sphere sampling, the generalized ODF is theoretically independent of the sparse sampling, and this eases inter-subject comparison.

+ The number of q-space samples required for the ODF estimation is moderate [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF][START_REF] Wu | Computation of Diffusion Function Measures in q-Space Using Magnetic Resonance Hybrid Diffusion Imaging[END_REF].

-The multi-Gaussian prior on the radial profile of the diffusion signal prevents the accurate approximation of a profile with a "diffusion-diffraction" pattern, as reported in [START_REF] Callaghan | Diffractionlike effects in nmr diffusion studies of fluids in porous solids[END_REF][START_REF] Cory | Measurement of translational displacement probabilities by nmr: an indicator of compartmentation[END_REF].

-The proposed method can extract only one feature of the diffusion propagator, namely the ODF.

-The fusion ODF is an ad hoc extension of the QBI method for multi-sphere sampling.

Generalized DTI

Local diffusion modeling [START_REF] Liu | Characterizing non-gaussian diffusion by using generalized diffusion tensors[END_REF] propose to extend the second-order tensor modeling of the diffusion signal by taking advantage of a sparse sampling instead of HARDI. This approach, called generalized diffusion tensor imaging (GDTI), differs from the HOT method presented in Section 3.3, for which the higher order tensors capture only the angular profile while the radial profile is assumed to be Gaussian. In the GDTI formalism, the coefficients of the fourth-order tensor or higher capture the difference with a purely Gaussian profile (captured by the second-order DTI). The diffusion signal E is expressed as:

E(q) = exp        L l=2 (i) l (2π) l ∆ - l -1 l + 1 δ D i 1 ,...,i l q i 1 q i 2 . . . q i l        . (88)
As noted by [START_REF] Liu | Characterizing non-gaussian diffusion by using generalized diffusion tensors[END_REF], the diffusion signal E can also be interpreted as a characteristic function, i.e., the Fourier transform of the water displacement PDF. Such a characteristic function is usually expressed by a cumulant-generating function Q i 1 ,...,i l of order l:

E(q) = exp        L l=1 (-i2π) l l! Q i 1 ,...,i l q i 1 . . . q i l        , (89)
where the cumulants Q i 1 ,...,i l of order l are the Taylor expansion coefficients of the characteristic function logarithm. These cumulants are proportional to the diffusion tensor of the same order D i 1 ,...,i l defined in equation ( 88), so that

Q i 1 ,...,i l = (-1) l l!D i 1 ,...,i l ∆ - l -1 l + 1 δ . (90)
We point out that for the special case of second-order tensors L = 2, as used in the DTI method (c.f . Section 2), equation (90) simplifies to the well known Einstein-Smoluchowski equation [START_REF] Einstein | On the movement of small particles suspended in stationary liquids required by the molecular-kinetic theory of heat[END_REF][START_REF] Smoluchowski | Drei vortrage uber diffusion, brownsche bewegung und koagulation von kolloidteilchen[END_REF]:

Q i 1 ,i 2 = pp T = 2D i 1 ,i 2 (∆ -δ/3). ( 91
)
Model estimation from the data As for the higher order tensor used in HARDI sampling (c.f . Section 3.3), [START_REF] Liu | Characterizing non-gaussian diffusion by using generalized diffusion tensors[END_REF] estimates the tensors D i 1 ,...,i l of order l with a matrix formulation of the n q system of equations ( 88), where n q is the number of samples of the q-space. The basis matrix M, for the n q samples set located at {q k } 1≤k≤n q in the q-space is:

M = (i) l (2π) l ∆ - l -1 l + 1 δ q k i 1 q k i 2 . . . q k i l 1≤k≤n q ,1≤l≤L
.

Liu et al. ( 2004) estimate the vector made of tensor components D i 1 ,...,i l of order l with Singular Value Decomposition (SVD) [START_REF] Stewart | Introduction to matrix computations[END_REF]. The cumulants Q i 1 ,...,i n are then computed using equation ( 90).

Processing and extraction of diffusion features

Diffusion Propagator: The estimated cumulants Q i 1 ,...,i n are analytically related to the diffusion propagator by a Gram-Charlier A series [START_REF] Kendall | The advanced theory of statistics[END_REF][START_REF] Liu | Characterizing non-gaussian diffusion by using generalized diffusion tensors[END_REF][START_REF] Mccullagh | Tensor methods in statistics[END_REF]:

P(p) = 1 2πQ i 1 ,i 2 exp - p • p 2Q i 1 ,i 2        1 + L l=3 Q i 1 ,...,i l l! H i 1 ,...,i l (p)        , (93)
where H i 1 ,...,i n (p) stands for the Grad Hermite of order n defined as [START_REF] Grad | Note on N-dimensional Hermite polynomials[END_REF]:

H i 1 ,...,i l (p) = (-1) l exp p 2 2 ∂ ∂p i 1 ∂ ∂p i 2 . . . ∂ ∂p i l exp p 2 2 . (94)
GDTI can also be interpreted as a multidimensional generalization of the simple harmonic oscillators presented in Section 4.3, which approximates the radial diffusion signal in one dimension in the basis of Hermite functions. This technique enables one to retrieve the complete diffusion propagator P without any prior assumptions on the diffusion signal E. From P, [START_REF] Liu | Characterizing non-gaussian diffusion by using generalized diffusion tensors[END_REF] compute the iso-radius P(p k), k ∈ S 2 of radius p , the skewness and the kurtosis map.

Advantages and limitations

+ The diffusion signal is estimated without prior assumptions about its radial or angular profile.

+ An analytical relationship exists between the diffusion signal approximation in the tensor basis and the diffusion propagator.

+ The number of samples in the q-space required for the estimation of the diffusion signal is moderate [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF][START_REF] Wu | Computation of Diffusion Function Measures in q-Space Using Magnetic Resonance Hybrid Diffusion Imaging[END_REF].

-The cumulant expansion of the diffusion signal is not guaranteed to converge for large values of the diffusion, and may in fact diverge in practice [START_REF] Cramér | Mathematical Methods of Statistics[END_REF][START_REF] Frøhlich | Effect of impermeable boundaries on diffusion-attenuated MR signal[END_REF][START_REF] Kendall | The advanced theory of statistics[END_REF].

Diffusion Propagator Imaging (DPI)

Local diffusion modeling The spherical harmonics functions, which are often used by the methods based on HARDI sampling, are the solution of the angular part of the Laplace equation, expressed in spherical coordinates (see Section 7.3). The solution of the full Laplace equation (i.e., with both angular and radial parts) results in the solid harmonics. Descoteaux et al. (2009b) proposed Diffusion Propagator Imaging (DPI), which approximates the diffusion signal E in the solid harmonics basis:

E(q) = ∞ l=0 l m=-l c lm q l+1 + d lm q l . y m l (u) ( 95
)
Model estimation from the data The signal approximation in the solid harmonics basis defined in equation ( 95) requires the computation of the coefficients c lm and d lm :

c lm = q∈R 3
E(q) y m l (u) q l+1 dq and d lm =

q∈R 3 E(q)q l y m l (u)dq.

(96)

For a set of n q samples in the q-space, the computation of the coefficients c lm and d lm leads to a system of n q linear equations, which can be solved using linear weighted least squares methods (Descoteaux et al., 2009b).

Processing and extraction of diffusion features

Diffusion Propagator: Once the diffusion signal is estimated in the solid harmonics basis, resulting in the coefficients c lm and d lm , it is possible to analytically compute, in the same basis, the diffusion propagator P for any arbitrary value of water displacement p (iso-radius) (Descoteaux et al., 2009b):

P(p) = 2 ∞ l=0 l m=-l (-1) l/2 2π l-1 p l-2 (2l -1)!! c lm y m l (p). ( 97
)
As mentioned in (Descoteaux et al., 2009b), the expression of P in equation ( 97) does not depend on the coefficients d lm , but these are still required for the signal estimation in the solid harmonics basis. Furthermore, the radial functions of equation ( 95) 1/q l+1 and q l , do not seem to be well adapted to the multi-Gaussian decay empirical observations of the signal.

Advantages and limitations

+ The diffusion signal is estimated without prior assumptions about its radial or angular profile.

+ An analytical relationship exists between the diffusion signal approximation in the solid harmonics basis and the diffusion propagator.

+ The number of samples in the q-space required for the estimation of the diffusion signal is moderate [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF][START_REF] Wu | Computation of Diffusion Function Measures in q-Space Using Magnetic Resonance Hybrid Diffusion Imaging[END_REF].

-It is unclear how the radial diffusion basis functions is adapted to the multi-Gaussian decay of the signal reported by empirical observations [START_REF] Niendorf | Biexponential diffusion attenuation in various states of brain tissue: implications for diffusion-weighted imaging[END_REF] or to the "diffusion-diffraction" pattern [START_REF] Callaghan | Diffractionlike effects in nmr diffusion studies of fluids in porous solids[END_REF][START_REF] Cory | Measurement of translational displacement probabilities by nmr: an indicator of compartmentation[END_REF].

Generalized DOT

Local diffusion modeling Özarslan et al. (2006) propose an extension of the DOT method for the case of multi-spherical sampling of the q-space, and generalize the original method to the multi-Gaussian case. The modeling of the diffusion signal E by n Gaussian functions is expressed as:

E(q) = n i=1 f i (u) exp -4π 2 q 2 τD i (u) , (98)
where f i stands for the weight of the i th Gaussian function while the sum satisfies n i=1 f i (u) = 1.

Model estimation from the data The generalized DOT requires an estimation of the parameters f i and D i , but the authors of (Özarslan et al., 2006) do not provide details of the method they use. It seems conceivable to use the non-linear least squares Levenberg-Marquardt technique [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF].

Once the parameters f i and D i are estimated, the computation of the original function I l of the mono-Gaussian case (c.f . equation ( 65)) is generalized to the multi-Gaussian case:

I l (u) = n i=1 f i (u)I l i (u) = n i=1 f i (u) 4π ∞ q=0 E i (q) j l (2πqp )q 2 dq . (99)
The function I l is then approximated in the spherical harmonics basis, resulting in the set of coefficients {a lm } 0≤l≤L , which is similarly to the original DOT method (see Section 3.5 for further details):

a lm = I l , y m l = u∈S 2 I l (u)y m l (u).du (100)
Processing and extraction of diffusion features The isoradius spherical feature is extracted directly from the coefficients a lm of equation ( 100), as in the case of the original DOT method presented in Section 3.5.

Advantages and limitations

+ The diffusion signal is estimated without prior assumptions about its angular profile.

+ An analytical relationship exists between the diffusion signal model and the diffusion propagator.

+ The number of samples in the q-space required for the estimation of the diffusion signal is moderate [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF][START_REF] Wu | Computation of Diffusion Function Measures in q-Space Using Magnetic Resonance Hybrid Diffusion Imaging[END_REF].

-The multi-Gaussian prior on the radial diffusion prevents the accurate approximation of the observed "diffusiondiffraction" pattern [START_REF] Callaghan | Diffractionlike effects in nmr diffusion studies of fluids in porous solids[END_REF][START_REF] Cory | Measurement of translational displacement probabilities by nmr: an indicator of compartmentation[END_REF]). Therefore it is not clear that this relatively constrained multi-Gaussian model results in accurate diffusion features.

-The proposed method can extract analytically only spherical features of the diffusion propagator, namely, the isoradius and the ODF.

-The non-linear pre-estimation of the parameters f i and D i may causes stability and accuracy issues since they depend on the initial values, which vary with data [START_REF] Aubert | Mathematical problems in image processing: partial differential equations and the calculus of variations[END_REF].

Tomographic Reconstruction

Local diffusion modeling [START_REF] Pickalov | 3d tomographic reconstruction of the average propagator from mri data[END_REF] propose to consider the Fourier transform which relates the diffusion signal E to its propagator P by the projection-slice theorem:

E(q) = p∈R 3 P(p) exp -i2πq T p dp (101) = p q p ⊥q P(p ) exp -i2πq T p dp exp -i2πq T p dp. ( 102
)
One recognizes the Radon transform of the marginal density p q of the propagator P [START_REF] Radon | Uber die Bestimmung von Funktionen durch ihre Integralwerte langs gewisser Mannigfaltigkeiten[END_REF]:

E(p) = p q p q (p) exp -i2πq T p dp. ( 103
)
The sampling of the diffusion signal along radial lines of the q-space enables one to compute the diffusion propagator P using the inverse Radon transform of p q , which is explicitly expressed as:

p q (p) = q p E(q) exp i2πq T p dq. ( 104
)
Model estimation from the data The inversion of the Radon transform as described in equation ( 104) is not numerically usable due to instability issues. This is a problem that is common to many fields beyond dMRI [START_REF] Deans | The Radon transform and some of its applications[END_REF]. Therefore, numerous methods for a stable inversion of the Radon transform have been devised and [START_REF] Pickalov | 3d tomographic reconstruction of the average propagator from mri data[END_REF] propose to use the iterative algorithm of Gerchberg and Papoulis [START_REF] Defrise | A regularized iterative algorithm for limitedangle inverse Radon transform[END_REF][START_REF] Veretennikov | Emission tomography of a micropinch discharge plasma[END_REF]).

To ensure the numerical stability of the reconstruction process, the authors introduce regularity on the signal E and on the diffusion propagator P, based on Tikhonov regularization.

Processing and extraction of diffusion features

The full diffusion propagator P is directly estimated by the inverse Fourier transform. Nonetheless, to our knowledge there is no method based on this formalism for the extraction of features of P.

Advantages and limitations

+ The diffusion signal is estimated without prior assumptions about its radial or angular profile.

+ A large literature exists on the stable inversion of the Radon transform.

-The number of q-space samples required to estimate the propagator estimation is relatively high (496 according to [START_REF] Pickalov | 3d tomographic reconstruction of the average propagator from mri data[END_REF]).

-Tikhonov regularization of the diffusion signal might result in overly smooth reconstruction.

Spherical Polar Fourier expansion (SPF)

Local diffusion modeling To take advantage of whaterver the acquisition protocol may be (generally one or several spheres in the q-space), Assemlal et al. (2009a) define an orthonormal basis from a combination of angular and radial elementary functions expressed in spherical coordinates. The authors then expand the MR signal attenuation E as the following series in the orthonormal basis Spherical Polar Fourier (SPF):

E(q) = ∞ n=0 ∞ l=0 l m=-l Ψ nlm (q) with Ψ nlm (q) = R n (q) y m l (u), (105)
so that n ∈ N is the radial index, and l ∈ N, m ∈ Z, -l ≤ m ≤ l are the angular indexes. The symbols a n,l,m are the series coefficients, y m l are the real spherical harmonics (SH), and R n is an orthonormal radial basis function made of Gaussian-Laguerre (GL) functions.

The radial component of the MR signal is reconstructed by the functions R n . Since the number of radial samples is highly limited, the radial signal is expected to be sparse in R n in the basis, i.e., to capture the radial attenuation of E with a few orders n. Assemlal et al. (2009a) propose a general estimation method for the radial part of the signal attenuation E using a normalized basis of generalized Gaussian-Laguerre polynomials R n :

R n (q) = κ n (ζ) exp - q 2 2ζ L 1/2 n q 2 ζ , (106)
κ n (ζ) = 2 ζ 3/2 n! Γ(n + 3/2) 1/2 (107) L (α) n (x) = ((α + 1) n /n!) 1F1(-n, α + 1, x), (108)
where ζ denotes the scale factor and L (α) n (x) are the generalized Laguerre polynomials which are the eigenfunctions of the Fourier-Bessel transform [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]. A low order N assumes a radial Gaussian behaviour, as in [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF][START_REF] Özarslan | Resolution of complex tissue microarchitecture using the diffusion orientation transform (dot)[END_REF]. In contrast, a high order N leads to a model-free estimation. Figure 13 illustrates the actual adequacy of the first R n functions in an experimental MR signal from erythrocytes (c.f . Figure10). The scale factor ζ can be easily estimated from experimental values of decay on data samples. Model estimation from the data The signal approximation in the SPF basis defined in equation ( 105) requires the computation of the coefficients c lm and d lm :

a nlm = q∈R 3 E(q)Ψ nlm (q)dq. ( 109
)
For a set of n q samples in the q-space, the computation of the coefficients a lm leads to a system of n q linear equations, which can be solved using linear weighted least squares. In the case where the noise is prominent, Assemlal et al. (2009a) use a variational framework to simultaneously estimate and regularize the diffusion image. The estimation process is based on a maximum a posteriori Rician probability (MAP), which takes into account the empirical nature of the noise.

Processing and extraction of diffusion features Additionally, Assemlal et al. (2009a) introduce a generic framework for the extraction of various features of the diffusion, and a fortiori of brain microstructure at a local voxel scale. Since we have a continuous representation of E from the SPF coefficients, let G(k) = P(p)h k (p)dp be a characteristic G of the PDF at point k, where h k denotes a projection function at point k. Table 1 represents several popular characteristics G which can be evaluated using this computational scheme. As a result, any characteristic G defined from equation ( 110) can be computed directly from the SPF coefficients. Indeed, since the SPF are an orthonormal basis the following relation holds:

G(k) = p∈R 3 P(p)H k (p)dp = q∈R 3 E(q)h k (q)dq = ∞ n,l,m a n,l,m b k n,l,m . (110)
where H k is the inverse three dimensional Fourier transform of h k and a nlm , b k nlm respectively denote the SPF expansion of E and H k . Therefore, the computation of a propagator feature G(k) from equation ( 110) simply turns into a very fast dot product between two vectors of SPF coefficients.

In addition, to the above framework, Cheng et al. (2010a,b) have recently demonstrated analytical solutions for the diffusion propagator (PDF) and the ODF (See Table 1).

Advantages and limitations

+ The diffusion signal is estimated without any prior assumptions about its radial or angular profile.

+ The number of samples in the q-space required for the estimation of the full diffusion signal is moderate (approximately 120 samples as reported in (Assemlal et al., 2009a)).

+ In addition to the fast linear estimation, the variational framework enables a robust estimation which takes into account the Rician distribution of the noise.

+ The SPF expansion method enables the linear extraction of various propagator features directly from the SPF coefficients. It does not require the reconstruction of the full diffusion propagator at first to extract one feature. in contrast to the DSI, GDTI, DPI and Tomographic reconstruction methods, and thus avoid unnecessary computations which may compromise the numerical accuracy.

+ Albeit not required, the analytical formula of the diffusion propagator features given in table 1 improve the accuracy of the results.

-The robust estimation is based on a slow iterative process, approximately one hour for a 128 × 128 × 64 image with 65 samples in the q-space (Assemlal et al., 2009a).

-The diffusion is a positive process [START_REF] Wedeen | Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging[END_REF], however the linear estimation gives no guarantee that the approximated diffusion signal in the SPF basis remains positive.

Features G

Projection Functions hk

Zero displacement probability 4 √ π(-1) n γ 3/4 Γ(n+3/2) n! δ l0 δ m0 Funk-Radon Transform (FRT) q 2 R n (q ) 2πP l (0)y m l (k) 0 th order Orientation Density Function (ODF 0 ) 2πζ Z n j=0 κ n (ζ) j-1/2 j (-1) n-j P l (0)y m l (k) 2 nd order Orientation Density Function (ODF 2 ) 4π(-1) n γ 3/4 Γ(n+3/2)

n! P l (0)y m l (k) Diffusion Propagator (PDF) 4 (-1) l/2 ζ (l+3)/2 π l+3/2 k l Γ(l+3/2) κ n (ζ) n j=0 (-1) j j! n+1/2 n-j 2 l/2+ j-1/2 Γ(l/2 + j + 3/2)1F1 2 j+l+3 2 ; l + 3 2 ; -2π 2 k 2 ζ y m l (k)
Table 1: Some features G of the diffusion propagator P which can be computed in the SPF framework [START_REF] Assemlal | Traitement et analyse d'images IRM de diffusion pour l'estimation de l'architecture locale des tissus[END_REF]Assemlal et al., 2009a), with their respective projection functions H k .

Conclusion

Despite its many successes, diffusion MRI has one major weakness -the acquisition time for a diffusion image sequence of a brain at full resolution is approximately one hour. During this time the patient must lie motionless, which is infeasible for clinical applications. The methods described in this review of the literature suggest a possible way to overcome this constraint, namely, by sparse sampling of q-space combined with angular and radial reconstruction. The existing methods in the literature can be classified into three groups, based on the nature of sampling used:

• angular sampling only: These can be used for ODF reconstruction, for finding isoprobability surfaces, and for fiber tracking in brain white matter;

• radial sampling only: These can be used for measurement of the apparent kurtosis coefficient and of the "diffusiondiffraction" phenomena; and

• combined angular and radial sampling: These can be used for improved ODF reconstruction and for inferring the full diffusion propagator (PDF).

For each manner of q-space sampling, Table 2 summarizes the hardware requirements and the main features that have been proposed in the literature. Due to the limited number of samples of the diffusion signal in practice, a direct processing of the input data is typically not reliable. Therefore numerous reconstruction models have been proposed in the literature, of which a large set have been described in this review. These methods can be grouped into two categories:

1. Those for which the estimation of the signal with parametric models is based on strong prior assumptions. Such models can drastically reduce the number of required samples, provided their complexity remains manageable (i.e., only a few parameters characterize the model). Nonetheless, a simple model may not be adequate particularly in situations where the cause and nature of the diffusion process is under investigation; and 2. Those for which the reconstruction of the diffusion signal is carried out without prior assumptions, by factorizing using a infinite number of basis functions. Mathematically, such methods ensure a perfect reconstruction of the signal. In practice, however, the low number of samples of the diffusion signal leads to the constraint that the infinite summation must be truncated. As a consequence, these methods impose a weak, implicit prior which depends on the particular choice of basis functions used.

Looking at the issue of clinical feasibility, the DTI method is the most popular since it has several advantages. It is fast and easy to implement and leads to an acceptable and reproducible approximation when the underlying microstructure is relatively coherent, i.e., oriented in one single direction.

Although the error of the DTI method in modeling the diffusion signal at a voxel scale is modest, it can be amplified by iterative fiber-tracking algorithms, resulting in non-reproducible neural pathways. In addition, DTI can not detect subtle modifications in the diffusion pattern induced by various neuronal diseases. Therefore it is already apparent that there is a strong need to move from DTI to HARDI acquisitions. Yet, which modeling method to choose among those that are presently clinically feasible is non-trivial and the answer depends on several constraints including the actual application, the scan time and the scanner performance. The Q-Ball Imaging technique (section 3.4) was the first method to take advantage of HARDI acquisition and as such is still widely successful, since it exhibits relatively little variation in the results when compared with other HARDI modeling methods.

At first sight, it might seem that the development of highly precise modeling methods of the diffusion signal is inappropriate in contrast to DTI given the constant reduction of voxel size in the last two decades. Indeed, one might envision voxel size converging to the dimension of an axon in the future, thus eliminating the need for higher order models altogether. In reality, such voxel size reduction will not be feasible in the current state since it would dramatically decrease the SNR as fewer and fewer water molecules would generate the diffusion signal. Introducing higher magnetic field strength in MRI scanners while keeping within the limits that can be tolerated by human subjects, might compensate for this loss in SNR. The main q-space sampling methods in dMRI, their respective hardware requirements and the features which can be extracted using each.

Appendix

In this section, we review the mathematical notation and the tools required for understanding the estimation methods used in DMRI analysis and reconstruction.

Spherical Coordinates

Recall that q = (2π)γδg, this suggests that the spherical nature of q-space sampling arises from the gradient hardware limitations of MRI scanners. Whereas there are virtually no limit to the number of orientations q/||q|| of the q-space one can acquire, the radial sampling using different diffusion wave-vector norm ||q|| pushes the gradient hardware of MRI scanners to its limits. This leads naturally to the processing of the acquired data in a spherical coordinates system. We define the radius as q ∈ [0, ∞[, the colatitude angle θ ∈ [0, π] and the longitudinal angle φ ∈ [0, 2π[ in accordance with the conventions of physics. Fig. 14 illustrates the definition of spherical coordinates for a point q ∈ R 3 : q = q • u, with q = ||q|| and u = q ||q|| ∈ S 2 ..

The Cartesian coordinates are expressed as:

Figure 14: An illustration of spherical coordinates. The point q ∈ [0, ∞) is defined by its radius q, and its two angles θ, φ.

x = q cos(φ) sin(θ) q = x 2 + y 2 + z 2 (112)

y = q sin(φ) sin(θ) θ = arccos(z/q) (113) z = q cos(θ) φ = arctan(z/q) (114)
The infinitesimal differential elements of surface area du and volume dq are defined as:

du = sin(θ) dθ dφ, (115)
dq = q 2 sin(θ) dq dθ dφ.

In the following, we review two mathematical tools that are widely used in DMRI, since they ease the modeling of the diffusion image from experimental observations: tensors and spherical harmonics.

Tensors

Definition

A tensor is a mathematical object which generalizes the notion of scalars (order 0 tensors), vectors (order 1 tensors) and matrices (order 2 tensors). A tensor does not depend on a reference system, indeed the tensor object T should not confused with its components T i 1 ...i p j 1 ... j q which do depend on the choice of a specific coordinate system. The tensor object itself T, however, remains invariant under changes in the coordinate system. der 0 order 1 order 2 order 3

Figure 15: Tensors components of increasing order. Left to Right: a scalar, a vector, a matrix and a tensor of order 3.

Let V be a finite dimensional vector space. A tensor T of type (m, n) and of order m + n is defined by a multilinear application T (m,n) (V) → R [START_REF] Abraham | Manifolds, tensor analysis, and applications[END_REF][START_REF] Borisenko | Vector and Tensor Analysis with Applications[END_REF]:

T (m,n) (V) = V ⊗ . . . ⊗ V m ⊗ V * ⊗ . . . ⊗ V * n (117)
where ⊗ stands for the tensor product and V * is the symbol for the dual of V, defined as the vector space of the linear functional on V. Let {e 1 , . . . , e m } ∈ V and {e 1 , . . . , e n } ∈ V * be the respective basis in V and V * . The components of the tensor T ∈ T (m,n) (V) are then expressed as: T = T i 1 ...i p j 1 ... j n e i 1 ⊗ . . . ⊗ e i m ⊗ e j 1 ⊗ . . . e j n .

(118)

The notation for tensor components T i 1 ...i m j 1 ... j n is quite similar to that for matrices (T = i j )). A matrix is a tensor of type (1, 1) so that T ∈ T (1,1) and is written as T = T j i (e i ⊗ e j ) in tensor notation. However, unlike a matrix which has two indexes, a tensor has m + n indices. We refer the reader to [START_REF] Abraham | Manifolds, tensor analysis, and applications[END_REF][START_REF] Borisenko | Vector and Tensor Analysis with Applications[END_REF]. for a more precise definition of tensors and their properties.

The Second Order Tensor

In DMRI, the diffusion tensor refers to the generalization of the apparent diffusion coefficient D to a second-order tensor D, which is symmetric and with positive entries, i.e. ∀x ∈ R 3 , x T Dx ≥ 0. Under this assumption, the components D XYZ of the tensor in an XYZ frame of reference are expressed as a 3 × 3 symmetric and positive-definite matrix: The orientation of this tensor is often expressed by using a basis consisting of eigenvectors e = (e 1 , e 2 , e 3 ), so that the tensor components form a diagonal matrix D λ whose elements λ i are the eigenvalues. In this setting, the tensor D is often visualized as an ellipsoid with its major axis directions and lengths given by the eigenvectors and eigenvalues (see Fig. 16). As illustrated

D XYZ =          
D Λ =          λ 1 0 0 0 λ 2 0 0 0 λ 3          

A Tensor Field

A tensor field is a tensor valued function defined on a topological space. In other words, each point of the field is associated with a tensor T. In the same fashion that a tensor object can be represented in several ways (using an ellipsoid, a parallelepiped, etc.), there are several graphical representations of tensor fields, as shown in Fig. 18. Fig. 19 gives examples of first-order and second-order tensor fields. We refer the reader to [START_REF] Tschumperlé | Vector-valued image regularization with PDEs: A common framework for different applications[END_REF][START_REF] Weickert | Visualization and Image Processing of Tensor Fields[END_REF] for an introduction to tensor fields and their applications. The modeling of the diffusion image by a tensor field was proposed in the classical method of diffusion tensor imaging (DTI) [START_REF] Basser | Estimation of the effective selfdiffusion tensor from the nmr spin echo[END_REF], where diffusion in the q-space is modeled by a tensor of order two. We now discuss another model that is adapted to a spherical sampling of q-space, namely, spherical harmonics.

The Spherical Harmonics

There are many functions to approximate a function on the sphere. Among these possibilities spherical harmonics (SH) are particularly suited to the processing of diffusion images, because they have the following advantages:

1. This basis is adapted to a spherical sampling of the qspace commonly called high angular resolution imaging (HARDI) [START_REF] Tuch | High angular resolution diffusion imaging of the human brain[END_REF]. 2. The spherical harmonic (SH) representation can express most signals in spherical media that are sufficiently compact and smooth, with a small number of coefficients. 3. The SH are the eigenfunctions of the Laplace-Beltrami erator, which allows one to sharpen or smooth the spherical shape by exploiting a simple linear relationship. 4. The rotation of the SH is expressed by a simple formula, which offers particular advantages in the SPF framework (as discussed in section 5.8) .

We refer the reader to [START_REF] Arfken | Mathematical methods for physicists[END_REF][START_REF] Healy | Spherical deconvolution[END_REF]) for a comprehensive overview of the properties of spherical harmonics. In the following we present the definition of complex spherical harmonics, followed by a discussion of the subset of real symmetric spherical harmonics, with properties that are suitable for diffusion images.

Definitions

A function f is called harmonic when it satisfies the Laplace equation, whose expression in Cartesian coordinates (x, y, z) is given by:

∇ 2 f = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 + ∂ 2 ∂z 2 f = 0. ( 120
)
The Eq. ( 120) is rewritten in spherical coordinates (q, θ, φ) as:

1 q 2 ∂ ∂q q 2 ∂ ∂r + 1 sin θ ∂ ∂θ sin θ ∂ ∂θ + 1 sin 2 θ ∂ 2 ∂φ 2 f = 0 (121)
The method of separation of variables q, θ and φ is used to solve the spherical Laplace Eq. ( 121), so we assume that f is written as [START_REF] Courant | Methods of mathematical physics[END_REF]:

f (q, θ, φ) = q l Θ(θ)Φ(φ).

We are interested in solutions to the angular part of the Laplace equations ( 121). The combination of Eq. ( 121) and (122) gives:

1 sin θ ∂ ∂θ sin θ ∂ ∂θ + 1 sin 2 θ ∂ 2 ∂φ 2 + l(l + 1) f = 0. ( 123
)
Finally, the solutions of Eq. ( 123) are called spherical harmonics and are written for l ∈ N and m ∈ Z, m ∈ [-l, l] ( [START_REF] Arfken | Mathematical methods for physicists[END_REF] as:

Y m l (θ, φ) = 2l + 1 4π (l -m)! (l + m)! P m l (cos θ)e imφ , (124)
with P m l the generalized Legendre polynomials defined on the interval [-1, 1]: The complex spherical harmonics Y m l from Eq. ( 124) form a basis of orthonormal functions: with (z) the complex conjugate and δ j j the Kronecker symbol, which is equal to one if j = j and is zero otherwise.

P m l (x) =              1 2 l l! d l
To simplify the notation, we introduce the following equivalence relation:

Y m l (u) = Y m l (θ, φ).
(127)

The Laplacian Operator on the Sphere

The Laplace-Beltrami operator is denoted by ∆ b and generalizes the Laplace operator for functions defined on surfaces and in particular on spheres. This operator is defined by:

∆ b = 1 sin θ ∂ ∂θ sin θ ∂ ∂θ + 1 sin 2 θ ∂ 2 ∂φ 2 . ( 128
)
Since the Laplace-Beltrami operator is involved in the angular Laplace equation given by Eq. ( 123), the spherical harmonics satisfy the relation:

∆ b Y m l = -l(l + 1)Y m l . ( 129
)
This relationship is used in several approaches to the regularization and enhancement of spherical functions (Assemlal et al., 2009a;[START_REF] Descoteaux | Regularized, fast and robust analytical q-ball imaging[END_REF].

A Real and Symmetric Subset

The number of samples of the observed signal determines the number of coefficients in the reconstruction by a sum of basis functions, here the SH. For a fixed number of coefficients, one should not "waste" these coefficients for the estimation of impossible or uninteresting shapes, but rather should use them wisely to gain accuracy. For this, two prior assumptions are commonly used in DMRI, which restrict the definition domain of spherical harmonics to the ideal diffusion signal [START_REF] Alexander | Detection and modeling of nongaussian apparent diffusion coefficient profiles in human brain data[END_REF][START_REF] Frank | Characterization of anisotropy in high angular resolution diffusion-weighted mri[END_REF]: (i) The signal is antipodally symmetric on the unit sphere, and (ii) It is real and positive. Many approaches in the literature assume that the diffusion at a voxel-scale is symmetric on the sphere. For this, the SH has an interesting property when subjected to a transformation T :

T : (θ, φ) → (π -θ, π + φ) Y m l (T (θ, φ)) = Y m l (θ, φ) if l ∈ 2N -Y m l (θ, φ) if l 2N (130)
Eq. ( 130) implies that when the l is even, the spherical function is symmetric on the unit sphere. Also, a symmetric function on the sphere can be expressed as a series of SH with even order l.

In addition, although in theory the diffusion MRI signal E is captured in the complex domain C, in practice the use of the real part of E allows one to overcome the sensitivity due to phase

Figure 1 :

 1 Figure 1: Our mathematical definition of the diffusion MRI image as expressed in the equation (1). Left: the position x-space Ω x of the image, which results as the Fourier transform of the acquired k-space. Right: the diffusion q-space Ω q for one voxel of the x-space (depicted here as a yellow rectangle).

Figure 2 :

 2 Figure 2: Acquisition of a slice of the MRI image. (a) Acquisition is done in the k-space on a slice-by-slice basis. (b) Module of the Fourier transform of the image (a).

Figure 4 :

 4 Figure 4: Simulated random walk as measured by the PGSE sequence. (a) Ideal case (δ ≈ 0), (b) Real case (δ 0). Adapted from[START_REF] Hagmann | From diffusion mri to brain connectomics[END_REF].

Figure 5 :

 5 Figure 5: Examples of local diffusion profiles E observed in the brain matter measured by dMRI. The data are represented here as volumetric images 64 × 64 × 64, where the center of the q-space q = 0 is the center of each image. (a) Free isotropic Gaussian diffusion. (b) Restricted diffusion due to the presence of a single fiber bundle. (c) Restricted diffusion in the presence of two fiber bundles in a crossing configuration. (d) Restricted diffusion which is isotropic in direction, but has a multi-Gaussian profile radially.

Figure 6 :

 6 Figure 6: The analysis of the diffusion signal is closely related to the sampling of the q-space. (a) Full sampling of the q-space is currently impractical in vivo due to the significant acquisition time it would imply. (b) Low angular resolution sampling used in DTI. (c) High angular resolution sampling (HARDI). (d) Radial only sampling used in diffusion NMR. (e) Sparse sampling which combines radial and angular measurements.

Figure 7 :

 7 Figure 7: DTI visualization using ellipsoids: from a global to a local vue. a) Visualization of the corpus callosum by a second order tensor field. b) Zoom on the bottom part of the corpus callosum. The privileged diffusion direction is the left-right axis, this is coherent since the corpus callosum links both brain hemispheres. The colours indicate the orientation of the principal eigenvector of each tensor (more details in figure 8).

Figure 8 :

 8 Figure 8: Some vectorial features based on the principal eigenvector of the DTI (second order tensor field). (a) The colour encodes the orientation of the vector and its norm is encoded by the brightness (red for left-right, blue for superior-anterior, green for anterior-posterior).(b) The curvature encodes the direction of each fiber and the brightness the smoothness strength. (c) Fibertracking: A visualization of nerve fiber bundles as three dimensional curves overlaid on the diffusion ellipsoids. These images are adapted from[START_REF] Tschumperlé | Pde's based regularization of multivalued images and applications[END_REF].

Figure 9 :

 9 Figure 9: HARDI sampling of the diffusion signal: a schematic overview of the modeling of an intravoxel crossing of bundles nerve fibers. (a) Two fiber bundles crossing at 90 • inside a voxel generate a diffusion signal. (b) ADC modeling using a DTI reconstruction.(c) Generalization of the ADC modeling with HARDI sampling. The maxima of the angular profile do not match the underlying fiber bundle directions given in (a). (d) Angular feature estimation of the diffusion propagator with HARDI sampling. Adapted from[START_REF] Descoteaux | High angular resolution diffusion mri: from local estimation to segmentation and tractography[END_REF].

  ) is fast and precise since it analytically links the EAP isoradius to the diffusion signal coefficients in the spherical harmonics basis. More recently, Aganj et al. (2009a); Canales-Rodríguez et al. (2010); Tristan-Vega et al. (

  (a) First row: The raw diffusion data E[q] visualized as consecutive 2D planes of the q-space. Second row: The iso-contours in two dimensions. Right: An iso-surface in three dimensions.(b) First row: The diffusion propagator P[p] reconstructed with the Fourier transform of E and visualized as consecutive 2D planes of the q-space. Second row: Iso-contours in two dimensions. Right: An iso-surface in three dimensions.SequenceImage size Voxel size (mm 3 ) TE/TR (ms) q-space samples Single-shot, single slice 64 × 64 × 1 3.6 × 3.6 × 3.6 156/3000 512 ∆/δ (ms) g max (mT/m) b max (s/mm 2 ) p min (µm) FOV (µm) Time (min)

Figure 13 :

 13 Figure13: The first few basis functions of the radial orthonormal basis R n , with scaling factor ζ = 100. Low order n functions exhibit Gaussian attenuation behaviour whereas higher order n are helpful to capture oscillating components of the MR signal.

  D xx D xy D xz D xy D yy D yz . D xz D yz D zz

Figure 16 :Figure 17 :

 1617 Figure16: The second-order diffusion tensor represented by an ellipsoid whose components are defined by a 3×3 symmetric positive-definite matrix. When the diffusion tensor is expressed in the basis of eigenvectors e i and eigenvalues λ i , the matrix components D Λ are diagonal.

Figure 18 :Figure 19 :

 1819 Figure 18: Some possible representations of a first-order tensor field: (a) Representation by a set of icons (arrows). (b) An image whose intensity varies with the strength and direction of smoothing of a random image using line integral convolution (LIC)[START_REF] Cabral | Imaging vector fields using line integral convolution[END_REF][START_REF] Stalling | Fast and resolution independent line integral convolution[END_REF][START_REF] Tschumperlé | Vector-valued image regularization with PDEs: A common framework for different applications[END_REF].

  dx l x 2 -1 l si m = 0 (-1) m (1x 2 ) m/2 d m dx m P l (x) si m > 0 (-1) |m| (l-|m|)! (l+|m|)! P |m| l (x) si m < 0 (125)

  φ)Y m l (θ, φ) sin(θ)dθdφ = δ ll δ mm (126)

  Parameters of the acquisition sequence of the CHARMED method, proposed by[START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF].

	Sequence	Image size	Voxel size (mm 3 ) TE/TR (ms) q-space samples
	DW-EPI	64 × 64 × 10	3 × 3 × 3		133/2700	169
	∆/δ (ms) g max (mT/m) b max (s/mm 2 ) p min (µm) FOV (µm) Time (min)
	53/47	34	10 000	10.8	-	∼ 17
	(a) Sequence	Image size	Voxel size (mm 3 ) TE/TR (ms) q-space samples
	SS-SE-EPI 128 × 128 × 30	2 × 2 × 3		122/11700	102
	∆/δ (ms) g max (mT/m) b max (s/mm 2 ) p min (µm) FOV (µm) Time (min)
	56/45	40	9375	13.2	65.8	∼ 30

Table 2 :

 2

artifacts, which can be generally regarded as noise [START_REF] Lebihan | Artifacts and pitfalls in diffusion MRI[END_REF][START_REF] Tuch | Diffusion mri of complex tissue structure[END_REF].

Using these two prior assumptions we can define a new basis of real and symmetric SH [START_REF] Anderson | Measurement of fiber orientation distributions using high angular resolution diffusion imaging[END_REF][START_REF] Descoteaux | Regularized, fast and robust analytical q-ball imaging[END_REF][START_REF] Frank | Characterization of anisotropy in high angular resolution diffusion-weighted mri[END_REF][START_REF] Hess | Q-ball reconstruction of multimodal fiber orientations using the spherical harmonic basis[END_REF]: Fig. 20 illustrates the SH y m l for the first few orders of l, which capture the low frequencies on the unit sphere. In the following, we introduce the equivalent notation to facilitate the presentation:

The real and symmetric spherical harmonics form an orthonormal basis on the unit sphere [START_REF] Blanco | Evaluation of the rotation matrices in the basis of real spherical harmonics[END_REF]:

We have seen that the basis of spherical harmonics {y j } breaks down most spherical continuous functions into a Fourier series expansion on the sphere.

Rotation of Spherical Harmonics

We now review the mathematical tools that are necessary for understanding the rotation of real spherical harmonics, namely the Euler angles and the Wigner rotation matrix.

Euler Angles The Euler angles have been developed to describe the rotation of a rigid object in three-dimensional Euclidean space R 3 . The specific orientation of an object is obtained by a sequence of three successive rotations given by the Euler angles, as illustrated in Fig. 21: 1. the rotation by angle α around the axis z which transforms the frame of reference from xyz to x y z;

2. the rotation by angle β around the axis y which transforms the frame of reference from x y z to x y z ; 3. the rotation by angle γ around the axis z which transforms the frame of reference x y z to x y z

Remark. There are twelve variants of the rotation scheme of Euler angles, which differ in the axis of rotation chosen from each of the three previous steps. In this paper, we adopt, without loss of generality, the configuration zyz [START_REF] Blanco | Evaluation of the rotation matrices in the basis of real spherical harmonics[END_REF]. Wigner rotation matrix When subject to a rotation expressed in terms of Euler angles (α, β, γ), complex spherical harmonics Y m l of order l can be expressed as a linear combination of spherical harmonics of even order Y ( l m ):

with the d-small Wigner matrix defined as [START_REF] Su | Rotation of real spherical harmonics[END_REF]:

We are interested in the real Wigner rotation matrix

We saw in Eq. ( 131) that real spherical harmonics are defined by a linear combination of complex spherical harmonics. Then the d-small Wigner matrix of Eq. ( 136), for real spherical harmonics is expressed as [START_REF] Ritchie | Parametric protein shape recognition[END_REF] Appendix B.4: