Haz-Edine Assemlal

David Tschumperlé

Luc Brun Greyc

Robust Variational Estimation of PDF functions from Diffusion MR Signal

We address the problem of robust estimation of tissue microstructure from Diffusion Magnetic Resonance Imaging (dMRI). On one hand, recent hardware improvements enable the acquisition of more detailed images, on the other hand, this comes along with a low Signal to Noise (SNR) ratio. In such a context, the approximation of the Rician acquisition noise as Gaussian is not accurate. We propose to estimate the volume of PDF-based characteristics from data samples by minimizing a nonlinear energy functional which considers Rician MR acquisition noise as well as additional spatial regularity constraints. This approach relies on the approximation of the MR signal by a series expansion based on Spherical Harmonics and Laguerre-Gaussian functions. Results are presented to depict the performance of this PDE-based approach on synthetic data and human brain data sets respectively.

Introduction

Water molecules exhibit Brownian motion which might be constraint by internal micro-structure of the brain white matter. Diffusion-Weighted Imaging (DWI) measures this local displacement using the pulse gradient spin echo sequence [START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF] in each voxel and thus provides images of the architecture of the brain. These images provide valuable information to diagnose early stages of stroke, brain diseases or neurological disorders [START_REF] Lebihan | Mr imaging of intravoxel incoherent motions: Application to diffusion and perfusion in neurologic disorders[END_REF]. However, this molecular displacement is not directly measured. Indeed, when the diffusion gradient pulse duration δ is negligible compared to diffusion time ∆, the MR signal E defined in Q-Space is related to the average displacement probability P by the Fourier transform [START_REF] Callaghan | Principles of Nuclear Magnetic Resonance Microscopy[END_REF] P (p) = q∈R 3 E(q) exp(-2πiqp T )dq, with E(q) = S(q) S 0 ,

where p is the displacement vector and q stands for the diffusion wave-vector of the Q-Space. The symbols S(q) and S 0 respectively denote the diffusion signal at gradient q and the baseline image at q = 0. Eq.( 1) naturally suggests to use the Fourier transform to numerically estimate the PDF (Probability Density Function). This technique known as Diffusion Spectrum Imaging (DSI) [START_REF] Wedeen | Mapping fiber orientation spectra in cerebral white matter with fourier transform diffusion mri[END_REF] is not clinically feasible due to the huge acquisition time required to retrieve the whole Q-Space coefficients. As a result of DSI ⋆ We are thankful to Cyceron for providing data and the fruitful technical discussions.

constraints, High Angular Resolution Diffusion Imaging (HARDI) [START_REF] Tuch | High angular resolution diffusion imaging of the human brain[END_REF] comes as an interesting alternative and suggests to sample the signal on a single sphere of the Q-Space.

Prior methods of the literature based on HARD images [START_REF] Tuch | Q-ball imaging[END_REF][START_REF] Yablonskiy | Statistical model for diffusion attenuated mr signal[END_REF][START_REF] Jian | A novel tensor distribution model for the diffusion-weighted mr signal[END_REF][START_REF] Liu | Characterizing non-gaussian diffusion by using generalized diffusion tensors[END_REF][START_REF] Özarslan | Resolution of complex tissue microarchitecture using the diffusion orientation transform (dot)[END_REF] use a single shell acquisition and have thus to assume strong priors on the radial behavior of the signal, i.e. mono-exponential decay. Nonetheless, sampling schemes on several spheres in the Q-Space have been proposed recently [START_REF] Özarslan | Resolution of complex tissue microarchitecture using the diffusion orientation transform (dot)[END_REF][START_REF] Cercignani | Optimal acquisition schemes for in vivo quantitative magnetization transfer mri[END_REF][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][START_REF] Pickalov | 3d tomographic reconstruction of the average propagator from mri data[END_REF]. Since the number of samples still remains too low to allow a precise Fourier transform, some methods rather consider computed tomography methods [START_REF] Pickalov | 3d tomographic reconstruction of the average propagator from mri data[END_REF] or an approximation of the MR signal radial attenuation by a multi-exponential function [START_REF] Özarslan | Resolution of complex tissue microarchitecture using the diffusion orientation transform (dot)[END_REF][START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF]. Note that these methods use a larger set of data but are still based on a priori models of the radial behavior of the input signal.

