larmor.dwi#
Diffusion-weighted imaging data, brain masks, denoising, and the DTI and DKI fits.
Result of a DKI fit, as maps on the image grid (zero outside the fitted mask). |
|
Diffusion kurtosis model: its settings, and |
|
Result of a DTI fit, as maps on the image grid (zero outside the fitted mask). |
|
Diffusion tensor model: its settings, and |
|
Result of |
|
Compute a brain mask from a b0 volume. |
|
Denoises a diffusion-weighted volume by MP-PCA and estimates its noise level. |
|
Returns the signal without diffusion weighting, S0, of every voxel of a mask. |
|
Returns the mean of a DWI's b0 volumes. |
|
Removes Gibbs ringing from a 3D or 4D volume, slice by slice. |
|
Select the volumes whose b-value is within |
- class larmor.dwi.DKIFit(dt, kt, s0, fa, md, ad, rd, pe, eigen_vals, mk, ak, rk, ka)#
Result of a DKI fit, as maps on the image grid (zero outside the fitted mask).
(X, Y, Z) is the spatial shape of the DWI. Diffusivities are in the inverse units of the b-values: mm^2/s for b-values in s/mm^2. Every map is float32.
The apparent kurtosis along a unit vector n is K(n) = MD^2 W(n) / D(n)^2, with D(n) = n^T D n and W(n) the quartic form of the kurtosis tensor.
- dt#
Diffusion tensor, shape (X, Y, Z, 6): its 6 unique elements in the order xx, xy, xz, yy, yz, zz, in mm^2/s.
- Type:
- kt#
Kurtosis tensor W, shape (X, Y, Z, 15): its 15 unique elements in the order W1111, W2222, W3333, W1112, W1113, W2221, W2223, W3331, W3332, W1122, W1133, W2233, W1123, W2213, W3312; dimensionless.
- Type:
- s0#
Fitted signal at b = 0, shape (X, Y, Z), in the units of the DWI.
- Type:
- fa#
Fractional anisotropy, shape (X, Y, Z), unitless, in [0, 1].
- Type:
- md#
Mean diffusivity (l1 + l2 + l3) / 3, shape (X, Y, Z), in mm^2/s.
- Type:
- ad#
Axial diffusivity l1, shape (X, Y, Z), in mm^2/s.
- Type:
- rd#
Radial diffusivity (l2 + l3) / 2, shape (X, Y, Z), in mm^2/s.
- Type:
- pe#
Principal eigenvector, the unit vector of l1, shape (X, Y, Z, 3), in the frame of the b-vectors. Its sign is arbitrary.
- Type:
- eigen_vals#
Eigenvalues l1 >= l2 >= l3, shape (X, Y, Z, 3), in mm^2/s. Negative eigenvalues are set to 0.
- Type:
- mk#
Mean kurtosis, shape (X, Y, Z); dimensionless.
- Type:
- ak#
Axial kurtosis, shape (X, Y, Z); dimensionless.
- Type:
- rk#
Radial kurtosis, shape (X, Y, Z); dimensionless.
- Type:
- ka#
Kurtosis anisotropy, the standard deviation of K(n) over the sphere, shape (X, Y, Z); dimensionless.
- Type:
- class larmor.dwi.DKIModel(*, fit_method='wls', constrained=True)#
Diffusion kurtosis model: its settings, and
fitto fit it to a DWI.The signal model of Jensen et al. (2005), log(S) = log(S0) - b g^T D g + b^2 MD^2 W(g) / 6, is linear in log(S0), the 6 elements of D and the 15 elements of MD^2 W, and is fitted to the log signal of each voxel of the mask by ordinary least squares (“ols”) or by weighted least squares with weights from the predicted signal (“wls”, see
larmor.dwi.fitting.fit_log_signal). W is then the fitted MD^2 W over MD^2.References
Jensen JH, Helpern JA, Ramani A, Lu H, Kaczynski K. Diffusional kurtosis imaging: the quantification of non-Gaussian water diffusion by means of magnetic resonance imaging. Magn Reson Med 53 (2005) 1432-1440. Veraart J, Sijbers J, Sunaert S, Leemans A, Jeurissen B. Weighted linear least squares estimation of diffusion MRI parameters: strengths, limitations, and pitfalls. NeuroImage 81 (2013) 335-346.
- fit(volume, gtab, mask=None, *, n_jobs=-1)#
Fits the diffusion and kurtosis tensors in every voxel of the mask.
Before the fit, each voxel’s S0 is that of
larmor.dwi.estimate_s0: where it differs from the mean of the voxel’s own b0s (implausible b0s), it replaces the b0 measurements, and every diffusion-weighted measurement is capped at it. A voxel without a finite, positive S0 is not fitted and is 0 in every map, as outside the mask (a warning gives their number). A measurement that is still not finite and positive is set to 1e-3 S0 so that its log is defined. Without a b0 in the table there is no S0 to compare with: nothing is replaced or capped, and such measurements are set to 1e-3 times the largest value of the data. MK, AK, RK and KA are those ofcompute_kurtosis_metrics.The constrained fit (Tabesh et al., 2011) minimises the same (weighted) sum of squares subject to D(n) >= 0, K(n) >= 0 and K(n) <= 3 / (b_max D(n)) along 256 directions n of a half sphere (the constraints are even in n), b_max the largest b-value of the table, all shells fitted, as in the paper. With the fitted MD^2 W, V(n) = MD^2 W(n), and K(n) = V(n) / D(n)^2, these are linear: D(n) >= 0, V(n) >= 0 and 3 D(n) / b_max - V(n) >= 0 (the last multiplied by D(n)^2). The problem is convex, so a voxel whose unconstrained solution meets every constraint keeps it; the others are solved as a least-squares problem with linear inequality constraints, reduced to a non-negative least-squares problem (Lawson and Hanson, 1974). For “wls”, the weights are those of the signal the unconstrained weighted fit predicts, with no further reweighting. Where the constrained solve finds no solution or does not converge, the voxel keeps its unconstrained fit. The bound grows as D(n) goes to 0, so voxels with near-zero diffusivity (at the edge of the mask, for example) can give very large kurtosis values. Voxels are fitted in parallel, each on its own, so the result does not depend on n_jobs.
- Parameters:
volume (Volume) – Diffusion-weighted volume, shape (X, Y, Z, N).
gtab (GradientTable) – Its gradient table, with at least two non-zero shells; the tensors come out in its frame, world coordinates.
mask (Volume | ndarray | None) – Voxels to fit (nonzero inside), on the volume’s grid (a mask volume in another orientation is refused). Defaults to None, which computes a brain mask from the mean b0 volume.
n_jobs (int) – Number of threads for the fit, at most the number of cores; -1 uses every core. Defaults to -1.
- Returns:
DKIFit – The tensors and their metrics as maps of shape (X, Y, Z, …), zero outside the mask.
- Raises:
ValueError – If the table or the mask does not match the volume, the measurements cannot determine the tensors (fewer than two non-zero shells, or too few directions), if the mask is computed and there is no b0, or if n_jobs is 0 or below -1.
- Return type:
References
Tabesh A, Jensen JH, Ardekani BA, Helpern JA. Estimation of tensors and tensor-derived measures in diffusional kurtosis imaging. Magn Reson Med 65 (2011) 823-836. Lawson CL, Hanson RJ. Solving Least Squares Problems. Prentice-Hall, 1974 (chapter 23).
- class larmor.dwi.DTIFit(dt, fa, md, ad, rd, pe, eigen_vals)#
Result of a DTI fit, as maps on the image grid (zero outside the fitted mask).
(X, Y, Z) is the spatial shape of the DWI. Diffusivities are in the inverse units of the b-values: mm^2/s for b-values in s/mm^2. Every map is float32.
- dt#
Diffusion tensor, shape (X, Y, Z, 6): its 6 unique elements in the order xx, xy, xz, yy, yz, zz, in mm^2/s.
- Type:
- fa#
Fractional anisotropy, shape (X, Y, Z), unitless, in [0, 1].
- Type:
- md#
Mean diffusivity (l1 + l2 + l3) / 3, shape (X, Y, Z), in mm^2/s.
- Type:
- ad#
Axial diffusivity l1, shape (X, Y, Z), in mm^2/s.
- Type:
- rd#
Radial diffusivity (l2 + l3) / 2, shape (X, Y, Z), in mm^2/s.
- Type:
- pe#
Principal eigenvector, the unit vector of l1, shape (X, Y, Z, 3), in the frame of the b-vectors. Its sign is arbitrary.
- Type:
- eigen_vals#
Eigenvalues l1 >= l2 >= l3, shape (X, Y, Z, 3), in mm^2/s. Negative eigenvalues are set to 0.
- Type:
- class larmor.dwi.DTIModel(*, fit_method='ols')#
Diffusion tensor model: its settings, and
fitto fit it to a DWI.The signal model log(S) = log(S0) - b g^T D g is linear in log(S0) and the 6 elements of D, and is fitted to the log signal of each voxel of the mask by ordinary least squares (“ols”) or by weighted least squares with weights from the predicted signal (“wls”, see
larmor.dwi.fitting.fit_log_signal).References
Basser PJ, Mattiello J, LeBihan D. Estimation of the effective self-diffusion tensor from the NMR spin echo. J Magn Reson B 103 (1994) 247-254. Veraart J, Sijbers J, Sunaert S, Leemans A, Jeurissen B. Weighted linear least squares estimation of diffusion MRI parameters: strengths, limitations, and pitfalls. NeuroImage 81 (2013) 335-346.
- fit(volume, gtab, mask=None, *, n_jobs=-1)#
Fits the tensor in every voxel of the mask.
Voxels are fitted in parallel, each on its own, so the result does not depend on n_jobs.
- Parameters:
volume (Volume) – Diffusion-weighted volume, shape (X, Y, Z, N).
gtab (GradientTable) – Its gradient table, with at least two distinct b-values (a b0 is not needed when a mask is given); the tensor comes out in its frame, world coordinates.
mask (Volume | ndarray | None) – Voxels to fit (nonzero inside), on the volume’s grid (a mask volume in another orientation is refused). Defaults to None, which computes a brain mask from the mean b0 volume.
n_jobs (int) – Number of threads for the fit, at most the number of cores; -1 uses every core. Defaults to -1.
- Returns:
DTIFit – The tensor and its metrics as maps of shape (X, Y, Z, …), zero outside the mask. Signal values that are not positive or finite are replaced by 1 before taking the log.
- Raises:
ValueError – If the table or the mask does not match the volume, the measurements cannot determine the tensor and S0 (for example one shell without b0, or fewer than six non-coplanar directions), if the mask is computed and there is no b0, or if n_jobs is 0 or below -1.
- Return type:
- class larmor.dwi.DenoisedDWI(volume, noise, num_signal)#
Result of
denoise_mppca: the denoised volume and the noise map.- volume#
The denoised volume, float32, on the input’s grid; zero outside the mask.
- noise#
The noise level sigma of each voxel, in the units of the signal, shape (X, Y, Z), float32; zero outside the mask.
- Type:
- num_signal#
The number of signal components P kept in each voxel’s patch, shape (X, Y, Z), int64; zero outside the mask.
- Type:
- larmor.dwi.create_brain_mask(data_b0, *, threshold=0.5, morph=7)#
Compute a brain mask from a b0 volume.
- Parameters:
- Returns:
Boolean mask, shape (X, Y, Z).
- Return type:
NDArray[bool]
- larmor.dwi.denoise_mppca(volume, mask=None, *, patch_size=None, n_jobs=-1, progress=False)#
Denoises a diffusion-weighted volume by MP-PCA and estimates its noise level.
Every voxel is denoised in its own patch, a box of patch_size voxels centred on it; near the edge of the image the box is shifted to stay inside it. The patch’s signals form a matrix X with one row per volume (M) and one column per voxel (N), with no mean removed. Pure noise of level sigma gives eigenvalues of X X^T / max(M, N) that follow the Marchenko-Pastur distribution, of mean sigma^2 and width 4 sqrt(gamma) sigma^2, with gamma the ratio of the number of noise components to max(M, N). The number of signal components P is the smallest p for which the mean of the smallest min(M, N) - p eigenvalues is at least the sigma^2 their range implies, and sigma^2 is that mean (Veraart et al. 2016, equations 10-12). The voxel’s signal is projected onto the P principal components; only the centre voxel of each patch is kept, so the result does not depend on the number of threads. Patches with fewer voxels than volumes are handled with the roles of M and N swapped.
The noise is assumed Gaussian, constant within a patch and uncorrelated: denoise first, before any interpolation, and note that the Rician bias of magnitude images is not corrected. A patch holding a NaN or infinite value gives NaN for its voxel.
- Parameters:
volume (Volume) – The diffusion-weighted volume, (X, Y, Z, M), with all its volumes (b0s included), as acquired.
mask (Volume | ndarray | None) – Voxels to denoise (nonzero inside), on the volume’s grid. Patches still use the voxels outside it. Defaults to None, every voxel.
patch_size (int | integer | Sequence[int | integer] | None) – The patch width in voxels, odd (Python or NumPy integers), the same along every axis or one per axis; along an axis shorter than that, the whole axis. Defaults to None, the smallest odd width w with w^3 >= M (5 for up to 125 volumes, 7 for up to 343).
n_jobs (int) – Number of threads, at most the cores; -1 is all cores. Defaults to -1.
progress (bool) – Whether to show a progress bar. Defaults to False.
- Returns:
DenoisedDWI – The denoised volume, the noise map sigma and the number of signal components per voxel.
- Raises:
ValueError – If the volume is not 4D, the mask does not match its grid, patch_size is not odd and positive, or n_jobs is 0 or below -1.
TypeError – If patch_size is not an integer (Python or NumPy) or a sequence of them.
- Return type:
References
Veraart J, Novikov DS, Christiaens D, Ades-aron B, Sijbers J, Fieremans E. Denoising of diffusion MRI using random matrix theory. NeuroImage 142 (2016) 394-406. Tournier JD et al. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage 202 (2019) 116137 (
dwidenoise -estimator Exp1, the reference in the tests).
- larmor.dwi.estimate_s0(volume, gtab, mask)#
Returns the signal without diffusion weighting, S0, of every voxel of a mask.
S0 is the mean of a voxel’s b0 measurements. Diffusion weighting can only lower the signal, so a voxel whose mean b0 is at or below the mean of its diffusion-weighted measurements has implausible b0s (dark, or 0, as at the edge of the brain). Such a voxel takes the mean S0 of the plausible mask voxels in its 3x3x3 neighbourhood, or in its 5x5x5 neighbourhood when there is none, and NaN when there is none there either. Voxels with a NaN or infinite value keep their mean b0 and never serve as neighbours.
- Parameters:
volume (Volume) – Diffusion-weighted volume, shape (X, Y, Z, N).
gtab (GradientTable) – Its gradient table; its
b0_maskselects the b0s.mask (ndarray) – Voxels to estimate S0 in, shape (X, Y, Z) bool.
- Returns:
S0, shape (X, Y, Z), float64; the mean b0 outside the mask.
- Raises:
ValueError – If the table does not match the volume, or has no b0.
- Return type:
- larmor.dwi.mean_b0(volume, gtab)#
Returns the mean of a DWI’s b0 volumes.
- Parameters:
volume (Volume) – Diffusion-weighted volume, shape (X, Y, Z, N).
gtab (GradientTable) – Its gradient table; its
b0_maskselects the b0s.
- Returns:
The mean b0 volume, shape (X, Y, Z), float32.
- Raises:
ValueError – If the table does not match the volume, or has no b0.
- Return type:
- larmor.dwi.remove_gibbs_ringing(volume, *, slice_axis=2, num_shifts=20, window=(1, 3), n_jobs=-1)#
Removes Gibbs ringing from a 3D or 4D volume, slice by slice.
Truncating k-space at the edge of the acquired frequencies makes sharp edges, such as the border between CSF and tissue, ring: the image oscillates next to them. The local sub-voxel shift method (Kellner et al., 2016) re-interpolates each line of a slice at sub-voxel shifts of j / (2 num_shifts) voxels,
|j| <= num_shifts, by multiplying its spectrum with a phase. For each voxel it keeps the shift whose oscillation next to the voxel is smallest, measured as the total variation between the neighbours window[0] to window[1] + 1 voxels away, on the left and on the right (an edge can be on either side), and interpolates that shifted line linearly back onto the voxel. Each slice is first split in k-space into two images: the one unringed along the first in-plane axis keeps the frequencies (k_1, k_2) (radians per voxel) with weight (1 + cos k_2) / (2 + cos k_1 + cos k_2), the one unringed along the second with (1 + cos k_1) / (2 + cos k_1 + cos k_2), and the two results are added. Lines wrap around at the edge of the slice, as the Fourier transform does.Run it on the images as the scanner reconstructed them, before any interpolation (motion or distortion correction), and after denoising: Gibbs removal changes the noise, which denoising relies on. The method assumes full k-space sampling; along a partial Fourier direction some ringing may remain. Each slice is computed on its own, so the result does not depend on the number of threads. A slice with a NaN or infinite value comes out all NaN.
- Parameters:
volume (Volume) – The image, (X, Y, Z) or (X, Y, Z, N); its scaling is applied.
slice_axis (int) – The voxel axis across the acquired slices (0, 1 or 2); the other two span each slice. Defaults to 2.
num_shifts (int) – Shifts tried on each side of zero. Defaults to 20.
window (tuple[int, int]) – The neighbours, in voxels from the voxel, over which the oscillation is measured, (first, last). Defaults to (1, 3).
n_jobs (int) – Number of threads, -1 for all cores. Defaults to -1.
- Returns:
Volume – The image without ringing, float32, on the same grid.
- Raises:
ValueError – If slice_axis is not 0, 1 or 2, num_shifts is below 1, the window is not 0 <= first <= last, or n_jobs is 0 or below -1.
- Return type:
References
Kellner E, Dhital B, Kiselev VG, Reisert M. Gibbs-ringing artifact removal based on local subvoxel-shifts. Magn Reson Med 76 (2016) 1574-1581. Tournier JD et al. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage 202 (2019) 116137.
- larmor.dwi.shell_mask(bvals, shell, *, tol=100.0)#
Select the volumes whose b-value is within
tolofshell(s/mm^2).