Constrained Spherical Deconvolution (CSD) is a well-established and popular model for estimating white matter fibre orientations from diffusion MRI data¹. Rather than estimating a finite number of discrete fibre populations per voxel, CSD instead provides a continuous representation of fibre density over all orientations - the Fibre Orientation Distribution (FOD). However, the spherical harmonic (SH) basis in which this estimate is provided is an indirect representation of the fibre density information, and does not fully exploit the expected sparse nature of the FOD. A common approach is to extract the peak orientations of the FOD; however this fails to capture valuable information regarding its precise shape. We propose a sparse parametrization of the FOD, which accurately captures all features of the FOD using a discrete set of basis functions, and is more amenable to subsequent processing such as tractography.

In spherical deconvolution-based models, the empirical diffusion signal is considered to be the spherical convolution of the FOD with the 'response function' $$$RF$$$ (the expected diffusion signal profile for a single coherent fibre population oriented along the z-axis):

$$DWI=\int_{0}^{\pi/2}\int_{0}^{2\pi}FOD(\theta,\phi)\otimes RF(\theta)d\phi d\theta$$

The FOD is estimated from the data through a deconvolution operation; in CSD, this is a continuous function on $$$S^2$$$ estimated in the SH basis. A small number of finite fibre bundles is expected to be found in each voxel, but the SH basis does not exploit this; furthermore, each may individually demonstrate fibre orientation dispersion, bending or fanning within the voxel, but this is only indirectly reconstructed.

Our proposed re-parametrization of the FOD is as follows:

$$FOD(\widehat{v})=\sum_{f}sz_f\times\delta_{lmax}(0.5cos^{-1} (1-2\Theta_{\widehat{v}f}))$$

$$\Theta_{\widehat{v}f}=\widehat{v}\cdot(B_{f}(\theta,\phi,\psi,\kappa_1,\kappa_2)\times\widehat{v})$$

Each discrete fibre bundle $$$f$$$ in the voxel is parametrized by a scaling factor $$$sz_{f}$$$, and matrix $$$B_{f}(\theta,\phi,\psi,\kappa_1,\kappa_2)$$$, which encapsulates peak orientation and anisotropic dispersion (as demonstrated in the Ball-and-Rackets model²). To reconstruct the FOD amplitude in any given unit direction $$$\widehat{v}$$$ on the sphere, an 'effective' angle from each bundle is derived based on matrices $$$B_{f}$$$, and used to access a lookup table containing values of an approximate delta function of the appropriate spherical harmonic order $$$\delta_{lmax}(\theta)$$$ (to the first zero crossing only). The anisotropic dispersion parameters $$$\kappa_1$$$ and $$$\kappa_2$$$ incorporated within $$$B_{f}$$$ permit angular dilation of the delta function by scaling this effective angle.

Following an initial FOD segmentation based on separation at troughs³, each FOD lobe is represented using a finite number of fibre bundle elements using the Levenberg-Marquadt algorithm, with heuristic model complexity selection; this process is shown for an example FOD in Figure 1.

Figure 2 shows how the re-parametrization algorithm is able to accurately reconstruct the FOD profile across a range of crossing and fanning angles. In fact, as the fanning angle reaches 90 degrees, the FOD displays two discrete peaks, yet the re-parametrization correctly identifies this as a single fibre population (as a two-fibre model does not properly fit the precise FOD shape).

Figure 3 shows how the parameter $$$\kappa_1$$$, which parametrizes the extent of dilation along the main fanning axis, varies as a function of the simulated fanning angle. This graph indicates the potential for using the dilation coefficients extracted from this parametrization to infer the fanning angle in the underlying fibre bundles.

Our method shares features with a number of diffusion models that aim to estimate within-voxel fibre orientation dispersion^2,4,5. An important advantage of our approach, compared to discrete mixture models operating on the diffusion data itself, is that if the model complexity is either reduced (by discarding very small fibre populations) or erroneously under-estimated, parameter estimates for the remaining fibre bundles are not correspondingly biased. Furthermore, our method should be less sensitive to parameter initialisation, as the FOD itself provides appropriate initial conditions.

A similar parametrization of the FOD has been proposed previously⁶. Unlike their approach, we use the entire FOD to drive the model fit (rather than a small neighbourhood of samples around each peak), and use the approximate delta function of the SH basis as the mixture model basis function to drive a highly accurate fit (unlike the Bingham distribution, which is non-zero over the whole half-sphere and therefore not consistent with the FOD).

We have demonstrated that our approach extracts per-fibre-bundle orientation and dispersion parameters from the FOD. Additionally, we have shown that the fibre model dispersion parameters may give a direct estimate of the underlying fibre fanning angle. One potential application for this result is the design of more advanced tractography algorithms that are informed by fibre bundle fanning angles even in the presence of complex crossing fibre configurations.