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 :

In the above equation, SO(3) is the set of three dimensional rotations and S² the three dimensional 2-sphere {ω∈ℝ3;||ω|| 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):

where ∀n∈ℕ,∀m∈[|-n,n|],gnm are the coefficients of the expansion of g over the SHs Ynm and ∀n∈ℕ,∀m∈[|-n,n|],∀j∈[|-n,n|],hnm,j are the coefficients of the expansion of h over the D-Wigner functions Dnm,j.

Let C a signal axisymmetric around uz→ 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 ∀m≠0,∀n,nodd,Cnm=0. Therefore, we have the SH expansion:

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 Dnm,j (Healy et al., 1998):

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

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

In the above equation, Dnm,0¯(ω) denotes the complex conjugate of the complex number Dnm,0(ω). Although the Dnm,0 are not defined over the unit sphere, consider that the Dnm,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:

where αn=4π2n+1.

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

where U_u is defined by ∀ω∈S2,Uu(ω)=∑n=0∞∑m=-nnαnYnm¯(u)Ynm(ω).

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:

where both Cun and Uun are vectors of length 2n+1. It is to be noted that Uun 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 Uun and Cn0 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.

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 uz→, i.e.:

where the ν_k are the mix coefficients and verify ∑k=1Kνk=1.

Equation 3 gives the SH expansion formulation:

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 Uukn. The per-direction responses are given by the νkCk,n0, while the sparse ODF, with exactly K directions for which it is non-zero, is given by the Uukn.

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 Ck,n0 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 Sn0 from prior knowledge of a precomputed dictionary

of convolution filters associated with rotations V = {v₁, v₂, ..., v_NR} at degree n₀. The OMP was chosen because the correlation between Uv1n0 and Uv2n0 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

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 Rk(ω)=∑n=1NνkCk,n0Yn0(ω). Note that for n = 0, the matrix (Uv00...UvK0) 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.

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 uz→ with gamma distributed radii, with all possible combinations of parameters values intra-axonal diffusivity Din∈{1.5,2,2.25,2.5,3}×10-9m2.s-1, extra-axonal diffusivity Dex∈{1,1.5,2,2.5,3}×10-9m2.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:

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λ||z2)exp(-Bλ⊥(1-z2)), 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.

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 = 3000s.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

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 λ^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).

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 000s.mm⁻², and up to order 10 for B = 10 000s.mm⁻² for SBSD and CSD with Dipy. For the MSMT-CSD and CSD run with MRtrix3 at a b-value B = 5 000s.mm⁻² on in-vivo data, default parameters were used.

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).

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.

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 000s.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 000s.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 000s.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.

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.

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 Uvkn 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.

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.

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).