Sparse Blind Spherical Deconvolution of diffusion weighted MRI

Diffusion-weighted magnetic resonance imaging provides invaluable insights into in-vivo neurological pathways. However, accurate and robust characterization of white matter fibers microstructure remains challenging. Widely used spherical deconvolution algorithms retrieve the fiber Orientation Distribution Function (ODF) by using an estimation of a response function, i.e., the signal arising from individual fascicles within a voxel. In this paper, an algorithm of blind spherical deconvolution is proposed, which only assumes the axial symmetry of the response function instead of its exact knowledge. This algorithm provides a method for estimating the peaks of the ODF in a voxel without any explicit response function, as well as a method for estimating signals associated with the peaks of the ODF, regardless of how those peaks were obtained. The two stages of the algorithm are tested on Monte Carlo simulations, as well as compared to state-of-the-art methods on real in-vivo data for the orientation retrieval task. Although the proposed algorithm was shown to attain lower angular errors than the state-of-the-art constrained spherical deconvolution algorithm on synthetic data, it was outperformed by state-of-the-art spherical deconvolution algorithms on in-vivo data. In conjunction with state-of-the art methods for axon bundles direction estimation, the proposed method showed its potential for the derivation of per-voxel per-direction metrics on synthetic as well as in-vivo data.

1 Introduction

Diffusion-weighted magnetic resonance imaging (DW-MRI) (Merboldt et al., 1985; Taylor and Bushell, 1985) is widely used for the study of tissue microstructure and fiber orientation. However, accurate in-vivo characterization of the complex and heterogeneous composition of nervous tissue remains challenging. Diffusion tensor imaging (DTI) (Basser et al., 1994) is one of the first and most commonly used diffusion model, mainly due to its ease of implementation, which only requires the acquisition of six DW-MRI images along with an additional non-weighted image. However, its underlying assumptions, particularly the presumption of Gaussian diffusion, restrict its capacity to characterize complex microstructures accurately. Following, more advanced models have been developed, trying to overcome some of the limitations of DTI, such as high angular resolution diffusion imaging (HARDI) together with Q-ball imaging (Tuch, 2004; Anderson, 2005). Q-ball imaging uses a reconstruction technique of the diffusion Orientation Distribution Function (dODF) based on the Funk-Radon Transform. As another possibility, Spherical Deconvolution (SD) (Healy et al., 1998; Tournier et al., 2004) uses a response kernel estimated from the DW-MRI data to perform the reconstruction of the ODF, later improved with a positivity constraint (Tournier et al., 2007) as well as a generalization to multiple shells and tissues (Jeurissen et al., 2014). A general formulation of the deconvolution problem was presented in Jian and Vemuri (2007) together with a comparison of various methods for solving it. Numerous additional methods were proposed and compared in Canales-Rodrguez et al. (2019). Especially, a sparse Bayesian learning framework, in which sparsity of the reconstructed ODF is promoted by assigning different variances to each entry of the function together with a zero mean Gaussian prior, has been proposed (Pisharady et al., 2018). It has the particularity of using a ball-and-stick based dictionary in order to retrieve diffusivity parameters in addition to orientations. More recently, ODF fingerprinting (Baete et al., 2019) as well as deep learning techniques based on rotationally equivariant layers were proposed (Elaldi et al., 2021, 2024) in order to improve the separability power of SD. Although Elaldi et al. (2021) highlights the formulation of the ODF reconstruction problem as blind deconvolution, a single response function is still used to train the deep neural networks in an unsupervised manner.

Consequently, SD algorithms have been somewhat limited in extracting additional signal features useful for estimating microstructural parameters such as diffusivity (in- or extra-axonal), although the ODF amplitude does contain information about parameters such as the fiber volume fraction (Raffelt et al., 2012). In order to recover per-voxel impulse responses as well as longitudinal and transverse diffusion coefficients, the Spherical Means Technique (SMT) was introduced in Kaden et al. (2016b) and generalized to multiple compartments in Kaden et al. (2016a). It uses an analytical model of the diffusion signal arising from an axon segment and leverages the invariance to the ODF of the mean signal of a voxel, provided that the response function is the same for all directions. Under similar assumptions, Anderson (2005) estimated per-voxel impulse responses in a constrained spherical deconvolution framework. It was later improved in Schultz and Groeschel (2013) with an ℓ_1/2 regularization and better estimation of the mean diffusivity of the per-voxel impulse responses. Comparably, other methods, either based on Monte-Carlo simulations or on compartment models were proposed to compute microstructural metrics or diffusion properties. Notable contributions in this domain include works by Zhang et al. (2012); Panagiotaki et al. (2014); Scherrer et al. (2016); Nedjati-Gilani et al. (2017); Rensonnet et al. (2019); Palombo et al. (2020), and De Almeida Martins et al. (2021).

In this paper, an algorithm of blind spherical deconvolution is proposed, which does not need an explicit response function nor a generative model of the signal arising from the white matter, and retrieves information about both the orientations as well as the per-voxel per-direction impulse responses of the axon bundles, under the assumption that the observed DW-MRI signal is a sparse sum of axially symmetric signals each with its own orientation. This is performed by exploiting the relationships between Spherical Harmonics (SH) and D-Wigner functions and reformulating the deconvolution problem: the rotation of a signal axially symmetric around a basis axis can be expressed as a convolution with a rotation filter which is not dependent on the observed signal but only on mathematical properties of the SHs and D-Wigner functions. First, the theory underlying the proposed algorithm is explained. Then, two experiences, one on synthetic data obtained with Monte Carlo simulation, and the other on in-vivo data, are presented. Finally, the results and contribution of this work to existing spherical deconvolution algorithms are discussed.

2 Theory

For a complete definition of the notations and conventions used in this work, see Appendix A. In this section, we derive a reformulated deconvolution problem under a sparsity assumption using SHs and D-Wigner functions. For this, we study the relationship between a convolution operator on the sphere and rotations of signals axisymmetric around $\overrightarrow{u_z}$.

2.1 Convolution on the sphere

There is a conceptual difficulty in the definition of a convolution operator on the sphere. Indeed, on the 2D plane the value of the convolution at point M can be described as the correlation of the flipped convolution kernel with the underlying function values around M. In order to define a convolution on the sphere, the straightforward analogy is to consider rotations as a substitute for translations. Therefore, the result is defined over the set of rotations SO(3), elements of which can be represented by the three Euler angles, instead of the unit sphere S₂. This conceptual difficulty led to numerous convolution operators on the sphere to be introduced in the literature, depending on the application (McEwen et al., 2007; Kennedy et al., 2011; Wei et al., 2011; Roddy and McEwen, 2021).

The convolution operator at the heart of the SD (Tournier et al., 2004), CSD (Tournier et al., 2007) and MSMT-CSD (Jeurissen et al., 2014) algorithms was introduced in Healy et al. (1998) and is defined as :

$\begin{array}{l}\forall h\in L^2(SO(3)),\forall g\in L^2(S^2),\forall \omega \in S^2,\\ (h*g)(\omega )={\displaystyle {\int }_{u\in S{O}_3}}h(u)g(u^{-1}\omega )du.\end{array}$

In the above equation, SO(3) is the set of three dimensional rotations and S² the three dimensional 2-sphere $\{\omega \in {ℝ}^3;||\omega ||{}_2=1$ }. L²(S²) and L²(SO(3)) are the set of square-integrable functions respectively defined over S² and SO(3).

A theorem similar to the classical convolution theorem holds (Healy et al., 1998):

$\begin{array}{ll}\forall h\in L^2(S{O}_3),\forall g\in L^2(S^2),{(h*g)}_n^m={\displaystyle \sum _{j=-n}^n}{h}_n^{m,j}{g}_n^j,& (1)\end{array}$

where $\forall n\in \mathbb{N},\forall m\in \left[|-n,n|\right],g_n^m$ are the coefficients of the expansion of g over the SHs $Y_n^m$ and $\forall n\in \mathbb{N},\forall m\in \left[|-n,n|\right],\forall j\in \left[|-n,n|\right],h_n^{m,j}$ are the coefficients of the expansion of h over the D-Wigner functions $D_n^{m,j}$.

2.2 Rotation of an axially symmetric signal

Let C a signal axisymmetric around $\overrightarrow{u_z}$ as well as antipodally symmetric. C can be interpreted as the response function of a parallel bundle of fibers, such as the response function of the CSD algorithm. This implies that $\forall m\ne 0,\forall n,n\mathrm{\text{odd}},C_n^m=0$. Therefore, we have the SH expansion:

$C(\omega )={\displaystyle \sum _{n=0}^{\infty }}{C}_n^0{Y}_n^0(\omega ).$

A fundamental property of the SHs gives the expression of a rotated SH as a linear combination of the SHs of the same degree, where the coefficients of the combination can explicitly be derived with the D-Wigner functions $D_n^{m,j}$ (Healy et al., 1998):

$\begin{array}{l}\forall n\in \mathbb{N},\forall j\in \left[|-n,n|\right],\forall u\in SO(3),\forall \omega \in S^2,\\ \hspace{1em}{Y}_n^j(u^{-1}\omega )=\sum _{m=-n}^n{D}_n^{m,j}(u)Y_n^m(\omega ).\end{array}$

Therefore, the rotation of C by u ∈ SO(3), noted C_u can be computed as:

$\begin{array}{l}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\forall u\in SO(3),\forall \omega \in S^2,\\ C_u(\omega )=C(u^{-1}\omega )=\sum _{n=0}^{\infty }\sum _{m=-n}^n{C}_n^0{D}_n^{m,0}(u)Y_n^m(\omega ).& (2)\end{array}$

The SHs are related to the D-Wigner functions by Healy et al. (1998)

$Y_n^m(\omega )=\sqrt{\frac{2n+1}{4\pi }}\overline{D_n^{m,0}}(\omega ).$

In the above equation, $\overline{D_n^{m,0}}(\omega )$ denotes the complex conjugate of the complex number $D_n^{m,0}(\omega )$. Although the $D_n^{m,0}$ are not defined over the unit sphere, consider that the $D_n^{m,0}$ is not dependent on the third Euler angle γ and that for a given u ∈ SO(3) represented by (α, β, 0), u can be seen as a point of S² with coordinates θ = β and Φ = α.

Thus, Equation 2 becomes:

$\forall u\in SO(3),\forall \omega \in S^2,C(u^{-1}\omega )={\displaystyle \sum _{n=0}^{\infty }}{C}_n^0{\displaystyle \sum _{m=-n}^n}{\alpha }_n\overline{Y_n^m}(u)Y_n^m(\omega ),$

where ${\alpha }_n=\sqrt{\frac{4\pi }{2n+1}}$.

Consequently, C_u can be rewritten with the convolution operator defined in Healy et al. (1998) (see Equation 1) as:

$C_u=(U_u*C),$

where U_u is defined by $\forall \omega \in S^2,U_{\mathbf{\text{u}}}(\omega )=\sum _{n=0}^{\infty }\sum _{m=-n}^n{\alpha }_n\overline{Y_n^m}(\mathbf{\text{u}})Y_n^m(\omega )$.

As in the foundational paper of Spherical Deconvolution (Tournier et al., 2004), this relationship is better visualized by regrouping coefficients of a same degree n in a vector:

$\begin{array}{ll}{C}_u^n=C_0^n{U}_u^n,& (3)\end{array}$

where both $C_{\mathbf{\text{u}}}^n$ and $U_{\mathbf{\text{u}}}^n$ are vectors of length 2n+1. It is to be noted that $U_{\mathbf{\text{u}}}^n$ can be explicitly computed given the rotation u: this is a core idea that will allow us to perform Blind Spherical Deconvolution, i.e. perform the retrieval of both $U_{\mathbf{\text{u}}}^n$ and $C_n^0$ at each voxel for each direction. This means that contrary to the CSD algorithm, in this work a different impulse response will be estimated for each axon bundle directions in every voxels. Other methods, such as the Spherical Means Technique (SMT) (Kaden et al., 2016b), have been proposed to estimate per-voxel impulse responses.

2.3 Linear superposition

The linear superposition of signals arising from intertwined axon bundles is assumed, as is common in DW-MRI (Tournier et al., 2013; Rensonnet, 2019).

Therefore, the DW-MRI signal S at a given voxel is considered to be a linear superposition of K signals C_k,uk, each being obtained by rotation of a canonical signal C_k axially symmetric around $\overrightarrow{u_z}$, i.e.:

$S(\omega )={\displaystyle \sum _{k=1}^K}{\nu }_k{C}_{u_k}(\omega )={\displaystyle \sum _{k=1}^K}{\nu }_k{C}_k(u_k^{-1}\omega ),$

where the ν_k are the mix coefficients and verify $\sum _{k=1}^K{\nu }_k=1$.

Equation 3 gives the SH expansion formulation:

$\begin{array}{ll}{S}^n={\displaystyle \sum _{k=1}^K}{\nu }_k{C}_{k,n}^0{U}_{u_k}^n.& (4)\end{array}$

Equation 4 shows that at each degree n, the vector of coefficients of the SH expansion of the observed signal is a linear mix of K vectors $U_{{\mathbf{\text{u}}}_k}^n$. The per-direction responses are given by the ${\nu }_k{C}_{k,n}^0$, while the sparse ODF, with exactly K directions for which it is non-zero, is given by the $U_{{\mathbf{\text{u}}}_k}^n$.

2.4 Summary of the proposed algorithm

The algorithm presented in this paper, and summarized in Figure 1, is the retrieval of those vectors through a complex Orthogonal Matching Pursuit (OMP) algorithm (Tropp and Gilbert, 2007; Fan et al., 2012) for a chosen degree n₀, and then the retrieval of the coefficients $C_{k,n}^0$ for all degrees n using the orientations found with the OMP. This approach will be refered to as Sparse Blind Spherical Deconvolution (SBSD) in the remainder of the work. The OMP retrieves a sparse representation of $S^{n_0}$ from prior knowledge of a precomputed dictionary

$D^{n_0}=(U_{{\text{v}}_1}^{n_0}{U}_{{\text{v}}_2}^{n_0}...U_{{\text{v}}_{N_R}}^{n_0})\in M_{2n_0+1,N_R}(ℂ),$

of convolution filters associated with rotations V = {v₁, v₂, ..., v_NR} at degree n₀. The OMP was chosen because the correlation between $U_{{\mathbf{\text{v}}}_1}^{n_0}$ and $U_{{\mathbf{\text{v}}}_2}^{n_0}$ quickly diminishes as the angular distance between v₁ and v₂ grows. The latter can be performed given any estimation {_{vk}1 ≤ k ≤ K} of the orientations, such as the peaks of the MSMT-CSD FOD, by solving

$\begin{array}{ll}{({\nu }_k{C}_{k,n}^0)}_k=arg{\displaystyle \underset{x}{min}}||S^n-{\displaystyle \sum _{k=1}^K}{x}_k{U}_{{\text{v}}_k}^n||{}^2+\lambda {(n(n+1))}^2||x||{}^2.& (5)\end{array}$

Equation 5 makes use of the regularization proposed in Descoteaux et al. (2007), with the n^th eigenvalue of the Laplace-Beltrami operator −n(n+1). The estimated coefficients are then used to reconstruct the per-voxel per-direction impulse responses $R_k(\omega )=\sum _{n=1}^N{\nu }_k{C}_{k,n}^0{Y}_n^0(\omega )$. Note that for n = 0, the matrix $(\begin{array}{c}{U}_{{\mathbf{\text{v}}}_0}^0...U_{{\mathbf{\text{v}}}_K}^0\end{array})$ is singular, due to the infinite possible combinations for attributing mean values to the estimated impulses while exactly reconstructing the observed signal. Therefore, the mean of the diffusion signal is first subtracted from itself and this zero-mean signal fitted with coefficients of degree >0. Then, for each estimated impulse, the coefficient of degree 0 is given the value of the opposite of its minimum so that the response has positive values.

Figure 1

Figure 1. Graphical summary of the proposed algorithm, named Sparse Blind Spherical Deconvolution (SBSD). The first stage estimates the orientations of the fiber bundles, and the second stage the SH expansion of each bundle. Both stages are independent, and any orientation estimate can be used as an input of the second stage. The main outputs are highlighted in red and are the bundles orientations ({_{uk}1 ≤ k ≤ K}) as well as their SHs expansion (${\left\{{\nu }_k{C}_{k,n}^0\right\}}_{\begin{array}{c}1\le k\le K\\ 0\le n\le n_{max},n\mathrm{\text{even}}\end{array}}$).

3 Methods

3.1 Experience A: validation on Monte Carlo simulations

A set C of 50 canonical signals was generated using Monte-Carlo simulations with the open-source MC-DC simulator (Rafael-Patino et al., 2020), using a substrate of cylinders oriented along $\overrightarrow{u_z}$ with gamma distributed radii, with all possible combinations of parameters values intra-axonal diffusivity $D_{in}\in \left\{1.5,2,2.25,2.5,3\right\}\times 10^{-9}{\mathrm{\text{m}}}^2{\mathrm{\text{.s}}}^{-1}$, extra-axonal diffusivity $D_{ex}\in \left\{1,1.5,2,2.5,3\right\}\times 10^{-9}{\mathrm{\text{m}}}^2{\mathrm{\text{.s}}}^{-1}$ and fiber volume fraction fvf ∈ {0.7, 0.8}. For all simulations, the gamma law was parametrized with α = 1.5 μm and β = 0.5. The simulated acquisition schemes have sampling of the sphere obtained with a Coulomb repulsion algorithm , with 100, 120, 130, 150or300 directions on the sphere for each b-value B ∈ {2000, 3000, 4000, 5000, 10000}s.mm⁻².

