FODs#

Spherical deconvolution estimates a fibre orientation distribution (FOD) in each voxel: how much of the voxel’s fibres runs in each direction. It models the signal as the sum of the signals of single fibres, each the white matter response turned to point along the fibre, and finds the distribution that gives back the measured signal. FODs are stored as spherical harmonic (SH) coefficients. Citing lists the papers of the methods, and Conventions gives the SH bases.

The data on this page

The examples on this page run on a small synthetic image, written to the current folder by the code below: dwi.nii.gz with its gradient files dwi.bval and dwi.bvec, and a mask, mask.nii.gz. It has 12 x 12 x 2 voxels of 2 mm, two b=0 images and 30 directions at each of b=1000 and 2000 s/mm^2. Two bundles cross in its middle, one along x and one along y; the other voxels hold isotropic diffusion. To use your own data, use your files instead.

import numpy as np

from larmor.io import save_volume, write_bvals_bvecs
from larmor.sd import AnisotropicRF
from larmor.types import GradientTable, Volume

rng = np.random.default_rng(0)
directions = rng.normal(size=(60, 3))
directions /= np.linalg.norm(directions, axis=1, keepdims=True)
bvals = np.r_[0, 0, [1000] * 30, [2000] * 30]
bvecs = np.r_[np.zeros((2, 3)), directions]
gtab = GradientTable(bvals, bvecs)

fibre = AnisotropicRF([1.7e-3, 0.3e-3, 0.3e-3], mean_kurtosis=0.5)
along_x, along_y = fibre.signal(bvals, gtab.bvecs, np.eye(3)[:2]).T
isotropic = np.exp(-bvals * 0.8e-3)

signal = np.broadcast_to(isotropic, (12, 12, 2, len(bvals))).copy()
signal[:, 4:8] = along_x
signal[4:8, :] = along_y
signal[4:8, 4:8] = (along_x + along_y) / 2
data = 1000 * signal + rng.normal(0, 10, signal.shape)

affine = np.diag([2.0, 2.0, 2.0, 1.0])
save_volume(Volume(data.astype(np.float32), affine), "dwi.nii.gz")
write_bvals_bvecs(gtab, "dwi.bval", "dwi.bvec", ref_path="dwi.nii.gz")
save_volume(Volume(np.ones((12, 12, 2), np.uint8), affine), "mask.nii.gz")

Fitting FODs#

The fit needs response functions: the white matter response, estimated from the voxels with FA > 0.7 unless its eigenvalues are given, and optionally isotropic responses for grey matter and CSF. With isotropic responses the fit is multi-tissue and also gives each tissue’s fraction per voxel.

from larmor.io import load_dwi, load_volume, save_volume
from larmor.sd import SDModel, build_response_functions

dwi, gtab = load_dwi("dwi.nii.gz", "dwi.bval", "dwi.bvec")
mask = load_volume("mask.nii.gz")
responses = build_response_functions(
    dwi, gtab, mask, gm_diffusivity=0.7e-3, csf_diffusivity=3e-3
)
fit = SDModel(responses, lmax=8).fit(dwi, gtab, mask=mask)

save_volume(dwi.with_data(fit.fods[0]), "fod.nii.gz")
save_volume(dwi.with_data(fit.tissue_fractions), "fractions.nii.gz")
print(f"lmax {fit.lmax}, basis {fit.basis}, shape {fit.fods[0].shape}")

build_response_functions returns the responses in the order white matter, grey matter, CSF; without gm_diffusivity and csf_diffusivity it returns the white matter response alone, for a single-tissue fit. wm_eigen_vals=[1.7e-3, 0.2e-3, 0.2e-3] gives the white matter response instead of estimating it.

lar fod dwi.nii.gz mask.nii.gz fod --bvals dwi.bval --bvecs dwi.bvec

The command writes the FOD to fod.nii.gz and the tissue fractions to fod_tissue_fraction.nii.gz. By default it fits white matter, grey matter (isotropic, 0.7e-3 mm^2/s, with GRL) and CSF (3e-3 mm^2/s); --gm-rf none and --csf-rf none leave them out, and --wm-rf 1.7e-3 0.2e-3 0.2e-3 gives the white matter eigenvalues instead of estimating them.

