Base-classes for reconstruction models and reconstruction fits.
All the models in the reconst module follow the same template: a Model object
is used to represent the abstract properties of the model, that are independent
of the specifics of the data . These properties are reused whenever fitting a
particular set of data (different voxels, for example).
Tools for fitting Bingham distributions to orientation distribution
functions (ODF), as described in [1]. The resulting
distributions can further be used to compute ODF-lobe-specific measures such as
the fiber density (FD) and fiber spread (FS) [1] and the
orientation dispersion index (ODI) [2].
Extract the diffusion tensor eigenvalues, the diffusion tensor eigenvector matrix, and the 15 independent elements of the kurtosis tensor from the model parameters estimated from the DKI model
Selector function to switch between the 2-stage Trust-Region Reflective based NLLS fitting method (also containing the linear fit): trr and the Variable Projections based fitting method: varpro.
The q:math:tau-dMRI model to analytically and continuously represent the q:math:tau diffusion signal attenuation over diffusion sensitization q and diffusion time \(\tau\).
This is an implementation of the sparse fascicle model described in
Rokem et al.[4]. The multi b-value version of this model is described
in Rokem et al.[5].
Note about the Transpose:
In the literature the matrix representation of these methods is often written
as Y = Bx where B is some design matrix and Y and x are column vectors. In our
case the input data, a dwi stored as a nifti file for example, is stored as row
vectors (ndarrays) of the form (x, y, z, n), where n is the number of diffusion
directions. We could transpose and reshape the data to be (n, x*y*z), so that
we could directly plug it into the above equation. However, I have chosen to
keep the data as is and implement the relevant equations rewritten in the
following form: Y.T = x.T B.T, or in python syntax data = np.dot(sh_coef, B.T)
where data is Y.T and sh_coef is x.T.
Overall Orientation Dispersion Index (ODI) for an ODF lobe.
Overall Orientation Dispersion Index (ODI) computed for an
ODF lobe from the overall concentration parameter (k_total).
Defined by equation 20 in [10].
Reconstruct ODFs from fitted Bingham parameters on multiple voxels.
Parameters:
bingham_paramsndarray (…., nl, 12)
ndarray containing the model parameters of Binghams fitted to ODFs in
the following order:
- Maximum value of the Bingham function (f0, index 0);
- concentration parameters k1 and k2 (indexes 1 and 2);
- elements of Bingham’s main direction (indexes 3-5);
- elements of Bingham’s dispersion major axis (indexes 6-8);
- elements of Bingham’s dispersion minor axis (indexes 9-11).
sphere: `Sphere` class instance
The Sphere providing the odf’s discrete directions
mask: ndarray, optional
Map marking the coordinates in the data that should be analyzed.
Default (None) means all voxels in the volume will be analyzed.
Compute fiber density for each lobe for a given Bingham ODF.
Measured in the unit 1/mm^3.
Parameters:
bingham_paramsndarray (…., nl, 12)
ndarray containing the model parameters of Bingham’s fitted to ODFs in
the following order:
- Maximum value of the Bingham function (f0, index 0);
- Concentration parameters k1 and k2 (indexes 1 and 2);
- Elements of Bingham’s main direction (indexes 3-5);
- Elements of Bingham’s dispersion major axis (indexes 6-8);
- Elements of Bingham’s dispersion minor axis (indexes 9-11).
subdivide: int >= 0, optional
Number of times the unit icosahedron used for integration
should be subdivided. The higher this value the more precise the
approximation will be, at the cost of longer execution times. The
default results in a sphere of 10242 points.
mask: ndarray, optional
Map marking the coordinates in the data that should be analyzed.
Default (None) means all voxels in the volume will be analyzed.
Returns:
fd: ndarray (…., nl)
Fiber density for each Bingham function.
Notes
Fiber density (fd) is given by the surface integral of the
Bingham function [1].
Compute density-weighted scalar maps for metrics of Bingham functions
fitted to multiple ODF lobes. The metric is computed as the
weighted average of a given metric across the multiple ODF lobes
(weights are defined by the Bingham fiber density estimates).
Parameters:
bmetric: ndarray(…, nl)
Any metric with values for nl ODF lobes.
bfd: ndarray(…, nl)
Bingham’s fiber density estimates for the nl ODF lobes
The coefficient of determination is calculated as:
\[R^2 = 100 * (1 - \frac{SSE}{SSD})\]
where SSE is the sum of the squared error between the model and the data
(sum of the squared residuals) and SSD is the sum of the squares of the
deviations of the data from the mean of the data (variance * N).
The type of the model to use for prediction. The corresponding Fit
object must have a predict function implemented One of the following:
reconst.dti.TensorModel or
reconst.csdeconv.ConstrainedSphericalDeconvModel.
datandarray
Diffusion MRI data acquired with the GradientTable of the model. Shape
will typically be (x, y, z, b) where xyz are spatial dimensions and
b is the number of bvals/bvecs in the GradientTable.
foldsint
The number of divisions to apply to the data
model_argslist
Additional arguments to the model initialization
model_kwargsdict
Additional key-word arguments to the model initialization. If contains
the kwarg mask, this will be used as a key-word argument to the fit
method of the model object, rather than being used in the
initialization of the model object
Notes
This function assumes that a prediction API is implemented in the Model
class for which prediction is conducted. That is, the Fit object that gets
generated upon fitting the model needs to have a predict method, which
receives a GradientTable class instance as input and produces a predicted
signal as output.
It also assumes that the model object has bval and bvec attributes
holding b-values and corresponding unit vectors.
ratio = \(\frac{\lambda_2}{\lambda_1}\) of the single fiber response
function
l_valuesndarray (N,)
The order (\(l\)) of spherical harmonic function associated with each row
of the deconvolution matrix. Only even orders are allowed.
r2_termbool
True if ODF comes from an ODF computed from a model using the \(r^2\)
term in the integral. For example, DSI, GQI, SHORE, CSA, Tensor,
Multi-tensor ODFs. This results in using the proper analytical response
function solution solving from the single-fiber ODF with the r^2 term.
This derivation is not published anywhere but is very similar to
[12].
Deconvolves the axially symmetric single fiber response function r_rh in
rotational harmonics coefficients from the diffusion weighted signal in
dwsignal [13].
Parameters:
dwsignalarray
Diffusion weighted signals to be deconvolved.
Xarray
Prediction matrix which estimates diffusion weighted signals from FOD
coefficients.
B_regarray (N, B)
SH basis matrix which maps FOD coefficients to FOD values on the
surface of the sphere. B_reg should be scaled to account for lambda.
taufloat
Threshold controlling the amplitude below which the corresponding fODF
is assumed to be zero. Ideally, tau should be set to zero. However, to
improve the stability of the algorithm, tau is set to tau*100 % of the
max fODF amplitude (here, 10% by default). This is similar to peak
detection where peaks below 0.1 amplitude are usually considered noise
peaks. Because SDT is based on a q-ball ODF deconvolution, and not
signal deconvolution, using the max instead of mean (as in CSD), is
more stable.
convergenceint
Maximum number of iterations to allow the deconvolution to converge.
Pndarray
This is an optimization to avoid computing dot(X.T,X) many times.
If the same X is used many times, P can be precomputed and
passed to this function.
Where \(X\) maps current FOD SH coefficients \(f_n\) to DW signals \(s\) and
\(H_{n-1}\) maps FOD SH coefficients \(f_n\) to amplitudes along set of
negative directions identified in previous iteration, i.e. the matrix
formed by the rows of \(B_{reg}\) for which \(Hf_{n-1}<0\) where \(B_{reg}\)
maps \(f_n\) to FOD amplitude on a sphere.
Define \(Q = X^TX + \lambda^2 H_{n-1}^TH_{n-1}\) , which by construction is a
square positive definite symmetric matrix of size \(n_{SH} by n_{SH}\). If
needed, positive definiteness can be enforced with a small minimum norm
regulariser (helps a lot with poorly conditioned direction sets and/or
superresolution):
ndarray of SH coefficients for the ODF spherical function to be
deconvolved
Rndarray ((sh_order_max+1)(sh_order_max+2)/2,
(sh_order_max+1)(sh_order_max+2)/2)
SDT matrix in SH basis
B_regndarray ((sh_order_max+1)(sh_order_max+2)/2,
(sh_order_max+1)(sh_order_max+2)/2)
SH basis matrix used for deconvolution
lambda_float, optional
lambda parameter in minimization equation
taufloat, optional
threshold (tau*max(fODF)) controlling the amplitude below
which the corresponding fODF is assumed to be zero.
r2_termbool, optional
True if ODF is computed from model that uses the \(r^2\) term in the
integral. Recall that Tuch’s ODF (used in Q-ball Imaging
[14]) and the true normalized ODF definition differ
from a \(r^2\) term in the ODF integral. The original Sharpening
Deconvolution Transform (SDT) technique [15]
is expecting Tuch’s ODF without the \(r^2\) (see
[12] for the mathematical details). Now, this
function supports ODF that have been computed using the \(r^2\) term
because the proper analytical response function has be derived. For
example, models such as DSI, GQI, SHORE, CSA, Tensor, Multi-tensor ODFs,
should now be deconvolved with the r2_term=True.
Sharpen odfs using the sharpening deconvolution transform.
This function can be used to sharpen any smooth ODF spherical function. In
theory, this should only be used to sharpen QballModel ODFs, but in
practice, one can play with the deconvolution ratio and sharpen almost any
ODF-like spherical function. The constrained-regularization is stable and
will not only sharpen the ODF peaks but also regularize the noisy peaks.
different spherical harmonic basis:
None for the default DIPY basis,
tournier07 for the Tournier 2007 [13] basis,
and descoteaux07 for the Descoteaux 2007
[7] basis (None defaults to
descoteaux07).
ratiofloat, optional
ratio of the smallest vs the largest eigenvalue of the single prolate
tensor response function (\(\frac{\lambda_2}{\lambda_1}\))
sh_order_maxint, optional
maximal SH order (\(l\)) of the SH representation
lambda_float, optional
lambda parameter (see odfdeconv)
taufloat, optional
tau parameter in the L matrix construction (see odfdeconv)
r2_termbool, optional
True if ODF is computed from model that uses the \(r^2\) term in the
integral. Recall that Tuch’s ODF (used in Q-ball Imaging
[14]) and the true normalized ODF definition differ
from a \(r^2\) term in the ODF integral. The original Sharpening
Deconvolution Transform (SDT) technique [15] is
expecting Tuch’s ODF without the \(r^2\) (see [7]
for the mathematical details). Now, this function supports ODF that
have been computed using the \(r^2\) term because the proper analytical
response function has be derived. For example, models such as DSI,
GQI, SHORE, CSA, Tensor, Multi-tensor ODFs, should now be deconvolved
with the r2_term=True.
Returns:
fodf_shndarray
sharpened odf expressed as spherical harmonics coefficients
Computation of mask for single-shell single-tissue (ssst) response
function using FA.
Parameters:
gtabGradientTable
Gradient table.
datandarray
diffusion data (4D)
roi_centerarray-like, (3,)
Center of ROI in data. If center is None, it is assumed that it is
the center of the volume with shape data.shape[:3].
roi_radiiint or array-like, (3,)
radii of cuboid ROI
fa_thrfloat
FA threshold
Returns:
maskndarray
Mask of voxels within the ROI and with FA above the FA threshold.
Notes
In CSD there is an important pre-processing step: the estimation of the
fiber response function. In order to do this, we look for voxels with very
anisotropic configurations. This function aims to accomplish that by
returning a mask of voxels within a ROI, that have a FA value above a
given threshold. For example we can use a ROI (20x20x20) at
the center of the volume and store the signal values for the voxels with
FA values higher than 0.7 (see [16]).
Computation of single-shell single-tissue (ssst) response
function from a given mask.
Parameters:
gtabGradientTable
Gradient table.
datandarray
diffusion data
maskndarray
mask from where to compute the response function
Returns:
responsetuple, (2,)
(evals, S0)
ratiofloat
The ratio between smallest versus largest eigenvalue of the response.
Notes
In CSD there is an important pre-processing step: the estimation of the
fiber response function. In order to do this, we look for voxels with very
anisotropic configurations. This information can be obtained by using
csdeconv.mask_for_response_ssst() through a mask of selected voxels
(see [16]). The present function uses such a mask to
compute the ssst response function.
For the response we also need to find the average S0 in the ROI. This is
possible using gtab.b0s_mask() we can find all the S0 volumes (which
correspond to b-values equal 0) in the dataset.
The response consists always of a prolate tensor created by averaging
the highest and second highest eigenvalues in the ROI with FA higher than
threshold. We also include the average S0s.
We also return the ratio which is used for the SDT models.
Automatic estimation of single-shell single-tissue (ssst) response
function using FA.
Parameters:
gtabGradientTable
Gradient table.
datandarray
diffusion data
roi_centerarray-like, (3,)
Center of ROI in data. If center is None, it is assumed that it is
the center of the volume with shape data.shape[:3].
roi_radiiint or array-like, (3,)
radii of cuboid ROI
fa_thrfloat
FA threshold
Returns:
responsetuple, (2,)
(evals, S0)
ratiofloat
The ratio between smallest versus largest eigenvalue of the response.
Notes
In CSD there is an important pre-processing step: the estimation of the
fiber response function. In order to do this, we look for voxels with very
anisotropic configurations. We get this information from
csdeconv.mask_for_response_ssst(), which returns a mask of selected voxels
(more details are available in the description of the function).
With the mask, we compute the response function by using
csdeconv.response_from_mask_ssst(), which returns the response and the
ratio (more details are available in the description of the function).
mask for recursive calibration, for example a white matter mask. It has
shape data.shape[0:3] and dtype=bool. Default: use the entire data
array.
sh_order_maxint, optional
maximal spherical harmonics order (l).
peak_thrfloat, optional
peak threshold, how large the second peak can be relative to the first
peak in order to call it a single fiber population
[17].
init_fafloat, optional
FA of the initial ‘fat’ response function (tensor).
init_tracefloat, optional
trace of the initial ‘fat’ response function (tensor).
iterint, optional
maximum number of iterations for calibration.
convergencefloat, optional
convergence criterion, maximum relative change of SH
coefficients.
parallelbool, optional
Whether to use parallelization in peak-finding during the calibration
procedure.
num_processesint, optional
If parallel is True, the number of subprocesses to use
(default multiprocessing.cpu_count()). If < 0 the maximal number of
cores minus num_processes+1 is used (enter -1 to use as many
cores as possible). 0 raises an error.
sphereSphere, optional.
The sphere used for peak finding.
Returns:
responsendarray
response function in SH coefficients
Notes
In CSD there is an important pre-processing step: the estimation of the
fiber response function. Using an FA threshold is not a very robust method.
It is dependent on the dataset (non-informed used subjectivity), and still
depends on the diffusion tensor (FA and first eigenvector),
which has low accuracy at high b-value. This function recursively
calibrates the response function, for more information see
[17].
where \(K_{aniso}\) is the anisotropic kurtosis,
\(\langle V_{\lambda}(D_c) \rangle\) represents the mean of the variance
of eigenvalues of the diffusion tensor, \(\overline{D}\) is the mean of
the diffusion tensor.
where: \(K_{\text{iso}}\) is the isotropic kurtosis, \(V({\overline{D}^c})\)
represents the variance of the diffusion tensor raised to the power \(c\),
\(\overline{D}\) is the mean of the diffusion tensor.
where \(\Psi\) is a variable representing a part of the total
excess kurtosis, \(D_{ij}\) are elements of the diffusion tensor,
\(\overline{D}\) is the mean of the diffusion tensor. \(\overline{W}\)
is the mean kurtosis, \(W_{ijkl}\) are elements of the kurtosis
tensor.
Extract the diffusion tensor eigenvalues, the diffusion tensor
eigenvector matrix, and the 21 independent elements of the covariance
tensor, and the 15 independent elements of the kurtosis tensor from the
model parameters estimated from the CTI model
Parameters:
cti_params: numpy.ndarray (…, 48)
All parameters estimated from the correlation tensor model.
Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the
first, second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
Twenty-One elements of the covariance tensor
Returns:
evalsarray (…, 3)
Eigenvalues from eigen decomposition of the tensor.
evecsarray (…, 3)
Associated eigenvectors from eigen decomposition of the tensor.
Eigenvectors are columnar (e.g. evecs[:,j] is associated with
evals[j])
Compute the diffusion kurtosis and covariance tensors using an
ordinary or weighted linear least squares approach
Parameters:
design_matrixarray (g, 43)
Design matrix holding the covariants used to solve for the regression
coefficients.
dataarray (g)
Data or response variables holding the data.
inverse_design_matrixarray (43, g)
Inverse of the design matrix.
weightsbool, optional
Parameter indicating whether weights are used.
min_diffusivityfloat, optional
Because negative eigenvalues are not physical and small eigenvalues,
much smaller than the diffusion weighting, cause quite a lot of noise
in metrics such as fa, diffusivity values smaller than min_diffusivity
are replaced with min_diffusivity.
Returns:
cti_paramsarray (48)
All parameters estimated from the diffusion kurtosis model for all N
voxels. Parameters are ordered as follows:
Three diffusion tensor eigenvalues.
Three blocks of three elements, containing the first second and
third coordinates of the diffusion tensor eigenvectors.
A boolean array used to mark the coordinates in the data that
should be analyzed that has the shape data.shape[-1]
num_iterint, optional
Number of times to iterate.
weights_methodcallable, optional
A function with args and returns as follows:
(weights, robust) =
weights_method(data, pred_sig,
design_matrix, leverages,
idx, num_iter,
robust)
Notes
On the first iteration, (C)WLS fit is performed. Subsequent iterations
will be weighted according to weights_method. Outlier rejection should
be handled within weights_method by setting the corresponding weights
to zero (if weights_method implements outlier rejection).
Compute axial kurtosis (AK) of a diffusion kurtosis tensor.
See [21] and [20] for
further details about the method.
Parameters:
min_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range, axial
kurtosis values that are smaller than min_kurtosis are replaced
with -3./7 (theoretical kurtosis limit
for regions that consist of water confined to spherical pores
[19]).
max_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range, axial
kurtosis values that are larger than max_kurtosis are replaced
with max_kurtosis.
analyticalbool, optional
If True, AK is calculated from rotated diffusion kurtosis tensor,
otherwise it will be computed from the apparent diffusion kurtosis
values along the principal axis of the diffusion tensor
(see notes).
Returns:
akarray
Calculated AK.
Notes
AK is defined as the directional kurtosis parallel to the fiber’s main
direction e1 [21], [20].
You can compute AK using to approaches:
AK is calculated from rotated diffusion kurtosis tensor, i.e.:
AK can be sampled from the principal axis of the diffusion tensor:
\[AK = K(\mathbf{e}_1)\]
Although both approaches leads to an exact calculation of AK, the
first approach will be referred to as the analytical method while the
second approach will be referred to as the numerical method based on
their analogy to the estimation strategies for MK and RK.
where \(W\) is the kurtosis tensor, MKT the kurtosis tensor mean, \(I^{(4)}\)
is the fully symmetric rank 2 isotropic tensor and \(||...||_F\) is the
tensor’s Frobenius norm [22].
The sphere providing sample directions for the initial search of
the maximum value of kurtosis.
gtolfloat, optional
This input is to refine kurtosis maximum under the precision of the
directions sampled on the sphere class instance. The gradient of
the convergence procedure must be less than gtol before successful
termination. If gtol is None, fiber direction is directly taken
from the initial sampled directions of the given sphere object
Compute mean kurtosis (MK) from the kurtosis tensor.
See Tabesh et al.[21] and [20] for
further details about the method.
Parameters:
min_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range, mean
kurtosis values that are smaller than min_kurtosis are replaced
with min_kurtosis. Default = -3./7 (theoretical kurtosis limit
for regions that consist of water confined to spherical pores
[19]).
max_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range, mean
kurtosis values that are larger than max_kurtosis are replaced
with max_kurtosis.
analyticalbool, optional
If True, MK is calculated using its analytical solution, otherwise
an exact numerical estimator is used (see Notes).
Returns:
mkarray
Calculated MK.
Notes
The MK is defined as the average of directional kurtosis coefficients
across all spatial directions, which can be formulated by the following
surface integral [21],
[20]:
where \(\hat{W}_{ijkl}\) are the components of the \(W\) tensor in the
coordinates system defined by the eigenvectors of the diffusion tensor
\(\mathbf{D}\) and
To keep kurtosis values within a plausible biophysical range, mean
kurtosis values that are smaller than min_kurtosis are replaced
with min_kurtosis. Default = -3./7 (theoretical kurtosis limit
for regions that consist of water confined to spherical pores
[19]).
max_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range, mean
kurtosis values that are larger than max_kurtosis are replaced
with max_kurtosis.
where \(D_{ij}\) and \(W_{ijkl}\) are the elements of the second-order DT
and the fourth-order KT tensors, respectively, and \(MD\) is the mean
diffusivity.
To keep kurtosis values within a plausible biophysical range,
radial kurtosis values that are smaller than min_kurtosis are
replaced with min_kurtosis. Default = -3./7 (theoretical kurtosis
limit for regions that consist of water confined to spherical pores
[19]).
max_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range,
radial kurtosis values that are larger than max_kurtosis are
replaced with max_kurtosis.
analyticalbool, optional
If True, RK is calculated using its analytical solution, otherwise
an exact numerical estimator is used (see Notes).
Returns:
rkarray
Calculated RK.
Notes
RK is defined as the average of the directional kurtosis perpendicular
to the fiber’s main direction e1 [21],
[20]:
To keep kurtosis values within a plausible biophysical range,
radial kurtosis values that are smaller than min_kurtosis are
replaced with min_kurtosis.
max_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range,
radial kurtosis values that are larger than max_kurtosis are
replaced with max_kurtosis.
Returns:
rtkarray
Calculated escaled radial tensor kurtosis (RTK).
Notes
Rescaled radial tensor kurtosis (RTK) is defined as
[25]:
where W is the kurtosis tensor rotated to a coordinate system in
which the 3 orthonormal eigenvectors of DT are the base coordinate,
MD is the mean diffusivity, and RD is the radial diffusivity.
Calculate the apparent diffusion coefficient (adc) in each direction of
a sphere for a single voxel Neto Henriques et al.[27].
Parameters:
dtarray (6,)
elements of the diffusion tensor of the voxel.
Varray (g, 3)
g directions of a Sphere in Cartesian coordinates
min_diffusivityfloat, optional
Because negative eigenvalues are not physical and small eigenvalues
cause quite a lot of noise in diffusion-based metrics, diffusivity
values smaller than min_diffusivity are replaced with
min_diffusivity. Default = 0
Returns:
adcndarray (g,)
Apparent diffusion coefficient (ADC) in all g directions of a sphere
for a single voxel.
Calculate the apparent kurtosis coefficient (akc) in each direction of
a sphere for a single voxel.
See [27] and [20] for
further details about the method.
Parameters:
dtarray (6,)
elements of the diffusion tensor of the voxel.
mdfloat
mean diffusivity of the voxel
ktarray (15,)
elements of the kurtosis tensor of the voxel.
Varray (g, 3)
g directions of a Sphere in Cartesian coordinates
min_diffusivityfloat, optional
Because negative eigenvalues are not physical and small eigenvalues
cause quite a lot of noise in diffusion-based metrics, diffusivity
values smaller than min_diffusivity are replaced with
min_diffusivity. Default = 0
min_kurtosisfloat, optional
Because high-amplitude negative values of kurtosis are not physically
and biologicaly pluasible, and these cause artefacts in
kurtosis-based measures, directional kurtosis values smaller than
min_kurtosis are replaced with min_kurtosis.
(theoretical kurtosis limit for regions that consist of water confined
to spherical pores [19]).
adcndarray(g,), optional
Apparent diffusion coefficient (ADC) in all g directions of a sphere
for a single voxel.
advndarray(g,), optional
Apparent diffusion variance (advc) in all g directions of a sphere for
a single voxel.
Returns:
akcndarray (g,)
Apparent kurtosis coefficient (AKC) in all g directions of a sphere for
a single voxel.
Calculate the apparent kurtosis coefficient (AKC) in each direction
of a sphere.
See [27] and [20] for
further details about the method.
Parameters:
dki_paramsndarray (x, y, z, 27) or (n, 27)
All parameters estimated from the diffusion kurtosis model.
Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the first,
second and third coordinates of the eigenvectors respectively
Fifteen elements of the kurtosis tensor
spherea Sphere class instance
The AKC will be calculated for each of the vertices in the sphere
min_diffusivityfloat, optional
Because negative eigenvalues are not physical and small eigenvalues
cause quite a lot of noise in diffusion-based metrics, diffusivity
values smaller than min_diffusivity are replaced with
min_diffusivity.
min_kurtosisfloat, optional
Because high-amplitude negative values of kurtosis are not physically
and biologicaly pluasible, and these cause artefacts in
kurtosis-based measures, directional kurtosis values smaller than
min_kurtosis are replaced with min_kurtosis. Default = -3./7
(theoretical kurtosis limit for regions that consist of water confined
to spherical pores [19]).
Returns:
akcndarray (x, y, z, g) or (n, g)
Apparent kurtosis coefficient (AKC) for all g directions of a sphere.
Notes
For each sphere direction with coordinates \((n_{1}, n_{2}, n_{3})\), the
calculation of AKC is done using formula [27]:
Compute mean kurtosis (MK) from the kurtosis tensor.
See [21] and [20] for
further details about the method.
Parameters:
dki_paramsndarray (x, y, z, 27) or (n, 27)
All parameters estimated from the diffusion kurtosis model.
Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the first,
second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
min_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range, mean
kurtosis values that are smaller than min_kurtosis are replaced with
min_kurtosis. Default = -3./7 (theoretical kurtosis limit for regions
that consist of water confined to spherical pores
[19]).
max_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range, mean
kurtosis values that are larger than max_kurtosis are replaced with
max_kurtosis.
analyticalbool, optional
If True, MK is calculated using its analytical solution, otherwise an
exact numerical estimator is used (see Notes).
Returns:
mkarray
Calculated MK.
Notes
The MK is defined as the average of directional kurtosis coefficients
across all spatial directions, which can be formulated by the following
surface integral [21], [20]:
where \(\hat{W}_{ijkl}\) are the components of the \(W\) tensor in the
coordinates system defined by the eigenvectors of the diffusion tensor
\(\mathbf{D}\) and
Compute radial kurtosis (RK) of a diffusion kurtosis tensor.
See [21] and [20] for
further details about the method.
Parameters:
dki_paramsndarray (x, y, z, 27) or (n, 27)
All parameters estimated from the diffusion kurtosis model.
Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the first,
second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
min_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range, radial
kurtosis values that are smaller than min_kurtosis are replaced with
min_kurtosis. Default = -3./7 (theoretical kurtosis limit for regions
that consist of water confined to spherical pores
[19]).
max_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range, radial
kurtosis values that are larger than max_kurtosis are replaced with
max_kurtosis.
analyticalbool, optional
If True, RK is calculated using its analytical solution, otherwise an
exact numerical estimator is used (see Notes). Default is set to True.
Returns:
rkarray
Calculated RK.
Notes
RK is defined as the average of the directional kurtosis perpendicular
to the fiber’s main direction e1 [21],
[20]:
Compute axial kurtosis (AK) from the kurtosis tensor.
See [21] and [20] for further
details about the method.
Parameters:
dki_paramsndarray (x, y, z, 27) or (n, 27)
All parameters estimated from the diffusion kurtosis model.
Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the first,
second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
min_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range, axial
kurtosis values that are smaller than min_kurtosis are replaced with
min_kurtosis. Default = -3./7 (theoretical kurtosis limit for regions
that consist of water confined to spherical pores
[19]).
max_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range, axial
kurtosis values that are larger than max_kurtosis are replaced with
max_kurtosis.
analyticalbool, optional
If True, AK is calculated from rotated diffusion kurtosis tensor,
otherwise it will be computed from the apparent diffusion kurtosis
values along the principal axis of the diffusion tensor (see notes).
Returns:
akarray
Calculated AK.
Notes
AK is defined as the directional kurtosis parallel to the fiber’s main
direction e1 [21], [20]. You
can compute AK using to approaches:
AK is calculated from rotated diffusion kurtosis tensor
[21], i.e.:
AK can be sampled from the principal axis of the diffusion tensor
[20]:
\[AK = K(\mathbf{e}_1)\]
Although both approaches leads to an exact calculation of AK, the first
approach will be referred to as the analytical method while the second
approach will be referred to as the numerical method based on their analogy
to the estimation strategies for MK and RK.
All parameters estimated from the diffusion kurtosis model.
Parameters are ordered as follows:
Three diffusion tensor’s eingenvalues
Three lines of the eigenvector matrix each containing the first,
second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
sphereSphere class instance, optional
The sphere providing sample directions for the initial search of the
maximal value of kurtosis.
gtolfloat, optional
This input is to refine kurtosis maximum under the precision of the
directions sampled on the sphere class instance. The gradient of the
convergence procedure must be less than gtol before successful
termination. If gtol is None, fiber direction is directly taken from
the initial sampled directions of the given sphere object
maskndarray
A boolean array used to mark the coordinates in the data that should be
analyzed that has the shape dki_params.shape[:-1]
Returns:
max_valuefloat
kurtosis tensor maximum value
max_dirarray (3,)
Cartesian coordinates of the direction of the maximal kurtosis value
All parameters estimated from the diffusion kurtosis model.
Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the first,
second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
min_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range, mean
kurtosis values that are smaller than min_kurtosis are replaced with
min_kurtosis. Default = -3./7 (theoretical kurtosis limit for regions
that consist of water confined to spherical pores
[19]).
max_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range, mean
kurtosis values that are larger than max_kurtosis are replaced with
max_kurtosis.
All parameters estimated from the diffusion kurtosis model.
Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the first,
second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
min_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range, radial
kurtosis values that are smaller than min_kurtosis are replaced with
min_kurtosis.
max_kurtosisfloat, optional
To keep kurtosis values within a plausible biophysical range, radial
kurtosis values that are larger than max_kurtosis are replaced with
max_kurtosis.
Returns:
rtkarray
Calculated rescaled radial tensor kurtosis (RTK).
Notes
Rescaled radial tensor kurtosis (RTK) is defined as
[24]:
where W is the kurtosis tensor rotated to a coordinate system in which the
3 orthonormal eigenvectors of DT are the base coordinate, MD is the mean
diffusivity, and RD is the radial diffusivity.
where \(W\) is the kurtosis tensor, MKT the kurtosis tensor mean, \(I^{(4)}\) is
the fully symmetric rank 2 isotropic tensor and \(||...||_F\) is the tensor’s
Frobenius norm [22].
Convert the 21 unique elements of the diffusion and kurtosis tensors
to the parameter format adopted in DIPY
Parameters:
resultsarray (21)
Unique elements of the diffusion and kurtosis tensors in the following
order: 1) six unique lower triangular DT elements; and 2) Fifteen
unique elements of the kurtosis tensor.
min_diffusivityfloat, optional
Because negative eigenvalues are not physical and small eigenvalues,
much smaller than the diffusion weighting, cause quite a lot of noise
in metrics such as fa, diffusivity values smaller than
min_diffusivity are replaced with min_diffusivity.
Returns:
dki_paramsarray (27)
All parameters estimated from the diffusion kurtosis model for all N
voxels. Parameters are ordered as follows:
Three diffusion tensor eigenvalues.
Three blocks of three elements, containing the first second and
third coordinates of the diffusion tensor eigenvectors.
Design matrix holding the covariants used to solve for the regression
coefficients.
dataarray (g)
Data or response variables holding the data.
inverse_design_matrixarray (22, g)
Inverse of the design matrix.
return_S0_hatbool, optional
Boolean to return (True) or not (False) the S0 values for the fit.
weightsarray ([X, Y, Z, …], g), optional
Weights to apply for fitting. These weights must correspond to the
squared residuals such that \(S = \sum_i w_i r_i^2\). If not provided,
weights are estimated as the squared predicted signal from an initial
OLS fit.
min_diffusivityfloat, optional
Because negative eigenvalues are not physical and small eigenvalues,
much smaller than the diffusion weighting, cause quite a lot of noise
in metrics such as fa, diffusivity values smaller than
min_diffusivity are replaced with min_diffusivity.
return_lower_triangularbool, optional
Boolean to return (True) or not (False) the coefficients of the fit.
return_leveragesbool, optional
Boolean to return (True) or not (False) the fitting leverages.
Returns:
dki_paramsarray (27)
All parameters estimated from the diffusion kurtosis model for all N
voxels. Parameters are ordered as follows:
Three diffusion tensor eigenvalues.
Three blocks of three elements, containing the first second and
third coordinates of the diffusion tensor eigenvectors.
Fifteen elements of the kurtosis tensor.
leveragesarray (g)
Leverages of the fitting problem (if return_leverages is True)
Design matrix holding the covariants used to solve for the regression
coefficients.
dataarray (g)
Data or response variables holding the data.
inverse_design_matrixarray (22, g)
Inverse of the design matrix.
sdpPositiveDefiniteLeastSquares instance
A CVXPY representation of a regularized least squares optimization
problem.
return_S0_hatbool, optional
Boolean to return (True) or not (False) the S0 values for the fit.
weightsarray ([X, Y, Z, …], g), optional
Weights to apply for fitting. These weights must correspond to the
squared residuals such that \(S = \sum_i w_i r_i^2\). If not provided,
weights are estimated as the squared predicted signal from an initial
OLS fit.
weightsbool, optional
Parameter indicating whether weights are used.
min_diffusivityfloat, optional
Because negative eigenvalues are not physical and small eigenvalues,
much smaller than the diffusion weighting, cause quite a lot of noise
in metrics such as fa, diffusivity values smaller than
min_diffusivity are replaced with min_diffusivity.
return_lower_triangularbool, optional
Boolean to return (True) or not (False) the coefficients of the fit.
return_leveragesbool, optional
Boolean to return (True) or not (False) the fitting leverages.
cvxpy_solverstr, optional
cvxpy solver name. Optionally optimize the positivity constraint with a
particular cvxpy solver. See https://www.cvxpy.org/ for details.
Default: None (cvxpy chooses its own solver).
Returns:
dki_paramsarray (27)
All parameters estimated from the diffusion kurtosis model for all N
voxels. Parameters are ordered as follows:
Three diffusion tensor eigenvalues.
Three blocks of three elements, containing the first second and
third coordinates of the diffusion tensor eigenvectors.
Fifteen elements of the kurtosis tensor.
leveragesarray (g)
Leverages of the fitting problem (if return_leverages is True)
Vector with the 15 independent elements of the kurtosis tensor
Basisarray (3, 3)
Vectors of the basis column-wise oriented
indsarray(m, 4), optional
Array of vectors containing the four indexes of m specific elements of
the rotated kurtosis tensor. If not specified all 15 elements of the
rotated kurtosis tensor are computed.
Returns:
Wrotarray (m,) or (15,)
Vector with the m independent elements of the rotated kurtosis tensor.
If ‘indices’ is not specified all 15 elements of the rotated kurtosis
tensor are computed.
Notes
The kurtosis tensor elements are assumed to be ordered as follows:
Extract the diffusion tensor eigenvalues, the diffusion tensor
eigenvector matrix, and the 15 independent elements of the kurtosis tensor
from the model parameters estimated from the DKI model
Parameters:
dki_paramsndarray (x, y, z, 27) or (n, 27)
All parameters estimated from the diffusion kurtosis model.
Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the first,
second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
Returns:
eigvalsarray (x, y, z, 3) or (n, 3)
Eigenvalues from eigen decomposition of the tensor.
eigvecsarray (x, y, z, 3, 3) or (n, 3, 3)
Associated eigenvectors from eigen decomposition of the tensor.
Eigenvectors are columnar (e.g. eigvecs[:,j] is associated with
eigvals[j])
Fit method of the Diffusion Kurtosis Microstructural Model
Parameters:
dataarray
An 4D matrix containing the diffusion-weighted data.
maskarray
A boolean array used to mark the coordinates in the data that
should be analyzed that has the shape data.shape[-1]
sphereSphere class instance, optional
The sphere providing sample directions for the initial search of
the maximal value of kurtosis.
gtolfloat, optional
This input is to refine kurtosis maxima under the precision of the
directions sampled on the sphere class instance. The gradient of
the convergence procedure must be less than gtol before successful
termination. If gtol is None, fiber direction is directly taken
from the initial sampled directions of the given sphere object
awf_onlybool, optiomal
If set to true only the axonal volume fraction is computed from
the kurtosis tensor. Default = False
Predict a signal for the DKI microstructural model class instance
given parameters.
Parameters:
paramsndarray (x, y, z, 40) or (n, 40)
All parameters estimated from the diffusion kurtosis
microstructural model.
Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the
first, second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
Six elements of the hindered diffusion tensor
Six elements of the restricted diffusion tensor
Axonal water fraction
S0float or ndarray, optional
The non diffusion-weighted signal in every voxel, or across all
voxels.
Notes
In the original article of DKI microstructural model
[31], the hindered and restricted tensors were
defined as the intra-cellular and extra-cellular diffusion compartments
respectively.
Returns the tortuosity of the hindered diffusion which is defined by ADe / RDe, where ADe and RDe are the axial and radial diffusivities of the hindered compartment.
Methods
ad()
Axial diffusivity (AD) calculated from cached eigenvalues.
adc(sphere)
Calculate the apparent diffusion coefficient (ADC) in each direction on the sphere for each voxel in the data
ak(*[, min_kurtosis, max_kurtosis, analytical])
Compute axial kurtosis (AK) of a diffusion kurtosis tensor.
akc(sphere)
Calculate the apparent kurtosis coefficient (AKC) in each direction on the sphere for each voxel in the data
color_fa()
Color fractional anisotropy of diffusion tensor
fa()
Fractional anisotropy (FA) calculated from cached eigenvalues.
ga()
Geodesic anisotropy (GA) calculated from cached eigenvalues.
kmax(*[, sphere, gtol, mask])
Compute the maximum value of a single voxel kurtosis tensor
linearity()
md()
Mean diffusivity (MD) calculated from cached eigenvalues.
mk(*[, min_kurtosis, max_kurtosis, analytical])
Compute mean kurtosis (MK) from the kurtosis tensor.
mkt(*[, min_kurtosis, max_kurtosis])
Compute mean of the kurtosis tensor (MKT).
mode()
Tensor mode calculated from cached eigenvalues.
odf(sphere)
The diffusion orientation distribution function (dODF).
where \(ADC_{r}\) and \(ADC_{h}\) are the apparent diffusion coefficients
of the diffusion hindered and restricted compartment for a given
direction \(\theta\), \(b\) is the b value provided in the GradientTable
input for that direction, \(f\) is the volume fraction of the restricted
diffusion compartment (also known as the axonal water fraction).
Returns the tortuosity of the hindered diffusion which is defined
by ADe / RDe, where ADe and RDe are the axial and radial diffusivities
of the hindered compartment.
All parameters estimated from the diffusion kurtosis model.
Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the first,
second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
sphereSphere class instance, optional
The sphere providing sample directions for the initial search of the
maximal value of kurtosis.
gtolfloat, optional
This input is to refine kurtosis maxima under the precision of the
directions sampled on the sphere class instance. The gradient of the
convergence procedure must be less than gtol before successful
termination. If gtol is None, fiber direction is directly taken from
the initial sampled directions of the given sphere object
maskndarray
A boolean array used to mark the coordinates in the data that should be
analyzed that has the shape dki_params.shape[:-1]
All parameters estimated from the diffusion kurtosis model.
Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the first,
second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
sphereSphere class instance, optional
The sphere providing sample directions to sample the restricted and
hindered cellular diffusion tensors. For more details see
[31].
awfndarray, optional
Array containing values of the axonal water fraction that has the shape
dki_params.shape[:-1]. If not given this will be automatically computed
using axonal_water_fraction() with function’s default precision.
maskndarray, optional
A boolean array used to mark the coordinates in the data that should be
analyzed that has the shape dki_params.shape[:-1]
Returns:
edtndarray (x, y, z, 6) or (n, 6)
Parameters of the hindered diffusion tensor.
idtndarray (x, y, z, 6) or (n, 6)
Parameters of the restricted diffusion tensor.
Notes
In the original article of DKI microstructural model
[31], the hindered and restricted tensors were
defined as the intra-cellular and extra-cellular diffusion compartments
respectively.
where \(ADC_{r}\) and \(ADC_{h}\) are the apparent diffusion coefficients of
the diffusion hindered and restricted compartment for a given direction
\(\theta\), \(b\) is the b value provided in the GradientTable input for that
direction, \(f\) is the volume fraction of the restricted diffusion
compartment (also known as the axonal water fraction).
In the original article of DKI microstructural model
[31], the hindered and restricted tensors were
defined as the intra-cellular and extra-cellular diffusion compartments
respectively.
Applies the 3D FFT in the q-space grid to generate the DSI diffusion propagator, remove the background noise with a hard threshold and then deconvolve the propagator with the Lucy-Richardson deconvolution algorithm
rtop_pdf(*[, normalized])
Calculates the return to origin probability from the propagator, which is the propagator evaluated at zero.
rtop_signal(*[, filtering])
Calculates the return to origin probability (rtop) from the signal rtop equals to the sum of all signal values
Applies the 3D FFT in the q-space grid to generate
the DSI diffusion propagator, remove the background noise with a
hard threshold and then deconvolve the propagator with the
Lucy-Richardson deconvolution algorithm
Useful when dMRI data are provided in one qspace hemisphere as
DiffusionSpectrum expects data to be in full qspace.
Parameters:
dataarray, shape (X, Y, Z, W)
where (X, Y, Z) volume size and W number of gradient directions
gtabGradientTable
container for b-values and b-vectors (gradient directions)
Returns:
new_dataarray, shape (X, Y, Z, 2 * W -1)
DWI data across the full Cartesian space.
new_gtabGradientTable
Gradient table.
Notes
We assume here that only on b0 is provided with the initial data. If that
is not the case then you will need to write your own preparation function
before providing the gradients and the data to the DiffusionSpectrumModel
class.
The diffusion orientation distribution function (dODF). This is an
estimate of the diffusion distance in each direction
Parameters:
sphereSphere class instance.
The dODF is calculated in the vertices of this input.
Returns:
odfndarray
The diffusion distance in every direction of the sphere in every
voxel in the input data.
Notes
This is based on equation 3 in [38]. To re-derive it
from scratch, follow steps in [12], Section 7.9
Equation 7.24 but with an \(r^2\) term in the integral.
Given a model fit, predict the signal on the vertices of a sphere
Parameters:
gtaba GradientTable class instance
This encodes the directions for which a prediction is made
S0float array, optional
The mean non-diffusion weighted signal in each voxel. Default:
The fitted S0 value in all voxels if it was fitted. Otherwise 1 in
all voxels.
stepint, optional
The chunk size as a number of voxels. Optional parameter with
default value 10,000.
In order to increase speed of processing, tensor fitting is done
simultaneously over many voxels. This parameter sets the number of
voxels that will be fit at once in each iteration. A larger step
value should speed things up, but it will also take up more memory.
It is advisable to keep an eye on memory consumption as this value
is increased.
Notes
The predicted signal is given by:
\[S(\theta, b) = S_0 * e^{-b ADC}\]
Where:
.. math:
ADC= \thetaQ \theta^T
\(\theta\) is a unit vector pointing at any direction on the sphere for
which a signal is to be predicted and \(b\) is the b value provided in
the GradientTable input for that direction
Note that the notation, \(<D>\), is often used as the mean diffusivity (MD)
of the diffusion tensor and can lead to confusions in the literature
(see [39] versus [40] versus
[41] for example). [40] defines
geodesic anisotropy (GA) with \(<D>\) as the MD in the denominator of the
sum. This is wrong. The original paper [39] defines
GA with \(<D> = det(D)^{1/3}\), as the isotropic part of the distance. This
might be an explanation for the confusion. The isotropic part of the
diffusion tensor in Euclidean space is the MD whereas the isotropic part of
the tensor in log-Euclidean space is \(det(D)^{1/3}\). The Appendix of
[39] and log-Euclidean derivations from
[41] are clear on this. Hence, all that to say that
\(<D> = det(D)^{1/3}\) here for the GA definition and not MD.
The quadratic form of a tensor, or an array with quadratic forms of
tensors. Should be of shape (x, y, z, 3, 3) or (n, 3, 3) or (3, 3).
Returns:
modearray
Calculated tensor mode in each spatial coordinate.
Notes
Mode ranges between -1 (planar anisotropy) and +1 (linear anisotropy)
with 0 representing isotropy. Mode is calculated with the following
equation (equation 9 in [43]):
Tensor parameters. The last dimension should have 12 tensor
parameters: 3 eigenvalues, followed by the 3 corresponding
eigenvectors.
gtaba GradientTable class instance
The gradient table for this prediction
S0float or ndarray
The non diffusion-weighted signal in every voxel, or across all
voxels. Default: 1
Notes
The predicted signal is given by:
\[S(\theta, b) = S_0 * e^{-b ADC}\]
where \(ADC = \theta Q \theta^T\), \(\theta\) is a unit vector pointing at any
direction on the sphere for which a signal is to be predicted, \(b\) is the b
value provided in the GradientTable input for that direction, \(Q\) is the
quadratic form of the tensor determined by the input parameters.
Wrap a fit_tensor func and iterate over chunks of data with given length
Splits data into a number of chunks of specified size and iterates the
decorated fit_tensor function over them. This is useful to counteract the
temporary but significant memory usage increase in fit_tensor functions
that use vectorized operations and need to store large temporary arrays for
their vectorized operations.
Parameters:
stepint, optional
The chunk size as a number of voxels. Optional parameter with default
value 10,000.
In order to increase speed of processing, tensor fitting is done
simultaneously over many voxels. This parameter sets the number of
voxels that will be fit at once in each iteration. A larger step value
should speed things up, but it will also take up more memory. It is
advisable to keep an eye on memory consumption as this value is
increased.
Design matrix holding the covariants used to solve for the regression
coefficients.
dataarray ([X, Y, Z, …], g)
Data or response variables holding the data. Note that the last
dimension should contain the data. It makes no copies of data.
weightsarray ([X, Y, Z, …], g), optional
Weights to apply for fitting. These weights must correspond to the
squared residuals such that \(S = \sum_i w_i r_i^2\).
If not provided, weights are estimated as the squared predicted signal
from an initial OLS fit [44].
return_S0_hatbool, optional
Boolean to return (True) or not (False) the S0 values for the fit.
return_lower_triangularbool, optional
Boolean to return (True) or not (False) the coefficients of the fit.
return_leveragesbool, optional
Boolean to return (True) or not (False) the fitting leverages.
Returns:
eigvalsarray (…, 3)
Eigenvalues from eigen decomposition of the tensor.
eigvecsarray (…, 3, 3)
Associated eigenvectors from eigen decomposition of the tensor.
Eigenvectors are columnar (e.g. eigvecs[:,j] is associated with
eigvals[j])
leveragesarray (g)
Leverages of the fitting problem (if return_leverages is True)
In Chung, et al. 2006, the regression of the WLS fit needed an unbiased
preliminary estimate of the weights and therefore the ordinary least
squares (OLS) estimates were used. A “two pass” method was implemented:
calculate OLS estimates of the data
apply the OLS estimates as weights to the WLS fit of the data
This ensured heteroscedasticity could be properly modeled for various
types of bootstrap resampling (namely residual bootstrap).
\[\begin{split}y = \mathrm{data} \\
X = \mathrm{design matrix} \\
\hat{\beta}_\mathrm{WLS} =
\mathrm{desired regression coefficients (e.g. tensor)}\\
\\
\hat{\beta}_\mathrm{WLS} = (X^T W X)^{-1} X^T W y \\
\\
W = \mathrm{diag}((X \hat{\beta}_\mathrm{OLS})^2),
\mathrm{where} \hat{\beta}_\mathrm{OLS} = (X^T X)^{-1} X^T y\end{split}\]
Fit the cumulant expansion params (e.g. DTI, DKI) using non-linear
least-squares.
Parameters:
design_matrixarray (g, Npar)
Design matrix holding the covariants used to solve for the regression
coefficients. First six parameters of design matrix should correspond
to the six unique diffusion tensor elements in the lower triangular
order (Dxx, Dxy, Dyy, Dxz, Dyz, Dzz), while last parameter to -log(S0)
dataarray ([X, Y, Z, …], g)
Data or response variables holding the data. Note that the last
dimension should contain the data. It makes no copies of data.
weightsarray ([X, Y, Z, …], g), optional
Weights to apply for fitting. These weights must correspond to the
squared residuals such that \(S = \sum_i w_i r_i^2\).
jacbool, optional
Use the Jacobian?
return_S0_hatbool, optional
Boolean to return (True) or not (False) the S0 values for the fit.
fail_is_nanbool, optional
Boolean to set failed NL fitting to NaN (True) or LS (False, default).
return_lower_triangularbool, optional
Boolean to return (True) or not (False) the coefficients of the fit.
return_leveragesbool, optional
Boolean to return (True) or not (False) the fitting leverages.
init_paramsarray ([X, Y, Z, …], Npar), optional
Parameters in lower triangular form as initial optimization guess.
Returns:
nlls_params: the eigen-values and eigen-vectors of the tensor in each
Compute a robust tensor fit using the RESTORE algorithm.
Note that the RESTORE algorithm defined in [45] does not define
Geman–McClure M-estimator weights as claimed (instead, Cauchy M-estimator
weights are defined), but this function does define correct Geman–McClure
M-estimator weights.
Parameters:
design_matrixarray of shape (g, 7)
Design matrix holding the covariants used to solve for the regression
coefficients.
dataarray of shape ([X, Y, Z, n_directions], g)
Data or response variables holding the data. Note that the last
dimension should contain the data. It makes no copies of data.
sigmafloat, optional
An estimate of the variance. [45] recommend to use
1.5267 * std(background_noise), where background_noise is estimated
from some part of the image known to contain no signal (only noise).
If not provided, will be estimated per voxel as:
sigma = 1.4826 * sqrt(N / (N - p)) * MAD(residuals)
as in [46] but with the additional correction factor
1.4826 required to link standard deviation to MAD.
jacbool, optional
Whether to use the Jacobian of the tensor to speed the non-linear
optimization procedure used to fit the tensor parameters (see also
nlls_fit_tensor()).
return_S0_hatbool, optional
Boolean to return (True) or not (False) the S0 values for the fit.
fail_is_nanbool, optional
Boolean to set failed NL fitting to NaN (True) or LS (False).
Returns:
restore_paramsan estimate of the tensor parameters in each voxel.
Take care to supply an appropriate weights_method for the fit_type.
It is possible to use NLLS fitting with weights designed for WLS fitting,
but this is a user error.
Returns eigenvalues and eigenvectors given a diffusion tensor
Computes tensor eigen decomposition to calculate eigenvalues and
eigenvectors (Basser et al., 1994a).
Parameters:
tensorarray (…, 3, 3)
Hermitian matrix representing a diffusion tensor.
min_diffusivityfloat, optional
Because negative eigenvalues are not physical and small eigenvalues,
much smaller than the diffusion weighting, cause quite a lot of noise
in metrics such as fa, diffusivity values smaller than
min_diffusivity are replaced with min_diffusivity.
Returns:
eigvalsarray (…, 3)
Eigenvalues from eigen decomposition of the tensor. Negative
eigenvalues are replaced by zero. Sorted from largest to smallest.
eigvecsarray (…, 3, 3)
Associated eigenvectors from eigen decomposition of the tensor.
Eigenvectors are columnar (e.g. eigvecs[…, :, j] is associated with
eigvals[…, j])
Fiber ORientation Estimated using Continuous Axially Symmetric Tensors
(FORECAST).
FORECAST [47], [48],
[49] is a Spherical Deconvolution reconstruction
model for multi-shell diffusion data which enables the calculation of a
voxel adaptive response function using the Spherical Mean Technique (SMT)
[48], [49].
With FORECAST it is possible to calculate crossing invariant parallel
diffusivity, perpendicular diffusivity, mean diffusivity, and fractional
anisotropy [48].
Predict a signal for this TensorModel class instance given
parameters.
Parameters:
fwdti_params(…, 13) ndarray
The last dimension should have 13 parameters: the 12 tensor
parameters (3 eigenvalues, followed by the 3 corresponding
eigenvectors) and the free water volume fraction.
S0float or ndarray
The non diffusion-weighted signal in every voxel, or across all
voxels. Default: 1
Returns:
S(…, N) ndarray
Simulated signal based on the free water DTI model
Signal prediction given the free water DTI model parameters.
Parameters:
params(…, 13) ndarray
Model parameters. The last dimension should have the 12 tensor
parameters (3 eigenvalues, followed by the 3 corresponding
eigenvectors) and the volume fraction of the free water compartment.
gtaba GradientTable class instance
The gradient table for this prediction
S0float or ndarray
The non diffusion-weighted signal in every voxel, or across all
voxels. Default: 1
Disofloat, optional
Value of the free water isotropic diffusion. Default is set to 3e-3
\(mm^{2}.s^{-1}\). Please adjust this value if you are assuming different
units of diffusion.
Returns:
S(…, N) ndarray
Simulated signal based on the free water DTI model
Notes
The predicted signal is given by:
\(S(\theta, b) = S_0 * [(1-f) * e^{-b ADC} + f * e^{-b D_{iso}]\), where
\(ADC = \theta Q \theta^T\), \(\theta\) is a unit vector pointing at any
direction on the sphere for which a signal is to be predicted, \(b\) is the b
value provided in the GradientTable input for that direction, \(Q\) is the
quadratic form of the tensor determined by the input parameters, \(f\) is the
free water diffusion compartment, \(D_{iso}\) is the free water diffusivity
which is equal to $3 * 10^{-3} mm^{2}s^{-1} [50].
Applies weighted linear least squares fit of the water free elimination
model to single voxel signals.
Parameters:
design_matrixarray (g, 7)
Design matrix holding the covariants used to solve for the regression
coefficients.
sigarray (g, )
Diffusion-weighted signal for a single voxel data.
S0float
Non diffusion weighted signal (i.e. signal for b-value=0).
Disofloat, optional
Value of the free water isotropic diffusion. Default is set to 3e-3
\(mm^{2}.s^{-1}\). Please adjust this value if you are assuming different
units of diffusion.
mdregfloat, optimal
DTI’s mean diffusivity regularization threshold. If standard DTI
diffusion tensor’s mean diffusivity is almost near the free water
diffusion value, the diffusion signal is assumed to be only free water
diffusion (i.e. volume fraction will be set to 1 and tissue’s diffusion
parameters are set to zero). Default md_reg is 2.7e-3 \(mm^{2}.s^{-1}\)
(corresponding to 90% of the free water diffusion value).
min_signalfloat
The minimum signal value. Needs to be a strictly positive
number. Default: minimal signal in the data provided to fit.
piterationsinter, optional
Number of iterations used to refine the precision of f. Default is set
to 3 corresponding to a precision of 0.01.
Returns:
fw_paramsndarray
All parameters estimated from the free water tensor model. Parameters
are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the
first, second and third coordinates of the eigenvector
The gradient table containing diffusion acquisition parameters.
datandarray ([X, Y, Z, …], g)
Data or response variables holding the data. Note that the last
dimension should contain the data. It makes no copies of data.
Disofloat, optional
Value of the free water isotropic diffusion. Default is set to 3e-3
\(mm^{2}.s^{-1}\). Please adjust this value if you are assuming different
units of diffusion.
maskarray, optional
A boolean array used to mark the coordinates in the data that should
be analyzed that has the shape data.shape[:-1]
min_signalfloat
The minimum signal value. Needs to be a strictly positive
number.
piterationsinter, optional
Number of iterations used to refine the precision of f. Default is set
to 3 corresponding to a precision of 0.01.
mdregfloat, optimal
DTI’s mean diffusivity regularization threshold. If standard DTI
diffusion tensor’s mean diffusivity is almost near the free water
diffusion value, the diffusion signal is assumed to be only free water
diffusion (i.e. volume fraction will be set to 1 and tissue’s diffusion
parameters are set to zero). Default md_reg is 2.7e-3 \(mm^{2}.s^{-1}\)
(corresponding to 90% of the free water diffusion value).
Returns:
fw_paramsndarray (x, y, z, 13)
Matrix containing in the last dimension the free water model parameters
in the following order:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the
first, second and third coordinates of the eigenvector
The volume fraction of the free water compartment.
Applies non linear least squares fit of the water free elimination
model to single voxel signals.
Parameters:
design_matrixarray (g, 7)
Design matrix holding the covariants used to solve for the regression
coefficients.
sigarray (g, )
Diffusion-weighted signal for a single voxel data.
S0float
Non diffusion weighted signal (i.e. signal for b-value=0).
Disofloat, optional
Value of the free water isotropic diffusion. Default is set to 3e-3
\(mm^{2}.s^{-1}\). Please adjust this value if you are assuming different
units of diffusion.
mdregfloat, optimal
DTI’s mean diffusivity regularization threshold. If standard DTI
diffusion tensor’s mean diffusivity is almost near the free water
diffusion value, the diffusion signal is assumed to be only free water
diffusion (i.e. volume fraction will be set to 1 and tissue’s diffusion
parameters are set to zero). Default md_reg is 2.7e-3 \(mm^{2}.s^{-1}\)
(corresponding to 90% of the free water diffusion value).
min_signalfloat, optional
The minimum signal value. Needs to be a strictly positive
number.
choleskybool, optional
If true it uses Cholesky decomposition to ensure that diffusion tensor
is positive define.
f_transformbool, optional
If true, the water volume fractions is converted during the convergence
procedure to ft = arcsin(2*f - 1) + pi/2, insuring f estimates between
0 and 1.
jacbool, optional
True to use the Jacobian.
weighting: str, optional
the weighting scheme to use in considering the
squared-error. Default behavior is to use uniform weighting. Other
options: ‘sigma’ ‘gmm’
sigma: float, optional
If the ‘sigma’ weighting scheme is used, a value of sigma needs to be
provided here. According to Chang et al.[45], a good value to use
is 1.5267 * std(background_noise), where background_noise is estimated
from some part of the image known to contain no signal (only noise).
Returns:
All parameters estimated from the free water tensor model.
Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the
first, second and third coordinates of the eigenvector
The volume fraction of the free water compartment.
Fit the water elimination tensor model using the non-linear least-squares.
Parameters:
gtaba GradientTable class instance
The gradient table containing diffusion acquisition parameters.
datandarray ([X, Y, Z, …], g)
Data or response variables holding the data. Note that the last
dimension should contain the data. It makes no copies of data.
maskarray, optional
A boolean array used to mark the coordinates in the data that should
be analyzed that has the shape data.shape[:-1]
Disofloat, optional
Value of the free water isotropic diffusion. Default is set to 3e-3
\(mm^{2}.s^{-1}\). Please adjust this value if you are assuming different
units of diffusion.
mdregfloat, optimal
DTI’s mean diffusivity regularization threshold. If standard DTI
diffusion tensor’s mean diffusivity is almost near the free water
diffusion value, the diffusion signal is assumed to be only free water
diffusion (i.e. volume fraction will be set to 1 and tissue’s diffusion
parameters are set to zero). Default md_reg is 2.7e-3 \(mm^{2}.s^{-1}\)
(corresponding to 90% of the free water diffusion value).
min_signalfloat, optional
The minimum signal value. Needs to be a strictly positive
number.
f_transformbool, optional
If true, the water volume fractions is converted during the convergence
procedure to ft = arcsin(2*f - 1) + pi/2, insuring f estimates between
0 and 1.
choleskybool, optional
If true it uses Cholesky decomposition to ensure that diffusion tensor
is positive define.
jacbool, optional
True to use the Jacobian.
weighting: str, optional
the weighting scheme to use in considering the
squared-error. Default behavior is to use uniform weighting. Other
options: ‘sigma’ ‘gmm’
sigma: float, optional
If the ‘sigma’ weighting scheme is used, a value of sigma needs to be
provided here. According to Chang et al.[45], a good value to use
is 1.5267 * std(background_noise), where background_noise is estimated
from some part of the image known to contain no signal (only noise).
Returns:
fw_paramsndarray (x, y, z, 13)
Matrix containing in the dimension the free water model parameters in
the following order:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the
first, second and third coordinates of the eigenvector
maximum qa value. Usually found in the CSF (corticospinal fluid).
Returns:
nqaarray, shape (x, Y, Z, N)
normalized quantitative anisotropy
Notes
Normalized quantitative anisotropy has the very useful property
to be very small near gray matter and background areas. Therefore,
it can be used to mask out white matter areas.
Predict a signal for this IvimModel class instance given parameters.
Parameters:
ivim_paramsarray
The ivim parameters as an array [S0, f, D_star and D]
gtabGradientTable class instance
Gradient directions and bvalues.
S0float, optional
This has been added just for consistency with the existing
API. Unlike other models, IVIM predicts S0 and this is over written
by the S0 value in params.
Creates a structure for the combining the diffusion and pseudo- diffusion by multiplying with the bvals and then exponentiating each of the two components for fitting as per the IVIM- two compartment model.
Performs the constrained search for the linear parameters f after
the estimation of x is done. Estimation of the linear parameters f
is a constrained linear least-squares optimization problem solved by
using a convex optimizer from cvxpy. The IVIM equation contains two
parameters that depend on the same volume fraction. Both are estimated
as separately in the convex optimizer.
Parameters:
phiarray
Returns an array calculated from :func: phi.
signalarray
The signal values measured for this model.
Returns:
f1, f2 (volume fractions)
Notes
cost function for differential evolution algorithm:
Constructs the objective for the :func: stoc_search_cost.
First calculates the Moore-Penrose inverse of the input phi and takes
a dot product with the measured signal. The result obtained is again
multiplied with phi to complete the projection of the variable into
a transformed space. (see [53] and
[54] for thorough discussion on Variable Projections
and relevant cost functions).
Parameters:
phiarray
Returns an array calculated from :func: Phi.
signalarray
The signal values measured for this model.
Returns:
(signal - S)^T(signal - S)
Notes
to make cost function for Differential Evolution algorithm:
.. math:
Cost function for the least square problem. The cost function is used
in the Least Squares function of SciPy in :func: fit. It guarantees
that stopping point of the algorithm is at least a stationary point
with reduction in the number of iterations required by the
differential evolution optimizer.
Parameters:
x_farray
Contains the parameters ‘x’ and ‘f’ combines in the same array.
Creates a structure for the combining the diffusion and pseudo-
diffusion by multiplying with the bvals and then exponentiating each of
the two components for fitting as per the IVIM- two compartment model.
Parameters:
xarray
input from the Differential Evolution optimizer.
Returns:
exp_phi1array
Combined array of parameters perfusion/pseudo-diffusion
and diffusion parameters.
Cost function for differential evolution algorithm. Performs a
stochastic search for the non-linear parameters ‘x’. The objective
function is calculated in the :func: ivim_mix_cost_one. The function
constructs the parameters using :func: phi.
S0 value here is not necessary and will not be used to predict the
signal. It has been added to conform to the structure of the
predict method in multi_voxel which requires a keyword argument S0.
Returns:
signalarray
The signal values predicted for this model using its parameters.
The Intravoxel incoherent motion (IVIM) model function.
Parameters:
paramsarray
An array of IVIM parameters - [S0, f, D_star, D].
gtabGradientTable class instance
Gradient directions and bvalues.
S0float, optional
This has been added just for consistency with the existing
API. Unlike other models, IVIM predicts S0 and this is over written
by the S0 value in params.
Returns:
Sarray
An array containing the IVIM signal estimated using given parameters.
Selector function to switch between the 2-stage Trust-Region Reflective
based NLLS fitting method (also containing the linear fit): trr and the
Variable Projections based fitting method: varpro.
Parameters:
fit_methodstring, optional
The value fit_method can either be ‘trr’ or ‘varpro’.
default : trr
Mean Apparent Propagator MRI (MAPMRI) of the diffusion signal.
The main idea in MAPMRI footcite:p:Ozarslan2013 is to model the diffusion
signal as a linear combination of the continuous functions presented in
footcite:p:Ozarslan2008 but extending it in three dimensions.
The main difference with the SHORE proposed in footcite:p:Merlet2013 is
that MAPMRI 3D extension is provided using a set of three basis functions
for the radial part, one for the signal along x, one for y and one for z,
while footcite:p:Merlet2013 uses one basis function to model the radial
part and real Spherical Harmonics to model the angular part.
From the MAPMRI coefficients is possible to use the analytical formulae
to estimate the ODF.
See [55] for additional tissue microstructure insights
provided by MAPMRI.
See also footcite:p:Fick2016b, footcite:p:Cheng2012,
footcite:p:Hosseinbor2013, footcite:p:Craven1979, and
footcite:p:DelaHaije2020 for additional insight into to the model.
Calculates the analytical Mean Squared Displacement (MSD).
It is defined as the Laplacian of the origin of the estimated signal
[56]. The analytical formula for the MAP-MRI basis
was derived in [57] eq. (C13, D1).
For the NG to be meaningful the mapmri scale factors must be estimated
only on data representing Gaussian diffusion of spins, i.e., bvals
smaller than about 2000 s/mm^2 [55].
Calculates the analytical parallel non-Gaussiannity (NG).
For the NG to be meaningful the mapmri scale factors must be estimated
only on data representing Gaussian diffusion of spins, i.e., bvals
smaller than about 2000 s/mm^2 [55].
Calculates the analytical perpendicular non-Gaussiannity (NG)
For the NG to be meaningful the mapmri scale factors must be estimated
only on data representing Gaussian diffusion of spins, i.e., bvals
smaller than about 2000 s/mm^2 [55].
Calculates the norm of the laplacian of the fitted signal.
This information could be useful to assess if the extrapolation of the
fitted signal contains spurious oscillations. A high laplacian may
indicate that these are present, and any q-space indices that use
integrals of the signal may be corrupted (e.g. RTOP, RTAP, RTPP, QIV).
Calculates the real analytical odf for a given discrete sphere.
Computes the design matrix of the ODF for the given sphere vertices
and radial moment [58] eq. (32). The radial moment
s acts as a sharpening method. The analytical equation for the spherical
ODF basis is given in [57] eq. (C8).
Diffusion propagator on a given set of real points.
if the array r_points is non writeable, then intermediate
results are cached for faster recalculation
Calculates the analytical Q-space Inverse Variance (QIV).
It is defined as the inverse of the Laplacian of the origin of the
estimated propagator [59] eq. (22). The
analytical formula for the MAP-MRI basis was derived in
[57] eq. (C14, D2).
The computation follows [58] Eq. 32, but it is done
for the isotropic propagator in footcite:p:Ozarslan2013 eq. (60).
Analytical derivation in [57] eq. (C8).
Parameters:
radial_orderunsigned int,
an even integer that represent the order of the basis
mufloat,
isotropic scale factor of the isotropic MAP-MRI basis
sunsigned int
radial moment of the ODF
verticesarray, shape (N,3)
points of the sphere shell in the r-space in which evaluate the ODF
The computation follows [58] Eq. 32, but it is done
for the isotropic propagator in [58] eq. (60). Here
we do not compute the sphere function but the spherical harmonics by only
integrating the radial part of the propagator. We use the same derivation of
the ODF in the isotropic implementation as in [57] eq.
(C8).
Parameters:
radial_orderunsigned int,
an even integer that represent the order of the basis
mufloat,
isotropic scale factor of the isotropic MAP-MRI basis
Here weights_array is a numpy array with all values that should be
considered in the GCV. It will run through the weights until the cost
function starts to increase, then stop and take the last value as the
optimum weight.
Helper function to set up and solve the Quadratic Program (QP) in CVXPY.
A QP problem has the following form:
minimize 1/2 x’ P x + Q’ x
subject to G x <= H
Here the QP solver is based on CVXPY and uses OSQP.
Parameters:
Pndarray
n x n matrix for the primal QP objective function.
Qndarray
n x 1 matrix for the primal QP objective function.
Fiber response function estimation for multi-shell data.
Parameters:
sh_order_maxint
Maximum spherical harmonics order (l).
bvalsndarray
Array containing the b-values. Must be unique b-values, like outputted
by dipy.core.gradients.unique_bvals_tolerance.
wm_rf(N-1, 4) ndarray
Response function of the WM tissue, for each bvals,
where N is the number of unique b-values including the b0.
gm_rf(N-1, 4) ndarray
Response function of the GM tissue, for each bvals.
csf_rf(N-1, 4) ndarray
Response function of the CSF tissue, for each bvals.
spheredipy.core.Sphere instance, optional
Sphere where the signal will be evaluated.
tolint, optional
Tolerance gap for b-values clustering.
btenscan be any of two options, optional
an array of strings of shape (N,) specifying
encoding tensor shape associated with all unique b-values
separately. N corresponds to the number of unique b-values,
including the b0. Options for elements in array: ‘LTE’,
‘PTE’, ‘STE’, ‘CTE’ corresponding to linear, planar, spherical, and
“cigar-shaped” tensor encoding.
an array of shape (N,3,3) specifying the b-tensor of each unique
b-values exactly. N corresponds to the number of unique b-values,
including the b0.
Computation of masks for multi-shell multi-tissue (msmt) response
function using FA and MD.
Parameters:
gtabGradientTable
Gradient table.
datandarray
diffusion data (4D)
roi_centerarray-like, (3,)
Center of ROI in data. If center is None, it is assumed that it is
the center of the volume with shape data.shape[:3].
roi_radiiint or array-like, (3,)
radii of cuboid ROI
wm_fa_thrfloat
FA threshold for WM.
gm_fa_thrfloat
FA threshold for GM.
csf_fa_thrfloat
FA threshold for CSF.
gm_md_thrfloat
MD threshold for GM.
csf_md_thrfloat
MD threshold for CSF.
Returns:
mask_wmndarray
Mask of voxels within the ROI and with FA above the FA threshold
for WM.
mask_gmndarray
Mask of voxels within the ROI and with FA below the FA threshold
for GM and with MD below the MD threshold for GM.
mask_csfndarray
Mask of voxels within the ROI and with FA below the FA threshold
for CSF and with MD below the MD threshold for CSF.
Notes
In msmt-CSD there is an important pre-processing step: the estimation of
every tissue’s response function. In order to do this, we look for voxels
corresponding to WM, GM and CSF. This function aims to accomplish that by
returning a mask of voxels within a ROI and who respect some threshold
constraints, for each tissue. More precisely, the WM mask must have a FA
value above a given threshold. The GM mask and CSF mask must have a FA
below given thresholds and a MD below other thresholds. To get the FA and
MD, we need to fit a Tensor model to the datasets.
(evals, S0) for CSF for each unique bvalues (except b0).
Notes
In msmt-CSD there is an important pre-processing step: the estimation of
every tissue’s response function. In order to do this, we look for voxels
corresponding to WM, GM and CSF. This information can be obtained by using
mcsd.mask_for_response_msmt() through masks of selected voxels. The present
function uses such masks to compute the msmt response functions.
For the responses, we base our approach on the function
csdeconv.response_from_mask_ssst(), with the added layers of multishell and
multi-tissue (see the ssst function for more information about the
computation of the ssst response function). This means that for each tissue
we use the previously found masks and loop on them. For each mask, we loop
on the b-values (clustered using the tolerance gap) to get many responses
and then average them to get one response per tissue.
(evals, S0) for CSF for each unique bvalues (except b0).
Notes
In msmt-CSD there is an important pre-processing step: the estimation of
every tissue’s response function. In order to do this, we look for voxels
corresponding to WM, GM and CSF. We get this information from
mcsd.mask_for_response_msmt(), which returns masks of selected voxels
(more details are available in the description of the function).
With the masks, we compute the response functions by using
mcsd.response_from_mask_msmt(), which returns the response for each
tissue (more details are available in the description of the function).
Computes the microscopic fractional anisotropy from the mean signal diffusional kurtosis parameters assuming the 2-compartmental spherical mean technique model.
This encodes the directions for which a prediction is made
S0float array
The mean non-diffusion weighted signal in each voxel. Default:
The fitted S0 value in all voxels if it was fitted. Otherwise 1 in
all voxels.
Returns:
S(…, N) ndarray
Simulated mean signal based on the mean signal kurtosis model
Notes
The predicted signal is given by:
\(MS(b) = S_0 * exp(-bD + 1/6 b^{2} D^{2} K)\), where \(D\) and \(k\) are the
mean signal diffusivity and mean signal kurtosis.
Computes the microscopic fractional anisotropy from the mean signal
diffusional kurtosis parameters assuming the 2-compartmental spherical
mean technique model.
See [61] and [62] for
further details about the method.
Returns:
smt2uFAndarray
Microscopic fractional anisotropy computed by converting MSDKI to
SMT2.
Notes
Computes the intrinsic diffusivity using equation 10 of
[62].
Computes the average signal across different diffusion directions
for each unique b-value
Parameters:
datandarray ([X, Y, Z, …], g)
ndarray containing the data signals in its last dimension.
gtaba GradientTable class instance
The gradient table containing diffusion acquisition parameters.
bmagThe order of magnitude that the bvalues have to differ to be
considered an unique b-value. Default: derive this value from the
maximal b-value provided: \(bmag=log_{10}(max(bvals)) - 1\).
Returns:
msignalndarray ([X, Y, Z, …, nub])
Mean signal along all gradient directions for each unique b-value
Note that the last dimension contains the signal means and nub is the
number of unique b-values.
ngndarray(nub)
Number of gradient directions used to compute the mean signal for
all unique b-values
Notes
This function assumes that directions are evenly sampled on the sphere or
on the hemisphere
Design matrix holding the covariants used to solve for the regression
coefficients of the mean signal diffusion kurtosis model. Note that
nub is the number of unique b-values
msignalndarray ([X, Y, Z, …, nub])
Mean signal along all gradient directions for each unique b-value
Note that the last dimension should contain the signal means and nub
is the number of unique b-values.
ngndarray(nub)
Number of gradient directions used to compute the mean signal for
all unique b-values
maskarray
A boolean array used to mark the coordinates in the data that
should be analyzed that has the shape data.shape[:-1]
min_signalfloat, optional
Voxel with mean signal intensities lower than the min positive signal
are not processed. Default: 0.0001
return_S0_hatbool
If True, also return S0 values for the fit.
Returns:
paramsarray (…, 2)
Containing the mean signal diffusivity and mean signal kurtosis
Constructs design matrix for the mean signal diffusion kurtosis model
Parameters:
ubvalsarray
Containing the unique b-values of the data.
Returns:
design_matrixarray (nb, 3)
Design matrix or B matrix for the mean signal diffusion kurtosis
model assuming that parameters are in the following order:
design_matrix[j, :] = (msd, msk, S0)
An object to simplify the interaction of the array with the ctypes module.
data
Python buffer object pointing to the start of the array’s data.
device
dtype
Data-type of the array’s elements.
flags
Information about the memory layout of the array.
flat
A 1-D iterator over the array.
imag
The imaginary part of the array.
itemset
itemsize
Length of one array element in bytes.
mT
View of the matrix transposed array.
nbytes
Total bytes consumed by the elements of the array.
ndim
Number of array dimensions.
newbyteorder
ptp
real
The real part of the array.
shape
Tuple of array dimensions.
size
Number of elements in the array.
strides
Tuple of bytes to step in each dimension when traversing an array.
Methods
__call__(*args, **kwargs)
Call self as a function.
all([axis, out, keepdims, where])
Returns True if all elements evaluate to True.
any([axis, out, keepdims, where])
Returns True if any of the elements of a evaluate to True.
argmax([axis, out, keepdims])
Return indices of the maximum values along the given axis.
argmin([axis, out, keepdims])
Return indices of the minimum values along the given axis.
argpartition(kth[, axis, kind, order])
Returns the indices that would partition this array.
argsort([axis, kind, order])
Returns the indices that would sort this array.
astype(dtype[, order, casting, subok, copy])
Copy of the array, cast to a specified type.
byteswap([inplace])
Swap the bytes of the array elements
choose(choices[, out, mode])
Use an index array to construct a new array from a set of choices.
clip([min, max, out])
Return an array whose values are limited to [min,max].
compress(condition[, axis, out])
Return selected slices of this array along given axis.
conj()
Complex-conjugate all elements.
conjugate()
Return the complex conjugate, element-wise.
copy([order])
Return a copy of the array.
cumprod([axis, dtype, out])
Return the cumulative product of the elements along the given axis.
cumsum([axis, dtype, out])
Return the cumulative sum of the elements along the given axis.
diagonal([offset, axis1, axis2])
Return specified diagonals.
dump(file)
Dump a pickle of the array to the specified file.
dumps()
Returns the pickle of the array as a string.
fill(value)
Fill the array with a scalar value.
flatten([order])
Return a copy of the array collapsed into one dimension.
getfield(dtype[, offset])
Returns a field of the given array as a certain type.
item(*args)
Copy an element of an array to a standard Python scalar and return it.
max([axis, out, keepdims, initial, where])
Return the maximum along a given axis.
mean([axis, dtype, out, keepdims, where])
Returns the average of the array elements along given axis.
min([axis, out, keepdims, initial, where])
Return the minimum along a given axis.
nonzero()
Return the indices of the elements that are non-zero.
partition(kth[, axis, kind, order])
Partially sorts the elements in the array in such a way that the value of the element in k-th position is in the position it would be in a sorted array.
prod([axis, dtype, out, keepdims, initial, ...])
Return the product of the array elements over the given axis
put(indices, values[, mode])
Set a.flat[n]=values[n] for all n in indices.
ravel([order])
Return a flattened array.
repeat(repeats[, axis])
Repeat elements of an array.
reshape(shape, /, *[, order])
Returns an array containing the same data with a new shape.
resize(new_shape[, refcheck])
Change shape and size of array in-place.
round([decimals, out])
Return a with each element rounded to the given number of decimals.
searchsorted(v[, side, sorter])
Find indices where elements of v should be inserted in a to maintain order.
setfield(val, dtype[, offset])
Put a value into a specified place in a field defined by a data-type.
setflags([write, align, uic])
Set array flags WRITEABLE, ALIGNED, WRITEBACKIFCOPY, respectively.
sort([axis, kind, order])
Sort an array in-place.
squeeze([axis])
Remove axes of length one from a.
std([axis, dtype, out, ddof, keepdims, where])
Returns the standard deviation of the array elements along given axis.
sum([axis, dtype, out, keepdims, initial, where])
Return the sum of the array elements over the given axis.
swapaxes(axis1, axis2)
Return a view of the array with axis1 and axis2 interchanged.
take(indices[, axis, out, mode])
Return an array formed from the elements of a at the given indices.
tobytes([order])
Construct Python bytes containing the raw data bytes in the array.
tofile(fid[, sep, format])
Write array to a file as text or binary (default).
tolist()
Return the array as an a.ndim-levels deep nested list of Python scalars.
tostring([order])
A compatibility alias for ~ndarray.tobytes, with exactly the same behavior.
trace([offset, axis1, axis2, dtype, out])
Return the sum along diagonals of the array.
transpose(*axes)
Returns a view of the array with axes transposed.
var([axis, dtype, out, ddof, keepdims, where])
Returns the variance of the array elements, along given axis.
Where \(\Psi\) is an orientation distribution function sampled discretely on
the unit sphere and angle brackets denote average over the samples on the
sphere.
The q:math:tau-dMRI model to analytically and continuously represent the
q:math:tau diffusion signal attenuation over diffusion sensitization q and
diffusion time \(\tau\).
The model [64] can be seen as an extension of the MAP-MRI
basis [58] towards different diffusion times.
The main idea is to model the diffusion signal over time and space as
a linear combination of continuous functions,
where \(\Phi\) and \(T\) are the spatial and temporal basis functions,
\(N_{\textbf{q}}\) and \(N_\tau\) are the maximum spatial and temporal
order, and \(i,k\) are basis order iterators.
The estimation of the coefficients \(c_i\) can be regularized using
either analytic Laplacian regularization, sparsity regularization using
the l1-norm, or both to do a type of elastic net regularization.
From the coefficients, there exists an analytical formula to estimate
the ODF, RTOP, RTAP, RTPP, QIV and MSD, for any diffusion time.
Parameters:
gtabGradientTable,
gradient directions and bvalues container class. The bvalues
should be in the normal s/mm^2. big_delta and small_delta need to be
given in seconds.
radial_orderunsigned int,
an even integer representing the spatial/radial order of the basis.
time_orderunsigned int,
an integer larger or equal than zero representing the time order
of the basis.
laplacian_regularizationbool,
Regularize using the Laplacian of the qt-dMRI basis.
laplacian_weighting: string or scalar,
The string ‘GCV’ makes it use generalized cross-validation to find
the regularization weight [60]. A scalar sets the
regularization weight to that value.
l1_regularizationbool,
Regularize by imposing sparsity in the coefficients using the
l1-norm.
l1_weighting‘CV’ or scalar,
The string ‘CV’ makes it use five-fold cross-validation to find
the regularization weight. A scalar sets the regularization weight
to that value.
cartesianbool
Whether to use the Cartesian or spherical implementation of the
qt-dMRI basis, which we first explored in [65].
anisotropic_scalingbool
Whether to use anisotropic scaling or isotropic scaling. This
option can be used to test if the Cartesian implementation is
equivalent with the spherical one when using the same scaling.
normalizationbool
Whether to normalize the basis functions such that their inner
product is equal to one. Normalization is only necessary when
imposing sparsity in the spherical basis if cartesian=False.
constrain_q0bool
whether to constrain the q0 point to unity along the tau-space.
This is necessary to ensure that \(E(0,\tau)=1\).
bval_thresholdfloat
the threshold b-value to be used, such that only data points below
that threshold are used when estimating the scale factors.
eigenvalue_thresholdfloat,
Sets the minimum of the tensor eigenvalues in order to avoid
stability problem.
cvxpy_solverstr, optional
cvxpy solver name. Optionally optimize the positivity constraint
with a particular cvxpy solver. See https://www.cvxpy.org/ for
details.
Default: None (cvxpy chooses its own solver)
Calculates the analytical Mean Squared Displacement (MSD) for a
given diffusion time tau.
It is defined as the Laplacian of the origin of the estimated signal
Cheng[56]. The analytical formula for the MAP-MRI basis
was derived in [57] eq. (C13, D1). The qtdmri
coefficients are first converted to mapmri coefficients following
[64].
Parameters:
taufloat
diffusion time (big_delta - small_delta / 3.) in seconds
Calculates the norm of the laplacian of the fitted signal.
This information could be useful to assess if the extrapolation of the
fitted signal contains spurious oscillations. A high laplacian norm may
indicate that these are present, and any q-space indices that
use integrals of the signal may be corrupted (e.g. RTOP, RTAP, RTPP,
QIV). In contrast to Fick et al.[57], the Laplacian now
describes oscillations in the 4-dimensional qt-signal
[64].
Calculates the real analytical odf for a given discrete sphere.
Computes the design matrix of the ODF for the given sphere vertices
and radial moment [58] eq. (32). The radial moment
s acts as a sharpening method. The analytical equation for the spherical
ODF basis is given in [57] eq. (C8). The qtdmri
coefficients are first converted to mapmri coefficients following
[64].
Parameters:
taufloat
diffusion time (big_delta - small_delta / 3.) in seconds
Diffusion propagator on a given set of real points.
if the array r_points is non writeable, then intermediate
results are cached for faster recalculation
Calculates the analytical Q-space Inverse Variance (QIV) for given
diffusion time tau.
It is defined as the inverse of the Laplacian of the origin of the
estimated propagator [59] eq. (22). The
analytical formula for the MAP-MRI basis was derived in
[57] eq. (C14, D2). The qtdmri coefficients are first
converted to mapmri coefficients following Fick et al.[64].
Parameters:
taufloat
diffusion time (big_delta - small_delta / 3.) in seconds
This function converts the qtdmri coefficients to mapmri
coefficients for a given tau.
Defined in [64], the conversion is performed by a
matrix multiplication that evaluates the time-depenent part of the basis
and multiplies it with the coefficients, after which coefficients with
the same spatial orders are summed up, resulting in mapmri coefficients.
Parameters:
taufloat
diffusion time (big_delta - small_delta / 3.) in seconds
Calculates the analytical return to the axis probability (RTAP)
for a given diffusion time tau.
See [58] eq. (40, 44a). The analytical formula for
the isotropic MAP-MRI basis was derived in [57] eq.
(C11). The qtdmri coefficients are first converted to mapmri
coefficients following [64].
Parameters:
taufloat
diffusion time (big_delta - small_delta / 3.) in seconds
Calculates the analytical return to the origin probability (RTOP)
for a given diffusion time tau.
See [58] eq. (36, 43). The analytical formula for
the isotropic MAP-MRI basis was derived in [57] eq.
(C11). The qtdmri coefficients are first converted to mapmri
coefficients following [64].
Parameters:
taufloat
diffusion time (big_delta - small_delta / 3.) in seconds
Calculates the analytical return to the plane probability (RTPP)
for a given diffusion time tau.
See [58] eq. (42). The analytical formula for the
isotropic MAP-MRI basis was derived in [57] eq. (C11).
The qtdmri coefficients are first converted to mapmri coefficients
following [64].
Parameters:
taufloat
diffusion time (big_delta - small_delta / 3.) in seconds
Normalization factor for Spherical MAP-MRI basis. The normalization
for a basis function with orders [j,l,m] depends only on orders j,l and
the isotropic scale factor.
dipy.reconst.qtdmri.qtdmri_signal_matrix(radial_order, time_order, us, ut, q, tau)[source]#
Constructs the design matrix as a product of 3 separated radial,
angular and temporal design matrices. It precomputes the relevant basis
orders for each one and finally puts them together according to the index
matrix
Constructs the design matrix as a product of 3 separated radial,
angular and temporal design matrices. It precomputes the relevant basis
orders for each one and finally puts them together according to the index
matrix
Constructs the design matrix as a product of 3 separated radial,
angular and temporal design matrices. It precomputes the relevant basis
orders for each one and finally puts them together according to the index
matrix
Computes the cartesian qt-dMRI Laplacian regularization matrix.
If given, uses precomputed matrices for temporal and spatial regularization
matrices to speed up computation. Follows the formulation of Appendix B
in [64].
Computes the spherical qt-dMRI Laplacian regularization matrix.
If given, uses precomputed matrices for temporal and spatial regularization
matrices to speed up computation. Follows the formulation of Appendix C
in [64].
This function visualizes a q-tau acquisition scheme as a function of
gradient strength and pulse separation (big_delta). It represents every
measurements at its G and big_delta position regardless of b-vector, with a
background of b-value isolines for reference. It assumes there is only one
unique pulse length (small_delta) in the acquisition scheme.
Parameters:
gtabGradientTable object
constructed gradient table with big_delta and small_delta given as
inputs.
big_delta_startfloat,
optional minimum big_delta that is plotted in seconds
big_delta_endfloat,
optional maximum big_delta that is plotted in seconds
G_startfloat,
optional minimum gradient strength that is plotted in T/m
G_endfloat,
optional maximum gradient strength that is plotted in T/m
bval_isolinesarray,
optional array of bvalue isolines that are plotted in the background
alpha_shadingfloat between [0-1]
optional shading of the bvalue colors in the background
Generate signals from this model class instance and given parameters.
Parameters:
paramsnumpy.ndarray
Array of shape (…, 28) containing the model parameters. Element 0
is the natural logarithm of the signal without diffusion-weighting,
elements 1-6 are the diffusion tensor elements in Voigt notation,
and elements 7-27 are the covariance tensor elements in Voigt
notation.
Try and make a standard array from an object array
This function takes an object array and attempts to convert it to a more
useful dtype. If array can be converted to a better dtype, Nones are
replaced by fill. To make the behaviour of this function more clear, here
are the most common cases:
obj_arr is an array of scalars of type T. Returns an array like
obj_arr.astype(T)
obj_arr is an array of arrays. All items in obj_arr have the same
shape S. Returns an array with shape obj_arr.shape+S
obj_arr is an array of arrays of different shapes. Returns obj_arr.
Items in obj_arr are not ndarrays or scalars. Returns obj_arr.
Apply a function of two arguments cumulatively to the items of an iterable, from left to right.
This effectively reduces the iterable to a single value. If initial is present,
it is placed before the items of the iterable in the calculation, and serves as
a default when the iterable is empty.
Indices of local maxima from vals given adjacent points
Parameters:
vals(N,) array, dtype np.float64
values at all vertices referred to in either of vertex_inds or
adj_inds’
vertex_inds(V,) array
indices into vals giving vertices that may be local maxima.
adj_indssequence
For every vertex in vertex_inds, the indices (into vals) of
the neighboring points
Returns:
inds(M,) array
Indices into vals giving local maxima of vals, given topology
from adj_inds, and restrictions from vertex_inds. Inds are
returned sorted by value at that index - i.e. smallest value (at
index) first.
Indices of local maxima from vals from count, array neighbors
Parameters:
vals(N,) array, dtype float
values at all vertices referred to in either of vertex_inds or
adj_inds’
vertinds(V,) array, dtype uint32
indices into vals giving vertices that may be local maxima.
adj_counts(V,) array, dtype uint32
For every vertex i in vertex_inds, the number of
neighbors for vertex i
adj_inds(P,) array, dtype uint32
Indices for neighbors for each point. P=sum(adj_counts)
Returns:
inds(M,) array
Indices into vals giving local maxima of vals, given topology
from adj_counts and adj_inds, and restrictions from
vertex_inds. Inds are returned sorted by value at that index -
i.e. smallest value (at index) first.
A point is a local maximum if it is > at least one neighbor and >= all
neighbors. If no points meet the above criteria, 1 maximum is returned such
that odf[maximum] == max(odf).
Remove vertices that are less than theta degrees from any other
Returns vertices that are at least theta degrees from any other vertex.
Vertex v and -v are considered the same so if v and -v are both in
vertices only one is kept. Also if v and w are both in vertices, w must
be separated by theta degrees from both v and -v to be unique. To disable
this, set remove_antipodal to False to keep both directions.
Parameters:
vertices(N, 3) ndarray
N unit vectors.
thetafloat
The minimum separation between vertices in degrees.
return_mapping{False, True}, optional
If True, return mapping as well as vertices and maybe indices
(see below).
return_indices{False, True}, optional
If True, return indices as well as vertices and maybe mapping
(see below).
remove_antipodal{False, True}, optional
If True, v and -v are considered equal, and only one will be kept.
Returns:
unique_vertices(M, 3) ndarray
Vertices sufficiently separated from one another.
mapping(N,) ndarray
For each element vertices[i] (\(i \in 0..N-1\)), the index \(j\) to a
vertex in unique_vertices that is less than theta degrees from
vertices[i]. Only returned if return_mapping is True.
indices(N,) ndarray
indices gives the reverse of mapping. For each element
unique_vertices[j] (\(j \in 0..M-1\)), the index \(i\) to a vertex in
vertices that is less than theta degrees from
unique_vertices[j]. If there is more than one element of
vertices that is less than theta degrees from unique_vertices[j],
return the first (lowest index) matching value. Only return if
return_indices is True.
i in descending array a so a[i] < a[0] * relative_threshold
Call T=a[0]*relative_threshold. Return value i will be the
smallest index in the descending array a such that a[i]<T.
Equivalently, i will be the largest index such that all(a[:i]>=T).
If all values in a are >= T, return the length of array a.
Parameters:
andarray, ndim=1, c-contiguous
Array to be searched. We assume a is in descending order.
relative_thresholdfloat
Applied threshold will be T with T=a[0]*relative_threshold.
Returns:
inp.intp
If T=a[0]*relative_threshold then i will be the largest index
such that all(a[:i]>=T). If all values in a are >= T then
i will be len(a).
Compute signal prediction on model gradient table given given fODF
and GM/CSF volume fractions for each voxel.
Parameters:
gtabGradientTable, optional
The gradients for which the signal will be predicted. Use the
model’s gradient table if None. Default: None
S0ndarray ([x, y, z]) or float, optional
The non diffusion-weighted signal value for each voxel. If a float,
the same value is used for each voxel. If None, 1 is used for
each voxel. Default: None
Fit fODF and GM/CSF volume fractions for a voxel using RUMBA-SD.
Deconvolves the kernel from the diffusion-weighted signal by computing a
maximum likelihood estimation of the fODF
[66]. Minimizes the negative log-likelihood of
the data under Rician or Noncentral Chi noise distributions by adapting the
iterative technique developed in Richardson-Lucy deconvolution.
Parameters:
data1d ndarray (N,)
Signal values for a single voxel.
kernel2d ndarray (N, M)
Deconvolution kernel mapping volume fractions of the M compartments to
N-length signal. Last two columns should be for GM and CSF.
n_iterint, optional
Number of iterations for fODF estimation. Must be a positive int.
Default: 600
recon_type{‘smf’, ‘sos’}, optional
MRI reconstruction method: spatial matched filter (SMF) or
sum-of-squares (SoS). SMF reconstruction generates Rician noise while
SoS reconstruction generates Noncentral Chi noise. Default: ‘smf’
n_coilsint, optional
Number of coils in MRI scanner – only relevant in SoS reconstruction.
Must be a positive int. Default: 1
Returns:
fit_vec1d ndarray (M,)
Vector containing fODF and GM/CSF volume fractions. First M-2
components are fODF while last two are GM and CSF respectively.
Notes
The diffusion MRI signal measured at a given voxel is a sum of
contributions from each intra-voxel compartment, including parallel white
matter (WM) fiber populations in a given orientation as well as effects
from GM and CSF. The equation governing these contributions is:
Where \(S_i\) is the resultant signal along the diffusion-sensitizing
gradient unit vector \(\textbf{v_i}; i = 1, ..., N\) with a b-value of \(b_i\).
\(f_j; j = 1, ..., M\) is the volume fraction of the \(j^{th}\) fiber
population with an anisotropic diffusion tensor \(\textbf{D_j}\).
\(f_{GM}\) and \(f_{CSF}\) are the volume fractions and \(D_{GM}\) and \(D_{CSF}\)
are the mean diffusivities of GM and CSF respectively.
This equation is linear in \(f_j, f_{GM}, f_{CSF}\) and can be simplified to
a single matrix multiplication:
\(\textbf{S} = \textbf{Hf}\)
Where \(\textbf{S}\) is the signal vector at a certain voxel, \(\textbf{H}\) is
the deconvolution kernel, and \(\textbf{f}\) is the vector of volume
fractions for each compartment.
Modern MRI scanners produce noise following a Rician or Noncentral Chi
distribution, depending on their signal reconstruction technique
footcite:p:`Constantinides1997. Using this linear model, it can be shown
that the likelihood of a signal under a Noncentral Chi noise model is:
Where \(S_i\) and \(\bar{S}_i = \textbf{Hf}\) are the measured and expected
signals respectively, and \(n\) is the number of coils in the scanner, and
\(I_{n-1}\) is the modified Bessel function of first kind of order \(n-1\).
This gives the likelihood under a Rician distribution when \(n\) is set to 1.
By taking the negative log of this with respect to \(\textbf{f}\) and setting
the derivative to 0, the \(\textbf{f}\) maximizing likelihood is found to be:
In order to apply this, a reasonable estimate of \(\sigma^2\) is required.
To find this, a separate iterative scheme is found using the derivative
of the negative log with respect to \(\sigma^2\), and is run in parallel.
This is shown here:
Compute kernel mapping orientation densities of white matter fiber
populations (along each vertex of the sphere) and isotropic volume
fractions to a diffusion weighted signal.
Parameters:
gtabGradientTable
Gradient table.
sphereSphere
Sphere with which to sample discrete fiber orientations in order to
construct kernel
wm_response1d ndarray or 2d ndarray or AxSymShResponse, optional
Tensor eigenvalues as a (3,) ndarray, multishell eigenvalues as
a (len(unique_bvals_tolerance(gtab.bvals))-1, 3) ndarray in
order of smallest to largest b-value, or an AxSymShResponse.
gm_responsefloat, optional
Mean diffusivity for GM compartment. If None, then grey
matter compartment set to all zeros.
csf_responsefloat, optional
Mean diffusivity for CSF compartment. If None, then CSF
compartment set to all zeros.
Returns:
kernel2d ndarray (N, M)
Computed kernel; can be multiplied with a vector consisting of volume
fractions for each of M-2 fiber populations as well as GM and CSF
fractions to produce a diffusion weighted signal.
Fit fODF for all voxels simultaneously using RUMBA-SD.
Deconvolves the kernel from the diffusion-weighted signal at each voxel by
computing a maximum likelihood estimation of the fODF
[66]. Global fitting also permits the use of
total variation regularization (RUMBA-SD + TV). The spatial dependence
introduced by TV promotes smoother solutions (i.e. prevents oscillations),
while still allowing for sharp discontinuities [68]. This
promotes smoothness and continuity along individual tracts while preventing
smoothing of adjacent tracts.
Generally, global_fit will proceed more quickly than the voxelwise fit
provided that the computer has adequate RAM (>= 16 GB should be more than
sufficient).
Parameters:
data4d ndarray (x, y, z, N)
Signal values for entire brain. None of the volume dimensions x, y, z
can be 1 if TV regularization is required.
kernel2d ndarray (N, M)
Deconvolution kernel mapping volume fractions of the M compartments to
N-length signal. Last two columns should be for GM and CSF.
mask3d ndarray(x, y, z)
Binary mask specifying voxels of interest with 1; fODF will only be
fit at these voxels (0 elsewhere).
n_iterint, optional
Number of iterations for fODF estimation. Must be a positive int.
recon_type{‘smf’, ‘sos’}, optional
MRI reconstruction method: spatial matched filter (SMF) or
sum-of-squares (SoS). SMF reconstruction generates Rician noise while
SoS reconstruction generates Noncentral Chi noise.
n_coilsint, optional
Number of coils in MRI scanner – only relevant in SoS reconstruction.
Must be a positive int.
use_tvbool, optional
If true, applies total variation regularization. This requires a brain
volume with no singleton dimensions.
verbosebool, optional
If true, logs updates on estimated signal-to-noise ratio after each
iteration.
Returns:
fit_array4d ndarray (x, y, z, M)
fODF and GM/CSF volume fractions computed for each voxel. First M-2
components are fODF, while last two are GM and CSf respectively.
where the first term is the negative log likelihood described in the notes
of rumba_deconv, and the second term is the TV energy, or the sum of
gradient absolute values for the fODF across the entire brain. This results
in a new multiplicative factor in the iterative scheme, now becoming:
(textbf{R}^k)_j = frac{1}{1 - alpha_{TV}divleft(frac{triangledown[
textbf{f}^k_{3D}]_j}{lverttriangledown[textbf{f}^k_{3D}]_j rvert}
right)biggrrvert_{x, y, z}}
Here, \(\triangledown\) is the symbol for the 3D gradient at any voxel.
The regularization strength, \(\alpha_{TV}\) is updated after each iteration
by the discrepancy principle – specifically, it is selected to match the
estimated variance after each iteration [69].
A base-class for the representation of isotropic signals.
The default behavior, suitable for single b-value data is to calculate the
mean in each voxel as an estimate of the signal that does not depend on
direction.
A boolean array used to mark the coordinates in the data that
should be analyzed. Has the shape data.shape[:-1]. Default: None,
which implies that all points should be analyzed.
Predict the isotropic signal, based on a gradient table. In this case,
the prediction will be for an exponential decay with the mean
diffusivity derived from the data that was fit.
Parameters:
gtaba GradientTable class instance, optional
Defaults to use the gtab from the IsotropicModel from which this
fit was derived.
A boolean array used to mark the coordinates in the data that
should be analyzed. Has the shape data.shape[:-1]. Default: None,
which implies that all points should be analyzed.
num_processesint, optional
Split the fit calculation to a pool of children processes using
joblib. This only applies to 4D data arrays. Default is 1,
which does not require joblib and will run fit serially.
If < 0 the maximal number of cores minus num_processes+1
is used (enter -1 to use as many cores as possible).
0 raises an error.
parallel_backend: str, ParallelBackendBase instance or None
Specify the parallelization backend implementation.
Supported backends are:
“loky” used by default, can induce some
communication and memory overhead when exchanging input and
output data with the worker Python processes.
“multiprocessing” previous process-based backend based on
multiprocessing.Pool. Less robust than loky.
“threading” is a very low-overhead backend but it suffers
from the Python Global Interpreter Lock if the called function
relies a lot on Python objects. “threading” is mostly useful
when the execution bottleneck is a compiled extension that
explicitly releases the GIL (for instance a Cython loop wrapped
in a “with nogil” block or an expensive call to a library such
as NumPy).
Sets the rows of the matrix, if the mode is ‘signal’, this should be a
GradientTable. If mode is ‘odf’ this should be a Sphere.
sphereSphere
Sets the columns of the matrix
responselist of 3 elements
The eigenvalues of a tensor which will serve as a kernel
function.
modestr {‘signal’ | ‘odf’}, optional
Choose the (default) ‘signal’ for a design matrix containing predicted
signal in the measurements defined by the gradient table for putative
fascicles oriented along the vertices of the sphere. Otherwise, choose
‘odf’ for an odf convolution matrix, with values of the odf calculated
from a tensor with the provided response eigenvalues, evaluated at the
b-vectors in the gradient table, for the tensors with principal
diffusion directions along the vertices of the sphere.
Returns:
matndarray
A design matrix that can be used for one of the following operations:
when the ‘signal’ mode is used, each column contains the putative
signal in each of the bvectors of the gtab if a fascicle is oriented
in the direction encoded by the sphere vertex corresponding to this
column. This is used for deconvolution with a measured DWI signal. If
the ‘odf’ mode is chosen, each column instead contains the values of
the tensor ODF for a tensor with a principal diffusion direction
corresponding to this vertex. This is used to generate odfs from the
fits of the SFM for the purpose of tracking.
Returns a residual bootstrap sample of the signal_object when indexed
Wraps a signal_object, this signal object can be an interpolator. When
indexed, the wrapper indexes the signal_object to get the signal.
There wrapper than samples the residual bootstrap distribution of signal and
returns that sample.
Spherical harmonics (SH) to rotational harmonics (RH)
Calculate the rotational harmonic decomposition up to
harmonic phase factor m, order l for an axially and antipodally
symmetric function. Note that all m!=0 coefficients
will be ignored as axial symmetry is assumed. Hence, there
will be (sh_order/2+1) non-zero coefficients.
ndarray of SH coefficients for the single fiber response function.
These coefficients must correspond to the real spherical harmonic
functions produced by shm.real_sh_descoteaux_from_index.
m_valuesndarray (N,)
The phase factors (\(m\)) of the spherical harmonic function associated with
each coefficient.
l_valuesndarray (N,)
The orders (\(l\)) of the spherical harmonic function associated with each
coefficient.
Returns:
r_rhndarray ((sh_order+1)*(sh_order+2)/2,)
Rotational harmonics coefficients representing the input r_sh
Generate Dirac delta function orientated in (theta, phi) on the sphere
The spherical harmonics (SH) representation of this Dirac is returned as
coefficients to spherical harmonic functions produced from descoteaux07
basis.
Parameters:
m_valuesndarray (N,)
The phase factors of the spherical harmonic function associated with
each coefficient.
l_valuesndarray (N,)
The order (\(l\)) of the spherical harmonic function associated with each
coefficient.
thetafloat [0, pi]
The polar (colatitudinal) coordinate.
phifloat [0, 2*pi]
The azimuthal (longitudinal) coordinate.
legacy: bool, optional
If true, uses DIPY’s legacy descoteaux07 implementation (where \(|m|\)
is used for m < 0). Else, implements the basis as defined in
Descoteaux et al. 2007 (without the absolute value).
Returns:
diracndarray
SH coefficients representing the Dirac function. The shape of this is
(m + 2) * (m + 1) / 2.
This may take scalar or array arguments. The inputs will be broadcast
against each other.
Parameters:
m_valuesarray of int |m|<=l
The phase factors (\(m\)) of the harmonics.
l_valuesarray of int l>=0
The orders (\(l\)) of the harmonics.
thetafloat [0, 2*pi]
The azimuthal (longitudinal) coordinate.
phifloat [0, pi]
The polar (colatitudinal) coordinate.
use_scipybool, optional
If True, use scipy implementation.
Returns:
y_mncomplex float
The harmonic \(Y_l^m\) sampled at theta and phi.
Notes
This is a faster implementation of scipy.special.sph_harm for
scipy version < 0.15.0. For scipy 0.15 and onwards, we use the scipy
implementation of the function.
The usual definitions for theta`and`phi used in DIPY are interchanged
in the method definition to agree with the definitions in
scipy.special.sph_harm, where theta represents the azimuthal coordinate
and phi represents the polar coordinate.
Although scipy uses a naming convention where m is the order and n
is the degree of the SH, the opposite of DIPY’s, their definition for
both parameters is the same as ours, with l>=0 and |m|<=l.
The definition adopted here follows [7], where the
real harmonic \(Y_l^m\) is defined to be:
Y_l^m =
\begin{cases}
\sqrt{2} * \Im(Y_l^m) \; if m > 0 \\
Y^0_l \; if m = 0 \\
\sqrt{2} * \Re(Y_l^m) \; if m < 0 \\
\end{cases}
This may take scalar or array arguments. The inputs will be broadcast
against each other.
Parameters:
m_valuesarray of int |m|<=l
The phase factors (\(m\)) of the harmonics.
l_valuesarray of int l>=0
The orders (\(l\)) of the harmonics.
thetafloat [0, pi]
The polar (colatitudinal) coordinate.
phifloat [0, 2*pi]
The azimuthal (longitudinal) coordinate.
legacy: bool, optional
If true, uses DIPY’s legacy descoteaux07 implementation (where \(|m|\)
is used for m < 0). Else, implements the basis as defined in
Descoteaux et al. 2007 (without the absolute value).
Returns:
real_shreal float
The real harmonic \(Y_l^m\) sampled at theta and phi.
Samples a real symmetric spherical harmonic basis at point on the sphere
dipy.reconst.shm.real_sym_sh_basis is deprecated, Please use dipy.reconst.shm.real_sh_descoteaux instead
deprecated from version: 1.3
Will raise <class ‘dipy.utils.deprecator.ExpiredDeprecationError’> as of version: 2.0
Samples the basis functions up to order sh_order_max at points on the
sphere given by theta and phi. The basis functions are defined here the
same way as in [7] where the real harmonic \(Y_l^m\)
is defined to be:
Y_l^m =
\begin{cases}
\sqrt{2} * \Im(Y_l^m) \; if m > 0 \\
Y^0_l \; if m = 0 \\
\sqrt{2} * \Im(Y_l^{|m|}) \; if m < 0 \\
\end{cases}
This may take scalar or array arguments. The inputs will be broadcast
against each other.
Parameters:
sh_order_maxint
The maximum order (\(l\)) of the spherical harmonic basis. Even int > 0,
max spherical harmonic order
thetafloat [0, 2*pi]
The azimuthal (longitudinal) coordinate.
phifloat [0, pi]
The polar (colatitudinal) coordinate.
Returns:
y_mnreal float
The real harmonic \(Y_l^m\) sampled at theta and phi
Returns the order (l) and phase_factor (m) of all the symmetric
spherical harmonics of order less then or equal to sh_order_max.
The results, m_list and l_list are kx1 arrays, where k depends on
sh_order_max.
They can be passed to real_sh_descoteaux_from_index() and
:func:real_sh_tournier_from_index.
Parameters:
sh_order_maxint
The maximum order (\(l\)) of the spherical harmonic basis.
Even int > 0, max order to return
In the literature this inverse is often written \((B^{T}B+L^{2})^{-1}B^{T}\).
However here this inverse is implemented using the pseudo-inverse because
it is more numerically stable than the direct implementation of the matrix
product.
Applies the Residual Bootstraps to the data given H and R
data must be normalized, ie 0 < data <= 1
This function, and the bootstrap_data_voxel function, calculate
residual-bootstrap samples given a Hat matrix and a Residual matrix. These
samples can be used for non-parametric statistics or for bootstrap
probabilistic tractography,
Maximum SH order (l) in the SH fit. For sh_order_max, there will be
(sh_order_max+1)*(sh_order_max+2)/2 SH coefficients for a
symmetric basis and (sh_order_max+1)*(sh_order_max+1)
coefficients for a full SH basis.
None for the default DIPY basis,
tournier07 for the Tournier 2007 [13],
[70] basis,
descoteaux07 for the Descoteaux 2007 [7]
basis,
(None defaults to descoteaux07).
full_basis: bool, optional
True for using a SH basis containing even and odd order SH functions.
False for using a SH basis consisting only of even order SH functions.
legacy: bool, optional
True to use a legacy basis definition for backward compatibility
with previous tournier07 and descoteaux07 implementations.
Spherical harmonics (SH) to spherical function (SF).
Parameters:
shndarray
SH coefficients representing a spherical function.
sphereSphere
The points on which to sample the spherical function.
sh_order_maxint, optional
Maximum SH order (l) in the SH fit. For sh_order_max, there will be
(sh_order_max+1)*(sh_order_max+2)/2 SH coefficients for a
symmetric basis and (sh_order_max+1)*(sh_order_max+1)
coefficients for a full SH basis.
None for the default DIPY basis,
tournier07 for the Tournier 2007 [13],
[70] basis,
descoteaux07 for the Descoteaux 2007 [7]
basis,
(None defaults to descoteaux07).
full_basis: bool, optional
True to use a SH basis containing even and odd order SH functions.
Else, use a SH basis consisting only of even order SH functions.
legacy: bool, optional
True to use a legacy basis definition for backward compatibility
with previous tournier07 and descoteaux07 implementations.
Matrix that transforms Spherical harmonics (SH) to spherical
function (SF).
Parameters:
sphereSphere
The points on which to sample the spherical function.
sh_order_maxint, optional
Maximum SH order in the SH fit. For sh_order_max, there will be
(sh_order_max+1)*(sh_order_max+2)/2 SH coefficients for a
symmetric basis and (sh_order_max+1)*(sh_order_max+1)
coefficients for a full SH basis.
None for the default DIPY basis,
tournier07 for the Tournier 2007 [13],
[70] basis,
descoteaux07 for the Descoteaux 2007 [7]
basis,
(None defaults to descoteaux07).
full_basis: bool, optional
If True, uses a SH basis containing even and odd order SH functions.
Else, uses a SH basis consisting only of even order SH functions.
legacy: bool, optional
True to use a legacy basis definition for backward compatibility
with previous tournier07 and descoteaux07 implementations.
return_invbool, optional
If True then the inverse of the matrix is also returned.
smoothfloat, optional
Lambda-regularization in the SH fit.
Returns:
Bndarray
Matrix that transforms spherical harmonics to spherical function
sf=np.dot(sh,B).
Where the last dimension, C, is made of a flattened array of \(l`x:math:`m\)
coefficients, where \(l\) are the SH orders, and \(m = 2l+1\),
So l=1 has 1 coefficient, l=2 has 5, … l=8 has 17 and so on.
A l=2 SH coefficient matrix will then be composed of a IxJxKx6 volume.
The power, \(n\) is usually set to \(n=2\).
The final AP image is then shifted by -log(norm_factor), to be strictly
non-negative. Remaining values < 0 are discarded (set to 0), per default,
and this option is controlled through the non_negative keyword argument.
Convert SH coefficients between legacy-descoteaux07 and tournier07.
Convert SH coefficients between the legacy descoteaux07 SH basis and
the non-legacy tournier07 SH basis. Because this conversion is equal to
its own inverse, it can be used to convert in either direction:
legacy-descoteaux to non-legacy-tournier or non-legacy-tournier to
legacy-descoteaux.
This can be used to convert SH representations between DIPY and MRtrix3.
See [7] and [70] for the
origin of these SH bases.
See [mrtrixbasis] for a description of the basis used in MRtrix3.
See [mrtrixdipybases] for more details on the conversion.
Parameters:
sh_coeffs: ndarray
A ndarray where the last dimension is the
SH coefficients estimates for that voxel.
Returns:
out_sh_coeffs: ndarray
The array of coefficients expressed in the “other” SH basis. If the
input was in the legacy-descoteaux basis then the output will be in the
non-legacy-tournier basis, and vice versa.
where \(\mathbf{q}\) is the wave vector which corresponds to different
gradient directions. Numerous continuous functions \(\phi_i\) can be used to
model \(S\). Some are presented in [76],
[77], and [78].
From the \(c_i\) coefficients, there exist analytical formulae to estimate
the ODF, the return to the origin probability (RTOP), the mean square
displacement (MSD), amongst others [58].
Diffusion propagator on a given set of real points.
if the array r_points is non writeable, then intermediate
results are cached for faster recalculation
Applies the analytical FFT on \(S\) to generate the diffusion
propagator. This is calculated on a discrete 3D grid in order to
obtain an EAP similar to that which is obtained with DSI.
The predicted signal, given a previous fit.
Has the same shape as data.
design_matrixndarray (g, …)
Design matrix holding the covariants used to solve for the
regression coefficients.
leveragesndarray
The leverages (diagonal of the ‘hat matrix’) of the fit.
idxint
The current iteration number.
total_idxint
The total number of iterations.
last_robustndarray
True for inlier indices and False for outlier indices. Must have the
same shape as data.
m_eststr, optional.
M-estimator weighting scheme to use. Currently,
‘gm’ (Geman-McClure) and ‘cauchy’ are provided.
cutofffloat, optional
Cut-off value for defining outliers based on fitting residuals.
Will be passed to the outlier_condition_func.
Typical example: |residuals|>cut_offxstandard_deviation
Robust weights are calculated specifically for the WLS problem, i.e. the
usual form of the WLS problem is accounted for when defining these new
weights, see [79]. On the second-to-last iteration,
OLS is performed without outliers. On the last iteration, WLS is performed
without outliers.
The predicted signal, given a previous fit.
Has the same shape as data.
design_matrixndarray (g, …)
Design matrix holding the covariants used to solve for the
regression coefficients.
leveragesndarray
The leverages (diagonal of the ‘hat matrix’) of the fit.
idxint
The current iteration number.
total_idxint
The total number of iterations.
last_robustndarray
True for inlier indices and False for outlier indices. Must have the
same shape as data.
m_eststr, optional.
M-estimator weighting scheme to use. Currently,
‘gm’ (Geman-McClure) and ‘cauchy’ are provided.
cutofffloat, optional
Cut-off value for defining outliers based on fitting residuals.
Will be passed to the outlier_condition_func.
Typical example: |residuals|>cut_offxstandard_deviation