Then, for each Signal to Noise Ratio (SNR) ∈{100, 50, 40, 30, 20, 10}, 20 000 samples were generated using the following procedure, which is also summarized in Figure 2:

1. Two canonical signals c_i1 and c_i2 are chosen in C.

2. Two directions $\overrightarrow{u_1}$ and $\overrightarrow{u_2}$ are uniformly chosen, with a minimal angle of 25° between them, and c_i1 (resp. c_i2) is rotated along $\overrightarrow{u_1}$ (resp. $\overrightarrow{u_2}$ ).

3. The synthetic signal is computed: y = ν₁c_i1+(1−ν₁)c_i2, where ν₁ is randomly selected in [0.5, 0.85] with uniform probability.

4. y is contaminated by noise following a Rice distribution (Rice, 1944; Gudbjartsson and Patz, 1995; Alexander, 2009) with the given SNR.

Figure 2

Figure 2. Graphical summary of the method used to generate synthetic data from Monte Carlo simulations computed with the open source MC-DC simulator. First, two atoms i₁ and i₂ are randomly selected among 50 canonical signals simulated with MC-DC. They are rotated along two random directions $\overrightarrow{{\mathbf{\text{u}}}_1}$ and $\overrightarrow{{\mathbf{\text{u}}}_2}$ with a minimal separation angle of 25°. Finally, the synthetic sample is obtained by linearly mixing the two rotated atoms and adding noise.

As commonly done in DW-MRI (Canales-Rodrguez et al., 2019), SNR was defined with respect to the signal value at B = 0s.mm⁻², i.e. SNR = s(B = 0)/σ. However, this means that for higher b-values, for which the signal decreases, noise becomes stronger compared to the noiseless signal.

The CSD algorithm was applied to the synthetic samples and peaks were extracted from the reconstructed ODF. The response function used was an average of the 50 simulated signals. The estimations of the directions were also computed using SBSD for n₀∈{4, 6, 8}. Moreover, estimation of the per-voxel per-direction impulse responses was performed using ground truth bundle orientations together with Equation 5, with λ = 10⁻⁴ for all settings explored. The results were then compared both to the ground truth as well as an estimation obtained via a SMT based algorithm. For the SMT algorithm, the per-voxel response function was reconstructed from the estimated λ_|| and λ_⊥ using $R(B,z)=exp(-B{\lambda }_{||}{z}^2)exp(-B{\lambda }_{\perp }(1-z^2))$, which is the analytical model hypothesized by SMT (Kaden et al., 2016b). Both parameters (λ_|| and λ_⊥) were estimated by using all shells investigated in this work, i.e. B∈{2000, 3000, 4000, 5000, 10000}s.mm⁻². For fair comparison, the latter estimation (one per voxel) was weighted with the ground truth mix coefficient ν_k for each axon bundle to obtain per-direction impulse responses.

3.2 Experience B: application on Human Connectome Project data

The Human Connectome Project (HCP) (Setsompop et al., 2013; Van Essen et al., 2013) provides data with the aim of improving the understanding of the human brain structure, function and connectivity. In order to assess the behavior of SBSD on real in-vivo data, the algorithm was run on 5 participants of the HCP Young Adult diffusion preprocessed dataset (Feinberg et al., 2010; Moeller et al., 2010; Setsompop et al., 2012; Xu et al., 2012; Milchenko and Marcus, 2013; Sotiropoulos et al., 2013). In order to ease the reproducibility of the results presented in this paper, no additional preprocessing was performed compared to the general (Jenkinson et al., 2002; Glasser et al., 2013) as well as diffusion specific (Andersson et al., 2003; Fischl, 2012; Jenkinson et al., 2012; Andersson and Sotiropoulos, 2015, 2016) preprocessings proposed by the HCP. Results for the orientation retrieval tasks were compared to results from the CSD and MSMT-CSD algorithms. Moreover, estimation of the per-voxel per-directions impulse responses were made using the MSMT-CSD directions together with Equation 5 for the b-value B = 3000 s.mm⁻² with 64 directions on the sphere. The MSMT-CSD directions were thresholded with a relative peak amplitude of 0.1. Once the per-directions impulses are estimated, any method can be used in order to estimate per-voxel per-directions descriptors. In order to present an example of derivation of such properties from the reconstructed impulses, a descriptor of the shape of the reconstructed impulses was computed using the following procedure: a diffusion tensor was fitted to the per-direction impulse responses estimated using Equation 5 with λ = 10⁻⁴. Then, similarly to the usual FA, an Anisotropy Index (AI) is computed as

$A{I}_k=\sqrt{\frac{3}{2}\frac{{({\lambda }_{k,1}-{\widehat{\lambda }}_k)}^2+{({\lambda }_{k,2}-{\widehat{\lambda }}_k)}^2+{({\lambda }_{k,3}-{\widehat{\lambda }}_k)}^2}{{\lambda }_{k,1}^2+{\lambda }_{k,2}^2+{\lambda }_{k,3}^2}}$

for each direction k at each voxel, where λ_{k, 1}, λ_{k, 2}andλ_{k, 3} are the eigenvalues of the diffusion tensor fitted to the estimated per-direction impulse for direction k and ${\widehat{\lambda }}_k$ their mean. The resulting maps were compared to maps of standard FA. Along-tract AI was computed for the corpus callosum and the frontal aslant tract using the per-direction AI maps. This was then compared to along-tract FA obtained from the conventional FA map. The along-tract metrics were computed using the publicly available UNRAVEL Python package (Delinte, 2023), which assigns microstructural properties along a tract by examining the nearest angular peaks (Delinte et al., 2023), similarly to earlier works such as Chandio et al. (2020).

3.3 SH expansion computation

The SHs expansion of the signals were obtained via LSI with the regularization proposed in Descoteaux et al. (2007) for all experiences. The LSI was performed with a truncation including SHs up to order 8 for b-values B ≤ 5 000 s.mm⁻², and up to order 10 for B = 10 000 s.mm⁻² for SBSD and CSD with Dipy. For the MSMT-CSD and CSD run with MRtrix3 at a b-value B = 5 000 s.mm⁻² on in-vivo data, default parameters were used.

3.4 Algorithms implementations

The SBSD algorithm was implemented using Python 3.11 and makes ample use of the Quaternionic (Boyle, 2023a) and Spherical (Boyle, 2023b) packages. The SMT algorithm was performed with the implementation publicly available at https://github.com/ekaden/smt.

For applications on in-vivo data, the Python code was parallelized using the dask package (Rocklin, 2015). The source code of the implementation is available at https://github.com/cfuchs2023/SBSD_public (https://doi.org/10.5281/zenodo.10866169) together with the scripts used to generate the results presented in this work.

For Experience A (on synthetic data), CSD and ODF peaks extraction were realized with the Dipy package (Garyfallidis et al., 2014), in order to compare SBSD to an implementation of CSD in the same language. For Experience B (on in-vivo data), CSD, MSMT-CSD and ODF peaks extraction were performed with the C++ implementation available in MRtrix3 (Tournier et al., 2019), and the CSD response function was obtained using the “Tournier” algorithm (Tournier et al., 2013) while the MSMT-CSD response functions were computed using the “dHollander” algorithm (Dhollander et al., 2016, 2019).

4 Results

As the algorithm comprises two stages (see Figure 1), the performance of each stage is evaluated individually.

4.1 Experience A: validation on Monte Carlo simulations

4.1.1 Fascicle orientation estimation

Figure 3 shows the proportion of fascicles for which an error of less than 10° was achieved, which will be referred to as the accuracy of the algorithm. Performing the orientation estimation with n₀ = 4 fails to achieve accuracy higher than 0.8 even for b-value B = 10 000 s.mm⁻², SNR = 100 and 300 directions. However, it is also the choice most robust to noise, consistently achieving accuracy greater than 0.7 for SNR ≥20. The highest degree (n₀ = 8) is the choice that attains the highest accuracy of 1, for settings of high SNR = 20 together with 250 directions at b-value B = 10 000 s.mm⁻². Although such acquisition schemes are not usual in clinical settings, it is notable that the acquisition scheme used by the Human Connectome Project (HCP) in the MGH HCP Adult Diffusion dataset has a shell of 256 directions at b-value B = 10 000 s.mm⁻². However, the accuracy at n₀ = 8 quickly diminishes at SNR 10, and reliably outperforms n₀ = 6 only for b-values B≥4 000 s.mm⁻². Degree 6 outperforms degree 8 at SNR 10 for all number of directions, and for some directions at SNR 20 and 30. It also significantly outperforms degree 8 at clinically usual b-value B = 3 000 s.mm⁻², achieving accuracy between 0.8 at SNR 20 and 0.95 at SNR 30.