Nonetheless, a recent method [START_REF] Assemlal | Efficient computation of pdf-based characteristics from diffusion mr signal[END_REF] tackles this problem with a continuous representation of data from multiple shells and a fast method for computation of functions of the PDF. It involves a damped least square estimation of the best-fitting coefficients in the Spherical Polar Fourier basis. However magnitude MR data are corrupted by a Rician noise not a Gaussian one, and consequently introduces a bias in intensity at low Signal-to-Noise Ratio (SNR) which reduces the tissue contrast. This arises from complex Gaussian noise in the original frequency domain measurements named K-Space [START_REF] Granlund | Signal Processing for Computer Vision[END_REF][START_REF] Gudbjartsson | The rician distribution of noisy mri data[END_REF][START_REF] Henkelman | Measurement of signal intensities in the presence of noise in mr images[END_REF][START_REF] Sijbers | Estimation of the noise in magnitude mr images[END_REF]. In this study we propose to extend the previous approach to robustness to Rician noise within a variational framework.

In section 2, we overview the mathematical background of the method introduced in [START_REF] Assemlal | Efficient computation of pdf-based characteristics from diffusion mr signal[END_REF]. Then, we present an variational framework extension of previous method for a robust estimation in section 3. Section 4 shows validation results on both numerical and real human data-sets. Finally, we draw conclusions of the proposed approach in section 5.

Mathematical background

To be as self-contained as possible, we briefly overview the method introduced in [START_REF] Assemlal | Efficient computation of pdf-based characteristics from diffusion mr signal[END_REF] using the Spherical Polar Fourier (SPF) expansion. In order to be able to reconstruct the PDF from Eq.( 1) even with few samples, we seek a basis in which the acquired signal is sparse. For convenience, a list of common notations used in this paper is reminded in Table 1.

Spherical Polar Fourier Expansion

Let E be the MR signal attenuation, we propose to express it as a serie in a spherical orthonormal basis named Spherical Polar Fourier (SPF) [START_REF] Ritchie | High-order analytic translation matrix elements for real-space sixdimensional polar fourier correlations[END_REF]

E(q) = S(q) S(0) = ∞ n=0 ∞ l=0 l m=-l a nlm R n (||q||)y m l q ||q|| , (2)
where a nlm are the expansion coefficients, y m l are real Spherical Harmonics (SH), and R n is an orthonormal radial basis function. The angular part of the signal E is reconstructed by the complex SH Y m l which form an orthonormal basis for functions defined on the single sphere. For this reason, they have been widely used in diffusion MRI [START_REF] Alexander | Detection and modeling of non-gaussian apparent diffusion coefficient profiles in human brain data[END_REF][START_REF] Frank | Characterization of anisotropy in high angular resolution diffusionweighted mri[END_REF]. Indeed, as the diffusion signal exhibits real and symmetry properties, the use of a subset of the complex basis made of real and symmetric SH y m l strenghten the robustness of the reconstruction to signal noise and reduces the number of required coefficients [START_REF] Alexander | Detection and modeling of non-gaussian apparent diffusion coefficient profiles in human brain data[END_REF][START_REF] Frank | Characterization of anisotropy in high angular resolution diffusionweighted mri[END_REF].

The radial part of the signal E is reconstructed by the elementary radial functions R n . A sparse representation of the radial signal should approximate it in a few radial order N . Based on these observations, we propose to estimate E using the normalized generalized Gaussian-Laguerre polynomials basis R n [START_REF] Ritchie | High-order analytic translation matrix elements for real-space sixdimensional polar fourier correlations[END_REF][START_REF] Biedenharn | Angular momentum in quantum physics[END_REF]:

R n (||q||) = 2 γ 3/2 n! Γ (n + 3/2) 1/2 exp - ||q|| 2 2γ L 1/2 n ||q|| 2 γ , (3)
where γ denotes the scale factor and L (α)

n are the generalized Laguerre polynomials. The Gaussian decay arises from the normalization of the Laguerre polynomials in spherical coordinates.