A fit with even degrees up to lmax has (lmax + 1)(lmax + 2) / 2 coefficients: 28 for degree 6, 45 for degree 8. The degree is at most lmax (8 by default), lowered until the coefficients are no more than the directions of the highest shell: these data, with 30 directions per shell, give degree 6. The FODs are in world coordinates.

Methods and solvers#

Method

Responses

Default solver

Weight of the inner shells in a multi-tissue fit

GRL, the default (Guo et al. 2020)

one anisotropic (white matter), any isotropic

damped Richardson-Lucy, DRL (Dell’Acqua et al. 2010)

0.2

mFOD (De Luca et al. 2020)

one or more anisotropic, one FOD each, and any isotropic

non-negative least squares with a ridge penalty, L2NNLS

1

Each method can use any solver: DRL, L2NNLS (reg 0.1 by default) or Lasso, an elastic net (reg 0.1, alpha 0.5; Tibshirani 1996; Zou and Hastie 2005). mFOD can model grey matter with the anisotropic NODDI response (NoddiRF; Zhang et al. 2012), which gives a second FOD:

from larmor.sd import L2NNLS, SDMethod

responses = build_response_functions(
    dwi, gtab, mask, gm_noddi=True, csf_diffusivity=3e-3
)
model = SDModel(responses, method=SDMethod.MFOD, solver=L2NNLS(reg=0.1))
fit = model.fit(dwi, gtab, mask=mask)
print(len(fit.fods), "FODs")

lar fod dwi.nii.gz mask.nii.gz mfod --bvals dwi.bval --bvecs dwi.bvec \
    --method mfod

With mFOD, grey matter defaults to the NODDI response (--gm-rf noddi) and the solver to L2NNLS. This writes one FOD per anisotropic response, mfod_fod1.nii.gz (white matter) and mfod_fod2.nii.gz (grey matter), and mfod_tissue_fraction.nii.gz.

inner_shell_weight (--inner-shell-weight) changes the weight of the measurements outside the highest shell in a multi-tissue fit.

SH bases#

Larmor writes and reads three real SH bases, with the same functions in a different order, sign or frame: "tournier", the default, as MRtrix3 uses it; "descoteaux"; and "leemans", as ExploreDTI and MRIToolkit use it. An image file does not record the basis, so it must be given again whenever an FOD is read: basis= in Python, --basis on the command line (lar fod, lar peaks, lar track). The fit writes its FODs in another basis with SDModel(..., basis="leemans") or lar fod --basis leemans, and existing FODs are converted with:

from larmor.io import load_volume, save_volume
from larmor.sd import convert_sh_basis, sh_layout

fod = load_volume("fod.nii.gz")
leemans = convert_sh_basis(fod.data, "tournier", "leemans")
save_volume(fod.with_data(leemans), "fod_leemans.nii.gz")
print(sh_layout(fod.data.shape[-1]))

sh_layout() tells the degree from the number of coefficients, and whether all degrees are stored, as for an asymmetric FOD.

lar sh-convert fod.nii.gz fod_leemans.nii.gz --from tournier --to leemans

See also#

References#

  • Guo F, Leemans A, Viergever MA, Dell’Acqua F, De Luca A. Generalized Richardson-Lucy (GRL) for analyzing multi-shell diffusion MRI data. NeuroImage 218 (2020) 116948.

  • De Luca A, Guo F, Froeling M, Leemans A. Spherical deconvolution with tissue-specific response functions and multi-shell diffusion MRI to estimate multiple fiber orientation distributions (mFODs). NeuroImage 222 (2020) 117206.

  • Dell’Acqua F, Scifo P, Rizzo G, et al. A modified damped Richardson-Lucy algorithm to reduce isotropic background effects in spherical deconvolution. NeuroImage 49 (2010) 1446-1458.

  • Tibshirani R. Regression shrinkage and selection via the lasso. J R Stat Soc B 58 (1996) 267-288.

  • Zou H, Hastie T. Regularization and variable selection via the elastic net. J R Stat Soc B 67 (2005) 301-320.

  • Zhang H, Schneider T, Wheeler-Kingshott CA, Alexander DC. NODDI: practical in vivo neurite orientation dispersion and density imaging of the human brain. NeuroImage 61 (2012) 1000-1016.

  • The methods and solvers are reimplementations of those in MRIToolkit (De Luca, github.com/delucaal/MRIToolkit).