Figure 3

Figure 3. Proportion of fascicles for which the angular error was less than 10°, referred to as accuracy. The algorithm was run on the same synthetic data with different choices for n₀, i.e., the degree of the coefficients used for direction estimation. For each of this choice, the accuracy is plotted for different values of SNR and number of directions of the acquisition schemes.

Figure 4 shows the average angular error over all fascicles for various choices of n₀ as well as values of SNR and number of directions on the sphere. Less than 10° average error is consistently achieved for n₀ = 6, while the lowest average error is achieved for n₀ = 8 for combinations of sufficiently high SNR, b-value and number of directions. Notably, n₀ = 8 achieves significantly lower average angular error than CSD in favorable settings, while n₀ = 6 achieves a lower improvement over CSD but on a larger set of SNR, b-values and number of directions. n₀ = 4 achieves lower or similar average angular error compared to CSD at b-values B = 3 000 and 2 000 s.mm⁻².

Figure 4

Figure 4. Average angular error (in degrees) over all fascicles plotted for various choices of n₀ for the SBSD algorithm as well as CSD, for various values of SNR and number of directions on the sphere.

In in-vivo data a deviation is expected from the formulation given in Equation 4. Specifically, it is expected that the performance will drop more significantly the higher the degree, since they are more sensitive to noise and errors in the SH expansion computation.

4.1.2 Reconstruction of per-voxel per-direction impulse responses

Figure 5 show the Mean Absolute Error (MAE) of the estimated impulses using either Equation 5 or SMT compared to the ground truth, while Figure 6 shows examples of estimated impulse responses. It is notable that the proposed method outperforms SMT on all settings of SNR and number of directions for b-values B∈{2 000, 3 000, 4 000, 5 000}s.mm⁻² in terms of MAE. Both methods exhibit little sensitivity to the number of sampling points of the signal in the range studied (from 100 to 300 sampling points, i.e., from 50 to 150 acquisition gradients) for b-values B ≤ 5 000 s.mm⁻². as well as a breakdown of performance at B = 10 000 s.mm⁻². This might seem counter-intuitive and is investigated thoroughly in the Discussion section. Figure 6 shows that estimating the impulses using Equation 5 leads to ringing at both ends of the estimation (i.e., close to the poles of the 2-sphere), while SMT does not provide accurate estimation over the wide range of b-values used.

Figure 5

Figure 5. Mean Absolute Error (MAE) of the estimated impulse responses compared to the ground truth. For each synthetic voxel and each fascicle, the error vivo normalized by the mean of the ground truth impulse response in order to compare results over a wide range of b-values, for which the signal changes scale. As SMT only estimates one impulse response per voxel, it was weighted with the ground truth mix coefficient associated with each fascicle for a fair comparison with the ground truth impulse responses.

Figure 6

Figure 6. Examples of impulses estimated by using either Equation 5 (in green) or SMT (in blue) on synthetic data generated following the procedure described on Figure 2 and in Section 3.1. The SMT response was weighted by the ground truth ν_k coefficients to obtain per-direction impulse responses. The Monte-Carlo simulated ground truth is shown in black. Below, the impulse response used by CSD, which is an average of all 50 canonical signals used to generate the synthetic data, is shown in red.

4.2 Experience B: application to HCP data

4.2.1 Fascicle orientation estimation

The SBSD algorithm was run on five patients (1007, 1010 and 1016, 1019 and 1031) of the MGH HCP Adult Diffusion dataset and results for the first three are presented on Figure 7 while Figure 8 shows additional details for the first two patients.

Figure 7

Figure 7. Results of the estimation of fibers orientation performed on three patients of the HCP at b-value B = 5, 000 s.mm⁻² with 128 directions on the sphere. The estimated directions are weighted by the modulus of the associated ${\nu }_k{C}_{k,0}^{n_0}$ (green highlight in Figure 1) estimated by SBSD. The results are visualized with the usual RGB convention where red corresponds to the left-right axis, green to the anterior-posterior axis and blue to the inferior-superior axis.

Figure 8

Figure 8. Zoom on an arbitrary brain region for the two first HCP patients shown on Figure 7, showing estimated fascicles orientations in an arbitrary region including parts of the cingulum, corpus callosum, corticospinal tract and arcuate fasciculus. The estimated orientations are peaks of the ODF for CSD and MSMT-CSD. The background is a fractional anisotropy map obtained with the DTI implementation available in Mrtrix3.

Figure 7 shows that on in-vivo data, SBSD is likely able to estimate the orientation of single bundles for n₀∈{2, 4} : the corpus callosum is red-colored and the cingulum green-colored. Although the corpus callosum and corticospinal tract can still be perceived for n₀ = 6, the results are largely contaminated by noise. For n₀ = 8, the ventricles are still visible but even large structures like the cingulum and corpus callosum are indistinguishable from the gray matter. Figure 8 shows a zoom for the two first patients (1007 and 1010) in the coronal view on a region containing parts of the corpus callosum, cingulum, cerebro spinal tracts and arcuate fasciculus. The results from the proposed SBSD algorithm are compared to those obtained with CSD and MSMT-CSD. In regions occupied by a single tract, such as in the corpus callosum or cingulum, results of SBSD are broadly comparable to CSD and in agreement with MSMT-CSD. However, SBSD often only estimates one of the orientation correctly and fails to retrieve the other, such as in the area where the arcuate fasciculus crosses projections from the corpus callosum.

4.2.2 Reconstruction of per-voxel per-direction impulse responses

Figure 9 shows examples of estimated impulse responses for two patients of the HCP at different locations in the brain. Similarly to the results on synthetic data, ringing is often observed at the ends of the reconstructed impulses. Figure 10 shows the standard FA computed from the diffusion signal, as well as per-directions AI maps computed from the estimated per-direction impulses. The per-directions AI maps exhibited at least one direction with high AI in the regions where the arcuate fasciculus and corticospinal tract cross projections from the corpus callosum, where the standard FA drops significantly. Figure 11 shows along-tract FA using maps of standard FA as well as along-tract AI using per-direction AI maps. Similarly to Figure 10, the per-directions AI are generally lower than the standard FA. Moreover, for the corpus callosum, the standard FA drops significantly near and within the crossing with the frontal aslant tracts, while the per-direction AI remains close to constant. Both metrics show a significant decrease when approaching the gray matter, i.e. at the ends of the tract.

Figure 9

Figure 9. Examples of per-voxel per-directions estimated impulse responses using MSMT-CSD ODF peaks at a b-value B = 3, 000 s.mm⁻² with 64 directions on the sphere following Equation 5. Below or above each 1D plot, a 3D visualization of the estimated impulse response is shown. Together with the impulses, the directions of the MSMT-CSD peaks projected in the plane of the slices are shown as crossing sticks. As the peaks were thresholded according to their amplitudes, there can be one, two or three directions per voxel. These ℓ₂ normalized directions were used to estimate the impulse responses.

Figure 10

Figure 10. FA map obtained by running DTI with MRtrix on the diffusion signal (left column) as well as AI maps derived from impulse responses estimated using Equation 5 together with the first (resp. second) peaks of the MSMT-CSD ODF (middle (resp. right) columns) at a b-value B = 3, 000 s.mm⁻² with 64 directions on the sphere. For the latter two, no responses were estimated in voxels where no MSMT-CSD peaks exceeded the threshold, resulting in an AI of 0 in those voxels.

Figure 11

Figure 11. Along-tract AI (first row of graphs) for the corpus callosum and frontal aslant tract for patient 1007, derived from impulse responses estimated following Equation 5. The second row of graphs shows the results using the standard FA map.

5 Discussion

5.1 Orientation retrieval on Monte Carlo simulations derived data

The proposed algorithm was first tested on synthetic data obtained with Monte Carlo simulations in Experience A. SBSD achieved lower average angular errors for direction estimation in a simplified synthetic setting at high SNR (≥20) compared to CSD, as shown in Figures 3, 4. This was likely due to the use of a fine grained dictionary in the SBSD algorithm, i.e., using a high angular resolution sampling of the sphere, which was made possible thanks to the low computational complexity of the OMP. However, the relaxation of hypothesis about the response function also induced a greater sensitivity to errors in the computation of the SH expansion of the signal, mainly driven by noise and number of directions on the sphere (Figures 3, 4). Recently, efforts have been made to develop spherical Fourier transform algorithms dedicated to axially symmetric signals (Bates et al., 2016) or methods more resistant to noise for SH expansion computation (Tarar and Khalid, 2021), but no clear state-of-the-art method significantly improving over the least squares inversion proposed in Descoteaux et al. (2007) has emerged yet. Overall, this experience validated the ability of the SBSD algorithm to accurately estimate the directions of two crossing fascicles of parallel fibers, provided that both the number of directions on the shell and SNR are sufficiently high.

5.2 Impulse responses estimation on Monte Carlo simulations derived data

As SMT hypothesizes that all axon bundles have the same response regardless of their orientation, the statistical independency of the diffusion properties of the two crossing bundles in the synthetic data made it a challenging task for SMT. Moreover, SMT uses an analytical model for the response R of axon fascicles such that ln(R) is linearly dependent on B. This is not the case in in-vivo data (Jensen et al., 2005) nor in the Monte-Carlo simulated data used for Experience A, which explains why SMT cannot accurately estimate the impulse responses over a wide range of b-values, as well as its significant drop in performance for B = 10 000 s.mm⁻². Conversely, since the proposed method does not need to assume a generative model for the impulse response, it was able to recover per-direction impulse responses accurately for all b-values B ≤ 5 000 s.mm. Ringing was often observed (see Figure 6) at both ends of the impulse responses estimated following Equation 5, which is linked to core properties of the SHs. Indeed, SHs of order 0 are proportional to Legendre polynomials. Therefore, the estimations are in fact polynomials of the same degree as the maximum degree used for SHs coefficients computation. However, the ground truth is exponential in nature and peaked around z = 0. Consequently, the truncature at degree 8 or 10 induced ringing in the reconstructed impulses. This is particularly challenging for settings of high b-values, because the relative increase in noise complicates the estimation of SHs expansion while the degree needed to accurately represents the signal increases. It explains the significant drop in performance observed for the b-value B = 10 000 s.mm⁻². Moreover, the estimations were found to systematically undershoot the peak at z = 0. Depending on the method used to link those impulses to diffusion properties, this could lead to systemic bias in the estimation of parameters such as longitudinal diffusivity. Overall, the proposed method for estimating per-voxel per-direction impulse responses was found to be robust and accurate on synthetic data generated following the procedure described on Figure 2 and in Section 3.1., significantly outperforming SMT in terms of MAE (Figure 5) for all settings of SNR and number of directions for a clinically realistic range of b-values B∈{2 000, 3 000, 4 000, 5 000}s.mm⁻². However, one should note that the ODF in the synthetic data were indeed sum of Dirac impulses on the sphere, with one peak per direction.

5.3 Axon bundles orientation prediction on in-vivo data

SBSD was then run on the reduced sampling of the HCP Young Adult diffusion dataset comprising five patients (1007, 1010, 1016, 1019 and 1031) (Experience B) and results were shown for the first three patients with more details being provided for the first two. Visually, results of orientation estimation of SBSD were broadly similar to those of CSD in regions where a single axon bundle was present, with both methods agreeing with MSMT-CSD (Figures 7, 8). Although SBSD, compared to CSD, exhibited much smaller magnitude for peaks in regions where no fibers were expected such as the ventricles, SBSD performed noticeably worse in regions comprising multiple crossing fiber bundles, failing to retrieve the expected orientations in a significant number of such voxels. This was coherent with the findings of greater noise sensitivity highlighted during Experience A. Moreover, it is likely that deviations from the expected model summarized in Equation 4, such as tortuosity of axon bundles which invalidates the hypothesis of axial symmetry of the signal, also played a role, although CSD also hypothesizes the axial symmetry of its response function.

5.4 Per-direction impulse responses estimation

Since the estimation of the per-direction impulse responses corresponding to axon bundles can be performed with any orientations estimate (Equation 5), such an estimation was computed using the peaks of the MSMT-CSD ODF. Although in-vivo ODF are often assumed to be sparse and peaked (Canales-Rodrguez et al., 2019), they are not sum of Dirac impulses because the axons are not perfectly parallel and have tortuosity. Therefore, the estimated impulses are likely flatter than the ground truth because they need to account for the spread of the ODF around the Dirac peaks modeled by the $U_{{\mathbf{\text{v}}}_k}^n$ in Equation 5. Both those factors mean that depending on the method used to compute descriptors from the impulse responses, bias could be introduced by the proposed method. In Figure 11, the standard FA decreases sharply in regions where multiple tracts crosses, which is a well-known property of the DTI model which measures the anisotropy of the diffusion signal on a voxel wise scale without considering crossings of axon bundles. Conversely, for the per-direction AI maps, at least one peak was associated with a high AI value in those regions and the along-tract AI in Figure 11 was close to constant for the corpus callosum, except at the ends of the tracts where the white matter reaches the gray matter. This is an argument showing that the proposed method at least partially disentangles diffusion properties of crossing axon bundles, although the estimated impulse responses include information about the ODF. This also partially explains the decreases of the per-direction AI close to the gray matter, because the ODF becomes more spread out which flattens the estimated impulse responses. This also means that linking the estimated impulse responses to properties such as axonal density is not straightforward.

5.5 Summary of the contribution

Given the results obtained on in-vivo data, the use of SBSD for orientation estimation in clinical settings is not recommended without further improvements. However, the work presented in this paper could open the door to a new family of spherical deconvolution algorithms, allowing for more detailed characterization of axon bundles, and is of interest for people developing methods for this particular problem. Indeed, although blind or ℓ₀ penalized methods for solving spherical deconvolution problems have been proposed (Canales-Rodrguez et al., 2019), our contribution is twofold. First, the blind framework proposed does not need prior information about a generative model of the white matter response, but only an assumption about its axial symmetry. Second, our algorithm effectively leverages ℓ₀ penalization in a fast and efficient way. Both of these contributions arise from the theoretical formulation of Equation 4, which gives a sparse representation of the SH expansion of the DW-MRI signals using a dictionary of highly uncorrelated atoms which are not dependent on the observations, but only on mathematical properties of the SHs and D-Wigner functions. This equation was derived under the explicit assumption that the observed DW-MRI signal is a sparse sum of axially symmetric signals, as is hypothesized in Daducci et al. (2014) or the NNLS-BSS-EBIC method presented in Canales-Rodrguez et al. (2019). Additionally, using the peaks of MSMT-CSD, maps of per-direction AI were derived from the per-direction impulse responses estimated using Equation 5, which showed resemblance to the brain anatomy (Figures 10, 11). Therefore, the second stage of the proposed algorithm (Figure 1) in conjunction with MSMT-CSD or another robust direction estimation method could be of interest to clinicians in order to study along-tract variation of metrics between populations.

5.6 Future works

We believe SBSD to be the most straightforward approach to exploit the core ideas exposed in the Theory section. This means there is likely room for building upon the algorithm exposed in this work; future works could include leveraging information from multiple shells instead of one, and using SBSD in a multi tissues framework. Moreover, as OMP uses a least squares problem during its greedy reconstruction of the sparse representation, it was natural to also solve a least squares problem for SH expansion estimation. However, taking into account the non-Gaussian nature of noise in magnitude reconstructed DW-MRI images would lead to other formulations of the optimization problem used to estimate the SH expansion associated with each bundle (Alexander, 2009).

6 Conclusion

This work proposed an algorithm of blind spherical deconvolution called Sparse Blind Spherical Deconvolution (SBSD), which, instead of necessitating exact knowledge of a response function, only assumes its axial symmetry. This relaxation allows for the retrieval of the individual signals arising from each bundle in addition to the orientations of the bundles. On the orientation retrieval task, SBSD outperformed the CSD algorithm in some settings of high SNR (≥20) in a synthetic and simplified test case, but was more sensitive to noise and errors in the computation of the SH expansion of the signal. On real in-vivo data, SBSD performed similarly to CSD, in voxels where one fiber population associated with a given orientation was dominant. However, in voxels where roughly similar fiber populations crossed, SBSD was found to be more prone to errors or straight failures than CSD, which was likely caused by SBSD's greater sensitivity to noise and SH expansion computation errors as well as deviation from the expected axial symmetry of impulse responses due to axonal tortuosity and fiber fanning. Both methods were outperformed by MSMT-CSD, which necessitates multi-shell data. Since the retrieval of signals arising from each bundle is agnostic to the method used to estimate the relevant orientations, per-direction impulse responses were estimated using peaks from the MSMT-CSD ODF. Maps of per-direction metrics computed using an approach similar to the standard fractional anisotropy derived from the estimated per-direction impulse responses showed resemblance to brain anatomy and remained close to constant in the studied regions of crossing fascicles along the studied tracts. Therefore, this method proved its potential for the derivation of per-voxel per-direction metrics, and potential improvements to spherical deconvolution algorithms by removing the need for an explicit response function for the white matter, thus enabling the retrieval of additional information about axon bundles.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: https://db.humanconnectome.org/.