The SPF forms an orthonormal basis on Spherical Harmonics (SH) and Gaussian-Laguerre polynomials. Consequently a low order truncation assumes a radial Gaussian behavior as in [START_REF] Özarslan | Resolution of complex tissue microarchitecture using the diffusion orientation transform (dot)[END_REF][START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF] and a high order truncation provides a model-free estimation. Besides, the square error between a function and its expansion in SPF to order n <= N and l <= L converges to zero as N and L go to infinity.

Fast Computation of Characteristics on the PDF

As 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 2 represents several popular characteristics G which can be evaluated using this computational scheme. A natural way to retrieve G would be to reconstruct E from the SPF coefficients, compute a Fast Fourier Transform (FFT) and finally calculate G on the PDF Table 2. A non-exhaustive list of some PDF characteristics G and their projection function h k at point k. FRT stands for the Funk-Radon Transform used in QBI, where J0 is the Bessel function of the first kind and ||q ′ || is the radius of the q-ball shell. ISO stands for isoprobability profiles. SD and FD respectively stands for slow and fast diffusion, where ||p ′ || is the radius limit between intra and extra cellular diffusion.

G ODF FRT ISO SD FD h k (p) δ(1 -|p•k| ||p||||k|| ) J0(2π||q ′ ||||p||) δ(k) ODF if ||p|| < ||p ′ || 0 if ||p|| > ||p ′ || 0 if ||p|| < ||p ′ || ODF if ||p|| > ||p ′ ||
volume; however such a scheme would induce cumbersome computations and raise numerical accuracy issues. So, any characteristic G defined from Eq.( 4) can alternatively 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 = ∞ nlm a nlm b k nlm (4)
where H k is the inverse Fourier transform of h k and a nlm , b k nlm respectively denote the SPF expansion of E and H k . Therefore, the numerical computation of G(k) which is an integration over an entire volume simply turns into a very fast dot product between two vectors of SPF coefficients.

We seek the SPF coefficients that represent the best the MR data samples E. However, E is strongly corrupted by noise and may lead to distortion of computed characteristics on the PDF.

Robustness to noise

Since the acquisition noise on the MR signal is not Gaussian, a least square fit is definitely not the best choice for such an estimation process. This issue arises especially when dealing with low SNR data as this is the case for very high-q values. Furthermore, independent voxel estimation does not reflect the spatial regularity of the diffusion function. We propose to tackle these issues with a variational framework which is adaptable to noise distribution and is able to use valuable information given by the neighbour voxels.

Variational Framework

The key idea is to estimate and regularize the whole volume of voxels at the same time. Indeed, it enables to take into account correlation between all parts of the processing pipeline instead of doing the different parts separately. Let E be the acquired dMRI volume corrupted by Rician noise, we seek the SPF coefficients A of the filtered dMRI volume Ê = M A, where the symbol M = (R n (||q j ||)y m l ( qj ||qj || )) nlm×j∈N 3 ×N denotes the SPF basis matrix. We propose to robustly estimate and regularize the SPF coefficients field from the dataset volume simultaneously by minimizing the following nonlinear functional energy:

min A Ω E ns k ψ( Êk ) + α r ϕ(||∇A||)dΩ E (5)
where Ω E ⊂ R 3 is the domain of datasets voxels. The likelihood term ψ(E k ) measures the dissimilitudes at voxel x ∈ Ω E between E and its reconstruction Ê at gradient direction k, ψ : R → R + and ϕ : R → R + are real and positive functions, α r ∈ R is the regularization weight and ||∇A|| the gradient norm defined as

||∇A|| = nlm ||∇A nlm || (6)
Note that if ψ(s) = s 2 and α r = 0 in Eq.( 5), we minimize the least square criterion. As the minimization cannot be computed straightforwardly, the gradient descent coming from the Euler-Lagrange derivation of Eq.5 leads to a set of multi-valued Partial Differential Equation (PDE) as described in Eq. [START_REF] Yablonskiy | Statistical model for diffusion attenuated mr signal[END_REF]. In practice, we first set A (t=0) to U 0 , an initial estimate of SPF coefficients. In order to estimate a solution, SPF coefficients velocity ∂A ∂t giving the direction from the current A to a solution is computed. The latter is done several times until convergence (typically when

ε ∈ R + , ε → 0, ∂A ∂t < ε),    A t=0 = U 0 ∂Aj ∂t = ns k M k,j ψ ′ ( Êk ) + α r div(ϕ(||∇A||)) (7)
The initial estimate U 0 is computed either by considering a random field or a more structured one. A good choice is to start from an initial set which is not so far from the global minimum; so the linear least square estimation seems to be an adequate alternative. Indeed, least square minimization is the global minimum when ψ(s) = s 2 and α r = 0. One can expect the minimum to be close enough to the least square minimum through variations of ψ and ϕ; and should consequently bring down the number of iterations required to converge.

