reconstreconst.basereconst.benchmarksreconst.benchmarks.bench_bounding_boxreconst.benchmarks.bench_csdreconst.benchmarks.bench_peaksreconst.benchmarks.bench_squashreconst.benchmarks.bench_vec_val_sumreconst.cachereconst.cross_validationreconst.csdeconvreconst.dkireconst.dki_microreconst.dsireconst.dtireconst.forecastreconst.fwdtireconst.gqireconst.ivimreconst.mapmrireconst.mcsdreconst.msdkireconst.multi_voxelreconst.odfreconst.peak_direction_getterreconst.qtdmrireconst.sfmreconst.shm
reconst.shorereconst.utilsReconstFitReconstModelConstrainedSphericalDeconvModelGradientTableCacheAxSymShResponseConstrainedSDTModelConstrainedSphericalDeconvModelSphHarmFitSphHarmModelTensorModelDiffusionKurtosisFitDiffusionKurtosisModelReconstModelTensorFitDiffusionKurtosisFitDiffusionKurtosisModelKurtosisMicrostructuralFitKurtosisMicrostructureModelCacheDiffusionSpectrumDeconvFitDiffusionSpectrumDeconvModelDiffusionSpectrumFitDiffusionSpectrumModelOdfFitOdfModelReconstModelTensorFitTensorModelCacheForecastFitForecastModelLooseVersionOdfFitOdfModelFreeWaterTensorFitFreeWaterTensorModelReconstModelTensorFitCacheGeneralizedQSamplingFitGeneralizedQSamplingModelOdfFitOdfModelIvimFitIvimModelTRRIvimModelVPLooseVersionReconstModelCacheLooseVersionMapmriFitMapmriModelOptimizerReconstFitReconstModelGradientTableLooseVersionMSDeconvFitMultiShellDeconvModelMultiShellResponseQpFitterTensorModelMeanDiffusionKurtosisFitMeanDiffusionKurtosisModelReconstModelCallableArrayMultiVoxelFitReconstFitOdfFitOdfModelReconstFitReconstModelCacheLooseVersionQtdmriFitQtdmriModelCacheExponentialIsotropicFitExponentialIsotropicModelIsotropicFitIsotropicModelReconstFitReconstModelSparseFascicleFitSparseFascicleModelCacheCsaOdfModelOdfFitOdfModelOpdtModelQballBaseModelQballModelResidualBootstrapWrapperSphHarmFitSphHarmModelCacheLooseVersionShoreFitShoreModelreconst
|
Run benchmarks for module using nose. |
|
Run tests for module using nose. |
reconst.baseBase-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 whenver fitting a particular set of data (different voxels, for example).
|
Abstract class which holds the fit result of ReconstModel |
|
Abstract class for signal reconstruction models |
reconst.benchmarks.bench_bounding_boxBenchmarks 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
|
Compute the bounding box of nonzero intensity voxels in the volume. |
|
Return elapsed time for executing code in the namespace of the caller. |
reconst.benchmarks.bench_csd
|
Methods |
|
Diffusion gradient information |
|
|
|
Load only the data array from a nifti file. |
|
|
Read stanford hardi data and label map. |
reconst.benchmarks.bench_peaksBenchmarks 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
Local maxima of a function evaluated on a discrete set of points. |
|
|
Return elapsed time for executing code in the namespace of the caller. |
|
Extract all unique edges from given triangular faces. |
reconst.benchmarks.bench_squashBenchmarks 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
|
Return elapsed time for executing code in the namespace of the caller. |
|
An N-dimensional iterator object to index arrays. |
|
Try and make a standard array from an object array |
Try and make a standard array from an object array |
|
|
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. |
reconst.benchmarks.bench_vec_val_sumBenchmarks 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
|
Return elapsed time for executing code in the namespace of the caller. |
|
Return a sample (or samples) from the “standard normal” distribution. |
Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs |
reconst.cacheCache values based on a key object (such as a sphere or gradient table). |
|
|
Decorator to create OneTimeProperty attributes. |
reconst.cross_validationCross-validation analysis of diffusion models.
|
Calculate the coefficient of determination for a model prediction, |
|
Perform k-fold cross-validation. |
reconst.csdeconv
|
A simple wrapper for response functions represented using only axially symmetric, even spherical harmonic functions (ie, m == 0 and n even). |
|
Methods |
|
Methods |
|
Diffusion data fit to a spherical harmonic model |
|
To be subclassed by all models that return a SphHarmFit when fit. |
|
Diffusion Tensor |
|
Automatic estimation of ssst response function using FA. |
|
Automatic estimation of single-shell single-tissue (ssst) response |
|
Return angles for Cartesian 3D coordinates x, y, and z |
|
Constrained-regularized spherical deconvolution (CSD) [1] |
|
Return decorator function function for deprecation warning / error. |
|
Estimate single fiber response function |
|
|
|
Build forward spherical deconvolution matrix |
|
Build forward sharpening deconvolution transform (SDT) matrix |
|
Return Fractional anisotropy (FA) of a diffusion tensor. |
|
provide triangulated spheres |
|
Produces a lazy index |
|
Legendre function of the first kind. |
|
Computation of mask for single-shell single-tissue (ssst) response |
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
|
An N-dimensional iterator object to index arrays. |
|
ODF constrained-regularized spherical deconvolution using the Sharpening Deconvolution Transform (SDT) [1], [2]. |
|
Sharpen odfs using the sharpening deconvolution transform [2] |
|
Fit the model to data and computes peaks and metrics |
|
Compute a definite integral. |
|
Compute real spherical harmonics. |
|
Samples a real symmetric spherical harmonic basis at point on the sphere |
|
Recursive calibration of response function using peak threshold |
|
Computation of single-shell single-tissue (ssst) response |
|
Computation of single-shell single-tissue (ssst) response |
|
Spherical harmonics (SH) to rotational harmonics (RH) |
|
Simulate diffusion-weighted signals with a single tensor. |
|
Returns the degree ( |
|
rotation matrix from 2 unit vectors |
reconst.dkiClasses and functions for fitting the diffusion kurtosis model
|
Class for fitting the Diffusion Kurtosis Model |
|
Class for the Diffusion Kurtosis Model |
|
Abstract class for signal reconstruction models |
|
|
|
Construct the full 4D kurtosis tensors from its 15 independent elements |
|
Rotate a kurtosis tensor from the standard Cartesian coordinate system to another coordinate system basis |
|
Computes the the specified index element of a kurtosis tensor rotated to the coordinate system basis B. |
|
Calculates the apparent kurtosis coefficient (AKC) in each direction of a sphere [1]. |
|
Computes axial Kurtosis (AK) from the kurtosis tensor [1], [2]. |
|
Computes the Carlson’s incomplete elliptic integral of the second kind defined as: |
|
Computes the Carlson’s incomplete elliptic integral of the first kind defined as: |
|
Return angles for Cartesian 3D coordinates x, y, and z |
|
Check if you have enough different b-values in your gradient table |
|
Returns eigenvalues and eigenvectors given a diffusion tensor |
|
Construct B design matrix for DKI. |
|
Calculates the apparent diffusion coefficient (adc) in each direction of a sphere for a single voxel [1]. |
|
Calculates the apparent diffusion variance (adv) in each direction of a sphere for a single voxel [1]. |
|
Calculates the apparent kurtosis coefficient (akc) in each direction of a sphere for a single voxel [1]. |
|
Predict a signal given diffusion kurtosis imaging parameters. |
Returns a tensor given the six unique tensor elements |
|
|
Provide full paths to example or test datasets. |
|
provide triangulated spheres |
|
Computes the anisotropy of the kurtosis tensor (KFA) [1]. |
|
Computes kurtosis maximum value |
Local maxima of a function evaluated on a discrete set of points. |
|
|
Returns the six lower triangular values of the tensor and a dummy variable if b0 is not None |
|
Mean Diffusivity (MD) of a diffusion tensor. |
|
Computes mean Kurtosis (MK) from the kurtosis tensor. |
|
Computes mean of the kurtosis tensor (MKT) [1]. |
|
An N-dimensional iterator object to index arrays. |
|
Fit the cumulant expansion params (e.g. |
|
Computes the diffusion and kurtosis tensors using an ordinary linear least squares (OLS) approach 1. |
|
Computes n evenly spaced perpendicular directions relative to a given vector v |
|
Radial Kurtosis (RK) of a diffusion kurtosis tensor [1], [2]. |
|
Use the RESTORE algorithm [1] to calculate a robust tensor fit |
|
Spherical to Cartesian coordinates |
|
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 |
Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs |
|
|
Computes the diffusion and kurtosis tensors using a weighted linear least squares (WLS) approach 1. |
reconst.dki_microClasses and functions for fitting the DKI-based microstructural model
|
Class for fitting the Diffusion Kurtosis Model |
|
Class for the Diffusion Kurtosis Model |
|
Class for fitting the Diffusion Kurtosis Microstructural Model |
|
Class for the Diffusion Kurtosis Microstructural Model |
|
Axial Diffusivity (AD) of a diffusion tensor. |
|
Computes the axonal water fraction from DKI [1]. |
|
Returns eigenvalues and eigenvectors given a diffusion tensor |
|
Extracts the restricted and hindered diffusion tensors of well aligned fibers from diffusion kurtosis imaging parameters [1]. |
|
Calculates the apparent diffusion coefficient (adc) in each direction of a sphere for a single voxel [1]. |
|
Calculates the apparent kurtosis coefficient (akc) in each direction of a sphere for a single voxel [1]. |
|
Signal prediction given the DKI microstructure model parameters. |
|
Constructs design matrix for DTI weighted least squares or least squares fitting. |
Returns a tensor given the six unique tensor elements |
|
|
provide triangulated spheres |
|
Computes kurtosis maximum value |
|
Returns the six lower triangular values of the tensor and a dummy variable if b0 is not None |
|
Mean Diffusivity (MD) of a diffusion tensor. |
|
An N-dimensional iterator object to index arrays. |
|
Radial Diffusivity (RD) of a diffusion tensor. |
|
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 |
|
Computes the tortuosity of the hindered diffusion compartment given its axial and radial diffusivities |
|
Trace of a diffusion tensor. |
Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs |
reconst.dsiCache values based on a key object (such as a sphere or gradient table). |
|
|
Methods |
|
Methods |
|
Methods |
|
Methods |
|
Methods |
|
An abstract class to be sub-classed by specific odf models |
|
Perform Lucy-Richardson deconvolution algorithm on a 3D array. |
|
create the 3D grid which holds the signal values (q-space) |
|
create a normalized version of gradients |
|
Return multidimensional discrete Fourier transform. |
|
Shift the zero-frequency component to the center of the spectrum. |
|
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 Cartesian grid mapping |
|
create a hanning window |
|
The inverse of fftshift. |
|
Map the input array to new coordinates by interpolation. |
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
|
Precompute coordinates for ODF calculation from the PDF |
|
Calculates the real ODF from the diffusion propagator(PDF) Pr |
|
Project any near identical bvecs to the other hemisphere |
|
Applies hard threshold on the propagator to remove background noise for the deconvolution. |
reconst.dtiClasses and functions for fitting tensors
|
Abstract class for signal reconstruction models |
|
|
|
Diffusion Tensor |
|
Calculate the apparent diffusion coefficient (ADC) in each direction of a sphere. |
|
Decorator to create OneTimeProperty attributes. |
|
Axial Diffusivity (AD) of a diffusion tensor. |
|
Color fractional anisotropy of diffusion tensor |
|
Returns eigenvalues and eigenvectors given a diffusion tensor |
|
Constructs design matrix for DTI weighted least squares or least squares fitting. |
|
The determinant of a tensor, given in quadratic form |
|
Calculate the deviatoric (anisotropic) part of the tensor [1]. |
|
Calculates tensor eigenvalues/eigenvectors from an array containing the lower diagonal form of the six unique tensor elements. |
|
Return Fractional anisotropy (FA) of a diffusion tensor. |
Returns a tensor given the six unique tensor elements |
|
|
Geodesic anisotropy (GA) of a diffusion tensor. |
|
provide triangulated spheres |
|
A general function for creating diffusion MR gradients. |
|
Calculate the isotropic part of the tensor [1]. |
|
Wrap a fit_tensor func and iterate over chunks of data with given length |
|
The linearity of the tensor [1] |
|
Returns the six lower triangular values of the tensor and a dummy variable if b0 is not None |
|
Mean Diffusivity (MD) of a diffusion tensor. |
|
Mode (MO) of a diffusion tensor [1]. |
|
Fit the cumulant expansion params (e.g. |
|
Calculate the Frobenius norm of a tensor quadratic form |
|
Computes ordinary least squares (OLS) fit to calculate self-diffusion tensor using a linear regression model [Rd310240b4eed-1]. |
|
Vectorized version of numpy.linalg.pinv |
|
The planarity of the tensor [1] |
|
Find the closest orientation of an evenly distributed sphere |
|
Radial Diffusivity (RD) of a diffusion tensor. |
|
Use the RESTORE algorithm [1] to calculate a robust tensor fit |
|
The sphericity of the tensor [1] |
|
Predict a signal given tensor parameters. |
|
Trace of a diffusion tensor. |
Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs |
|
|
Return vector Euclidean (L2) norm |
|
Computes weighted least squares (WLS) fit to calculate self-diffusion tensor using a linear regression model [1]. |
reconst.forecastCache values based on a key object (such as a sphere or gradient table). |
|
|
|
|
Fiber ORientation Estimated using Continuous Axially Symmetric Tensors (FORECAST) [1,2,3]_. |
|
Version numbering for anarchists and software realists. |
|
Methods |
|
An abstract class to be sub-classed by specific odf models |
|
Return angles for Cartesian 3D coordinates x, y, and z |
|
Constrained-regularized spherical deconvolution (CSD) [1] |
|
Calculate the mean signal for each shell. |
|
Calculates the difference between the mean signal calculated using the parameter vector x and the average signal E using FORECAST and SMT |
|
Compute the FORECAST radial matrix |
|
Returns the Laplace-Beltrami regularization matrix for FORECAST |
|
Minimize the sum of squares of a set of equations. |
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
|
Return package-like thing and module setup for package name |
|
|
|
Compute real spherical harmonics. |
|
Compute the SH matrix \(\rho\) |
|
Issue a warning, or maybe ignore it or raise an exception. |
reconst.fwdtiClasses and functions for fitting tensors without free water contamination
|
Class for fitting the Free Water Tensor Model |
|
Class for the Free Water Elimination Diffusion Tensor Model |
|
Abstract class for signal reconstruction models |
|
|
|
Check if you have enough different b-values in your gradient table |
Convert Cholesky decompostion elements to the diffusion tensor elements |
|
|
Returns eigenvalues and eigenvectors given a diffusion tensor |
|
Constructs design matrix for DTI weighted least squares or least squares fitting. |
Returns a tensor given the six unique tensor elements |
|
|
Signal prediction given the free water DTI model parameters. |
|
Returns the six lower triangular values of the tensor and a dummy variable if b0 is not None |
|
Performs Cholesky decomposition of the diffusion tensor |
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
|
An N-dimensional iterator object to index arrays. |
|
Fit the water elimination tensor model using the non-linear least-squares. |
|
Applies non linear least squares fit of the water free elimination model to single voxel signals. |
Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs |
|
|
Computes weighted least squares (WLS) fit to calculate self-diffusion tensor using a linear regression model [1]. |
|
Applies weighted linear least squares fit of the water free elimination model to single voxel signals. |
reconst.gqiClasses and functions for generalized q-sampling
Cache values based on a key object (such as a sphere or gradient table). |
|
|
Methods |
|
Methods |
|
Methods |
|
An abstract class to be sub-classed by specific odf models |
|
|
|
finds the ‘vertices’ in the equatorial zone conjugate to ‘pole’ with width half ‘width’ degrees |
|
The general fractional anisotropy of a function evaluated on the unit sphere |
Local maxima of a function evaluated on a discrete set of points. |
|
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
|
Normalize quantitative anisotropy. |
|
non-parametric anisotropy |
|
|
|
|
|
|
|
find ‘vertices’ within the cone of ‘width’ degrees around ‘pole’ |
|
finds the ‘vertices’ in the equatorial band around the ‘pole’ of radius ‘width’ degrees |
Remove vertices that are less than theta degrees from any other |
|
|
Part of the GQI2 integral |
|
|
maps a 3-vector into the z-upper hemisphere |
reconst.ivimClasses and functions for fitting ivim model
|
|
|
Ivim model |
|
Methods |
|
Version numbering for anarchists and software realists. |
|
Abstract class for signal reconstruction models |
|
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. |
|
Finds the global minimum of a multivariate function. |
|
Error function used to fit f and D_star keeping S0 and D fixed |
|
Function used to predict IVIM signal when S0 and D are known by considering f and D_star as the unknown parameters. |
|
Selector function to switch between the 2-stage Trust-Region Reflective based NLLS fitting method (also containing the linear fit): trr and the Variable Projections based fitting method: varpro. |
|
The Intravoxel incoherent motion (IVIM) model function. |
|
Solve a nonlinear least-squares problem with bounds on the variables. |
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
|
Return package-like thing and module setup for package name |
reconst.mapmriCache values based on a key object (such as a sphere or gradient table). |
|
|
Version numbering for anarchists and software realists. |
|
|
|
Mean Apparent Propagator MRI (MAPMRI) [1] of the diffusion signal. |
|
|
|
Abstract class which holds the fit result of ReconstModel |
|
Abstract class for signal reconstruction models |
|
Calculates the B coefficients from [1] Eq. |
|
Calculates the isotropic B coefficients from [1] Fig 8. |
|
Custom Binomial function |
|
Return angles for Cartesian 3D coordinates x, y, and z |
|
Create the real space table, that contains the points in which to compute the pdf. |
|
|
|
Double factorial. |
|
The GCV cost function that is iterated [4]. |
|
Generalized Cross Validation Function [Rb690cd738504-1] eq. |
|
Generalized Cross Validation Function 1 eq. |
|
Generalized (associated) Laguerre polynomial. |
|
A general function for creating diffusion MR gradients. |
|
Physicist’s Hermite polynomial. |
|
Estimated isotropic scaling factor _[1] Eq. |
|
R(m,n) static matrix for Laplacian regularization [R932dd40ca52e-1] eq. |
|
L(m, n) static matrix for Laplacian regularization [Reb78d789d6c4-1] eq. |
|
S(n, m) static matrix for Laplacian regularization [Rb93dd9dab8c9-1] eq. |
|
Generate the static portions of the Laplacian regularization matrix according to [R1d585103467a-1] eq. |
|
Calculates the indices for the MAPMRI [1] basis in x, y and z. |
Computes mu dependent part of M. |
|
Computes mu independent part of K. |
|
Computed the mu dependent part of the signal design matrix. |
|
Computed the mu independent part of the signal design matrix. |
|
|
Calculates the indices for the isotropic MAPMRI basis [1] Fig 8. |
Computes the Laplacian regularization matrix for MAP-MRI’s isotropic implementation [R156f27ca005f-1] eq. |
|
Computes the Laplacian regularization matrix for MAP-MRI’s isotropic implementation [Rdcc29394f577-1] eq. |
|
|
Compute the isotropic MAPMRI ODF matrix [1] Eq. |
|
Compute the isotropic MAPMRI ODF matrix [1] Eq. |
|
Three dimensional isotropic MAPMRI signal basis function from [1] Eq. |
|
Three dimensional isotropic MAPMRI propagator basis function from [1] Eq. |
|
Radial part of the isotropic 1D-SHORE propagator basis [1] eq. |
|
Radial part of the isotropic 1D-SHORE signal basis [1] eq. |
|
Put the Laplacian regularization matrix together [Rc66aaccd07c1-1] eq. |
|
Compute the MAPMRI ODF matrix [1] Eq. |
|
One dimensional MAPMRI basis function from [1] Eq. |
|
Compute the MAPMRI phi matrix for the signal [1] eq. |
|
One dimensional MAPMRI propagator basis function from [1] Eq. |
|
Compute the MAPMRI psi matrix for the propagator [1] eq. |
|
Find x!. |
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
|
Return package-like thing and module setup for package name |
|
Compute real spherical harmonics. |
|
The factorial of a number or array of numbers. |
|
Returns the degree ( |
|
Issue a warning, or maybe ignore it or raise an exception. |
reconst.mcsd
|
Diffusion gradient information |
|
Version numbering for anarchists and software realists. |
|
|
|
Methods |
|
|
|
Methods |
|
Diffusion Tensor |
|
Automatic estimation of multi-shell multi-tissue (msmt) response |
|
Return Fractional anisotropy (FA) of a diffusion tensor. |
|
Get indices where the b-value is bval |
|
A general function for creating diffusion MR gradients. |
|
Computation of masks for multi-shell multi-tissue (msmt) response |
|
Mean Diffusivity (MD) of a diffusion tensor. |
|
Fiber response function estimation for multi-shell data. |
|
Builds a basis for multi-shell multi-tissue CSD model. |
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
|
Return package-like thing and module setup for package name |
|
Computation of multi-shell multi-tissue (msmt) response |
|
Computation of single-shell single-tissue (ssst) response |
|
Simulate diffusion-weighted signals with a single tensor. |
|
Helper function to set up and solve the Quadratic Program (QP) in CVXPY. |
|
Gives the unique b-values of the data, within a tolerance gap |
reconst.msdkiClasses and functions for fitting the mean signal diffusion kurtosis model
|
|
|
Mean signal Diffusion Kurtosis Model |
|
Abstract class for signal reconstruction models |
|
Decorator to create OneTimeProperty attributes. |
|
Computes the axonal water fraction from the mean signal kurtosis assuming the 2-compartmental spherical mean technique model [1], [2] |
|
Check if you have enough different b-values in your gradient table |
|
Constructs design matrix for the mean signal diffusion kurtosis model |
|
Computes the average signal across different diffusion directions for each unique b-value |
|
Predict the mean signal given the parameters of the mean signal DKI, an GradientTable object and S0 signal. |
|
Computes mean signal kurtosis from axonal water fraction estimates of the SMT2 model |
|
An N-dimensional iterator object to index arrays. |
|
“This function rounds the b-values |
|
This function gives the unique rounded b-values of the data |
|
Fits the mean signal diffusion kurtosis imaging based on a weighted least square solution [1]. |
reconst.multi_voxelTools to easily make multi voxel models
An array which can be called like a function |
|
|
Holds an array of fits and allows access to their attributes and methods |
|
Abstract class which holds the fit result of ReconstModel |
|
Create a view into the array with the given shape and strides. |
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
|
An N-dimensional iterator object to index arrays. |
reconst.odf
|
Methods |
|
An abstract class to be sub-classed by specific odf models |
|
Abstract class which holds the fit result of ReconstModel |
|
Abstract class for signal reconstruction models |
|
The general fractional anisotropy of a function evaluated on the unit sphere |
|
Min-max normalization of a function evaluated on the unit sphere |
reconst.peak_direction_getter
|
|
|
Issue a warning, or maybe ignore it or raise an exception. |
reconst.qtdmriCache values based on a key object (such as a sphere or gradient table). |
|
|
Version numbering for anarchists and software realists. |
|
Methods |
|
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\). |
|
Generalized Cross Validation Function that is iterated [1]. |
|
Step function of H(x)=1 if x>=0 and zero otherwise. |
|
|
|
Angular basis independent of spatial scaling factor us. |
|
Return angles for Cartesian 3D coordinates x, y, and z |
|
Generates EAP grid (for potential positivity constraint). |
|
Constructs design matrix for DTI weighted least squares or least squares fitting. |
|
cross-validation function to find the optimal weight of alpha for sparsity regularization when also Laplacian regularization is used. |
|
The factorial of a number or array of numbers. |
|
Double factorial. |
|
Minimize a function func using the L-BFGS-B algorithm. |
|
Generalized Cross Validation Function [1]. |
|
Generalized (associated) Laguerre polynomial. |
A general function for creating diffusion MR gradients. |
|
|
cross-validation function to find the optimal weight of alpha for sparsity regularization |
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
|
Return package-like thing and module setup for package name |
|
Partial temporal Laplacian regularization matrix following Appendix B in [1]. |
|
Partial spherical spatial Laplacian regularization matrix following the equation below Eq. |
|
Partial cartesian spatial Laplacian regularization matrix following second line of Eq. |
|
Partial temporal Laplacian regularization matrix following Appendix B in [1]. |
|
Partial spherical spatial Laplacian regularization matrix following the equation below Eq. |
|
Partial cartesian spatial Laplacian regularization matrix following equation Eq. |
|
Partial temporal Laplacian regularization matrix following Appendix B in [1]. |
|
Constructs design matrix for fitting an exponential to the diffusion time points. |
|
Constructs the design matrix as a product of 3 separated radial, angular and temporal design matrices. |
|
|
|
Computes the SHORE basis order indices according to [1]. |
|
Constructs the design matrix as a product of 3 separated radial, angular and temporal design matrices. |
|
|
|
Computes the SHORE basis order indices according to [1]. |
Computes the spherical qt-dMRI Laplacian regularization matrix. |
|
|
Constructs design matrix for fitting an exponential to the diffusion time points. |
|
|
|
|
Generates the matrix that maps the spherical qtdmri coefficients to MAP-MRI coefficients. |
|
|
Computes the cartesian qt-dMRI Laplacian regularization matrix. |
|
Normalization factor for Spherical MAP-MRI basis. |
Normalization factor for Cartesian MAP-MRI basis. |
|
|
Computes the total number of coefficients of the qtdmri basis given a radial and temporal order. |
|
Constructs the design matrix as a product of 3 separated radial, angular and temporal design matrices. |
|
Function to generate the qtdmri signal basis. |
Normalization factor for the temporal basis |
|
|
Generates the matrix that maps the qtdmri coefficients to MAP-MRI coefficients. |
|
|
|
Spatial basis dependent on spatial scaling factor us |
|
Compute real spherical harmonics. |
|
Temporal basis dependent on temporal scaling factor ut |
This function visualizes a q-tau acquisition scheme as a function of gradient strength and pulse separation (big_delta). |
|
|
Issue a warning, or maybe ignore it or raise an exception. |
reconst.sfmThe 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].
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
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 values based on a key object (such as a sphere or gradient table). |
|
|
A fit to the ExponentialIsotropicModel object, based on data. |
Representing the isotropic signal as a fit to an exponential decay function with b-values |
|
|
A fit object for representing the isotropic signal as the mean of the diffusion-weighted signal. |
|
A base-class for the representation of isotropic signals. |
|
Abstract class which holds the fit result of ReconstModel |
|
Abstract class for signal reconstruction models |
|
Methods |
|
Methods |
|
Decorator to create OneTimeProperty attributes. |
|
Compute the arithmetic mean along the specified axis, ignoring NaNs. |
|
Return package-like thing and module setup for package name |
|
Construct the SFM design matrix |
reconst.shmTools for using spherical harmonic models to fit diffusion data.
Angle Consideration.
Q-ball imaging.
populations in white matter of the brain with Funk-Radon transform.
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 values based on a key object (such as a sphere or gradient table). |
|
|
Implementation of Constant Solid Angle reconstruction method. |
|
Methods |
|
An abstract class to be sub-classed by specific odf models |
|
Implementation of Orientation Probability Density Transform reconstruction method. |
|
To be subclassed by Qball type models. |
|
Implementation of regularized Qball reconstruction method. |
|
Returns a residual bootstrap sample of the signal_object when indexed |
|
Diffusion data fit to a spherical harmonic model |
|
To be subclassed by all models that return a SphHarmFit when fit. |
|
Calculate anisotropic power map with a given SH coefficient matrix. |
|
Decorator to create OneTimeProperty attributes. |
|
Applies the Residual Bootstraps to the data given H and R |
|
Like bootstrap_data_array but faster when for a single voxel |
|
Calculate the maximal harmonic order, given that you know the number of parameters that were estimated. |
|
Return angles for Cartesian 3D coordinates x, y, and z |
|
Convert SH coefficients in legacy SH basis to SH coefficients of the new SH basis for |
|
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 coefficients in new SH basis to SH coefficients for the legacy SH basis for |
|
Return decorator function function for deprecation warning / error. |
|
Build forward spherical deconvolution matrix |
|
Generate Dirac delta function orientated in (theta, phi) on the sphere |
|
Returns the hat matrix for the design matrix B |
|
Produces a lazy index |
|
Returns a matrix for computing leveraged, centered residuals from data |
|
Normalizes the data with respect to the mean b0 |
|
Given a number |
|
Return random integers from low (inclusive) to high (exclusive). |
|
Compute real spherical harmonics as in Descoteaux et al. |
|
Compute real spherical harmonics as in Descoteaux et al. |
|
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: |
|
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: |
|
Compute real spherical harmonics. |
|
Samples a real symmetric spherical harmonic basis at point on the sphere |
|
dipy.reconst.shm.real_sym_sh_mrtix is deprecated, Please use dipy.reconst.shm.real_sh_tournier instead |
|
Spherical function to spherical harmonics (SH). |
|
Spherical harmonics (SH) to rotational harmonics (RH) |
|
Spherical harmonics (SH) to spherical function (SF). |
|
Matrix that transforms Spherical harmonics (SH) to spherical function (SF). |
|
Regularized pseudo-inverse |
|
Returns the degree ( |
|
Compute spherical harmonics. |
|
Issue a warning, or maybe ignore it or raise an exception. |
reconst.shoreCache values based on a key object (such as a sphere or gradient table). |
|
|
Version numbering for anarchists and software realists. |
|
|
|
Simple Harmonic Oscillator based Reconstruction and Estimation (SHORE) [1] of the diffusion signal. |
|
Return angles for Cartesian 3D coordinates x, y, and z |
|
Create the real space table, that contains the points in which |
|
Find x!. |
|
Generalized (associated) Laguerre polynomial. |
|
Returns the angular regularisation matrix for SHORE basis |
|
Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition |
|
Returns the angular regularisation matrix for SHORE basis |
|
Return package-like thing and module setup for package name |
|
Compute real spherical harmonics. |
|
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} |
|
Compute the SHORE matrix for modified Merlet’s 3D-SHORE [1] |
|
Compute the SHORE ODF matrix [1]” |
|
Compute the SHORE propagator matrix [1]” |
|
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. |
|
Issue a warning, or maybe ignore it or raise an exception. |
reconst.utils
|
Construct B design matrix for DKI. |
dipy.reconst.bench(label='fast', verbose=1, extra_argv=None)Run benchmarks for module using nose.
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’.
Verbosity value for benchmark outputs, in the range 1-10. Default is 1.
List with any extra arguments to pass to nosetests.
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
>>> success
True
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.
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’.
Verbosity value for test outputs, in the range 1-10. Default is 1.
List with any extra arguments to pass to nosetests.
If True, run doctests in module. Default is False.
If True, report coverage of NumPy code. Default is False. (This requires the coverage module).
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.
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 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:
>>> np.lib.test()
Examples
>>> result = np.lib.test()
Running unit tests for numpy.lib
...
Ran 976 tests in 3.933s
OK
>>> result.errors
[]
>>> result.knownfail
[]
ReconstFitdipy.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.
ReconstModeldipy.reconst.base.ReconstModel(gtab)Bases: object
Abstract class for signal reconstruction models
Methods
fit |
dipy.reconst.benchmarks.bench_bounding_box.bounding_box(vol)Compute the bounding box of nonzero intensity voxels in the volume.
Volume to compute bounding box on.
Array containg minimum index of each dimension
Array containg maximum index of each dimension
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.
The code to be timed.
The number of times the code is executed. Default is 1. The code is only compiled once.
A label to identify code_str with. This is passed into compile
as the second argument (for run-time error messages).
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
ConstrainedSphericalDeconvModeldipy.reconst.benchmarks.bench_csd.ConstrainedSphericalDeconvModel(gtab, response, reg_sphere=None, sh_order=8, lambda_=1, tau=0.1, convergence=50)Bases: dipy.reconst.shm.SphHarmModel
Methods
|
Clear the cache. |
|
Retrieve a value from the cache. |
|
Store a value in the cache. |
|
Fit method for every voxel in data |
|
Compute a signal prediction given spherical harmonic coefficients for the provided GradientTable class instance. |
|
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].
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]).
sphere used to build the regularization B matrix. Default: ‘symmetric362’.
maximal spherical harmonics order. Default: 8
weight given to the constrained-positivity regularization part of the deconvolution equation (see [1]). Default: 1
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
Maximum number of iterations to allow the deconvolution to converge.
References
Tournier, J.D., et al. NeuroImage 2007. Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution
Descoteaux, M., et al. IEEE TMI 2009. Deterministic and Probabilistic Tractography Based on Complex Fibre Orientation Distributions
Côté, M-A., et al. Medical Image Analysis 2013. Tractometer: Towards validation of tractography pipelines
Tournier, J.D, et al. Imaging Systems and Technology 2012. MRtrix: Diffusion Tractography in Crossing Fiber Regions
predict(sh_coeff, gtab=None, S0=1.0)Compute a signal prediction given spherical harmonic coefficients for the provided GradientTable class instance.
The spherical harmonic representation of the FOD from which to make the signal prediction.
The gradients for which the signal will be predicted. Use the model’s gradient table by default.
The non diffusion-weighted signal value.
The predicted signal.
GradientTabledipy.reconst.benchmarks.bench_csd.GradientTable(gradients, big_delta=None, small_delta=None, b0_threshold=50, btens=None)Bases: object
Diffusion gradient information
Diffusion gradients. The direction of each of these vectors corresponds to the b-vector, and the length corresponds to the b-value.
Gradients with b-value less than or equal to b0_threshold are considered as b0s i.e. without diffusion weighting.
See also
gradient_tableNotes
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
diffusion gradients
The b-value, or magnitude, of each gradient direction.
The q-value for each gradient direction. Needs big and small delta.
The direction, represented as a unit vector, of each gradient.
Boolean array indicating which gradients have no diffusion weighting, ie b-value is close to 0.
Gradients with b-value less than or equal to b0_threshold are considered to not have diffusion weighting.
The b-tensor of each gradient direction.
Methods
b0s_mask |
|
bvals |
|
bvecs |
|
gradient_strength |
|
qvals |
|
tau |
dipy.reconst.benchmarks.bench_csd.load_nifti_data(fname, as_ndarray=True)Load only the data array from a nifti file.
Full path to the file.
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)
See also
load_niftidipy.reconst.benchmarks.bench_peaks.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.
The function evaluated on a set of discrete points.
The set of neighbor relations between the points. Every edge, ie edges[i, :], is a pair of neighboring points.
Value of odf at a maximum point. Peak values is sorted in descending order.
Indices of maximum points. Sorted in the same order as peak_values so odf[peak_indices[i]] == peak_values[i].
See also
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).
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.
The code to be timed.
The number of times the code is executed. Default is 1. The code is only compiled once.
A label to identify code_str with. This is passed into compile
as the second argument (for run-time error messages).
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
dipy.reconst.benchmarks.bench_peaks.unique_edges(faces, return_mapping=False)Extract all unique edges from given triangular faces.
Vertex indices forming triangular faces.
If true, a mapping to the edges of each face is returned.
Unique edges.
For each face, [x, y, z], a mapping to it’s edges [a, b, c].
y
/ / a/
/ / /__________ x c z
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.
The code to be timed.
The number of times the code is executed. Default is 1. The code is only compiled once.
A label to identify code_str with. This is passed into compile
as the second argument (for run-time error messages).
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
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.
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)
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.
The array to be converted.
Where arr has Nones.
Nones are replaced by fill.
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')
dipy.reconst.benchmarks.bench_squash.quick_squash()Try and make a standard array from an object array
This function takes an object array and attempts to convert it to a more useful dtype. If array can be converted to a better dtype, Nones are replaced by fill. To make the behaviour of this function more clear, here are the most common cases:
obj_arr is an array of scalars of type T. Returns an array like obj_arr.astype(T)
obj_arr is an array of arrays. All items in obj_arr have the same
shape S. Returns an array with shape obj_arr.shape + S
obj_arr is an array of arrays of different shapes. Returns obj_arr.
Items in obj_arr are not ndarrays or scalars. Returns obj_arr.
The array to be converted.
mask is nonzero where obj_arr has Nones.
Nones are replaced by fill.
Examples
>>> arr = np.empty(3, dtype=object)
>>> arr.fill(2)
>>> quick_squash(arr)
array([2, 2, 2])
>>> arr[0] = None
>>> quick_squash(arr)
array([0, 2, 2])
>>> arr.fill(np.ones(2))
>>> r = quick_squash(arr)
>>> r.shape
(3, 2)
>>> r.dtype
dtype('float64')
dipy.reconst.benchmarks.bench_squash.reduce(function, sequence[, initial]) → valueApply 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.
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.
The code to be timed.
The number of times the code is executed. Default is 1. The code is only compiled once.
A label to identify code_str with. This is passed into compile
as the second argument (for run-time error messages).
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
dipy.reconst.benchmarks.bench_vec_val_sum.randn(d0, d1, ..., dn)Return a sample (or samples) from the “standard normal” distribution.
Note
This is a convenience function for users porting code from Matlab, and wraps standard_normal. That function takes a tuple to specify the size of the output, which is consistent with other NumPy functions like numpy.zeros and numpy.ones.
Note
New code should use the standard_normal method of a default_rng()
instance instead; please see the random-quick-start.
If positive int_like arguments are provided, randn generates an array
of shape (d0, d1, ..., dn), filled
with random floats sampled from a univariate “normal” (Gaussian)
distribution of mean 0 and variance 1. A single float randomly sampled
from the distribution is returned if no argument is provided.
The dimensions of the returned array, must be non-negative. If no argument is given a single Python float is returned.
A (d0, d1, ..., dn)-shaped array of floating-point samples from
the standard normal distribution, or a single such float if
no parameters were supplied.
See also
standard_normalSimilar, but takes a tuple as its argument.
normalAlso accepts mu and sigma arguments.
Generator.standard_normalwhich should be used for new code.
Notes
For random samples from \(N(\mu, \sigma^2)\), use:
sigma * np.random.randn(...) + mu
Examples
>>> np.random.randn()
2.1923875335537315 # random
Two-by-four array of samples from N(3, 6.25):
>>> 3 + 2.5 * np.random.randn(2, 4)
array([[-4.49401501, 4.00950034, -1.81814867, 7.29718677], # random
[ 0.39924804, 4.68456316, 4.99394529, 4.84057254]]) # random
dipy.reconst.benchmarks.bench_vec_val_sum.vec_val_vect()Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs
containing tensor in last two dimensions; M, N usually equal to (3, 3)
diagonal values carried in last dimension, ... shape above must
match that for vecs
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).
... dimensions of vecs, valsN dimensions of vecs, valsExamples
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]])
Cachedipy.reconst.cache.CacheBases: 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
Clear the cache. |
|
|
Retrieve a value from the cache. |
|
Store a value in the cache. |
cache_get(tag, key, default=None)Retrieve a value from the cache.
Description of the cached value.
Key object used to look up the cached value.
Value to be returned if no cached entry is found.
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.
Description of the cached value.
Key object used to look up the cached value.
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)
>>> c = Cache()
>>> parameters = (1, 2, 3)
>>> X1 = compute_expensive_matrix(parameters)
>>> c.cache_set('expensive_matrix', parameters, X1)
>>> X2 = c.cache_get('expensive_matrix', parameters)
>>> X1 is X2
True
dipy.reconst.cache.auto_attr(func)Decorator to create OneTimeProperty attributes.
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
dipy.reconst.cross_validation.coeff_of_determination(data, model, axis=-1)relative to data.
The data
The predictions of a model for this data. Same shape as the data.
The axis along which different samples are laid out (default: -1).
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).
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.
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.
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.
The number of divisions to apply to the data
Additional arguments to the model initialization
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
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.
AxSymShResponsedipy.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).
Signal with no diffusion weighting.
Response function signal as coefficients to axially symmetric, even spherical harmonic.
Methods
|
A basis that maps the response coefficients onto a sphere. |
|
Evaluates the response function on sphere. |
ConstrainedSDTModeldipy.reconst.csdeconv.ConstrainedSDTModel(gtab, ratio, reg_sphere=None, sh_order=8, lambda_=1.0, tau=0.1)Bases: dipy.reconst.shm.SphHarmModel
Methods
|
Clear the cache. |
|
Retrieve a value from the cache. |
|
Store a value in the cache. |
|
Fit method for every voxel in data |
|
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.
ratio of the smallest vs the largest eigenvalue of the single prolate tensor response function
sphere used to build the regularization B matrix
maximal spherical harmonics order
weight given to the constrained-positivity regularization part of the deconvolution equation
threshold (tau *mean(fODF)) controlling the amplitude below which the corresponding fODF is assumed to be zero.
References
ConstrainedSphericalDeconvModeldipy.reconst.csdeconv.ConstrainedSphericalDeconvModel(gtab, response, reg_sphere=None, sh_order=8, lambda_=1, tau=0.1, convergence=50)Bases: dipy.reconst.shm.SphHarmModel
Methods
|
Clear the cache. |
|
Retrieve a value from the cache. |
|
Store a value in the cache. |
|
Fit method for every voxel in data |
|
Compute a signal prediction given spherical harmonic coefficients for the provided GradientTable class instance. |
|
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].
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]).
sphere used to build the regularization B matrix. Default: ‘symmetric362’.
maximal spherical harmonics order. Default: 8
weight given to the constrained-positivity regularization part of the deconvolution equation (see [1]). Default: 1
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
Maximum number of iterations to allow the deconvolution to converge.
References
Tournier, J.D., et al. NeuroImage 2007. Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution
Descoteaux, M., et al. IEEE TMI 2009. Deterministic and Probabilistic Tractography Based on Complex Fibre Orientation Distributions
Côté, M-A., et al. Medical Image Analysis 2013. Tractometer: Towards validation of tractography pipelines
Tournier, J.D, et al. Imaging Systems and Technology 2012. MRtrix: Diffusion Tractography in Crossing Fiber Regions
predict(sh_coeff, gtab=None, S0=1.0)Compute a signal prediction given spherical harmonic coefficients for the provided GradientTable class instance.
The spherical harmonic representation of the FOD from which to make the signal prediction.
The gradients for which the signal will be predicted. Use the model’s gradient table by default.
The non diffusion-weighted signal value.
The predicted signal.
SphHarmFitdipy.reconst.csdeconv.SphHarmFit(model, shm_coef, mask)Bases: dipy.reconst.odf.OdfFit
Diffusion data fit to a spherical harmonic model
shm_coeffThe spherical harmonic coefficients of the odf
Methods
|
Samples the odf function on the points of a sphere |
|
Predict the diffusion signal from the model coefficients. |
gfa |
odf(sphere)Samples the odf function on the points of a sphere
The points on which to sample the odf.
The value of the odf on each point of sphere.
SphHarmModeldipy.reconst.csdeconv.SphHarmModel(gtab)Bases: dipy.reconst.odf.OdfModel, dipy.reconst.cache.Cache
To be subclassed by all models that return a SphHarmFit when fit.
Methods
|
Clear the cache. |
|
Retrieve a value from the cache. |
|
Store a value in the cache. |
|
To be implemented by specific odf models |
|
The matrix needed to sample ODFs from coefficients of the model. |
__init__(gtab)Initialization of the abstract class for signal reconstruction models
sampling_matrix(sphere)The matrix needed to sample ODFs from coefficients of the model.
Points used to sample ODF.
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.
TensorModeldipy.reconst.csdeconv.TensorModel(gtab, fit_method='WLS', return_S0_hat=False, *args, **kwargs)Bases: dipy.reconst.base.ReconstModel
Diffusion Tensor
Methods
|
Fit method of the DTI model class |
|
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].
str can be one of the following:
dti.wls_fit_tensor()
dti.ols_fit_tensor()
dti.nlls_fit_tensor()
fitting [3]
dti.restore_fit_tensor()
Boolean to return (True) or not (False) the S0 values for the fit.
fit_method. See dti.wls_fit_tensor, dti.ols_fit_tensor for details
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
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.
Basser, P., Pierpaoli, C., 1996. Microstructural and physiological features of tissues elucidated by quantitative diffusion-tensor MRI. Journal of Magnetic Resonance 111, 209-219.
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
The measured signal from one voxel.
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.
The last dimension should have 12 tensor parameters: 3 eigenvalues, followed by the 3 eigenvectors
The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1
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
deprecated from version: 1.2
Will raise <class ‘dipy.utils.deprecator.ExpiredDeprecationError’> as of version: 1.4
diffusion data
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].
radius of cubic ROI
FA threshold
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.
If True, returns the number of voxels used for estimating the response function
(evals, S0)
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).
dipy.reconst.csdeconv.auto_response_ssst(gtab, data, roi_center=None, roi_radii=10, fa_thr=0.7)function using FA.
diffusion data
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].
radii of cuboid ROI
FA threshold
(evals, S0)
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).
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\)
x coordinate in Cartesian space
y coordinate in Cartesian space
z coordinate
radius
inclination (polar) angle
azimuth angle
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.
Diffusion weighted signals to be deconvolved.
Prediction matrix which estimates diffusion weighted signals from FOD coefficients.
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.
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.
Maximum number of iterations to allow the deconvolution to converge.
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.
(sh_order + 1)*(sh_order + 2)/2,)Spherical harmonics coefficients of the constrained-regularized fiber ODF.
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
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).
Message explaining deprecation, giving possible alternatives.
Released version at which object was first deprecated.
Last released version at which this function will still raise a deprecation warning. Versions higher than this will raise an error.
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.
Class of warning to generate for deprecation.
Class of error to generate when version_comparator returns 1 for a
given argument of until.
Function returning a decorator.
dipy.reconst.csdeconv.forward_sdeconv_mat(r_rh, n)Build forward spherical deconvolution matrix
Rotational harmonics coefficients for the single fiber response
function. Each element rh[i] is associated with spherical harmonics
of degree 2*i.
The order of spherical harmonic function associated with each row of the deconvolution matrix. Only even orders are allowed
Deconvolution matrix with shape (N, N)
dipy.reconst.csdeconv.forward_sdt_deconv_mat(ratio, n, r2_term=False)Build forward sharpening deconvolution transform (SDT) matrix
ratio = \(\frac{\lambda_2}{\lambda_1}\) of the single fiber response function
The degree of spherical harmonic function associated with each row of the deconvolution matrix. Only even degrees are allowed.
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].
SDT deconvolution matrix
Funk-Radon Transform (FRT) matrix
References
Descoteaux, M. PhD Thesis. INRIA Sophia-Antipolis. 2008.
dipy.reconst.csdeconv.fractional_anisotropy(evals, axis=-1)Return Fractional anisotropy (FA) of a diffusion tensor.
Eigenvalues of a diffusion tensor.
Axis of evals which contains 3 eigenvalues.
Calculated FA. Range is 0 <= FA <= 1.
Notes
FA is calculated using the following equation:
dipy.reconst.csdeconv.get_sphere(name='symmetric362')provide triangulated spheres
which sphere - one of: * ‘symmetric362’ * ‘symmetric642’ * ‘symmetric724’ * ‘repulsion724’ * ‘repulsion100’ * ‘repulsion200’
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"
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
Zhang, Shanjie and Jin, Jianming. “Computation of Special Functions”, John Wiley and Sons, 1996. https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
dipy.reconst.csdeconv.mask_for_response_ssst(gtab, data, roi_center=None, roi_radii=10, fa_thr=0.7)function using FA.
diffusion data (4D)
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].
radii of cuboid ROI
FA threshold
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
Tournier, J.D., et al. NeuroImage 2004. Direct estimation of the
fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution
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.
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)
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].
(sh_order + 1)*(sh_order + 2)/2,)ndarray of SH coefficients for the ODF spherical function to be deconvolved
(sh_order + 1)(sh_order + 2)/2,(sh_order + 1)(sh_order + 2)/2)
SDT matrix in SH basis
(sh_order + 1)(sh_order + 2)/2,(sh_order + 1)(sh_order + 2)/2)
SH basis matrix used for deconvolution
lambda parameter in minimization equation (default 1.0)
threshold (tau *max(fODF)) controlling the amplitude below which the corresponding fODF is assumed to be zero.
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.
(sh_order + 1)(sh_order + 2)/2,)Spherical harmonics coefficients of the constrained-regularized fiber ODF
Number of iterations in the constrained-regularization used for convergence
References
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.
(sh_order + 1)*(sh_order + 2)/2, )array of odfs expressed as spherical harmonics coefficients
sphere used to build the regularization matrix
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).
ratio of the smallest vs the largest eigenvalue of the single prolate tensor response function (\(\frac{\lambda_2}{\lambda_1}\))
maximal SH order of the SH representation
lambda parameter (see odfdeconv) (default 1.0)
tau parameter in the L matrix construction (see odfdeconv) (default 0.1)
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.
sharpened odf expressed as spherical harmonics coefficients
References
Tuch, D. MRM 2004. Q-Ball Imaging.
Descoteaux, M., et al. IEEE TMI 2009. Deterministic and Probabilistic Tractography Based on Complex Fibre Orientation Distributions
Descoteaux, M., Angelino, E., Fitzgibbons, S. and Deriche, R. Regularized, Fast, and Robust Analytical Q-ball Imaging. Magn. Reson. Med. 2007;58:497-510.
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.
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, nbr_processes=None)Fit the model to data and computes peaks and metrics
model will be used to fit the data.
The Sphere providing discrete directions for evaluation.
Only return peaks greater than relative_peak_threshold * m where m
is the largest peak.
directions. If two peaks are too close only the larger of the two is returned.
If mask is provided, voxels that are False in mask are skipped and no peaks are returned.
If True, the odfs are returned.
If True, the odf as spherical harmonics coefficients is returned
Voxels with gfa less than gfa_thr are skipped, no peaks are returned.
If true, all peak values are calculated relative to max(odf).
Maximum SH order in the SH fit. For sh_order, there will be
(sh_order + 1) * (sh_order + 2) / 2 SH coefficients (default 8).
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).
Lambda-regularization in the SH fit (default 0.0).
Maximum number of peaks found (default 5 peaks).
Matrix that transforms spherical harmonics to spherical function
sf = np.dot(sh, B).
Inverse of B.
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'.
If parallel is True, the number of subprocesses to use (default multiprocessing.cpu_count()).
An object with gfa, peak_directions, peak_values,
peak_indices, odf, shm_coeffs as attributes
References
Descoteaux, M., Angelino, E., Fitzgibbons, S. and Deriche, R. Regularized, Fast, and Robust Analytical Q-ball Imaging. Magn. Reson. Med. 2007;58:497-510.
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.
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)Compute a definite integral.
Integrate func from a to b (possibly infinite interval) using a technique from the Fortran library QUADPACK.
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.
Lower limit of integration (use -numpy.inf for -infinity).
Upper limit of integration (use numpy.inf for +infinity).
Extra arguments to pass to func.
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.
The integral of func from a to b.
An estimate of the absolute error in the result.
A dictionary containing additional information. Run scipy.integrate.quad_explain() for more information.
A convergence message.
Appended only with ‘cos’ or ‘sin’ weighting and infinite integration limits, it contains an explanation of the codes in infodict[‘ierlst’]
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.
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.
An upper bound on the number of subintervals used in the adaptive algorithm.
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.
String indicating weighting function. Full explanation for this and the remaining arguments can be found below.
Variables for use with weighting functions.
Optional input for reusing Chebyshev moments.
An upper bound on the number of Chebyshev moments.
Upper bound on the number of cycles (>=3) for use with a sinusoidal weighting and an infinite end-point.
See also
dblquaddouble integral
tplquadtriple integral
nquadn-dimensional integrals (uses quad recursively)
fixed_quadfixed-order Gaussian quadrature
quadratureadaptive Gaussian quadrature
odeintODE integrator
odeODE integrator
simpsintegrator for sampled data
rombintegrator for sampled data
scipy.specialfor 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:
The number of function evaluations.
The number, K, of subintervals produced in the subdivision process.
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.
A rank-1 array of length M, the first K elements of which are the right end points of the subintervals.
A rank-1 array of length M, the first K elements of which are the integral approximations on the subintervals.
A rank-1 array of length M, the first K elements of which are the moduli of the absolute error estimates on the subintervals.
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.
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.
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).
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 function used |
|
|---|---|---|
‘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:
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.
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).
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’:
The number of subintervals needed for the integration (call it K_f).
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.
A rank-1 array of length M_f containing the error estimate
corresponding to the interval in the same position in
infodict['rslist'].
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.
Examples
Calculate \(\int^4_0 x^2 dx\) and compare with an analytic result
>>> from scipy import integrate
>>> 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)
>>> 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)
dipy.reconst.csdeconv.real_sph_harm(m, n, theta, phi)Compute real spherical harmonics.
dipy.reconst.shm.real_sph_harm is deprecated, Please use dipy.reconst.shm.real_sh_descoteaux_from_index instead
deprecated from version: 1.3
Will raise <class ‘dipy.utils.deprecator.ExpiredDeprecationError’> as of version: 2.0
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 broadcasted against each other.
|m| <= nThe degree of the harmonic.
>= 0The order of the harmonic.
The polar (colatitudinal) coordinate.
The azimuthal (longitudinal) coordinate.
The real harmonic \(Y^m_n\) sampled at theta and phi.
See also
scipy.special.sph_harmdipy.reconst.csdeconv.real_sym_sh_basis(sh_order, theta, phi)Samples a real symmetric spherical harmonic basis at point on the sphere
dipy.reconst.shm.real_sym_sh_basis is deprecated, Please use dipy.reconst.shm.real_sh_descoteaux instead
deprecated from version: 1.3
Will raise <class ‘dipy.utils.deprecator.ExpiredDeprecationError’> as of version: 2.0
Samples the basis functions up to order sh_order at points on the sphere given by theta and phi. The basis functions are defined here the same way as in Descoteaux et al. 2007 [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 broadcasted against each other.
even int > 0, max spherical harmonic order
The azimuthal (longitudinal) coordinate.
The polar (colatitudinal) coordinate.
The real harmonic \(Y^m_n\) sampled at theta and phi
The degree of the harmonics.
The order of the harmonics.
References
Descoteaux, M., Angelino, E., Fitzgibbons, S. and Deriche, R. Regularized, Fast, and Robust Analytical Q-ball Imaging. Magn. Reson. Med. 2007;58:497-510.
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=True, nbr_processes=None, sphere=<dipy.core.sphere.HemiSphere object>)Recursive calibration of response function using peak threshold
diffusion data
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.
maximal spherical harmonics order. Default: 8
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
FA of the initial ‘fat’ response function (tensor). Default: 0.08
trace of the initial ‘fat’ response function (tensor). Default: 0.0021
maximum number of iterations for calibration. Default: 8.
convergence criterion, maximum relative change of SH coefficients. Default: 0.001.
Whether to use parallelization in peak-finding during the calibration procedure. Default: True
If parallel is True, the number of subprocesses to use (default multiprocessing.cpu_count()).
The sphere used for peak finding. Default: default_sphere.
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
Tax, C.M.W., et al. NeuroImage 2014. Recursive calibration of the fiber response function for spherical deconvolution of diffusion MRI data.
dipy.reconst.csdeconv.response_from_mask(gtab, data, mask)function from a given mask.
dipy.reconst.csdeconv.response_from_mask is deprecated, Please use dipy.reconst.csdeconv.response_from_mask_ssst instead
deprecated from version: 1.2
Will raise <class ‘dipy.utils.deprecator.ExpiredDeprecationError’> as of version: 1.4
diffusion data
mask from where to compute the response function
(evals, S0)
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
Tournier, J.D., et al. NeuroImage 2004. Direct estimation of the
fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution
dipy.reconst.csdeconv.response_from_mask_ssst(gtab, data, mask)function from a given mask.
diffusion data
mask from where to compute the response function
(evals, S0)
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
Tournier, J.D., et al. NeuroImage 2004. Direct estimation of the
fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution
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.
ndarray of SH coefficients for the single fiber response function. These coefficients must correspond to the real spherical harmonic functions produced by shm.real_sph_harm.
The degree of the spherical harmonic function associated with each coefficient.
The order of the spherical harmonic function associated with each coefficient.
(sh_order + 1)*(sh_order + 2)/2,)Rotational harmonics coefficients representing the input r_sh
See also
shm.real_sph_harm, shm.real_sym_sh_basisReferences
Tournier, J.D., et al. NeuroImage 2007. Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution
dipy.reconst.csdeconv.single_tensor(gtab, S0=1, evals=None, evecs=None, snr=None)Simulate diffusion-weighted signals with a single tensor.
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.
Strength of signal in the presence of no diffusion gradient (also
called the b=0 value).
Eigenvalues of the diffusion tensor. By default, values typical for prolate white matter are used.
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.
Signal to noise ratio, assuming Rician noise. None implies no noise.
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
M. Descoteaux, “High Angular Resolution Diffusion MRI: from Local Estimation to Segmentation and Tractography”, PhD thesis, University of Nice-Sophia Antipolis, p. 42, 2008.
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.
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.
even int > 0, max order to return
True for SH basis with even and odd order terms
degrees of even spherical harmonics
orders of even spherical harmonics
See also
shm.real_sh_descoteaux_from_index, shm.real_sh_tournier_from_indexdipy.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.
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.])
DiffusionKurtosisFitdipy.reconst.dki.DiffusionKurtosisFit(model, model_params)Bases: dipy.reconst.dti.TensorFit
Class for fitting the Diffusion Kurtosis Model
directionsFor tracking - return the primary direction in each voxel
evalsReturns the eigenvalues of the tensor as an array
evecsReturns the eigenvectors of the tensor as an array, columnwise
kfaReturns the kurtosis tensor (KFA) 1.
ktReturns the 15 independent elements of the kurtosis tensor as an array
quadratic_formCalculates the 3x3 diffusion tensor for each voxel
Methods
|
Axial diffusivity (AD) calculated from cached eigenvalues. |
|
Calculate the apparent diffusion coefficient (ADC) in each direction on |
|
Axial Kurtosis (AK) of a diffusion kurtosis tensor [1]. |
|
Calculates the apparent kurtosis coefficient (AKC) in each direction on the sphere for each voxel in the data |
|
Color fractional anisotropy of diffusion tensor |
|
Fractional anisotropy (FA) calculated from cached eigenvalues. |
|
Geodesic anisotropy (GA) calculated from cached eigenvalues. |
|
Computes the maximum value of a single voxel kurtosis tensor |
|
|
|
Mean diffusivity (MD) calculated from cached eigenvalues. |
|
Computes mean Kurtosis (MK) from the kurtosis tensor. |
|
Computes mean of the kurtosis tensor (MKT) [1]. |
|
Tensor mode calculated from cached eigenvalues. |
|
The diffusion orientation distribution function (dODF). |
|
|
|
Given a DKI model fit, predict the signal on the vertices of a gradient table |
|
Radial diffusivity (RD) calculated from cached eigenvalues. |
|
Radial Kurtosis (RK) of a diffusion kurtosis tensor [1]. |
|
|
|
Trace of the tensor calculated from cached eigenvalues. |
lower_triangular |
__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
Class instance containing the Diffusion Kurtosis Model for the fit
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].
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])
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
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.
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 can be sampled from the principal axis of the diffusion tensor:
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
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of non-Gaussian water diffusion by kurtosis analysis. NMR in Biomedicine 23(7): 698-710
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
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
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:
where \(W_{ijkl}\) are the elements of the kurtosis tensor, MD the mean diffusivity and ADC the apparent diffusion coefficent computed as:
where \(D_{ij}\) are the elements of the diffusion tensor.
kfaReturns the kurtosis tensor (KFA) [1].
Notes
The KFA is defined as [1]:
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
kmax(sphere='repulsion100', gtol=1e-05, mask=None)Computes the maximum value of a single voxel kurtosis tensor
The sphere providing sample directions for the initial search of the maximum value of kurtosis.
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
kurtosis tensor maximum value
mk(min_kurtosis=-0.42857142857142855, max_kurtosis=10, analytical=True)Computes mean Kurtosis (MK) from the kurtosis tensor.
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])
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
If True, MK is calculated using its analytical solution, otherwise an exact numerical estimator is used (see Notes). Default is set to True.
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]_:
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:
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
where \(R_f\) and \(R_d\) are the Carlson’s elliptic integrals.
References
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of non-Gaussian water diffusion by kurtosis analysis. NMR in Biomedicine 23(7): 698-710
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.
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
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].
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])
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
Calculated mean kurtosis tensor.
Notes
The MKT is defined as [1]:
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
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
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
The gradient table for this prediction
The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1
Notes
The predicted signal is given by:
\(\mathbf{D(n)}\) and \(\mathbf{K(n)}\) can be computed from the DT and KT using the following equations:
and
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].
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])
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
If True, RK is calculated using its analytical solution, otherwise an exact numerical estimator is used (see Notes). Default is set to True
Calculated RK.
Notes
RK is defined as the average of the directional kurtosis perpendicular to the fiber’s main direction e1 [1], [2]:
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]:
where:
and
References
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of non-Gaussian water diffusion by kurtosis analysis. NMR in Biomedicine 23(7): 698-710
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
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
DiffusionKurtosisModeldipy.reconst.dki.DiffusionKurtosisModel(gtab, fit_method='WLS', *args, **kwargs)Bases: dipy.reconst.base.ReconstModel
Class for the Diffusion Kurtosis Model
Methods
|
Fit method of the DKI model class |
|
Predict a signal for this DKI model class instance given parameters. |
__init__(gtab, fit_method='WLS', *args, **kwargs)Diffusion Kurtosis Tensor Model [1]
str can be one of the following: ‘OLS’ or ‘ULLS’ for ordinary least squares
dki.ols_fit_dki
fit_method. See dki.ols_fit_dki, dki.wls_fit_dki for details
References
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
The measured signal from one voxel.
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.
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
The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1
ReconstModeldipy.reconst.dki.ReconstModel(gtab)Bases: object
Abstract class for signal reconstruction models
Methods
fit |
TensorFitdipy.reconst.dki.TensorFit(model, model_params, model_S0=None)Bases: object
directionsFor tracking - return the primary direction in each voxel
evalsReturns the eigenvalues of the tensor as an array
evecsReturns the eigenvectors of the tensor as an array, columnwise
quadratic_formCalculates the 3x3 diffusion tensor for each voxel
Methods
|
Axial diffusivity (AD) calculated from cached eigenvalues. |
|
Calculate the apparent diffusion coefficient (ADC) in each direction on |
|
Color fractional anisotropy of diffusion tensor |
|
Fractional anisotropy (FA) calculated from cached eigenvalues. |
|
Geodesic anisotropy (GA) calculated from cached eigenvalues. |
|
|
|
Mean diffusivity (MD) calculated from cached eigenvalues. |
|
Tensor mode calculated from cached eigenvalues. |
|
The diffusion orientation distribution function (dODF). |
|
|
|
Given a model fit, predict the signal on the vertices of a sphere |
|
Radial diffusivity (RD) calculated from cached eigenvalues. |
|
|
|
Trace of the tensor calculated from cached eigenvalues. |
lower_triangular |
ad()Axial diffusivity (AD) calculated from cached eigenvalues.
Calculated AD.
Notes
RD is calculated with the following equation:
adc(sphere)Calculate the apparent diffusion coefficient (ADC) in each direction on the sphere for each voxel in the data
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.
linearity()Calculated linearity of the diffusion tensor [1].
Notes
Linearity is calculated with the following equation:
References
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.
md()Mean diffusivity (MD) calculated from cached eigenvalues.
Calculated MD.
Notes
MD is calculated with the following equation:
odf(sphere)The diffusion orientation distribution function (dODF). This is an estimate of the diffusion distance in each direction
The dODF is calculated in the vertices of this input.
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
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
Descoteaux, M. (2008). PhD Thesis: High Angular Resolution Diffusion MRI: from Local Estimation to Segmentation and Tractography. ftp://ftp-sop.inria.fr/athena/Publications/PhDs/descoteaux_thesis.pdf
planarity()Calculated sphericity of the diffusion tensor [1].
Notes
Sphericity is calculated with the following equation:
References
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
This encodes the directions for which a prediction is made
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.
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:
Where: .. math
ADC = heta Q heta^T
: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
rd()Radial diffusivity (RD) calculated from cached eigenvalues.
Calculated RD.
Notes
RD is calculated with the following equation:
sphericity()Calculated sphericity of the diffusion tensor [1].
Notes
Sphericity is calculated with the following equation:
References
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.
dipy.reconst.dki.Wcons(k_elements)Construct the full 4D kurtosis tensors from its 15 independent elements
elements of the kurtosis tensor in the following order:
& … \ & W_{xyyy} & W_{yyyz} & W_{xzzz} & W_{yzzz} & W_{xxyy} & … \ & W_{xxzz} & W_{yyzz} & W_{xxyz} & W_{xyyz} & W_{xyzz} & & )end{matrix}
Full 4D kurtosis tensor
dipy.reconst.dki.Wrotate(kt, Basis)Rotate a kurtosis tensor from the standard Cartesian coordinate system to another coordinate system basis
Vector with the 15 independent elements of the kurtosis tensor
Vectors of the basis column-wise oriented
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.
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:
& … \ & 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
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.
Array containing the 15 independent elements of the kurtosis tensor
Rotated kurtosis tensor element index i (0 for x, 1 for y, 2 for z)
Rotated kurtosis tensor element index j (0 for x, 1 for y, 2 for z)
Rotated kurtosis tensor element index k (0 for x, 1 for y, 2 for z)
Rotated kurtosis tensor element index l (0 for x, 1 for y, 2 for z)
Vectors of the basis column-wise oriented
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
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].
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
The AKC will be calculated for each of the vertices in the sphere
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
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])
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]:
where \(W_{ijkl}\) are the elements of the kurtosis tensor, MD the mean diffusivity and ADC the apparent diffusion coefficent computed as:
where \(D_{ij}\) are the elements of the diffusion tensor.
References
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
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
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].
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
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])
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
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.
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 can be sampled from the principal axis of the diffusion tensor:
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
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of non-Gaussian water diffusion by kurtosis analysis. NMR in Biomedicine 23(7): 698-710
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
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
dipy.reconst.dki.carlson_rd(x, y, z, errtol=0.0001)Computes the Carlson’s incomplete elliptic integral of the second kind defined as:
First independent variable of the integral.
Second independent variable of the integral.
Third independent variable of the integral.
Error tolerance. Integral is computed with relative error less in magnitude than the defined value
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.
dipy.reconst.dki.carlson_rf(x, y, z, errtol=0.0003)Computes the Carlson’s incomplete elliptic integral of the first kind defined as:
First independent variable of the integral.
Second independent variable of the integral.
Third independent variable of the integral.
Error tolerance. Integral is computed with relative error less in magnitude than the defined value
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
Carlson, B.C., 1994. Numerical computation of real or complex elliptic integrals. arXiv:math/9409227 [math.CA]
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\)
x coordinate in Cartesian space
y coordinate in Cartesian space
z coordinate
radius
inclination (polar) angle
azimuth angle
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
The number of different b-values you are checking for.
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)
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\).
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).
Hermitian matrix representing a diffusion tensor.
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.
Eigenvalues from eigen decomposition of the tensor. Negative eigenvalues are replaced by zero. Sorted from largest to smallest.
Associated eigenvectors from eigen decomposition of the tensor. Eigenvectors are columnar (e.g. eigvecs[…, :, j] is associated with eigvals[…, j])
dipy.reconst.dki.design_matrix(gtab)Construct B design matrix for DKI.
Measurement directions.
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)
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].
elements of the diffusion tensor of the voxel.
g directions of a Sphere in Cartesian coordinates
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
Apparent diffusion coefficient (adc) in all g directions of a sphere for a single voxel.
References
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].
elements of the diffusion tensor of the voxel.
elements of the kurtosis tensor of the voxel.
g directions of a Sphere in Cartesian coordinates
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]_)
Apparent diffusion coefficient (adc) in all g directions of a sphere for a single voxel.
Apparent diffusion variance coefficient (advc) in all g directions of a sphere for a single voxel.
Apparent diffusion variance (adv) in all g directions of a sphere for a single voxel.
References
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].
elements of the diffusion tensor of the voxel.
mean diffusivity of the voxel
elements of the kurtosis tensor of the voxel.
g directions of a Sphere in Cartesian coordinates
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
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])
Apparent diffusion coefficient (adc) in all g directions of a sphere for a single voxel.
Apparent diffusion variance (advc) in all g directions of a sphere for a single voxel.
Apparent kurtosis coefficient (AKC) in all g directions of a sphere for a single voxel.
References
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
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
dipy.reconst.dki.dki_prediction(dki_params, gtab, S0=1.0)Predict a signal given diffusion kurtosis imaging parameters.
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
The gradient table for this prediction
The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1
Simulated signal based on the DKI model:
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.
Unique elements of the tensors
3 by 3 tensors
dipy.reconst.dki.get_fnames(name='small_64D')Provide full paths to example or test datasets.
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
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
dipy.reconst.dki.get_sphere(name='symmetric362')provide triangulated spheres
which sphere - one of: * ‘symmetric362’ * ‘symmetric642’ * ‘symmetric724’ * ‘repulsion724’ * ‘repulsion100’ * ‘repulsion200’
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"
dipy.reconst.dki.kurtosis_fractional_anisotropy(dki_params)Computes the anisotropy of the kurtosis tensor (KFA) [1].
All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
Calculated mean kurtosis tensor.
Notes
The KFA is defined as [1]:
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
dipy.reconst.dki.kurtosis_maximum(dki_params, sphere='repulsion100', gtol=0.01, mask=None)Computes kurtosis maximum value
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
The sphere providing sample directions for the initial search of the maximal value of kurtosis.
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
A boolean array used to mark the coordinates in the data that should be analyzed that has the shape dki_params.shape[:-1]
kurtosis tensor maximum value
Cartesian coordinates of the direction of the maximal kurtosis value
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.
The function evaluated on a set of discrete points.
The set of neighbor relations between the points. Every edge, ie edges[i, :], is a pair of neighboring points.
Value of odf at a maximum point. Peak values is sorted in descending order.
Indices of maximum points. Sorted in the same order as peak_values so odf[peak_indices[i]] == peak_values[i].
See also
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).
dipy.reconst.dki.lower_triangular(tensor, b0=None)Returns the six lower triangular values of the tensor and a dummy variable if b0 is not None
a collection of 3, 3 diffusion tensors
if b0 is not none log(b0) is returned as the dummy variable
If b0 is none, then the shape will be (…, 6) otherwise (…, 7)
dipy.reconst.dki.mean_diffusivity(evals, axis=-1)Mean Diffusivity (MD) of a diffusion tensor.
Eigenvalues of a diffusion tensor.
Axis of evals which contains 3 eigenvalues.
Calculated MD.
Notes
MD is calculated with the following equation:
dipy.reconst.dki.mean_kurtosis(dki_params, min_kurtosis=-0.42857142857142855, max_kurtosis=3, analytical=True)Computes mean Kurtosis (MK) from the kurtosis tensor.
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
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])
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
If True, MK is calculated using its analytical solution, otherwise an exact numerical estimator is used (see Notes). Default is set to True
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]_:
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:
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
where \(R_f\) and \(R_d\) are the Carlson’s elliptic integrals.
References
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of non-Gaussian water diffusion by kurtosis analysis. NMR in Biomedicine 23(7): 698-710
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.
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
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
dipy.reconst.dki.mean_kurtosis_tensor(dki_params, min_kurtosis=-0.42857142857142855, max_kurtosis=10)Computes mean of the kurtosis tensor (MKT) [1].
All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
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])
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
Calculated mean kurtosis tensor.
Notes
The MKT is defined as [1]:
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
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
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
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.
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)
dipy.reconst.dki.nlls_fit_tensor(design_matrix, data, weighting=None, sigma=None, jac=True, return_S0_hat=False)Fit the cumulant expansion params (e.g. DTI, DKI) using non-linear least-squares.
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)
Data or response variables holding the data. Note that the last dimension should contain the data. It makes no copies of data.
the weighting scheme to use in considering the squared-error. Default behavior is to use uniform weighting. Other options: ‘sigma’ ‘gmm’
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).
Use the Jacobian? Default: True
Boolean to return (True) or not (False) the S0 values for the fit.
voxel.
dipy.reconst.dki.ols_fit_dki(design_matrix, data)Computes the diffusion and kurtosis tensors using an ordinary linear least squares (OLS) approach 1.
Design matrix holding the covariants used to solve for the regression coefficients.
Data or response variables holding the data. Note that the last dimension should contain the data. It makes no copies of data.
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
See also
wls_fit_dki, nls_fit_dkiReferences
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
dipy.reconst.dki.perpendicular_directions(v, num=30, half=False)Computes n evenly spaced perpendicular directions relative to a given vector v
Array containing the three cartesian coordinates of vector v
Number of perpendicular directions to generate
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
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:
The first axis is vector v
The second axis is defined as the normalized vector given by the
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:
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>`_.
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].
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
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])
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
If True, RK is calculated using its analytical solution, otherwise an exact numerical estimator is used (see Notes). Default is set to True.
Calculated RK.
Notes
RK is defined as the average of the directional kurtosis perpendicular to the fiber’s main direction e1 [1], [2]:
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]:
where:
and
References
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of non-Gaussian water diffusion by kurtosis analysis. NMR in Biomedicine 23(7): 698-710
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
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
dipy.reconst.dki.restore_fit_tensor(design_matrix, data, sigma=None, jac=True, return_S0_hat=False)Use the RESTORE algorithm [1] to calculate a robust tensor fit
Design matrix holding the covariants used to solve for the regression coefficients.
Data or response variables holding the data. Note that the last dimension should contain the data. It makes no copies of data.
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).
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
Boolean to return (True) or not (False) the S0 values for the fit.
References
estimation of tensors by outlier rejection. MRM, 53: 1088-95.
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’
radius
inclination or polar angle
azimuth angle
x coordinate(s) in Cartesion space
y coordinate(s) in Cartesian space
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.
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
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
Eigenvalues from eigen decomposition of the tensor.
Associated eigenvectors from eigen decomposition of the tensor. Eigenvectors are columnar (e.g. eigvecs[:,j] is associated with eigvals[j])
Fifteen elements of the kurtosis tensor
dipy.reconst.dki.vec_val_vect()Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs
containing tensor in last two dimensions; M, N usually equal to (3, 3)
diagonal values carried in last dimension, ... shape above must
match that for vecs
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).
... dimensions of vecs, valsN dimensions of vecs, valsExamples
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]])
dipy.reconst.dki.wls_fit_dki(design_matrix, data)Computes the diffusion and kurtosis tensors using a weighted linear least squares (WLS) approach 1.
Design matrix holding the covariants used to solve for the regression coefficients.
Data or response variables holding the data. Note that the last dimension should contain the data. It makes no copies of data.
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
2013. Weighted linear least squares estimation of diffusion MRI parameters: Strengths, limitations, and pitfalls. Magn Reson Med 81, 335-346.
DiffusionKurtosisFitdipy.reconst.dki_micro.DiffusionKurtosisFit(model, model_params)Bases: dipy.reconst.dti.TensorFit
Class for fitting the Diffusion Kurtosis Model
directionsFor tracking - return the primary direction in each voxel
evalsReturns the eigenvalues of the tensor as an array
evecsReturns the eigenvectors of the tensor as an array, columnwise
kfaReturns the kurtosis tensor (KFA) 1.
ktReturns the 15 independent elements of the kurtosis tensor as an array
quadratic_formCalculates the 3x3 diffusion tensor for each voxel
Methods
|
Axial diffusivity (AD) calculated from cached eigenvalues. |
|
Calculate the apparent diffusion coefficient (ADC) in each direction on |
|
Axial Kurtosis (AK) of a diffusion kurtosis tensor [1]. |
|
Calculates the apparent kurtosis coefficient (AKC) in each direction on the sphere for each voxel in the data |
|
Color fractional anisotropy of diffusion tensor |
|
Fractional anisotropy (FA) calculated from cached eigenvalues. |
|
Geodesic anisotropy (GA) calculated from cached eigenvalues. |
|
Computes the maximum value of a single voxel kurtosis tensor |
|
|
|
Mean diffusivity (MD) calculated from cached eigenvalues. |
|
Computes mean Kurtosis (MK) from the kurtosis tensor. |
|
Computes mean of the kurtosis tensor (MKT) [1]. |
|
Tensor mode calculated from cached eigenvalues. |
|
The diffusion orientation distribution function (dODF). |
|
|
|
Given a DKI model fit, predict the signal on the vertices of a gradient table |
|
Radial diffusivity (RD) calculated from cached eigenvalues. |
|
Radial Kurtosis (RK) of a diffusion kurtosis tensor [1]. |
|
|
|
Trace of the tensor calculated from cached eigenvalues. |
lower_triangular |
__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
Class instance containing the Diffusion Kurtosis Model for the fit
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].
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])
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
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.
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 can be sampled from the principal axis of the diffusion tensor:
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
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of non-Gaussian water diffusion by kurtosis analysis. NMR in Biomedicine 23(7): 698-710
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
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
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:
where \(W_{ijkl}\) are the elements of the kurtosis tensor, MD the mean diffusivity and ADC the apparent diffusion coefficent computed as:
where \(D_{ij}\) are the elements of the diffusion tensor.
kfaReturns the kurtosis tensor (KFA) [1].
Notes
The KFA is defined as [1]:
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
kmax(sphere='repulsion100', gtol=1e-05, mask=None)Computes the maximum value of a single voxel kurtosis tensor
The sphere providing sample directions for the initial search of the maximum value of kurtosis.
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
kurtosis tensor maximum value
mk(min_kurtosis=-0.42857142857142855, max_kurtosis=10, analytical=True)Computes mean Kurtosis (MK) from the kurtosis tensor.
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])
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
If True, MK is calculated using its analytical solution, otherwise an exact numerical estimator is used (see Notes). Default is set to True.
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]_:
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:
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
where \(R_f\) and \(R_d\) are the Carlson’s elliptic integrals.
References
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of non-Gaussian water diffusion by kurtosis analysis. NMR in Biomedicine 23(7): 698-710
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.
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
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].
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])
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
Calculated mean kurtosis tensor.
Notes
The MKT is defined as [1]:
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
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
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
The gradient table for this prediction
The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1
Notes
The predicted signal is given by:
\(\mathbf{D(n)}\) and \(\mathbf{K(n)}\) can be computed from the DT and KT using the following equations:
and
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].
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])
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
If True, RK is calculated using its analytical solution, otherwise an exact numerical estimator is used (see Notes). Default is set to True
Calculated RK.
Notes
RK is defined as the average of the directional kurtosis perpendicular to the fiber’s main direction e1 [1], [2]:
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]:
where:
and
References
Jensen, J.H., Helpern, J.A., 2010. MRI quantification of non-Gaussian water diffusion by kurtosis analysis. NMR in Biomedicine 23(7): 698-710
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
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
DiffusionKurtosisModeldipy.reconst.dki_micro.DiffusionKurtosisModel(gtab, fit_method='WLS', *args, **kwargs)Bases: dipy.reconst.base.ReconstModel
Class for the Diffusion Kurtosis Model
Methods
|
Fit method of the DKI model class |
|
Predict a signal for this DKI model class instance given parameters. |
__init__(gtab, fit_method='WLS', *args, **kwargs)Diffusion Kurtosis Tensor Model [1]
str can be one of the following: ‘OLS’ or ‘ULLS’ for ordinary least squares
dki.ols_fit_dki
fit_method. See dki.ols_fit_dki, dki.wls_fit_dki for details
References
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
The measured signal from one voxel.
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.
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
The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1
KurtosisMicrostructuralFitdipy.reconst.dki_micro.KurtosisMicrostructuralFit(model, model_params)Bases: dipy.reconst.dki.DiffusionKurtosisFit
Class for fitting the Diffusion Kurtosis Microstructural Model
awfReturns the volume fraction of the restricted diffusion compartment also known as axonal water fraction.
axonal_diffusivityReturns the axonal diffusivity defined as the restricted diffusion tensor trace 1.
directionsFor tracking - return the primary direction in each voxel
evalsReturns the eigenvalues of the tensor as an array
evecsReturns the eigenvectors of the tensor as an array, columnwise
hindered_adReturns the axial diffusivity of the hindered compartment.
hindered_evalsReturns the eigenvalues of the hindered diffusion compartment.
hindered_rdReturns the radial diffusivity of the hindered compartment.
kfaReturns the kurtosis tensor (KFA) 1.
ktReturns the 15 independent elements of the kurtosis tensor as an array
quadratic_formCalculates the 3x3 diffusion tensor for each voxel
restricted_evalsReturns the eigenvalues of the restricted diffusion compartment.
tortuosityReturns the tortuosity of the hindered diffusion which is defined by ADe / RDe, where ADe and RDe are the axial and radial diffusivities of the hindered compartment 1.
Methods
|
Axial diffusivity (AD) calculated from cached eigenvalues. |
|
Calculate the apparent diffusion coefficient (ADC) in each direction on |
|
Axial Kurtosis (AK) of a diffusion kurtosis tensor [R0b1a747e81c9-1]. |
|
Calculates the apparent kurtosis coefficient (AKC) in each direction on the sphere for each voxel in the data |
|
Color fractional anisotropy of diffusion tensor |
|
Fractional anisotropy (FA) calculated from cached eigenvalues. |
|
Geodesic anisotropy (GA) calculated from cached eigenvalues. |
|
Computes the maximum value of a single voxel kurtosis tensor |
|
|
|
Mean diffusivity (MD) calculated from cached eigenvalues. |
|
Computes mean Kurtosis (MK) from the kurtosis tensor. |
|
Computes mean of the kurtosis tensor (MKT) [R8b3dd90f2e0d-1]. |
|
Tensor mode calculated from cached eigenvalues. |
|
The diffusion orientation distribution function (dODF). |
|
|
|
Given a DKI microstructural model fit, predict the signal on the vertices of a gradient table |
|
Radial diffusivity (RD) calculated from cached eigenvalues. |
|
Radial Kurtosis (RK) of a diffusion kurtosis tensor [Rc4101656d30e-1]. |
|
|
|
Trace of the tensor calculated from cached eigenvalues. |
lower_triangular |
__init__(model, model_params)Initialize a KurtosisMicrostructural Fit class instance.
Class instance containing the Diffusion Kurtosis Model for the fit
All parameters estimated from the diffusion kurtosis microstructural model. Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
Six elements of the hindered diffusion tensor
Six elements of the restricted diffusion tensor
Axonal water fraction
Notes
In the original article of DKI microstructural model [1], the hindered and restricted tensors were definde as the intra-cellular and extra-cellular diffusion compartments respectively.
References
Fieremans, E., Jensen, J.H., Helpern, J.A., 2011. White Matter Characterization with Diffusion Kurtosis Imaging. Neuroimage 58(1): 177-188. doi:10.1016/j.neuroimage.2011.06.006
awfReturns the volume fraction of the restricted diffusion compartment also known as axonal water fraction.
Notes
The volume fraction of the restricted diffusion compartment can be seem as the volume fraction of the intra-cellular compartment [1].
References
Fieremans, E., Jensen, J.H., Helpern, J.A., 2011. White Matter Characterization with Diffusion Kurtosis Imaging. Neuroimage 58(1): 177-188. doi:10.1016/j.neuroimage.2011.06.006
axonal_diffusivityReturns the axonal diffusivity defined as the restricted diffusion tensor trace [1].
References
Fieremans, E., Jensen, J.H., Helpern, J.A., 2011. White Matter Characterization with Diffusion Kurtosis Imaging. Neuroimage 58(1): 177-188. doi:10.1016/j.neuroimage.2011.06.006
hindered_adReturns the axial diffusivity of the hindered compartment.
Notes
The hindered diffusion tensor can be seem as the tissue’s extra-cellular diffusion compartment [1].
References
Fieremans, E., Jensen, J.H., Helpern, J.A., 2011. White Matter Characterization with Diffusion Kurtosis Imaging. Neuroimage 58(1): 177-188. doi:10.1016/j.neuroimage.2011.06.006
hindered_evalsReturns the eigenvalues of the hindered diffusion compartment.
Notes
The hindered diffusion tensor can be seem as the tissue’s extra-cellular diffusion compartment [1].
References
Fieremans, E., Jensen, J.H., Helpern, J.A., 2011. White Matter Characterization with Diffusion Kurtosis Imaging. Neuroimage 58(1): 177-188. doi:10.1016/j.neuroimage.2011.06.006
hindered_rdReturns the radial diffusivity of the hindered compartment.
Notes
The hindered diffusion tensor can be seem as the tissue’s extra-cellular diffusion compartment [1].
References
Fieremans, E., Jensen, J.H., Helpern, J.A., 2011. White Matter Characterization with Diffusion Kurtosis Imaging. Neuroimage 58(1): 177-188. doi:10.1016/j.neuroimage.2011.06.006
predict(gtab, S0=1.0)Given a DKI microstructural model fit, predict the signal on the vertices of a gradient table
The gradient table for this prediction
The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1
Notes
The predicted signal is given by:
\(S(\theta, b) = S_0 * [f * e^{-b ADC_{r}} + (1-f) * e^{-b ADC_{h}]\), where \(ADC_{r}\) and \(ADC_{h}\) are the apparent diffusion coefficients of the diffusion hindered and restricted compartment for a given direction \(\theta\), \(b\) is the b value provided in the GradientTable input for that direction, \(f\) is the volume fraction of the restricted diffusion compartment (also known as the axonal water fraction).
restricted_evalsReturns the eigenvalues of the restricted diffusion compartment.
Notes
The restricted diffusion tensor can be seem as the tissue’s intra-cellular diffusion compartment [1].
References
Fieremans, E., Jensen, J.H., Helpern, J.A., 2011. White Matter Characterization with Diffusion Kurtosis Imaging. Neuroimage 58(1): 177-188. doi:10.1016/j.neuroimage.2011.06.006
tortuosityReturns the tortuosity of the hindered diffusion which is defined by ADe / RDe, where ADe and RDe are the axial and radial diffusivities of the hindered compartment [1].
Notes
The hindered diffusion tensor can be seem as the tissue’s extra-cellular diffusion compartment [1].
References
KurtosisMicrostructureModeldipy.reconst.dki_micro.KurtosisMicrostructureModel(gtab, fit_method='WLS', *args, **kwargs)Bases: dipy.reconst.dki.DiffusionKurtosisModel
Class for the Diffusion Kurtosis Microstructural Model
Methods
|
Fit method of the Diffusion Kurtosis Microstructural Model |
|
Predict a signal for the DKI microstructural model class instance given parameters. |
__init__(gtab, fit_method='WLS', *args, **kwargs)Initialize a KurtosisMicrostrutureModel class instance [1].
str can be one of the following: ‘OLS’ or ‘ULLS’ to fit the diffusion tensor and kurtosis tensor using the ordinary linear least squares solution
dki.ols_fit_dki
‘WLS’ or ‘UWLLS’ to fit the diffusion tensor and kurtosis tensor using the ordinary linear least squares solution
dki.wls_fit_dki
fit_method. See dki.ols_fit_dki, dki.wls_fit_dki for details
References
Fieremans, E., Jensen, J.H., Helpern, J.A., 2011. White Matter Characterization with Diffusion Kurtosis Imaging. Neuroimage 58(1): 177-188. doi:10.1016/j.neuroimage.2011.06.006
fit(data, mask=None, sphere='repulsion100', gtol=0.01, awf_only=False)Fit method of the Diffusion Kurtosis Microstructural Model
An 4D matrix containing the diffusion-weighted data.
A boolean array used to mark the coordinates in the data that should be analyzed that has the shape data.shape[-1]
The sphere providing sample directions for the initial search of the maximal value of kurtosis.
This input is to refine kurtosis maxima under the precision of the directions sampled on the sphere class instance. The gradient of the convergence procedure must be less than gtol before successful termination. If gtol is None, fiber direction is directly taken from the initial sampled directions of the given sphere object
If set to true only the axonal volume fraction is computed from the kurtosis tensor. Default = False
predict(params, S0=1.0)Predict a signal for the DKI microstructural model class instance given parameters.
All parameters estimated from the diffusion kurtosis microstructural model. Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
Six elements of the hindered diffusion tensor
Six elements of the restricted diffusion tensor
Axonal water fraction
The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1
Notes
In the original article of DKI microstructural model [1], the hindered and restricted tensors were definde as the intra-cellular and extra-cellular diffusion compartments respectively.
References
Fieremans, E., Jensen, J.H., Helpern, J.A., 2011. White Matter Characterization with Diffusion Kurtosis Imaging. Neuroimage 58(1): 177-188. doi:10.1016/j.neuroimage.2011.06.006
dipy.reconst.dki_micro.axial_diffusivity(evals, axis=-1)Axial Diffusivity (AD) of a diffusion tensor. Also called parallel diffusivity.
Eigenvalues of a diffusion tensor, must be sorted in descending order along axis.
Axis of evals which contains 3 eigenvalues.
Calculated AD.
Notes
AD is calculated with the following equation:
dipy.reconst.dki_micro.axonal_water_fraction(dki_params, sphere='repulsion100', gtol=0.01, mask=None)Computes the axonal water fraction from DKI [1].
All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
The sphere providing sample directions for the initial search of the maximal value of kurtosis.
This input is to refine kurtosis maxima under the precision of the directions sampled on the sphere class instance. The gradient of the convergence procedure must be less than gtol before successful termination. If gtol is None, fiber direction is directly taken from the initial sampled directions of the given sphere object
A boolean array used to mark the coordinates in the data that should be analyzed that has the shape dki_params.shape[:-1]
Axonal Water Fraction
References
dipy.reconst.dki_micro.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).
Hermitian matrix representing a diffusion tensor.
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.
Eigenvalues from eigen decomposition of the tensor. Negative eigenvalues are replaced by zero. Sorted from largest to smallest.
Associated eigenvectors from eigen decomposition of the tensor. Eigenvectors are columnar (e.g. eigvecs[…, :, j] is associated with eigvals[…, j])
dipy.reconst.dki_micro.diffusion_components(dki_params, sphere='repulsion100', awf=None, mask=None)Extracts the restricted and hindered diffusion tensors of well aligned fibers from diffusion kurtosis imaging parameters [1].
All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows:
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
The sphere providing sample directions to sample the restricted and hindered cellular diffusion tensors. For more details see Fieremans et al., 2011.
Array containing values of the axonal water fraction that has the shape
dki_params.shape[:-1]. If not given this will be automatically computed
using axonal_water_fraction()” with function’s default precision.
A boolean array used to mark the coordinates in the data that should be analyzed that has the shape dki_params.shape[:-1]
Parameters of the hindered diffusion tensor.
Parameters of the restricted diffusion tensor.
Notes
In the original article of DKI microstructural model [1], the hindered and restricted tensors were definde as the intra-cellular and extra-cellular diffusion compartments respectively.
References
dipy.reconst.dki_micro.directional_diffusion(dt, V, min_diffusivity=0)Calculates the apparent diffusion coefficient (adc) in each direction of a sphere for a single voxel [1].
elements of the diffusion tensor of the voxel.
g directions of a Sphere in Cartesian coordinates
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
Apparent diffusion coefficient (adc) in all g directions of a sphere for a single voxel.
References
dipy.reconst.dki_micro.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].
elements of the diffusion tensor of the voxel.
mean diffusivity of the voxel
elements of the kurtosis tensor of the voxel.
g directions of a Sphere in Cartesian coordinates
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
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])
Apparent diffusion coefficient (adc) in all g directions of a sphere for a single voxel.
Apparent diffusion variance (advc) in all g directions of a sphere for a single voxel.
Apparent kurtosis coefficient (AKC) in all g directions of a sphere for a single voxel.
References
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
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
dipy.reconst.dki_micro.dkimicro_prediction(params, gtab, S0=1)Signal prediction given the DKI microstructure model parameters.
Three diffusion tensor’s eigenvalues
Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector
Fifteen elements of the kurtosis tensor
Six elements of the hindered diffusion tensor
Six elements of the restricted diffusion tensor
Axonal water fraction
The gradient table for this prediction
The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1
Simulated signal based on the DKI microstructure model
Notes
1) The predicted signal is given by: \(S(\theta, b) = S_0 * [f * e^{-b ADC_{r}} + (1-f) * e^{-b ADC_{h}]\), where :math:` ADC_{r} and ADC_{h} are the apparent diffusion coefficients of the diffusion hindered and restricted compartment for a given direction theta:math:, b:math: is the b value provided in the GradientTable input for that direction, `f$ is the volume fraction of the restricted diffusion compartment (also known as the axonal water fraction).
2) In the original article of DKI microstructural model 1, the hindered and restricted tensors were definde as the intra-cellular and extra-cellular diffusion compartments respectively.
dipy.reconst.dki_micro.dti_design_matrix(gtab, dtype=None)Constructs design matrix for DTI weighted least squares or least squares fitting. (Basser et al., 1994a)
Parameter to control the dtype of returned designed matrix
Design matrix or B matrix assuming Gaussian distributed tensor model design_matrix[j, :] = (Bxx, Byy, Bzz, Bxy, Bxz, Byz, dummy)
dipy.reconst.dki_micro.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.
Unique elements of the tensors
3 by 3 tensors
dipy.reconst.dki_micro.get_sphere(name='symmetric362')provide triangulated spheres
which sphere - one of: * ‘symmetric362’ * ‘symmetric642’ * ‘symmetric724’ * ‘repulsion724’ * ‘repulsion100’ * ‘repulsion200’
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"
dipy.reconst.dki_micro.kurtosis_maximum(dki_params, sphere='repulsion100', gtol=0.01, mask=None)Computes kurtosis maximum value
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
The sphere providing sample directions for the initial search of the maximal value of kurtosis.
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
A boolean array used to mark the coordinates in the data that should be analyzed that has the shape dki_params.shape[:-1]
kurtosis tensor maximum value
Cartesian coordinates of the direction of the maximal kurtosis value
dipy.reconst.dki_micro.lower_triangular(tensor, b0=None)Returns the six lower triangular values of the tensor and a dummy variable if b0 is not None
a collection of 3, 3 diffusion tensors
if b0 is not none log(b0) is returned as the dummy variable
If b0 is none, then the shape will be (…, 6) otherwise (…, 7)
dipy.reconst.dki_micro.mean_diffusivity(evals, axis=-1)Mean Diffusivity (MD) of a diffusion tensor.
Eigenvalues of a diffusion tensor.
Axis of evals which contains 3 eigenvalues.
Calculated MD.
Notes
MD is calculated with the following equation:
dipy.reconst.dki_micro.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.
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)
dipy.reconst.dki_micro.radial_diffusivity(evals, axis=-1)Radial Diffusivity (RD) of a diffusion tensor. Also called perpendicular diffusivity.
Eigenvalues of a diffusion tensor, must be sorted in descending order along axis.
Axis of evals which contains 3 eigenvalues.
Calculated RD.
Notes
RD is calculated with the following equation:
dipy.reconst.dki_micro.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
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
Eigenvalues from eigen decomposition of the tensor.
Associated eigenvectors from eigen decomposition of the tensor. Eigenvectors are columnar (e.g. eigvecs[:,j] is associated with eigvals[j])
Fifteen elements of the kurtosis tensor
dipy.reconst.dki_micro.tortuosity(hindered_ad, hindered_rd)Computes the tortuosity of the hindered diffusion compartment given its axial and radial diffusivities
Array containing the values of the hindered axial diffusivity.
Array containing the values of the hindered radial diffusivity.
dipy.reconst.dki_micro.trace(evals, axis=-1)Trace of a diffusion tensor.
Eigenvalues of a diffusion tensor.
Axis of evals which contains 3 eigenvalues.
Calculated trace of the diffusion tensor.
Notes
Trace is calculated with the following equation:
dipy.reconst.dki_micro.vec_val_vect()Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs
containing tensor in last two dimensions; M, N usually equal to (3, 3)
diagonal values carried in last dimension, ... shape above must
match that for vecs
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).
... dimensions of vecs, valsN dimensions of vecs, valsExamples
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]])
Cachedipy.reconst.dsi.CacheBases: 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
Clear the cache. |
|
|
Retrieve a value from the cache. |
|
Store a value in the cache. |
cache_get(tag, key, default=None)Retrieve a value from the cache.
Description of the cached value.
Key object used to look up the cached value.
Value to be returned if no cached entry is found.
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.
Description of the cached value.
Key object used to look up the cached value.
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)
>>> c = Cache()
>>> parameters = (1, 2, 3)
>>> X1 = compute_expensive_matrix(parameters)
>>> c.cache_set('expensive_matrix', parameters, X1)
>>> X2 = c.cache_get('expensive_matrix', parameters)
>>> X1 is X2
True
DiffusionSpectrumDeconvFitdipy.reconst.dsi.DiffusionSpectrumDeconvFit(model, data)Bases: dipy.reconst.dsi.DiffusionSpectrumFit
Methods
|
Calculates the mean squared displacement on the discrete propagator |
|
Calculates the real discrete odf for a given discrete sphere |
|
Applies the 3D FFT in the q-space grid to generate the DSI diffusion propagator, remove the background noise with a hard threshold and then deconvolve the propagator with the Lucy-Richardson deconvolution algorithm |
|
Calculates the return to origin probability from the propagator, which is the propagator evaluated at zero (see Descoteaux et Al. |
|
Calculates the return to origin probability (rtop) from the signal rtop equals to the sum of all signal values |
DiffusionSpectrumDeconvModeldipy.reconst.dsi.DiffusionSpectrumDeconvModel(gtab, qgrid_size=35, r_start=4.1, r_end=13.0, r_step=0.4, filter_width=inf, normalize_peaks=False)Bases: dipy.reconst.dsi.DiffusionSpectrumModel
Methods
|
Clear the cache. |
|
Retrieve a value from the cache. |
|
Store a value in the cache. |
|
Fit method for every voxel in data |
__init__(gtab, qgrid_size=35, r_start=4.1, r_end=13.0, r_step=0.4, filter_width=inf, normalize_peaks=False)Diffusion Spectrum Deconvolution
The idea is to remove the convolution on the DSI propagator that is caused by the truncation of the q-space in the DSI sampling.
P_{dsi}(mathbf{r}) & = & S_{0}^{-1}iiintlimits_{| mathbf{q} | le mathbf{q_{max}}} S(mathbf{q})exp(-i2pimathbf{q}cdotmathbf{r})dmathbf{q} \ & = & S_{0}^{-1}iiintlimits_{mathbf{q}} left( S(mathbf{q}) cdot M(mathbf{q}) right) exp(-i2pimathbf{q}cdotmathbf{r})dmathbf{q} \ & = & P(mathbf{r}) otimes left( S_{0}^{-1}iiintlimits_{mathbf{q}} M(mathbf{q}) exp(-i2pimathbf{q}cdotmathbf{r})dmathbf{q} right) \
end{eqnarray*}
where \(\mathbf{r}\) is the displacement vector and \(\mathbf{q}\) is the wave vector which corresponds to different gradient directions, \(M(\mathbf{q})\) is a mask corresponding to your q-space sampling and \(\otimes\) is the convolution operator [1].
Gradient directions and bvalues container class
has to be an odd number. Sets the size of the q_space grid.
For example if qgrid_size is 35 then the shape of the grid will be
(35, 35, 35).
ODF is sampled radially in the PDF. This parameters shows where the sampling should start.
Radial endpoint of ODF sampling
Step size of the ODf sampling from r_start to r_end
Strength of the hanning filter
References
Canales-Rodriguez E.J et al., “Deconvolution in Diffusion
Spectrum Imaging”, Neuroimage, 2010.
Biggs David S.C. et al., “Acceleration of Iterative Image
Restoration Algorithms”, Applied Optics, vol. 36, No. 8, p. 1766-1775, 1997.
DiffusionSpectrumFitdipy.reconst.dsi.DiffusionSpectrumFit(model, data)Bases: dipy.reconst.odf.OdfFit
Methods
|
Calculates the mean squared displacement on the discrete propagator |
|
Calculates the real discrete odf for a given discrete sphere |
|
Applies the 3D FFT in the q-space grid to generate the diffusion propagator |
|
Calculates the return to origin probability from the propagator, which is the propagator evaluated at zero (see Descoteaux et Al. |
|
Calculates the return to origin probability (rtop) from the signal rtop equals to the sum of all signal values |
__init__(model, data)Calculates PDF and ODF and other properties for a single voxel
DiffusionSpectrumModel
signal values
msd_discrete(normalized=True)Calculates the mean squared displacement on the discrete propagator
MSD:{DSI}=int_{-infty}^{infty}int_{-infty}^{infty}int_{-infty}^{infty} P(hat{mathbf{r}}) cdot hat{mathbf{r}}^{2} dr_x dr_y dr_z
end{equation}
where \(\hat{\mathbf{r}}\) is a point in the 3D Propagator space (see Wu et al. [1]).
Whether to normalize the propagator by its sum in order to obtain a pdf. Default: True
the mean square displacement
References
Wu Y. et al., “Hybrid diffusion imaging”, NeuroImage, vol 36,
617-629, 2007.
odf(sphere)Calculates the real discrete odf for a given discrete sphere
psi_{DSI}(hat{mathbf{u}})=int_{0}^{infty}P(rhat{mathbf{u}})r^{2}dr
end{equation}
where \(\hat{\mathbf{u}}\) is the unit vector which corresponds to a sphere point.
rtop_pdf(normalized=True)Calculates the return to origin probability from the propagator, which is the propagator evaluated at zero (see Descoteaux et Al. [1], Tuch [2], Wu et al. [3]) rtop = P(0)
Whether to normalize the propagator by its sum in order to obtain a pdf. Default: True.
the return to origin probability
References
Descoteaux M. et al., “Multiple q-shell diffusion propagator
imaging”, Medical Image Analysis, vol 15, No. 4, p. 603-621, 2011.
Tuch D.S., “Diffusion MRI of Complex Tissue Structure”, PhD Thesis, 2002.
Wu Y. et al., “Computation of Diffusion Function Measures
in q -Space Using Magnetic Resonance Hybrid Diffusion Imaging”, IEEE TRANSACTIONS ON MEDICAL IMAGING, vol. 27, No. 6, p. 858-865, 2008
DiffusionSpectrumModeldipy.reconst.dsi.DiffusionSpectrumModel(gtab, qgrid_size=17, r_start=2.1, r_end=6.0, r_step=0.2, filter_width=32, normalize_peaks=False)Bases: dipy.reconst.odf.OdfModel, dipy.reconst.cache.Cache
Methods
|
Clear the cache. |
|
Retrieve a value from the cache. |
|
Store a value in the cache. |
|
Fit method for every voxel in data |
__init__(gtab, qgrid_size=17, r_start=2.1, r_end=6.0, r_step=0.2, filter_width=32, normalize_peaks=False)Diffusion Spectrum Imaging
The theoretical idea underlying this method is that the diffusion propagator \(P(\mathbf{r})\) (probability density function of the average spin displacements) can be estimated by applying 3D FFT to the signal values \(S(\mathbf{q})\)
P(mathbf{r}) & = & S_{0}^{-1}int S(mathbf{q})exp(-i2pimathbf{q}cdotmathbf{r})dmathbf{r}
end{eqnarray}
where \(\mathbf{r}\) is the displacement vector and \(\mathbf{q}\) is the wave vector which corresponds to different gradient directions. Method used to calculate the ODFs. Here we implement the method proposed by Wedeen et al. [1].
The main assumption for this model is fast gradient switching and that the acquisition gradients will sit on a keyhole Cartesian grid in q_space [3].
Gradient directions and bvalues container class
has to be an odd number. Sets the size of the q_space grid.
For example if qgrid_size is 17 then the shape of the grid will be
(17, 17, 17).
ODF is sampled radially in the PDF. This parameters shows where the sampling should start.
Radial endpoint of ODF sampling
Step size of the ODf sampling from r_start to r_end
Strength of the hanning filter
See also
dipy.reconst.gqi.GeneralizedQSamplingNotes
A. Have in mind that DSI expects gradients on both hemispheres. If your gradients span only one hemisphere you need to duplicate the data and project them to the other hemisphere before calling this class. The function dipy.reconst.dsi.half_to_full_qspace can be used for this purpose.
B. If you increase the size of the grid (parameter qgrid_size) you will most likely also need to update the r_* parameters. This is because the added zero padding from the increase of gqrid_size also introduces a scaling of the PDF.
We assume that data only one b0 volume is provided.
References
Wedeen V.J et al., “Mapping Complex Tissue Architecture With
Diffusion Spectrum Magnetic Resonance Imaging”, MRM 2005.
Canales-Rodriguez E.J et al., “Deconvolution in Diffusion
Spectrum Imaging”, Neuroimage, 2010.
Garyfallidis E, “Towards an accurate brain tractography”, PhD
thesis, University of Cambridge, 2012.
Examples
In this example where we provide the data, a gradient table and a reconstruction sphere, we calculate generalized FA for the first voxel in the data with the reconstruction performed using DSI.
>>> import warnings
>>> from dipy.data import dsi_voxels, default_sphere
>>> data, gtab = dsi_voxels()
>>> from dipy.reconst.dsi import DiffusionSpectrumModel
>>> ds = DiffusionSpectrumModel(gtab)
>>> dsfit = ds.fit(data)
>>> from dipy.reconst.odf import gfa
>>> np.round(gfa(dsfit.odf(default_sphere))[0, 0, 0], 2)
0.11
OdfFitdipy.reconst.dsi.OdfFit(model, data)Bases: dipy.reconst.base.ReconstFit
Methods
|
To be implemented but specific odf models |
OdfModeldipy.reconst.dsi.OdfModel(gtab)Bases: dipy.reconst.base.ReconstModel
An abstract class to be sub-classed by specific odf models
All odf models should provide a fit method which may take data as it’s first and only argument.
Methods
|
To be implemented by specific odf models |
dipy.reconst.dsi.LR_deconv(prop, psf, numit=5, acc_factor=1)Perform Lucy-Richardson deconvolution algorithm on a 3D array.
The 3D volume to be deconvolve
The filter that will be used for the deconvolution.
Number of Lucy-Richardson iteration to perform.
Exponential acceleration factor as in [1].
References
Biggs David S.C. et al., “Acceleration of Iterative Image Restoration Algorithms”, Applied Optics, vol. 36, No. 8, p. 1766-1775, 1997.
dipy.reconst.dsi.fftn(x, shape=None, axes=None, overwrite_x=False)Return multidimensional discrete Fourier transform.
The returned array contains:
y[j_1,..,j_d] = sum[k_1=0..n_1-1, ..., k_d=0..n_d-1]
x[k_1,..,k_d] * prod[i=1..d] exp(-sqrt(-1)*2*pi/n_i * j_i * k_i)
where d = len(x.shape) and n = x.shape.
The (N-D) array to transform.
The shape of the result. If both shape and axes (see below) are
None, shape is x.shape; if shape is None but axes is
not None, then shape is scipy.take(x.shape, axes, axis=0).
If shape[i] > x.shape[i], the ith dimension is padded with zeros.
If shape[i] < x.shape[i], the ith dimension is truncated to
length shape[i].
If any element of shape is -1, the size of the corresponding
dimension of x is used.
The axes of x (y if shape is not None) along which the transform is applied. The default is over all axes.
If True, the contents of x can be destroyed. Default is False.
The (N-D) DFT of the input array.
See also
ifftnNotes
If x is real-valued, then
y[..., j_i, ...] == y[..., n_i-j_i, ...].conjugate().
Both single and double precision routines are implemented. Half precision inputs will be converted to single precision. Non-floating-point inputs will be converted to double precision. Long-double precision inputs are not supported.
Examples
>>> from scipy.fftpack import fftn, ifftn
>>> y = (-np.arange(16), 8 - np.arange(16), np.arange(16))
>>> np.allclose(y, fftn(ifftn(y)))
True
dipy.reconst.dsi.fftshift(x, axes=None)Shift the zero-frequency component to the center of the spectrum.
This function swaps half-spaces for all axes listed (defaults to all).
Note that y[0] is the Nyquist component only if len(x) is even.
Input array.
Axes over which to shift. Default is None, which shifts all axes.
The shifted array.
See also
ifftshiftThe inverse of fftshift.
Examples
>>> freqs = np.fft.fftfreq(10, 0.1)
>>> freqs
array([ 0., 1., 2., ..., -3., -2., -1.])
>>> np.fft.fftshift(freqs)
array([-5., -4., -3., -2., -1., 0., 1., 2., 3., 4.])
Shift the zero-frequency component only along the second axis:
>>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
>>> freqs
array([[ 0., 1., 2.],
[ 3., 4., -4.],
[-3., -2., -1.]])
>>> np.fft.fftshift(freqs, axes=(1,))
array([[ 2., 0., 1.],
[-4., 3., 4.],
[-1., -3., -2.]])
dipy.reconst.dsi.half_to_full_qspace(data, gtab)Half to full Cartesian grid mapping
Useful when dMRI data are provided in one qspace hemisphere as DiffusionSpectrum expects data to be in full qspace.
where (X, Y, Z) volume size and W number of gradient directions
container for b-values and b-vectors (gradient directions)
Notes
We assume here that only on b0 is provided with the initial data. If that is not the case then you will need to write your own preparation function before providing the gradients and the data to the DiffusionSpectrumModel class.
dipy.reconst.dsi.hanning_filter(gtab, filter_width, origin)create a hanning window
The signal is premultiplied by a Hanning window before Fourier transform in order to ensure a smooth attenuation of the signal at high q values.
center of qspace
where N is the number of non-b0 gradient directions
dipy.reconst.dsi.ifftshift(x, axes=None)The inverse of fftshift. Although identical for even-length x, the functions differ by one sample for odd-length x.
Input array.
Axes over which to calculate. Defaults to None, which shifts all axes.
The shifted array.
See also
fftshiftShift zero-frequency component to the center of the spectrum.
Examples
>>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
>>> freqs
array([[ 0., 1., 2.],
[ 3., 4., -4.],
[-3., -2., -1.]])
>>> np.fft.ifftshift(np.fft.fftshift(freqs))
array([[ 0., 1., 2.],
[ 3., 4., -4.],
[-3., -2., -1.]])
dipy.reconst.dsi.map_coordinates(input, coordinates, output=None, order=3, mode='constant', cval=0.0, prefilter=True)Map the input array to new coordinates by interpolation.
The array of coordinates is used to find, for each point in the output, the corresponding coordinates in the input. The value of the input at those coordinates is determined by spline interpolation of the requested order.
The shape of the output is derived from that of the coordinate array by dropping the first axis. The values of the array along the first axis are the coordinates in the input array at which the output value is found.
The input array.
The coordinates at which input is evaluated.
The array in which to place the output, or the dtype of the returned array. By default an array of the same dtype as input will be created.
The order of the spline interpolation, default is 3. The order has to be in the range 0-5.
The mode parameter determines how the input array is extended beyond its boundaries. Default is ‘constant’. Behavior for each valid value is as follows:
The input is extended by reflecting about the edge of the last pixel.
The input is extended by filling all values beyond the edge with the same constant value, defined by the cval parameter.
The input is extended by replicating the last pixel.
The input is extended by reflecting about the center of the last pixel.
The input is extended by wrapping around to the opposite edge.
Value to fill past edges of input if mode is ‘constant’. Default is 0.0.
Determines if the input array is prefiltered with spline_filter before interpolation. The default is True, which will create a temporary float64 array of filtered values if order > 1. If setting this to False, the output will be slightly blurred if order > 1, unless the input is prefiltered, i.e. it is the result of calling spline_filter on the original input.
The result of transforming the input. The shape of the output is derived from that of coordinates by dropping the first axis.
See also
spline_filter, geometric_transform, scipy.interpolateExamples
>>> from scipy import ndimage
>>> a = np.arange(12.).reshape((4, 3))
>>> a
array([[ 0., 1., 2.],
[ 3., 4., 5.],
[ 6., 7., 8.],
[ 9., 10., 11.]])
>>> ndimage.map_coordinates(a, [[0.5, 2], [0.5, 1]], order=1)
array([ 2., 7.])
Above, the interpolated value of a[0.5, 0.5] gives output[0], while a[2, 1] is output[1].
>>> inds = np.array([[0.5, 2], [0.5, 4]])
>>> ndimage.map_coordinates(a, inds, order=1, cval=-33.3)
array([ 2. , -33.3])
>>> ndimage.map_coordinates(a, inds, order=1, mode='nearest')
array([ 2., 8.])
>>> ndimage.map_coordinates(a, inds, order=1, cval=0, output=bool)
array([ True, False], dtype=bool)
dipy.reconst.dsi.pdf_odf(Pr, rradius, interp_coords)Calculates the real ODF from the diffusion propagator(PDF) Pr
probability density function
interpolation range on the radius
coordinates in the pdf for interpolating the odf
ReconstModeldipy.reconst.dti.ReconstModel(gtab)Bases: object
Abstract class for signal reconstruction models
Methods
fit |
TensorFitdipy.reconst.dti.TensorFit(model, model_params, model_S0=None)Bases: object
directionsFor tracking - return the primary direction in each voxel
evalsReturns the eigenvalues of the tensor as an array
evecsReturns the eigenvectors of the tensor as an array, columnwise
quadratic_formCalculates the 3x3 diffusion tensor for each voxel
Methods
|
Axial diffusivity (AD) calculated from cached eigenvalues. |
|
Calculate the apparent diffusion coefficient (ADC) in each direction on |
|
Color fractional anisotropy of diffusion tensor |
|
Fractional anisotropy (FA) calculated from cached eigenvalues. |
|
Geodesic anisotropy (GA) calculated from cached eigenvalues. |
|
|
|
Mean diffusivity (MD) calculated from cached eigenvalues. |
|
Tensor mode calculated from cached eigenvalues. |
|
The diffusion orientation distribution function (dODF). |
|
|
|
Given a model fit, predict the signal on the vertices of a sphere |
|
Radial diffusivity (RD) calculated from cached eigenvalues. |
|
|
|
Trace of the tensor calculated from cached eigenvalues. |
lower_triangular |
ad()Axial diffusivity (AD) calculated from cached eigenvalues.
Calculated AD.
Notes
RD is calculated with the following equation:
adc(sphere)Calculate the apparent diffusion coefficient (ADC) in each direction on the sphere for each voxel in the data
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.
linearity()Calculated linearity of the diffusion tensor [1].
Notes
Linearity is calculated with the following equation:
References
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.
md()Mean diffusivity (MD) calculated from cached eigenvalues.
Calculated MD.
Notes
MD is calculated with the following equation:
odf(sphere)The diffusion orientation distribution function (dODF). This is an estimate of the diffusion distance in each direction
The dODF is calculated in the vertices of this input.
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
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
Descoteaux, M. (2008). PhD Thesis: High Angular Resolution Diffusion MRI: from Local Estimation to Segmentation and Tractography. ftp://ftp-sop.inria.fr/athena/Publications/PhDs/descoteaux_thesis.pdf
planarity()Calculated sphericity of the diffusion tensor [1].
Notes
Sphericity is calculated with the following equation:
References
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
This encodes the directions for which a prediction is made
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.
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:
Where: .. math
ADC = heta Q heta^T
: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
rd()Radial diffusivity (RD) calculated from cached eigenvalues.
Calculated RD.
Notes
RD is calculated with the following equation:
sphericity()Calculated sphericity of the diffusion tensor [1].
Notes
Sphericity is calculated with the following equation:
References
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.
TensorModeldipy.reconst.dti.TensorModel(gtab, fit_method='WLS', return_S0_hat=False, *args, **kwargs)Bases: dipy.reconst.base.ReconstModel
Diffusion Tensor
Methods
|
Fit method of the DTI model class |
|
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].
str can be one of the following:
dti.wls_fit_tensor()
dti.ols_fit_tensor()
dti.nlls_fit_tensor()
fitting [3]
dti.restore_fit_tensor()
Boolean to return (True) or not (False) the S0 values for the fit.
fit_method. See dti.wls_fit_tensor, dti.ols_fit_tensor for details
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
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.
Basser, P., Pierpaoli, C., 1996. Microstructural and physiological features of tissues elucidated by quantitative diffusion-tensor MRI. Journal of Magnetic Resonance 111, 209-219.
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
The measured signal from one voxel.
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.
The last dimension should have 12 tensor parameters: 3 eigenvalues, followed by the 3 eigenvectors
The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1
dipy.reconst.dti.apparent_diffusion_coef(q_form, sphere)Calculate the apparent diffusion coefficient (ADC) in each direction of a sphere.
The quadratic form of a tensor, or an array with quadratic forms of tensors. Should be of shape (…, 3, 3)
The ADC will be calculated for each of the vertices in the sphere
Notes
The calculation of ADC, relies on the following relationship:
Where Q is the quadratic form of the tensor.
dipy.reconst.dti.auto_attr(func)Decorator to create OneTimeProperty attributes.
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
dipy.reconst.dti.axial_diffusivity(evals, axis=-1)Axial Diffusivity (AD) of a diffusion tensor. Also called parallel diffusivity.
Eigenvalues of a diffusion tensor, must be sorted in descending order along axis.
Axis of evals which contains 3 eigenvalues.
Calculated AD.
Notes
AD is calculated with the following equation:
dipy.reconst.dti.color_fa(fa, evecs)Color fractional anisotropy of diffusion tensor
Array of the fractional anisotropy (can be 1D, 2D or 3D)
eigen vectors from the tensor model
Colormap of the FA with red for the x value, y for the green value and z for the blue value.
Notes
It is computed from the clipped FA between 0 and 1 using the following formula
dipy.reconst.dti.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).
Hermitian matrix representing a diffusion tensor.
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.
Eigenvalues from eigen decomposition of the tensor. Negative eigenvalues are replaced by zero. Sorted from largest to smallest.
Associated eigenvectors from eigen decomposition of the tensor. Eigenvectors are columnar (e.g. eigvecs[…, :, j] is associated with eigvals[…, j])
dipy.reconst.dti.design_matrix(gtab, dtype=None)Constructs design matrix for DTI weighted least squares or least squares fitting. (Basser et al., 1994a)
Parameter to control the dtype of returned designed matrix
Design matrix or B matrix assuming Gaussian distributed tensor model design_matrix[j, :] = (Bxx, Byy, Bzz, Bxy, Bxz, Byz, dummy)
dipy.reconst.dti.determinant(q_form)The determinant of a tensor, given in quadratic form
The quadratic form of a tensor, or an array with quadratic forms of tensors. Should be of shape (x, y, z, 3, 3) or (n, 3, 3) or (3, 3).
The determinant of the tensor in each spatial coordinate
dipy.reconst.dti.deviatoric(q_form)Calculate the deviatoric (anisotropic) part of the tensor [1].
The quadratic form of a tensor, or an array with quadratic forms of tensors. Should be of shape (x,y,z,3,3) or (n, 3, 3) or (3,3).
The deviatoric part of the tensor in each spatial coordinate.
Notes
The deviatoric part of the tensor is defined as (equations 3-5 in [1]):
Where \(A\) is the tensor quadratic form and \(\bar{A}\) is the anisotropic part of the tensor.
References
dipy.reconst.dti.eig_from_lo_tri(data, min_diffusivity=0)Calculates tensor eigenvalues/eigenvectors from an array containing the lower diagonal form of the six unique tensor elements.
diffusion tensors elements stored in lower triangular order
See decompose_tensor()
Eigen-values and eigen-vectors of the same array.
dipy.reconst.dti.fractional_anisotropy(evals, axis=-1)Return Fractional anisotropy (FA) of a diffusion tensor.
Eigenvalues of a diffusion tensor.
Axis of evals which contains 3 eigenvalues.
Calculated FA. Range is 0 <= FA <= 1.
Notes
FA is calculated using the following equation:
dipy.reconst.dti.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.
Unique elements of the tensors
3 by 3 tensors
dipy.reconst.dti.geodesic_anisotropy(evals, axis=-1)Geodesic anisotropy (GA) of a diffusion tensor.
Eigenvalues of a diffusion tensor.
Axis of evals which contains 3 eigenvalues.
Calculated GA. In the range 0 to +infinity
Notes
GA is calculated using the following equation given in [1]:
Note that the notation, \(<D>\), is often used as the mean diffusivity (MD) of the diffusion tensor and can lead to confusions in the literature (see [1] versus [2] versus [3] for example). Reference [2] defines geodesic anisotropy (GA) with \(<D>\) as the MD in the denominator of the sum. This is wrong. The original paper [1] defines GA with \(<D> = det(D)^{1/3}\), as the isotropic part of the distance. This might be an explanation for the confusion. The isotropic part of the diffusion tensor in Euclidean space is the MD whereas the isotropic part of the tensor in log-Euclidean space is \(det(D)^{1/3}\). The Appendix of [1] and log-Euclidean derivations from [3] are clear on this. Hence, all that to say that \(<D> = det(D)^{1/3}\) here for the GA definition and not MD.
References
P. G. Batchelor, M. Moakher, D. Atkinson, F. Calamante, A. Connelly, “A rigorous framework for diffusion tensor calculus”, Magnetic Resonance in Medicine, vol. 53, pp. 221-225, 2005.
M. M. Correia, V. F. Newcombe, G.B. Williams. “Contrast-to-noise ratios for indices of anisotropy obtained from diffusion MRI: a study with standard clinical b-values at 3T”. NeuroImage, vol. 57, pp. 1103-1115, 2011.
A. D. Lee, etal, P. M. Thompson. “Comparison of fractional and geodesic anisotropy in diffusion tensor images of 90 monozygotic and dizygotic twins”. 5th IEEE International Symposium on Biomedical Imaging (ISBI), pp. 943-946, May 2008.
V. Arsigny, P. Fillard, X. Pennec, N. Ayache. “Log-Euclidean metrics for fast and simple calculus on diffusion tensors.” Magnetic Resonance in Medecine, vol 56, pp. 411-421, 2006.
dipy.reconst.dti.get_sphere(name='symmetric362')provide triangulated spheres
which sphere - one of: * ‘symmetric362’ * ‘symmetric642’ * ‘symmetric724’ * ‘repulsion724’ * ‘repulsion100’ * ‘repulsion200’
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"
dipy.reconst.dti.gradient_table(bvals, bvecs=None, big_delta=None, small_delta=None, b0_threshold=50, atol=0.01, btens=None)A general function for creating diffusion MR gradients.
It reads, loads and prepares scanner parameters like the b-values and b-vectors so that they can be useful during the reconstruction process.
an array of shape (N,) or (1, N) or (N, 1) with the b-values.
a path for the file which contains an array like the above (1).
an array of shape (N, 4) or (4, N). Then this parameter is considered to be a b-table which contains both bvals and bvecs. In this case the next parameter is skipped.
a path for the file which contains an array like the one at (3).
an array of shape (N, 3) or (3, N) with the b-vectors.
a path for the file which contains an array like the previous.
acquisition pulse separation time in seconds (default None)
acquisition pulse duration time in seconds (default None)
All b-values with values less than or equal to bo_threshold are considered as b0s i.e. without diffusion weighting.
All b-vectors need to be unit vectors up to a tolerance.
a string specifying the shape of the encoding tensor for all volumes in data. Options: ‘LTE’, ‘PTE’, ‘STE’, ‘CTE’ corresponding to linear, planar, spherical, and “cigar-shaped” tensor encoding. Tensors are rotated so that linear and cigar tensors are aligned with the corresponding gradient direction and the planar tensor’s normal is aligned with the corresponding gradient direction. Magnitude is scaled to match the b-value.
an array of strings of shape (N,), (N, 1), or (1, N) specifying encoding tensor shape for each volume separately. N corresponds to the number volumes in data. Options for elements in array: ‘LTE’, ‘PTE’, ‘STE’, ‘CTE’ corresponding to linear, planar, spherical, and “cigar-shaped” tensor encoding. Tensors are rotated so that linear and cigar tensors are aligned with the corresponding gradient direction and the planar tensor’s normal is aligned with the corresponding gradient direction. Magnitude is scaled to match the b-value.
an array of shape (N,3,3) specifying the b-tensor of each volume exactly. N corresponds to the number volumes in data. No rotation or scaling is performed.
A GradientTable with all the gradient information.
Notes
Often b0s (b-values which correspond to images without diffusion weighting) have 0 values however in some cases the scanner cannot provide b0s of an exact 0 value and it gives a bit higher values e.g. 6 or 12. This is the purpose of the b0_threshold in the __init__.
We assume that the minimum number of b-values is 7.
B-vectors should be unit vectors.
Examples
>>> from dipy.core.gradients import gradient_table
>>> bvals = 1500 * np.ones(7)
>>> bvals[0] = 0
>>> sq2 = np.sqrt(2) / 2
>>> bvecs = np.array([[0, 0, 0],
... [1, 0, 0],
... [0, 1, 0],
... [0, 0, 1],
... [sq2, sq2, 0],
... [sq2, 0, sq2],
... [0, sq2, sq2]])
>>> gt = gradient_table(bvals, bvecs)
>>> gt.bvecs.shape == bvecs.shape
True
>>> gt = gradient_table(bvals, bvecs.T)
>>> gt.bvecs.shape == bvecs.T.shape
False
dipy.reconst.dti.isotropic(q_form)Calculate the isotropic part of the tensor [1].
The quadratic form of a tensor, or an array with quadratic forms of tensors. Should be of shape (x,y,z,3,3) or (n, 3, 3) or (3,3).
The isotropic part of the tensor in each spatial coordinate
Notes
The isotropic part of a tensor is defined as (equations 3-5 of [1]):
References
dipy.reconst.dti.iter_fit_tensor(step=10000.0)Wrap a fit_tensor func and iterate over chunks of data with given length
Splits data into a number of chunks of specified size and iterates the decorated fit_tensor function over them. This is useful to counteract the temporary but significant memory usage increase in fit_tensor functions that use vectorized operations and need to store large temporary arrays for their vectorized operations.
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.
dipy.reconst.dti.linearity(evals, axis=-1)The linearity of the tensor [1]
Eigenvalues of a diffusion tensor.
Axis of evals which contains 3 eigenvalues.
Calculated linearity of the diffusion tensor.
Notes
Linearity is calculated with the following equation:
References
dipy.reconst.dti.lower_triangular(tensor, b0=None)Returns the six lower triangular values of the tensor and a dummy variable if b0 is not None
a collection of 3, 3 diffusion tensors
if b0 is not none log(b0) is returned as the dummy variable
If b0 is none, then the shape will be (…, 6) otherwise (…, 7)
dipy.reconst.dti.mean_diffusivity(evals, axis=-1)Mean Diffusivity (MD) of a diffusion tensor.
Eigenvalues of a diffusion tensor.
Axis of evals which contains 3 eigenvalues.
Calculated MD.
Notes
MD is calculated with the following equation:
dipy.reconst.dti.mode(q_form)Mode (MO) of a diffusion tensor [1].
The quadratic form of a tensor, or an array with quadratic forms of tensors. Should be of shape (x, y, z, 3, 3) or (n, 3, 3) or (3, 3).
Calculated tensor mode in each spatial coordinate.
Notes
Mode ranges between -1 (planar anisotropy) and +1 (linear anisotropy) with 0 representing orthotropy. Mode is calculated with the following equation (equation 9 in [1]):
Where \(\widetilde{A}\) is the deviatoric part of the tensor quadratic form.
References
dipy.reconst.dti.nlls_fit_tensor(design_matrix, data, weighting=None, sigma=None, jac=True, return_S0_hat=False)Fit the cumulant expansion params (e.g. DTI, DKI) using non-linear least-squares.
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)
Data or response variables holding the data. Note that the last dimension should contain the data. It makes no copies of data.
the weighting scheme to use in considering the squared-error. Default behavior is to use uniform weighting. Other options: ‘sigma’ ‘gmm’
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).
Use the Jacobian? Default: True
Boolean to return (True) or not (False) the S0 values for the fit.
voxel.
dipy.reconst.dti.norm(q_form)Calculate the Frobenius norm of a tensor quadratic form
The quadratic form of a tensor, or an array with quadratic forms of tensors. Should be of shape (x,y,z,3,3) or (n, 3, 3) or (3,3).
The Frobenius norm of the 3,3 tensor q_form in each spatial coordinate.
See also
np.linalg.normNotes
The Frobenius norm is defined as:
||A||_F = [sum_{i,j} abs(a_{i,j})^2]^{1/2}
dipy.reconst.dti.ols_fit_tensor(design_matrix, data, return_S0_hat=False)Computes ordinary least squares (OLS) fit to calculate self-diffusion tensor using a linear regression model [Rd310240b4eed-1].
Design matrix holding the covariants used to solve for the regression coefficients.
Data or response variables holding the data. Note that the last dimension should contain the data. It makes no copies of data.
Boolean to return (True) or not (False) the S0 values for the fit.
Eigenvalues from eigen decomposition of the tensor.
Associated eigenvectors from eigen decomposition of the tensor. Eigenvectors are columnar (e.g. eigvecs[:,j] is associated with eigvals[j])
See also
WLS_fit_tensor, decompose_tensor, design_matrixNotes
References
dipy.reconst.dti.pinv(a, rcond=1e-15)Vectorized version of numpy.linalg.pinv
If numpy version is less than 1.8, it falls back to iterating over np.linalg.pinv since there isn’t a vectorized version of np.linalg.svd available.
Matrix to be pseudo-inverted.
Cutoff for small singular values.
The pseudo-inverse of a.
If the SVD computation does not converge.
See also
np.linalg.pinvdipy.reconst.dti.planarity(evals, axis=-1)The planarity of the tensor [1]
Eigenvalues of a diffusion tensor.
Axis of evals which contains 3 eigenvalues.
Calculated linearity of the diffusion tensor.
Notes
Planarity is calculated with the following equation:
References
dipy.reconst.dti.radial_diffusivity(evals, axis=-1)Radial Diffusivity (RD) of a diffusion tensor. Also called perpendicular diffusivity.
Eigenvalues of a diffusion tensor, must be sorted in descending order along axis.
Axis of evals which contains 3 eigenvalues.
Calculated RD.
Notes
RD is calculated with the following equation:
dipy.reconst.dti.restore_fit_tensor(design_matrix, data, sigma=None, jac=True, return_S0_hat=False)Use the RESTORE algorithm [1] to calculate a robust tensor fit
Design matrix holding the covariants used to solve for the regression coefficients.
Data or response variables holding the data. Note that the last dimension should contain the data. It makes no copies of data.
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).
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
Boolean to return (True) or not (False) the S0 values for the fit.
References
estimation of tensors by outlier rejection. MRM, 53: 1088-95.
dipy.reconst.dti.sphericity(evals, axis=-1)The sphericity of the tensor [1]
Eigenvalues of a diffusion tensor.
Axis of evals which contains 3 eigenvalues.
Calculated sphericity of the diffusion tensor.
Notes
Sphericity is calculated with the following equation:
References
dipy.reconst.dti.tensor_prediction(dti_params, gtab, S0)Predict a signal given tensor parameters.
Tensor parameters. The last dimension should have 12 tensor parameters: 3 eigenvalues, followed by the 3 corresponding eigenvectors.
The gradient table for this prediction
The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1
Notes
The predicted signal is given by: \(S( heta, b) = S_0 * e^{-b ADC}\), where \(ADC = heta Q heta^T\), :math:` heta` is a unit vector pointing at any direction on the sphere for which a signal is to be predicted, \(b\) is the b value provided in the GradientTable input for that direction, \(Q\) is the quadratic form of the tensor determined by the input parameters.
dipy.reconst.dti.trace(evals, axis=-1)Trace of a diffusion tensor.
Eigenvalues of a diffusion tensor.
Axis of evals which contains 3 eigenvalues.
Calculated trace of the diffusion tensor.
Notes
Trace is calculated with the following equation:
dipy.reconst.dti.vec_val_vect()Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs
containing tensor in last two dimensions; M, N usually equal to (3, 3)
diagonal values carried in last dimension, ... shape above must
match that for vecs
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).
... dimensions of vecs, valsN dimensions of vecs, valsExamples
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]])
dipy.reconst.dti.vector_norm(vec, axis=-1, keepdims=False)Return vector Euclidean (L2) norm
See unit vector and Euclidean norm
Vectors to norm.
Axis over which to norm. By default norm over last axis. If axis is None, vec is flattened then normed.
If True, the output will have the same number of dimensions as vec, with shape 1 on axis.
Euclidean norms of vectors.
Examples
>>> import numpy as np
>>> vec = [[8, 15, 0], [0, 36, 77]]
>>> vector_norm(vec)
array([ 17., 85.])
>>> vector_norm(vec, keepdims=True)
array([[ 17.],
[ 85.]])
>>> vector_norm(vec, axis=0)
array([ 8., 39., 77.])
dipy.reconst.dti.wls_fit_tensor(design_matrix, data, return_S0_hat=False)Computes weighted least squares (WLS) fit to calculate self-diffusion tensor using a linear regression model [1].
Design matrix holding the covariants used to solve for the regression coefficients.
Data or response variables holding the data. Note that the last dimension should contain the data. It makes no copies of data.
Boolean to return (True) or not (False) the S0 values for the fit.
Eigenvalues from eigen decomposition of the tensor.
Associated eigenvectors from eigen decomposition of the tensor. Eigenvectors are columnar (e.g. eigvecs[:,j] is associated with eigvals[j])
See also
Notes
In Chung, et al. 2006, the regression of the WLS fit needed an unbiased preliminary estimate of the weights and therefore the ordinary least squares (OLS) estimates were used. A “two pass” method was implemented:
calculate OLS estimates of the data
apply the OLS estimates as weights to the WLS fit of the data
This ensured heteroscedasticity could be properly modeled for various types of bootstrap resampling (namely residual bootstrap).
References
Cachedipy.reconst.forecast.CacheBases: 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
Clear the cache. |
|
|
Retrieve a value from the cache. |
|
Store a value in the cache. |
cache_get(tag, key, default=None)Retrieve a value from the cache.
Description of the cached value.
Key object used to look up the cached value.
Value to be returned if no cached entry is found.
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.
Description of the cached value.
Key object used to look up the cached value.
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)
>>> c = Cache()
>>> parameters = (1, 2, 3)
>>> X1 = compute_expensive_matrix(parameters)
>>> c.cache_set('expensive_matrix', parameters, X1)
>>> X2 = c.cache_get('expensive_matrix', parameters)
>>> X1 is X2
True
ForecastFitdipy.reconst.forecast.ForecastFit(model, data, sh_coef, d_par, d_perp)Bases: dipy.reconst.odf.OdfFit
Methods
Calculates the fractional anisotropy. |
|
Calculates the mean diffusivity. |
|
|
Calculates the fODF for a given discrete sphere. |
|
Calculates the fODF for a given discrete sphere. |
__init__(model, data, sh_coef, d_par, d_perp)Calculates diffusion properties for a single voxel
AnalyticalModel
fitted data
forecast sh coefficients
parallel diffusivity
perpendicular diffusivity
odf(sphere, clip_negative=True)Calculates the fODF for a given discrete sphere.
the odf sphere
if True clip the negative odf values to 0, default True
ForecastModeldipy.reconst.forecast.ForecastModel(gtab, sh_order=8, lambda_lb=0.001, dec_alg='CSD', sphere=None, lambda_csd=1.0)Bases: dipy.reconst.odf.OdfModel, dipy.reconst.cache.Cache
Fiber ORientation Estimated using Continuous Axially Symmetric Tensors (FORECAST) [1,2,3]_. FORECAST is a Spherical Deconvolution reconstruction model for multi-shell diffusion data which enables the calculation of a voxel adaptive response function using the Spherical Mean Tecnique (SMT) [2,3]_.
With FORECAST it is possible to calculate crossing invariant parallel diffusivity, perpendicular diffusivity, mean diffusivity, and fractional anisotropy [2]
Notes
The implementation of FORECAST may require CVXPY (http://www.cvxpy.org/).
References
Anderson A. W., “Measurement of Fiber Orientation Distributions Using High Angular Resolution Diffusion Imaging”, Magnetic Resonance in Medicine, 2005.
Kaden E. et al., “Quantitative Mapping of the Per-Axon Diffusion Coefficients in Brain White Matter”, Magnetic Resonance in Medicine, 2016.
Zucchelli E. et al., “A generalized SMT-based framework for Diffusion MRI microstructural model estimation”, MICCAI Workshop on Computational DIFFUSION MRI (CDMRI), 2017.
Methods
|
Clear the cache. |
|
Retrieve a value from the cache. |
|
Store a value in the cache. |
|
Fit method for every voxel in data |
__init__(gtab, sh_order=8, lambda_lb=0.001, dec_alg='CSD', sphere=None, lambda_csd=1.0)Analytical and continuous modeling of the diffusion signal with respect to the FORECAST basis [1,2,3]_. This implementation is a modification of the original FORECAST model presented in [1] adapted for multi-shell data as in [2,3]_ .
The main idea is to model the diffusion signal as the combination of a single fiber response function \(F(\mathbf{b})\) times the fODF \(\rho(\mathbf{v})\)
E(mathbf{b}) = int_{mathbf{v} in mathcal{S}^2} rho(mathbf{v}) F({mathbf{b}} | mathbf{v}) d mathbf{v}
end{equation}
where \(\mathbf{b}\) is the b-vector (b-value times gradient direction) and \(\mathbf{v}\) is an unit vector representing a fiber direction.
In FORECAST \(\rho\) is modeled using real symmetric Spherical Harmonics (SH) and \(F(\mathbf(b))\) is an axially symmetric tensor.
gradient directions and bvalues container class.
an even integer that represent the SH order of the basis (max 12)
Laplace-Beltrami regularization weight.
Spherical deconvolution algorithm. The possible values are Weighted Least Squares (‘WLS’), Positivity Constraints using CVXPY (‘POS’) and the Constraint Spherical Deconvolution algorithm (‘CSD’). Default is ‘CSD’.
sphere points where to enforce positivity when ‘POS’ or ‘CSD’ dec_alg are selected.
CSD regularization weight.
References
Anderson A. W., “Measurement of Fiber Orientation Distributions Using High Angular Resolution Diffusion Imaging”, Magnetic Resonance in Medicine, 2005.
Kaden E. et al., “Quantitative Mapping of the Per-Axon Diffusion Coefficients in Brain White Matter”, Magnetic Resonance in Medicine, 2016.
Zucchelli M. et al., “A generalized SMT-based framework for Diffusion MRI microstructural model estimation”, MICCAI Workshop on Computational DIFFUSION MRI (CDMRI), 2017.
Examples
In this example, where the data, gradient table and sphere tessellation used for reconstruction are provided, we model the diffusion signal with respect to the FORECAST and compute the fODF, parallel and perpendicular diffusivity.
>>> from dipy.data import default_sphere, get_3shell_gtab
>>> gtab = get_3shell_gtab()
>>> from dipy.sims.voxel import multi_tensor
>>> mevals = np.array(([0.0017, 0.0003, 0.0003],
... [0.0017, 0.0003, 0.0003]))
>>> angl = [(0, 0), (60, 0)]
>>> data, sticks = multi_tensor(gtab,
... mevals,
... S0=100.0,
... angles=angl,
... fractions=[50, 50],
... snr=None)
>>> from dipy.reconst.forecast import ForecastModel
>>> fm = ForecastModel(gtab, sh_order=6)
>>> f_fit = fm.fit(data)
>>> d_par = f_fit.dpar
>>> d_perp = f_fit.dperp
>>> fodf = f_fit.odf(default_sphere)
LooseVersiondipy.reconst.forecast.LooseVersion(vstring=None)Bases: distutils.version.Version
Version numbering for anarchists and software realists. Implements the standard interface for version number classes as described above. A version number consists of a series of numbers, separated by either periods or strings of letters. When comparing version numbers, the numeric components will be compared numerically, and the alphabetic components lexically. The following are all valid version numbers, in no particular order:
1.5.1 1.5.2b2 161 3.10a 8.02 3.4j 1996.07.12 3.2.pl0 3.1.1.6 2g6 11g 0.960923 2.2beta29 1.13++ 5.5.kw 2.0b1pl0
In fact, there is no such thing as an invalid version number under this scheme; the rules for comparison are simple and predictable, but may not always give the results you want (for some definition of “want”).
Methods
parse |
OdfFitdipy.reconst.forecast.OdfFit(model, data)Bases: dipy.reconst.base.ReconstFit
Methods
|
To be implemented but specific odf models |
OdfModeldipy.reconst.forecast.OdfModel(gtab)Bases: dipy.reconst.base.ReconstModel
An abstract class to be sub-classed by specific odf models
All odf models should provide a fit method which may take data as it’s first and only argument.
Methods
|
To be implemented by specific odf models |
dipy.reconst.forecast.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\)
x coordinate in Cartesian space
y coordinate in Cartesian space
z coordinate
radius
inclination (polar) angle
azimuth angle
dipy.reconst.forecast.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.
Diffusion weighted signals to be deconvolved.
Prediction matrix which estimates diffusion weighted signals from FOD coefficients.
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.
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.
Maximum number of iterations to allow the deconvolution to converge.
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.
(sh_order + 1)*(sh_order + 2)/2,)Spherical harmonics coefficients of the constrained-regularized fiber ODF.
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