Author contributions

CF: Conceptualization, Methodology, Writing – original draft, Writing – review & editing, Data curation, Formal analysis, Software. QD: Conceptualization, Data curation, Investigation, Writing – review & editing. ND: Conceptualization, Writing – review & editing. MD: Investigation, Writing – review & editing. BM: Funding acquisition, Project administration, Supervision, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was funded in part by the MedReSyst project, supported by funds from FEDER and the Walloon region. Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11 and by the Walloon Region. QD was a research fellow of the Fonds de la Recherche Scientifique - FNRS of Belgium.

Acknowledgments

The authors thank the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research and by the McDonnell Center for Systems Neuroscience at Washington University, which provided data for Experience B. Code written by PhD, Gaëtan Rensonnet was used for the rotation of signals in the data generation framework highlighted in Figure 2.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Alexander, D. C. (2009). “Modelling, fitting and sampling in diffusion MRI,” in Visualization and Processing of Tensor Fields, eds. D. Laidlaw, J. Weickert, G. Farin, H. C. Hege, D. Hoffman, C. R. Johnson, K. Polthier (Berlin, Heidelberg: Springer Berlin Heidelberg), 3–20.

Google Scholar

Anderson, A. W. (2005). Measurement of fiber orientation distributions using high angular resolution diffusion imaging. Magn. Reson. Med. 54, 1194–1206. doi: 10.1002/mrm.20667

PubMed Abstract | Crossref Full Text | Google Scholar

Andersson, J. L., Skare, S., and Ashburner, J. (2003). How to correct susceptibility distortions in spin-echo-planar images: application to diffusion tensor imaging. Neuroimage 20, 870–888. doi: 10.1016/S1053-8119(03)00336-7

PubMed Abstract | Crossref Full Text | Google Scholar

Andersson, J. L., and Sotiropoulos, S. N. (2015). Non-parametric representation and prediction of single- and multi-shell diffusion-weighted MRI data using Gaussian processes. Neuroimage 122, 166–176. doi: 10.1016/j.neuroimage.2015.07.067

PubMed Abstract | Crossref Full Text | Google Scholar

Andersson, J. L., and Sotiropoulos, S. N. (2016). An integrated approach to correction for off-resonance effects and subject movement in diffusion MR imaging. Neuroimage 125, 1063–1078. doi: 10.1016/j.neuroimage.2015.10.019

PubMed Abstract | Crossref Full Text | Google Scholar

Baete, S. H., Cloos, M. A., Lin, Y.-C., Placantonakis, D. G., Shepherd, T., and Boada, F. E. (2019). Fingerprinting orientation distribution functions in diffusion MRI detects smaller crossing angles. Neuroimage 198, 231–241. doi: 10.1016/j.neuroimage.2019.05.024

PubMed Abstract | Crossref Full Text | Google Scholar

Basser, P., Mattiello, J., and LeBihan, D. (1994). MR diffusion tensor spectroscopy and imaging. Biophys. J. 66, 259–267. doi: 10.1016/S0006-3495(94)80775-1

PubMed Abstract | Crossref Full Text | Google Scholar

Bates, A. P., Khalid, Z., and Kennedy, R. A. (2016). An optimal dimensionality sampling scheme on the sphere with accurate and efficient spherical harmonic transform for diffusion MRI. IEEE Signal Proc. Lett. 1, 15–19. doi: 10.1109/LSP.2015.2498162

Crossref Full Text | Google Scholar

Boyle, M. (2023a). The Quaternionic Package. Zenodo. doi: 10.5281/ZENODO.10214723

Crossref Full Text | Google Scholar

Boyle, M. (2023b). The Spherical Package. Zenodo. doi: 10.5281/ZENODO.10214833

Crossref Full Text | Google Scholar

Canales-Rodrguez, E. J., Legarreta, J. H., Pizzolato, M., Rensonnet, G., Girard, G., Patino, J. R., et al. (2019). Sparse wars: a survey and comparative study of spherical deconvolution algorithms for diffusion MRI. Neuroimage 184, 140–160. doi: 10.1016/j.neuroimage.2018.08.071

PubMed Abstract | Crossref Full Text | Google Scholar

Chandio, B. Q., Risacher, S. L., Pestilli, F., Bullock, D., Yeh, F.-C., Koudoro, S., et al. (2020). Bundle analytics, a computational framework for investigating the shapes and profiles of brain pathways across populations. Sci. Rep. 10:17149. doi: 10.1038/s41598-020-74054-4

PubMed Abstract | Crossref Full Text | Google Scholar

Daducci, A., Van De Ville, D., Thiran, J.-P., and Wiaux, Y. (2014). Sparse regularization for fiber ODF reconstruction: from the suboptimality of ℓ₂ and ℓ₁ priors to ℓ₀. Med. Image Anal. 18, 820–833. doi: 10.1016/j.media.2014.01.011

PubMed Abstract | Crossref Full Text | Google Scholar

De Almeida Martins, J. P., Nilsson, M., Lampinen, B., Palombo, M., While, P. T., Westin, C.-F., et al. (2021). Neural networks for parameter estimation in microstructural MRI: application to a diffusion-relaxation model of white matter. Neuroimage 244:118601. doi: 10.1016/j.neuroimage.2021.118601

PubMed Abstract | Crossref Full Text | Google Scholar

Delinte, N. (2023). DelinteNicolas/UNRAVEL: v1.4.14. doi: 10.5281/ZENODO.10259897

Crossref Full Text | Google Scholar

Delinte, N., Dricot, L., Macq, B., Gosse, C., Van Reybroeck, M., and Rensonnet, G. (2023). Unraveling multi-fixel microstructure with tractography and angular weighting. Front. Neurosci. 17:1199568. doi: 10.3389/fnins.2023.1199568

PubMed Abstract | Crossref Full Text | Google Scholar

Descoteaux, M., Angelino, E., Fitzgibbons, S., and Deriche, R. (2007). Regularized, fast, and robust analytical Q-ball imaging. Magn. Reson. Med. 58, 497–510. doi: 10.1002/mrm.21277

PubMed Abstract | Crossref Full Text | Google Scholar

Dhollander, T., Mito, R., Raffelt, D., and Connelly, A. (2019). “Improved white matter response function estimation for 3-tissue constrained spherical deconvolution,” in Proceedings of the 27th Annual Meeting of the International Society of Magnetic Resonance in Medicine (Berkeley, CA).

Google Scholar

Dhollander, T., Raffelt, D., and Connelly, A. (2016). “Unsupervised 3-tissue response function estimation from single-shell or multi-shell diffusion MR data without a co-registered T1 image,” in ISMRM Workshop on Breaking the Barriers of Diffusion MRI.

Google Scholar

Elaldi, A., Dey, N., Kim, H., and Gerig, G. (2021). “Equivariant spherical deconvolution: Learning sparse orientation distribution functions from spherical data,” in Information Processing in Medical Imaging, eds. A. Feragen, S. Sommer, J. Schnabel, and M. Nielsen (Cham: Springer International Publishing), 267–278.

PubMed Abstract | Google Scholar

Elaldi, A., Gerig, G., and Dey, N. (2024). “(E(3) × SO(3))-equivariant networks for spherical deconvolution in diffusio MRI. Medical imaging with deep learning,” in Proceedings of Machine Learning Research, 301–319. Available online at: https://proceedings.mlr.press/v227/elaldi24a.html

PubMed Abstract | Google Scholar

Fan, R., Wan, Q., Liu, Y., Chen, H., and Zhang, X. (2012). Complex Orthogonal Matching Pursuit and Its Exact Recovery Conditions. Cornell University.

Google Scholar

Feinberg, D. A., Moeller, S., Smith, S. M., Auerbach, E., Ramanna, S., Glasser, M. F., et al. (2010). Multiplexed Echo planar imaging for sub-second whole brain FMRI and fast diffusion imaging. PLoS ONE 5:e15710. doi: 10.1371/journal.pone.0015710

PubMed Abstract | Crossref Full Text | Google Scholar

Fischl, B. (2012). FreeSurfer. Neuroimage 62, 774–781. doi: 10.1016/j.neuroimage.2012.01.021

PubMed Abstract | Crossref Full Text | Google Scholar

Garyfallidis, E., Brett, M., Amirbekian, B., Rokem, A., Van Der Walt, S., Descoteaux, M., et al. (2014). Dipy, a library for the analysis of diffusion MRI data. Front. Neuroinform. 8:8. doi: 10.3389/fninf.2014.00008

