Gradient tables and shells#

A GradientTable holds the diffusion weighting of each volume of a diffusion image, in the order of the volumes: a b-value in s/mm^2 and a unit direction in world coordinates (RAS+). On construction, each b-value is scaled by the squared length of its vector and each vector is made a unit vector, so a vector of length 0.5 at b=2800 becomes a unit vector at b=700.

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")

Reading and writing#

FSL bvecs files store the directions in the frame of the image’s voxel axes, with x flipped when the affine’s determinant is positive; Larmor calls this frame "fsl", the default. Reading or writing them therefore needs the image they belong to. Files that hold the directions in world coordinates are read with frame="world". MRtrix gradient files (.b, one row x y z b per volume) are always in world coordinates.

from larmor.io import read_bvals_bvecs, read_mrtrix_gradients, write_mrtrix_gradients

gtab = read_bvals_bvecs("dwi.bval", "dwi.bvec", ref_path="dwi.nii.gz")
write_mrtrix_gradients(gtab, "dwi.b")
same = read_mrtrix_gradients("dwi.b")

The commands read bvals and bvecs files with --bvals and --bvecs, and --bvec-frame world for b-vectors in world coordinates. lar convert writes the table into the header of a MIF image, or to new bvals and bvecs files:

lar convert dwi.nii.gz dwi.mif --bvals dwi.bval --bvecs dwi.bvec
lar convert dwi.mif dwi_copy.nii.gz --bvals-out dwi_copy.bval --bvecs-out dwi_copy.bvec

write_bvals_bvecs(), read_b_matrix() and write_b_matrix() read and write the other formats; File formats lists them.

Shells#

Volumes at or below b0_threshold (50 s/mm^2 by default) count as b=0. The other b-values are sorted, and a gap of more than shell_epsilon (80 s/mm^2 by default) between two of them starts a new shell, so b-values spread around a shell (995, 1000, 1005) count as one.

print(gtab.shells)
print(gtab.count_shells(), gtab.b0_mask.sum())
print(gtab.shell_index[:5])

shells gives the b-value of each shell, the mean of its volumes, and shell_index the shell of each volume, -1 for a b=0 volume.

lar info prints the number of b=0 volumes and each shell, with its number of volumes, from a MIF image’s table:

lar info dwi.mif

Keeping one shell#

shell_mask() selects the volumes within a tolerance of a b-value (100 s/mm^2 by default). To keep the b=0 volumes and the b=1000 shell, for a DTI fit on that shell alone:

import numpy as np

from larmor.dwi import shell_mask
from larmor.io import load_dwi, save_volume
from larmor.types import GradientTable

dwi, gtab = load_dwi("dwi.nii.gz", "dwi.bval", "dwi.bvec")
keep = gtab.b0_mask | shell_mask(gtab.bvals, 1000)
b1000 = dwi.with_data(dwi.data[..., keep])
b1000_gtab = GradientTable(gtab.bvals[keep], gtab.bvecs[keep])
save_volume(b1000, "dwi_b1000.mif", gradients=b1000_gtab)

Correcting a table#

When the directions do not fit the data, for example with the sign of y reversed by the scanner’s export, flip() negates one component of every direction and permute() reorders the components. Reorienting an image never needs them, since the table is in world coordinates. The direction colours of a DTI fit (see DTI) show whether a table is right: the corpus callosum runs left to right.

See also#