# Spherical Harmonic bases

Spherical Harmonics (SH) are functions defined on the sphere. A collection of SH can be used as a basis function to represent and reconstruct any function on the surface of a unit sphere.

Spherical harmonics are orthonormal functions defined by:

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

where \(l\) is the order, \(m\) is the phase factor, \(P_l^m\) is an associated \(l\)-th order, \(m\)-th phase factor Legendre polynomial, and \((\theta, \phi)\) is the representation of the direction vector in spherical coordinates. The relation between \(Y_l^m\) and \(Y_l^{-m}\) is given by:

\[Y_l^{-m}(\theta, \phi) = (-1)^m \overline{Y_l^m}\]

where \(\overline{Y_l^m}\) is the complex conjugate of \(Y_l^m\) defined as \(\overline{Y_l^m} = \Re(Y_l^m) - \Im(Y_l^m)\).

A function \(f(\theta, \phi)\) can be represented using a spherical harmonics basis using the spherical harmonics coefficients \(a_l^m\), which can be computed using the expression:

\[a_l^m = \int_S f(\theta, \phi) Y_l^m(\theta, \phi) ds\]

Once the coefficients are computed, the function \(f(\theta, \phi)\) can be computed as:

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

In HARDI, the Orientation Distribution Function (ODF) is a function on the sphere. Therefore, SH functions offer the ideal framework for reconstructing the ODF. Descoteaux et al.[1] use the Q-Ball Imaging (QBI) formalization to recover the ODF, while Tournier et al.[2] use the Spherical Deconvolution (SD) framework.

Several modified SH bases have been proposed in the diffusion imaging literature for the computation of the ODF. DIPY implements two of these in the `shm` module. Below are the formal definitions taken directly from the literature.

- The basis proposed by Descoteaux et al.[1]:

\[\begin{split}Y_i(\theta, \phi) = \begin{cases} \sqrt{2} * \Re(Y_l^m(\theta, \phi)) & -l \leq m < 0, \\ Y_l^0(\theta, \phi) & m = 0, \\ \sqrt{2} * \Im(Y_l^m(\theta, \phi)) & 0 < m \leq l \end{cases}\end{split}\]

- The basis proposed by Tournier et al.[2]:

\[\begin{split}Y_i(\theta, \phi) = \begin{cases} \Im(Y_l^m(\theta, \phi)) & -l \leq m < 0, \\ Y_l^0(\theta, \phi) & m = 0, \\ \Re(Y_l^m(\theta, \phi)) & 0 < m \leq l \end{cases}\end{split}\]

In both cases, \(\Re\) denotes the real part of the spherical harmonic basis, and \(\Im\) denotes the imaginary part. The SH bases are both orthogonal and real. Moreover, the descoteaux07 basis is orthonormal.

By alternately selecting the real or imaginary part of the original SH basis, the modified SH bases have the properties of being both orthogonal and real. Moreover, due to the presence of the \(\sqrt{2}\) factor, the basis proposed by Descoteaux et al. is orthonormal.

The SH bases implemented in DIPY for versions 1.2 and below differ slightly from the literature. Their implementation is given below.

- The `descoteaux07` basis is based on the one proposed by Descoteaux et al.[1] and is given by:

\[\begin{split}Y_i(\theta, \phi) = \begin{cases} \sqrt{2} * \Re(Y_l^{|m|}(\theta, \phi)) & -l \leq m < 0, \\ Y_l^0(\theta, \phi) & m = 0, \\ \sqrt{2} * \Im(Y_l^m(\theta, \phi)) & 0 < m \leq l \end{cases}\end{split}\]

- The `tournier07` basis is based on the one proposed by Tournier et al.[2] and is given by:

\[\begin{split}Y_i(\theta, \phi) = \begin{cases} \Im(Y_l^{|m|}(\theta, \phi)) & -l \leq m < 0, \\ Y_l^0(\theta, \phi) & m = 0, \\ \Re(Y_l^m(\theta, \phi)) & 0 < m \leq l \end{cases}\end{split}\]