Likelihood function ψ

The diffusion MR magnitude images are corrupted by noise and the best ψ function is the one specific to the MR scanners, that is to say the Rice distribution whose probability density function is:

p(E| Ê, σ) = E σ 2 exp -(E 2 + Ê2 ) 2σ 2 I 0 E • Ê σ 2 ( 8
)
where σ is the standard deviation of the noise and I 0 is the modified zerothorder Bessel function of the first kind. We adapt the Rician bias correction filter introduced in [START_REF] Basu | Rician noise removal in diffusion tensor mri[END_REF] from 2nd-order DTI to the SPF basis. It is based on a maximum a posteriori approach so we construct the filtered volume Ê that maximizes the log-posterior probability:

log p( Ê|E) = log p(E| Ê) + log p( Ê) -log p(E) (9)
where p(E| Ê) is the likelihood term, p( Ê) is the prior or the regularization term and p(E) is the normalizing constant. We are interested in the likelihood term, thus combining Eq.( 8) and Eq.( 9) the pointwise log-likelihood becomes

log p(E| Ê, σ) = log E σ 2 - (E 2 + Ê2 ) 2σ 2 + log I 0 E • Ê σ 2 = ψ( Ê) (10)
Fig. 3.2 illustrates variation of the opposite function with scalar values of Ê when E = 1 and σ = 0.5. The energy is low when E ≈ Ê and increases with their dissimilitudes. Note that σ has to be known a priori and can be either retrieved as a parameter specific to the MR scanner, or can be either computed from a uniform area as described in [START_REF] Sijbers | Estimation of the noise in magnitude mr images[END_REF]. Combining Eq.( 7) and the derivative of Eq.( 10) with respect to A j gives the PDE adapted to Rician noise,

∂A j ∂t = ns k M k,j σ 2   -Êk + E k   I 1 E k • Êk σ 2 I 0 E k • Êk σ 2     + α r div(ϕ(||∇A||)) (11)

Regularization function ϕ

Regarding the spatial regularization, various functions ϕ of the image processing literature can be proposed as long as it preserves important features of the image. Indeed, regularization should be strong on homogeneous area (low ||∇A||), and preserve contours not only between isotropic and anisotropic regions but also among voxels with different number of fibers (large ||∇A||). We drop the angular and radial regularization step of Eq.4 in [START_REF] Assemlal | Efficient computation of pdf-based characteristics from diffusion mr signal[END_REF] since the spatial regularization puts sufficient constraints on the diffusion signal to be estimated. Our experiments have confirmed that combining these regularization to the spatial one is useless.

Experiments

In this section, we present results of our method on both synthetic simulations and real human brain data-sets. We first focus on a comparison between the damped least square estimation introduced in previous work and the robust variational framework introduced in this paper. Then we illustrate the flexibility of the proposed approach with a comparison of ODF computed with QBI method as proposed in [START_REF] Descoteaux | Regularized, fast and robust analytical q-ball imaging[END_REF] and with our method [START_REF] Assemlal | Efficient computation of pdf-based characteristics from diffusion mr signal[END_REF] on an in-vivo dataset.

Numerical Simulations

We have applied the above scheme to the simulations of a single fiber and crossing fiber configurations. The following synthetic multi-exponential model was used to generate data, In order to assess the robustness to noise of our proposed variational framework, we produced a synthetic phantom of crossing fibers (horizontal and vertical networks) surrounded by water regions (upper left area) (c.f . Fig. 3a). To simulate dMRI acquisitions, we added Rician noise of variance σ to the signal E which was then sampled using medium resolution protocol as described above.

