reconst
bench ([label, verbose, extra_argv])
|
Run benchmarks for module using nose. |
test ([label, verbose, extra_argv, doctests, ...])
|
Run tests for module using nose. |
Module: reconst.base
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).
ReconstFit (model, data)
|
Abstract class which holds the fit result of ReconstModel |
ReconstModel (gtab)
|
Abstract class for signal reconstruction models |
Module: reconst.benchmarks
Module: reconst.benchmarks.bench_bounding_box
Benchmarks for bounding_box
Run all benchmarks with:
import dipy.reconst as dire
dire.bench()
With Pytest, Run this benchmark with:
pytest -svv -c bench.ini /path/to/bench_bounding_box.py
bench_bounding_box ()
|
|
measure (code_str[, times, label])
|
Return elapsed time for executing code in the namespace of the caller. |
Module: reconst.benchmarks.bench_csd
Module: reconst.benchmarks.bench_peaks
Benchmarks for peak finding
Run all benchmarks with:
import dipy.reconst as dire
dire.bench()
With Pytest, Run this benchmark with:
pytest -svv -c bench.ini /path/to/bench_peaks.py
bench_local_maxima ()
|
|
measure (code_str[, times, label])
|
Return elapsed time for executing code in the namespace of the caller. |
Module: reconst.benchmarks.bench_squash
Benchmarks for fast squashing
Run all benchmarks with:
import dipy.reconst as dire
dire.bench()
With Pytest, Run this benchmark with:
pytest -svv -c bench.ini /path/to/bench_squash.py
bench_quick_squash ()
|
|
measure (code_str[, times, label])
|
Return elapsed time for executing code in the namespace of the caller. |
ndindex (shape)
|
An N-dimensional iterator object to index arrays. |
old_squash (arr[, mask, fill])
|
Try and make a standard array from an object array |
reduce (function, sequence[, initial])
|
Apply a function of two arguments cumulatively to the items of a sequence, from left to right, so as to reduce the sequence to a single value. |
Module: reconst.benchmarks.bench_vec_val_sum
Benchmarks for vec / val summation routine
Run benchmarks with:
import dipy.reconst as dire
dire.bench()
With Pytest, Run this benchmark with:
pytest -svv -c bench.ini /path/to/bench_vec_val_sum.py
bench_vec_val_vect ()
|
|
measure (code_str[, times, label])
|
Return elapsed time for executing code in the namespace of the caller. |
Module: reconst.cache
Cache ()
|
Cache values based on a key object (such as a sphere or gradient table). |
auto_attr (func)
|
Decorator to create OneTimeProperty attributes. |
Module: reconst.cross_validation
Cross-validation analysis of diffusion models.
coeff_of_determination (data, model[, axis])
|
Calculate the coefficient of determination for a model prediction, |
kfold_xval (model, data, folds, *model_args, ...)
|
Perform k-fold cross-validation. |
Module: reconst.csdeconv
AxSymShResponse (S0, dwi_response[, bvalue])
|
A simple wrapper for response functions represented using only axially symmetric, even spherical harmonic functions (ie, m == 0 and n even). |
ConstrainedSDTModel (gtab, ratio[, ...])
|
Methods
|
ConstrainedSphericalDeconvModel (gtab, response)
|
Methods
|
SphHarmFit (model, shm_coef, mask)
|
Diffusion data fit to a spherical harmonic model |
SphHarmModel (gtab)
|
To be subclassed by all models that return a SphHarmFit when fit. |
TensorModel (gtab[, fit_method, return_S0_hat])
|
Diffusion Tensor |
auto_response (gtab, data[, roi_center, ...])
|
Automatic estimation of ssst response function using FA. |
auto_response_ssst (gtab, data[, roi_center, ...])
|
Automatic estimation of single-shell single-tissue (ssst) response |
cart2sphere (x, y, z)
|
Return angles for Cartesian 3D coordinates x, y, and z |
csdeconv (dwsignal, X, B_reg[, tau, ...])
|
Constrained-regularized spherical deconvolution (CSD) [1] |
deprecate_with_version (message[, since, ...])
|
Return decorator function function for deprecation warning / error. |
deprecated_params (old_name[, new_name, ...])
|
Deprecate a renamed or removed function argument. |
estimate_response (gtab, evals, S0)
|
Estimate single fiber response function |
fa_trace_to_lambdas ([fa, trace])
|
|
forward_sdeconv_mat (r_rh, n)
|
Build forward spherical deconvolution matrix |
forward_sdt_deconv_mat (ratio, n[, r2_term])
|
Build forward sharpening deconvolution transform (SDT) matrix |
fractional_anisotropy (evals[, axis])
|
Return Fractional anisotropy (FA) of a diffusion tensor. |
get_sphere ([name])
|
provide triangulated spheres |
lazy_index (index)
|
Produces a lazy index |
lpn (n, z)
|
Legendre function of the first kind. |
mask_for_response_ssst (gtab, data[, ...])
|
Computation of mask for single-shell single-tissue (ssst) response |
multi_voxel_fit (single_voxel_fit)
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
ndindex (shape)
|
An N-dimensional iterator object to index arrays. |
odf_deconv (odf_sh, R, B_reg[, lambda_, tau, ...])
|
ODF constrained-regularized spherical deconvolution using the Sharpening Deconvolution Transform (SDT) [1], [2]. |
odf_sh_to_sharp (odfs_sh, sphere[, basis, ...])
|
Sharpen odfs using the sharpening deconvolution transform [2] |
peaks_from_model (model, data, sphere, ...[, ...])
|
Fit the model to data and computes peaks and metrics |
quad (func, a, b[, args, full_output, ...])
|
Compute a definite integral. |
real_sh_descoteaux (sh_order, theta, phi[, ...])
|
Compute real spherical harmonics as in Descoteaux et al. 2007 [Ra6d8f6cd2652-1], where the real harmonic \(Y^m_n\) is defined to be:. |
real_sh_descoteaux_from_index (m, n, theta, phi)
|
Compute real spherical harmonics as in Descoteaux et al. 2007 [R2763a8caef3c-1], where the real harmonic \(Y^m_n\) is defined to be:. |
recursive_response (gtab, data[, mask, ...])
|
Recursive calibration of response function using peak threshold |
response_from_mask (gtab, data, mask)
|
Computation of single-shell single-tissue (ssst) response |
response_from_mask_ssst (gtab, data, mask)
|
Computation of single-shell single-tissue (ssst) response |
sh_to_rh (r_sh, m, n)
|
Spherical harmonics (SH) to rotational harmonics (RH) |
single_tensor (gtab[, S0, evals, evecs, snr])
|
Simulate diffusion-weighted signals with a single tensor. |
sph_harm_ind_list (sh_order[, full_basis])
|
Returns the degree (m ) and order (n ) of all the symmetric spherical harmonics of degree less then or equal to sh_order . |
vec2vec_rotmat (u, v)
|
rotation matrix from 2 unit vectors |
Module: reconst.dki
Classes and functions for fitting the diffusion kurtosis model
DiffusionKurtosisFit (model, model_params)
|
Class for fitting the Diffusion Kurtosis Model |
DiffusionKurtosisModel (gtab[, fit_method])
|
Class for the Diffusion Kurtosis Model |
ReconstModel (gtab)
|
Abstract class for signal reconstruction models |
TensorFit (model, model_params[, model_S0])
|
- Attributes:
|
Wcons (k_elements)
|
Construct the full 4D kurtosis tensors from its 15 independent elements |
Wrotate (kt, Basis)
|
Rotate a kurtosis tensor from the standard Cartesian coordinate system to another coordinate system basis |
Wrotate_element (kt, indi, indj, indk, indl, B)
|
Computes the the specified index element of a kurtosis tensor rotated to the coordinate system basis B. |
apparent_kurtosis_coef (dki_params, sphere[, ...])
|
Calculates the apparent kurtosis coefficient (AKC) in each direction of a sphere [1]. |
axial_kurtosis (dki_params[, min_kurtosis, ...])
|
Computes axial Kurtosis (AK) from the kurtosis tensor [1], [2]. |
carlson_rd (x, y, z[, errtol])
|
Computes the Carlson's incomplete elliptic integral of the second kind defined as: |
carlson_rf (x, y, z[, errtol])
|
Computes the Carlson's incomplete elliptic integral of the first kind defined as: |
cart2sphere (x, y, z)
|
Return angles for Cartesian 3D coordinates x, y, and z |
check_multi_b (gtab, n_bvals[, non_zero, bmag])
|
Check if you have enough different b-values in your gradient table |
decompose_tensor (tensor[, min_diffusivity])
|
Returns eigenvalues and eigenvectors given a diffusion tensor |
design_matrix (gtab)
|
Construct B design matrix for DKI. |
directional_diffusion (dt, V[, min_diffusivity])
|
Calculates the apparent diffusion coefficient (adc) in each direction of a sphere for a single voxel [1]. |
directional_diffusion_variance (kt, V[, ...])
|
Calculates the apparent diffusion variance (adv) in each direction of a sphere for a single voxel [1]. |
directional_kurtosis (dt, md, kt, V[, ...])
|
Calculates the apparent kurtosis coefficient (akc) in each direction of a sphere for a single voxel [1]. |
dki_prediction (dki_params, gtab[, S0])
|
Predict a signal given diffusion kurtosis imaging parameters. |
from_lower_triangular (D)
|
Returns a tensor given the six unique tensor elements |
get_fnames ([name])
|
Provide full paths to example or test datasets. |
get_sphere ([name])
|
provide triangulated spheres |
kurtosis_fractional_anisotropy (dki_params)
|
Computes the anisotropy of the kurtosis tensor (KFA) [1]. |
kurtosis_maximum (dki_params[, sphere, gtol, ...])
|
Computes kurtosis maximum value |
local_maxima
|
Local maxima of a function evaluated on a discrete set of points. |
lower_triangular (tensor[, b0])
|
Returns the six lower triangular values of the tensor ordered as (Dxx, Dxy, Dyy, Dxz, Dyz, Dzz) and a dummy variable if b0 is not None. |
mean_diffusivity (evals[, axis])
|
Mean Diffusivity (MD) of a diffusion tensor. |
mean_kurtosis (dki_params[, min_kurtosis, ...])
|
Computes mean Kurtosis (MK) from the kurtosis tensor. |
mean_kurtosis_tensor (dki_params[, ...])
|
Computes mean of the kurtosis tensor (MKT) [1]. |
ndindex (shape)
|
An N-dimensional iterator object to index arrays. |
nlls_fit_tensor (design_matrix, data[, ...])
|
Fit the cumulant expansion params (e.g. |
ols_fit_dki (design_matrix, data)
|
Computes the diffusion and kurtosis tensors using an ordinary linear least squares (OLS) approach . |
perpendicular_directions (v[, num, half])
|
Computes n evenly spaced perpendicular directions relative to a given vector v |
radial_kurtosis (dki_params[, min_kurtosis, ...])
|
Radial Kurtosis (RK) of a diffusion kurtosis tensor [1], [2]. |
restore_fit_tensor (design_matrix, data[, ...])
|
Use the RESTORE algorithm [1] to calculate a robust tensor fit |
sphere2cart (r, theta, phi)
|
Spherical to Cartesian coordinates |
split_dki_param (dki_params)
|
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 |
vec_val_vect
|
Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs |
wls_fit_dki (design_matrix, data)
|
Computes the diffusion and kurtosis tensors using a weighted linear least squares (WLS) approach . |
Module: reconst.dki_micro
Classes and functions for fitting the DKI-based microstructural model
DiffusionKurtosisFit (model, model_params)
|
Class for fitting the Diffusion Kurtosis Model |
DiffusionKurtosisModel (gtab[, fit_method])
|
Class for the Diffusion Kurtosis Model |
KurtosisMicrostructuralFit (model, model_params)
|
Class for fitting the Diffusion Kurtosis Microstructural Model |
KurtosisMicrostructureModel (gtab[, fit_method])
|
Class for the Diffusion Kurtosis Microstructural Model |
axial_diffusivity (evals[, axis])
|
Axial Diffusivity (AD) of a diffusion tensor. |
axonal_water_fraction (dki_params[, sphere, ...])
|
Computes the axonal water fraction from DKI [1]. |
decompose_tensor (tensor[, min_diffusivity])
|
Returns eigenvalues and eigenvectors given a diffusion tensor |
diffusion_components (dki_params[, sphere, ...])
|
Extracts the restricted and hindered diffusion tensors of well aligned fibers from diffusion kurtosis imaging parameters [1]. |
directional_diffusion (dt, V[, min_diffusivity])
|
Calculates the apparent diffusion coefficient (adc) in each direction of a sphere for a single voxel [1]. |
directional_kurtosis (dt, md, kt, V[, ...])
|
Calculates the apparent kurtosis coefficient (akc) in each direction of a sphere for a single voxel [1]. |
dkimicro_prediction (params, gtab[, S0])
|
Signal prediction given the DKI microstructure model parameters. |
dti_design_matrix (gtab[, dtype])
|
Constructs design matrix for DTI weighted least squares or least squares fitting. |
from_lower_triangular (D)
|
Returns a tensor given the six unique tensor elements |
get_sphere ([name])
|
provide triangulated spheres |
kurtosis_maximum (dki_params[, sphere, gtol, ...])
|
Computes kurtosis maximum value |
lower_triangular (tensor[, b0])
|
Returns the six lower triangular values of the tensor ordered as (Dxx, Dxy, Dyy, Dxz, Dyz, Dzz) and a dummy variable if b0 is not None. |
mean_diffusivity (evals[, axis])
|
Mean Diffusivity (MD) of a diffusion tensor. |
ndindex (shape)
|
An N-dimensional iterator object to index arrays. |
radial_diffusivity (evals[, axis])
|
Radial Diffusivity (RD) of a diffusion tensor. |
split_dki_param (dki_params)
|
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 |
tortuosity (hindered_ad, hindered_rd)
|
Computes the tortuosity of the hindered diffusion compartment given its axial and radial diffusivities |
trace (evals[, axis])
|
Trace of a diffusion tensor. |
vec_val_vect
|
Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs |
Module: reconst.dsi
Cache ()
|
Cache values based on a key object (such as a sphere or gradient table). |
DiffusionSpectrumDeconvFit (model, data)
|
Methods
|
DiffusionSpectrumDeconvModel (gtab[, ...])
|
Methods
|
DiffusionSpectrumFit (model, data)
|
Methods
|
DiffusionSpectrumModel (gtab[, qgrid_size, ...])
|
Methods
|
OdfFit (model, data)
|
Methods
|
OdfModel (gtab)
|
An abstract class to be sub-classed by specific odf models |
LR_deconv (prop, psf[, numit, acc_factor])
|
Perform Lucy-Richardson deconvolution algorithm on a 3D array. |
create_qspace (gtab, origin)
|
create the 3D grid which holds the signal values (q-space) |
create_qtable (gtab, origin)
|
create a normalized version of gradients |
fftn (x[, shape, axes, overwrite_x])
|
Return multidimensional discrete Fourier transform. |
fftshift (x[, axes])
|
Shift the zero-frequency component to the center of the spectrum. |
gen_PSF (qgrid_sampling, siz_x, siz_y, siz_z)
|
Generate a PSF for DSI Deconvolution by taking the ifft of the binary q-space sampling mask and truncating it to keep only the center. |
half_to_full_qspace (data, gtab)
|
Half to full Cartesian grid mapping |
hanning_filter (gtab, filter_width, origin)
|
create a hanning window |
ifftshift (x[, axes])
|
The inverse of fftshift. |
map_coordinates (input, coordinates[, ...])
|
Map the input array to new coordinates by interpolation. |
multi_voxel_fit (single_voxel_fit)
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
pdf_interp_coords (sphere, rradius, origin)
|
Precompute coordinates for ODF calculation from the PDF |
pdf_odf (Pr, rradius, interp_coords)
|
Calculates the real ODF from the diffusion propagator(PDF) Pr |
project_hemisph_bvecs (gtab)
|
Project any near identical bvecs to the other hemisphere |
threshold_propagator (P[, estimated_snr])
|
Applies hard threshold on the propagator to remove background noise for the deconvolution. |
Module: reconst.dti
Classes and functions for fitting tensors
ReconstModel (gtab)
|
Abstract class for signal reconstruction models |
TensorFit (model, model_params[, model_S0])
|
- Attributes:
|
TensorModel (gtab[, fit_method, return_S0_hat])
|
Diffusion Tensor |
Version (version)
|
This class abstracts handling of a project's versions. |
apparent_diffusion_coef (q_form, sphere)
|
Calculate the apparent diffusion coefficient (ADC) in each direction of a sphere. |
auto_attr (func)
|
Decorator to create OneTimeProperty attributes. |
axial_diffusivity (evals[, axis])
|
Axial Diffusivity (AD) of a diffusion tensor. |
color_fa (fa, evecs)
|
Color fractional anisotropy of diffusion tensor |
decompose_tensor (tensor[, min_diffusivity])
|
Returns eigenvalues and eigenvectors given a diffusion tensor |
design_matrix (gtab[, dtype])
|
Constructs design matrix for DTI weighted least squares or least squares fitting. |
determinant (q_form)
|
The determinant of a tensor, given in quadratic form |
deviatoric (q_form)
|
Calculate the deviatoric (anisotropic) part of the tensor [1]. |
eig_from_lo_tri (data[, min_diffusivity])
|
Calculates tensor eigenvalues/eigenvectors from an array containing the lower diagonal form of the six unique tensor elements. |
fractional_anisotropy (evals[, axis])
|
Return Fractional anisotropy (FA) of a diffusion tensor. |
from_lower_triangular (D)
|
Returns a tensor given the six unique tensor elements |
geodesic_anisotropy (evals[, axis])
|
Geodesic anisotropy (GA) of a diffusion tensor. |
get_sphere ([name])
|
provide triangulated spheres |
gradient_table (bvals[, bvecs, big_delta, ...])
|
A general function for creating diffusion MR gradients. |
isotropic (q_form)
|
Calculate the isotropic part of the tensor [1]. |
iter_fit_tensor ([step])
|
Wrap a fit_tensor func and iterate over chunks of data with given length |
linearity (evals[, axis])
|
The linearity of the tensor [1] |
lower_triangular (tensor[, b0])
|
Returns the six lower triangular values of the tensor ordered as (Dxx, Dxy, Dyy, Dxz, Dyz, Dzz) and a dummy variable if b0 is not None. |
mean_diffusivity (evals[, axis])
|
Mean Diffusivity (MD) of a diffusion tensor. |
mode (q_form)
|
Mode (MO) of a diffusion tensor [1]. |
nlls_fit_tensor (design_matrix, data[, ...])
|
Fit the cumulant expansion params (e.g. |
norm (q_form)
|
Calculate the Frobenius norm of a tensor quadratic form |
ols_fit_tensor (design_matrix, data[, ...])
|
Computes ordinary least squares (OLS) fit to calculate self-diffusion tensor using a linear regression model [Rd310240b4eed-1]. |
pinv (a[, rcond])
|
Vectorized version of numpy.linalg.pinv |
planarity (evals[, axis])
|
The planarity of the tensor [1] |
quantize_evecs (evecs[, odf_vertices])
|
Find the closest orientation of an evenly distributed sphere |
radial_diffusivity (evals[, axis])
|
Radial Diffusivity (RD) of a diffusion tensor. |
restore_fit_tensor (design_matrix, data[, ...])
|
Use the RESTORE algorithm [1] to calculate a robust tensor fit |
sphericity (evals[, axis])
|
The sphericity of the tensor [1] |
tensor_prediction (dti_params, gtab, S0)
|
Predict a signal given tensor parameters. |
trace (evals[, axis])
|
Trace of a diffusion tensor. |
vec_val_vect
|
Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs |
vector_norm (vec[, axis, keepdims])
|
Return vector Euclidean (L2) norm |
wls_fit_tensor (design_matrix, data[, ...])
|
Computes weighted least squares (WLS) fit to calculate self-diffusion tensor using a linear regression model [1]. |
Module: reconst.eudx_direction_getter
Module: reconst.forecast
Cache ()
|
Cache values based on a key object (such as a sphere or gradient table). |
ForecastFit (model, data, sh_coef, d_par, d_perp)
|
- Attributes:
|
ForecastModel (gtab[, sh_order, lambda_lb, ...])
|
Fiber ORientation Estimated using Continuous Axially Symmetric Tensors (FORECAST) [1,2,3]_. |
OdfFit (model, data)
|
Methods
|
OdfModel (gtab)
|
An abstract class to be sub-classed by specific odf models |
Version (version)
|
This class abstracts handling of a project's versions. |
cart2sphere (x, y, z)
|
Return angles for Cartesian 3D coordinates x, y, and z |
csdeconv (dwsignal, X, B_reg[, tau, ...])
|
Constrained-regularized spherical deconvolution (CSD) [1] |
find_signal_means (b_unique, data_norm, ...)
|
Calculate the mean signal for each shell. |
forecast_error_func (x, b_unique, E)
|
Calculates the difference between the mean signal calculated using the parameter vector x and the average signal E using FORECAST and SMT |
forecast_matrix (sh_order, d_par, d_perp, bvals)
|
Compute the FORECAST radial matrix |
lb_forecast (sh_order)
|
Returns the Laplace-Beltrami regularization matrix for FORECAST |
leastsq (func, x0[, args, Dfun, full_output, ...])
|
Minimize the sum of squares of a set of equations. |
multi_voxel_fit (single_voxel_fit)
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
optional_package (name[, trip_msg])
|
Return package-like thing and module setup for package name |
psi_l (l, b)
|
|
real_sh_descoteaux_from_index (m, n, theta, phi)
|
Compute real spherical harmonics as in Descoteaux et al. 2007 [Rded9a56186a6-1], where the real harmonic \(Y^m_n\) is defined to be:. |
rho_matrix (sh_order, vecs)
|
Compute the SH matrix \(\rho\) |
warn (/, message[, category, stacklevel, source])
|
Issue a warning, or maybe ignore it or raise an exception. |
Module: reconst.fwdti
Classes and functions for fitting tensors without free water
contamination
FreeWaterTensorFit (model, model_params)
|
Class for fitting the Free Water Tensor Model |
FreeWaterTensorModel (gtab[, fit_method])
|
Class for the Free Water Elimination Diffusion Tensor Model |
ReconstModel (gtab)
|
Abstract class for signal reconstruction models |
TensorFit (model, model_params[, model_S0])
|
- Attributes:
|
check_multi_b (gtab, n_bvals[, non_zero, bmag])
|
Check if you have enough different b-values in your gradient table |
cholesky_to_lower_triangular (R)
|
Convert Cholesky decompostion elements to the diffusion tensor elements |
decompose_tensor (tensor[, min_diffusivity])
|
Returns eigenvalues and eigenvectors given a diffusion tensor |
design_matrix (gtab[, dtype])
|
Constructs design matrix for DTI weighted least squares or least squares fitting. |
from_lower_triangular (D)
|
Returns a tensor given the six unique tensor elements |
fwdti_prediction (params, gtab[, S0, Diso])
|
Signal prediction given the free water DTI model parameters. |
lower_triangular (tensor[, b0])
|
Returns the six lower triangular values of the tensor ordered as (Dxx, Dxy, Dyy, Dxz, Dyz, Dzz) and a dummy variable if b0 is not None. |
lower_triangular_to_cholesky (tensor_elements)
|
Performs Cholesky decomposition of the diffusion tensor |
multi_voxel_fit (single_voxel_fit)
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
ndindex (shape)
|
An N-dimensional iterator object to index arrays. |
nls_fit_tensor (gtab, data[, mask, Diso, ...])
|
Fit the water elimination tensor model using the non-linear least-squares. |
nls_iter (design_matrix, sig, S0[, Diso, ...])
|
Applies non linear least squares fit of the water free elimination model to single voxel signals. |
vec_val_vect
|
Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs |
wls_fit_tensor (gtab, data[, Diso, mask, ...])
|
Computes weighted least squares (WLS) fit to calculate self-diffusion tensor using a linear regression model [1]. |
wls_iter (design_matrix, sig, S0[, Diso, ...])
|
Applies weighted linear least squares fit of the water free elimination model to single voxel signals. |
Module: reconst.gqi
Classes and functions for generalized q-sampling
Cache ()
|
Cache values based on a key object (such as a sphere or gradient table). |
GeneralizedQSamplingFit (model, data)
|
Methods
|
GeneralizedQSamplingModel (gtab[, method, ...])
|
Methods
|
OdfFit (model, data)
|
Methods
|
OdfModel (gtab)
|
An abstract class to be sub-classed by specific odf models |
equatorial_maximum (vertices, odf, pole, width)
|
|
equatorial_zone_vertices (vertices, pole[, width])
|
finds the 'vertices' in the equatorial zone conjugate to 'pole' with width half 'width' degrees |
multi_voxel_fit (single_voxel_fit)
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
normalize_qa (qa[, max_qa])
|
Normalize quantitative anisotropy. |
npa (self, odf[, width])
|
non-parametric anisotropy |
odf_sum (odf)
|
|
patch_maximum (vertices, odf, pole, width)
|
|
patch_sum (vertices, odf, pole, width)
|
|
patch_vertices (vertices, pole, width)
|
find 'vertices' within the cone of 'width' degrees around 'pole' |
polar_zone_vertices (vertices, pole[, width])
|
finds the 'vertices' in the equatorial band around the 'pole' of radius 'width' degrees |
squared_radial_component (x[, tol])
|
Part of the GQI2 integral |
triple_odf_maxima (vertices, odf, width)
|
|
upper_hemi_map (v)
|
maps a 3-vector into the z-upper hemisphere |
Module: reconst.ivim
Classes and functions for fitting ivim model
IvimFit (model, model_params)
|
- Attributes:
|
IvimModelTRR (gtab[, split_b_D, split_b_S0, ...])
|
Ivim model |
IvimModelVP (gtab[, bounds, maxiter, xtol])
|
Methods
|
ReconstModel (gtab)
|
Abstract class for signal reconstruction models |
Version (version)
|
This class abstracts handling of a project's versions. |
IvimModel (gtab[, fit_method])
|
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. |
differential_evolution (func, bounds[, args, ...])
|
Finds the global minimum of a multivariate function. |
f_D_star_error (params, gtab, signal, S0, D)
|
Error function used to fit f and D_star keeping S0 and D fixed |
f_D_star_prediction (params, gtab, S0, D)
|
Function used to predict IVIM signal when S0 and D are known by considering f and D_star as the unknown parameters. |
ivim_model_selector (gtab[, fit_method])
|
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. |
ivim_prediction (params, gtab)
|
The Intravoxel incoherent motion (IVIM) model function. |
least_squares (fun, x0[, jac, bounds, ...])
|
Solve a nonlinear least-squares problem with bounds on the variables. |
multi_voxel_fit (single_voxel_fit)
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
optional_package (name[, trip_msg])
|
Return package-like thing and module setup for package name |
Module: reconst.mapmri
Cache ()
|
Cache values based on a key object (such as a sphere or gradient table). |
MapmriFit (model, mapmri_coef, mu, R, lopt[, ...])
|
- Attributes:
|
MapmriModel (gtab[, radial_order, ...])
|
Mean Apparent Propagator MRI (MAPMRI) [1] of the diffusion signal. |
Optimizer (fun, x0[, args, method, jac, ...])
|
- Attributes:
|
PositiveDefiniteLeastSquares (m[, A, L])
|
Methods
|
ReconstFit (model, data)
|
Abstract class which holds the fit result of ReconstModel |
ReconstModel (gtab)
|
Abstract class for signal reconstruction models |
Version (version)
|
This class abstracts handling of a project's versions. |
b_mat (index_matrix)
|
Calculates the B coefficients from [1] Eq. |
b_mat_isotropic (index_matrix)
|
Calculates the isotropic B coefficients from [1] Fig 8. |
binomialfloat (n, k)
|
Custom Binomial function |
cart2sphere (x, y, z)
|
Return angles for Cartesian 3D coordinates x, y, and z |
create_rspace (gridsize, radius_max)
|
Create the real space table, that contains the points in which to compute the pdf. |
delta (n, m)
|
|
factorial2 (n[, exact])
|
Double factorial. |
gcv_cost_function (weight, args)
|
The GCV cost function that is iterated [4]. |
generalized_crossvalidation (data, M, LR[, ...])
|
Generalized Cross Validation Function [1] eq. |
generalized_crossvalidation_array (data, M, LR)
|
Generalized Cross Validation Function eq. |
genlaguerre (n, alpha[, monic])
|
Generalized (associated) Laguerre polynomial. |
gradient_table (bvals[, bvecs, big_delta, ...])
|
A general function for creating diffusion MR gradients. |
hermite (n[, monic])
|
Physicist's Hermite polynomial. |
isotropic_scale_factor (mu_squared)
|
Estimated isotropic scaling factor _[1] Eq. |
load_sdp_constraints (model_name[, order])
|
Import semidefinite programming constraint matrices for different models, generated as described for example in [1]. |
map_laplace_s (n, m)
|
R(m,n) static matrix for Laplacian regularization [1] eq. |
map_laplace_t (n, m)
|
L(m, n) static matrix for Laplacian regularization [1] eq. |
map_laplace_u (n, m)
|
S(n, m) static matrix for Laplacian regularization [1] eq. |
mapmri_STU_reg_matrices (radial_order)
|
Generate the static portions of the Laplacian regularization matrix according to [1] eq. |
mapmri_index_matrix (radial_order)
|
Calculates the indices for the MAPMRI [1] basis in x, y and z. |
mapmri_isotropic_K_mu_dependent (...)
|
Computes mu dependent part of M. |
mapmri_isotropic_K_mu_independent (...)
|
Computes mu independent part of K. |
mapmri_isotropic_M_mu_dependent (...)
|
Computed the mu dependent part of the signal design matrix. |
mapmri_isotropic_M_mu_independent (...)
|
Computed the mu independent part of the signal design matrix. |
mapmri_isotropic_index_matrix (radial_order)
|
Calculates the indices for the isotropic MAPMRI basis [1] Fig 8. |
mapmri_isotropic_laplacian_reg_matrix (...)
|
Computes the Laplacian regularization matrix for MAP-MRI's isotropic implementation [1] eq. |
mapmri_isotropic_laplacian_reg_matrix_from_index_matrix (...)
|
Computes the Laplacian regularization matrix for MAP-MRI's isotropic implementation [1] eq. |
mapmri_isotropic_odf_matrix (radial_order, ...)
|
Compute the isotropic MAPMRI ODF matrix [1] Eq. |
mapmri_isotropic_odf_sh_matrix (radial_order, ...)
|
Compute the isotropic MAPMRI ODF matrix [1] Eq. |
mapmri_isotropic_phi_matrix (radial_order, mu, q)
|
Three dimensional isotropic MAPMRI signal basis function from [1] Eq. |
mapmri_isotropic_psi_matrix (radial_order, ...)
|
Three dimensional isotropic MAPMRI propagator basis function from [1] Eq. |
mapmri_isotropic_radial_pdf_basis (j, l, mu, r)
|
Radial part of the isotropic 1D-SHORE propagator basis [1] eq. |
mapmri_isotropic_radial_signal_basis (j, l, ...)
|
Radial part of the isotropic 1D-SHORE signal basis [1] eq. |
mapmri_laplacian_reg_matrix (ind_mat, mu, ...)
|
Put the Laplacian regularization matrix together [1] eq. |
mapmri_odf_matrix (radial_order, mu, s, vertices)
|
Compute the MAPMRI ODF matrix [1] Eq. |
mapmri_phi_1d (n, q, mu)
|
One dimensional MAPMRI basis function from [1] Eq. |
mapmri_phi_matrix (radial_order, mu, q_gradients)
|
Compute the MAPMRI phi matrix for the signal [1] eq. |
mapmri_psi_1d (n, x, mu)
|
One dimensional MAPMRI propagator basis function from [1] Eq. |
mapmri_psi_matrix (radial_order, mu, rgrad)
|
Compute the MAPMRI psi matrix for the propagator [1] eq. |
mfactorial (x, /)
|
Find x!. |
multi_voxel_fit (single_voxel_fit)
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
optional_package (name[, trip_msg])
|
Return package-like thing and module setup for package name |
real_sh_descoteaux_from_index (m, n, theta, phi)
|
Compute real spherical harmonics as in Descoteaux et al. 2007 [R85c991cbd5ae-1], where the real harmonic \(Y^m_n\) is defined to be:. |
sfactorial (n[, exact])
|
The factorial of a number or array of numbers. |
sph_harm_ind_list (sh_order[, full_basis])
|
Returns the degree (m ) and order (n ) of all the symmetric spherical harmonics of degree less then or equal to sh_order . |
warn (/, message[, category, stacklevel, source])
|
Issue a warning, or maybe ignore it or raise an exception. |
Module: reconst.mcsd
GradientTable (gradients[, big_delta, ...])
|
Diffusion gradient information |
MSDeconvFit (model, coeff, mask)
|
- Attributes:
|
MultiShellDeconvModel (gtab, response[, ...])
|
Methods
|
MultiShellResponse (response, sh_order, shells)
|
- Attributes:
|
QpFitter (X, reg)
|
Methods
|
TensorModel (gtab[, fit_method, return_S0_hat])
|
Diffusion Tensor |
Version (version)
|
This class abstracts handling of a project's versions. |
auto_response_msmt (gtab, data[, tol, ...])
|
Automatic estimation of multi-shell multi-tissue (msmt) response |
fractional_anisotropy (evals[, axis])
|
Return Fractional anisotropy (FA) of a diffusion tensor. |
get_bval_indices (bvals, bval[, tol])
|
Get indices where the b-value is bval |
gradient_table (bvals[, bvecs, big_delta, ...])
|
A general function for creating diffusion MR gradients. |
mask_for_response_msmt (gtab, data[, ...])
|
Computation of masks for multi-shell multi-tissue (msmt) response |
mean_diffusivity (evals[, axis])
|
Mean Diffusivity (MD) of a diffusion tensor. |
multi_shell_fiber_response (sh_order, bvals, ...)
|
Fiber response function estimation for multi-shell data. |
multi_tissue_basis (gtab, sh_order, iso_comp)
|
Builds a basis for multi-shell multi-tissue CSD model. |
multi_voxel_fit (single_voxel_fit)
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
optional_package (name[, trip_msg])
|
Return package-like thing and module setup for package name |
response_from_mask_msmt (gtab, data, mask_wm, ...)
|
Computation of multi-shell multi-tissue (msmt) response |
response_from_mask_ssst (gtab, data, mask)
|
Computation of single-shell single-tissue (ssst) response |
single_tensor (gtab[, S0, evals, evecs, snr])
|
Simulate diffusion-weighted signals with a single tensor. |
solve_qp (P, Q, G, H)
|
Helper function to set up and solve the Quadratic Program (QP) in CVXPY. |
unique_bvals_tolerance (bvals[, tol])
|
Gives the unique b-values of the data, within a tolerance gap |
Module: reconst.msdki
Classes and functions for fitting the mean signal diffusion kurtosis
model
MeanDiffusionKurtosisFit (model, model_params)
|
- Attributes:
|
MeanDiffusionKurtosisModel (gtab[, bmag, ...])
|
Mean signal Diffusion Kurtosis Model |
ReconstModel (gtab)
|
Abstract class for signal reconstruction models |
auto_attr (func)
|
Decorator to create OneTimeProperty attributes. |
awf_from_msk (msk[, mask])
|
Computes the axonal water fraction from the mean signal kurtosis assuming the 2-compartmental spherical mean technique model [1], [2] |
check_multi_b (gtab, n_bvals[, non_zero, bmag])
|
Check if you have enough different b-values in your gradient table |
design_matrix (ubvals)
|
Constructs design matrix for the mean signal diffusion kurtosis model |
mean_signal_bvalue (data, gtab[, bmag])
|
Computes the average signal across different diffusion directions for each unique b-value |
msdki_prediction (msdki_params, gtab[, S0])
|
Predict the mean signal given the parameters of the mean signal DKI, an GradientTable object and S0 signal. |
msk_from_awf (f)
|
Computes mean signal kurtosis from axonal water fraction estimates of the SMT2 model |
ndindex (shape)
|
An N-dimensional iterator object to index arrays. |
round_bvals (bvals[, bmag])
|
"This function rounds the b-values |
unique_bvals_magnitude (bvals[, bmag, rbvals])
|
This function gives the unique rounded b-values of the data |
wls_fit_msdki (design_matrix, msignal, ng[, ...])
|
Fits the mean signal diffusion kurtosis imaging based on a weighted least square solution [1]. |
Module: reconst.multi_voxel
Tools to easily make multi voxel models
CallableArray
|
An array which can be called like a function |
MultiVoxelFit (model, fit_array, mask)
|
Holds an array of fits and allows access to their attributes and methods |
ReconstFit (model, data)
|
Abstract class which holds the fit result of ReconstModel |
tqdm (*_, **__)
|
Decorate an iterable object, returning an iterator which acts exactly like the original iterable, but prints a dynamically updating progressbar every time a value is requested. |
as_strided (x[, shape, strides, subok, writeable])
|
Create a view into the array with the given shape and strides. |
multi_voxel_fit (single_voxel_fit)
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
ndindex (shape)
|
An N-dimensional iterator object to index arrays. |
Module: reconst.odf
OdfFit (model, data)
|
Methods
|
OdfModel (gtab)
|
An abstract class to be sub-classed by specific odf models |
ReconstFit (model, data)
|
Abstract class which holds the fit result of ReconstModel |
ReconstModel (gtab)
|
Abstract class for signal reconstruction models |
gfa (samples)
|
The general fractional anisotropy of a function evaluated on the unit sphere |
minmax_normalize (samples[, out])
|
Min-max normalization of a function evaluated on the unit sphere |
Module: reconst.qtdmri
Cache ()
|
Cache values based on a key object (such as a sphere or gradient table). |
QtdmriFit (model, qtdmri_coef, us, ut, ...)
|
Methods
|
QtdmriModel (gtab[, radial_order, ...])
|
The q:math:tau-dMRI model [1] to analytically and continuously represent the q:math:tau diffusion signal attenuation over diffusion sensitization q and diffusion time \(\tau\). |
Version (version)
|
This class abstracts handling of a project's versions. |
GCV_cost_function (weight, arguments)
|
Generalized Cross Validation Function that is iterated [1]. |
H (value)
|
Step function of H(x)=1 if x>=0 and zero otherwise. |
angular_basis_EAP_opt (j, l, m, r, theta, phi)
|
|
angular_basis_opt (l, m, q, theta, phi)
|
Angular basis independent of spatial scaling factor us. |
cart2sphere (x, y, z)
|
Return angles for Cartesian 3D coordinates x, y, and z |
create_rt_space_grid (grid_size_r, ...)
|
Generates EAP grid (for potential positivity constraint). |
design_matrix_spatial (bvecs, qvals)
|
Constructs design matrix for DTI weighted least squares or least squares fitting. |
elastic_crossvalidation (b0s_mask, E, M, L, lopt)
|
cross-validation function to find the optimal weight of alpha for sparsity regularization when also Laplacian regularization is used. |
factorial (n[, exact])
|
The factorial of a number or array of numbers. |
factorial2 (n[, exact])
|
Double factorial. |
fmin_l_bfgs_b (func, x0[, fprime, args, ...])
|
Minimize a function func using the L-BFGS-B algorithm. |
generalized_crossvalidation (data, M, LR[, ...])
|
Generalized Cross Validation Function [1]. |
genlaguerre (n, alpha[, monic])
|
Generalized (associated) Laguerre polynomial. |
gradient_table_from_gradient_strength_bvecs (...)
|
A general function for creating diffusion MR gradients. |
l1_crossvalidation (b0s_mask, E, M[, ...])
|
cross-validation function to find the optimal weight of alpha for sparsity regularization |
multi_voxel_fit (single_voxel_fit)
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
optional_package (name[, trip_msg])
|
Return package-like thing and module setup for package name |
part1_reg_matrix_tau (ind_mat, ut)
|
Partial temporal Laplacian regularization matrix following Appendix B in [1]. |
part23_iso_reg_matrix_q (ind_mat, us)
|
Partial spherical spatial Laplacian regularization matrix following the equation below Eq. |
part23_reg_matrix_q (ind_mat, U_mat, T_mat, us)
|
Partial cartesian spatial Laplacian regularization matrix following second line of Eq. |
part23_reg_matrix_tau (ind_mat, ut)
|
Partial temporal Laplacian regularization matrix following Appendix B in [1]. |
part4_iso_reg_matrix_q (ind_mat, us)
|
Partial spherical spatial Laplacian regularization matrix following the equation below Eq. |
part4_reg_matrix_q (ind_mat, U_mat, us)
|
Partial cartesian spatial Laplacian regularization matrix following equation Eq. |
part4_reg_matrix_tau (ind_mat, ut)
|
Partial temporal Laplacian regularization matrix following Appendix B in [1]. |
qtdmri_anisotropic_scaling (data, q, bvecs, tau)
|
Constructs design matrix for fitting an exponential to the diffusion time points. |
qtdmri_eap_matrix (radial_order, time_order, ...)
|
Constructs the design matrix as a product of 3 separated radial, angular and temporal design matrices. |
qtdmri_eap_matrix_ (radial_order, time_order, ...)
|
|
qtdmri_index_matrix (radial_order, time_order)
|
Computes the SHORE basis order indices according to [1]. |
qtdmri_isotropic_eap_matrix (radial_order, ...)
|
Constructs the design matrix as a product of 3 separated radial, angular and temporal design matrices. |
qtdmri_isotropic_eap_matrix_ (radial_order, ...)
|
|
qtdmri_isotropic_index_matrix (radial_order, ...)
|
Computes the SHORE basis order indices according to [1]. |
qtdmri_isotropic_laplacian_reg_matrix (...[, ...])
|
Computes the spherical qt-dMRI Laplacian regularization matrix. |
qtdmri_isotropic_scaling (data, q, tau)
|
Constructs design matrix for fitting an exponential to the diffusion time points. |
qtdmri_isotropic_signal_matrix (radial_order, ...)
|
|
qtdmri_isotropic_signal_matrix_ (...[, ...])
|
|
qtdmri_isotropic_to_mapmri_matrix (...)
|
Generates the matrix that maps the spherical qtdmri coefficients to MAP-MRI coefficients. |
qtdmri_laplacian_reg_matrix (ind_mat, us, ut)
|
Computes the cartesian qt-dMRI Laplacian regularization matrix. |
qtdmri_mapmri_isotropic_normalization (j, l, u0)
|
Normalization factor for Spherical MAP-MRI basis. |
qtdmri_mapmri_normalization (mu)
|
Normalization factor for Cartesian MAP-MRI basis. |
qtdmri_number_of_coefficients (radial_order, ...)
|
Computes the total number of coefficients of the qtdmri basis given a radial and temporal order. |
qtdmri_signal_matrix (radial_order, ...)
|
Constructs the design matrix as a product of 3 separated radial, angular and temporal design matrices. |
qtdmri_signal_matrix_ (radial_order, ...[, ...])
|
Function to generate the qtdmri signal basis. |
qtdmri_temporal_normalization (ut)
|
Normalization factor for the temporal basis |
qtdmri_to_mapmri_matrix (radial_order, ...)
|
Generates the matrix that maps the qtdmri coefficients to MAP-MRI coefficients. |
radial_basis_EAP_opt (j, l, us, r)
|
|
radial_basis_opt (j, l, us, q)
|
Spatial basis dependent on spatial scaling factor us |
real_sh_descoteaux_from_index (m, n, theta, phi)
|
Compute real spherical harmonics as in Descoteaux et al. 2007 [R6e399bd6d05a-1], where the real harmonic \(Y^m_n\) is defined to be:. |
temporal_basis (o, ut, tau)
|
Temporal basis dependent on temporal scaling factor ut |
visualise_gradient_table_G_Delta_rainbow (gtab)
|
This function visualizes a q-tau acquisition scheme as a function of gradient strength and pulse separation (big_delta). |
warn (/, message[, category, stacklevel, source])
|
Issue a warning, or maybe ignore it or raise an exception. |
Module: reconst.qti
Classes and functions for fitting the covariance tensor model of q-space
trajectory imaging (QTI) by Westin et al. as presented in “Q-space trajectory
imaging for multidimensional diffusion MRI of the human brain” NeuroImage vol.
135 (2016): 345-62. https://doi.org/10.1016/j.neuroimage.2016.02.039
QtiFit (params)
|
Methods
|
QtiModel (gtab[, fit_method, cvxpy_solver])
|
Methods
|
ReconstModel (gtab)
|
Abstract class for signal reconstruction models |
Version (version)
|
This class abstracts handling of a project's versions. |
auto_attr (func)
|
Decorator to create OneTimeProperty attributes. |
cvxpy_1x21_to_6x6 (V)
|
Convert 1 x 21 vector into a symmetric 6 x 6 matrix. |
cvxpy_1x6_to_3x3 (V)
|
Convert a 1 x 6 vector into a symmetric 3 x 3 matrix. |
design_matrix (btens)
|
Calculate the design matrix from the b-tensors. |
dtd_covariance (DTD)
|
Calculate covariance of a diffusion tensor distribution (DTD). |
from_21x1_to_6x6 (V)
|
Convert 21 x 1 vectors into symmetric 6 x 6 matrices. |
from_3x3_to_6x1 (T)
|
Convert symmetric 3 x 3 matrices into 6 x 1 vectors. |
from_6x1_to_3x3 (V)
|
Convert 6 x 1 vectors into symmetric 3 x 3 matrices. |
from_6x6_to_21x1 (T)
|
Convert symmetric 6 x 6 matrices into 21 x 1 vectors. |
ndindex (shape)
|
An N-dimensional iterator object to index arrays. |
optional_package (name[, trip_msg])
|
Return package-like thing and module setup for package name |
qti_signal (gtab, D, C[, S0])
|
Generate signals using the covariance tensor signal representation. |
warn (/, message[, category, stacklevel, source])
|
Issue a warning, or maybe ignore it or raise an exception. |
Module: reconst.quick_squash
Detect common dtype across object array
quick_squash
|
Try and make a standard array from an object array |
reduce (function, sequence[, initial])
|
Apply a function of two arguments cumulatively to the items of a sequence, from left to right, so as to reduce the sequence to a single value. |
Module: reconst.recspeed
Optimized routines for creating voxel diffusion models
Module: reconst.rumba
Robust and Unbiased Model-BAsed Spherical Deconvolution (RUMBA-SD)
AxSymShResponse (S0, dwi_response[, bvalue])
|
A simple wrapper for response functions represented using only axially symmetric, even spherical harmonic functions (ie, m == 0 and n even). |
OdfFit (model, data)
|
Methods
|
OdfModel (gtab)
|
An abstract class to be sub-classed by specific odf models |
RumbaFit (model, model_params)
|
Methods
|
RumbaSDModel (gtab[, wm_response, ...])
|
Methods
|
Sphere ([x, y, z, theta, phi, xyz, faces, edges])
|
Points on the unit sphere. |
all_tensor_evecs (e0)
|
Given the principle tensor axis, return the array of all eigenvectors column-wise (or, the rotation matrix that orientates the tensor). |
auto_attr (func)
|
Decorator to create OneTimeProperty attributes. |
bounding_box (vol)
|
Compute the bounding box of nonzero intensity voxels in the volume. |
crop (vol, mins, maxs)
|
Crops the input volume. |
generate_kernel (gtab, sphere, wm_response, ...)
|
Generate deconvolution kernel |
get_bval_indices (bvals, bval[, tol])
|
Get indices where the b-value is bval |
get_sphere ([name])
|
provide triangulated spheres |
gradient_table (bvals[, bvecs, big_delta, ...])
|
A general function for creating diffusion MR gradients. |
lazy_index (index)
|
Produces a lazy index |
mbessel_ratio (n, x)
|
Fast computation of modified Bessel function ratio (first kind). |
normalize_data (data, where_b0[, min_signal, out])
|
Normalizes the data with respect to the mean b0 |
rumba_deconv (data, kernel[, n_iter, ...])
|
Fit fODF and GM/CSF volume fractions for a voxel using RUMBA-SD [1]. |
rumba_deconv_global (data, kernel, mask[, ...])
|
Fit fODF for all voxels simultaneously using RUMBA-SD. |
single_tensor (gtab[, S0, evals, evecs, snr])
|
Simulate diffusion-weighted signals with a single tensor. |
unique_bvals_tolerance (bvals[, tol])
|
Gives the unique b-values of the data, within a tolerance gap |
vec2vec_rotmat (u, v)
|
rotation matrix from 2 unit vectors |
Module: reconst.sfm
The Sparse Fascicle Model.
This is an implementation of the sparse fascicle model described in
[Rokem2015]. The multi b-value version of this model is described in
[Rokem2014].
[Rokem2015]
Ariel Rokem, Jason D. Yeatman, Franco Pestilli, Kendrick
N. Kay, Aviv Mezer, Stefan van der Walt, Brian A. Wandell
(2015). Evaluating the accuracy of diffusion MRI models in white
matter. PLoS ONE 10(4): e0123272. doi:10.1371/journal.pone.0123272
[Rokem2014]
Ariel Rokem, Kimberly L. Chan, Jason D. Yeatman, Franco
Pestilli, Brian A. Wandell (2014). Evaluating the accuracy of diffusion
models at multiple b-values with cross-validation. ISMRM 2014.
Cache ()
|
Cache values based on a key object (such as a sphere or gradient table). |
ExponentialIsotropicFit (model, params)
|
A fit to the ExponentialIsotropicModel object, based on data. |
ExponentialIsotropicModel (gtab)
|
Representing the isotropic signal as a fit to an exponential decay function with b-values |
IsotropicFit (model, params)
|
A fit object for representing the isotropic signal as the mean of the diffusion-weighted signal. |
IsotropicModel (gtab)
|
A base-class for the representation of isotropic signals. |
OrderedDict
|
Dictionary that remembers insertion order |
ReconstFit (model, data)
|
Abstract class which holds the fit result of ReconstModel |
ReconstModel (gtab)
|
Abstract class for signal reconstruction models |
SparseFascicleFit (model, beta, S0, iso)
|
Methods
|
SparseFascicleModel (gtab[, sphere, ...])
|
Methods
|
auto_attr (func)
|
Decorator to create OneTimeProperty attributes. |
determine_num_processes (num_processes)
|
Determine the effective number of processes for parallelization. |
nanmean (a[, axis, dtype, out, keepdims, where])
|
Compute the arithmetic mean along the specified axis, ignoring NaNs. |
optional_package (name[, trip_msg])
|
Return package-like thing and module setup for package name |
sfm_design_matrix (gtab, sphere, response[, mode])
|
Construct the SFM design matrix |
Module: reconst.shm
Tools for using spherical harmonic models to fit diffusion data.
References
[1]
Aganj, I., et al. 2009. ODF Reconstruction in Q-Ball Imaging With Solid
Angle Consideration.
[2]
Descoteaux, M., et al. 2007. Regularized, fast, and robust analytical
Q-ball imaging.
[3]
Tristan-Vega, A., et al. 2010. A new methodology for estimation of fiber
populations in white matter of the brain with Funk-Radon transform.
[4]
Tristan-Vega, A., et al. 2009. Estimation of fiber orientation
probability density functions in high angular resolution diffusion
imaging.
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.
Cache ()
|
Cache values based on a key object (such as a sphere or gradient table). |
CsaOdfModel (gtab, sh_order[, smooth, ...])
|
Implementation of Constant Solid Angle reconstruction method. |
OdfFit (model, data)
|
Methods
|
OdfModel (gtab)
|
An abstract class to be sub-classed by specific odf models |
OpdtModel (gtab, sh_order[, smooth, ...])
|
Implementation of Orientation Probability Density Transform reconstruction method. |
QballBaseModel (gtab, sh_order[, smooth, ...])
|
To be subclassed by Qball type models. |
QballModel (gtab, sh_order[, smooth, ...])
|
Implementation of regularized Qball reconstruction method. |
ResidualBootstrapWrapper (signal_object, B, ...)
|
Returns a residual bootstrap sample of the signal_object when indexed |
SphHarmFit (model, shm_coef, mask)
|
Diffusion data fit to a spherical harmonic model |
SphHarmModel (gtab)
|
To be subclassed by all models that return a SphHarmFit when fit. |
anisotropic_power (sh_coeffs[, norm_factor, ...])
|
Calculate anisotropic power map with a given SH coefficient matrix. |
auto_attr (func)
|
Decorator to create OneTimeProperty attributes. |
bootstrap_data_array (data, H, R[, permute])
|
Applies the Residual Bootstraps to the data given H and R |
bootstrap_data_voxel (data, H, R[, permute])
|
Like bootstrap_data_array but faster when for a single voxel |
calculate_max_order (n_coeffs[, full_basis])
|
Calculate the maximal harmonic order, given that you know the number of parameters that were estimated. |
cart2sphere (x, y, z)
|
Return angles for Cartesian 3D coordinates x, y, and z |
convert_sh_from_legacy (sh_coeffs, sh_basis)
|
Convert SH coefficients in legacy SH basis to SH coefficients of the new SH basis for descoteaux07 [1] or tournier07 [R8020d68d5dcd-2]_[R8020d68d5dcd-3]_ bases. |
convert_sh_to_full_basis (sh_coeffs)
|
Given an array of SH coeffs from a symmetric basis, returns the coefficients for the full SH basis by filling odd order SH coefficients with zeros |
convert_sh_to_legacy (sh_coeffs, sh_basis[, ...])
|
Convert SH coefficients in new SH basis to SH coefficients for the legacy SH basis for descoteaux07 [1] or tournier07 [R2032a14b59d6-2]_[R2032a14b59d6-3]_ bases. |
deprecate_with_version (message[, since, ...])
|
Return decorator function function for deprecation warning / error. |
forward_sdeconv_mat (r_rh, n)
|
Build forward spherical deconvolution matrix |
gen_dirac (m, n, theta, phi[, legacy])
|
Generate Dirac delta function orientated in (theta, phi) on the sphere |
hat (B)
|
Returns the hat matrix for the design matrix B |
lazy_index (index)
|
Produces a lazy index |
lcr_matrix (H)
|
Returns a matrix for computing leveraged, centered residuals from data |
normalize_data (data, where_b0[, min_signal, out])
|
Normalizes the data with respect to the mean b0 |
order_from_ncoef (ncoef[, full_basis])
|
Given a number n of coefficients, calculate back the sh_order |
randint (low[, high, size, dtype])
|
Return random integers from low (inclusive) to high (exclusive). |
real_sh_descoteaux (sh_order, theta, phi[, ...])
|
Compute real spherical harmonics as in Descoteaux et al. 2007 [R191e35d27b5b-1], where the real harmonic \(Y^m_n\) is defined to be:. |
real_sh_descoteaux_from_index (m, n, theta, phi)
|
Compute real spherical harmonics as in Descoteaux et al. 2007 [R700bd48b1688-1], where the real harmonic \(Y^m_n\) is defined to be:. |
real_sh_tournier (sh_order, theta, phi[, ...])
|
Compute real spherical harmonics as initially defined in Tournier 2007 [1] then updated in MRtrix3 [2], where the real harmonic \(Y^m_n\) is defined to be: |
real_sh_tournier_from_index (m, n, theta, phi)
|
Compute real spherical harmonics as initially defined in Tournier 2007 [1] then updated in MRtrix3 [2], where the real harmonic \(Y^m_n\) is defined to be: |
real_sph_harm (m, n, theta, phi)
|
Compute real spherical harmonics. |
real_sym_sh_basis (sh_order, theta, phi)
|
Samples a real symmetric spherical harmonic basis at point on the sphere |
real_sym_sh_mrtrix (sh_order, theta, phi)
|
dipy.reconst.shm.real_sym_sh_mrtrix is deprecated, Please use dipy.reconst.shm.real_sh_tournier instead |
sf_to_sh (sf, sphere[, sh_order, basis_type, ...])
|
Spherical function to spherical harmonics (SH). |
sh_to_rh (r_sh, m, n)
|
Spherical harmonics (SH) to rotational harmonics (RH) |
sh_to_sf (sh, sphere[, sh_order, basis_type, ...])
|
Spherical harmonics (SH) to spherical function (SF). |
sh_to_sf_matrix (sphere[, sh_order, ...])
|
Matrix that transforms Spherical harmonics (SH) to spherical function (SF). |
smooth_pinv (B, L)
|
Regularized pseudo-inverse |
sph_harm_ind_list (sh_order[, full_basis])
|
Returns the degree (m ) and order (n ) of all the symmetric spherical harmonics of degree less then or equal to sh_order . |
spherical_harmonics (m, n, theta, phi[, ...])
|
Compute spherical harmonics. |
warn (/, message[, category, stacklevel, source])
|
Issue a warning, or maybe ignore it or raise an exception. |
Module: reconst.shore
Cache ()
|
Cache values based on a key object (such as a sphere or gradient table). |
ShoreFit (model, shore_coef)
|
- Attributes:
|
ShoreModel (gtab[, radial_order, zeta, ...])
|
Simple Harmonic Oscillator based Reconstruction and Estimation (SHORE) [1] of the diffusion signal. |
Version (version)
|
This class abstracts handling of a project's versions. |
cart2sphere (x, y, z)
|
Return angles for Cartesian 3D coordinates x, y, and z |
create_rspace (gridsize, radius_max)
|
Create the real space table, that contains the points in which |
factorial (x, /)
|
Find x!. |
genlaguerre (n, alpha[, monic])
|
Generalized (associated) Laguerre polynomial. |
l_shore (radial_order)
|
Returns the angular regularisation matrix for SHORE basis |
multi_voxel_fit (single_voxel_fit)
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
n_shore (radial_order)
|
Returns the angular regularisation matrix for SHORE basis |
optional_package (name[, trip_msg])
|
Return package-like thing and module setup for package name |
real_sh_descoteaux_from_index (m, n, theta, phi)
|
Compute real spherical harmonics as in Descoteaux et al. 2007 [Rf0cd5775fa0d-1], where the real harmonic \(Y^m_n\) is defined to be:. |
shore_indices (radial_order, index)
|
Given the basis order and the index, return the shore indices n, l, m for modified Merlet's 3D-SHORE ..math:: :nowrap: begin{equation} textbf{E}(qtextbf{u})=sum_{l=0, even}^{N_{max}} sum_{n=l}^{(N_{max}+l)/2} sum_{m=-l}^l c_{nlm} phi_{nlm}(qtextbf{u}) end{equation} |
shore_matrix (radial_order, zeta, gtab[, tau])
|
Compute the SHORE matrix for modified Merlet's 3D-SHORE [1] |
shore_matrix_odf (radial_order, zeta, ...)
|
Compute the SHORE ODF matrix [1]" |
shore_matrix_pdf (radial_order, zeta, rtab)
|
Compute the SHORE propagator matrix [1]" |
shore_order (n, l, m)
|
Given the indices (n,l,m) of the basis, return the minimum order for those indices and their index for modified Merlet's 3D-SHORE. |
warn (/, message[, category, stacklevel, source])
|
Issue a warning, or maybe ignore it or raise an exception. |
Module: reconst.vec_val_sum
vec_val_vect
|
Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs |
bench
-
dipy.reconst.bench(label='fast', verbose=1, extra_argv=None)
Run benchmarks for module using nose.
- Parameters:
- label{‘fast’, ‘full’, ‘’, attribute identifier}, optional
Identifies the benchmarks to run. This can be a string to pass to
the nosetests executable with the ‘-A’ option, or one of several
special values. Special values are:
‘fast’ - the default - which corresponds to the nosetests -A
option of ‘not slow’.
‘full’ - fast (as above) and slow benchmarks as in the
‘no -A’ option to nosetests - this is the same as ‘’.
None or ‘’ - run all tests.
attribute_identifier - string passed directly to nosetests as ‘-A’.
- verboseint, optional
Verbosity value for benchmark outputs, in the range 1-10. Default is 1.
- extra_argvlist, optional
List with any extra arguments to pass to nosetests.
- Returns:
- successbool
Returns True if running the benchmarks works, False if an error
occurred.
Notes
Benchmarks are like tests, but have names starting with “bench” instead
of “test”, and can be found under the “benchmarks” sub-directory of the
module.
Each NumPy module exposes bench in its namespace to run all benchmarks
for it.
Examples
>>> success = np.lib.bench()
Running benchmarks for numpy.lib
...
using 562341 items:
unique:
0.11
unique1d:
0.11
ratio: 1.0
nUnique: 56230 == 56230
...
OK
test
-
dipy.reconst.test(label='fast', verbose=1, extra_argv=None, doctests=False, coverage=False, raise_warnings=None, timer=False)
Run tests for module using nose.
- Parameters:
- label{‘fast’, ‘full’, ‘’, attribute identifier}, optional
Identifies the tests to run. This can be a string to pass to
the nosetests executable with the ‘-A’ option, or one of several
special values. Special values are:
‘fast’ - the default - which corresponds to the nosetests -A
option of ‘not slow’.
‘full’ - fast (as above) and slow tests as in the
‘no -A’ option to nosetests - this is the same as ‘’.
None or ‘’ - run all tests.
attribute_identifier - string passed directly to nosetests as ‘-A’.
- verboseint, optional
Verbosity value for test outputs, in the range 1-10. Default is 1.
- extra_argvlist, optional
List with any extra arguments to pass to nosetests.
- doctestsbool, optional
If True, run doctests in module. Default is False.
- coveragebool, optional
If True, report coverage of NumPy code. Default is False.
(This requires the
coverage module).
- raise_warningsNone, str or sequence of warnings, optional
This specifies which warnings to configure as ‘raise’ instead
of being shown once during the test execution. Valid strings are:
“develop” : equals (Warning,)
“release” : equals ()
, do not raise on any warnings.
- timerbool or int, optional
Timing of individual tests with nose-timer
(which needs to be
installed). If True, time tests and report on all of them.
If an integer (say N
), report timing results for N
slowest
tests.
- Returns:
- resultobject
Returns the result of running the tests as a
nose.result.TextTestResult
object.
Notes
Each NumPy module exposes test in its namespace to run all tests for it.
For example, to run all tests for numpy.lib:
Examples
>>> result = np.lib.test()
Running unit tests for numpy.lib
...
Ran 976 tests in 3.933s
OK
>>> result.errors
[]
>>> result.knownfail
[]
-
class dipy.reconst.base.ReconstFit(model, data)
Bases: object
Abstract class which holds the fit result of ReconstModel
For example that could be holding FA or GFA etc.
-
__init__(model, data)
-
class dipy.reconst.base.ReconstModel(gtab)
Bases: object
Abstract class for signal reconstruction models
Methods
-
__init__(gtab)
Initialization of the abstract class for signal reconstruction models
- Parameters:
- gtabGradientTable class instance
-
fit(data, mask=None, **kwargs)
bench_bounding_box
-
dipy.reconst.benchmarks.bench_bounding_box.bench_bounding_box()
measure
-
dipy.reconst.benchmarks.bench_bounding_box.measure(code_str, times=1, label=None)
Return elapsed time for executing code in the namespace of the caller.
The supplied code string is compiled with the Python builtin compile
.
The precision of the timing is 10 milli-seconds. If the code will execute
fast on this timescale, it can be executed many times to get reasonable
timing accuracy.
- Parameters:
- code_strstr
The code to be timed.
- timesint, optional
The number of times the code is executed. Default is 1. The code is
only compiled once.
- labelstr, optional
A label to identify code_str with. This is passed into compile
as the second argument (for run-time error messages).
- Returns:
- elapsedfloat
Total elapsed time in seconds for executing code_str times times.
Examples
>>> times = 10
>>> etime = np.testing.measure('for i in range(1000): np.sqrt(i**2)', times=times)
>>> print("Time for a single execution : ", etime / times, "s")
Time for a single execution : 0.005 s
-
class dipy.reconst.benchmarks.bench_csd.ConstrainedSphericalDeconvModel(gtab, response, reg_sphere=None, sh_order=8, lambda_=1, tau=0.1, convergence=50)
Bases: SphHarmModel
Methods
cache_clear ()
|
Clear the cache. |
cache_get (tag, key[, default])
|
Retrieve a value from the cache. |
cache_set (tag, key, value)
|
Store a value in the cache. |
fit (data[, mask])
|
Fit method for every voxel in data |
predict (sh_coeff[, gtab, S0])
|
Compute a signal prediction given spherical harmonic coefficients for the provided GradientTable class instance. |
sampling_matrix (sphere)
|
The matrix needed to sample ODFs from coefficients of the model. |
-
__init__(gtab, response, reg_sphere=None, sh_order=8, lambda_=1, tau=0.1, convergence=50)
Constrained Spherical Deconvolution (CSD) [1].
Spherical deconvolution computes a fiber orientation distribution
(FOD), also called fiber ODF (fODF) [2], as opposed to a diffusion ODF
as the QballModel or the CsaOdfModel. This results in a sharper angular
profile with better angular resolution that is the best object to be
used for later deterministic and probabilistic tractography [3].
A sharp fODF is obtained because a single fiber response function is
injected as a priori knowledge. The response function is often
data-driven and is thus provided as input to the
ConstrainedSphericalDeconvModel. It will be used as deconvolution
kernel, as described in [1].
- Parameters:
- gtabGradientTable
- responsetuple or AxSymShResponse object
A tuple with two elements. The first is the eigen-values as an (3,)
ndarray and the second is the signal value for the response
function without diffusion weighting (i.e. S0). This is to be able
to generate a single fiber synthetic signal. The response function
will be used as deconvolution kernel ([1]).
- reg_sphereSphere (optional)
sphere used to build the regularization B matrix.
Default: ‘symmetric362’.
- sh_orderint (optional)
maximal spherical harmonics order. Default: 8
- lambda_float (optional)
weight given to the constrained-positivity regularization part of
the deconvolution equation (see [1]). Default: 1
- taufloat (optional)
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 mean fODF amplitude (here, 10% by default)
(see [1]). Default: 0.1
- convergenceint
Maximum number of iterations to allow the deconvolution to
converge.
References
[1]
(1,2,3,4,5)
Tournier, J.D., et al. NeuroImage 2007. Robust determination of
the fibre orientation distribution in diffusion MRI:
Non-negativity constrained super-resolved spherical
deconvolution
[2]
Descoteaux, M., et al. IEEE TMI 2009. Deterministic and
Probabilistic Tractography Based on Complex Fibre Orientation
Distributions
[3]
Côté, M-A., et al. Medical Image Analysis 2013. Tractometer:
Towards validation of tractography pipelines
[4]
Tournier, J.D, et al. Imaging Systems and Technology
2012. MRtrix: Diffusion Tractography in Crossing Fiber Regions
-
fit(data, mask=None)
Fit method for every voxel in data
-
predict(sh_coeff, gtab=None, S0=1.0)
Compute a signal prediction given spherical harmonic coefficients
for the provided GradientTable class instance.
- Parameters:
- sh_coeffndarray
The spherical harmonic representation of the FOD from which to make
the signal prediction.
- gtabGradientTable
The gradients for which the signal will be predicted. Uses the
model’s gradient table by default.
- S0ndarray or float
The non diffusion-weighted signal value.
- Returns:
- pred_signdarray
The predicted signal.
-
class dipy.reconst.benchmarks.bench_csd.GradientTable(gradients, big_delta=None, small_delta=None, b0_threshold=50, btens=None)
Bases: object
Diffusion gradient information
- Parameters:
- gradientsarray_like (N, 3)
Diffusion gradients. The direction of each of these vectors corresponds
to the b-vector, and the length corresponds to the b-value.
- b0_thresholdfloat
Gradients with b-value less than or equal to b0_threshold are
considered as b0s i.e. without diffusion weighting.
Notes
The GradientTable object is immutable. Do NOT assign attributes.
If you have your gradient table in a bval & bvec format, we recommend
using the factory function gradient_table
- Attributes:
- gradients(N,3) ndarray
diffusion gradients
- bvals(N,) ndarray
The b-value, or magnitude, of each gradient direction.
- qvals: (N,) ndarray
The q-value for each gradient direction. Needs big and small
delta.
- bvecs(N,3) ndarray
The direction, represented as a unit vector, of each gradient.
- b0s_mask(N,) ndarray
Boolean array indicating which gradients have no diffusion
weighting, ie b-value is close to 0.
- b0_thresholdfloat
Gradients with b-value less than or equal to b0_threshold are
considered to not have diffusion weighting.
- btens(N,3,3) ndarray
The b-tensor of each gradient direction.
Methods
b0s_mask |
|
bvals |
|
bvecs |
|
gradient_strength |
|
qvals |
|
tau |
|
-
__init__(gradients, big_delta=None, small_delta=None, b0_threshold=50, btens=None)
Constructor for GradientTable class
-
b0s_mask()
-
bvals()
-
bvecs()
-
gradient_strength()
-
property info
-
qvals()
-
tau()
bench_csdeconv
-
dipy.reconst.benchmarks.bench_csd.bench_csdeconv(center=(50, 40, 40), width=12)
load_nifti_data
-
dipy.reconst.benchmarks.bench_csd.load_nifti_data(fname, as_ndarray=True)
Load only the data array from a nifti file.
- Parameters:
- fnamestr
Full path to the file.
- as_ndarray: bool, optional
convert nibabel ArrayProxy to a numpy.ndarray.
If you want to save memory and delay this casting, just turn this
option to False (default: True)
- Returns:
- data: np.ndarray or nib.ArrayProxy
num_grad
-
dipy.reconst.benchmarks.bench_csd.num_grad(gtab)
read_stanford_labels
-
dipy.reconst.benchmarks.bench_csd.read_stanford_labels()
Read stanford hardi data and label map.
bench_local_maxima
-
dipy.reconst.benchmarks.bench_peaks.bench_local_maxima()
measure
-
dipy.reconst.benchmarks.bench_peaks.measure(code_str, times=1, label=None)
Return elapsed time for executing code in the namespace of the caller.
The supplied code string is compiled with the Python builtin compile
.
The precision of the timing is 10 milli-seconds. If the code will execute
fast on this timescale, it can be executed many times to get reasonable
timing accuracy.
- Parameters:
- code_strstr
The code to be timed.
- timesint, optional
The number of times the code is executed. Default is 1. The code is
only compiled once.
- labelstr, optional
A label to identify code_str with. This is passed into compile
as the second argument (for run-time error messages).
- Returns:
- elapsedfloat
Total elapsed time in seconds for executing code_str times times.
Examples
>>> times = 10
>>> etime = np.testing.measure('for i in range(1000): np.sqrt(i**2)', times=times)
>>> print("Time for a single execution : ", etime / times, "s")
Time for a single execution : 0.005 s
bench_quick_squash
-
dipy.reconst.benchmarks.bench_squash.bench_quick_squash()
measure
-
dipy.reconst.benchmarks.bench_squash.measure(code_str, times=1, label=None)
Return elapsed time for executing code in the namespace of the caller.
The supplied code string is compiled with the Python builtin compile
.
The precision of the timing is 10 milli-seconds. If the code will execute
fast on this timescale, it can be executed many times to get reasonable
timing accuracy.
- Parameters:
- code_strstr
The code to be timed.
- timesint, optional
The number of times the code is executed. Default is 1. The code is
only compiled once.
- labelstr, optional
A label to identify code_str with. This is passed into compile
as the second argument (for run-time error messages).
- Returns:
- elapsedfloat
Total elapsed time in seconds for executing code_str times times.
Examples
>>> times = 10
>>> etime = np.testing.measure('for i in range(1000): np.sqrt(i**2)', times=times)
>>> print("Time for a single execution : ", etime / times, "s")
Time for a single execution : 0.005 s
ndindex
-
dipy.reconst.benchmarks.bench_squash.ndindex(shape)
An N-dimensional iterator object to index arrays.
Given the shape of an array, an ndindex instance iterates over
the N-dimensional index of the array. At each iteration a tuple
of indices is returned; the last dimension is iterated over first.
- Parameters:
- shapetuple of ints
The dimensions of the array.
Examples
>>> from dipy.core.ndindex import ndindex
>>> shape = (3, 2, 1)
>>> for index in ndindex(shape):
... print(index)
(0, 0, 0)
(0, 1, 0)
(1, 0, 0)
(1, 1, 0)
(2, 0, 0)
(2, 1, 0)
old_squash
-
dipy.reconst.benchmarks.bench_squash.old_squash(arr, mask=None, fill=0)
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:
arr is an array of scalars of type T. Returns an array like
arr.astype(T)
arr is an array of arrays. All items in arr have the same shape
S. Returns an array with shape arr.shape + S.
arr is an array of arrays of different shapes. Returns arr.
Items in arr are not ndarrys or scalars. Returns arr.
- Parameters:
- arrarray, dtype=object
The array to be converted.
- maskarray, dtype=bool, optional
Where arr has Nones.
- fillnumber, optional
Nones are replaced by fill.
- Returns:
- resultarray
Examples
>>> arr = np.empty(3, dtype=object)
>>> arr.fill(2)
>>> old_squash(arr)
array([2, 2, 2])
>>> arr[0] = None
>>> old_squash(arr)
array([0, 2, 2])
>>> arr.fill(np.ones(2))
>>> r = old_squash(arr)
>>> r.shape == (3, 2)
True
>>> r.dtype
dtype('float64')
reduce
-
dipy.reconst.benchmarks.bench_squash.reduce(function, sequence[, initial]) → value
Apply a function of two arguments cumulatively to the items of a sequence,
from left to right, so as to reduce the sequence to a single value.
For example, reduce(lambda x, y: x+y, [1, 2, 3, 4, 5]) calculates
((((1+2)+3)+4)+5). If initial is present, it is placed before the items
of the sequence in the calculation, and serves as a default when the
sequence is empty.
bench_vec_val_vect
-
dipy.reconst.benchmarks.bench_vec_val_sum.bench_vec_val_vect()
measure
-
dipy.reconst.benchmarks.bench_vec_val_sum.measure(code_str, times=1, label=None)
Return elapsed time for executing code in the namespace of the caller.
The supplied code string is compiled with the Python builtin compile
.
The precision of the timing is 10 milli-seconds. If the code will execute
fast on this timescale, it can be executed many times to get reasonable
timing accuracy.
- Parameters:
- code_strstr
The code to be timed.
- timesint, optional
The number of times the code is executed. Default is 1. The code is
only compiled once.
- labelstr, optional
A label to identify code_str with. This is passed into compile
as the second argument (for run-time error messages).
- Returns:
- elapsedfloat
Total elapsed time in seconds for executing code_str times times.
Examples
>>> times = 10
>>> etime = np.testing.measure('for i in range(1000): np.sqrt(i**2)', times=times)
>>> print("Time for a single execution : ", etime / times, "s")
Time for a single execution : 0.005 s
-
class dipy.reconst.cache.Cache
Bases: object
Cache values based on a key object (such as a sphere or gradient table).
Notes
This class is meant to be used as a mix-in:
class MyModel(Model, Cache):
pass
class MyModelFit(Fit):
pass
Inside a method on the fit, typical usage would be:
def odf(sphere):
M = self.model.cache_get('odf_basis_matrix', key=sphere)
if M is None:
M = self._compute_basis_matrix(sphere)
self.model.cache_set('odf_basis_matrix', key=sphere, value=M)
Methods
cache_clear ()
|
Clear the cache. |
cache_get (tag, key[, default])
|
Retrieve a value from the cache. |
cache_set (tag, key, value)
|
Store a value in the cache. |
-
__init__(*args, **kwargs)
-
cache_clear()
Clear the cache.
-
cache_get(tag, key, default=None)
Retrieve a value from the cache.
- Parameters:
- tagstr
Description of the cached value.
- keyobject
Key object used to look up the cached value.
- defaultobject
Value to be returned if no cached entry is found.
- Returns:
- vobject
Value from the cache associated with (tag, key)
. Returns
default if no cached entry is found.
-
cache_set(tag, key, value)
Store a value in the cache.
- Parameters:
- tagstr
Description of the cached value.
- keyobject
Key object used to look up the cached value.
- valueobject
Value stored in the cache for each unique combination
of (tag, key)
.
Examples
>>> def compute_expensive_matrix(parameters):
... # Imagine the following computation is very expensive
... return (p**2 for p in parameters)
>>> parameters = (1, 2, 3)
>>> X1 = compute_expensive_matrix(parameters)
>>> c.cache_set('expensive_matrix', parameters, X1)
>>> X2 = c.cache_get('expensive_matrix', parameters)
auto_attr
-
dipy.reconst.cache.auto_attr(func)
Decorator to create OneTimeProperty attributes.
- Parameters:
- funcmethod
The method that will be called the first time to compute a value.
Afterwards, the method’s name will be a standard attribute holding the
value of this computation.
Examples
>>> class MagicProp(object):
... @auto_attr
... def a(self):
... return 99
...
>>> x = MagicProp()
>>> 'a' in x.__dict__
False
>>> x.a
99
>>> 'a' in x.__dict__
True
coeff_of_determination
-
dipy.reconst.cross_validation.coeff_of_determination(data, model, axis=-1)
- Calculate the coefficient of determination for a model prediction,
relative to data.
- Parameters:
- datandarray
The data
- modelndarray
The predictions of a model for this data. Same shape as the data.
- axis: int, optional
The axis along which different samples are laid out (default: -1).
- Returns:
- CODndarray
The coefficient of determination. This has shape data.shape[:-1]
rac{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).
kfold_xval
-
dipy.reconst.cross_validation.kfold_xval(model, data, folds, *model_args, **model_kwargs)
Perform k-fold cross-validation.
It generate out-of-sample predictions for each measurement.
- Parameters:
- modelModel class instance
The type of the model to use for prediction. The corresponding Fit
object must have a predict function implementd 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.
References
[1]
Rokem, A., Chan, K.L. Yeatman, J.D., Pestilli, F., Mezer, A.,
Wandell, B.A., 2014. Evaluating the accuracy of diffusion models at
multiple b-values with cross-validation. ISMRM 2014.
-
class dipy.reconst.csdeconv.AxSymShResponse(S0, dwi_response, bvalue=None)
Bases: object
A simple wrapper for response functions represented using only axially
symmetric, even spherical harmonic functions (ie, m == 0 and n even).
- Parameters:
- S0float
Signal with no diffusion weighting.
- dwi_responsearray
Response function signal as coefficients to axially symmetric, even
spherical harmonic.
Methods
basis (sphere)
|
A basis that maps the response coefficients onto a sphere. |
on_sphere (sphere)
|
Evaluates the response function on sphere. |
-
__init__(S0, dwi_response, bvalue=None)
-
basis(sphere)
A basis that maps the response coefficients onto a sphere.
-
on_sphere(sphere)
Evaluates the response function on sphere.
-
class dipy.reconst.csdeconv.ConstrainedSDTModel(gtab, ratio, reg_sphere=None, sh_order=8, lambda_=1.0, tau=0.1)
Bases: SphHarmModel
Methods
cache_clear ()
|
Clear the cache. |
cache_get (tag, key[, default])
|
Retrieve a value from the cache. |
cache_set (tag, key, value)
|
Store a value in the cache. |
fit (data[, mask])
|
Fit method for every voxel in data |
sampling_matrix (sphere)
|
The matrix needed to sample ODFs from coefficients of the model. |
-
__init__(gtab, ratio, reg_sphere=None, sh_order=8, lambda_=1.0, tau=0.1)
Spherical Deconvolution Transform (SDT) [1].
The SDT computes a fiber orientation distribution (FOD) as opposed to a
diffusion ODF as the QballModel or the CsaOdfModel. This results in a
sharper angular profile with better angular resolution. The Constrained
SDTModel is similar to the Constrained CSDModel but mathematically it
deconvolves the q-ball ODF as oppposed to the HARDI signal (see [1]
for a comparison and a through discussion).
A sharp fODF is obtained because a single fiber response function is
injected as a priori knowledge. In the SDTModel, this response is a
single fiber q-ball ODF as opposed to a single fiber signal function
for the CSDModel. The response function will be used as deconvolution
kernel.
- Parameters:
- gtabGradientTable
- ratiofloat
ratio of the smallest vs the largest eigenvalue of the single
prolate tensor response function
- reg_sphereSphere
sphere used to build the regularization B matrix
- sh_orderint
maximal spherical harmonics order
- lambda_float
weight given to the constrained-positivity regularization part of
the deconvolution equation
- taufloat
threshold (tau *mean(fODF)) controlling the amplitude below
which the corresponding fODF is assumed to be zero.
References
[1]
(1,2)
Descoteaux, M., et al. IEEE TMI 2009. Deterministic and
Probabilistic Tractography Based on Complex Fibre Orientation
Distributions.
-
fit(data, mask=None)
Fit method for every voxel in data
-
class dipy.reconst.csdeconv.ConstrainedSphericalDeconvModel(gtab, response, reg_sphere=None, sh_order=8, lambda_=1, tau=0.1, convergence=50)
Bases: SphHarmModel
Methods
cache_clear ()
|
Clear the cache. |
cache_get (tag, key[, default])
|
Retrieve a value from the cache. |
cache_set (tag, key, value)
|
Store a value in the cache. |
fit (data[, mask])
|
Fit method for every voxel in data |
predict (sh_coeff[, gtab, S0])
|
Compute a signal prediction given spherical harmonic coefficients for the provided GradientTable class instance. |
sampling_matrix (sphere)
|
The matrix needed to sample ODFs from coefficients of the model. |
-
__init__(gtab, response, reg_sphere=None, sh_order=8, lambda_=1, tau=0.1, convergence=50)
Constrained Spherical Deconvolution (CSD) [1].
Spherical deconvolution computes a fiber orientation distribution
(FOD), also called fiber ODF (fODF) [2], as opposed to a diffusion ODF
as the QballModel or the CsaOdfModel. This results in a sharper angular
profile with better angular resolution that is the best object to be
used for later deterministic and probabilistic tractography [3].
A sharp fODF is obtained because a single fiber response function is
injected as a priori knowledge. The response function is often
data-driven and is thus provided as input to the
ConstrainedSphericalDeconvModel. It will be used as deconvolution
kernel, as described in [1].
- Parameters:
- gtabGradientTable
- responsetuple or AxSymShResponse object
A tuple with two elements. The first is the eigen-values as an (3,)
ndarray and the second is the signal value for the response
function without diffusion weighting (i.e. S0). This is to be able
to generate a single fiber synthetic signal. The response function
will be used as deconvolution kernel ([1]).
- reg_sphereSphere (optional)
sphere used to build the regularization B matrix.
Default: ‘symmetric362’.
- sh_orderint (optional)
maximal spherical harmonics order. Default: 8
- lambda_float (optional)
weight given to the constrained-positivity regularization part of
the deconvolution equation (see [1]). Default: 1
- taufloat (optional)
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 mean fODF amplitude (here, 10% by default)
(see [1]). Default: 0.1
- convergenceint
Maximum number of iterations to allow the deconvolution to
converge.
References
[1]
(1,2,3,4,5)
Tournier, J.D., et al. NeuroImage 2007. Robust determination of
the fibre orientation distribution in diffusion MRI:
Non-negativity constrained super-resolved spherical
deconvolution
[2]
Descoteaux, M., et al. IEEE TMI 2009. Deterministic and
Probabilistic Tractography Based on Complex Fibre Orientation
Distributions
[3]
Côté, M-A., et al. Medical Image Analysis 2013. Tractometer:
Towards validation of tractography pipelines
[4]
Tournier, J.D, et al. Imaging Systems and Technology
2012. MRtrix: Diffusion Tractography in Crossing Fiber Regions
-
fit(data, mask=None)
Fit method for every voxel in data
-
predict(sh_coeff, gtab=None, S0=1.0)
Compute a signal prediction given spherical harmonic coefficients
for the provided GradientTable class instance.
- Parameters:
- sh_coeffndarray
The spherical harmonic representation of the FOD from which to make
the signal prediction.
- gtabGradientTable
The gradients for which the signal will be predicted. Uses the
model’s gradient table by default.
- S0ndarray or float
The non diffusion-weighted signal value.
- Returns:
- pred_signdarray
The predicted signal.
-
class dipy.reconst.csdeconv.SphHarmFit(model, shm_coef, mask)
Bases: OdfFit
Diffusion data fit to a spherical harmonic model
- Attributes:
- shape
shm_coeff
The spherical harmonic coefficients of the odf
Methods
odf (sphere)
|
Samples the odf function on the points of a sphere |
predict ([gtab, S0])
|
Predict the diffusion signal from the model coefficients. |
-
__init__(model, shm_coef, mask)
-
gfa()
-
odf(sphere)
Samples the odf function on the points of a sphere
- Parameters:
- sphereSphere
The points on which to sample the odf.
- Returns:
- valuesndarray
The value of the odf on each point of sphere.
-
predict(gtab=None, S0=1.0)
Predict the diffusion signal from the model coefficients.
- Parameters:
- gtaba GradientTable class instance
The directions and bvalues on which prediction is desired
- S0float array
The mean non-diffusion-weighted signal in each voxel.
-
property shape
-
property shm_coeff
The spherical harmonic coefficients of the odf
Make this a property for now, if there is a use case for modifying
the coefficients we can add a setter or expose the coefficients more
directly
-
class dipy.reconst.csdeconv.SphHarmModel(gtab)
Bases: OdfModel
, Cache
To be subclassed by all models that return a SphHarmFit when fit.
Methods
cache_clear ()
|
Clear the cache. |
cache_get (tag, key[, default])
|
Retrieve a value from the cache. |
cache_set (tag, key, value)
|
Store a value in the cache. |
fit (data)
|
To be implemented by specific odf models |
sampling_matrix (sphere)
|
The matrix needed to sample ODFs from coefficients of the model. |
-
__init__(gtab)
Initialization of the abstract class for signal reconstruction models
- Parameters:
- gtabGradientTable class instance
-
sampling_matrix(sphere)
The matrix needed to sample ODFs from coefficients of the model.
- Parameters:
- sphereSphere
Points used to sample ODF.
- Returns:
- sampling_matrixarray
The size of the matrix will be (N, M) where N is the number of
vertices on sphere and M is the number of coefficients needed by
the model.
-
class dipy.reconst.csdeconv.TensorModel(gtab, fit_method='WLS', return_S0_hat=False, *args, **kwargs)
Bases: ReconstModel
Diffusion Tensor
Methods
fit (data[, mask])
|
Fit method of the DTI model class |
predict (dti_params[, S0])
|
Predict a signal for this TensorModel class instance given parameters. |
-
__init__(gtab, fit_method='WLS', return_S0_hat=False, *args, **kwargs)
A Diffusion Tensor Model [1], [2].
- Parameters:
- gtabGradientTable class instance
- fit_methodstr or callable
str can be one of the following:
- ‘WLS’ for weighted least squares
dti.wls_fit_tensor()
- ‘LS’ or ‘OLS’ for ordinary least squares
dti.ols_fit_tensor()
- ‘NLLS’ for non-linear least-squares
dti.nlls_fit_tensor()
- ‘RT’ or ‘restore’ or ‘RESTORE’ for RESTORE robust tensor
fitting [3]
dti.restore_fit_tensor()
- callable has to have the signature:
fit_method(design_matrix, data, *args, **kwargs)
- return_S0_hatbool
Boolean to return (True) or not (False) the S0 values for the fit.
- args, kwargsarguments and key-word arguments passed to the
fit_method. See dti.wls_fit_tensor, dti.ols_fit_tensor for details
- min_signalfloat
The minimum signal value. Needs to be a strictly positive
number. Default: minimal signal in the data provided to fit.
Notes
In order to increase speed of processing, tensor fitting is done
simultaneously over many voxels. Many fit_methods use the ‘step’
parameter to set the number of voxels that will be fit at once in each
iteration. This is the chunk size as a number of voxels. 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.
E.g., in iter_fit_tensor()
we have a default step value of
1e4
References
[1]
Basser, P.J., Mattiello, J., LeBihan, D., 1994. Estimation of
the effective self-diffusion tensor from the NMR spin echo. J Magn
Reson B 103, 247-254.
[2]
Basser, P., Pierpaoli, C., 1996. Microstructural and
physiological features of tissues elucidated by quantitative
diffusion-tensor MRI. Journal of Magnetic Resonance 111, 209-219.
[3]
Lin-Ching C., Jones D.K., Pierpaoli, C. 2005. RESTORE: Robust
estimation of tensors by outlier rejection. MRM 53: 1088-1095
-
fit(data, mask=None)
Fit method of the DTI model class
- Parameters:
- dataarray
The measured signal from one voxel.
- maskarray
A boolean array used to mark the coordinates in the data that
should be analyzed that has the shape data.shape[:-1]
-
predict(dti_params, S0=1.0)
Predict a signal for this TensorModel class instance given parameters.
- Parameters:
- dti_paramsndarray
The last dimension should have 12 tensor parameters: 3
eigenvalues, followed by the 3 eigenvectors
- S0float or ndarray
The non diffusion-weighted signal in every voxel, or across all
voxels. Default: 1
auto_response
-
dipy.reconst.csdeconv.auto_response(gtab, data, roi_center=None, roi_radius=10, fa_thr=0.7, fa_callable=None, return_number_of_voxels=None)
Automatic estimation of ssst response function using FA.
dipy.reconst.csdeconv.auto_response is deprecated, Please use dipy.reconst.csdeconv.auto_response_ssst instead
- Parameters:
- gtabGradientTable
- 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_radiusint
radius of cubic ROI
- fa_thrfloat
FA threshold
- fa_callablecallable
A callable that defines an operation that compares FA with the fa_thr.
The operator should have two positional arguments
(e.g., fa_operator(FA, fa_thr)) and it should return a bool array.
- return_number_of_voxelsbool
If True, returns the number of voxels used for estimating 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. 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).
auto_response_ssst
-
dipy.reconst.csdeconv.auto_response_ssst(gtab, data, roi_center=None, roi_radii=10, fa_thr=0.7)
- Automatic estimation of single-shell single-tissue (ssst) response
function using FA.
- Parameters:
- gtabGradientTable
- 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).
cart2sphere
-
dipy.reconst.csdeconv.cart2sphere(x, y, z)
Return angles for Cartesian 3D coordinates x, y, and z
See doc for sphere2cart
for angle conventions and derivation
of the formulae.
\(0\le\theta\mathrm{(theta)}\le\pi\) and \(-\pi\le\phi\mathrm{(phi)}\le\pi\)
- Parameters:
- xarray_like
x coordinate in Cartesian space
- yarray_like
y coordinate in Cartesian space
- zarray_like
z coordinate
- Returns:
- rarray
radius
- thetaarray
inclination (polar) angle
- phiarray
azimuth angle
csdeconv
-
dipy.reconst.csdeconv.csdeconv(dwsignal, X, B_reg, tau=0.1, convergence=50, P=None)
Constrained-regularized spherical deconvolution (CSD) [1]
Deconvolves the axially symmetric single fiber response function r_rh in
rotational harmonics coefficients from the diffusion weighted signal in
dwsignal.
- 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.
- Returns:
- fodf_shndarray (
(sh_order + 1)*(sh_order + 2)/2
,) Spherical harmonics coefficients of the constrained-regularized fiber
ODF.
- num_itint
Number of iterations in the constrained-regularization used for
convergence.
Notes
This section describes how the fitting of the SH coefficients is done.
Problem is to minimise per iteration:
\(F(f_n) = ||Xf_n - S||^2 + \lambda^2 ||H_{n-1} f_n||^2\)
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.
Solve by differentiating and setting to zero:
\(\Rightarrow \frac{\delta F}{\delta f_n} = 2X^T(Xf_n - S) + 2 \lambda^2
H_{n-1}^TH_{n-1}f_n=0\)
Or:
\((X^TX + \lambda^2 H_{n-1}^TH_{n-1})f_n = X^Ts\)
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):
\(Q = X^TX + (\lambda H_{n-1}^T) (\lambda H_{n-1}) + \mu I\)
Solve \(Qf_n = X^Ts\) using Cholesky decomposition:
\(Q = LL^T\)
where \(L\) is lower triangular. Then problem can be solved by
back-substitution:
\(L_y = X^Ts\)
\(L^Tf_n = y\)
To speeds things up further, form \(P = X^TX + \mu I\), and update to form
\(Q\) by rankn update with \(H_{n-1}\). The dipy implementation looks like:
form initially \(P = X^T X + \mu I\) and \(\lambda B_{reg}\)
for each voxel: form \(z = X^Ts\)
estimate \(f_0\) by solving \(Pf_0=z\). We use a simplified \(l_{max}=4\)
solution here, but it might not make a big difference.
Then iterate until no change in rows of \(H\) used in \(H_n\)
form \(H_{n}\) given \(f_{n-1}\)
form \(Q = P + (\lambda H_{n-1}^T) (\lambda H_{n-1}\)) (this can
be done by rankn update, but we currently do not use rankn
update).
solve \(Qf_n = z\) using Cholesky decomposition
We’d like to thanks Donald Tournier for his help with describing and
implementing this algorithm.
References
[1]
(1,2)
Tournier, J.D., et al. NeuroImage 2007. Robust determination of the
fibre orientation distribution in diffusion MRI: Non-negativity
constrained super-resolved spherical deconvolution.
deprecate_with_version
-
dipy.reconst.csdeconv.deprecate_with_version(message, since='', until='', version_comparator=<function cmp_pkg_version>, warn_class=<class 'DeprecationWarning'>, error_class=<class 'dipy.utils.deprecator.ExpiredDeprecationError'>)
Return decorator function function for deprecation warning / error.
The decorated function / method will:
Raise the given warning_class warning when the function / method gets
called, up to (and including) version until (if specified);
Raise the given error_class error when the function / method gets
called, when the package version is greater than version until (if
specified).
- Parameters:
- messagestr
Message explaining deprecation, giving possible alternatives.
- sincestr, optional
Released version at which object was first deprecated.
- untilstr, optional
Last released version at which this function will still raise a
deprecation warning. Versions higher than this will raise an
error.
- version_comparatorcallable
Callable accepting string as argument, and return 1 if string
represents a higher version than encoded in the version_comparator, 0
if the version is equal, and -1 if the version is lower. For example,
the version_comparator may compare the input version string to the
current package version string.
- warn_classclass, optional
Class of warning to generate for deprecation.
- error_classclass, optional
Class of error to generate when version_comparator returns 1 for a
given argument of until
.
- Returns:
- deprecatorfunc
Function returning a decorator.
deprecated_params
-
dipy.reconst.csdeconv.deprecated_params(old_name, new_name=None, since='', until='', version_comparator=<function cmp_pkg_version>, arg_in_kwargs=False, warn_class=<class 'dipy.utils.deprecator.ArgsDeprecationWarning'>, error_class=<class 'dipy.utils.deprecator.ExpiredDeprecationError'>, alternative='')
Deprecate a renamed or removed function argument.
The decorator assumes that the argument with the old_name
was removed
from the function signature and the new_name
replaced it at the
same position in the signature. If the old_name
argument is
given when calling the decorated function the decorator will catch it and
issue a deprecation warning and pass it on as new_name
argument.
- Parameters:
- old_namestr or list/tuple thereof
The old name of the argument.
- new_namestr or list/tuple thereof or
None
, optional The new name of the argument. Set this to None to remove the
argument old_name
instead of renaming it.
- sincestr or number or list/tuple thereof, optional
The release at which the old argument became deprecated.
- untilstr or number or list/tuple thereof, optional
Last released version at which this function will still raise a
deprecation warning. Versions higher than this will raise an
error.
- version_comparatorcallable
Callable accepting string as argument, and return 1 if string
represents a higher version than encoded in the version_comparator
,
0 if the version is equal, and -1 if the version is lower. For example,
the version_comparator
may compare the input version string to the
current package version string.
- arg_in_kwargsbool or list/tuple thereof, optional
If the argument is not a named argument (for example it
was meant to be consumed by **kwargs
) set this to
True
. Otherwise the decorator will throw an Exception
if the new_name
cannot be found in the signature of
the decorated function.
Default is False
.
- warn_classwarning, optional
Warning to be issued.
- error_classException, optional
Error to be issued
- alternativestr, optional
An alternative function or class name that the user may use in
place of the deprecated object if new_name
is None. The deprecation
warning will tell the user about this alternative if provided.
- Raises:
- TypeError
If the new argument name cannot be found in the function
signature and arg_in_kwargs was False or if it is used to
deprecate the name of the *args
-, **kwargs
-like arguments.
At runtime such an Error is raised if both the new_name
and old_name were specified when calling the function and
“relax=False”.
Notes
This function is based on the Astropy (major modification).
https://github.com/astropy/astropy. See COPYING file distributed along with
the astropy package for the copyright and license terms.
Examples
The deprecation warnings are not shown in the following examples.
To deprecate a positional or keyword argument::
>>> from dipy.utils.deprecator import deprecated_params
>>> @deprecated_params(‘sig’, ‘sigma’, ‘0.3’)
… def test(sigma):
… return sigma
>>> test(2)
2
>>> test(sigma=2)
2
>>> test(sig=2) # doctest: +SKIP
2
It is also possible to replace multiple arguments. The old_name
,
new_name
and since
have to be tuple or list and contain the
same number of entries::
>>> @deprecated_params([‘a’, ‘b’], [‘alpha’, ‘beta’],
… [‘0.2’, 0.4])
… def test(alpha, beta):
… return alpha, beta
>>> test(a=2, b=3) # doctest: +SKIP
(2, 3)
estimate_response
-
dipy.reconst.csdeconv.estimate_response(gtab, evals, S0)
Estimate single fiber response function
- Parameters:
- gtabGradientTable
- evalsndarray
- S0float
non diffusion weighted
- Returns:
- Sestimated signal
fa_trace_to_lambdas
-
dipy.reconst.csdeconv.fa_trace_to_lambdas(fa=0.08, trace=0.0021)
forward_sdeconv_mat
-
dipy.reconst.csdeconv.forward_sdeconv_mat(r_rh, n)
Build forward spherical deconvolution matrix
- Parameters:
- r_rhndarray
Rotational harmonics coefficients for the single fiber response
function. Each element rh[i]
is associated with spherical harmonics
of degree 2*i
.
- nndarray
The order of spherical harmonic function associated with each row of
the deconvolution matrix. Only even orders are allowed
- Returns:
- Rndarray (N, N)
Deconvolution matrix with shape (N, N)
forward_sdt_deconv_mat
-
dipy.reconst.csdeconv.forward_sdt_deconv_mat(ratio, n, r2_term=False)
Build forward sharpening deconvolution transform (SDT) matrix
- Parameters:
- ratiofloat
ratio = \(\frac{\lambda_2}{\lambda_1}\) of the single fiber response
function
- nndarray (N,)
The degree of spherical harmonic function associated with each row of
the deconvolution matrix. Only even degrees 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 [1].
- Returns:
- Rndarray (N, N)
SDT deconvolution matrix
- Pndarray (N, N)
Funk-Radon Transform (FRT) matrix
References
[1]
Descoteaux, M. PhD Thesis. INRIA Sophia-Antipolis. 2008.
fractional_anisotropy
-
dipy.reconst.csdeconv.fractional_anisotropy(evals, axis=-1)
Return Fractional anisotropy (FA) of a diffusion tensor.
- Parameters:
- evalsarray-like
Eigenvalues of a diffusion tensor.
- axisint
Axis of evals which contains 3 eigenvalues.
- Returns:
- faarray
Calculated FA. Range is 0 <= FA <= 1.
Notes
FA is calculated using the following equation:
\[FA = \sqrt{\frac{1}{2}\frac{(\lambda_1-\lambda_2)^2+(\lambda_1-
\lambda_3)^2+(\lambda_2-\lambda_3)^2}{\lambda_1^2+
\lambda_2^2+\lambda_3^2}}\]
get_sphere
-
dipy.reconst.csdeconv.get_sphere(name='symmetric362')
provide triangulated spheres
- Parameters:
- namestr
which sphere - one of:
* ‘symmetric362’
* ‘symmetric642’
* ‘symmetric724’
* ‘repulsion724’
* ‘repulsion100’
* ‘repulsion200’
- Returns:
- spherea dipy.core.sphere.Sphere class instance
Examples
>>> import numpy as np
>>> from dipy.data import get_sphere
>>> sphere = get_sphere('symmetric362')
>>> verts, faces = sphere.vertices, sphere.faces
>>> verts.shape == (362, 3)
True
>>> faces.shape == (720, 3)
True
>>> verts, faces = get_sphere('not a sphere name')
Traceback (most recent call last):
...
DataError: No sphere called "not a sphere name"
lazy_index
-
dipy.reconst.csdeconv.lazy_index(index)
Produces a lazy index
Returns a slice that can be used for indexing an array, if no slice can be
made index is returned as is.
lpn
-
dipy.reconst.csdeconv.lpn(n, z)
Legendre function of the first kind.
Compute sequence of Legendre functions of the first kind (polynomials),
Pn(z) and derivatives for all degrees from 0 to n (inclusive).
See also special.legendre for polynomial class.
References
mask_for_response_ssst
-
dipy.reconst.csdeconv.mask_for_response_ssst(gtab, data, roi_center=None, roi_radii=10, fa_thr=0.7)
- Computation of mask for single-shell single-tissue (ssst) response
function using FA.
- Parameters:
- gtabGradientTable
- 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 [1]).
References
[1]
Tournier, J.D., et al. NeuroImage 2004. Direct estimation of the
fiber orientation density function from diffusion-weighted MRI
data using spherical deconvolution
multi_voxel_fit
-
dipy.reconst.csdeconv.multi_voxel_fit(single_voxel_fit)
Method decorator to turn a single voxel model fit
definition into a multi voxel model fit definition
ndindex
-
dipy.reconst.csdeconv.ndindex(shape)
An N-dimensional iterator object to index arrays.
Given the shape of an array, an ndindex instance iterates over
the N-dimensional index of the array. At each iteration a tuple
of indices is returned; the last dimension is iterated over first.
- Parameters:
- shapetuple of ints
The dimensions of the array.
Examples
>>> from dipy.core.ndindex import ndindex
>>> shape = (3, 2, 1)
>>> for index in ndindex(shape):
... print(index)
(0, 0, 0)
(0, 1, 0)
(1, 0, 0)
(1, 1, 0)
(2, 0, 0)
(2, 1, 0)
odf_deconv
-
dipy.reconst.csdeconv.odf_deconv(odf_sh, R, B_reg, lambda_=1.0, tau=0.1, r2_term=False)
ODF constrained-regularized spherical deconvolution using
the Sharpening Deconvolution Transform (SDT) [1], [2].
- Parameters:
- odf_shndarray (
(sh_order + 1)*(sh_order + 2)/2
,) ndarray of SH coefficients for the ODF spherical function to be
deconvolved
- Rndarray (
(sh_order + 1)(sh_order + 2)/2
, (sh_order + 1)(sh_order + 2)/2
)
SDT matrix in SH basis
- B_regndarray (
(sh_order + 1)(sh_order + 2)/2
, (sh_order + 1)(sh_order + 2)/2
)
SH basis matrix used for deconvolution
- lambda_float
lambda parameter in minimization equation (default 1.0)
- taufloat
threshold (tau *max(fODF)) controlling the amplitude below
which the corresponding fODF is assumed to be zero.
- r2_termbool
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 [1]) and
the true normalized ODF definition differ from a \(r^2\) term in the ODF
integral. The original Sharpening Deconvolution Transform (SDT)
technique [2] is expecting Tuch’s ODF without the \(r^2\) (see [3] 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 (
(sh_order + 1)(sh_order + 2)/2
,) Spherical harmonics coefficients of the constrained-regularized fiber
ODF
- num_itint
Number of iterations in the constrained-regularization used for
convergence
References
[1]
(1,2,3)
Tuch, D. MRM 2004. Q-Ball Imaging.
[2]
(1,2,3)
Descoteaux, M., et al. IEEE TMI 2009. Deterministic and
Probabilistic Tractography Based on Complex Fibre Orientation
Distributions
[3]
Descoteaux, M, PhD thesis, INRIA Sophia-Antipolis, 2008.
odf_sh_to_sharp
-
dipy.reconst.csdeconv.odf_sh_to_sharp(odfs_sh, sphere, basis=None, ratio=0.2, sh_order=8, lambda_=1.0, tau=0.1, r2_term=False)
Sharpen odfs using the sharpening deconvolution transform [2]
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.
- Parameters:
- odfs_shndarray (
(sh_order + 1)*(sh_order + 2)/2
, ) array of odfs expressed as spherical harmonics coefficients
- sphereSphere
sphere used to build the regularization matrix
- basis{None, ‘tournier07’, ‘descoteaux07’}
different spherical harmonic basis:
None
for the default DIPY basis,
tournier07
for the Tournier 2007 [4] basis, and
descoteaux07
for the Descoteaux 2007 [3] basis
(None
defaults to descoteaux07
).
- ratiofloat,
ratio of the smallest vs the largest eigenvalue of the single prolate
tensor response function (\(\frac{\lambda_2}{\lambda_1}\))
- sh_orderint
maximal SH order of the SH representation
- lambda_float
lambda parameter (see odfdeconv) (default 1.0)
- taufloat
tau parameter in the L matrix construction (see odfdeconv)
(default 0.1)
- r2_termbool
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 [1]) and
the true normalized ODF definition differ from a \(r^2\) term in the ODF
integral. The original Sharpening Deconvolution Transform (SDT)
technique [2] is expecting Tuch’s ODF without the \(r^2\) (see [3] 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
References
[1]
Tuch, D. MRM 2004. Q-Ball Imaging.
[2]
(1,2,3)
Descoteaux, M., et al. IEEE TMI 2009. Deterministic and
Probabilistic Tractography Based on Complex Fibre Orientation
Distributions
[3]
(1,2)
Descoteaux, M., Angelino, E., Fitzgibbons, S. and Deriche, R.
Regularized, Fast, and Robust Analytical Q-ball Imaging.
Magn. Reson. Med. 2007;58:497-510.
[4]
Tournier J.D., Calamante F. and Connelly A. Robust determination
of the fibre orientation distribution in diffusion MRI:
Non-negativity constrained super-resolved spherical deconvolution.
NeuroImage. 2007;35(4):1459-1472.
peaks_from_model
-
dipy.reconst.csdeconv.peaks_from_model(model, data, sphere, relative_peak_threshold, min_separation_angle, mask=None, return_odf=False, return_sh=True, gfa_thr=0, normalize_peaks=False, sh_order=8, sh_basis_type=None, npeaks=5, B=None, invB=None, parallel=False, num_processes=None)
Fit the model to data and computes peaks and metrics
- Parameters:
- modela model instance
model will be used to fit the data.
- datandarray
Diffusion data.
- sphereSphere
The Sphere providing discrete directions for evaluation.
- relative_peak_thresholdfloat
Only return peaks greater than relative_peak_threshold * m
where m
is the largest peak.
- min_separation_anglefloat in [0, 90] The minimum distance between
directions. If two peaks are too close only the larger of the two is
returned.
- maskarray, optional
If mask is provided, voxels that are False in mask are skipped and
no peaks are returned.
- return_odfbool
If True, the odfs are returned.
- return_shbool
If True, the odf as spherical harmonics coefficients is returned
- gfa_thrfloat
Voxels with gfa less than gfa_thr are skipped, no peaks are returned.
- normalize_peaksbool
If true, all peak values are calculated relative to max(odf).
- sh_orderint, optional
Maximum SH order in the SH fit. For sh_order, there will be
(sh_order + 1) * (sh_order + 2) / 2
SH coefficients (default 8).
- sh_basis_type{None, ‘tournier07’, ‘descoteaux07’}
None
for the default DIPY basis,
tournier07
for the Tournier 2007 [2] basis, and
descoteaux07
for the Descoteaux 2007 [1] basis
(None
defaults to descoteaux07
).
- npeaksint
Maximum number of peaks found (default 5 peaks).
- Bndarray, optional
Matrix that transforms spherical harmonics to spherical function
sf = np.dot(sh, B)
.
- invBndarray, optional
Inverse of B.
- parallel: bool
If True, use multiprocessing to compute peaks and metric
(default False). Temporary files are saved in the default temporary
directory of the system. It can be changed using import tempfile
and tempfile.tempdir = '/path/to/tempdir'
.
- num_processes: int, 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.
- Returns:
- pamPeaksAndMetrics
An object with gfa
, peak_directions
, peak_values
,
peak_indices
, odf
, shm_coeffs
as attributes
References
[1]
Descoteaux, M., Angelino, E., Fitzgibbons, S. and Deriche, R.
Regularized, Fast, and Robust Analytical Q-ball Imaging.
Magn. Reson. Med. 2007;58:497-510.
[2]
Tournier J.D., Calamante F. and Connelly A. Robust determination
of the fibre orientation distribution in diffusion MRI:
Non-negativity constrained super-resolved spherical deconvolution.
NeuroImage. 2007;35(4):1459-1472.
quad
-
dipy.reconst.csdeconv.quad(func, a, b, args=(), full_output=0, epsabs=1.49e-08, epsrel=1.49e-08, limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50, limlst=50, complex_func=False)
Compute a definite integral.
Integrate func from a to b (possibly infinite interval) using a
technique from the Fortran library QUADPACK.
- Parameters:
- func{function, scipy.LowLevelCallable}
A Python function or method to integrate. If func takes many
arguments, it is integrated along the axis corresponding to the
first argument.
If the user desires improved integration performance, then f may
be a scipy.LowLevelCallable with one of the signatures:
double func(double x)
double func(double x, void *user_data)
double func(int n, double *xx)
double func(int n, double *xx, void *user_data)
The user_data
is the data contained in the scipy.LowLevelCallable.
In the call forms with xx
, n
is the length of the xx
array which contains xx[0] == x
and the rest of the items are
numbers contained in the args
argument of quad.
In addition, certain ctypes call signatures are supported for
backward compatibility, but those should not be used in new code.
- afloat
Lower limit of integration (use -numpy.inf for -infinity).
- bfloat
Upper limit of integration (use numpy.inf for +infinity).
- argstuple, optional
Extra arguments to pass to func.
- full_outputint, optional
Non-zero to return a dictionary of integration information.
If non-zero, warning messages are also suppressed and the
message is appended to the output tuple.
- complex_funcbool, optional
Indicate if the function’s (func) return type is real
(complex_func=False
: default) or complex (complex_func=True
).
In both cases, the function’s argument is real.
If full_output is also non-zero, the infodict, message, and
explain for the real and complex components are returned in
a dictionary with keys “real output” and “imag output”.
- Returns:
- yfloat
The integral of func from a to b.
- abserrfloat
An estimate of the absolute error in the result.
- infodictdict
A dictionary containing additional information.
- message
A convergence message.
- explain
Appended only with ‘cos’ or ‘sin’ weighting and infinite
integration limits, it contains an explanation of the codes in
infodict[‘ierlst’]
- Other Parameters:
- epsabsfloat or int, optional
Absolute error tolerance. Default is 1.49e-8. quad tries to obtain
an accuracy of abs(i-result) <= max(epsabs, epsrel*abs(i))
where i
= integral of func from a to b, and result
is the
numerical approximation. See epsrel below.
- epsrelfloat or int, optional
Relative error tolerance. Default is 1.49e-8.
If epsabs <= 0
, epsrel must be greater than both 5e-29
and 50 * (machine epsilon)
. See epsabs above.
- limitfloat or int, optional
An upper bound on the number of subintervals used in the adaptive
algorithm.
- points(sequence of floats,ints), optional
A sequence of break points in the bounded integration interval
where local difficulties of the integrand may occur (e.g.,
singularities, discontinuities). The sequence does not have
to be sorted. Note that this option cannot be used in conjunction
with weight
.
- weightfloat or int, optional
String indicating weighting function. Full explanation for this
and the remaining arguments can be found below.
- wvaroptional
Variables for use with weighting functions.
- woptsoptional
Optional input for reusing Chebyshev moments.
- maxp1float or int, optional
An upper bound on the number of Chebyshev moments.
- limlstint, optional
Upper bound on the number of cycles (>=3) for use with a sinusoidal
weighting and an infinite end-point.
See also
dblquad
double integral
tplquad
triple integral
nquad
n-dimensional integrals (uses quad recursively)
fixed_quad
fixed-order Gaussian quadrature
quadrature
adaptive Gaussian quadrature
odeint
ODE integrator
ode
ODE integrator
simpson
integrator for sampled data
romb
integrator for sampled data
scipy.special
for coefficients and roots of orthogonal polynomials
Notes
Extra information for quad() inputs and outputs
If full_output is non-zero, then the third output argument
(infodict) is a dictionary with entries as tabulated below. For
infinite limits, the range is transformed to (0,1) and the
optional outputs are given with respect to this transformed range.
Let M be the input argument limit and let K be infodict[‘last’].
The entries are:
- ‘neval’
The number of function evaluations.
- ‘last’
The number, K, of subintervals produced in the subdivision process.
- ‘alist’
A rank-1 array of length M, the first K elements of which are the
left end points of the subintervals in the partition of the
integration range.
- ‘blist’
A rank-1 array of length M, the first K elements of which are the
right end points of the subintervals.
- ‘rlist’
A rank-1 array of length M, the first K elements of which are the
integral approximations on the subintervals.
- ‘elist’
A rank-1 array of length M, the first K elements of which are the
moduli of the absolute error estimates on the subintervals.
- ‘iord’
A rank-1 integer array of length M, the first L elements of
which are pointers to the error estimates over the subintervals
with L=K
if K<=M/2+2
or L=M+1-K
otherwise. Let I be the
sequence infodict['iord']
and let E be the sequence
infodict['elist']
. Then E[I[1]], ..., E[I[L]]
forms a
decreasing sequence.
If the input argument points is provided (i.e., it is not None),
the following additional outputs are placed in the output
dictionary. Assume the points sequence is of length P.
- ‘pts’
A rank-1 array of length P+2 containing the integration limits
and the break points of the intervals in ascending order.
This is an array giving the subintervals over which integration
will occur.
- ‘level’
A rank-1 integer array of length M (=limit), containing the
subdivision levels of the subintervals, i.e., if (aa,bb) is a
subinterval of (pts[1], pts[2])
where pts[0]
and pts[2]
are adjacent elements of infodict['pts']
, then (aa,bb) has level l
if |bb-aa| = |pts[2]-pts[1]| * 2**(-l)
.
- ‘ndin’
A rank-1 integer array of length P+2. After the first integration
over the intervals (pts[1], pts[2]), the error estimates over some
of the intervals may have been increased artificially in order to
put their subdivision forward. This array has ones in slots
corresponding to the subintervals for which this happens.
Weighting the integrand
The input variables, weight and wvar, are used to weight the
integrand by a select list of functions. Different integration
methods are used to compute the integral with these weighting
functions, and these do not support specifying break points. The
possible values of weight and the corresponding weighting functions are.
weight
|
Weight function used |
wvar
|
‘cos’ |
cos(w*x) |
wvar = w |
‘sin’ |
sin(w*x) |
wvar = w |
‘alg’ |
g(x) = ((x-a)**alpha)*((b-x)**beta) |
wvar = (alpha, beta) |
‘alg-loga’ |
g(x)*log(x-a) |
wvar = (alpha, beta) |
‘alg-logb’ |
g(x)*log(b-x) |
wvar = (alpha, beta) |
‘alg-log’ |
g(x)*log(x-a)*log(b-x) |
wvar = (alpha, beta) |
‘cauchy’ |
1/(x-c) |
wvar = c |
wvar holds the parameter w, (alpha, beta), or c depending on the weight
selected. In these expressions, a and b are the integration limits.
For the ‘cos’ and ‘sin’ weighting, additional inputs and outputs are
available.
For finite integration limits, the integration is performed using a
Clenshaw-Curtis method which uses Chebyshev moments. For repeated
calculations, these moments are saved in the output dictionary:
- ‘momcom’
The maximum level of Chebyshev moments that have been computed,
i.e., if M_c
is infodict['momcom']
then the moments have been
computed for intervals of length |b-a| * 2**(-l)
,
l=0,1,...,M_c
.
- ‘nnlog’
A rank-1 integer array of length M(=limit), containing the
subdivision levels of the subintervals, i.e., an element of this
array is equal to l if the corresponding subinterval is
|b-a|* 2**(-l)
.
- ‘chebmo’
A rank-2 array of shape (25, maxp1) containing the computed
Chebyshev moments. These can be passed on to an integration
over the same interval by passing this array as the second
element of the sequence wopts and passing infodict[‘momcom’] as
the first element.
If one of the integration limits is infinite, then a Fourier integral is
computed (assuming w neq 0). If full_output is 1 and a numerical error
is encountered, besides the error message attached to the output tuple,
a dictionary is also appended to the output tuple which translates the
error codes in the array info['ierlst']
to English messages. The
output information dictionary contains the following entries instead of
‘last’, ‘alist’, ‘blist’, ‘rlist’, and ‘elist’:
- ‘lst’
The number of subintervals needed for the integration (call it K_f
).
- ‘rslst’
A rank-1 array of length M_f=limlst, whose first K_f
elements
contain the integral contribution over the interval
(a+(k-1)c, a+kc)
where c = (2*floor(|w|) + 1) * pi / |w|
and k=1,2,...,K_f
.
- ‘erlst’
A rank-1 array of length M_f
containing the error estimate
corresponding to the interval in the same position in
infodict['rslist']
.
- ‘ierlst’
A rank-1 integer array of length M_f
containing an error flag
corresponding to the interval in the same position in
infodict['rslist']
. See the explanation dictionary (last entry
in the output tuple) for the meaning of the codes.
Details of QUADPACK level routines
quad calls routines from the FORTRAN library QUADPACK. This section
provides details on the conditions for each routine to be called and a
short description of each routine. The routine called depends on
weight, points and the integration limits a and b.
QUADPACK routine |
weight |
points |
infinite bounds |
qagse |
None |
No |
No |
qagie |
None |
No |
Yes |
qagpe |
None |
Yes |
No |
qawoe |
‘sin’, ‘cos’ |
No |
No |
qawfe |
‘sin’, ‘cos’ |
No |
either a or b |
qawse |
‘alg*’ |
No |
No |
qawce |
‘cauchy’ |
No |
No |
The following provides a short desciption from [1] for each
routine.
- qagse
is an integrator based on globally adaptive interval
subdivision in connection with extrapolation, which will
eliminate the effects of integrand singularities of
several types.
- qagie
handles integration over infinite intervals. The infinite range is
mapped onto a finite interval and subsequently the same strategy as
in QAGS
is applied.
- qagpe
serves the same purposes as QAGS, but also allows the
user to provide explicit information about the location
and type of trouble-spots i.e. the abscissae of internal
singularities, discontinuities and other difficulties of
the integrand function.
- qawoe
is an integrator for the evaluation of
\(\int^b_a \cos(\omega x)f(x)dx\) or
\(\int^b_a \sin(\omega x)f(x)dx\)
over a finite interval [a,b], where \(\omega\) and \(f\)
are specified by the user. The rule evaluation component is based
on the modified Clenshaw-Curtis technique
An adaptive subdivision scheme is used in connection
with an extrapolation procedure, which is a modification
of that in QAGS
and allows the algorithm to deal with
singularities in \(f(x)\).
- qawfe
calculates the Fourier transform
\(\int^\infty_a \cos(\omega x)f(x)dx\) or
\(\int^\infty_a \sin(\omega x)f(x)dx\)
for user-provided \(\omega\) and \(f\). The procedure of
QAWO
is applied on successive finite intervals, and convergence
acceleration by means of the \(\varepsilon\)-algorithm is applied
to the series of integral approximations.
- qawse
approximate \(\int^b_a w(x)f(x)dx\), with \(a < b\) where
\(w(x) = (x-a)^{\alpha}(b-x)^{\beta}v(x)\) with
\(\alpha,\beta > -1\), where \(v(x)\) may be one of the
following functions: \(1\), \(\log(x-a)\), \(\log(b-x)\),
\(\log(x-a)\log(b-x)\).
The user specifies \(\alpha\), \(\beta\) and the type of the
function \(v\). A globally adaptive subdivision strategy is
applied, with modified Clenshaw-Curtis integration on those
subintervals which contain a or b.
- qawce
compute \(\int^b_a f(x) / (x-c)dx\) where the integral must be
interpreted as a Cauchy principal value integral, for user specified
\(c\) and \(f\). The strategy is globally adaptive. Modified
Clenshaw-Curtis integration is used on those intervals containing the
point \(x = c\).
Integration of Complex Function of a Real Variable
A complex valued function, \(f\), of a real variable can be written as
\(f = g + ih\). Similarly, the integral of \(f\) can be
written as
\[\int_a^b f(x) dx = \int_a^b g(x) dx + i\int_a^b h(x) dx\]
assuming that the integrals of \(g\) and \(h\) exist
over the inteval \([a,b]\) [2]. Therefore, quad
integrates
complex-valued functions by integrating the real and imaginary components
separately.
References
[1]
Piessens, Robert; de Doncker-Kapenga, Elise;
Überhuber, Christoph W.; Kahaner, David (1983).
QUADPACK: A subroutine package for automatic integration.
Springer-Verlag.
ISBN 978-3-540-12553-2.
[2]
McCullough, Thomas; Phillips, Keith (1973).
Foundations of Analysis in the Complex Plane.
Holt Rinehart Winston.
ISBN 0-03-086370-8
Examples
Calculate \(\int^4_0 x^2 dx\) and compare with an analytic result
>>> from scipy import integrate
>>> import numpy as np
>>> x2 = lambda x: x**2
>>> integrate.quad(x2, 0, 4)
(21.333333333333332, 2.3684757858670003e-13)
>>> print(4**3 / 3.) # analytical result
21.3333333333
Calculate \(\int^\infty_0 e^{-x} dx\)
>>> invexp = lambda x: np.exp(-x)
>>> integrate.quad(invexp, 0, np.inf)
(1.0, 5.842605999138044e-11)
Calculate \(\int^1_0 a x \,dx\) for \(a = 1, 3\)
>>> f = lambda x, a: a*x
>>> y, err = integrate.quad(f, 0, 1, args=(1,))
>>> y
0.5
>>> y, err = integrate.quad(f, 0, 1, args=(3,))
>>> y
1.5
Calculate \(\int^1_0 x^2 + y^2 dx\) with ctypes, holding
y parameter as 1:
testlib.c =>
double func(int n, double args[n]){
return args[0]*args[0] + args[1]*args[1];}
compile to library testlib.*
from scipy import integrate
import ctypes
lib = ctypes.CDLL('/home/.../testlib.*') #use absolute path
lib.func.restype = ctypes.c_double
lib.func.argtypes = (ctypes.c_int,ctypes.c_double)
integrate.quad(lib.func,0,1,(1))
#(1.3333333333333333, 1.4802973661668752e-14)
print((1.0**3/3.0 + 1.0) - (0.0**3/3.0 + 0.0)) #Analytic result
# 1.3333333333333333
Be aware that pulse shapes and other sharp features as compared to the
size of the integration interval may not be integrated correctly using
this method. A simplified example of this limitation is integrating a
y-axis reflected step function with many zero values within the integrals
bounds.
>>> y = lambda x: 1 if x<=0 else 0
>>> integrate.quad(y, -1, 1)
(1.0, 1.1102230246251565e-14)
>>> integrate.quad(y, -1, 100)
(1.0000000002199108, 1.0189464580163188e-08)
>>> integrate.quad(y, -1, 10000)
(0.0, 0.0)
real_sh_descoteaux
-
dipy.reconst.csdeconv.real_sh_descoteaux(sh_order, theta, phi, full_basis=False, legacy=True)
Compute real spherical harmonics as in Descoteaux et al. 2007 [Ra6d8f6cd2652-1],
where the real harmonic \(Y^m_n\) is defined to be:
Imag(\(Y^m_n\)) * sqrt(2) if m > 0
\(Y^0_n\) if m = 0
Real(\(Y^m_n\)) * sqrt(2) if m < 0
This may take scalar or array arguments. The inputs will be broadcast
against each other.
- Parameters:
- sh_orderint
The maximum degree or the spherical harmonic basis.
- thetafloat [0, pi]
The polar (colatitudinal) coordinate.
- phifloat [0, 2*pi]
The azimuthal (longitudinal) coordinate.
- full_basis: bool, optional
If true, returns a basis including odd order SH functions as well as
even order SH functions. Otherwise returns only even order SH
functions.
- legacy: bool, optional
If true, uses DIPY’s legacy descoteaux07 implementation (where |m|
for m < 0). Else, implements the basis as defined in Descoteaux et al.
2007.
- Returns:
- real_shreal float
The real harmonic \(Y^m_n\) sampled at theta
and phi
.
- marray
The degree of the harmonics.
- narray
The order of the harmonics.
real_sh_descoteaux_from_index
-
dipy.reconst.csdeconv.real_sh_descoteaux_from_index(m, n, theta, phi, legacy=True)
Compute real spherical harmonics as in Descoteaux et al. 2007 [R2763a8caef3c-1],
where the real harmonic \(Y^m_n\) is defined to be:
Imag(\(Y^m_n\)) * sqrt(2) if m > 0
\(Y^0_n\) if m = 0
Real(\(Y^m_n\)) * sqrt(2) if m < 0
This may take scalar or array arguments. The inputs will be broadcast
against each other.
- Parameters:
- mint
|m| <= n
The degree of the harmonic.
- nint
>= 0
The order of the harmonic.
- 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^m_n\) sampled at theta
and phi
.
recursive_response
-
dipy.reconst.csdeconv.recursive_response(gtab, data, mask=None, sh_order=8, peak_thr=0.01, init_fa=0.08, init_trace=0.0021, iter=8, convergence=0.001, parallel=False, num_processes=None, sphere=<dipy.core.sphere.HemiSphere object>)
Recursive calibration of response function using peak threshold
- Parameters:
- gtabGradientTable
- datandarray
diffusion data
- maskndarray, optional
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_orderint, optional
maximal spherical harmonics order. Default: 8
- 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 [1]. Default: 0.01
- init_fafloat, optional
FA of the initial ‘fat’ response function (tensor). Default: 0.08
- init_tracefloat, optional
trace of the initial ‘fat’ response function (tensor). Default: 0.0021
- iterint, optional
maximum number of iterations for calibration. Default: 8.
- convergencefloat, optional
convergence criterion, maximum relative change of SH
coefficients. Default: 0.001.
- parallelbool, optional
Whether to use parallelization in peak-finding during the calibration
procedure. Default: True
- 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. Default: default_sphere.
- 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 [1].
References
[1]
Tax, C.M.W., et al. NeuroImage 2014. Recursive calibration of
the fiber response function for spherical deconvolution of
diffusion MRI data.
response_from_mask
-
dipy.reconst.csdeconv.response_from_mask(gtab, data, mask)
- Computation of single-shell single-tissue (ssst) response
function from a given mask.
dipy.reconst.csdeconv.response_from_mask is deprecated, Please use dipy.reconst.csdeconv.response_from_mask_ssst instead
- Parameters:
- gtabGradientTable
- 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[R64bc28ca561d-1]_). 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.
References
[1]
Tournier, J.D., et al. NeuroImage 2004. Direct estimation of the
fiber orientation density function from diffusion-weighted MRI
data using spherical deconvolution
response_from_mask_ssst
-
dipy.reconst.csdeconv.response_from_mask_ssst(gtab, data, mask)
- Computation of single-shell single-tissue (ssst) response
function from a given mask.
- Parameters:
- gtabGradientTable
- 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[Rb280cc58c30e-1]_). 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.
References
[1]
Tournier, J.D., et al. NeuroImage 2004. Direct estimation of the
fiber orientation density function from diffusion-weighted MRI
data using spherical deconvolution
sh_to_rh
-
dipy.reconst.csdeconv.sh_to_rh(r_sh, m, n)
Spherical harmonics (SH) to rotational harmonics (RH)
Calculate the rotational harmonic decomposition up to
harmonic order n
, degree m
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.
- Parameters:
- r_shndarray (N,)
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.
- mndarray (N,)
The degree of the spherical harmonic function associated with each
coefficient.
- nndarray (N,)
The order 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
See also
shm.real_sh_descoteaux_from_index
, shm.real_sh_descoteaux
References
[1]
Tournier, J.D., et al. NeuroImage 2007. Robust determination of the
fibre orientation distribution in diffusion MRI: Non-negativity
constrained super-resolved spherical deconvolution
single_tensor
-
dipy.reconst.csdeconv.single_tensor(gtab, S0=1, evals=None, evecs=None, snr=None)
Simulate diffusion-weighted signals with a single tensor.
- Parameters:
- gtabGradientTable
Table with information of b-values and gradient directions g.
Note that if gtab has a btens attribute, simulations will be performed
according to the given b-tensor B information.
- S0double, optional
Strength of signal in the presence of no diffusion gradient (also
called the b=0
value).
- evals(3,) ndarray, optional
Eigenvalues of the diffusion tensor. By default, values typical for
prolate white matter are used.
- evecs(3, 3) ndarray, optional
Eigenvectors of the tensor. You can also think of this as a rotation
matrix that transforms the direction of the tensor. The eigenvectors
need to be column wise.
- snrfloat, optional
Signal to noise ratio, assuming Rician noise. None implies no noise.
- Returns:
- S(N,) ndarray
- Simulated signal:
S(b, g) = S_0 e^(-b g^T R D R.T g)
, if gtab.tens=None
S(B) = S_0 e^(-B:D)
, if gtab.tens information is given
References
[1]
M. Descoteaux, “High Angular Resolution Diffusion MRI: from Local
Estimation to Segmentation and Tractography”, PhD thesis,
University of Nice-Sophia Antipolis, p. 42, 2008.
[2]
E. Stejskal and J. Tanner, “Spin diffusion measurements: spin echos
in the presence of a time-dependent field gradient”, Journal of
Chemical Physics, nr. 42, pp. 288–292, 1965.
sph_harm_ind_list
-
dipy.reconst.csdeconv.sph_harm_ind_list(sh_order, full_basis=False)
Returns the degree (m
) and order (n
) of all the symmetric spherical
harmonics of degree less then or equal to sh_order
. The results,
m_list
and n_list
are kx1 arrays, where k depends on sh_order
.
They can be passed to real_sh_descoteaux_from_index()
and
:func:real_sh_tournier_from_index
.
- Parameters:
- sh_orderint
even int > 0, max order to return
- full_basis: bool, optional
True for SH basis with even and odd order terms
- Returns:
- m_listarray
degrees of even spherical harmonics
- n_listarray
orders of even spherical harmonics
See also
shm.real_sh_descoteaux_from_index
, shm.real_sh_tournier_from_index
vec2vec_rotmat
-
dipy.reconst.csdeconv.vec2vec_rotmat(u, v)
rotation matrix from 2 unit vectors
u, v being unit 3d vectors return a 3x3 rotation matrix R than aligns u to
v.
In general there are many rotations that will map u to v. If S is any
rotation using v as an axis then R.S will also map u to v since (S.R)u =
S(Ru) = Sv = v. The rotation R returned by vec2vec_rotmat leaves fixed the
perpendicular to the plane spanned by u and v.
The transpose of R will align v to u.
- Parameters:
- uarray, shape(3,)
- varray, shape(3,)
- Returns:
- Rarray, shape(3,3)
Examples
>>> import numpy as np
>>> from dipy.core.geometry import vec2vec_rotmat
>>> u=np.array([1,0,0])
>>> v=np.array([0,1,0])
>>> R=vec2vec_rotmat(u,v)
>>> np.dot(R,u)
array([ 0., 1., 0.])
>>> np.dot(R.T,v)
array([ 1., 0., 0.])
-
class dipy.reconst.dki.DiffusionKurtosisFit(model, model_params)
Bases: TensorFit
Class for fitting the Diffusion Kurtosis Model
- Attributes:
- S0_hat
directions
For tracking - return the primary direction in each voxel
evals
Returns the eigenvalues of the tensor as an array
evecs
Returns the eigenvectors of the tensor as an array, columnwise
kfa
Returns the kurtosis tensor (KFA) .
kt
Returns the 15 independent elements of the kurtosis tensor as an array
quadratic_form
Calculates the 3x3 diffusion tensor for each voxel
- shape
Methods
ad ()
|
Axial diffusivity (AD) calculated from cached eigenvalues. |
adc (sphere)
|
Calculate the apparent diffusion coefficient (ADC) in each direction on |
ak ([min_kurtosis, max_kurtosis, analytical])
|
Axial Kurtosis (AK) of a diffusion kurtosis tensor [1]. |
akc (sphere)
|
Calculates 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])
|
Computes the maximum value of a single voxel kurtosis tensor |
linearity ()
|
- Returns:
|
md ()
|
Mean diffusivity (MD) calculated from cached eigenvalues. |
mk ([min_kurtosis, max_kurtosis, analytical])
|
Computes mean Kurtosis (MK) from the kurtosis tensor. |
mkt ([min_kurtosis, max_kurtosis])
|
Computes mean of the kurtosis tensor (MKT) [1]. |
mode ()
|
Tensor mode calculated from cached eigenvalues. |
odf (sphere)
|
The diffusion orientation distribution function (dODF). |
planarity ()
|
- Returns:
|
predict (gtab[, S0])
|
Given a DKI model fit, predict the signal on the vertices of a gradient table |
rd ()
|
Radial diffusivity (RD) calculated from cached eigenvalues. |
rk ([min_kurtosis, max_kurtosis, analytical])
|
Radial Kurtosis (RK) of a diffusion kurtosis tensor [1]. |
sphericity ()
|
- Returns:
|
trace ()
|
Trace of the tensor calculated from cached eigenvalues. |
-
__init__(model, model_params)
Initialize a DiffusionKurtosisFit class instance.
Since DKI is an extension of DTI, class instance is defined as subclass
of the TensorFit from dti.py
- Parameters:
- modelDiffusionKurtosisModel Class instance
Class instance containing the Diffusion Kurtosis Model for the fit
- model_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
-
ak(min_kurtosis=-0.42857142857142855, max_kurtosis=10, analytical=True)
Axial Kurtosis (AK) of a diffusion kurtosis tensor [1].
- 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 [2])
- 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. Default = 10
- 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). Default is set to True.
- Returns:
- akarray
Calculated AK.
Notes
AK is defined as the directional kurtosis parallel to the fiber’s main
direction e1 [1], [2]. You can compute AK using to approaches:
AK is calculated from rotated diffusion kurtosis tensor [2], i.e.:
\[AK = \hat{W}_{1111}
\frac{(\lambda_{1}+\lambda_{2}+\lambda_{3})^2}{(9 \lambda_{1}^2)}\]
AK can be sampled from the principal axis of the diffusion tensor:
\[AK = K(\mathbf{\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.
References
[1]
(1,2,3)
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of
non-Gaussian water diffusion by kurtosis analysis. NMR in
Biomedicine 23(7): 698-710
[2]
(1,2,3)
Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011.
Estimation of tensors and tensor-derived measures in diffusional
kurtosis imaging. Magn Reson Med. 65(3), 823-836
[3]
Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging:
Robust estimation from DW-MRI using homogeneous polynomials.
Proceedings of the 8th {IEEE} International Symposium on
Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265.
doi: 10.1109/ISBI.2011.5872402
-
akc(sphere)
Calculates the apparent kurtosis coefficient (AKC) in each
direction on the sphere for each voxel in the data
- Parameters:
- sphereSphere class instance
- Returns:
- akcndarray
The estimates of the apparent kurtosis coefficient in every
direction on the input sphere
Notes
For each sphere direction with coordinates \((n_{1}, n_{2}, n_{3})\), the
calculation of AKC is done using formula:
\[AKC(n)=\frac{MD^{2}}{ADC(n)^{2}}\sum_{i=1}^{3}\sum_{j=1}^{3}
\sum_{k=1}^{3}\sum_{l=1}^{3}n_{i}n_{j}n_{k}n_{l}W_{ijkl}\]
where \(W_{ijkl}\) are the elements of the kurtosis tensor, MD the mean
diffusivity and ADC the apparent diffusion coefficent computed as:
\[ADC(n)=\sum_{i=1}^{3}\sum_{j=1}^{3}n_{i}n_{j}D_{ij}\]
where \(D_{ij}\) are the elements of the diffusion tensor.
-
property kfa
Returns the kurtosis tensor (KFA) [1].
Notes
The KFA is defined as [1]:
\[KFA \equiv
\frac{||\mathbf{W} - MKT \mathbf{I}^{(4)}||_F}{||\mathbf{W}||_F}\]
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 [1].
References
[1]
(1,2,3)
Glenn, G. R., Helpern, J. A., Tabesh, A., and Jensen, J. H.
(2015). Quantitative assessment of diffusional kurtosis
anisotropy. NMR in Biomedicine 28, 448–459. doi:10.1002/nbm.3271
-
kmax(sphere='repulsion100', gtol=1e-05, mask=None)
Computes the maximum value of a single voxel kurtosis tensor
- Parameters:
- sphereSphere class instance, optional
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
- Returns:
- max_valuefloat
kurtosis tensor maximum value
-
property kt
Returns the 15 independent elements of the kurtosis tensor as an array
-
mk(min_kurtosis=-0.42857142857142855, max_kurtosis=10, analytical=True)
Computes mean Kurtosis (MK) from the kurtosis tensor.
- 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 [4])
- 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. Default = 10
- analyticalbool (optional)
If True, MK is calculated using its analytical solution, otherwise
an exact numerical estimator is used (see Notes). Default is set to
True.
- 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[R1a4c5980fd18-1]_:
\[MK \equiv \frac{1}{4\pi} \int d\Omega_\mathbf{n} K(\mathbf{n})\]
This integral can be numerically solved by averaging directional
kurtosis values sampled for directions of a spherical t-design [2].
Alternatively, MK can be solved from the analytical solution derived by
Tabesh et al. [3]. This solution is given by:
\[\begin{split}MK=F_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{1111}+
F_1(\lambda_2,\lambda_1,\lambda_3)\hat{W}_{2222}+
F_1(\lambda_3,\lambda_2,\lambda_1)\hat{W}_{3333}+ \\
F_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233}+
F_2(\lambda_2,\lambda_1,\lambda_3)\hat{W}_{1133}+
F_2(\lambda_3,\lambda_2,\lambda_1)\hat{W}_{1122}\end{split}\]
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
\[ \begin{align}\begin{aligned}\begin{split}F_1(\lambda_1,\lambda_2,\lambda_3)=
\frac{(\lambda_1+\lambda_2+\lambda_3)^2}
{18(\lambda_1-\lambda_2)(\lambda_1-\lambda_3)}
[\frac{\sqrt{\lambda_2\lambda_3}}{\lambda_1}
R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\
\frac{3\lambda_1^2-\lambda_1\lambda_2-\lambda_2\lambda_3-
\lambda_1\lambda_3}
{3\lambda_1 \sqrt{\lambda_2 \lambda_3}}
R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-1 ]\end{split}\\\begin{split}F_2(\lambda_1,\lambda_2,\lambda_3)=
\frac{(\lambda_1+\lambda_2+\lambda_3)^2}
{3(\lambda_2-\lambda_3)^2}
[\frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}}
R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\
\frac{2\lambda_1-\lambda_2-\lambda_3}{3\sqrt{\lambda_2 \lambda_3}}
R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-2]\end{split}\end{aligned}\end{align} \]
where \(R_f\) and \(R_d\) are the Carlson’s elliptic integrals.
References
[1]
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of
non-Gaussian water diffusion by kurtosis analysis. NMR in
Biomedicine 23(7): 698-710
[2]
Hardin, R.H., Sloane, N.J.A., 1996. McLaren’s Improved Snub Cube
and Other New Spherical Designs in Three Dimensions. Discrete
and Computational Geometry 15, 429-441.
[3]
Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011.
Estimation of tensors and tensor-derived measures in diffusional
kurtosis imaging. Magn Reson Med. 65(3), 823-836
[4]
Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging:
Robust estimation from DW-MRI using homogeneous polynomials.
Proceedings of the 8th {IEEE} International Symposium on
Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265.
doi: 10.1109/ISBI.2011.5872402
-
mkt(min_kurtosis=-0.42857142857142855, max_kurtosis=10)
Computes mean of the kurtosis tensor (MKT) [1].
- 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 [2])
- 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. Default = 10
- Returns:
- mktarray
Calculated mean kurtosis tensor.
Notes
The MKT is defined as [1]:
\[MKT \equiv \frac{1}{4\pi} \int d
\Omega_{\mathnbf{n}} n_i n_j n_k n_l W_{ijkl}\]
which can be directly computed from the trace of the kurtosis tensor:
\[\]
MKT = frac{1}{5} Tr(mathbf{W}) = frac{1}{5}
(W_{1111} + W_{2222} + W_{3333} + 2W_{1122} + 2W_{1133} + 2W_{2233})
References
[1]
(1,2,3)
Hansen, B., Lund, T. E., Sangill, R., and Jespersen, S. N. 2013.
Experimentally and computationally fast method for estimation
of a mean kurtosis. Magnetic Resonance in Medicine69, 1754–1760.
388. doi:10.1002/mrm.24743
[2]
Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging:
Robust estimation from DW-MRI using homogeneous polynomials.
Proceedings of the 8th {IEEE} International Symposium on
Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265.
doi: 10.1109/ISBI.2011.5872402
-
predict(gtab, S0=1.0)
Given a DKI model fit, predict the signal on the vertices of a
gradient table
- Parameters:
- gtaba GradientTable class instance
The gradient table for this prediction
- S0float or ndarray (optional)
The non diffusion-weighted signal in every voxel, or across all
voxels. Default: 1
Notes
The predicted signal is given by:
\[S(n,b)=S_{0}e^{-bD(n)+\frac{1}{6}b^{2}D(n)^{2}K(n)}\]
\(\mathbf{D(n)}\) and \(\mathbf{K(n)}\) can be computed from the DT and KT
using the following equations:
\[D(n)=\sum_{i=1}^{3}\sum_{j=1}^{3}n_{i}n_{j}D_{ij}\]
and
\[K(n)=\frac{MD^{2}}{D(n)^{2}}\sum_{i=1}^{3}\sum_{j=1}^{3}
\sum_{k=1}^{3}\sum_{l=1}^{3}n_{i}n_{j}n_{k}n_{l}W_{ijkl}\]
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.
-
rk(min_kurtosis=-0.42857142857142855, max_kurtosis=10, analytical=True)
Radial Kurtosis (RK) of a diffusion kurtosis tensor [1].
- Parameters:
- 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
[3])
- 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. Default = 10
- 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 [1], [2]:
\[\]
- RK equiv frac{1}{2pi} int dOmega _mathbf{theta}
K(mathbf{theta}) delta (mathbf{theta}cdot mathbf{e}_1)
This equation can be numerically computed by averaging apparent
directional kurtosis samples for directions perpendicular to e1.
Otherwise, RK can be calculated from its analytical solution [2]:
\[K_{\bot} = G_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2222} +
G_1(\lambda_1,\lambda_3,\lambda_2)\hat{W}_{3333} +
G_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233}\]
where:
\[G_1(\lambda_1,\lambda_2,\lambda_3)=
\frac{(\lambda_1+\lambda_2+\lambda_3)^2}{18\lambda_2(\lambda_2-
\lambda_3)} \left (2\lambda_2 +
\frac{\lambda_3^2-3\lambda_2\lambda_3}{\sqrt{\lambda_2\lambda_3}}
\right)\]
and
\[ G_2(\lambda_1,\lambda_2,\lambda_3)=
\frac{(\lambda_1+\lambda_2+\lambda_3)^2}{(\lambda_2-\lambda_3)^2}
\left ( \frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}}-
2\right )\]
References
[1]
(1,2,3)
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of
non-Gaussian water diffusion by kurtosis analysis. NMR in
Biomedicine 23(7): 698-710
[2]
(1,2)
Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011.
Estimation of tensors and tensor-derived measures in diffusional
kurtosis imaging. Magn Reson Med. 65(3), 823-836
[3]
Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging:
Robust estimation from DW-MRI using homogeneous polynomials.
Proceedings of the 8th {IEEE} International Symposium on
Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265.
doi: 10.1109/ISBI.2011.5872402
-
class dipy.reconst.dki.DiffusionKurtosisModel(gtab, fit_method='WLS', *args, **kwargs)
Bases: ReconstModel
Class for the Diffusion Kurtosis Model
Methods
fit (data[, mask])
|
Fit method of the DKI model class |
predict (dki_params[, S0])
|
Predict a signal for this DKI model class instance given parameters. |
-
__init__(gtab, fit_method='WLS', *args, **kwargs)
Diffusion Kurtosis Tensor Model [1]
- Parameters:
- gtabGradientTable class instance
- fit_methodstr or callable
str can be one of the following:
‘OLS’ or ‘ULLS’ for ordinary least squares
- ‘WLS’ or ‘UWLLS’ for weighted ordinary least squares
dki.wls_fit_dki
- callable has to have the signature:
fit_method(design_matrix, data, *args, **kwargs)
- args, kwargsarguments and key-word arguments passed to the
fit_method. See dki.ols_fit_dki, dki.wls_fit_dki for details
References
[1]
Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011.
Estimation of tensors and tensor-derived measures in diffusional
kurtosis imaging. Magn Reson Med. 65(3), 823-836
-
fit(data, mask=None)
Fit method of the DKI model class
- Parameters:
- dataarray
The measured signal from one voxel.
- maskarray
A boolean array used to mark the coordinates in the data that
should be analyzed that has the shape data.shape[-1]
-
predict(dki_params, S0=1.0)
Predict a signal for this DKI model class instance given parameters.
- 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
2. Three lines of the eigenvector matrix each containing the
first, second and third coordinates of the eigenvector
3. Fifteen elements of the kurtosis tensor
- S0float or ndarray (optional)
The non diffusion-weighted signal in every voxel, or across all
voxels. Default: 1
-
class dipy.reconst.dki.ReconstModel(gtab)
Bases: object
Abstract class for signal reconstruction models
Methods
-
__init__(gtab)
Initialization of the abstract class for signal reconstruction models
- Parameters:
- gtabGradientTable class instance
-
fit(data, mask=None, **kwargs)
-
class dipy.reconst.dki.TensorFit(model, model_params, model_S0=None)
Bases: object
- Attributes:
- S0_hat
directions
For tracking - return the primary direction in each voxel
evals
Returns the eigenvalues of the tensor as an array
evecs
Returns the eigenvectors of the tensor as an array, columnwise
quadratic_form
Calculates the 3x3 diffusion tensor for each voxel
- shape
Methods
ad ()
|
Axial diffusivity (AD) calculated from cached eigenvalues. |
adc (sphere)
|
Calculate the apparent diffusion coefficient (ADC) in each direction on |
color_fa ()
|
Color fractional anisotropy of diffusion tensor |
fa ()
|
Fractional anisotropy (FA) calculated from cached eigenvalues. |
ga ()
|
Geodesic anisotropy (GA) calculated from cached eigenvalues. |
linearity ()
|
- Returns:
|
md ()
|
Mean diffusivity (MD) calculated from cached eigenvalues. |
mode ()
|
Tensor mode calculated from cached eigenvalues. |
odf (sphere)
|
The diffusion orientation distribution function (dODF). |
planarity ()
|
- Returns:
|
predict (gtab[, S0, step])
|
Given a model fit, predict the signal on the vertices of a sphere |
rd ()
|
Radial diffusivity (RD) calculated from cached eigenvalues. |
sphericity ()
|
- Returns:
|
trace ()
|
Trace of the tensor calculated from cached eigenvalues. |
-
__init__(model, model_params, model_S0=None)
Initialize a TensorFit class instance.
-
property S0_hat
-
ad()
Axial diffusivity (AD) calculated from cached eigenvalues.
- Returns:
- adarray (V, 1)
Calculated AD.
Notes
RD is calculated with the following equation:
\[AD = \lambda_1\]
-
adc(sphere)
Calculate the apparent diffusion coefficient (ADC) in each direction on
the sphere for each voxel in the data
- Parameters:
- sphereSphere class instance
- Returns:
- adcndarray
The estimates of the apparent diffusion coefficient in every
direction on the input sphere
ec{b} Q
ec{b}^T
Where Q is the quadratic form of the tensor.
-
color_fa()
Color fractional anisotropy of diffusion tensor
-
property directions
For tracking - return the primary direction in each voxel
-
property evals
Returns the eigenvalues of the tensor as an array
-
property evecs
Returns the eigenvectors of the tensor as an array, columnwise
-
fa()
Fractional anisotropy (FA) calculated from cached eigenvalues.
-
ga()
Geodesic anisotropy (GA) calculated from cached eigenvalues.
-
linearity()
- Returns:
- linearityarray
Calculated linearity of the diffusion tensor [1].
Notes
Linearity is calculated with the following equation:
\[Linearity =
\frac{\lambda_1-\lambda_2}{\lambda_1+\lambda_2+\lambda_3}\]
References
[1]
Westin C.-F., Peled S., Gubjartsson H., Kikinis R., Jolesz
F., “Geometrical diffusion measures for MRI from tensor basis
analysis” in Proc. 5th Annual ISMRM, 1997.
-
lower_triangular(b0=None)
-
md()
Mean diffusivity (MD) calculated from cached eigenvalues.
- Returns:
- mdarray (V, 1)
Calculated MD.
Notes
MD is calculated with the following equation:
\[MD = \frac{\lambda_1+\lambda_2+\lambda_3}{3}\]
-
mode()
Tensor mode calculated from cached eigenvalues.
-
odf(sphere)
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 [1]. To re-derive it from
scratch, follow steps in [2], Section 7.9 Equation
7.24 but with an \(r^2\) term in the integral.
References
[1]
Aganj, I., Lenglet, C., Sapiro, G., Yacoub, E., Ugurbil,
K., & Harel, N. (2010). Reconstruction of the orientation
distribution function in single- and multiple-shell q-ball imaging
within constant solid angle. Magnetic Resonance in Medicine, 64(2),
554-566. doi:DOI: 10.1002/mrm.22365
-
planarity()
- Returns:
- sphericityarray
Calculated sphericity of the diffusion tensor [1].
Notes
Sphericity is calculated with the following equation:
\[Sphericity =
\frac{2 (\lambda_2 - \lambda_3)}{\lambda_1+\lambda_2+\lambda_3}\]
References
[1]
Westin C.-F., Peled S., Gubjartsson H., Kikinis R., Jolesz
F., “Geometrical diffusion measures for MRI from tensor basis
analysis” in Proc. 5th Annual ISMRM, 1997.
-
predict(gtab, S0=None, step=None)
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
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
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( heta, b) = S_0 * e^{-b ADC}\]
Where:
.. math
:math:` heta` 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
-
property quadratic_form
Calculates the 3x3 diffusion tensor for each voxel
-
rd()
Radial diffusivity (RD) calculated from cached eigenvalues.
- Returns:
- rdarray (V, 1)
Calculated RD.
Notes
RD is calculated with the following equation:
\[RD = \frac{\lambda_2 + \lambda_3}{2}\]
-
property shape
-
sphericity()
- Returns:
- sphericityarray
Calculated sphericity of the diffusion tensor [1].
Notes
Sphericity is calculated with the following equation:
\[Sphericity = \frac{3 \lambda_3}{\lambda_1+\lambda_2+\lambda_3}\]
References
[1]
Westin C.-F., Peled S., Gubjartsson H., Kikinis R., Jolesz
F., “Geometrical diffusion measures for MRI from tensor basis
analysis” in Proc. 5th Annual ISMRM, 1997.
-
trace()
Trace of the tensor calculated from cached eigenvalues.
- Returns:
- tracearray (V, 1)
Calculated trace.
Notes
The trace is calculated with the following equation:
\[trace = \lambda_1 + \lambda_2 + \lambda_3\]
Wcons
-
dipy.reconst.dki.Wcons(k_elements)
Construct the full 4D kurtosis tensors from its 15 independent
elements
- Parameters:
- k_elements(15,)
elements of the kurtosis tensor in the following order:
- .. math::
- begin{matrix} ( & W_{xxxx} & W_{yyyy} & W_{zzzz} & W_{xxxy} & W_{xxxz}
& … \
& W_{xyyy} & W_{yyyz} & W_{xzzz} & W_{yzzz} & W_{xxyy}
& … \
& W_{xxzz} & W_{yyzz} & W_{xxyz} & W_{xyyz} & W_{xyzz}
& & )end{matrix}
- Returns:
- Warray(3, 3, 3, 3)
Full 4D kurtosis tensor
Wrotate
-
dipy.reconst.dki.Wrotate(kt, Basis)
Rotate a kurtosis tensor from the standard Cartesian coordinate system
to another coordinate system basis
- Parameters:
- kt(15,)
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
KT elements are assumed to be ordered as follows:
\[\]
- begin{matrix} ( & W_{xxxx} & W_{yyyy} & W_{zzzz} & W_{xxxy} & W_{xxxz}
& … \
& W_{xyyy} & W_{yyyz} & W_{xzzz} & W_{yzzz} & W_{xxyy}
& … \
& W_{xxzz} & W_{yyzz} & W_{xxyz} & W_{xyyz} & W_{xyzz}
& & )end{matrix}
References
[1] Hui ES, Cheung MM, Qi L, Wu EX, 2008. Towards better MR
characterization of neural tissues using directional diffusion kurtosis
analysis. Neuroimage 42(1): 122-34
Wrotate_element
-
dipy.reconst.dki.Wrotate_element(kt, indi, indj, indk, indl, B)
Computes the the specified index element of a kurtosis tensor rotated
to the coordinate system basis B.
- Parameters:
- ktndarray (x, y, z, 15) or (n, 15)
Array containing the 15 independent elements of the kurtosis tensor
- indiint
Rotated kurtosis tensor element index i (0 for x, 1 for y, 2 for z)
- indjint
Rotated kurtosis tensor element index j (0 for x, 1 for y, 2 for z)
- indkint
Rotated kurtosis tensor element index k (0 for x, 1 for y, 2 for z)
- indl: int
Rotated kurtosis tensor element index l (0 for x, 1 for y, 2 for z)
- B: array (x, y, z, 3, 3) or (n, 15)
Vectors of the basis column-wise oriented
- Returns:
- Wrefloat
rotated kurtosis tensor element of index ind_i, ind_j, ind_k, ind_l
Notes
It is assumed that initial kurtosis tensor elementes are defined on the
Cartesian coordinate system.
References
[1] Hui ES, Cheung MM, Qi L, Wu EX, 2008. Towards better MR
characterization of neural tissues using directional diffusion kurtosis
analysis. Neuroimage 42(1): 122-34
apparent_kurtosis_coef
-
dipy.reconst.dki.apparent_kurtosis_coef(dki_params, sphere, min_diffusivity=0, min_kurtosis=-0.42857142857142855)
Calculates the apparent kurtosis coefficient (AKC) in each direction
of a sphere [1].
- 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. Default = 0
- min_kurtosisfloat (optional)
Because high-amplitude negative values of kurtosis are not physicaly
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 [2])
- 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 [1]:
\[AKC(n)=\frac{MD^{2}}{ADC(n)^{2}}\sum_{i=1}^{3}\sum_{j=1}^{3}
\sum_{k=1}^{3}\sum_{l=1}^{3}n_{i}n_{j}n_{k}n_{l}W_{ijkl}\]
where \(W_{ijkl}\) are the elements of the kurtosis tensor, MD the mean
diffusivity and ADC the apparent diffusion coefficent computed as:
\[ADC(n)=\sum_{i=1}^{3}\sum_{j=1}^{3}n_{i}n_{j}D_{ij}\]
where \(D_{ij}\) are the elements of the diffusion tensor.
References
[1]
(1,2,3)
Neto Henriques R, Correia MM, Nunes RG, Ferreira HA (2015).
Exploring the 3D geometry of the diffusion kurtosis tensor -
Impact on the development of robust tractography procedures and
novel biomarkers, NeuroImage 111: 85-99
[2]
Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging:
Robust estimation from DW-MRI using homogeneous polynomials.
Proceedings of the 8th {IEEE} International Symposium on
Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265.
doi: 10.1109/ISBI.2011.5872402
axial_kurtosis
-
dipy.reconst.dki.axial_kurtosis(dki_params, min_kurtosis=-0.42857142857142855, max_kurtosis=10, analytical=True)
Computes axial Kurtosis (AK) from the kurtosis tensor [1], [2].
- 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 [3])
- 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. Default = 10
- 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).
Default is set to True.
- Returns:
- akarray
Calculated AK.
Notes
AK is defined as the directional kurtosis parallel to the fiber’s main
direction e1 [1], [2]. You can compute AK using to approaches:
AK is calculated from rotated diffusion kurtosis tensor [2], i.e.:
\[AK = \hat{W}_{1111}
\frac{(\lambda_{1}+\lambda_{2}+\lambda_{3})^2}{(9 \lambda_{1}^2)}\]
AK can be sampled from the principal axis of the diffusion tensor:
\[AK = K(\mathbf{\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.
References
[1]
(1,2,3)
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of
non-Gaussian water diffusion by kurtosis analysis. NMR in
Biomedicine 23(7): 698-710
[2]
(1,2,3,4)
Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011.
Estimation of tensors and tensor-derived measures in diffusional
kurtosis imaging. Magn Reson Med. 65(3), 823-836
[3]
Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging:
Robust estimation from DW-MRI using homogeneous polynomials.
Proceedings of the 8th {IEEE} International Symposium on
Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265.
doi: 10.1109/ISBI.2011.5872402
carlson_rd
-
dipy.reconst.dki.carlson_rd(x, y, z, errtol=0.0001)
Computes the Carlson’s incomplete elliptic integral of the second kind
defined as:
\[R_D = \frac{3}{2} \int_{0}^{\infty} (t+x)^{-\frac{1}{2}}
(t+y)^{-\frac{1}{2}}(t+z) ^{-\frac{3}{2}}\]
- Parameters:
- xndarray
First independent variable of the integral.
- yndarray
Second independent variable of the integral.
- zndarray
Third independent variable of the integral.
- errtolfloat
Error tolerance. Integral is computed with relative error less in
magnitude than the defined value
- Returns:
- RDndarray
Value of the incomplete second order elliptic integral
Notes
x, y, and z have to be nonnegative and at most x or y is zero.
carlson_rf
-
dipy.reconst.dki.carlson_rf(x, y, z, errtol=0.0003)
Computes the Carlson’s incomplete elliptic integral of the first kind
defined as:
\[R_F = \frac{1}{2} \int_{0}^{\infty} \left [(t+x)(t+y)(t+z) \right ]
^{-\frac{1}{2}}dt\]
- Parameters:
- xndarray
First independent variable of the integral.
- yndarray
Second independent variable of the integral.
- zndarray
Third independent variable of the integral.
- errtolfloat
Error tolerance. Integral is computed with relative error less in
magnitude than the defined value
- Returns:
- RFndarray
Value of the incomplete first order elliptic integral
Notes
x, y, and z have to be nonnegative and at most one of them is zero.
References
[1]
Carlson, B.C., 1994. Numerical computation of real or complex
elliptic integrals. arXiv:math/9409227 [math.CA]
cart2sphere
-
dipy.reconst.dki.cart2sphere(x, y, z)
Return angles for Cartesian 3D coordinates x, y, and z
See doc for sphere2cart
for angle conventions and derivation
of the formulae.
\(0\le\theta\mathrm{(theta)}\le\pi\) and \(-\pi\le\phi\mathrm{(phi)}\le\pi\)
- Parameters:
- xarray_like
x coordinate in Cartesian space
- yarray_like
y coordinate in Cartesian space
- zarray_like
z coordinate
- Returns:
- rarray
radius
- thetaarray
inclination (polar) angle
- phiarray
azimuth angle
check_multi_b
-
dipy.reconst.dki.check_multi_b(gtab, n_bvals, non_zero=True, bmag=None)
Check if you have enough different b-values in your gradient table
- Parameters:
- gtabGradientTable class instance.
- n_bvalsint
The number of different b-values you are checking for.
- non_zerobool
Whether to check only non-zero bvalues. In this case, we will require
at least n_bvals non-zero b-values (where non-zero is defined
depending on the gtab object’s b0_threshold attribute)
- bmagint
The order of magnitude of the b-values used. The function will
normalize the b-values relative \(10^{bmag}\). Default: derive this
value from the maximal b-value provided:
\(bmag=log_{10}(max(bvals)) - 1\).
- Returns:
- boolWhether there are at least n_bvals different b-values in the
- gradient table used.
decompose_tensor
-
dipy.reconst.dki.decompose_tensor(tensor, min_diffusivity=0)
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
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])
design_matrix
-
dipy.reconst.dki.design_matrix(gtab)
Construct B design matrix for DKI.
- Parameters:
- gtabGradientTable
Measurement directions.
- Returns:
- Barray (N, 22)
Design matrix or B matrix for the DKI model
B[j, :] = (Bxx, Bxy, Bzz, Bxz, Byz, Bzz,
Bxxxx, Byyyy, Bzzzz, Bxxxy, Bxxxz,
Bxyyy, Byyyz, Bxzzz, Byzzz, Bxxyy,
Bxxzz, Byyzz, Bxxyz, Bxyyz, Bxyzz,
BlogS0)
directional_diffusion
-
dipy.reconst.dki.directional_diffusion(dt, V, min_diffusivity=0)
Calculates the apparent diffusion coefficient (adc) in each direction
of a sphere for a single voxel [1].
- 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.
References
[1]
(1,2)
Neto Henriques R, Correia MM, Nunes RG, Ferreira HA (2015).
Exploring the 3D geometry of the diffusion kurtosis tensor -
Impact on the development of robust tractography procedures and
novel biomarkers, NeuroImage 111: 85-99
directional_diffusion_variance
-
dipy.reconst.dki.directional_diffusion_variance(kt, V, min_kurtosis=-0.42857142857142855)
Calculates the apparent diffusion variance (adv) in each direction
of a sphere for a single voxel [1].
- Parameters:
- dtarray (6,)
elements of the diffusion tensor 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_kurtosisfloat (optional)
Because high-amplitude negative values of kurtosis are not physicaly
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 [2]_)
- adcndarray(g,) (optional)
Apparent diffusion coefficient (adc) in all g directions of a sphere
for a single voxel.
- advndarray(g,) (optional)
Apparent diffusion variance coefficient (advc) in all g directions of
a sphere for a single voxel.
- Returns:
- advndarray (g,)
Apparent diffusion variance (adv) in all g directions of a sphere for
a single voxel.
References
[1]
(1,2)
Neto Henriques R, Correia MM, Nunes RG, Ferreira HA (2015).
Exploring the 3D geometry of the diffusion kurtosis tensor -
Impact on the development of robust tractography procedures and
novel biomarkers, NeuroImage 111: 85-99
directional_kurtosis
-
dipy.reconst.dki.directional_kurtosis(dt, md, kt, V, min_diffusivity=0, min_kurtosis=-0.42857142857142855, adc=None, adv=None)
Calculates the apparent kurtosis coefficient (akc) in each direction
of a sphere for a single voxel [1].
- 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 physicaly
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 [2])
- 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.
References
[1]
(1,2)
Neto Henriques R, Correia MM, Nunes RG, Ferreira HA (2015).
Exploring the 3D geometry of the diffusion kurtosis tensor -
Impact on the development of robust tractography procedures and
novel biomarkers, NeuroImage 111: 85-99
[2]
Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging:
Robust estimation from DW-MRI using homogeneous polynomials.
Proceedings of the 8th {IEEE} International Symposium on
Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265.
doi: 10.1109/ISBI.2011.5872402
dki_prediction
-
dipy.reconst.dki.dki_prediction(dki_params, gtab, S0=1.0)
Predict a signal given diffusion kurtosis imaging parameters.
- 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
- gtaba GradientTable class instance
The gradient table for this prediction
- S0float or ndarray (optional)
The non diffusion-weighted signal in every voxel, or across all
voxels. Default: 1
- Returns:
- S(…, N) ndarray
Simulated signal based on the DKI model:
\[S=S_{0}e^{-bD+\]
- rac{1}{6}b^{2}D^{2}K}
from_lower_triangular
-
dipy.reconst.dki.from_lower_triangular(D)
Returns a tensor given the six unique tensor elements
Given the six unique tensor elements (in the order: Dxx, Dxy, Dyy, Dxz,
Dyz, Dzz) returns a 3 by 3 tensor. All elements after the sixth are
ignored.
- Parameters:
- Darray_like, (…, >6)
Unique elements of the tensors
- Returns:
- tensorndarray (…, 3, 3)
3 by 3 tensors
get_fnames
-
dipy.reconst.dki.get_fnames(name='small_64D')
Provide full paths to example or test datasets.
- Parameters:
- namestr
the filename/s of which dataset to return, one of:
‘small_64D’ small region of interest nifti,bvecs,bvals 64 directions
‘small_101D’ small region of interest nifti, bvecs, bvals
101 directions
‘aniso_vox’ volume with anisotropic voxel size as Nifti
‘fornix’ 300 tracks in Trackvis format (from Pittsburgh
Brain Competition)
‘gqi_vectors’ the scanner wave vectors needed for a GQI acquisitions
of 101 directions tested on Siemens 3T Trio
‘small_25’ small ROI (10x8x2) DTI data (b value 2000, 25 directions)
‘test_piesno’ slice of N=8, K=14 diffusion data
‘reg_c’ small 2D image used for validating registration
‘reg_o’ small 2D image used for validation registration
‘cb_2’ two vectorized cingulum bundles
- Returns:
- fnamestuple
filenames for dataset
Examples
>>> import numpy as np
>>> from dipy.io.image import load_nifti
>>> from dipy.data import get_fnames
>>> fimg, fbvals, fbvecs = get_fnames('small_101D')
>>> bvals=np.loadtxt(fbvals)
>>> bvecs=np.loadtxt(fbvecs).T
>>> data, affine = load_nifti(fimg)
>>> data.shape == (6, 10, 10, 102)
True
>>> bvals.shape == (102,)
True
>>> bvecs.shape == (102, 3)
True
get_sphere
-
dipy.reconst.dki.get_sphere(name='symmetric362')
provide triangulated spheres
- Parameters:
- namestr
which sphere - one of:
* ‘symmetric362’
* ‘symmetric642’
* ‘symmetric724’
* ‘repulsion724’
* ‘repulsion100’
* ‘repulsion200’
- Returns:
- spherea dipy.core.sphere.Sphere class instance
Examples
>>> import numpy as np
>>> from dipy.data import get_sphere
>>> sphere = get_sphere('symmetric362')
>>> verts, faces = sphere.vertices, sphere.faces
>>> verts.shape == (362, 3)
True
>>> faces.shape == (720, 3)
True
>>> verts, faces = get_sphere('not a sphere name')
Traceback (most recent call last):
...
DataError: No sphere called "not a sphere name"
kurtosis_fractional_anisotropy
-
dipy.reconst.dki.kurtosis_fractional_anisotropy(dki_params)
Computes the anisotropy of the kurtosis tensor (KFA) [1].
- 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
- ——-
- kfaarray
Calculated mean kurtosis tensor.
Notes
The KFA is defined as [1]:
\[KFA \equiv
\frac{||\mathbf{W} - MKT \mathbf{I}^{(4)}||_F}{||\mathbf{W}||_F}\]
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 [1].
References
[1]
(1,2,3,4)
Glenn, G. R., Helpern, J. A., Tabesh, A., and Jensen, J. H. (2015).
Quantitative assessment of diffusional kurtosis anisotropy.
NMR in Biomedicine 28, 448–459. doi:10.1002/nbm.3271
kurtosis_maximum
-
dipy.reconst.dki.kurtosis_maximum(dki_params, sphere='repulsion100', gtol=0.01, mask=None)
Computes kurtosis maximum value
- 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 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
local_maxima
-
dipy.reconst.dki.local_maxima()
Local maxima of a function evaluated on a discrete set of points.
If a function is evaluated on some set of points where each pair of
neighboring points is an edge in edges, find the local maxima.
- Parameters:
- odfarray, 1d, dtype=double
The function evaluated on a set of discrete points.
- edgesarray (N, 2)
The set of neighbor relations between the points. Every edge, ie
edges[i, :], is a pair of neighboring points.
- Returns:
- peak_valuesndarray
Value of odf at a maximum point. Peak values is sorted in descending
order.
- peak_indicesndarray
Indices of maximum points. Sorted in the same order as peak_values so
odf[peak_indices[i]] == peak_values[i].
Notes
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).
lower_triangular
-
dipy.reconst.dki.lower_triangular(tensor, b0=None)
Returns the six lower triangular values of the tensor ordered as
(Dxx, Dxy, Dyy, Dxz, Dyz, Dzz) and a dummy variable if b0 is not None.
- Parameters:
- tensorarray_like (…, 3, 3)
a collection of 3, 3 diffusion tensors
- b0float
if b0 is not none log(b0) is returned as the dummy variable
- Returns:
- Dndarray
If b0 is none, then the shape will be (…, 6) otherwise (…, 7)
mean_diffusivity
-
dipy.reconst.dki.mean_diffusivity(evals, axis=-1)
Mean Diffusivity (MD) of a diffusion tensor.
- Parameters:
- evalsarray-like
Eigenvalues of a diffusion tensor.
- axisint
Axis of evals which contains 3 eigenvalues.
- Returns:
- mdarray
Calculated MD.
Notes
MD is calculated with the following equation:
\[MD = \frac{\lambda_1 + \lambda_2 + \lambda_3}{3}\]
mean_kurtosis
-
dipy.reconst.dki.mean_kurtosis(dki_params, min_kurtosis=-0.42857142857142855, max_kurtosis=3, analytical=True)
Computes mean Kurtosis (MK) from the kurtosis tensor.
- 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 [4])
- 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. Default = 10
- analyticalbool (optional)
If True, MK is calculated using its analytical solution, otherwise an
exact numerical estimator is used (see Notes). Default is set to True
- 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[R953e26c55b6a-1]_:
\[MK \equiv \frac{1}{4\pi} \int d\Omega_\mathbf{n} K(\mathbf{n})\]
This integral can be numerically solved by averaging directional
kurtosis values sampled for directions of a spherical t-design [2].
Alternatively, MK can be solved from the analytical solution derived by
Tabesh et al. [3]. This solution is given by:
\[\begin{split}MK=F_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{1111}+
F_1(\lambda_2,\lambda_1,\lambda_3)\hat{W}_{2222}+
F_1(\lambda_3,\lambda_2,\lambda_1)\hat{W}_{3333}+ \\
F_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233}+
F_2(\lambda_2,\lambda_1,\lambda_3)\hat{W}_{1133}+
F_2(\lambda_3,\lambda_2,\lambda_1)\hat{W}_{1122}\end{split}\]
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
\[ \begin{align}\begin{aligned}\begin{split}F_1(\lambda_1,\lambda_2,\lambda_3)=
\frac{(\lambda_1+\lambda_2+\lambda_3)^2}
{18(\lambda_1-\lambda_2)(\lambda_1-\lambda_3)}
[\frac{\sqrt{\lambda_2\lambda_3}}{\lambda_1}
R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\
\frac{3\lambda_1^2-\lambda_1\lambda_2-\lambda_2\lambda_3-
\lambda_1\lambda_3}
{3\lambda_1 \sqrt{\lambda_2 \lambda_3}}
R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-1 ]\end{split}\\\begin{split}F_2(\lambda_1,\lambda_2,\lambda_3)=
\frac{(\lambda_1+\lambda_2+\lambda_3)^2}
{3(\lambda_2-\lambda_3)^2}
[\frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}}
R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\
\frac{2\lambda_1-\lambda_2-\lambda_3}{3\sqrt{\lambda_2 \lambda_3}}
R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-2]\end{split}\end{aligned}\end{align} \]
where \(R_f\) and \(R_d\) are the Carlson’s elliptic integrals.
References
[1]
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of
non-Gaussian water diffusion by kurtosis analysis. NMR in
Biomedicine 23(7): 698-710
[2]
Hardin, R.H., Sloane, N.J.A., 1996. McLaren’s Improved Snub Cube and
Other New Spherical Designs in Three Dimensions. Discrete and
Computational Geometry 15, 429-441.
[3]
Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011.
Estimation of tensors and tensor-derived measures in diffusional
kurtosis imaging. Magn Reson Med. 65(3), 823-836
[4]
Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging:
Robust estimation from DW-MRI using homogeneous polynomials.
Proceedings of the 8th {IEEE} International Symposium on
Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265.
doi: 10.1109/ISBI.2011.5872402
mean_kurtosis_tensor
-
dipy.reconst.dki.mean_kurtosis_tensor(dki_params, min_kurtosis=-0.42857142857142855, max_kurtosis=10)
Computes mean of the kurtosis tensor (MKT) [1].
- 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 [2])
- 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. Default = 10
- Returns
- ——-
- mktarray
Calculated mean kurtosis tensor.
Notes
The MKT is defined as [1]:
\[MKT \equiv \frac{1}{4\pi} \int d
\Omega_{\mathnbf{n}} n_i n_j n_k n_l W_{ijkl}\]
which can be directly computed from the trace of the kurtosis tensor:
\[\]
MKT = frac{1}{5} Tr(mathbf{W}) = frac{1}{5}
(W_{1111} + W_{2222} + W_{3333} + 2W_{1122} + 2W_{1133} + 2W_{2233})
References
[1]
(1,2,3)
Hansen, B., Lund, T. E., Sangill, R., and Jespersen, S. N. (2013).
Experimentally and computationally fast method for estimation of
a mean kurtosis.Magnetic Resonance in Medicine69, 1754–1760.388
doi:10.1002/mrm.24743
[2]
Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging:
Robust estimation from DW-MRI using homogeneous polynomials.
Proceedings of the 8th {IEEE} International Symposium on
Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265.
doi: 10.1109/ISBI.2011.5872402
ndindex
-
dipy.reconst.dki.ndindex(shape)
An N-dimensional iterator object to index arrays.
Given the shape of an array, an ndindex instance iterates over
the N-dimensional index of the array. At each iteration a tuple
of indices is returned; the last dimension is iterated over first.
- Parameters:
- shapetuple of ints
The dimensions of the array.
Examples
>>> from dipy.core.ndindex import ndindex
>>> shape = (3, 2, 1)
>>> for index in ndindex(shape):
... print(index)
(0, 0, 0)
(0, 1, 0)
(1, 0, 0)
(1, 1, 0)
(2, 0, 0)
(2, 1, 0)
nlls_fit_tensor
-
dipy.reconst.dki.nlls_fit_tensor(design_matrix, data, weighting=None, sigma=None, jac=True, return_S0_hat=False, fail_is_nan=False)
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.
- weighting: str
the weighting scheme to use in considering the
squared-error. Default behavior is to use uniform weighting. Other
options: ‘sigma’ ‘gmm’
- sigma: float
If the ‘sigma’ weighting scheme is used, a value of sigma needs to be
provided here. According to [Chang2005], 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).
- jacbool
Use the Jacobian? Default: True
- return_S0_hatbool
Boolean to return (True) or not (False) the S0 values for the fit.
- fail_is_nanbool
Boolean to set failed NL fitting to NaN (True) or LS (False, default).
- Returns:
- nlls_params: the eigen-values and eigen-vectors of the tensor in each
voxel.
ols_fit_dki
-
dipy.reconst.dki.ols_fit_dki(design_matrix, data)
Computes the diffusion and kurtosis tensors using an ordinary linear
least squares (OLS) approach .
- Parameters:
- design_matrixarray (g, 22)
Design matrix holding the covariants used to solve for the regression
coefficients.
- dataarray (N, g)
Data or response variables holding the data. Note that the last
dimension should contain the data. It makes no copies of data.
- Returns:
- dki_paramsarray (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
References
- [1] Lu, H., Jensen, J. H., Ramani, A., and Helpern, J. A. (2006).
Three-dimensional characterization of non-gaussian water diffusion in
humans using diffusion kurtosis imaging. NMR in Biomedicine 19,
236–247. doi:10.1002/nbm.1020
perpendicular_directions
-
dipy.reconst.dki.perpendicular_directions(v, num=30, half=False)
Computes n evenly spaced perpendicular directions relative to a given
vector v
- Parameters:
- varray (3,)
Array containing the three cartesian coordinates of vector v
- numint, optional
Number of perpendicular directions to generate
- halfbool, optional
If half is True, perpendicular directions are sampled on half of the
unit circumference perpendicular to v, otherwive perpendicular
directions are sampled on the full circumference. Default of half is
False
- Returns:
- psamplesarray (n, 3)
array of vectors perpendicular to v
Notes
Perpendicular directions are estimated using the following two step
procedure:
1) the perpendicular directions are first sampled in a unit
circumference parallel to the plane normal to the x-axis.
2) Samples are then rotated and aligned to the plane normal to vector
v. The rotational matrix for this rotation is constructed as reference
frame basis which axis are the following:
cross product between vector v and the unit vector aligned to the
x-axis
- The third axis is defined as the cross product between the
previous computed vector and vector v.
Following this two steps, coordinates of the final perpendicular directions
are given as:
\[\left [ -\sin(a_{i}) \sqrt{{v_{y}}^{2}+{v_{z}}^{2}}
\; , \;
\frac{v_{x}v_{y}\sin(a_{i})-v_{z}\cos(a_{i})}
{\sqrt{{v_{y}}^{2}+{v_{z}}^{2}}}
\; , \;
\frac{v_{x}v_{z}\sin(a_{i})-v_{y}\cos(a_{i})}
{\sqrt{{v_{y}}^{2}+{v_{z}}^{2}}} \right ]\]
This procedure has a singularity when vector v is aligned to the x-axis. To
solve this singularity, perpendicular directions in procedure’s step 1 are
defined in the plane normal to y-axis and the second axis of the rotated
frame of reference is computed as the normalized vector given by the cross
product between vector v and the unit vector aligned to the y-axis.
Following this, the coordinates of the perpendicular directions are given
as:
left [ -frac{left (v_{x}v_{y}sin(a_{i})+v_{z}cos(a_{i}) right )}
{sqrt{{v_{x}}^{2}+{v_{z}}^{2}}}
; , ;
sin(a_{i}) sqrt{{v_{x}}^{2}+{v_{z}}^{2}}
; , ;
frac{v_{y}v_{z}sin(a_{i})+v_{x}cos(a_{i})}
{sqrt{{v_{x}}^{2}+{v_{z}}^{2}}} right ]
For more details on this calculation, see ` here <http://gsoc2015dipydki.blogspot.it/2015/07/rnh-post-8-computing-perpendicular.html>`_.
radial_kurtosis
-
dipy.reconst.dki.radial_kurtosis(dki_params, min_kurtosis=-0.42857142857142855, max_kurtosis=10, analytical=True)
Radial Kurtosis (RK) of a diffusion kurtosis tensor [1], [2].
- 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 [3])
- 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. Default = 10
- 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 [1], [2]:
\[\]
- RK equiv frac{1}{2pi} int dOmega _mathbf{theta} K(mathbf{theta})
delta (mathbf{theta}cdot mathbf{e}_1)
This equation can be numerically computed by averaging apparent
directional kurtosis samples for directions perpendicular to e1.
Otherwise, RK can be calculated from its analytical solution [2]:
\[K_{\bot} = G_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2222} +
G_1(\lambda_1,\lambda_3,\lambda_2)\hat{W}_{3333} +
G_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233}\]
where:
\[G_1(\lambda_1,\lambda_2,\lambda_3)=
\frac{(\lambda_1+\lambda_2+\lambda_3)^2}{18\lambda_2(\lambda_2-
\lambda_3)} \left (2\lambda_2 +
\frac{\lambda_3^2-3\lambda_2\lambda_3}{\sqrt{\lambda_2\lambda_3}}
\right)\]
and
\[G_2(\lambda_1,\lambda_2,\lambda_3)=
\frac{(\lambda_1+\lambda_2+\lambda_3)^2}{(\lambda_2-\lambda_3)^2}
\left ( \frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}}-2\right )\]
References
[1]
(1,2,3)
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of
non-Gaussian water diffusion by kurtosis analysis. NMR in
Biomedicine 23(7): 698-710
[2]
(1,2,3,4)
Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011.
Estimation of tensors and tensor-derived measures in diffusional
kurtosis imaging. Magn Reson Med. 65(3), 823-836
[3]
Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging:
Robust estimation from DW-MRI using homogeneous polynomials.
Proceedings of the 8th {IEEE} International Symposium on Biomedical
Imaging: From Nano to Macro, ISBI 2011, 262-265.
doi: 10.1109/ISBI.2011.5872402
restore_fit_tensor
-
dipy.reconst.dki.restore_fit_tensor(design_matrix, data, sigma=None, jac=True, return_S0_hat=False, fail_is_nan=False)
Use the RESTORE algorithm [1] to calculate a robust tensor fit
- 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, array of shape [n_directions], array of shape [X, Y, Z]
An estimate of the variance. [1] 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).
Array with ndim > 1 corresponds to spatially varying sigma, so if
providing spatially-flattened data and spatially-varying sigma,
provide array with shape [num_vox, 1].
- 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()
). Default: True
- return_S0_hatbool
Boolean to return (True) or not (False) the S0 values for the fit.
- fail_is_nanbool
Boolean to set failed NL fitting to NaN (True) or LS (False, default).
- Returns:
- restore_paramsan estimate of the tensor parameters in each voxel.
References
[1]
(1,2,3)
Chang, L-C, Jones, DK and Pierpaoli, C (2005). RESTORE: robust
estimation of tensors by outlier rejection. MRM, 53: 1088-95.
sphere2cart
-
dipy.reconst.dki.sphere2cart(r, theta, phi)
Spherical to Cartesian coordinates
This is the standard physics convention where theta is the
inclination (polar) angle, and phi is the azimuth angle.
Imagine a sphere with center (0,0,0). Orient it with the z axis
running south-north, the y axis running west-east and the x axis
from posterior to anterior. theta (the inclination angle) is the
angle to rotate from the z-axis (the zenith) around the y-axis,
towards the x axis. Thus the rotation is counter-clockwise from the
point of view of positive y. phi (azimuth) gives the angle of
rotation around the z-axis towards the y axis. The rotation is
counter-clockwise from the point of view of positive z.
Equivalently, given a point P on the sphere, with coordinates x, y,
z, theta is the angle between P and the z-axis, and phi is
the angle between the projection of P onto the XY plane, and the X
axis.
Geographical nomenclature designates theta as ‘co-latitude’, and phi
as ‘longitude’
- Parameters:
- rarray_like
radius
- thetaarray_like
inclination or polar angle
- phiarray_like
azimuth angle
- Returns:
- xarray
x coordinate(s) in Cartesion space
- yarray
y coordinate(s) in Cartesian space
- zarray
z coordinate
Notes
See these pages:
for excellent discussion of the many different conventions
possible. Here we use the physics conventions, used in the
wikipedia page.
Derivations of the formulae are simple. Consider a vector x, y, z of
length r (norm of x, y, z). The inclination angle (theta) can be
found from: cos(theta) == z / r -> z == r * cos(theta). This gives
the hypotenuse of the projection onto the XY plane, which we will
call Q. Q == r*sin(theta). Now x / Q == cos(phi) -> x == r *
sin(theta) * cos(phi) and so on.
We have deliberately named this function sphere2cart
rather than
sph2cart
to distinguish it from the Matlab function of that
name, because the Matlab function uses an unusual convention for the
angles that we did not want to replicate. The Matlab function is
trivial to implement with the formulae given in the Matlab help.
split_dki_param
-
dipy.reconst.dki.split_dki_param(dki_params)
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])
- ktarray (x, y, z, 15) or (n, 15)
Fifteen elements of the kurtosis tensor
vec_val_vect
-
dipy.reconst.dki.vec_val_vect()
Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs
- Parameters:
- vecsshape (…, M, N) array
containing tensor in last two dimensions; M, N usually equal to (3, 3)
- valsshape (…, N) array
diagonal values carried in last dimension, ...
shape above must
match that for vecs
- Returns:
- resshape (…, M, M) array
For all the dimensions ellided by ...
, loops to get (M, N) vec
matrix, and (N,) vals
vector, and calculates
vec.dot(np.diag(val).dot(vec.T)
.
- Raises:
- ValueErrornon-matching
...
dimensions of vecs, vals
- ValueErrornon-matching
N
dimensions of vecs, vals
Examples
Make a 3D array where the first dimension is only 1
>>> vecs = np.arange(9).reshape((1, 3, 3))
>>> vals = np.arange(3).reshape((1, 3))
>>> vec_val_vect(vecs, vals)
array([[[ 9., 24., 39.],
[ 24., 66., 108.],
[ 39., 108., 177.]]])
That’s the same as the 2D case (apart from the float casting):
>>> vecs = np.arange(9).reshape((3, 3))
>>> vals = np.arange(3)
>>> np.dot(vecs, np.dot(np.diag(vals), vecs.T))
array([[ 9, 24, 39],
[ 24, 66, 108],
[ 39, 108, 177]])
wls_fit_dki
-
dipy.reconst.dki.wls_fit_dki(design_matrix, data)
Computes the diffusion and kurtosis tensors using a weighted linear
least squares (WLS) approach .
- Parameters:
- design_matrixarray (g, 22)
Design matrix holding the covariants used to solve for the regression
coefficients.
- dataarray (N, g)
Data or response variables holding the data. Note that the last
dimension should contain the data. It makes no copies of data.
- Returns:
- dki_paramsarray (N, 27)
All parameters estimated from the diffusion kurtosis model for all N
voxels.
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
References
- [1] Veraart, J., Sijbers, J., Sunaert, S., Leemans, A., Jeurissen, B.,
2013. Weighted linear least squares estimation of diffusion MRI
parameters: Strengths, limitations, and pitfalls. Magn Reson Med 81,
335-346.
-
class dipy.reconst.dki_micro.DiffusionKurtosisFit(model, model_params)
Bases: TensorFit
Class for fitting the Diffusion Kurtosis Model
- Attributes:
- S0_hat
directions
For tracking - return the primary direction in each voxel
evals
Returns the eigenvalues of the tensor as an array
evecs
Returns the eigenvectors of the tensor as an array, columnwise
kfa
Returns the kurtosis tensor (KFA) .
kt
Returns the 15 independent elements of the kurtosis tensor as an array
quadratic_form
Calculates the 3x3 diffusion tensor for each voxel
- shape
Methods
ad ()
|
Axial diffusivity (AD) calculated from cached eigenvalues. |
adc (sphere)
|
Calculate the apparent diffusion coefficient (ADC) in each direction on |
ak ([min_kurtosis, max_kurtosis, analytical])
|
Axial Kurtosis (AK) of a diffusion kurtosis tensor [1]. |
akc (sphere)
|
Calculates 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])
|
Computes the maximum value of a single voxel kurtosis tensor |
linearity ()
|
- Returns:
|
md ()
|
Mean diffusivity (MD) calculated from cached eigenvalues. |
mk ([min_kurtosis, max_kurtosis, analytical])
|
Computes mean Kurtosis (MK) from the kurtosis tensor. |
mkt ([min_kurtosis, max_kurtosis])
|
Computes mean of the kurtosis tensor (MKT) [1]. |
mode ()
|
Tensor mode calculated from cached eigenvalues. |
odf (sphere)
|
The diffusion orientation distribution function (dODF). |
planarity ()
|
- Returns:
|
predict (gtab[, S0])
|
Given a DKI model fit, predict the signal on the vertices of a gradient table |
rd ()
|
Radial diffusivity (RD) calculated from cached eigenvalues. |
rk ([min_kurtosis, max_kurtosis, analytical])
|
Radial Kurtosis (RK) of a diffusion kurtosis tensor [1]. |
sphericity ()
|
- Returns:
|
trace ()
|
Trace of the tensor calculated from cached eigenvalues. |
-
__init__(model, model_params)
Initialize a DiffusionKurtosisFit class instance.
Since DKI is an extension of DTI, class instance is defined as subclass
of the TensorFit from dti.py
- Parameters:
- modelDiffusionKurtosisModel Class instance
Class instance containing the Diffusion Kurtosis Model for the fit
- model_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
-
ak(min_kurtosis=-0.42857142857142855, max_kurtosis=10, analytical=True)
Axial Kurtosis (AK) of a diffusion kurtosis tensor [1].
- 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 [2])
- 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. Default = 10
- 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). Default is set to True.
- Returns:
- akarray
Calculated AK.
Notes
AK is defined as the directional kurtosis parallel to the fiber’s main
direction e1 [1], [2]. You can compute AK using to approaches:
AK is calculated from rotated diffusion kurtosis tensor [2], i.e.:
\[AK = \hat{W}_{1111}
\frac{(\lambda_{1}+\lambda_{2}+\lambda_{3})^2}{(9 \lambda_{1}^2)}\]
AK can be sampled from the principal axis of the diffusion tensor:
\[AK = K(\mathbf{\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.
References
[1]
(1,2,3)
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of
non-Gaussian water diffusion by kurtosis analysis. NMR in
Biomedicine 23(7): 698-710
[2]
(1,2,3)
Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011.
Estimation of tensors and tensor-derived measures in diffusional
kurtosis imaging. Magn Reson Med. 65(3), 823-836
[3]
Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging:
Robust estimation from DW-MRI using homogeneous polynomials.
Proceedings of the 8th {IEEE} International Symposium on
Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265.
doi: 10.1109/ISBI.2011.5872402
-
akc(sphere)
Calculates the apparent kurtosis coefficient (AKC) in each
direction on the sphere for each voxel in the data
- Parameters:
- sphereSphere class instance
- Returns:
- akcndarray
The estimates of the apparent kurtosis coefficient in every
direction on the input sphere
Notes
For each sphere direction with coordinates \((n_{1}, n_{2}, n_{3})\), the
calculation of AKC is done using formula:
\[AKC(n)=\frac{MD^{2}}{ADC(n)^{2}}\sum_{i=1}^{3}\sum_{j=1}^{3}
\sum_{k=1}^{3}\sum_{l=1}^{3}n_{i}n_{j}n_{k}n_{l}W_{ijkl}\]
where \(W_{ijkl}\) are the elements of the kurtosis tensor, MD the mean
diffusivity and ADC the apparent diffusion coefficent computed as:
\[ADC(n)=\sum_{i=1}^{3}\sum_{j=1}^{3}n_{i}n_{j}D_{ij}\]
where \(D_{ij}\) are the elements of the diffusion tensor.
-
property kfa
Returns the kurtosis tensor (KFA) [1].
Notes
The KFA is defined as [1]:
\[KFA \equiv
\frac{||\mathbf{W} - MKT \mathbf{I}^{(4)}||_F}{||\mathbf{W}||_F}\]
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 [1].
References
[1]
(1,2,3)
Glenn, G. R., Helpern, J. A., Tabesh, A., and Jensen, J. H.
(2015). Quantitative assessment of diffusional kurtosis
anisotropy. NMR in Biomedicine 28, 448–459. doi:10.1002/nbm.3271
-
kmax(sphere='repulsion100', gtol=1e-05, mask=None)
Computes the maximum value of a single voxel kurtosis tensor
- Parameters:
- sphereSphere class instance, optional
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
- Returns:
- max_valuefloat
kurtosis tensor maximum value
-
property kt
Returns the 15 independent elements of the kurtosis tensor as an array
-
mk(min_kurtosis=-0.42857142857142855, max_kurtosis=10, analytical=True)
Computes mean Kurtosis (MK) from the kurtosis tensor.
- 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 [4])
- 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. Default = 10
- analyticalbool (optional)
If True, MK is calculated using its analytical solution, otherwise
an exact numerical estimator is used (see Notes). Default is set to
True.
- 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[Rb657f27beb9e-1]_:
\[MK \equiv \frac{1}{4\pi} \int d\Omega_\mathbf{n} K(\mathbf{n})\]
This integral can be numerically solved by averaging directional
kurtosis values sampled for directions of a spherical t-design [2].
Alternatively, MK can be solved from the analytical solution derived by
Tabesh et al. [3]. This solution is given by:
\[\begin{split}MK=F_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{1111}+
F_1(\lambda_2,\lambda_1,\lambda_3)\hat{W}_{2222}+
F_1(\lambda_3,\lambda_2,\lambda_1)\hat{W}_{3333}+ \\
F_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233}+
F_2(\lambda_2,\lambda_1,\lambda_3)\hat{W}_{1133}+
F_2(\lambda_3,\lambda_2,\lambda_1)\hat{W}_{1122}\end{split}\]
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
\[ \begin{align}\begin{aligned}\begin{split}F_1(\lambda_1,\lambda_2,\lambda_3)=
\frac{(\lambda_1+\lambda_2+\lambda_3)^2}
{18(\lambda_1-\lambda_2)(\lambda_1-\lambda_3)}
[\frac{\sqrt{\lambda_2\lambda_3}}{\lambda_1}
R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\
\frac{3\lambda_1^2-\lambda_1\lambda_2-\lambda_2\lambda_3-
\lambda_1\lambda_3}
{3\lambda_1 \sqrt{\lambda_2 \lambda_3}}
R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-1 ]\end{split}\\\begin{split}F_2(\lambda_1,\lambda_2,\lambda_3)=
\frac{(\lambda_1+\lambda_2+\lambda_3)^2}
{3(\lambda_2-\lambda_3)^2}
[\frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}}
R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\
\frac{2\lambda_1-\lambda_2-\lambda_3}{3\sqrt{\lambda_2 \lambda_3}}
R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-2]\end{split}\end{aligned}\end{align} \]
where \(R_f\) and \(R_d\) are the Carlson’s elliptic integrals.
References
[1]
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of
non-Gaussian water diffusion by kurtosis analysis. NMR in
Biomedicine 23(7): 698-710
[2]
Hardin, R.H., Sloane, N.J.A., 1996. McLaren’s Improved Snub Cube
and Other New Spherical Designs in Three Dimensions. Discrete
and Computational Geometry 15, 429-441.
[3]
Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011.
Estimation of tensors and tensor-derived measures in diffusional
kurtosis imaging. Magn Reson Med. 65(3), 823-836
[4]
Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging:
Robust estimation from DW-MRI using homogeneous polynomials.
Proceedings of the 8th {IEEE} International Symposium on
Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265.
doi: 10.1109/ISBI.2011.5872402
-
mkt(min_kurtosis=-0.42857142857142855, max_kurtosis=10)
Computes mean of the kurtosis tensor (MKT) [1].
- 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 [2])
- 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. Default = 10
- Returns:
- mktarray
Calculated mean kurtosis tensor.
Notes
The MKT is defined as [1]:
\[MKT \equiv \frac{1}{4\pi} \int d
\Omega_{\mathnbf{n}} n_i n_j n_k n_l W_{ijkl}\]
which can be directly computed from the trace of the kurtosis tensor:
\[\]
MKT = frac{1}{5} Tr(mathbf{W}) = frac{1}{5}
(W_{1111} + W_{2222} + W_{3333} + 2W_{1122} + 2W_{1133} + 2W_{2233})
References
[1]
(1,2,3)
Hansen, B., Lund, T. E., Sangill, R., and Jespersen, S. N. 2013.
Experimentally and computationally fast method for estimation
of a mean kurtosis. Magnetic Resonance in Medicine69, 1754–1760.
388. doi:10.1002/mrm.24743
[2]
Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging:
Robust estimation from DW-MRI using homogeneous polynomials.
Proceedings of the 8th {IEEE} International Symposium on
Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265.
doi: 10.1109/ISBI.2011.5872402
-
predict(gtab, S0=1.0)
Given a DKI model fit, predict the signal on the vertices of a
gradient table
- Parameters:
- gtaba GradientTable class instance
The gradient table for this prediction
- S0float or ndarray (optional)
The non diffusion-weighted signal in every voxel, or across all
voxels. Default: 1
Notes
The predicted signal is given by:
\[S(n,b)=S_{0}e^{-bD(n)+\frac{1}{6}b^{2}D(n)^{2}K(n)}\]
\(\mathbf{D(n)}\) and \(\mathbf{K(n)}\) can be computed from the DT and KT
using the following equations:
\[D(n)=\sum_{i=1}^{3}\sum_{j=1}^{3}n_{i}n_{j}D_{ij}\]
and
\[K(n)=\frac{MD^{2}}{D(n)^{2}}\sum_{i=1}^{3}\sum_{j=1}^{3}
\sum_{k=1}^{3}\sum_{l=1}^{3}n_{i}n_{j}n_{k}n_{l}W_{ijkl}\]
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.
-
rk(min_kurtosis=-0.42857142857142855, max_kurtosis=10, analytical=True)
Radial Kurtosis (RK) of a diffusion kurtosis tensor [1].
- Parameters:
- 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
[3])
- 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. Default = 10
- 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 [1], [2]:
\[\]
- RK equiv frac{1}{2pi} int dOmega _mathbf{theta}
K(mathbf{theta}) delta (mathbf{theta}cdot mathbf{e}_1)
This equation can be numerically computed by averaging apparent
directional kurtosis samples for directions perpendicular to e1.
Otherwise, RK can be calculated from its analytical solution [2]:
\[K_{\bot} = G_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2222} +
G_1(\lambda_1,\lambda_3,\lambda_2)\hat{W}_{3333} +
G_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233}\]
where:
\[G_1(\lambda_1,\lambda_2,\lambda_3)=
\frac{(\lambda_1+\lambda_2+\lambda_3)^2}{18\lambda_2(\lambda_2-
\lambda_3)} \left (2\lambda_2 +
\frac{\lambda_3^2-3\lambda_2\lambda_3}{\sqrt{\lambda_2\lambda_3}}
\right)\]
and
\[ G_2(\lambda_1,\lambda_2,\lambda_3)=
\frac{(\lambda_1+\lambda_2+\lambda_3)^2}{(\lambda_2-\lambda_3)^2}
\left ( \frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}}-
2\right )\]
References
[1]
(1,2,3)
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of
non-Gaussian water diffusion by kurtosis analysis. NMR in
Biomedicine 23(7): 698-710
[2]
(1,2)
Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011.
Estimation of tensors and tensor-derived measures in diffusional
kurtosis imaging. Magn Reson Med. 65(3), 823-836
[3]
Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging:
Robust estimation from DW-MRI using homogeneous polynomials.
Proceedings of the 8th {IEEE} International Symposium on
Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265.
doi: 10.1109/ISBI.2011.5872402
-
class dipy.reconst.dki_micro.DiffusionKurtosisModel(gtab, fit_method='WLS', *args, **kwargs)
Bases: ReconstModel
Class for the Diffusion Kurtosis Model
Methods
fit (data[, mask])
|
Fit method of the DKI model class |
predict (dki_params[, S0])
|
Predict a signal for this DKI model class instance given parameters. |
-
__init__(gtab, fit_method='WLS', *args, **kwargs)
Diffusion Kurtosis Tensor Model [1]
- Parameters:
- gtabGradientTable class instance
- fit_methodstr or callable
str can be one of the following:
‘OLS’ or ‘ULLS’ for ordinary least squares
- ‘WLS’ or ‘UWLLS’ for weighted ordinary least squares
dki.wls_fit_dki
- callable has to have the signature:
fit_method(design_matrix, data, *args, **kwargs)
- args, kwargsarguments and key-word arguments passed to the
fit_method. See dki.ols_fit_dki, dki.wls_fit_dki for details
References
[1]
Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011.
Estimation of tensors and tensor-derived measures in diffusional
kurtosis imaging. Magn Reson Med. 65(3), 823-836
-
fit(data, mask=None)
Fit method of the DKI model class
- Parameters:
- dataarray
The measured signal from one voxel.
- maskarray
A boolean array used to mark the coordinates in the data that
should be analyzed that has the shape data.shape[-1]
-
predict(dki_params, S0=1.0)
Predict a signal for this DKI model class instance given parameters.
- 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
2. Three lines of the eigenvector matrix each containing the
first, second and third coordinates of the eigenvector
3. Fifteen elements of the kurtosis tensor
- S0float or ndarray (optional)
The non diffusion-weighted signal in every voxel, or across all
voxels. Default: 1
-
class dipy.reconst.dki_micro.KurtosisMicrostructuralFit(model, model_params)
Bases: DiffusionKurtosisFit
Class for fitting the Diffusion Kurtosis Microstructural Model
- Attributes:
- S0_hat
awf
Returns the volume fraction of the restricted diffusion compartment also known as axonal water fraction.
axonal_diffusivity
Returns the axonal diffusivity defined as the restricted diffusion tensor trace .
directions
For tracking - return the primary direction in each voxel
evals
Returns the eigenvalues of the tensor as an array
evecs
Returns the eigenvectors of the tensor as an array, columnwise
hindered_ad
Returns the axial diffusivity of the hindered compartment.