These bases differ from the literature by the presence of an absolute value around \(m\) when \(m < 0\). Due to relations \(-p = |p| ; \forall p < 0\) and \(Y_l^{-m}(\theta, \phi) = (-1)^m \overline{Y_l^m}\), the effect of this change is a sign flip for the SH functions of even degree \(m < 0\). This has no effect on the mathematical properties of each basis.

The `tournier07` SH basis defined above is the basis used in MRtrix 0.2 Tournier et al.[3]. However, the omission of the \(\sqrt{2}\) factor seen in the basis from Descoteaux et al.[1] makes it non-orthonormal. For this reason, the MRtrix3 Tournier et al.[4] SH basis uses a new basis including the normalization factor.

Since DIPY 1.3, the `descoteaux07` and `tournier07` SH bases have been updated in order to agree with the literature and the latest MRtrix3 implementation. While previous bases are still available as legacy bases, the `descoteaux07` and `tournier07` bases now default to:

\[\begin{split}Y_i(\theta, \phi) = \begin{cases} \sqrt{2} * \Re(Y_l^m(\theta, \phi)) & -l \leq m < 0, \\ Y_l^0(\theta, \phi) & m = 0, \\ \sqrt{2} * \Im(Y_l^m(\theta, \phi)) & 0 < m \leq l \end{cases}\end{split}\]

for the `descoteaux07` basis and

\[\begin{split}Y_i(\theta, \phi) = \begin{cases} \sqrt{2} * \Im(Y_l^{|m|}(\theta, \phi)) & -l \leq m < 0, \\ Y_l^0(\theta, \phi) & m = 0, \\ \sqrt{2} * \Re(Y_l^m(\theta, \phi)) & 0 < m \leq l \end{cases}\end{split}\]

for the `tournier07` basis. Both bases are very similar, with their only difference being the sign of \(m\) for which the imaginary and real parts of the spherical harmonic \(Y_{l}^m\) are used.

In practice, a maximum order \(k\) is used to truncate the SH series. By only taking into account even order SH functions, the above bases can be used to reconstruct symmetric spherical functions. The choice of an even order is motivated by the symmetry of the diffusion process around the origin.

Both bases are also available as full SH bases, where odd order SH functions are also taken into account when reconstructing a spherical function. These full bases can successfully reconstruct asymmetric signals as well as symmetric signals.

NOTE: The definition of spherical harmonics that DIPY utilizes does not match the one in Wikipedia and scipy. Instead, DIPY follows the dMRI literature conventions, like in `descoteaux07` and `tournier07`. The code in DIPY also follows the following convention: Let the SH be noted as \(Y_{l}^m\). Then, \(l\) is referred to as either order or l_value(s), and \(m\) is referred to as either phase factor or m_value(s). These decisions were made as a result of the PR in https://github.com/dipy/dipy/pull/3086

## References

[1] (1,2,3,4)

Maxime Descoteaux, Elaine Angelino, Shaun Fitzgibbons, and Rachid Deriche. Regularized, fast, and robust analytical Q-ball imaging. Magnetic Resonance in Medicine, 58(3):497–510, 2007. URL: https://doi.org/10.1002/mrm.21277, doi:10.1002/mrm.21277.

[2] (1,2,3)

Jacques-Donald Tournier, Fernando Calamante, and Alan Connelly. Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution. NeuroImage, 35(4):1459–1472, May 2007. URL: https://doi.org/10.1016/j.neuroimage.2012.06.005, doi:10.1016/j.neuroimage.2007.02.016.

[3]

Jacques-Donald Tournier, Fernando Calamante, and Alan Connelly. MRtrix: Diffusion tractography in crossing fiber regions. International Journal of Imaging Systems and Technology, 22(1):53–66, 2012. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/ima.22005, doi:https://doi.org/10.1002/ima.22005.

[4]

Jacques-Donald Tournier, Robert Smith, David Raffelt, Rami Tabbara, Thijs Dhollander, Maximilian Pietsch, Daan Christiaens, Ben Jeurissen, Chun-Hung Yeh, and Alan Connelly. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 202:116137, November 2019. URL: https://doi.org/10.1016/j.neuroimage.2019.116137, doi:10.1016/j.neuroimage.2019.116137.