E(q) = N b k=1 f k exp -q T D k q where N b k=1 f k = 1.
Fig. 2 shows the results of a comparison between the reconstruction of E with Gaussian and the Rician likelihood functions on a noisy dataset (Fig. 2b). A postprocessing contrast enhancement with the same parameters was applied to all images (a-f) to highlight artifacts. Although the Gaussian function is classically used in the least square minimization [START_REF] Assemlal | Efficient computation of pdf-based characteristics from diffusion mr signal[END_REF], it is not robust to noise and creates undesirable radial oscillations at high q values (Fig. 2c). On the contrary, the Rician likelihood function strongly attenuate this drawback and gives a correct estimation of E (Fig. 2d).

The Generalized Fractional Anisotropy (GFA) measure [START_REF] Tuch | High angular resolution diffusion imaging of the human brain[END_REF] in Fig. 3 is a generalization of the fractional anisotropy (FA) measure of DTI. Each image was normalized independently to enhance visualization contrast. This qualitative comparison highlights the need for spatial regularization within the estimation process. Indeed, GFA is an adequate measure to have insight on the global coherence of the dataset volume estimation since every voxel is summarized by a scalar value. When it comes to noisy input data, regularization greatly improves the spatial coherence of the volume estimation as illustrated in Fig. 3(c,d). It is worth noting that the gradient norm ||∇A|| is an adequate measure to set apart isotropic area from anisotropic area and subsequently, divergence div(ϕ||∇A||) performs well in regularizing homogeneous area without degrading the contours. We computed statistics on the performance of the PDE estimation with various likelihood and regularization functions ψ and ϕ. The PSNR (Peak Signal to Noise Ratio) between the reconstruction and the original data stands for the PSNR between the ground truth dataset volume E and its estimation Ê. Values of Ê were restricted to [0, 1] in order to reflect the signal attenuation properties. Fig. 4a, illustrates the PSNR(E, Ê) of the reconstruction versus the quality of input datasets. Out of the results, the Rician likelihood function outperforms the Gaussian function and greatly improves the PSNR of the estimation. This gap can be explained by a more robust radial fit thanks to the Rician likelihood function as shown in Fig. 2(c,d). Indeed, the number of radial sample in this experiment is very limited (2 shells) whereas the number of angular samples is sufficient (42 directions).

Fig. 4b shows influence of regularization function ϕ on the PSNR of the reconstruction Ê. Although this brings modest improvements, the spatial regularization ensures numerical stability of the estimation by adding constraints when there are only very few samples available. Besides, this also brings stability to the fiber-tracking algorithms and helps to better estimate the white matter nerve fibers tracks [START_REF] Assemlal | Fiber tracking on hardi data using robust odf fields[END_REF].

In vivo experiments

Diffusion-weighted images were acquired in two shells along 32 directions at b = 1000 s/mm 2 and b = 3000 s/mm 2 , and a single image at b = 0 s/mm 2 . Thus, there were a total of 65 images acquired in a sequence of 15 minutes on a 3 T Philips scanner. The SENSE parallel imaging protocol was used with a factor of acceleration set to 2; and only 80% of the k-space was acquired. Matrix size was 112 × 112 × 60 and the image resolution was 2 × 2 × 2 mm 3 . Repetition time was TR=11490 ms, echo time was TE=85 ms. Time between two pulses and time of diffusion gradients were respectively ∆ = 42.2 ms and δ = 26.3 ms. Terms up to N = 1 and L = 4 were used in the calculations. Computations were done in less than an hour on a 3 Ghz processor, and includes calculations of SPF coefficients and projections along the 642 directions for the whole data-set 112 × 112 × 60 volume.