PubMed Abstract | Crossref Full Text | Google Scholar

Glasser, M. F., Sotiropoulos, S. N., Wilson, J. A., Coalson, T. S., Fischl, B., Andersson, J. L., et al. (2013). The minimal preprocessing pipelines for the Human Connectome Project. Neuroimage 80, 105–124. doi: 10.1016/j.neuroimage.2013.04.127

PubMed Abstract | Crossref Full Text | Google Scholar

Gudbjartsson, H., and Patz, S. (1995). The rician distribution of noisy MRI data. Magn. Reson. Med. 34, 910–914. doi: 10.1002/mrm.1910340618

PubMed Abstract | Crossref Full Text | Google Scholar

Healy, D. M., Hendriks, H., and Kim, P. T. (1998). Spherical deconvolution. J. Multivar. Anal. 67, 1–22. doi: 10.1006/jmva.1998.1757

Crossref Full Text | Google Scholar

Jenkinson, M., Bannister, P., Brady, M., and Smith, S. (2002). Improved optimization for the robust and accurate linear registration and motion correction of brain images. Neuroimage 17, 825–841. doi: 10.1006/nimg.2002.1132

PubMed Abstract | Crossref Full Text | Google Scholar

Jenkinson, M., Beckmann, C. F., Behrens, T. E., Woolrich, M. W., and Smith, S. M. (2012). FSL. Neuroimage 62, 782–790. doi: 10.1016/j.neuroimage.2011.09.015

PubMed Abstract | Crossref Full Text | Google Scholar

Jensen, J. H., Helpern, J. A., Ramani, A., Lu, H., and Kaczynski, K. (2005). Diffusional kurtosis imaging: the quantification of non-gaussian water diffusion by means of magnetic resonance imaging. Magn. Reson. Med. 53, 1432–1440. doi: 10.1002/mrm.20508

PubMed Abstract | Crossref Full Text | Google Scholar

Jeurissen, B., Tournier, J.-D., Dhollander, T., Connelly, A., and Sijbers, J. (2014). Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data. Neuroimage 103, 411–426. doi: 10.1016/j.neuroimage.2014.07.061

PubMed Abstract | Crossref Full Text | Google Scholar

Jian, B., and Vemuri, B. (2007). A unified computational framework for deconvolution to reconstruct multiple fibers from diffusion weighted MRI. IEEE Trans. Med. Imaging 26, 1464–1471. doi: 10.1109/TMI.2007.907552

PubMed Abstract | Crossref Full Text | Google Scholar

Kaden, E., Kelm, N. D., Carson, R. P., Does, M. D., and Alexander, D. C. (2016a). Multi-compartment microscopic diffusion imaging. Neuroimage 139, 346–359. doi: 10.1016/j.neuroimage.2016.06.002

PubMed Abstract | Crossref Full Text | Google Scholar

Kaden, E., Kruggel, F., and Alexander, D. C. (2016b). Quantitative mapping of the per-axon diffusion coefficients in brain white matter. Magn. Reson. Med. 75, 1752–1763. doi: 10.1002/mrm.25734

PubMed Abstract | Crossref Full Text | Google Scholar

Kennedy, R. A., Lamahewa, T. A., and Wei, L. (2011). On azimuthally symmetric 2-sphere convolution. Digit. Signal Process. 21, 660–666. doi: 10.1016/j.dsp.2011.05.002

Crossref Full Text | Google Scholar

McEwen, J. D., Hobson, M. P., Mortlock, D. J., and Lasenby, A. N. (2007). Fast directional continuous spherical wavelet transform algorithms. IEEE Trans. Signal Proc. 55, 520–529. doi: 10.1109/TSP.2006.887148

Crossref Full Text | Google Scholar

Merboldt, K.-D., Hanicke, W., and Frahm, J. (1985). Self-diffusion NMR imaging using stimulated echoes. J. Magnet. Reson. 64, 479–486. doi: 10.1016/0022-2364(85)90111-8

Crossref Full Text | Google Scholar

Milchenko, M., and Marcus, D. (2013). Obscuring surface anatomy in volumetric imaging data. Neuroinformatics 11, 65–75. doi: 10.1007/s12021-012-9160-3

PubMed Abstract | Crossref Full Text | Google Scholar

Moeller, S., Yacoub, E., Olman, C. A., Auerbach, E., Strupp, J., Harel, N., et al. (2010). Multiband multislice GE-EPI at 7 tesla, with 16-fold acceleration using partial parallel imaging with application to high spatial and temporal whole-brain fMRI. Magn. Reson. Med. 63, 1144–1153. doi: 10.1002/mrm.22361

PubMed Abstract | Crossref Full Text | Google Scholar

Nedjati-Gilani, G. L., Schneider, T., Hall, M. G., Cawley, N., Hill, I., Ciccarelli, O., et al. (2017). Machine learning based compartment models with permeability for white matter microstructure imaging. Neuroimage 150, 119–135. doi: 10.1016/j.neuroimage.2017.02.013

PubMed Abstract | Crossref Full Text | Google Scholar

Palombo, M., Ianus, A., Guerreri, M., Nunes, D., Alexander, D. C., Shemesh, N., et al. (2020). SANDI: a compartment-based model for non-invasive apparent soma and neurite imaging by diffusion MRI. Neuroimage 215:116835. doi: 10.1016/j.neuroimage.2020.116835

PubMed Abstract | Crossref Full Text | Google Scholar

Panagiotaki, E., Walker-Samuel, S., Siow, B., Johnson, S. P., Rajkumar, V., Pedley, R. B., et al. (2014). Noninvasive quantification of solid tumor microstructure using VERDICT MRI. Cancer Res. 74, 1902–1912. doi: 10.1158/0008-5472.CAN-13-2511

PubMed Abstract | Crossref Full Text | Google Scholar

Pisharady, P. K., Sotiropoulos, S. N., Duarte-Carvajalino, J. M., Sapiro, G., and Lenglet, C. (2018). Estimation of white matter fiber parameters from compressed multiresolution diffusion MRI using sparse Bayesian learning. Neuroimage 167, 488–503. doi: 10.1016/j.neuroimage.2017.06.052

PubMed Abstract | Crossref Full Text | Google Scholar

Rafael-Patino, J., Romascano, D., Ramirez-Manzanares, A., Canales-Rodrguez, E. J., Girard, G., and Thiran, J.-P. (2020). Robust Monte-Carlo simulations in diffusion-MRI: effect of the substrate complexity and parameter choice on the reproducibility of results. Front. Neuroinform. 14:8. doi: 10.3389/fninf.2020.00008

PubMed Abstract | Crossref Full Text | Google Scholar

Raffelt, D., Tournier, J.-D., Rose, S., Ridgway, G. R., Henderson, R., Crozier, S., et al. (2012). Apparent fibre density: a novel measure for the analysis of diffusion-weighted magnetic resonance images. Neuroimage 59, 3976–3994. doi: 10.1016/j.neuroimage.2011.10.045

PubMed Abstract | Crossref Full Text | Google Scholar

Rensonnet, G. (2019). In vivo diffusion magnetic resonance imaging of the white matter microstructure from dictionaries generated by Monte Carlo simulations: development and validation. Vaud: Lausanne, EPFL.

Google Scholar

Rensonnet, G., Scherrer, B., Girard, G., Jankovski, A., Warfield, S. K., Macq, B., et al. (2019). Towards microstructure fingerprinting: estimation of tissue properties from a dictionary of Monte Carlo diffusion MRI simulations. Neuroimage 184, 964–980. doi: 10.1016/j.neuroimage.2018.09.076

PubMed Abstract | Crossref Full Text | Google Scholar

Rice, S. O. (1944). Mathematical analysis of random noise. Bell Syst. Techn. J. 23, 282–332. doi: 10.1002/j.1538-7305.1944.tb00874.x

Crossref Full Text | Google Scholar

Rocklin, M. (2015). “Dask: Parallel computation with blocked algorithms and task scheduling,” in Proceedings of the 14th Python in Science Conference, Number 130-136 (Princeton: Citeseer).

Google Scholar

Roddy, P. J., and McEwen, J. D. (2021). Sifting convolution on the sphere. IEEE Signal Process. Lett. 28, 304–308. doi: 10.1109/LSP.2021.3050961

Crossref Full Text | Google Scholar

Scherrer, B., Schwartzman, A., Taquet, M., Sahin, M., Prabhu, S. P., and Warfield, S. K. (2016). Characterizing brain tissue by assessment of the distribution of anisotropic microstructural environments in diffusion-compartment imaging (DIAMOND). Magn. Reson. Med. 76, 963–977. doi: 10.1002/mrm.25912

