Peaks#

The peaks of an FOD are its local maxima: the directions of the fibres in a voxel, each with an amplitude. find_peaks() finds them from the FOD’s SH coefficients, largest first, and drops those far smaller than the largest. Conventions gives the SH bases and the layout of the peaks image.

The data on this page

The examples on this page run on a small synthetic FOD image, written to the current folder by the code below: fod.nii.gz, SH coefficients up to degree 8 in the tournier basis, and mask.nii.gz, the voxels of its two bundles. It has 20 x 20 x 4 voxels of 2 mm. One bundle runs along x and one along y, and where they cross the FOD has both lobes. To use your own data, use your files instead.

import numpy as np

from larmor.io import save_volume
from larmor.sd import fit_sh
from larmor.types import Volume

n = 600
z = 1 - (2 * np.arange(n) + 1) / n
azimuth = np.pi * (3 - np.sqrt(5)) * np.arange(n)
radius = np.sqrt(1 - z**2)
sphere = np.stack([radius * np.cos(azimuth), radius * np.sin(azimuth), z], axis=1)


def fod(*fibres):
    """SH coefficients of a narrow lobe along each fibre."""
    cosines = sphere @ np.array(fibres).T
    return fit_sh(np.exp(-8 * (1 - cosines**2)).sum(axis=1), sphere, lmax=8)


x, y = np.eye(3)[0], np.eye(3)[1]
along_x = np.zeros((20, 20, 4), bool)
along_x[:, 8:12] = True
along_y = np.zeros((20, 20, 4), bool)
along_y[8:12, :] = True
coeffs = np.zeros((20, 20, 4, 45), np.float32)
coeffs[along_x] = fod(x)
coeffs[along_y] = fod(y)
coeffs[along_x & along_y] = fod(x, y)

affine = np.diag([2.0, 2.0, 2.0, 1.0])
save_volume(Volume(coeffs, affine), "fod.nii.gz")
save_volume(Volume((along_x | along_y).astype(np.uint8), affine), "mask.nii.gz")

Finding the peaks#

import numpy as np

from larmor.io import load_volume, save_volume
from larmor.sd import find_peaks

fod = load_volume("fod.nii.gz")
mask = load_volume("mask.nii.gz").data > 0
peaks = find_peaks(fod.data, mask=mask, num_peaks=3, relative_threshold=0.1)

crossing = np.isfinite(peaks.amplitudes[10, 10, 0])
print(np.round(peaks.directions[10, 10, 0][crossing], 2))

vectors = peaks.directions * peaks.amplitudes[..., None]
save_volume(fod.with_data(vectors.reshape(*fod.data.shape[:3], -1)), "peaks.nii.gz")

The result, a Peaks, holds the unit directions, shape (X, Y, Z, P, 3), and their amplitudes, shape (X, Y, Z, P), largest first; where a voxel has fewer than P peaks, the others are NaN, and so are the voxels outside the mask. The last lines write the peaks image lar peaks writes.

lar peaks fod.nii.gz peaks.nii.gz --mask mask.nii.gz --num 3 --relative-threshold 0.1

The peaks image has 3 volumes per peak, largest first: x, y and z of the peak’s direction scaled by its amplitude, in world coordinates, NaN where a voxel has fewer peaks.

Options#

  • Number of peaks. 3 by default (num_peaks, --num), or 6 for an asymmetric FOD.

  • Thresholds. A peak below relative_threshold (0.1 by default) times the voxel’s largest peak, or below the absolute threshold (0 by default), is dropped, and so are maxima with an amplitude of 0 or below.

  • Basis. The FOD’s SH basis, "tournier" by default (basis, --basis); see FODs.

  • Symmetric and asymmetric FODs. A symmetric FOD, with even degrees only, has a lobe at both ends of a fibre, so each fibre is one peak, pointing up (z >= 0). An asymmetric FOD, with all degrees, is recognised by its number of coefficients; a peak and its opposite are then two peaks and keep their sign, so a straight fibre gives two.

  • Close fibres. Where two fibres cross at a small angle, their lobes merge into one peak.

Viewing the peaks#

lar bv -p peaks.nii.gz -f fod.nii.gz

The viewer draws each peak as a line through its voxel, over the FOD glyphs.

See also#

References#

  • The search follows sh2peaks of MRtrix3: Tournier JD et al. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage 202 (2019) 116137.