Representative images of GFA on the data are presented in Fig. 5 and reveals microstructures around the genu of the corpus callosum. The first line of Fig. 5 are the results from previous work [START_REF] Assemlal | Efficient computation of pdf-based characteristics from diffusion mr signal[END_REF] and shows how it compares to standard dMRI methods. DTI performs well in corpus callosum but fails in voxel with orientational heterogeneity as shown in Fig. 5b. QBI can successfully retrieve multiple fibers orientations using the analytical stable ODF reconstruction de- scribed in [START_REF] Descoteaux | Regularized, fast and robust analytical q-ball imaging[END_REF] but is sensitive to noise, especially in region of cerebrospinal fluid (c.f . Fig. 5c). On the contrary, the ODF obtained by the SPF estimation approach does not fall into this pitfall (Fig. 5d), it successfully retrieves anisotropic shapes in brain white matter fibers regions and isotropic shape in cerebrospinal fluid area (c.f . Fig. 5d). The second line of Fig. 5 shows the performances of the proposed variational framework. Out of the results, Rician likelihood function does not modify much the GFA map computed on ODF (c.f . Fig. 5(d,e)). It was expected as artifacts on Ê are mostly radial distortion and GFA focuses on angular variations. However, spatial regularization strongly influence results, depending on the regularization strength α r as illustrated in Fig. 5(f-h).

Discussion and conclusion

In this paper, we proposed a variational approach which robustly estimates at a stretch the whole volume of PDF functions as a set of Spherical Polar Fourier (SPF) coefficients. This is done by minimizing an energy that simultaneously considers the Rician model of the MRI noise and the regularization on spatial constraints. Results demonstrate that the ability to reconstruct a voxel taking the whole neighborhood information into account strongly improve the spatial coherence of the reconstruction. Besides, fiber-tracking is unstable on noisy datasets and this last property may greatly improve the ability to recover reliable and accurate intra-voxel fibers distributions within the human brain.

Fig. 1 .

 1 Fig. 1. Energy associated to respectively Gaussian and Rician likelihood ψ functions. Note the bias introduced by the Rician function on low SNR data. E = 1 and σ = 0.5.

  The symbol N b stands for the number of fibers and D k is a 3 × 3 symmetric definite positive matrix defining the diffusion anisotropy. Diffusion images were synthesized following 3 sampling protocols: low resolution (1 shell b = 3000 s/mm 2 ), medium resolution (2 shells b = {1000 , 3000} s/mm 2 ) and high resolution (5 shells b = {500 , 1000 , 1700 , 2400 , 3000} s/mm 2 ) along with a single baseline image acquired at b = 0 s/mm 2 . Each shells is composed of 42 directions along the edges of a subdivided icosahedron. Estimation parameters were chosen empirically for each sampling protocol: low resolution {N = 0, L = 4, γ = 100, λ N = 0, λ L = 6 × 10 -5 }, medium resolution {N = 1, L = 4, γ = 70, λ N = λ L = 0}, high resolution {N = 4, L = 4, γ = 50, λ N = λ L = 0}.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Qualitative comparison between Gaussian and Rician likelihood functions influence on reconstruction of an noisy dataset. PSNR(noisy,original)=18.5. Graphics (e,f) are the absolute difference between respectively (c,d) and (a).

  a) Comparison of likelihood functions. b) Comparison of regularization functions.

Fig. 4 .

 4 Fig.[START_REF] Wedeen | Mapping fiber orientation spectra in cerebral white matter with fourier transform diffusion mri[END_REF]. Synthetic phantom of networks of crossing fibers (c.f . Fig.3a). a) Performances of likelihood functions on increasing levels of noise. b) Performances of regularization functions on increasing regularization strength αr.

  f) Soft Reg. g) Med. Reg. h) Strong Reg.

Fig. 5 .

 5 Fig. 5. Comparison of GFA [5] on region of corpus callosum and lateral ventricles. (a) Baseline image, q = 0. (b) DTI anisotropy map. (c) Q-Ball Imaging. (d) Previous work using damped least-square estimation. (e) Variational approach using Rician likelihood function. (f-g) Variational approach using Rician likelihood function + Hyper Surface regularization function.

Table 1 .

 1 A list of major notations used in this paper

	Symbol Description	Symbol Description
	PDF	Probability Density Function p, k	Displacement vectors in R 3
	ODF	Orientation Density Function P (p)	Average displacement probabilities
	FRT	Funk-Radon Transform	q	Diffusion space vector in R 3
	SH	Spherical Harmonics	S(q)	MR signal at diffusion gradient q
	SPF	Spherical Polar Fourier	E(q)	MR signal attenuation S(q)/S(0)
	a nlm	SPF expansion coefficient	G(k)	PDF characteristic at point k
		at order n, l and m	h k (p) projection function of G(k)