Spherical Harmonics

For Spherical Deconvolution (SD) as implemented in MRtrix, processing is done in the Spherical Harmonic (SH) basis; this mathematical formulation provides a smooth representation of data distributed on the sphere. When we do SD, the resulting Fibre Orientation Distributions (FODs) are written to an image. These FOD images contain coefficients in this SH basis, that when interpreted correctly, produce the FOD butterflies we all know and love. If you’ve ever looked at the raw image volumes from an FOD image, you’ll know that all but the first one are basically not interpretable.

What are Spherical Harmonics?

Spherical harmonics are special functions defined on the surface of a sphere. They form a complete orthonormal set and can therefore be used to represent any well-behaved spherical function. In many ways, they are the equivalent to the Fourier series for functions defined over spherical (rather than Cartesian) coordinates. They are defined as:

\[Y_l^m(\theta,\phi) = \sqrt{\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}} P_l^m(\cos \theta) e^{im\phi}\]

with integer order $l$ and phase $m$, where $l \geq 0$ and $-l \leq m \leq l$ (note that the terms degree and order are also commonly used to denote $l$ & $m$ respectively), and associated Legendre polynomials $P_l^m$. The harmonic order $l$ corresponds to the angular frequency of the basis function; for example, all $l=2$ SH basis functions feature 2 full oscillations around some equator on the sphere. The harmonic phase $m$ correspond to different orthogonal modes at this frequency; e.g. where the oscillations occur around a different plane.

Any well-behaved function on the sphere $f(\theta,\phi)$ can be expressed as its spherical harmonic expansion:

\[f(\theta,\phi) = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} c_l^m Y_l^m(\theta,\phi)\]

For smooth functions that have negligible high angular frequency content, the series can be truncated at some suitable maximum harmonic order $l_\text{max}$ with little to no loss of accuracy:

\[f(\theta,\phi) = \sum_{l=0}^{l_\text{max}} \sum_{m=-l}^{l} c_l^m Y_l^m(\theta,\phi)\]

The spherical harmonic series therefore provides a compact represention for smooth functions on the sphere. Moreover, due to its formulation, it has many compelling mathematical properties that simplify otherwise complex operations, including notably spherical (de)convolution.

Formulation used in MRtrix3

Due to the nature of the problems addressed in diffusion MRI, in MRtrix3 a simplified version of the SH series is used:

The data involved are real (the phase information is invariably discarded due to its instability to motion), so we can use a real basis with no imaginary components.
The problems involved all exhibit antipodal symmetry (i.e. symmetry about the origin, $f(\mathbf{x}) = f(-\mathbf{x})$), so we can ignore all odd order terms in the series (since these correspond to strictly antisymmetric terms).

The SH basis functions $Y_{lm}(\theta,\phi)$ used in MRtrix3 are therefore:

\[\begin{split}Y_{lm}(\theta,\phi) = \begin{cases} 0 & \text{if $l$ is odd}, \\ \sqrt{2} \: \text{Im} \left[ Y_l^{-m}(\theta,\phi) \right] & \text{if $m < 0$},\\ Y_l^0(\theta,\phi) & \text{if $m = 0$},\\ \sqrt{2} \: \text{Re} \left[ Y_l^m(\theta,\phi) \right] & \text{if $m > 0$},\\ \end{cases}\end{split}\]

Storage conventions

Images that contain spherical harmonic coefficients are stored as 4-dimensional images, with each voxel’s coefficients stored along the fourth axis. Only the even degree coefficients are stored (since odd $l$ coefficients are assumed to be zero).

The SH coefficients $c_{lm}$ corresponding to the basis functions $Y_{lm}(\theta,\phi)$ are stored in the corresponding image volume at index $V_{lm}$ according to the following equation:

$V_{lm} = \frac{1}{2} l(l+1) + m$

The first few volumes of the image therefore correspond to SH coefficients as follows:

Volume $V_{lm}$	$c_{lm}$
0	$l=0$, $m=0$
1	$l=2$, $m=-2$
2	$l=2$, $m=-1$
3	$l=2$, $m=0$
4	$l=2$, $m=1$
5	$l=2$, $m=2$
6	$l=4$, $m=-4$
7	$l=4$, $m=-3$
…	…

The total number of volumes N in the image depends on the highest angular frequency band included in the series, referred to as the maximal spherical harmonic order $l_\text{max}$:

$N= \frac{1}{2} (l_\text{max}+1) (l_\text{max}+2)$

$l_\text{max}$	$N$
0	1
2	6
4	15
6	28
8	45
10	66
12	91
…	…

Representation of response functions

The response functions required when performing spherical (de)convolution correspond to axially symmetric kernels (they typically represent the ideal signal for a coherently aligned bundle of fibres aligned with the $z$ axis). Due to this symmetry, all $m \neq 0$ coefficients can be assumed to be zero. Therefore, response functions can be fully represented using only their even order $l$, zero phase $m=0$ coefficients (the so-called zonal harmonics).

Response functions are stored in plain text files. A vector of values stored on one row of such a file are interpreted as these even $l$, $m=0$ terms. The number of coefficients stored for a response function of a given maximal harmonic order $l_\text{max}$ is $1+l_\text{max}/2$.

Response files can contain multiple rows, in which case they are assumed to represent a multi-shell response, with one set of coefficients per b-value, in order of increasing b-value (i.e. the first row would normally correspond to the $b=0$ ‘shell’, with all $l>0$ terms set to zero). The b-values themselves are not stored in the response file, but are assumed to match the values in the DW encoding of the diffusion MRI dataset to be processed.

Volume \(V_{lm}\)	\(c_{lm}\)
0	\(l=0\), \(m=0\)
1	\(l=2\), \(m=-2\)
2	\(l=2\), \(m=-1\)
3	\(l=2\), \(m=0\)
4	\(l=2\), \(m=1\)
5	\(l=2\), \(m=2\)
6	\(l=4\), \(m=-4\)
7	\(l=4\), \(m=-3\)
…	…

\(l_\text{max}\)	\(N\)
0	1
2	6
4	15
6	28
8	45
10	66
12	91
…	…