PubMed Abstract | Crossref Full Text | Google Scholar

Schultz, T., and Groeschel, S. (2013). “Auto-calibrating spherical deconvolution based on odf sparsity,” in Medical Image Computing and Computer-Assisted Intervention-MICCAI 2013, eds. K. Mori, I. Sakuma, Y. Sato, C. Barillot, and N. Navab (Berlin, Heidelberg: Springer Berlin Heidelberg), 663–670.

PubMed Abstract | Google Scholar

Setsompop, K., Gagoski, B. A., Polimeni, J. R., Witzel, T., Wedeen, V. J., and Wald, L. L. (2012). Blipped-controlled aliasing in parallel imaging for simultaneous multislice echo planar imaging with reduced factor penalty. Magn. Reson. Med. 67, 1210–1224. doi: 10.1002/mrm.23097

PubMed Abstract | Crossref Full Text | Google Scholar

Setsompop, K., Kimmlingen, R., Eberlein, E., Witzel, T., Cohen-Adad, J., McNab, J., et al. (2013). Pushing the limits of in vivo diffusion MRI for the Human Connectome Project. Neuroimage 80, 220–233. doi: 10.1016/j.neuroimage.2013.05.078

PubMed Abstract | Crossref Full Text | Google Scholar

Sotiropoulos, S. N., Moeller, S., Jbabdi, S., Xu, J., Andersson, J. L., Auerbach, E. J., et al. (2013). Effects of image reconstruction on fiber orientation mapping from multichannel diffusion MRI: Reducing the noise floor using SENSE. Magn. Reson. Med. 70, 1682–1689. doi: 10.1002/mrm.24623

PubMed Abstract | Crossref Full Text | Google Scholar

Tarar, M. O., and Khalid, Z. (2021). “Reconstruction of finite rate of innovation spherical signals in the presence of noise using deep learning architecture,” in 2020 28th European Signal Processing Conference (EUSIPCO) (Amsterdam: IEEE), 1487–1491.

Google Scholar

Taylor, D. G., and Bushell, M. C. (1985). The spatial mapping of translational diffusion coefficients by the NMR imaging technique. Phys. Med. Biol. 30, 345–349. doi: 10.1088/0031-9155/30/4/009

PubMed Abstract | Crossref Full Text | Google Scholar

Tournier, J., Calamante, F., and Connelly, A. (2013). Determination of the appropriate b value and number of gradient directions for high angular resolution diffusion weighted imaging. NMR Biomed. 26, 1775–1786. doi: 10.1002/nbm.3017

PubMed Abstract | Crossref Full Text | Google Scholar

Tournier, J.-D., Calamante, F., and Connelly, A. (2007). Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution. NeuroImage, 35, 1459–1472. doi: 10.1016/j.neuroimage.2007.02.016

PubMed Abstract | Crossref Full Text | Google Scholar

Tournier, J.-D., Calamante, F., Gadian, D. G., and Connelly, A. (2004). Direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution. Neuroimage 23, 1176–1185. doi: 10.1016/j.neuroimage.2004.07.037

PubMed Abstract | Crossref Full Text | Google Scholar

Tournier, J.-D., Smith, R., Raffelt, D., Tabbara, R., Dhollander, T., Pietsch, M., et al. (2019). MRtrix3: a fast, flexible and open software framework for medical image processing and visualisation. Neuroimage 202, 116137. doi: 10.1016/j.neuroimage.2019.116137

PubMed Abstract | Crossref Full Text | Google Scholar

Tropp, J. A., and Gilbert, A. C. (2007). Signal recovery from random measurements via orthogonal matching pursuit. IEEE Trans. Informat. Theory 53, 4655–4666. doi: 10.1109/TIT.2007.909108

Crossref Full Text | Google Scholar

Tuch, D. S. (2004). Q ball imaging. Magn. Reson. Med. 52, 1358–1372. doi: 10.1002/mrm.20279

PubMed Abstract | Crossref Full Text | Google Scholar

Van Essen, D. C., Smith, S. M., Barch, D. M., Behrens, T. E., Yacoub, E., and Ugurbil, K. (2013). The WU-Minn human connectome project: an overview. Neuroimage 80, 62–79. doi: 10.1016/j.neuroimage.2013.05.041

PubMed Abstract | Crossref Full Text | Google Scholar

Wei, L., Kennedy, R. A., and Lamahewa, T. A. (2011). Quadratic variational framework for signal design on the 2-sphere. IEEE Trans. Signal Proc. 59, 5243–5252. doi: 10.1109/TSP.2011.2162506

Crossref Full Text | Google Scholar

Xu, J., Moeller, S., Strupp, J., Auerbach, E. J., Chen, L., Feinberg, D. A., et al. (2012). “Highly accelerated whole brain imaging using aligned-blipped-controlled-aliasing multiband EPI,” in Proceedings of the International Society for Magnetic Resonance in Medicine (Berkeley, CA).

Google Scholar

Zhang, H., Schneider, T., Wheeler-Kingshott, C. A., and Alexander, D. C. (2012). NODDI: Practical in vivo neurite orientation dispersion and density imaging of the human brain. Neuroimage 61, 1000–1016. doi: 10.1016/j.neuroimage.2012.03.072

PubMed Abstract | Crossref Full Text | Google Scholar

Appendix A

Angles conventions

Spherical coordinates are defined according to the common mathematical convention where θ is the colatitude [i.e. the oriented angle $(\overrightarrow{u_Z},\overrightarrow{OM})$] and Φ the longitude [i.e., the oriented angle $(\overrightarrow{u_x},{\overrightarrow{OM}}_z)$ where ${\overrightarrow{OM}}_z$ is the projection of $\overrightarrow{OM}$ in the plane z = 0]. For Euler angles, the common zyz convention is used together with the notations α, β, γ : γ is the first rotation, around $\overrightarrow{u_z}$, then β around $\overrightarrow{u_y}$ and finally α around $\overrightarrow{u_z}$. Consequently, the rotation of a signal axysymmetric around $\overrightarrow{u_z}$ is only dependent on α and β.

Spherical harmonics and D-Wigner functions

Multiple conventions can be found in the literature for the definition of spherical harmonics (SHs) or D-Wigner functions. In this work, the following definition is used for spherical harmonics:

$\begin{array}{l}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\forall n\in \mathbb{N}\text{,}\forall m\in \left[|-n,n|\right],\\ Y_n^m(\theta ,\Phi )=\sqrt{\frac{(2n+1)(n-\left|m\right|)!}{4\pi (n+\left|m\right|)}}{P}_n^{\left|m\right|}(cos(\theta ))e^{im\Phi },\end{array}$

where

$P_n^m(x)={(-1)}^m{(1-x^2)}^{\frac{m}{2}}\frac{{\text{d}}^m{P}_n(x)}{\text{d}{x}^m}$

are the associated Legendre polynomials and

$P_n(x)=\frac{1}{2^{n}n!}\frac{{\text{d}}^n}{\text{d}{x}^n}{(x^2-1)}^n$

the Legendre polynomials. SHs form an orthonormal basis of the Hilbert space L²(S²) and therefore, any f ∈ L²(S²) can be defined with its expansion over the SHs $\forall \omega \in S^2,f(\omega )=\sum _{n=0}^{\infty }\sum _{m=-n}^n{f}_n^m{Y}_n^m(\omega )$. The coefficients $f_n^m$ are defined by $f_n^m={\int }_{\omega \in S^2}f(\omega )\overline{Y_n^m(\omega )}d\omega$. Note that the antipodal symmetry of diffusion signals implies that all $f_n^m=0$ for n odd.

In this work, the following definition is used for D-Wigner functions

$\begin{array}{l}\forall n\in \mathbb{N}\text{,}\forall m\in \left[|-n,n|\right]\text{,}\forall q\in \left[|-n,n|\right]\text{,}\\ \hspace{1em}{D}_n^{m,q}(\alpha ,\beta ,\gamma )=e^{-jm\alpha }{d}_n^{m,q}(\beta )e^{-jq\gamma },\end{array}$

with the d-Wigner functions $d_n^{m,q}$

$\begin{array}{l}\begin{array}{c}{d}_n^{m,q}(x)=j^{m-q}\frac{sin{(x)}^{q-m}{(1+cos(x))}^m}{2^n{((n+m)!(n-q)!)}^{1/2}}{\left[\frac{(n-q)!}{(n+q)!}\right]}^{1/2}\\ \frac{{\text{d}}^{n+q}}{\text{d}cos{(x)}^{n+q}}\left[{(cos(x)-1)}^{n+m}{(cos(x)+1)}^{n-m}\right].\end{array}\end{array}$