larmor.sd#
Spherical deconvolution: response functions and FOD fitting.
Damped Richardson-Lucy (dRL) deconvolution. |
|
Non-negative least squares with a ridge (L2) penalty on the FOD. |
|
Elastic-net regression of the signal on the kernel's columns (LASSO-type). |
|
Response of a diffusion tensor with a mean kurtosis term. |
|
Isotropic response, exp(-b D), the same along every direction. |
|
The FODs and fractions of an MRIToolkit GRL or mFOD fit. |
|
Response of the NODDI model: dispersed neurites, their surroundings, free water. |
|
Peaks extracted from FODs. |
|
Base class of the response functions that make up a kernel. |
|
Spherical deconvolution of DWI data into FODs. |
|
Result of a spherical deconvolution fit. |
|
Spherical deconvolution methods. |
|
Builds the response functions of an FOD fit. |
|
Converts given coefficients of shape (..., C) from one basis to another. |
|
Estimates the white matter response from the voxels with FA > 0.7. |
|
Extracts the peaks of FODs. |
|
Fits SH coefficients to amplitudes sampled on unit directions. |
|
Reads the struct MRTTrack's PerformDeconv returns ( |
|
The number of FODs (anisotropic responses) a GRL or mFOD result holds. |
|
Reads one FOD image of a MATLAB file in MRIToolkit's layout, as world SH. |
|
The lmax of |
- class larmor.sd.DRL(*, num_iter=200, hr_num_iter=600, target_adc=0.0007, nu_scale=0.2)#
Damped Richardson-Lucy (dRL) deconvolution.
Richardson-Lucy updates the FOD multiplicatively, so it stays non-negative, and sharpens it with every iteration. The damping slows the update of directions whose amplitude is below a threshold nu, set from the kernel (see calibrate_drl_nu), so that the isotropic background does not grow into spurious peaks. A voxel stops early when its FOD no longer changes.
- num_iter#
Most iterations of a single-tissue fit and of each solve of the multi-tissue loop.
- Type:
References
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. doi:10.1016/j.neuroimage.2009.09.033
Follows MRIToolkit’s implementation (MRTTrack).
- class larmor.sd.L2NNLS(*, reg=0.1, tol_x=1e-06)#
Non-negative least squares with a ridge (L2) penalty on the FOD.
Each voxel’s FOD f minimises ||K f - s||^2 + reg^2 ||f||^2 subject to f >= 0, solved by the Lawson-Hanson active-set method on K with reg times the identity stacked under it. The solves of the multi-tissue loop use plain NNLS (reg 0); the final and single-tissue solves use the penalty.
References
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. doi:10.1016/j.neuroimage.2020.117206
Lawson CL, Hanson RJ. Solving Least Squares Problems. Prentice-Hall, 1974.
Follows MRIToolkit’s L2LSQ method (MRTTrack).
- class larmor.sd.Lasso(*, reg=0.1, alpha=0.5, tol=0.0001, max_iter=100000)#
Elastic-net regression of the signal on the kernel’s columns (LASSO-type).
The kernel’s columns are centred and scaled to unit variance and the signal centred; each voxel’s coefficients b minimise (1 / 2m) ||s0 - K0 b||^2 + reg (alpha ||b||_1 + (1 - alpha) / 2 ||b||^2), found by coordinate descent, and the FOD is b divided by the column scales. The L1 part favours sparse FODs. The FOD is not constrained to be non-negative: in a multi-tissue fit the median threshold of the fraction fits keeps negative values out of the fractions, but the FOD itself, as in a single-tissue fit, can keep them. When reg is large enough that no direction explains enough of the signal, the FOD is zero. Every stage of a fit uses the same penalty.
- tol#
Coordinate descent stops when no coefficient changes by more than this, relative to 1 + its size, over a sweep; MATLAB’s lasso default.
- Type:
References
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. doi:10.1111/j.1467-9868.2005.00503.x
Follows MRIToolkit’s LASSO method (MRTTrack), which uses MATLAB’s lasso with Alpha 0.5.
- class larmor.sd.AnisotropicRF(eigen_vals, mean_kurtosis, *, signal_offset=0.0)#
Response of a diffusion tensor with a mean kurtosis term.
The signal along gradient g at b-value b is exp(-b g^T D g + b^2 MD^2 K / 6) + offset, with D the tensor turned so that its first eigenvector points along the fibre, MD its mean eigenvalue and K the mean kurtosis: the diffusion kurtosis model with an isotropic kurtosis tensor.
- eigen_vals#
Tensor eigenvalues in mm^2/s; the first is along the fibre.
- Type:
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. doi:10.1002/mrm.20508
The kernel follows MRIToolkit’s (MRTTrack).
- signal(bvals, bvecs, dirs)#
Signal with the fibre along each direction, shape (n_meas, n_dirs).
- Parameters:
- Returns:
The signal with the fibre along each direction, shape (n_meas, n_dirs).
- Return type:
- class larmor.sd.IsotropicRF(diffusivity)#
Isotropic response, exp(-b D), the same along every direction.
- Parameters:
diffusivity (float) – Isotropic diffusivity D in mm^2/s.
- class larmor.sd.MRIToolkitFOD(fods, tissue_fractions, rss)#
The FODs and fractions of an MRIToolkit GRL or mFOD fit.
- fods#
One FOD per anisotropic response, in the order the responses were added: SH coefficients (X, Y, Z, C) float32 in the MRtrix
tournierbasis, world coordinates, scaled as Larmor’s FODs (300 / (4 pi) times the amplitudes the fit gives on its 300 directions).- Type:
- tissue_fractions#
(X, Y, Z, K), one signal fraction per response, in the order they were added; the file does not name them.
- rss#
(X, Y, Z), the residual sum of squares, or None if not stored.
- Type:
larmor.types.volume.Volume | None
- class larmor.sd.NoddiRF(*, intra_fraction=0.4, parallel_diffusivity=0.0017, concentration=1.0, iso_fraction=0.0, iso_diffusivity=0.003, s0=1.0)#
Response of the NODDI model: dispersed neurites, their surroundings, free water.
The neurites are sticks with diffusivity d along them and none across, their orientations spread around the response’s axis by a Watson distribution of concentration kappa. The space around them is a tensor with the same parallel diffusivity and, by tortuosity, d (1 - f_in) across, averaged over the same distribution. With cos t the cosine between the gradient and the axis:
E_in = 1/2 sum_l L_l(b d) c_l(kappa) sqrt((4l + 1) / 4 pi) P_2l(cos t) E_ex = exp(-b ((d_par - d_perp) cos^2 t + d_perp)) E = S0 ((1 - f_iso) ((1 - f_in) E_ex + f_in E_in) + f_iso exp(-b d_iso))
with c_l the Watson distribution’s SH coefficients, L_l the Legendre-Gaussian integrals and d_par, d_perp the Watson-averaged diffusivities. These are computed with the toolbox’s formulas, accurate to about 1e-7 or better for kappa between 0.3 and 30 and for b d = 0 or above 0.2; near the formulas’ switch points (kappa 0.1, b d 0.05) and above kappa 30 the high-degree terms lose precision.
References
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. doi:10.1016/j.neuroimage.2012.03.072
Follows the NODDI toolbox (SynthMeasWatsonSHStickTortIsoV_B0), as MRIToolkit’s mFOD uses it for grey matter.
- signal(bvals, bvecs, dirs)#
Signal with the axis along each direction, shape (n_meas, n_dirs).
- Parameters:
- Returns:
The signal with the axis along each direction, shape (n_meas, n_dirs).
- Return type:
- class larmor.sd.Peaks(directions, amplitudes)#
Peaks extracted from FODs.
- directions#
Unit vectors along the peaks, shape (…, P, 3) where P is the number of peaks asked for. Sorted by amplitude, largest first. If fewer than P peaks are found, the remaining ones are NaN. A symmetric FOD’s peaks have z >= 0; an asymmetric FOD’s keep their sign, since a peak and its opposite are two lobes.
- Type:
numpy._typing._array_like.NDArray[numpy.float64]
- amplitudes#
Amplitude of each peak, shape (…, P), NaN where directions is NaN.
- Type:
numpy._typing._array_like.NDArray[numpy.float64]
- class larmor.sd.ResponseFunction#
Base class of the response functions that make up a kernel.
- abstractmethod signal(bvals, bvecs, dirs)#
Signal of this response function for each measurement.
- Parameters:
- Returns:
Signal with one column per orientation, shape (n_meas, n_cols).
- Return type:
- class larmor.sd.SDModel(rfs, *, method=SDMethod.GRL, solver=None, inner_shell_weight=None, lmax=8, basis='tournier', progress=False)#
Spherical deconvolution of DWI data into FODs.
Configure the model with its response functions, method and solver, then call
fiton the data. With a single anisotropic response function and no isotropic one, the fit is one FOD solve per voxel. Otherwise it is multi-tissue: per voxel, a low-resolution loop alternates FOD and tissue-fraction fits until the fractions settle, then a high-resolution FOD fit and the final fractions. b-values are in s/mm^2 and the response functions’ diffusivities in mm^2/s; FODs come out in the frame of the DWI’s b-vectors.References
GRL: 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. doi:10.1016/j.neuroimage.2020.116948
mFOD: 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. doi:10.1016/j.neuroimage.2020.117206
The solvers (dRL, L2NNLS, LASSO) cite their own papers: see
larmor.sd.DRL,larmor.sd.L2NNLSandlarmor.sd.Lasso.- fit(volume, gtab, mask=None, *, n_jobs=-1)#
Performs an FOD fit.
- Parameters:
volume (Volume) – Diffusion-weighted volume to fit, shape (X, Y, Z, N).
gtab (GradientTable) – Its gradient table; its
b0_maskselects the b0s the signal is normalised by.mask (Volume | ndarray | None) – Voxels to fit (nonzero inside), on the volume’s grid; a mask volume in another orientation of that grid is reoriented to the volume’s, with a warning. Defaults to None, which fits every voxel.
n_jobs (int) – Number of threads for the fit, at most the number of cores; -1 uses every core. Defaults to -1.
- Returns:
SDFit – The FODs as SH coefficients in the model’s
basis, in the frame of the gradient table (world coordinates), with the tissue fractions of a multi-tissue fit and the lmax actually used. Each voxel’s signal is divided by its S0 fromestimate_s0and capped at 1. Voxels whose signal has a NaN or infinite value, or whose S0 is NaN, are not fitted: their FODs and fractions are NaN.- Raises:
ValueError – If n_jobs is 0 or below -1.
- Return type:
- class larmor.sd.SDFit(fods, tissue_fractions, basis, lmax)#
Result of a spherical deconvolution fit.
(X, Y, Z) is the spatial shape of the DWI, C is the number of spherical harmonic (SH) coefficients up to lmax, and N is the number of response functions. The bases supported so far have even degrees only, for which C = (lmax + 1)(lmax + 2) / 2, e.g. 28 for lmax 6.
- fods#
One FOD per anisotropic response function, in the order given to
SDModel, each as SH coefficients inbasis, shape (X, Y, Z, C). Zero outside the mask.- Type:
- tissue_fractions#
Fraction of each response function per voxel, shape (X, Y, Z, N), summing to 1 inside the mask. The anisotropic response functions come first, then the isotropic ones, each group in the order given to
SDModel. None for a single-tissue fit (one anisotropic response function and no isotropic one).- Type:
numpy.ndarray | None
- basis#
SH basis of the FODs: “tournier” (MRtrix, readable by mrview), “leemans” (ExploreDTI / MRIToolkit) or “descoteaux”.
- Type:
- class larmor.sd.SDMethod#
Spherical deconvolution methods.
GRL fits one anisotropic response function (white matter) and any number of isotropic ones (such as grey matter and CSF); its FOD solver is damped Richardson-Lucy and its inner shells are weighted 0.2. mFOD fits one or more anisotropic response functions, one FOD each, with all shells weighted equally; its FOD solver is L2-regularised NNLS.
References
GRL: 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. doi:10.1016/j.neuroimage.2020.116948
mFOD: 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. doi:10.1016/j.neuroimage.2020.117206
Both follow MRIToolkit’s implementation (MRTTrack).
- GRL = 'GRL'#
- MFOD = 'mFOD'#
- default_solver()#
The FOD solver the method uses when none is given.
- larmor.sd.build_response_functions(volume, gtab, mask, *, wm_eigen_vals=None, mean_kurtosis=None, gm_diffusivity=None, csf_diffusivity=None, gm_noddi=False)#
Builds the response functions of an FOD fit.
They come in the order white matter, grey matter, CSF; grey matter only when its diffusivity is given or gm_noddi is set, CSF only when its diffusivity is given.
- Parameters:
volume (Volume) – Diffusion-weighted volume, shape (X, Y, Z, N).
gtab (GradientTable) – Its gradient table.
mask (Volume | ndarray) – Voxels used to estimate the white matter response (nonzero inside), on the volume’s grid. Only used when wm_eigen_vals is None.
wm_eigen_vals (Sequence[float] | None) – White matter tensor eigenvalues in mm^2/s, the first along the fibre. None estimates the response from the data (estimate_wm_rf).
mean_kurtosis (float | None) – Mean kurtosis of the white matter response. None takes the estimated kurtosis when the response is estimated, and 0 when wm_eigen_vals is given; a number replaces either.
gm_diffusivity (float | None) – Grey matter isotropic diffusivity in mm^2/s; None leaves grey matter out.
csf_diffusivity (float | None) – CSF isotropic diffusivity in mm^2/s; None leaves CSF out.
gm_noddi (bool) – Models grey matter with the anisotropic NODDI response (NoddiRF with its defaults) instead of an isotropic one, as mFOD does.
- Returns:
[white matter, grey matter, CSF], without the ones left out. With any isotropic response, or NODDI grey matter, the fit is multi-tissue.
- Raises:
ValueError – If the white matter response is estimated and no usable voxel in the mask has FA > 0.7, or if both gm_diffusivity and gm_noddi are given.
- Return type:
- larmor.sd.convert_sh_basis(coeffs, from_basis, to_basis)#
Converts given coefficients of shape (…, C) from one basis to another.
- Parameters:
- Returns:
npt.NDArray[np.floating] – coefficients in the new basis.
- Raises:
ValueError – If a basis is unknown, or the number of coefficients is neither a symmetric nor a full SH layout.
- Return type:
NDArray[floating]
- larmor.sd.estimate_wm_rf(volume, gtab, mask)#
Estimates the white matter response from the voxels with FA > 0.7.
In every voxel of the mask, a tensor is fitted to the log signal by iterated weighted least squares, together with an isotropic kurtosis term (b * 1e-3)^2 / 6 when the data has at least two shells above b = 500. Voxels with a signal that is not positive and finite are skipped. Over the voxels with FA > 0.7 (single-fibre white matter), the response’s eigenvalues are the median of each eigenvalue, the second and third then averaged so the response is axially symmetric, and its kurtosis is the mean of the voxel values 1e-6 * p_K / MD^2 that lie in (0, 4).
- Parameters:
volume (Volume) – Diffusion-weighted volume, shape (X, Y, Z, N).
gtab (GradientTable) – Its gradient table.
mask (Volume | ndarray) – Voxels to consider (nonzero inside), on the volume’s grid.
- Returns:
AnisotropicRF – Eigenvalues [l1, l23, l23] in mm^2/s and the mean kurtosis; the kurtosis is 0 with fewer than two shells above b = 500.
- Raises:
ValueError – If the mask does not match the volume, or no usable voxel in it has FA > 0.7.
- Return type:
References
Adapted from MRIToolkit’s Eigenval_IsoK_WM_FromData.
- larmor.sd.find_peaks(coeffs, mask=None, *, basis='tournier', num_peaks=None, relative_threshold=0.1, threshold=0.0, chunk_size=10000, n_jobs=-1)#
Extracts the peaks of FODs.
The maxima of the FOD’s amplitude are found by climbing from the vertices of an icosphere, and peaks closer than 0.99 in
|cos|are merged. Maxima with an amplitude of 0 or below are dropped, as are those underthresholdand those underrelative_thresholdtimes the largest peak.An asymmetric FOD (all SH degrees, recognised by the number of coefficients) is searched the same way, but a peak and its opposite are two peaks, merged only when closer than 0.99 in signed cos, and peaks keep their sign. A straight fibre then gives two peaks, one per lobe.
- Parameters:
coeffs (NDArray[floating]) – FOD coefficients, of shape (…, C): the even degrees of a symmetric FOD, or all degrees of an asymmetric one.
mask (NDArray[bool] | None) – Only computes peaks where mask is True, shape (…); the other voxels are NaN. Defaults to None (all voxels are considered).
basis (str) – Basis in which the FODs are estimated. Defaults to “tournier”.
num_peaks (int | None) – Number of peaks to extract. Defaults to None: 3 for a symmetric FOD, 6 for an asymmetric one, where each fibre is two peaks.
relative_threshold (float) – Threshold relative to the largest peak that considers a peak significant. Defaults to 0.1.
threshold (float) – Threshold under which peaks are not considered. Defaults to 0.0.
chunk_size (int) – Number of voxels whose amplitudes on the icosphere are held in memory at once. Defaults to 10_000.
n_jobs (int) – Number of threads, at most the number of cores, or -1 for every core. Defaults to -1.
- Returns:
Peaks – The peaks extracted.
- Raises:
ValueError – If num_peaks or chunk_size is below 1, if n_jobs is 0 or below -1, if relative_threshold is not between 0 and 1, if mask does not have the shape of coeffs without its last axis, if the coefficient count fits no SH layout or convert_sh_basis refuses the basis, or if lmax is above MAX_LMAX.
- Return type:
- larmor.sd.fit_sh(amplitudes, directions, lmax, *, basis='tournier', full=False)#
Fits SH coefficients to amplitudes sampled on unit directions.
For example, an FOD made from its amplitudes on many directions over the sphere.
- Parameters:
amplitudes (ArrayLike) – (…, N) amplitudes, one per direction along the last axis.
directions (ArrayLike) – (N, 3) unit vectors, in the frame of the coefficients.
lmax (int) – The maximum SH degree.
basis (str) – ‘tournier’ (MRtrix), ‘descoteaux’ or ‘leemans’.
full (bool) – If True, all degrees, for an asymmetric function: (lmax + 1)^2 coefficients. Otherwise the even degrees only: (lmax + 1)(lmax + 2) / 2.
- Returns:
(…, C) least-squares coefficients – evaluated at
directions, they give back the amplitudes as closely as degreelmaxallows.- Raises:
ValueError – If the amplitudes’ last axis does not have one value per direction.
- Return type:
NDArray[float64]
- larmor.sd.read_mritoolkit_fod(source, reference, *, variable=None, lmax=8)#
Reads the struct MRTTrack’s PerformDeconv returns (
FOD,fractions,RSS).The FOD is stored as 300 amplitudes per voxel and anisotropic response, on the 300 directions of
DIRS_300; they are fitted with SH up tolmaxat those directions taken to world coordinates. They are scaled as Larmor’s FODs (300 / (4 pi) times the amplitudes), so they compare with Larmor’s own fits of the same data; MRIToolkit’s exports multiply by its kernel’s scale instead, which the file does not store, so the two differ by one factor.- Parameters:
source (str | Path | MatlabFile) – The
matfile the struct was saved in, or it opened withread_matlab.reference (Grid | str | Path | MatlabFile) – The grid the fit was made on: a Grid, the MRIToolkit DTI or DKI fit (
mat) it came from, or a NIfTI or MIF file of that image. The struct stores no grid.variable (str | None) – The struct’s name in the file. Defaults to the only struct variable with a field
FOD.lmax (int) – The highest SH degree of the FODs (even).
- Returns:
The FODs, the fractions and the residuals on the reference grid.
- Raises:
ValueError – If the file has no such struct, or its arrays do not fit the reference grid.
- Return type:
- larmor.sd.read_mritoolkit_fod_count(source, *, variable=None)#
The number of FODs (anisotropic responses) a GRL or mFOD result holds.
- Parameters:
source (str | Path | MatlabFile) – The
matfile of a GRL or mFOD result, or the file opened withread_matlab.variable (str | None) – The struct holding the FOD field. Defaults to None, the only struct of the file with one.
- Returns:
The count – the FOD field’s last axis divided by the 300 directions of one FOD, at least 1.
- Raises:
ValueError – If
variableis None and the file does not hold exactly one struct with an FOD field.- Return type:
- larmor.sd.read_mritoolkit_sh(source, variable, reference=None, *, lmax=8)#
Reads one FOD image of a MATLAB file in MRIToolkit’s layout, as world SH.
The image’s first three axes are MRIToolkit’s (as
read_mritoolkit_volumereads them), and its directions are in MRIToolkit’s frame (the NIfTI voxel axes j, i, k). Its last axis holds either 300 amplitudes onDIRS_300, fitted as byread_mritoolkit_fod(lmaxlmax, Larmor’s scale), or SH coefficients in the MRtrixtournierbasis (even degrees, or all for an asymmetric FOD), which are evaluated on directions over the whole sphere in that frame and fitted again at those directions taken to world coordinates. A rotation or reflection keeps each degree’s span, so the result is that FOD, exactly up to rounding, in world coordinates.- Parameters:
source (str | Path | MatlabFile) – The
matfile, or it opened withread_matlab.variable (str) – The image, as
MatlabFile.readtakes it.reference (Grid | str | Path | MatlabFile | None) – The grid, as for
read_mritoolkit_volume.lmax (int) – The SH degree of FODs given as amplitudes.
- Returns:
SH coefficients (X, Y, Z, C) float32, MRtrix basis, world coordinates.
- Raises:
ValueError – If the last axis holds neither 300 amplitudes nor an SH layout.
- Return type:
- larmor.sd.sh_layout(num_coeffs)#
The lmax of
num_coeffsSH coefficients, and whether all degrees are stored.A symmetric FOD stores the even degrees, (lmax + 1)(lmax + 2) / 2 coefficients; an asymmetric one all degrees, (lmax + 1)^2. A count that fits both is read as symmetric, which makes 1 coefficient symmetric lmax 0.
- Parameters:
num_coeffs (int) – Number of coefficients.
- Returns:
tuple[int, bool] – lmax, and whether all degrees are stored (asymmetric).
- Raises:
ValueError – If the count fits neither layout.
- Return type: