# reconst

 bench([label, verbose, extra_argv]) Run benchmarks for module using nose. test([label, verbose, extra_argv, doctests, …]) Run tests for module using nose.

## Module: reconst.base

Base-classes for reconstruction models and reconstruction fits.

All the models in the reconst module follow the same template: a Model object is used to represent the abstract properties of the model, that are independent of the specifics of the data . These properties are reused whenver fitting a particular set of data (different voxels, for example).

 ReconstFit(model, data) Abstract class which holds the fit result of ReconstModel ReconstModel(gtab) Abstract class for signal reconstruction models

## Module: reconst.benchmarks.bench_bounding_box

 Compute the bounding box of nonzero intensity voxels in the volume. measure(code_str[, times, label]) Return elapsed time for executing code in the namespace of the caller.

## Module: reconst.benchmarks.bench_csd

 ConstrainedSphericalDeconvModel(gtab, response) Methods GradientTable(gradients[, big_delta, …]) Diffusion gradient information bench_csdeconv([center, width]) num_grad(gtab) Read stanford hardi data and label map

## Module: reconst.benchmarks.bench_peaks

Benchmarks for peak finding

Run all benchmarks with:

import dipy.reconst as dire
dire.bench()


With Pytest, Run this benchmark with:

pytest -svv -c bench.ini /path/to/bench_peaks.py

 local_maxima Local maxima of a function evaluated on a discrete set of points. measure(code_str[, times, label]) Return elapsed time for executing code in the namespace of the caller. unique_edges(faces[, return_mapping]) Extract all unique edges from given triangular faces.

## Module: reconst.benchmarks.bench_squash

Benchmarks for fast squashing

Run all benchmarks with:

import dipy.reconst as dire
dire.bench()


With Pytest, Run this benchmark with:

pytest -svv -c bench.ini /path/to/bench_squash.py

 measure(code_str[, times, label]) Return elapsed time for executing code in the namespace of the caller. ndindex(shape) An N-dimensional iterator object to index arrays. old_squash(arr[, mask, fill]) Try and make a standard array from an object array quick_squash Try and make a standard array from an object array reduce(function, sequence[, initial]) Apply a function of two arguments cumulatively to the items of a sequence, from left to right, so as to reduce the sequence to a single value.

## Module: reconst.benchmarks.bench_vec_val_sum

Benchmarks for vec / val summation routine

Run benchmarks with:

import dipy.reconst as dire
dire.bench()


With Pytest, Run this benchmark with:

pytest -svv -c bench.ini /path/to/bench_vec_val_sum.py

 measure(code_str[, times, label]) Return elapsed time for executing code in the namespace of the caller. randn(d0, d1, …, dn) Return a sample (or samples) from the “standard normal” distribution. vec_val_vect Vectorize vecs.diag(vals).vecs.T for last 2 dimensions of vecs

## Module: reconst.cache

 Cache Cache values based on a key object (such as a sphere or gradient table). auto_attr(func) Decorator to create OneTimeProperty attributes.

## Module: reconst.cross_validation

Cross-validation analysis of diffusion models

 coeff_of_determination(data, model[, axis]) Calculate the coefficient of determination for a model prediction, relative kfold_xval(model, data, folds, *model_args, …) Perform k-fold cross-validation to generate out-of-sample predictions for each measurement.

## Module: reconst.csdeconv

 AxSymShResponse(S0, dwi_response[, bvalue]) A simple wrapper for response functions represented using only axially symmetric, even spherical harmonic functions (ie, m == 0 and n even). ConstrainedSDTModel(gtab, ratio[, …]) Methods ConstrainedSphericalDeconvModel(gtab, response) Methods SphHarmFit(model, shm_coef, mask) Diffusion data fit to a spherical harmonic model SphHarmModel(gtab) To be subclassed by all models that return a SphHarmFit when fit. TensorModel(gtab[, fit_method, return_S0_hat]) Diffusion Tensor auto_response(gtab, data[, roi_center, …]) Automatic estimation of response function using FA. cart2sphere(x, y, z) Return angles for Cartesian 3D coordinates x, y, and z csdeconv(dwsignal, X, B_reg[, tau, …]) Constrained-regularized spherical deconvolution (CSD) [1] estimate_response(gtab, evals, S0) Estimate single fiber response function fa_inferior(FA, fa_thr) Check that the FA is lower than the FA threshold fa_superior(FA, fa_thr) Check that the FA is greater than the FA threshold fa_trace_to_lambdas([fa, trace]) forward_sdeconv_mat(r_rh, n) Build forward spherical deconvolution matrix forward_sdt_deconv_mat(ratio, n[, r2_term]) Build forward sharpening deconvolution transform (SDT) matrix fractional_anisotropy(evals[, axis]) Fractional anisotropy (FA) of a diffusion tensor. get_sphere([name]) provide triangulated spheres lazy_index(index) Produces a lazy index lpn(n, z) Legendre function of the first kind. multi_voxel_fit(single_voxel_fit) Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition ndindex(shape) An N-dimensional iterator object to index arrays. odf_deconv(odf_sh, R, B_reg[, lambda_, tau, …]) ODF constrained-regularized spherical deconvolution using the Sharpening Deconvolution Transform (SDT) [1], [2]. odf_sh_to_sharp(odfs_sh, sphere[, basis, …]) Sharpen odfs using the sharpening deconvolution transform [2] peaks_from_model(model, data, sphere, …[, …]) Fit the model to data and computes peaks and metrics quad(func, a, b[, args, full_output, …]) Compute a definite integral. real_sph_harm(m, n, theta, phi) Compute real spherical harmonics. real_sym_sh_basis(sh_order, theta, phi) Samples a real symmetric spherical harmonic basis at point on the sphere recursive_response(gtab, data[, mask, …]) Recursive calibration of response function using peak threshold response_from_mask(gtab, data, mask) Estimate the response function from a given mask. sh_to_rh(r_sh, m, n) Spherical harmonics (SH) to rotational harmonics (RH) single_tensor(gtab[, S0, evals, evecs, snr]) Simulated Q-space signal with a single tensor. sph_harm_ind_list(sh_order) Returns the degree (n) and order (m) of all the symmetric spherical harmonics of degree less then or equal to sh_order. vec2vec_rotmat(u, v) rotation matrix from 2 unit vectors

## Module: reconst.dki

Classes and functions for fitting the diffusion kurtosis model

 DiffusionKurtosisFit(model, model_params) Class for fitting the Diffusion Kurtosis Model DiffusionKurtosisModel(gtab[, fit_method]) Class for the Diffusion Kurtosis Model ReconstModel(gtab) Abstract class for signal reconstruction models TensorFit(model, model_params[, model_S0]) Attributes Wcons(k_elements) Construct the full 4D kurtosis tensors from its 15 independent elements Wrotate(kt, Basis) Rotate a kurtosis tensor from the standard Cartesian coordinate system to another coordinate system basis Wrotate_element(kt, indi, indj, indk, indl, B) Computes the the specified index element of a kurtosis tensor rotated to the coordinate system basis B. apparent_kurtosis_coef(dki_params, sphere[, …]) Calculates the apparent kurtosis coefficient (AKC) in each direction of a sphere [1]. axial_kurtosis(dki_params[, min_kurtosis, …]) Computes axial Kurtosis (AK) from the kurtosis tensor. carlson_rd(x, y, z[, errtol]) Computes the Carlson’s incomplete elliptic integral of the second kind defined as: carlson_rf(x, y, z[, errtol]) Computes the Carlson’s incomplete elliptic integral of the first kind defined as: cart2sphere(x, y, z) Return angles for Cartesian 3D coordinates x, y, and z check_multi_b(gtab, n_bvals[, non_zero, bmag]) Check if you have enough different b-values in your gradient table decompose_tensor(tensor[, min_diffusivity]) Returns eigenvalues and eigenvectors given a diffusion tensor design_matrix(gtab) Construct B design matrix for DKI. directional_diffusion(dt, V[, min_diffusivity]) Calculates the apparent diffusion coefficient (adc) in each direction of a sphere for a single voxel [1]. directional_diffusion_variance(kt, V[, …]) Calculates the apparent diffusion variance (adv) in each direction of a sphere for a single voxel [1]. directional_kurtosis(dt, md, kt, V[, …]) Calculates the apparent kurtosis coefficient (akc) in each direction of a sphere for a single voxel [1]. dki_prediction(dki_params, gtab[, S0]) Predict a signal given diffusion kurtosis imaging parameters. Returns a tensor given the six unique tensor elements get_sphere([name]) provide triangulated spheres kurtosis_maximum(dki_params[, sphere, gtol, …]) Computes kurtosis maximum value local_maxima Local maxima of a function evaluated on a discrete set of points. lower_triangular(tensor[, b0]) Returns the six lower triangular values of the tensor and a dummy variable if b0 is not None mean_diffusivity(evals[, axis]) Mean Diffusivity (MD) of a diffusion tensor. mean_kurtosis(dki_params[, min_kurtosis, …]) Computes mean Kurtosis (MK) from the kurtosis tensor [1]. ndindex(shape) An N-dimensional iterator object to index arrays. ols_fit_dki(design_matrix, data) Computes ordinary least squares (OLS) fit to calculate the diffusion tensor and kurtosis tensor using a linear regression diffusion kurtosis model 1. radial_kurtosis(dki_params[, min_kurtosis, …]) Radial Kurtosis (RK) of a diffusion kurtosis tensor [1]. sphere2cart(r, theta, phi) Spherical to Cartesian coordinates split_dki_param(dki_params) Extract the diffusion tensor eigenvalues, the diffusion tensor eigenvector matrix, and the 15 independent elements of the kurtosis tensor from the model parameters estimated from the DKI model vec_val_vect Vectorize vecs.diag(vals).vecs.T for last 2 dimensions of vecs wls_fit_dki(design_matrix, data) Computes weighted linear least squares (WLS) fit to calculate the diffusion tensor and kurtosis tensor using a weighted linear regression diffusion kurtosis model 1.

## Module: reconst.dki_micro

Classes and functions for fitting the DKI-based microstructural model

 DiffusionKurtosisFit(model, model_params) Class for fitting the Diffusion Kurtosis Model DiffusionKurtosisModel(gtab[, fit_method]) Class for the Diffusion Kurtosis Model KurtosisMicrostructuralFit(model, model_params) Class for fitting the Diffusion Kurtosis Microstructural Model KurtosisMicrostructureModel(gtab[, fit_method]) Class for the Diffusion Kurtosis Microstructural Model axial_diffusivity(evals[, axis]) Axial Diffusivity (AD) of a diffusion tensor. axonal_water_fraction(dki_params[, sphere, …]) Computes the axonal water fraction from DKI [1]. decompose_tensor(tensor[, min_diffusivity]) Returns eigenvalues and eigenvectors given a diffusion tensor diffusion_components(dki_params[, sphere, …]) Extracts the restricted and hindered diffusion tensors of well aligned fibers from diffusion kurtosis imaging parameters [1]. directional_diffusion(dt, V[, min_diffusivity]) Calculates the apparent diffusion coefficient (adc) in each direction of a sphere for a single voxel [1]. directional_kurtosis(dt, md, kt, V[, …]) Calculates the apparent kurtosis coefficient (akc) in each direction of a sphere for a single voxel [1]. dkimicro_prediction(params, gtab[, S0]) Signal prediction given the DKI microstructure model parameters. dti_design_matrix(gtab[, dtype]) Constructs design matrix for DTI weighted least squares or least squares fitting. Returns a tensor given the six unique tensor elements get_sphere([name]) provide triangulated spheres kurtosis_maximum(dki_params[, sphere, gtol, …]) Computes kurtosis maximum value lower_triangular(tensor[, b0]) Returns the six lower triangular values of the tensor and a dummy variable if b0 is not None mean_diffusivity(evals[, axis]) Mean Diffusivity (MD) of a diffusion tensor. ndindex(shape) An N-dimensional iterator object to index arrays. radial_diffusivity(evals[, axis]) Radial Diffusivity (RD) of a diffusion tensor. split_dki_param(dki_params) Extract the diffusion tensor eigenvalues, the diffusion tensor eigenvector matrix, and the 15 independent elements of the kurtosis tensor from the model parameters estimated from the DKI model tortuosity(hindered_ad, hindered_rd) Computes the tortuosity of the hindered diffusion compartment given its axial and radial diffusivities trace(evals[, axis]) Trace of a diffusion tensor. vec_val_vect Vectorize vecs.diag(vals).vecs.T for last 2 dimensions of vecs

## Module: reconst.dsi

 Cache Cache values based on a key object (such as a sphere or gradient table). DiffusionSpectrumDeconvFit(model, data) Methods DiffusionSpectrumDeconvModel(gtab[, …]) Methods DiffusionSpectrumFit(model, data) Methods DiffusionSpectrumModel(gtab[, qgrid_size, …]) Methods OdfFit(model, data) Methods OdfModel(gtab) An abstract class to be sub-classed by specific odf models LR_deconv(prop, psf[, numit, acc_factor]) Perform Lucy-Richardson deconvolution algorithm on a 3D array. create_qspace(gtab, origin) create the 3D grid which holds the signal values (q-space) create_qtable(gtab, origin) create a normalized version of gradients fftn(x[, shape, axes, overwrite_x]) Return multidimensional discrete Fourier transform. fftshift(x[, axes]) Shift the zero-frequency component to the center of the spectrum. gen_PSF(qgrid_sampling, siz_x, siz_y, siz_z) Generate a PSF for DSI Deconvolution by taking the ifft of the binary q-space sampling mask and truncating it to keep only the center. half_to_full_qspace(data, gtab) Half to full Cartesian grid mapping hanning_filter(gtab, filter_width, origin) create a hanning window ifftshift(x[, axes]) The inverse of fftshift. map_coordinates(input, coordinates[, …]) Map the input array to new coordinates by interpolation. multi_voxel_fit(single_voxel_fit) Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition pdf_interp_coords(sphere, rradius, origin) Precompute coordinates for ODF calculation from the PDF pdf_odf(Pr, rradius, interp_coords) Calculates the real ODF from the diffusion propagator(PDF) Pr Project any near identical bvecs to the other hemisphere threshold_propagator(P[, estimated_snr]) Applies hard threshold on the propagator to remove background noise for the deconvolution.

## Module: reconst.dti

Classes and functions for fitting tensors

 ReconstModel(gtab) Abstract class for signal reconstruction models TensorFit(model, model_params[, model_S0]) Attributes TensorModel(gtab[, fit_method, return_S0_hat]) Diffusion Tensor apparent_diffusion_coef(q_form, sphere) Calculate the apparent diffusion coefficient (ADC) in each direction of a auto_attr(func) Decorator to create OneTimeProperty attributes. axial_diffusivity(evals[, axis]) Axial Diffusivity (AD) of a diffusion tensor. color_fa(fa, evecs) Color fractional anisotropy of diffusion tensor decompose_tensor(tensor[, min_diffusivity]) Returns eigenvalues and eigenvectors given a diffusion tensor design_matrix(gtab[, dtype]) Constructs design matrix for DTI weighted least squares or least squares fitting. determinant(q_form) The determinant of a tensor, given in quadratic form deviatoric(q_form) Calculate the deviatoric (anisotropic) part of the tensor [1]. eig_from_lo_tri(data[, min_diffusivity]) Calculates tensor eigenvalues/eigenvectors from an array containing the lower diagonal form of the six unique tensor elements. fractional_anisotropy(evals[, axis]) Fractional anisotropy (FA) of a diffusion tensor. Returns a tensor given the six unique tensor elements geodesic_anisotropy(evals[, axis]) Geodesic anisotropy (GA) of a diffusion tensor. get_sphere([name]) provide triangulated spheres gradient_table(bvals[, bvecs, big_delta, …]) A general function for creating diffusion MR gradients. isotropic(q_form) Calculate the isotropic part of the tensor [Rd0568a744381-1]. iter_fit_tensor([step]) Wrap a fit_tensor func and iterate over chunks of data with given length linearity(evals[, axis]) The linearity of the tensor 1 lower_triangular(tensor[, b0]) Returns the six lower triangular values of the tensor and a dummy variable if b0 is not None mean_diffusivity(evals[, axis]) Mean Diffusivity (MD) of a diffusion tensor. mode(q_form) Mode (MO) of a diffusion tensor [1]. nlls_fit_tensor(design_matrix, data[, …]) Fit the tensor params using non-linear least-squares. norm(q_form) Calculate the Frobenius norm of a tensor quadratic form ols_fit_tensor(design_matrix, data[, …]) Computes ordinary least squares (OLS) fit to calculate self-diffusion tensor using a linear regression model 1. pinv(a[, rcond]) Vectorized version of numpy.linalg.pinv planarity(evals[, axis]) The planarity of the tensor 1 quantize_evecs(evecs[, odf_vertices]) Find the closest orientation of an evenly distributed sphere radial_diffusivity(evals[, axis]) Radial Diffusivity (RD) of a diffusion tensor. restore_fit_tensor(design_matrix, data[, …]) Use the RESTORE algorithm [Chang2005] to calculate a robust tensor fit sphericity(evals[, axis]) The sphericity of the tensor 1 tensor_prediction(dti_params, gtab, S0) Predict a signal given tensor parameters. trace(evals[, axis]) Trace of a diffusion tensor. vec_val_vect Vectorize vecs.diag(vals).vecs.T for last 2 dimensions of vecs vector_norm(vec[, axis, keepdims]) Return vector Euclidean (L2) norm wls_fit_tensor(design_matrix, data[, …]) Computes weighted least squares (WLS) fit to calculate self-diffusion tensor using a linear regression model [1].

## Module: reconst.forecast

 Cache Cache values based on a key object (such as a sphere or gradient table). ForecastFit(model, data, sh_coef, d_par, d_perp) Attributes ForecastModel(gtab[, sh_order, lambda_lb, …]) Fiber ORientation Estimated using Continuous Axially Symmetric Tensors (FORECAST) [1,2,3]_. OdfFit(model, data) Methods OdfModel(gtab) An abstract class to be sub-classed by specific odf models cart2sphere(x, y, z) Return angles for Cartesian 3D coordinates x, y, and z csdeconv(dwsignal, X, B_reg[, tau, …]) Constrained-regularized spherical deconvolution (CSD) [1] find_signal_means(b_unique, data_norm, …) Calculate the mean signal for each shell. forecast_error_func(x, b_unique, E) Calculates the difference between the mean signal calculated using the parameter vector x and the average signal E using FORECAST and SMT forecast_matrix(sh_order, d_par, d_perp, bvals) Compute the FORECAST radial matrix lb_forecast(sh_order) Returns the Laplace-Beltrami regularization matrix for FORECAST leastsq(func, x0[, args, Dfun, full_output, …]) Minimize the sum of squares of a set of equations. multi_voxel_fit(single_voxel_fit) Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition optional_package(name[, trip_msg]) Return package-like thing and module setup for package name psi_l(l, b) real_sph_harm(m, n, theta, phi) Compute real spherical harmonics. rho_matrix(sh_order, vecs) Compute the SH matrix $$\rho$$ warn Issue a warning, or maybe ignore it or raise an exception.

## Module: reconst.fwdti

Classes and functions for fitting tensors without free water contamination

 FreeWaterTensorFit(model, model_params) Class for fitting the Free Water Tensor Model FreeWaterTensorModel(gtab[, fit_method]) Class for the Free Water Elimination Diffusion Tensor Model ReconstModel(gtab) Abstract class for signal reconstruction models TensorFit(model, model_params[, model_S0]) Attributes check_multi_b(gtab, n_bvals[, non_zero, bmag]) Check if you have enough different b-values in your gradient table Convert Cholesky decompostion elements to the diffusion tensor elements decompose_tensor(tensor[, min_diffusivity]) Returns eigenvalues and eigenvectors given a diffusion tensor design_matrix(gtab[, dtype]) Constructs design matrix for DTI weighted least squares or least squares fitting. Returns a tensor given the six unique tensor elements fwdti_prediction(params, gtab[, S0, Diso]) Signal prediction given the free water DTI model parameters. lower_triangular(tensor[, b0]) Returns the six lower triangular values of the tensor and a dummy variable if b0 is not None lower_triangular_to_cholesky(tensor_elements) Performs Cholesky decomposition of the diffusion tensor multi_voxel_fit(single_voxel_fit) Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition ndindex(shape) An N-dimensional iterator object to index arrays. nls_fit_tensor(gtab, data[, mask, Diso, …]) Fit the water elimination tensor model using the non-linear least-squares. nls_iter(design_matrix, sig, S0[, Diso, …]) Applies non linear least squares fit of the water free elimination model to single voxel signals. vec_val_vect Vectorize vecs.diag(vals).vecs.T for last 2 dimensions of vecs wls_fit_tensor(gtab, data[, Diso, mask, …]) Computes weighted least squares (WLS) fit to calculate self-diffusion tensor using a linear regression model [1]. wls_iter(design_matrix, sig, S0[, Diso, …]) Applies weighted linear least squares fit of the water free elimination model to single voxel signals.

## Module: reconst.gqi

Classes and functions for generalized q-sampling

 Cache Cache values based on a key object (such as a sphere or gradient table). GeneralizedQSamplingFit(model, data) Methods GeneralizedQSamplingModel(gtab[, method, …]) Methods OdfFit(model, data) Methods OdfModel(gtab) An abstract class to be sub-classed by specific odf models equatorial_maximum(vertices, odf, pole, width) equatorial_zone_vertices(vertices, pole[, width]) finds the ‘vertices’ in the equatorial zone conjugate to ‘pole’ with width half ‘width’ degrees gfa(samples) The general fractional anisotropy of a function evaluated on the unit sphere local_maxima Local maxima of a function evaluated on a discrete set of points. multi_voxel_fit(single_voxel_fit) Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition normalize_qa(qa[, max_qa]) Normalize quantitative anisotropy. npa(self, odf[, width]) non-parametric anisotropy odf_sum(odf) patch_maximum(vertices, odf, pole, width) patch_sum(vertices, odf, pole, width) patch_vertices(vertices, pole, width) find ‘vertices’ within the cone of ‘width’ degrees around ‘pole’ polar_zone_vertices(vertices, pole[, width]) finds the ‘vertices’ in the equatorial band around the ‘pole’ of radius ‘width’ degrees remove_similar_vertices Remove vertices that are less than theta degrees from any other squared_radial_component(x[, tol]) Part of the GQI2 integral triple_odf_maxima(vertices, odf, width) maps a 3-vector into the z-upper hemisphere

## Module: reconst.ivim

Classes and functions for fitting ivim model

 IvimFit(model, model_params) Attributes IvimModelTRR(gtab[, split_b_D, split_b_S0, …]) Ivim model IvimModelVP(gtab[, bounds, maxiter, xtol]) Methods ReconstModel(gtab) Abstract class for signal reconstruction models IvimModel(gtab[, fit_method]) Selector function to switch between the 2-stage Trust-Region Reflective based NLLS fitting method (also containing the linear fit): trr and the Variable Projections based fitting method: varpro. differential_evolution(func, bounds[, args, …]) Finds the global minimum of a multivariate function. f_D_star_error(params, gtab, signal, S0, D) Error function used to fit f and D_star keeping S0 and D fixed f_D_star_prediction(params, gtab, S0, D) Function used to predict IVIM signal when S0 and D are known by considering f and D_star as the unknown parameters. ivim_model_selector(gtab[, fit_method]) Selector function to switch between the 2-stage Trust-Region Reflective based NLLS fitting method (also containing the linear fit): trr and the Variable Projections based fitting method: varpro. ivim_prediction(params, gtab) The Intravoxel incoherent motion (IVIM) model function. least_squares(fun, x0[, jac, bounds, …]) Solve a nonlinear least-squares problem with bounds on the variables. multi_voxel_fit(single_voxel_fit) Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition optional_package(name[, trip_msg]) Return package-like thing and module setup for package name

## Module: reconst.mapmri

 Cache Cache values based on a key object (such as a sphere or gradient table). MapmriFit(model, mapmri_coef, mu, R, lopt[, …]) Attributes MapmriModel(gtab[, radial_order, …]) Mean Apparent Propagator MRI (MAPMRI) [1] of the diffusion signal. Optimizer(fun, x0[, args, method, jac, …]) Attributes ReconstFit(model, data) Abstract class which holds the fit result of ReconstModel ReconstModel(gtab) Abstract class for signal reconstruction models b_mat(index_matrix) Calculates the B coefficients from [1] Eq. b_mat_isotropic(index_matrix) Calculates the isotropic B coefficients from [1] Fig 8. binomialfloat(n, k) Custom Binomial function cart2sphere(x, y, z) Return angles for Cartesian 3D coordinates x, y, and z create_rspace(gridsize, radius_max) Create the real space table, that contains the points in which to compute the pdf. delta(n, m) factorial2(n[, exact]) Double factorial. gcv_cost_function(weight, args) The GCV cost function that is iterated [4]. generalized_crossvalidation(data, M, LR[, …]) Generalized Cross Validation Function [Rb690cd738504-1] eq. generalized_crossvalidation_array(data, M, LR) Generalized Cross Validation Function 1 eq. genlaguerre(n, alpha[, monic]) Generalized (associated) Laguerre polynomial. gradient_table(bvals[, bvecs, big_delta, …]) A general function for creating diffusion MR gradients. hermite(n[, monic]) Physicist’s Hermite polynomial. isotropic_scale_factor(mu_squared) Estimated isotropic scaling factor _[1] Eq. map_laplace_s(n, m) R(m,n) static matrix for Laplacian regularization [R932dd40ca52e-1] eq. map_laplace_t(n, m) L(m, n) static matrix for Laplacian regularization [Reb78d789d6c4-1] eq. map_laplace_u(n, m) S(n, m) static matrix for Laplacian regularization [Rb93dd9dab8c9-1] eq. mapmri_STU_reg_matrices(radial_order) Generate the static portions of the Laplacian regularization matrix according to [R1d585103467a-1] eq. mapmri_index_matrix(radial_order) 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. mapmri_isotropic_index_matrix(radial_order) 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. mapmri_isotropic_odf_matrix(radial_order, …) Compute the isotropic MAPMRI ODF matrix [1] Eq. mapmri_isotropic_odf_sh_matrix(radial_order, …) Compute the isotropic MAPMRI ODF matrix [1] Eq. mapmri_isotropic_phi_matrix(radial_order, mu, q) Three dimensional isotropic MAPMRI signal basis function from [1] Eq. mapmri_isotropic_psi_matrix(radial_order, …) Three dimensional isotropic MAPMRI propagator basis function from [1] Eq. mapmri_isotropic_radial_pdf_basis(j, l, mu, r) Radial part of the isotropic 1D-SHORE propagator basis [1] eq. Radial part of the isotropic 1D-SHORE signal basis [1] eq. mapmri_laplacian_reg_matrix(ind_mat, mu, …) Put the Laplacian regularization matrix together [Rc66aaccd07c1-1] eq. mapmri_odf_matrix(radial_order, mu, s, vertices) Compute the MAPMRI ODF matrix [1] Eq. mapmri_phi_1d(n, q, mu) One dimensional MAPMRI basis function from [1] Eq. mapmri_phi_matrix(radial_order, mu, q_gradients) Compute the MAPMRI phi matrix for the signal [1] eq. mapmri_psi_1d(n, x, mu) One dimensional MAPMRI propagator basis function from [1] Eq. mapmri_psi_matrix(radial_order, mu, rgrad) Compute the MAPMRI psi matrix for the propagator [1] eq. mfactorial factorial(x) -> Integral multi_voxel_fit(single_voxel_fit) Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition optional_package(name[, trip_msg]) Return package-like thing and module setup for package name real_sph_harm(m, n, theta, phi) Compute real spherical harmonics. sfactorial(n[, exact]) The factorial of a number or array of numbers. sph_harm_ind_list(sh_order) Returns the degree (n) and order (m) of all the symmetric spherical harmonics of degree less then or equal to sh_order. warn Issue a warning, or maybe ignore it or raise an exception.

## Module: reconst.mcsd

 MSDeconvFit(model, coeff, mask) Attributes MultiShellDeconvModel(gtab, response[, …]) Methods MultiShellResponse(response, sh_order, shells) Attributes QpFitter(X, reg) Methods multi_tissue_basis(gtab, sh_order, iso_comp) Builds a basis for multi-shell multi-tissue CSD model. multi_voxel_fit(single_voxel_fit) Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition optional_package(name[, trip_msg]) Return package-like thing and module setup for package name solve_qp(P, Q, G, H) Helper function to set up and solve the Quadratic Program (QP) in CVXPY.

## Module: reconst.msdki

Classes and functions for fitting the mean signal diffusion kurtosis model

 MeanDiffusionKurtosisFit(model, model_params) Attributes MeanDiffusionKurtosisModel(gtab[, bmag, …]) Mean signal Diffusion Kurtosis Model ReconstModel(gtab) Abstract class for signal reconstruction models auto_attr(func) Decorator to create OneTimeProperty attributes. check_multi_b(gtab, n_bvals[, non_zero, bmag]) Check if you have enough different b-values in your gradient table design_matrix(ubvals) Constructs design matrix for the mean signal diffusion kurtosis model mean_signal_bvalue(data, gtab[, bmag]) Computes the average signal across different diffusion directions for each unique b-value msdki_prediction(msdki_params, gtab[, S0]) Predict the mean signal given the parameters of the mean signal DKI, an GradientTable object and S0 signal. ndindex(shape) An N-dimensional iterator object to index arrays. round_bvals(bvals[, bmag]) “This function rounds the b-values unique_bvals(bvals[, bmag, rbvals]) This function gives the unique rounded b-values of the data wls_fit_msdki(design_matrix, msignal, ng[, …]) Fits the mean signal diffusion kurtosis imaging based on a weighted least square solution [1].

## Module: reconst.multi_voxel

Tools to easily make multi voxel models

 CallableArray An array which can be called like a function MultiVoxelFit(model, fit_array, mask) Holds an array of fits and allows access to their attributes and methods ReconstFit(model, data) Abstract class which holds the fit result of ReconstModel as_strided(x[, shape, strides, subok, writeable]) Create a view into the array with the given shape and strides. multi_voxel_fit(single_voxel_fit) Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition ndindex(shape) An N-dimensional iterator object to index arrays.

## Module: reconst.odf

 OdfFit(model, data) Methods OdfModel(gtab) An abstract class to be sub-classed by specific odf models ReconstFit(model, data) Abstract class which holds the fit result of ReconstModel ReconstModel(gtab) Abstract class for signal reconstruction models gfa(samples) The general fractional anisotropy of a function evaluated on the unit sphere minmax_normalize(samples[, out]) Min-max normalization of a function evaluated on the unit sphere

## Module: reconst.qtdmri

 Cache Cache values based on a key object (such as a sphere or gradient table). QtdmriFit(model, qtdmri_coef, us, ut, …) Methods QtdmriModel(gtab[, radial_order, …]) The q:math:tau-dMRI model [1] to analytically and continuously represent the q:math:tau diffusion signal attenuation over diffusion sensitization q and diffusion time $$\tau$$. GCV_cost_function(weight, arguments) Generalized Cross Validation Function that is iterated [1]. H(value) Step function of H(x)=1 if x>=0 and zero otherwise. angular_basis_EAP_opt(j, l, m, r, theta, phi) angular_basis_opt(l, m, q, theta, phi) Angular basis independent of spatial scaling factor us. cart2sphere(x, y, z) Return angles for Cartesian 3D coordinates x, y, and z create_rt_space_grid(grid_size_r, …) Generates EAP grid (for potential positivity constraint). design_matrix_spatial(bvecs, qvals[, dtype]) Constructs design matrix for DTI weighted least squares or least squares fitting. elastic_crossvalidation(b0s_mask, E, M, L, lopt) cross-validation function to find the optimal weight of alpha for sparsity regularization when also Laplacian regularization is used. factorial(n[, exact]) The factorial of a number or array of numbers. factorial2(n[, exact]) Double factorial. fmin_l_bfgs_b(func, x0[, fprime, args, …]) Minimize a function func using the L-BFGS-B algorithm. generalized_crossvalidation(data, M, LR[, …]) Generalized Cross Validation Function [1]. genlaguerre(n, alpha[, monic]) Generalized (associated) Laguerre polynomial. A general function for creating diffusion MR gradients. l1_crossvalidation(b0s_mask, E, M[, …]) cross-validation function to find the optimal weight of alpha for sparsity regularization multi_voxel_fit(single_voxel_fit) Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition optional_package(name[, trip_msg]) Return package-like thing and module setup for package name part1_reg_matrix_tau(ind_mat, ut) Partial temporal Laplacian regularization matrix following Appendix B in [1]. part23_iso_reg_matrix_q(ind_mat, us) Partial spherical spatial Laplacian regularization matrix following the equation below Eq. part23_reg_matrix_q(ind_mat, U_mat, T_mat, us) Partial cartesian spatial Laplacian regularization matrix following second line of Eq. part23_reg_matrix_tau(ind_mat, ut) Partial temporal Laplacian regularization matrix following Appendix B in [1]. part4_iso_reg_matrix_q(ind_mat, us) Partial spherical spatial Laplacian regularization matrix following the equation below Eq. part4_reg_matrix_q(ind_mat, U_mat, us) Partial cartesian spatial Laplacian regularization matrix following equation Eq. part4_reg_matrix_tau(ind_mat, ut) Partial temporal Laplacian regularization matrix following Appendix B in [1]. qtdmri_anisotropic_scaling(data, q, bvecs, tau) Constructs design matrix for fitting an exponential to the diffusion time points. qtdmri_eap_matrix(radial_order, time_order, …) Constructs the design matrix as a product of 3 separated radial, angular and temporal design matrices. qtdmri_eap_matrix_(radial_order, time_order, …) qtdmri_index_matrix(radial_order, time_order) Computes the SHORE basis order indices according to [1]. qtdmri_isotropic_eap_matrix(radial_order, …) Constructs the design matrix as a product of 3 separated radial, angular and temporal design matrices. qtdmri_isotropic_eap_matrix_(radial_order, …) qtdmri_isotropic_index_matrix(radial_order, …) Computes the SHORE basis order indices according to [1]. Computes the spherical qt-dMRI Laplacian regularization matrix. qtdmri_isotropic_scaling(data, q, tau) Constructs design matrix for fitting an exponential to the diffusion time points. qtdmri_isotropic_signal_matrix(radial_order, …) Generates the matrix that maps the spherical qtdmri coefficients to MAP-MRI coefficients. qtdmri_laplacian_reg_matrix(ind_mat, us, ut) Computes the cartesian qt-dMRI Laplacian regularization matrix. Normalization factor for Spherical MAP-MRI basis. Normalization factor for Cartesian MAP-MRI basis. qtdmri_number_of_coefficients(radial_order, …) Computes the total number of coefficients of the qtdmri basis given a radial and temporal order. qtdmri_signal_matrix(radial_order, …) Constructs the design matrix as a product of 3 separated radial, angular and temporal design matrices. qtdmri_signal_matrix_(radial_order, …[, …]) Function to generate the qtdmri signal basis. Normalization factor for the temporal basis qtdmri_to_mapmri_matrix(radial_order, …) Generates the matrix that maps the qtdmri coefficients to MAP-MRI coefficients. radial_basis_EAP_opt(j, l, us, r) radial_basis_opt(j, l, us, q) Spatial basis dependent on spatial scaling factor us real_sph_harm(m, n, theta, phi) Compute real spherical harmonics. temporal_basis(o, ut, tau) 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). warn Issue a warning, or maybe ignore it or raise an exception.

## Module: reconst.sfm

The Sparse Fascicle Model.

This is an implementation of the sparse fascicle model described in [R204bb22f26e5-Rokem2015]. The multi b-value version of this model is described in [R204bb22f26e5-Rokem2014].

R204bb22f26e5-Rokem2015

Ariel Rokem, Jason D. Yeatman, Franco Pestilli, Kendrick N. Kay, Aviv Mezer, Stefan van der Walt, Brian A. Wandell (2015). Evaluating the accuracy of diffusion MRI models in white matter. PLoS ONE 10(4): e0123272. doi:10.1371/journal.pone.0123272

R204bb22f26e5-Rokem2014

Ariel Rokem, Kimberly L. Chan, Jason D. Yeatman, Franco Pestilli, Brian A. Wandell (2014). Evaluating the accuracy of diffusion models at multiple b-values with cross-validation. ISMRM 2014.

 Cache Cache values based on a key object (such as a sphere or gradient table). ExponentialIsotropicFit(model, params) A fit to the ExponentialIsotropicModel object, based on data. Representing the isotropic signal as a fit to an exponential decay function with b-values IsotropicFit(model, params) A fit object for representing the isotropic signal as the mean of the diffusion-weighted signal. IsotropicModel(gtab) A base-class for the representation of isotropic signals. ReconstFit(model, data) Abstract class which holds the fit result of ReconstModel ReconstModel(gtab) Abstract class for signal reconstruction models SparseFascicleFit(model, beta, S0, iso) Methods SparseFascicleModel(gtab[, sphere, …]) Methods auto_attr(func) Decorator to create OneTimeProperty attributes. nanmean(a[, axis, dtype, out, keepdims]) Compute the arithmetic mean along the specified axis, ignoring NaNs. optional_package(name[, trip_msg]) Return package-like thing and module setup for package name sfm_design_matrix(gtab, sphere, response[, mode]) Construct the SFM design matrix

## Module: reconst.shm

Tools for using spherical harmonic models to fit diffusion data

### References

Aganj, I., et al. 2009. ODF Reconstruction in Q-Ball Imaging With Solid

Angle Consideration.

Descoteaux, M., et al. 2007. Regularized, fast, and robust analytical

Q-ball imaging.

Tristan-Vega, A., et al. 2010. A new methodology for estimation of fiber

populations in white matter of the brain with Funk-Radon transform.

Tristan-Vega, A., et al. 2009. Estimation of fiber orientation probability

density functions in high angular resolution diffusion imaging.

Note about the Transpose: In the literature the matrix representation of these methods is often written as Y = Bx where B is some design matrix and Y and x are column vectors. In our case the input data, a dwi stored as a nifti file for example, is stored as row vectors (ndarrays) of the form (x, y, z, n), where n is the number of diffusion directions. We could transpose and reshape the data to be (n, x*y*z), so that we could directly plug it into the above equation. However, I have chosen to keep the data as is and implement the relevant equations rewritten in the following form: Y.T = x.T B.T, or in python syntax data = np.dot(sh_coef, B.T) where data is Y.T and sh_coef is x.T.

 Cache Cache values based on a key object (such as a sphere or gradient table). CsaOdfModel(gtab, sh_order[, smooth, …]) Implementation of Constant Solid Angle reconstruction method. LooseVersion([vstring]) Version numbering for anarchists and software realists. OdfFit(model, data) Methods OdfModel(gtab) An abstract class to be sub-classed by specific odf models OpdtModel(gtab, sh_order[, smooth, …]) Implementation of Orientation Probability Density Transform reconstruction method. QballBaseModel(gtab, sh_order[, smooth, …]) To be subclassed by Qball type models. QballModel(gtab, sh_order[, smooth, …]) Implementation of regularized Qball reconstruction method. ResidualBootstrapWrapper(signal_object, B, …) Returns a residual bootstrap sample of the signal_object when indexed SphHarmFit(model, shm_coef, mask) Diffusion data fit to a spherical harmonic model SphHarmModel(gtab) To be subclassed by all models that return a SphHarmFit when fit. anisotropic_power(sh_coeffs[, norm_factor, …]) Calculates anisotropic power map with a given SH coefficient matrix auto_attr(func) Decorator to create OneTimeProperty attributes. bootstrap_data_array(data, H, R[, permute]) Applies the Residual Bootstraps to the data given H and R bootstrap_data_voxel(data, H, R[, permute]) Like bootstrap_data_array but faster when for a single voxel calculate_max_order(n_coeffs) Calculate the maximal harmonic order, given that you know the cart2sphere(x, y, z) Return angles for Cartesian 3D coordinates x, y, and z concatenate([axis, out]) Join a sequence of arrays along an existing axis. diag(v[, k]) Extract a diagonal or construct a diagonal array. diff(a[, n, axis, prepend, append]) Calculate the n-th discrete difference along the given axis. dot(a, b[, out]) Dot product of two arrays. empty(shape[, dtype, order]) Return a new array of given shape and type, without initializing entries. eye(N[, M, k, dtype, order]) Return a 2-D array with ones on the diagonal and zeros elsewhere. forward_sdeconv_mat(r_rh, n) Build forward spherical deconvolution matrix gen_dirac(m, n, theta, phi) Generate Dirac delta function orientated in (theta, phi) on the sphere Returns the hat matrix for the design matrix B lazy_index(index) Produces a lazy index Returns a matrix for computing leveraged, centered residuals from data lpn(n, z) Legendre function of the first kind. normalize_data(data, where_b0[, min_signal, out]) Normalizes the data with respect to the mean b0 order_from_ncoef(ncoef) Given a number n of coefficients, calculate back the sh_order pinv(a[, rcond]) Compute the (Moore-Penrose) pseudo-inverse of a matrix. randint(low[, high, size, dtype]) Return random integers from low (inclusive) to high (exclusive). real_sph_harm(m, n, theta, phi) Compute real spherical harmonics. real_sym_sh_basis(sh_order, theta, phi) Samples a real symmetric spherical harmonic basis at point on the sphere real_sym_sh_mrtrix(sh_order, theta, phi) Compute real spherical harmonics as in Tournier 2007 [2], where the real harmonic $$Y^m_n$$ is defined to be. sf_to_sh(sf, sphere[, sh_order, basis_type, …]) Spherical function to spherical harmonics (SH). sh_to_rh(r_sh, m, n) Spherical harmonics (SH) to rotational harmonics (RH) sh_to_sf(sh, sphere, sh_order[, basis_type]) Spherical harmonics (SH) to spherical function (SF). sh_to_sf_matrix(sphere, sh_order[, …]) Matrix that transforms Spherical harmonics (SH) to spherical function (SF). smooth_pinv(B, L) Regularized pseudo-inverse sph_harm_ind_list(sh_order) Returns the degree (n) and order (m) of all the symmetric spherical harmonics of degree less then or equal to sh_order. spherical_harmonics(m, n, theta, phi) Compute spherical harmonics svd(a[, full_matrices, compute_uv]) Singular Value Decomposition. unique(ar[, return_index, return_inverse, …]) Find the unique elements of an array.

## Module: reconst.shore

 Cache Cache values based on a key object (such as a sphere or gradient table). ShoreFit(model, shore_coef) Attributes ShoreModel(gtab[, radial_order, zeta, …]) Simple Harmonic Oscillator based Reconstruction and Estimation (SHORE) [1] of the diffusion signal. cart2sphere(x, y, z) Return angles for Cartesian 3D coordinates x, y, and z create_rspace(gridsize, radius_max) Create the real space table, that contains the points in which Find x!. genlaguerre(n, alpha[, monic]) Generalized (associated) Laguerre polynomial. l_shore(radial_order) Returns the angular regularisation matrix for SHORE basis multi_voxel_fit(single_voxel_fit) Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition n_shore(radial_order) Returns the angular regularisation matrix for SHORE basis optional_package(name[, trip_msg]) Return package-like thing and module setup for package name real_sph_harm(m, n, theta, phi) Compute real spherical harmonics. shore_indices(radial_order, index) Given the basis order and the index, return the shore indices n, l, m for modified Merlet’s 3D-SHORE ..math:: :nowrap: begin{equation} textbf{E}(qtextbf{u})=sum_{l=0, even}^{N_{max}} sum_{n=l}^{(N_{max}+l)/2} sum_{m=-l}^l c_{nlm} phi_{nlm}(qtextbf{u}) end{equation} shore_matrix(radial_order, zeta, gtab[, tau]) Compute the SHORE matrix for modified Merlet’s 3D-SHORE [1] shore_matrix_odf(radial_order, zeta, …) Compute the SHORE ODF matrix [1]” shore_matrix_pdf(radial_order, zeta, rtab) Compute the SHORE propagator matrix [1]” shore_order(n, l, m) Given the indices (n,l,m) of the basis, return the minimum order for those indices and their index for modified Merlet’s 3D-SHORE. warn Issue a warning, or maybe ignore it or raise an exception.

## Module: reconst.utils

 Construct B design matrix for DKI.

### bench

dipy.reconst.bench(label='fast', verbose=1, extra_argv=None)

Run benchmarks for module using nose.

Parameters
label{‘fast’, ‘full’, ‘’, attribute identifier}, optional

Identifies the benchmarks to run. This can be a string to pass to the nosetests executable with the ‘-A’ option, or one of several special values. Special values are:

• ‘fast’ - the default - which corresponds to the nosetests -A option of ‘not slow’.

• ‘full’ - fast (as above) and slow benchmarks as in the ‘no -A’ option to nosetests - this is the same as ‘’.

• None or ‘’ - run all tests.

• attribute_identifier - string passed directly to nosetests as ‘-A’.

verboseint, optional

Verbosity value for benchmark outputs, in the range 1-10. Default is 1.

extra_argvlist, optional

List with any extra arguments to pass to nosetests.

Returns
successbool

Returns True if running the benchmarks works, False if an error occurred.

Notes

Benchmarks are like tests, but have names starting with “bench” instead of “test”, and can be found under the “benchmarks” sub-directory of the module.

Each NumPy module exposes bench in its namespace to run all benchmarks for it.

Examples

>>> success = np.lib.bench() #doctest: +SKIP
Running benchmarks for numpy.lib
...
using 562341 items:
unique:
0.11
unique1d:
0.11
ratio: 1.0
nUnique: 56230 == 56230
...
OK

>>> success #doctest: +SKIP
True


### test

dipy.reconst.test(label='fast', verbose=1, extra_argv=None, doctests=False, coverage=False, raise_warnings=None, timer=False)

Run tests for module using nose.

Parameters
label{‘fast’, ‘full’, ‘’, attribute identifier}, optional

Identifies the tests to run. This can be a string to pass to the nosetests executable with the ‘-A’ option, or one of several special values. Special values are:

• ‘fast’ - the default - which corresponds to the nosetests -A option of ‘not slow’.

• ‘full’ - fast (as above) and slow tests as in the ‘no -A’ option to nosetests - this is the same as ‘’.

• None or ‘’ - run all tests.

• attribute_identifier - string passed directly to nosetests as ‘-A’.

verboseint, optional

Verbosity value for test outputs, in the range 1-10. Default is 1.

extra_argvlist, optional

List with any extra arguments to pass to nosetests.

doctestsbool, optional

If True, run doctests in module. Default is False.

coveragebool, optional

If True, report coverage of NumPy code. Default is False. (This requires the coverage module).

raise_warningsNone, str or sequence of warnings, optional

This specifies which warnings to configure as ‘raise’ instead of being shown once during the test execution. Valid strings are:

• “develop” : equals (Warning,)

• “release” : equals (), do not raise on any warnings.

timerbool or int, optional

Timing of individual tests with nose-timer (which needs to be installed). If True, time tests and report on all of them. If an integer (say N), report timing results for N slowest tests.

Returns
resultobject

Returns the result of running the tests as a nose.result.TextTestResult object.

Notes

Each NumPy module exposes test in its namespace to run all tests for it. For example, to run all tests for numpy.lib:

>>> np.lib.test() #doctest: +SKIP


Examples

>>> result = np.lib.test() #doctest: +SKIP
Running unit tests for numpy.lib
...
Ran 976 tests in 3.933s


OK

>>> result.errors #doctest: +SKIP
[]
>>> result.knownfail #doctest: +SKIP
[]


### ReconstFit

class dipy.reconst.base.ReconstFit(model, data)

Bases: object

Abstract class which holds the fit result of ReconstModel

For example that could be holding FA or GFA etc.

__init__(model, data)

Initialize self. See help(type(self)) for accurate signature.

### ReconstModel

class dipy.reconst.base.ReconstModel(gtab)

Bases: object

Abstract class for signal reconstruction models

Methods

 fit
__init__(gtab)

Initialization of the abstract class for signal reconstruction models

Parameters
fit(data, mask=None, **kwargs)

### bench_bounding_box

dipy.reconst.benchmarks.bench_bounding_box.bench_bounding_box()

### bounding_box

dipy.reconst.benchmarks.bench_bounding_box.bounding_box(vol)

Compute the bounding box of nonzero intensity voxels in the volume.

Parameters
volndarray

Volume to compute bounding box on.

Returns
npminslist

Array containg minimum index of each dimension

npmaxslist

Array containg maximum index of each dimension

### measure

dipy.reconst.benchmarks.bench_bounding_box.measure(code_str, times=1, label=None)

Return elapsed time for executing code in the namespace of the caller.

The supplied code string is compiled with the Python builtin compile. The precision of the timing is 10 milli-seconds. If the code will execute fast on this timescale, it can be executed many times to get reasonable timing accuracy.

Parameters
code_strstr

The code to be timed.

timesint, optional

The number of times the code is executed. Default is 1. The code is only compiled once.

labelstr, optional

A label to identify code_str with. This is passed into compile as the second argument (for run-time error messages).

Returns
elapsedfloat

Total elapsed time in seconds for executing code_str times times.

Examples

>>> times = 10
>>> etime = np.testing.measure('for i in range(1000): np.sqrt(i**2)', times=times)
>>> print("Time for a single execution : ", etime / times, "s")  # doctest: +SKIP
Time for a single execution :  0.005 s


### ConstrainedSphericalDeconvModel

class dipy.reconst.benchmarks.bench_csd.ConstrainedSphericalDeconvModel(gtab, response, reg_sphere=None, sh_order=8, lambda_=1, tau=0.1, convergence=50)

Methods

 cache_clear() Clear the cache. cache_get(tag, key[, default]) Retrieve a value from the cache. cache_set(tag, key, value) Store a value in the cache. fit(data[, mask]) Fit method for every voxel in data predict(sh_coeff[, gtab, S0]) Compute a signal prediction given spherical harmonic coefficients for the provided GradientTable class instance. sampling_matrix(sphere) The matrix needed to sample ODFs from coefficients of the model.
__init__(gtab, response, reg_sphere=None, sh_order=8, lambda_=1, tau=0.1, convergence=50)

Constrained Spherical Deconvolution (CSD) [1].

Spherical deconvolution computes a fiber orientation distribution (FOD), also called fiber ODF (fODF) [2], as opposed to a diffusion ODF as the QballModel or the CsaOdfModel. This results in a sharper angular profile with better angular resolution that is the best object to be used for later deterministic and probabilistic tractography [3].

A sharp fODF is obtained because a single fiber response function is injected as a priori knowledge. The response function is often data-driven and is thus provided as input to the ConstrainedSphericalDeconvModel. It will be used as deconvolution kernel, as described in [1].

Parameters
responsetuple or AxSymShResponse object

A tuple with two elements. The first is the eigen-values as an (3,) ndarray and the second is the signal value for the response function without diffusion weighting (i.e. S0). This is to be able to generate a single fiber synthetic signal. The response function will be used as deconvolution kernel ([1]).

reg_sphereSphere (optional)

sphere used to build the regularization B matrix. Default: ‘symmetric362’.

sh_orderint (optional)

maximal spherical harmonics order. Default: 8

lambda_float (optional)

weight given to the constrained-positivity regularization part of the deconvolution equation (see [1]). Default: 1

taufloat (optional)

threshold controlling the amplitude below which the corresponding fODF is assumed to be zero. Ideally, tau should be set to zero. However, to improve the stability of the algorithm, tau is set to tau*100 % of the mean fODF amplitude (here, 10% by default) (see [1]). Default: 0.1

convergenceint

Maximum number of iterations to allow the deconvolution to converge.

References

1(1,2,3,4,5,6)

Tournier, J.D., et al. NeuroImage 2007. Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution

2(1,2)

Descoteaux, M., et al. IEEE TMI 2009. Deterministic and Probabilistic Tractography Based on Complex Fibre Orientation Distributions

3(1,2)

Côté, M-A., et al. Medical Image Analysis 2013. Tractometer: Towards validation of tractography pipelines

4

Tournier, J.D, et al. Imaging Systems and Technology 2012. MRtrix: Diffusion Tractography in Crossing Fiber Regions

fit(data, mask=None)

Fit method for every voxel in data

predict(sh_coeff, gtab=None, S0=1.0)

Compute a signal prediction given spherical harmonic coefficients for the provided GradientTable class instance.

Parameters
sh_coeffndarray

The spherical harmonic representation of the FOD from which to make the signal prediction.

The gradients for which the signal will be predicted. Use the model’s gradient table by default.

S0ndarray or float

The non diffusion-weighted signal value.

Returns
pred_signdarray

The predicted signal.

### GradientTable

class dipy.reconst.benchmarks.bench_csd.GradientTable(gradients, big_delta=None, small_delta=None, b0_threshold=50)

Bases: object

Parameters

Diffusion gradients. The direction of each of these vectors corresponds to the b-vector, and the length corresponds to the b-value.

b0_thresholdfloat

Gradients with b-value less than or equal to b0_threshold are considered as b0s i.e. without diffusion weighting.

gradient_table

Notes

The GradientTable object is immutable. Do NOT assign attributes. If you have your gradient table in a bval & bvec format, we recommend using the factory function gradient_table

Attributes

bvals(N,) ndarray

The b-value, or magnitude, of each gradient direction.

qvals: (N,) ndarray

The q-value for each gradient direction. Needs big and small delta.

bvecs(N,3) ndarray

The direction, represented as a unit vector, of each gradient.

Boolean array indicating which gradients have no diffusion weighting, ie b-value is close to 0.

b0_thresholdfloat

Gradients with b-value less than or equal to b0_threshold are considered to not have diffusion weighting.

Methods

__init__(gradients, big_delta=None, small_delta=None, b0_threshold=50)

Constructor for GradientTable class

b0s_mask()
bvals()
bvecs()
gradient_strength()
property info
qvals()
tau()

### bench_csdeconv

dipy.reconst.benchmarks.bench_csd.bench_csdeconv(center=(50, 40, 40), width=12)

dipy.reconst.benchmarks.bench_csd.num_grad(gtab)

dipy.reconst.benchmarks.bench_csd.read_stanford_labels()

Read stanford hardi data and label map

### bench_local_maxima

dipy.reconst.benchmarks.bench_peaks.bench_local_maxima()

### local_maxima

dipy.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.

Parameters
odfarray, 1d, dtype=double

The function evaluated on a set of discrete points.

edgesarray (N, 2)

The set of neighbor relations between the points. Every edge, ie edges[i, :], is a pair of neighboring points.

Returns
peak_valuesndarray

Value of odf at a maximum point. Peak values is sorted in descending order.

peak_indicesndarray

Indices of maximum points. Sorted in the same order as peak_values so odf[peak_indices[i]] == peak_values[i].

Notes

A point is a local maximum if it is > at least one neighbor and >= all neighbors. If no points meet the above criteria, 1 maximum is returned such that odf[maximum] == max(odf).

### measure

dipy.reconst.benchmarks.bench_peaks.measure(code_str, times=1, label=None)

Return elapsed time for executing code in the namespace of the caller.

The supplied code string is compiled with the Python builtin compile. The precision of the timing is 10 milli-seconds. If the code will execute fast on this timescale, it can be executed many times to get reasonable timing accuracy.

Parameters
code_strstr

The code to be timed.

timesint, optional

The number of times the code is executed. Default is 1. The code is only compiled once.

labelstr, optional

A label to identify code_str with. This is passed into compile as the second argument (for run-time error messages).

Returns
elapsedfloat

Total elapsed time in seconds for executing code_str times times.

Examples

>>> times = 10
>>> etime = np.testing.measure('for i in range(1000): np.sqrt(i**2)', times=times)
>>> print("Time for a single execution : ", etime / times, "s")  # doctest: +SKIP
Time for a single execution :  0.005 s


### unique_edges

dipy.reconst.benchmarks.bench_peaks.unique_edges(faces, return_mapping=False)

Extract all unique edges from given triangular faces.

Parameters
faces(N, 3) ndarray

Vertex indices forming triangular faces.

return_mappingbool

If true, a mapping to the edges of each face is returned.

Returns
edges(N, 2) ndarray

Unique edges.

mapping(N, 3)

For each face, [x, y, z], a mapping to it’s edges [a, b, c].

   y
/               /               a/    
/                  /                   /__________          x      c     z


### bench_quick_squash

dipy.reconst.benchmarks.bench_squash.bench_quick_squash()

### measure

dipy.reconst.benchmarks.bench_squash.measure(code_str, times=1, label=None)

Return elapsed time for executing code in the namespace of the caller.

The supplied code string is compiled with the Python builtin compile. The precision of the timing is 10 milli-seconds. If the code will execute fast on this timescale, it can be executed many times to get reasonable timing accuracy.

Parameters
code_strstr

The code to be timed.

timesint, optional

The number of times the code is executed. Default is 1. The code is only compiled once.

labelstr, optional

A label to identify code_str with. This is passed into compile as the second argument (for run-time error messages).

Returns
elapsedfloat

Total elapsed time in seconds for executing code_str times times.

Examples

>>> times = 10
>>> etime = np.testing.measure('for i in range(1000): np.sqrt(i**2)', times=times)
>>> print("Time for a single execution : ", etime / times, "s")  # doctest: +SKIP
Time for a single execution :  0.005 s


### ndindex

dipy.reconst.benchmarks.bench_squash.ndindex(shape)

An N-dimensional iterator object to index arrays.

Given the shape of an array, an ndindex instance iterates over the N-dimensional index of the array. At each iteration a tuple of indices is returned; the last dimension is iterated over first.

Parameters
shapetuple of ints

The dimensions of the array.

Examples

>>> from dipy.core.ndindex import ndindex
>>> shape = (3, 2, 1)
>>> for index in ndindex(shape):
...     print(index)
(0, 0, 0)
(0, 1, 0)
(1, 0, 0)
(1, 1, 0)
(2, 0, 0)
(2, 1, 0)


### old_squash

dipy.reconst.benchmarks.bench_squash.old_squash(arr, mask=None, fill=0)

Try and make a standard array from an object array

This function takes an object array and attempts to convert it to a more useful dtype. If array can be converted to a better dtype, Nones are replaced by fill. To make the behaviour of this function more clear, here are the most common cases:

1. arr is an array of scalars of type T. Returns an array like arr.astype(T)

2. arr is an array of arrays. All items in arr have the same shape S. Returns an array with shape arr.shape + S.

3. arr is an array of arrays of different shapes. Returns arr.

4. Items in arr are not ndarrys or scalars. Returns arr.

Parameters
arrarray, dtype=object

The array to be converted.

Where arr has Nones.

fillnumber, optional

Nones are replaced by fill.

Returns
resultarray

Examples

>>> arr = np.empty(3, dtype=object)
>>> arr.fill(2)
>>> old_squash(arr)
array([2, 2, 2])
>>> arr[0] = None
>>> old_squash(arr)
array([0, 2, 2])
>>> arr.fill(np.ones(2))
>>> r = old_squash(arr)
>>> r.shape == (3, 2)
True
>>> r.dtype
dtype('float64')


### quick_squash

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:

1. obj_arr is an array of scalars of type T. Returns an array like obj_arr.astype(T)

2. 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

3. obj_arr is an array of arrays of different shapes. Returns obj_arr.

4. Items in obj_arr are not ndarrays or scalars. Returns obj_arr.

Parameters
obj_arrarray, dtype=object

The array to be converted.

mask is nonzero where obj_arr has Nones.

fillnumber, optional

Nones are replaced by fill.

Returns
resultarray

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')


### reduce

dipy.reconst.benchmarks.bench_squash.reduce(function, sequence[, initial]) → value

Apply a function of two arguments cumulatively to the items of a sequence, from left to right, so as to reduce the sequence to a single value. For example, reduce(lambda x, y: x+y, [1, 2, 3, 4, 5]) calculates ((((1+2)+3)+4)+5). If initial is present, it is placed before the items of the sequence in the calculation, and serves as a default when the sequence is empty.

### bench_vec_val_vect

dipy.reconst.benchmarks.bench_vec_val_sum.bench_vec_val_vect()

### measure

dipy.reconst.benchmarks.bench_vec_val_sum.measure(code_str, times=1, label=None)

Return elapsed time for executing code in the namespace of the caller.

The supplied code string is compiled with the Python builtin compile. The precision of the timing is 10 milli-seconds. If the code will execute fast on this timescale, it can be executed many times to get reasonable timing accuracy.

Parameters
code_strstr

The code to be timed.

timesint, optional

The number of times the code is executed. Default is 1. The code is only compiled once.

labelstr, optional

A label to identify code_str with. This is passed into compile as the second argument (for run-time error messages).

Returns
elapsedfloat

Total elapsed time in seconds for executing code_str times times.

Examples

>>> times = 10
>>> etime = np.testing.measure('for i in range(1000): np.sqrt(i**2)', times=times)
>>> print("Time for a single execution : ", etime / times, "s")  # doctest: +SKIP
Time for a single execution :  0.005 s


### randn

dipy.reconst.benchmarks.bench_vec_val_sum.randn(d0, d1, ..., dn)

Return a sample (or samples) from the “standard normal” distribution.

If positive, int_like or int-convertible 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 (if any of the $$d_i$$ are floats, they are first converted to integers by truncation). A single float randomly sampled from the distribution is returned if no argument is provided.

This is a convenience function. If you want an interface that takes a tuple as the first argument, use numpy.random.standard_normal instead.

Parameters
d0, d1, …, dnint, optional

The dimensions of the returned array, should be all positive. If no argument is given a single Python float is returned.

Returns
Zndarray or float

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.

standard_normal

Similar, but takes a tuple as its argument.

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):

>>> 2.5 * np.random.randn(2, 4) + 3
array([[-4.49401501,  4.00950034, -1.81814867,  7.29718677],  #random
[ 0.39924804,  4.68456316,  4.99394529,  4.84057254]]) #random


### vec_val_vect

dipy.reconst.benchmarks.bench_vec_val_sum.vec_val_vect()

Vectorize vecs.diag(vals).vecs.T for last 2 dimensions of vecs

Parameters
vecsshape (…, M, N) array

containing tensor in last two dimensions; M, N usually equal to (3, 3)

valsshape (…, N) array

diagonal values carried in last dimension, ... shape above must match that for vecs

Returns
resshape (…, M, M) array

For all the dimensions ellided by ..., loops to get (M, N) vec matrix, and (N,) vals vector, and calculates vec.dot(np.diag(val).dot(vec.T).

Raises
ValueErrornon-matching ... dimensions of vecs, vals
ValueErrornon-matching N dimensions of vecs, vals

Examples

Make a 3D array where the first dimension is only 1

>>> vecs = np.arange(9).reshape((1, 3, 3))
>>> vals = np.arange(3).reshape((1, 3))
>>> vec_val_vect(vecs, vals)
array([[[   9.,   24.,   39.],
[  24.,   66.,  108.],
[  39.,  108.,  177.]]])


That’s the same as the 2D case (apart from the float casting):

>>> vecs = np.arange(9).reshape((3, 3))
>>> vals = np.arange(3)
>>> np.dot(vecs, np.dot(np.diag(vals), vecs.T))
array([[  9,  24,  39],
[ 24,  66, 108],
[ 39, 108, 177]])


### with_einsum

dipy.reconst.benchmarks.bench_vec_val_sum.with_einsum(f)

### Cache

class dipy.reconst.cache.Cache

Bases: object

Cache values based on a key object (such as a sphere or gradient table).

Notes

This class is meant to be used as a mix-in:

class MyModel(Model, Cache):
pass

class MyModelFit(Fit):
pass


Inside a method on the fit, typical usage would be:

def odf(sphere):
M = self.model.cache_get('odf_basis_matrix', key=sphere)

if M is None:
M = self._compute_basis_matrix(sphere)
self.model.cache_set('odf_basis_matrix', key=sphere, value=M)


Methods

 Clear the cache. cache_get(tag, key[, default]) Retrieve a value from the cache. cache_set(tag, key, value) Store a value in the cache.
__init__(self, /, *args, **kwargs) Initialize self. See help(type(self)) for accurate signature. cache_clear() Clear the cache. cache_get(tag, key, default=None) Retrieve a value from the cache. Parameters tagstr Description of the cached value. keyobject Key object used to look up the cached value. defaultobject Value to be returned if no cached entry is found. Returns vobject Value from the cache associated with (tag, key). Returns default if no cached entry is found. cache_set(tag, key, value) Store a value in the cache. Parameters tagstr Description of the cached value. keyobject Key object used to look up the cached value. valueobject Value stored in the cache for each unique combination of (tag, key). Examples >>> def compute_expensive_matrix(parameters): ... # Imagine the following computation is very expensive ... return (p**2 for p in parameters)  >>> 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  ### auto_attr dipy.reconst.cache.auto_attr(func) Decorator to create OneTimeProperty attributes. Parameters funcmethod The method that will be called the first time to compute a value. Afterwards, the method’s name will be a standard attribute holding the value of this computation. Examples >>> class MagicProp(object): ... @auto_attr ... def a(self): ... return 99 ... >>> x = MagicProp() >>> 'a' in x.__dict__ False >>> x.a 99 >>> 'a' in x.__dict__ True  ### coeff_of_determination dipy.reconst.cross_validation.coeff_of_determination(data, model, axis=-1) Calculate the coefficient of determination for a model prediction, relative to data. Parameters datandarray The data modelndarray The predictions of a model for this data. Same shape as the data. axis: int, optional The axis along which different samples are laid out (default: -1). Returns CODndarray The coefficient of determination. This has shape data.shape[:-1] rac{SSE}{SSD}) where SSE is the sum of the squared error between the model and the data (sum of the squared residuals) and SSD is the sum of the squares of the deviations of the data from the mean of the data (variance * N). ### kfold_xval dipy.reconst.cross_validation.kfold_xval(model, data, folds, *model_args, **model_kwargs) Perform k-fold cross-validation to generate out-of-sample predictions for each measurement. Parameters modelModel class instance The type of the model to use for prediction. The corresponding Fit object must have a predict function implementd One of the following: reconst.dti.TensorModel or reconst.csdeconv.ConstrainedSphericalDeconvModel. datandarray Diffusion MRI data acquired with the GradientTable of the model. Shape will typically be (x, y, z, b) where xyz are spatial dimensions and b is the number of bvals/bvecs in the GradientTable. foldsint The number of divisions to apply to the data model_argslist Additional arguments to the model initialization model_kwargsdict Additional key-word arguments to the model initialization. If contains the kwarg mask, this will be used as a key-word argument to the fit method of the model object, rather than being used in the initialization of the model object Notes This function assumes that a prediction API is implemented in the Model class for which prediction is conducted. That is, the Fit object that gets generated upon fitting the model needs to have a predict method, which receives a GradientTable class instance as input and produces a predicted signal as output. It also assumes that the model object has bval and bvec attributes holding b-values and corresponding unit vectors. References 1 Rokem, A., Chan, K.L. Yeatman, J.D., Pestilli, F., Mezer, A., Wandell, B.A., 2014. Evaluating the accuracy of diffusion models at multiple b-values with cross-validation. ISMRM 2014. ### AxSymShResponse class dipy.reconst.csdeconv.AxSymShResponse(S0, dwi_response, bvalue=None) Bases: object A simple wrapper for response functions represented using only axially symmetric, even spherical harmonic functions (ie, m == 0 and n even). Parameters S0float Signal with no diffusion weighting. dwi_responsearray Response function signal as coefficients to axially symmetric, even spherical harmonic. Methods  basis(sphere) A basis that maps the response coefficients onto a sphere. on_sphere(sphere) Evaluates the response function on sphere. __init__(S0, dwi_response, bvalue=None) Initialize self. See help(type(self)) for accurate signature. basis(sphere) A basis that maps the response coefficients onto a sphere. on_sphere(sphere) Evaluates the response function on sphere. ### ConstrainedSDTModel class dipy.reconst.csdeconv.ConstrainedSDTModel(gtab, ratio, reg_sphere=None, sh_order=8, lambda_=1.0, tau=0.1) Methods  cache_clear() Clear the cache. cache_get(tag, key[, default]) Retrieve a value from the cache. cache_set(tag, key, value) Store a value in the cache. fit(data[, mask]) Fit method for every voxel in data sampling_matrix(sphere) The matrix needed to sample ODFs from coefficients of the model. __init__(gtab, ratio, reg_sphere=None, sh_order=8, lambda_=1.0, tau=0.1) Spherical Deconvolution Transform (SDT) [1]. The SDT computes a fiber orientation distribution (FOD) as opposed to a diffusion ODF as the QballModel or the CsaOdfModel. This results in a sharper angular profile with better angular resolution. The Constrained SDTModel is similar to the Constrained CSDModel but mathematically it deconvolves the q-ball ODF as oppposed to the HARDI signal (see [1] for a comparison and a through discussion). A sharp fODF is obtained because a single fiber response function is injected as a priori knowledge. In the SDTModel, this response is a single fiber q-ball ODF as opposed to a single fiber signal function for the CSDModel. The response function will be used as deconvolution kernel. Parameters gtabGradientTable ratiofloat ratio of the smallest vs the largest eigenvalue of the single prolate tensor response function reg_sphereSphere sphere used to build the regularization B matrix sh_orderint maximal spherical harmonics order lambda_float weight given to the constrained-positivity regularization part of the deconvolution equation taufloat threshold (tau *mean(fODF)) controlling the amplitude below which the corresponding fODF is assumed to be zero. References 1(1,2,3) Descoteaux, M., et al. IEEE TMI 2009. Deterministic and Probabilistic Tractography Based on Complex Fibre Orientation Distributions. fit(data, mask=None) Fit method for every voxel in data ### ConstrainedSphericalDeconvModel class dipy.reconst.csdeconv.ConstrainedSphericalDeconvModel(gtab, response, reg_sphere=None, sh_order=8, lambda_=1, tau=0.1, convergence=50) Methods  cache_clear() Clear the cache. cache_get(tag, key[, default]) Retrieve a value from the cache. cache_set(tag, key, value) Store a value in the cache. fit(data[, mask]) Fit method for every voxel in data predict(sh_coeff[, gtab, S0]) Compute a signal prediction given spherical harmonic coefficients for the provided GradientTable class instance. sampling_matrix(sphere) The matrix needed to sample ODFs from coefficients of the model. __init__(gtab, response, reg_sphere=None, sh_order=8, lambda_=1, tau=0.1, convergence=50) Constrained Spherical Deconvolution (CSD) [1]. Spherical deconvolution computes a fiber orientation distribution (FOD), also called fiber ODF (fODF) [2], as opposed to a diffusion ODF as the QballModel or the CsaOdfModel. This results in a sharper angular profile with better angular resolution that is the best object to be used for later deterministic and probabilistic tractography [3]. A sharp fODF is obtained because a single fiber response function is injected as a priori knowledge. The response function is often data-driven and is thus provided as input to the ConstrainedSphericalDeconvModel. It will be used as deconvolution kernel, as described in [1]. Parameters gtabGradientTable responsetuple or AxSymShResponse object A tuple with two elements. The first is the eigen-values as an (3,) ndarray and the second is the signal value for the response function without diffusion weighting (i.e. S0). This is to be able to generate a single fiber synthetic signal. The response function will be used as deconvolution kernel ([1]). reg_sphereSphere (optional) sphere used to build the regularization B matrix. Default: ‘symmetric362’. sh_orderint (optional) maximal spherical harmonics order. Default: 8 lambda_float (optional) weight given to the constrained-positivity regularization part of the deconvolution equation (see [1]). Default: 1 taufloat (optional) threshold controlling the amplitude below which the corresponding fODF is assumed to be zero. Ideally, tau should be set to zero. However, to improve the stability of the algorithm, tau is set to tau*100 % of the mean fODF amplitude (here, 10% by default) (see [1]). Default: 0.1 convergenceint Maximum number of iterations to allow the deconvolution to converge. References 1(1,2,3,4,5,6) Tournier, J.D., et al. NeuroImage 2007. Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution 2(1,2) Descoteaux, M., et al. IEEE TMI 2009. Deterministic and Probabilistic Tractography Based on Complex Fibre Orientation Distributions 3(1,2) Côté, M-A., et al. Medical Image Analysis 2013. Tractometer: Towards validation of tractography pipelines 4 Tournier, J.D, et al. Imaging Systems and Technology 2012. MRtrix: Diffusion Tractography in Crossing Fiber Regions fit(data, mask=None) Fit method for every voxel in data predict(sh_coeff, gtab=None, S0=1.0) Compute a signal prediction given spherical harmonic coefficients for the provided GradientTable class instance. Parameters sh_coeffndarray The spherical harmonic representation of the FOD from which to make the signal prediction. gtabGradientTable The gradients for which the signal will be predicted. Use the model’s gradient table by default. S0ndarray or float The non diffusion-weighted signal value. Returns pred_signdarray The predicted signal. ### SphHarmFit class dipy.reconst.csdeconv.SphHarmFit(model, shm_coef, mask) Diffusion data fit to a spherical harmonic model Attributes shape shm_coeff The spherical harmonic coefficients of the odf Methods  odf(sphere) Samples the odf function on the points of a sphere predict([gtab, S0]) Predict the diffusion signal from the model coefficients.  gfa __init__(model, shm_coef, mask) Initialize self. See help(type(self)) for accurate signature. gfa() odf(sphere) Samples the odf function on the points of a sphere Parameters sphereSphere The points on which to sample the odf. Returns valuesndarray The value of the odf on each point of sphere. predict(gtab=None, S0=1.0) Predict the diffusion signal from the model coefficients. Parameters gtaba GradientTable class instance The directions and bvalues on which prediction is desired S0float array The mean non-diffusion-weighted signal in each voxel. Default: 1.0 in all voxels property shape property shm_coeff The spherical harmonic coefficients of the odf Make this a property for now, if there is a usecase for modifying the coefficients we can add a setter or expose the coefficients more directly ### SphHarmModel class dipy.reconst.csdeconv.SphHarmModel(gtab) To be subclassed by all models that return a SphHarmFit when fit. Methods  cache_clear() Clear the cache. cache_get(tag, key[, default]) Retrieve a value from the cache. cache_set(tag, key, value) Store a value in the cache. fit(data) To be implemented by specific odf models sampling_matrix(sphere) The matrix needed to sample ODFs from coefficients of the model. __init__(gtab) Initialization of the abstract class for signal reconstruction models Parameters gtabGradientTable class instance sampling_matrix(sphere) The matrix needed to sample ODFs from coefficients of the model. Parameters sphereSphere Points used to sample ODF. Returns sampling_matrixarray The size of the matrix will be (N, M) where N is the number of vertices on sphere and M is the number of coefficients needed by the model. ### TensorModel class dipy.reconst.csdeconv.TensorModel(gtab, fit_method='WLS', return_S0_hat=False, *args, **kwargs) Diffusion Tensor Methods  fit(data[, mask]) Fit method of the DTI model class predict(dti_params[, S0]) Predict a signal for this TensorModel class instance given parameters. __init__(gtab, fit_method='WLS', return_S0_hat=False, *args, **kwargs) A Diffusion Tensor Model [1], [2]. Parameters gtabGradientTable class instance fit_methodstr or callable str can be one of the following: ‘WLS’ for weighted least squares dti.wls_fit_tensor() ‘LS’ or ‘OLS’ for ordinary least squares dti.ols_fit_tensor() ‘NLLS’ for non-linear least-squares dti.nlls_fit_tensor() ‘RT’ or ‘restore’ or ‘RESTORE’ for RESTORE robust tensor fitting [3] dti.restore_fit_tensor() callable has to have the signature: fit_method(design_matrix, data, *args, **kwargs) return_S0_hatbool Boolean to return (True) or not (False) the S0 values for the fit. args, kwargsarguments and key-word arguments passed to the fit_method. See dti.wls_fit_tensor, dti.ols_fit_tensor for details min_signalfloat The minimum signal value. Needs to be a strictly positive number. Default: minimal signal in the data provided to fit. Notes In order to increase speed of processing, tensor fitting is done simultaneously over many voxels. Many fit_methods use the ‘step’ parameter to set the number of voxels that will be fit at once in each iteration. This is the chunk size as a number of voxels. A larger step value should speed things up, but it will also take up more memory. It is advisable to keep an eye on memory consumption as this value is increased. E.g., in iter_fit_tensor() we have a default step value of 1e4 References 1(1,2) Basser, P.J., Mattiello, J., LeBihan, D., 1994. Estimation of the effective self-diffusion tensor from the NMR spin echo. J Magn Reson B 103, 247-254. 2(1,2) Basser, P., Pierpaoli, C., 1996. Microstructural and physiological features of tissues elucidated by quantitative diffusion-tensor MRI. Journal of Magnetic Resonance 111, 209-219. 3(1,2) Lin-Ching C., Jones D.K., Pierpaoli, C. 2005. RESTORE: Robust estimation of tensors by outlier rejection. MRM 53: 1088-1095 fit(data, mask=None) Fit method of the DTI model class Parameters dataarray The measured signal from one voxel. maskarray A boolean array used to mark the coordinates in the data that should be analyzed that has the shape data.shape[:-1] predict(dti_params, S0=1.0) Predict a signal for this TensorModel class instance given parameters. Parameters dti_paramsndarray The last dimension should have 12 tensor parameters: 3 eigenvalues, followed by the 3 eigenvectors S0float or ndarray The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1 ### auto_response dipy.reconst.csdeconv.auto_response(gtab, data, roi_center=None, roi_radius=10, fa_thr=0.7, fa_callable=<function fa_superior>, return_number_of_voxels=False) Automatic estimation of response function using FA. Parameters gtabGradientTable datandarray diffusion data roi_centertuple, (3,) Center of ROI in data. If center is None, it is assumed that it is the center of the volume with shape data.shape[:3]. roi_radiusint radius of cubic ROI fa_thrfloat FA threshold fa_callablecallable A callable that defines an operation that compares FA with the fa_thr. The operator should have two positional arguments (e.g., fa_operator(FA, fa_thr)) and it should return a bool array. return_number_of_voxelsbool If True, returns the number of voxels used for estimating the response function. Returns responsetuple, (2,) (evals, S0) ratiofloat The ratio between smallest versus largest eigenvalue of the response. number of voxelsint (optional) The number of voxels used for estimating the response function. 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. For example we can use an ROI (20x20x20) at the center of the volume and store the signal values for the voxels with FA values higher than 0.7. Of course, if we haven’t precalculated FA we need to fit a Tensor model to the datasets. Which is what we do in this 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. If requested, the number of voxels used for estimating the response function is also returned, which can be used to judge the fidelity of the response function. As a rule of thumb, at least 300 voxels should be used to estimate a good response function (see [1]). References 1(1,2) Tournier, J.D., et al. NeuroImage 2004. Direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution ### cart2sphere dipy.reconst.csdeconv.cart2sphere(x, y, z) Return angles for Cartesian 3D coordinates x, y, and z See doc for sphere2cart for angle conventions and derivation of the formulae. $$0\le\theta\mathrm{(theta)}\le\pi$$ and $$-\pi\le\phi\mathrm{(phi)}\le\pi$$ Parameters xarray_like x coordinate in Cartesian space yarray_like y coordinate in Cartesian space zarray_like z coordinate Returns rarray radius thetaarray inclination (polar) angle phiarray azimuth angle ### csdeconv dipy.reconst.csdeconv.csdeconv(dwsignal, X, B_reg, tau=0.1, convergence=50, P=None) Constrained-regularized spherical deconvolution (CSD) [1] Deconvolves the axially symmetric single fiber response function r_rh in rotational harmonics coefficients from the diffusion weighted signal in dwsignal. Parameters dwsignalarray Diffusion weighted signals to be deconvolved. Xarray Prediction matrix which estimates diffusion weighted signals from FOD coefficients. B_regarray (N, B) SH basis matrix which maps FOD coefficients to FOD values on the surface of the sphere. B_reg should be scaled to account for lambda. taufloat Threshold controlling the amplitude below which the corresponding fODF is assumed to be zero. Ideally, tau should be set to zero. However, to improve the stability of the algorithm, tau is set to tau*100 % of the max fODF amplitude (here, 10% by default). This is similar to peak detection where peaks below 0.1 amplitude are usually considered noise peaks. Because SDT is based on a q-ball ODF deconvolution, and not signal deconvolution, using the max instead of mean (as in CSD), is more stable. convergenceint Maximum number of iterations to allow the deconvolution to converge. Pndarray This is an optimization to avoid computing dot(X.T, X) many times. If the same X is used many times, P can be precomputed and passed to this function. Returns fodf_shndarray ((sh_order + 1)*(sh_order + 2)/2,) Spherical harmonics coefficients of the constrained-regularized fiber ODF. num_itint Number of iterations in the constrained-regularization used for convergence. Notes This section describes how the fitting of the SH coefficients is done. Problem is to minimise per iteration: $$F(f_n) = ||Xf_n - S||^2 + \lambda^2 ||H_{n-1} f_n||^2$$ Where $$X$$ maps current FOD SH coefficients $$f_n$$ to DW signals $$s$$ and $$H_{n-1}$$ maps FOD SH coefficients $$f_n$$ to amplitudes along set of negative directions identified in previous iteration, i.e. the matrix formed by the rows of $$B_{reg}$$ for which $$Hf_{n-1}<0$$ where $$B_{reg}$$ maps $$f_n$$ to FOD amplitude on a sphere. Solve by differentiating and setting to zero: $$\Rightarrow \frac{\delta F}{\delta f_n} = 2X^T(Xf_n - S) + 2 \lambda^2 H_{n-1}^TH_{n-1}f_n=0$$ Or: $$(X^TX + \lambda^2 H_{n-1}^TH_{n-1})f_n = X^Ts$$ Define $$Q = X^TX + \lambda^2 H_{n-1}^TH_{n-1}$$ , which by construction is a square positive definite symmetric matrix of size $$n_{SH} by n_{SH}$$. If needed, positive definiteness can be enforced with a small minimum norm regulariser (helps a lot with poorly conditioned direction sets and/or superresolution): $$Q = X^TX + (\lambda H_{n-1}^T) (\lambda H_{n-1}) + \mu I$$ Solve $$Qf_n = X^Ts$$ using Cholesky decomposition: $$Q = LL^T$$ where $$L$$ is lower triangular. Then problem can be solved by back-substitution: $$L_y = X^Ts$$ $$L^Tf_n = y$$ To speeds things up further, form $$P = X^TX + \mu I$$, and update to form $$Q$$ by rankn update with $$H_{n-1}$$. The dipy implementation looks like: form initially $$P = X^T X + \mu I$$ and $$\lambda B_{reg}$$ for each voxel: form $$z = X^Ts$$ estimate $$f_0$$ by solving $$Pf_0=z$$. We use a simplified $$l_{max}=4$$ solution here, but it might not make a big difference. Then iterate until no change in rows of $$H$$ used in $$H_n$$ form $$H_{n}$$ given $$f_{n-1}$$ form $$Q = P + (\lambda H_{n-1}^T) (\lambda H_{n-1}$$) (this can be done by rankn update, but we currently do not use rankn update). solve $$Qf_n = z$$ using Cholesky decomposition We’d like to thanks Donald Tournier for his help with describing and implementing this algorithm. References 1(1,2,3) Tournier, J.D., et al. NeuroImage 2007. Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution. ### estimate_response dipy.reconst.csdeconv.estimate_response(gtab, evals, S0) Estimate single fiber response function Parameters gtabGradientTable evalsndarray S0float non diffusion weighted Returns Sestimated signal ### fa_inferior dipy.reconst.csdeconv.fa_inferior(FA, fa_thr) Check that the FA is lower than the FA threshold Parameters FAarray Fractional Anisotropy fa_thrint FA threshold Returns True when the FA value is lower than the FA threshold, otherwise False. ### fa_superior dipy.reconst.csdeconv.fa_superior(FA, fa_thr) Check that the FA is greater than the FA threshold Parameters FAarray Fractional Anisotropy fa_thrint FA threshold Returns True when the FA value is greater than the FA threshold, otherwise False. ### fa_trace_to_lambdas dipy.reconst.csdeconv.fa_trace_to_lambdas(fa=0.08, trace=0.0021) ### forward_sdeconv_mat dipy.reconst.csdeconv.forward_sdeconv_mat(r_rh, n) Build forward spherical deconvolution matrix Parameters r_rhndarray Rotational harmonics coefficients for the single fiber response function. Each element rh[i] is associated with spherical harmonics of degree 2*i. nndarray The degree of spherical harmonic function associated with each row of the deconvolution matrix. Only even degrees are allowed Returns Rndarray (N, N) Deconvolution matrix with shape (N, N) ### forward_sdt_deconv_mat dipy.reconst.csdeconv.forward_sdt_deconv_mat(ratio, n, r2_term=False) Build forward sharpening deconvolution transform (SDT) matrix Parameters ratiofloat ratio = :math: rac{lambda_2}{lambda_1} of the single fiber response function nndarray (N,) The degree of spherical harmonic function associated with each row of the deconvolution matrix. Only even degrees are allowed. r2_termbool True if ODF comes from an ODF computed from a model using the $$r^2$$ term in the integral. For example, DSI, GQI, SHORE, CSA, Tensor, Multi-tensor ODFs. This results in using the proper analytical response function solution solving from the single-fiber ODF with the r^2 term. This derivation is not published anywhere but is very similar to [Rf897e7c36096-1]. Returns Rndarray (N, N) SDT deconvolution matrix Pndarray (N, N) Funk-Radon Transform (FRT) matrix ### fractional_anisotropy dipy.reconst.csdeconv.fractional_anisotropy(evals, axis=-1) Fractional anisotropy (FA) of a diffusion tensor. Parameters evalsarray-like Eigenvalues of a diffusion tensor. axisint Axis of evals which contains 3 eigenvalues. Returns faarray Calculated FA. Range is 0 <= FA <= 1. Notes FA is calculated using the following equation: $FA = \sqrt{\frac{1}{2}\frac{(\lambda_1-\lambda_2)^2+(\lambda_1- \lambda_3)^2+(\lambda_2-\lambda_3)^2}{\lambda_1^2+ \lambda_2^2+\lambda_3^2}}$ ### get_sphere dipy.reconst.csdeconv.get_sphere(name='symmetric362') provide triangulated spheres Parameters namestr which sphere - one of: * ‘symmetric362’ * ‘symmetric642’ * ‘symmetric724’ * ‘repulsion724’ * ‘repulsion100’ * ‘repulsion200’ Returns spherea dipy.core.sphere.Sphere class instance Examples >>> import numpy as np >>> from dipy.data import get_sphere >>> sphere = get_sphere('symmetric362') >>> verts, faces = sphere.vertices, sphere.faces >>> verts.shape == (362, 3) True >>> faces.shape == (720, 3) True >>> verts, faces = get_sphere('not a sphere name') #doctest: +IGNORE_EXCEPTION_DETAIL Traceback (most recent call last): ... DataError: No sphere called "not a sphere name"  ### lazy_index dipy.reconst.csdeconv.lazy_index(index) Produces a lazy index Returns a slice that can be used for indexing an array, if no slice can be made index is returned as is. ### lpn dipy.reconst.csdeconv.lpn(n, z) Legendre function of the first kind. Compute sequence of Legendre functions of the first kind (polynomials), Pn(z) and derivatives for all degrees from 0 to n (inclusive). See also special.legendre for polynomial class. References 1 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 ### multi_voxel_fit dipy.reconst.csdeconv.multi_voxel_fit(single_voxel_fit) Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition ### ndindex dipy.reconst.csdeconv.ndindex(shape) An N-dimensional iterator object to index arrays. Given the shape of an array, an ndindex instance iterates over the N-dimensional index of the array. At each iteration a tuple of indices is returned; the last dimension is iterated over first. Parameters shapetuple of ints The dimensions of the array. Examples >>> from dipy.core.ndindex import ndindex >>> shape = (3, 2, 1) >>> for index in ndindex(shape): ... print(index) (0, 0, 0) (0, 1, 0) (1, 0, 0) (1, 1, 0) (2, 0, 0) (2, 1, 0)  ### odf_deconv dipy.reconst.csdeconv.odf_deconv(odf_sh, R, B_reg, lambda_=1.0, tau=0.1, r2_term=False) ODF constrained-regularized spherical deconvolution using the Sharpening Deconvolution Transform (SDT) [1], [2]. Parameters odf_shndarray ((sh_order + 1)*(sh_order + 2)/2,) ndarray of SH coefficients for the ODF spherical function to be deconvolved Rndarray ((sh_order + 1)(sh_order + 2)/2, (sh_order + 1)(sh_order + 2)/2) SDT matrix in SH basis B_regndarray ((sh_order + 1)(sh_order + 2)/2, (sh_order + 1)(sh_order + 2)/2) SH basis matrix used for deconvolution lambda_float lambda parameter in minimization equation (default 1.0) taufloat threshold (tau *max(fODF)) controlling the amplitude below which the corresponding fODF is assumed to be zero. r2_termbool True if ODF is computed from model that uses the $$r^2$$ term in the integral. Recall that Tuch’s ODF (used in Q-ball Imaging [1]) and the true normalized ODF definition differ from a $$r^2$$ term in the ODF integral. The original Sharpening Deconvolution Transform (SDT) technique [2] is expecting Tuch’s ODF without the $$r^2$$ (see [3] for the mathematical details). Now, this function supports ODF that have been computed using the $$r^2$$ term because the proper analytical response function has be derived. For example, models such as DSI, GQI, SHORE, CSA, Tensor, Multi-tensor ODFs, should now be deconvolved with the r2_term=True. Returns fodf_shndarray ((sh_order + 1)(sh_order + 2)/2,) Spherical harmonics coefficients of the constrained-regularized fiber ODF num_itint Number of iterations in the constrained-regularization used for convergence References 1(1,2,3,4) Tuch, D. MRM 2004. Q-Ball Imaging. 2(1,2,3,4) Descoteaux, M., et al. IEEE TMI 2009. Deterministic and Probabilistic Tractography Based on Complex Fibre Orientation Distributions 3(1,2) Descoteaux, M, PhD thesis, INRIA Sophia-Antipolis, 2008. ### odf_sh_to_sharp dipy.reconst.csdeconv.odf_sh_to_sharp(odfs_sh, sphere, basis=None, ratio=0.2, sh_order=8, lambda_=1.0, tau=0.1, r2_term=False) Sharpen odfs using the sharpening deconvolution transform [2] This function can be used to sharpen any smooth ODF spherical function. In theory, this should only be used to sharpen QballModel ODFs, but in practice, one can play with the deconvolution ratio and sharpen almost any ODF-like spherical function. The constrained-regularization is stable and will not only sharpen the ODF peaks but also regularize the noisy peaks. Parameters odfs_shndarray ((sh_order + 1)*(sh_order + 2)/2, ) array of odfs expressed as spherical harmonics coefficients sphereSphere sphere used to build the regularization matrix basis{None, ‘tournier07’, ‘descoteaux07’} different spherical harmonic basis: None for the default DIPY basis, tournier07 for the Tournier 2007 [4] basis, and descoteaux07 for the Descoteaux 2007 [3] basis (None defaults to descoteaux07). ratiofloat, ratio of the smallest vs the largest eigenvalue of the single prolate tensor response function ($$\frac{\lambda_2}{\lambda_1}$$) sh_orderint maximal SH order of the SH representation lambda_float lambda parameter (see odfdeconv) (default 1.0) taufloat tau parameter in the L matrix construction (see odfdeconv) (default 0.1) r2_termbool True if ODF is computed from model that uses the $$r^2$$ term in the integral. Recall that Tuch’s ODF (used in Q-ball Imaging [1]) and the true normalized ODF definition differ from a $$r^2$$ term in the ODF integral. The original Sharpening Deconvolution Transform (SDT) technique [2] is expecting Tuch’s ODF without the $$r^2$$ (see [3] for the mathematical details). Now, this function supports ODF that have been computed using the $$r^2$$ term because the proper analytical response function has be derived. For example, models such as DSI, GQI, SHORE, CSA, Tensor, Multi-tensor ODFs, should now be deconvolved with the r2_term=True. Returns fodf_shndarray sharpened odf expressed as spherical harmonics coefficients References 1(1,2) Tuch, D. MRM 2004. Q-Ball Imaging. 2(1,2,3,4) Descoteaux, M., et al. IEEE TMI 2009. Deterministic and Probabilistic Tractography Based on Complex Fibre Orientation Distributions 3(1,2,3) Descoteaux, M., Angelino, E., Fitzgibbons, S. and Deriche, R. Regularized, Fast, and Robust Analytical Q-ball Imaging. Magn. Reson. Med. 2007;58:497-510. 4(1,2) Tournier J.D., Calamante F. and Connelly A. Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution. NeuroImage. 2007;35(4):1459-1472. ### peaks_from_model dipy.reconst.csdeconv.peaks_from_model(model, data, sphere, relative_peak_threshold, min_separation_angle, mask=None, return_odf=False, return_sh=True, gfa_thr=0, normalize_peaks=False, sh_order=8, sh_basis_type=None, npeaks=5, B=None, invB=None, parallel=False, nbr_processes=None) Fit the model to data and computes peaks and metrics Parameters modela model instance model will be used to fit the data. sphereSphere The Sphere providing discrete directions for evaluation. relative_peak_thresholdfloat Only return peaks greater than relative_peak_threshold * m where m is the largest peak. min_separation_anglefloat in [0, 90] The minimum distance between directions. If two peaks are too close only the larger of the two is returned. maskarray, optional If mask is provided, voxels that are False in mask are skipped and no peaks are returned. return_odfbool If True, the odfs are returned. return_shbool If True, the odf as spherical harmonics coefficients is returned gfa_thrfloat Voxels with gfa less than gfa_thr are skipped, no peaks are returned. normalize_peaksbool If true, all peak values are calculated relative to max(odf). sh_orderint, optional Maximum SH order in the SH fit. For sh_order, there will be (sh_order + 1) * (sh_order + 2) / 2 SH coefficients (default 8). sh_basis_type{None, ‘tournier07’, ‘descoteaux07’} None for the default DIPY basis, tournier07 for the Tournier 2007 [2] basis, and descoteaux07 for the Descoteaux 2007 [1] basis (None defaults to descoteaux07). sh_smoothfloat, optional Lambda-regularization in the SH fit (default 0.0). npeaksint Maximum number of peaks found (default 5 peaks). Bndarray, optional Matrix that transforms spherical harmonics to spherical function sf = np.dot(sh, B). invBndarray, optional Inverse of B. parallel: bool If True, use multiprocessing to compute peaks and metric (default False). Temporary files are saved in the default temporary directory of the system. It can be changed using import tempfile and tempfile.tempdir = '/path/to/tempdir'. nbr_processes: int If parallel is True, the number of subprocesses to use (default multiprocessing.cpu_count()). Returns pamPeaksAndMetrics An object with gfa, peak_directions, peak_values, peak_indices, odf, shm_coeffs as attributes References 1(1,2) Descoteaux, M., Angelino, E., Fitzgibbons, S. and Deriche, R. Regularized, Fast, and Robust Analytical Q-ball Imaging. Magn. Reson. Med. 2007;58:497-510. 2(1,2) Tournier J.D., Calamante F. and Connelly A. Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution. NeuroImage. 2007;35(4):1459-1472. ### quad dipy.reconst.csdeconv.quad(func, a, b, args=(), full_output=0, epsabs=1.49e-08, epsrel=1.49e-08, limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50, limlst=50) Compute a definite integral. Integrate func from a to b (possibly infinite interval) using a technique from the Fortran library QUADPACK. Parameters func{function, scipy.LowLevelCallable} A Python function or method to integrate. If func takes many arguments, it is integrated along the axis corresponding to the first argument. If the user desires improved integration performance, then f may be a scipy.LowLevelCallable with one of the signatures: double func(double x) double func(double x, void *user_data) double func(int n, double *xx) double func(int n, double *xx, void *user_data)  The user_data is the data contained in the scipy.LowLevelCallable. In the call forms with xx, n is the length of the xx array which contains xx[0] == x and the rest of the items are numbers contained in the args argument of quad. In addition, certain ctypes call signatures are supported for backward compatibility, but those should not be used in new code. afloat Lower limit of integration (use -numpy.inf for -infinity). bfloat Upper limit of integration (use numpy.inf for +infinity). argstuple, optional Extra arguments to pass to func. full_outputint, optional Non-zero to return a dictionary of integration information. If non-zero, warning messages are also suppressed and the message is appended to the output tuple. Returns yfloat The integral of func from a to b. abserrfloat An estimate of the absolute error in the result. infodictdict A dictionary containing additional information. Run scipy.integrate.quad_explain() for more information. message A convergence message. explain Appended only with ‘cos’ or ‘sin’ weighting and infinite integration limits, it contains an explanation of the codes in infodict[‘ierlst’] Other Parameters epsabsfloat or int, optional Absolute error tolerance. epsrelfloat or int, optional Relative error tolerance. limitfloat or int, optional An upper bound on the number of subintervals used in the adaptive algorithm. points(sequence of floats,ints), optional A sequence of break points in the bounded integration interval where local difficulties of the integrand may occur (e.g., singularities, discontinuities). The sequence does not have to be sorted. weightfloat or int, optional String indicating weighting function. Full explanation for this and the remaining arguments can be found below. wvaroptional Variables for use with weighting functions. woptsoptional Optional input for reusing Chebyshev moments. maxp1float or int, optional An upper bound on the number of Chebyshev moments. limlstint, optional Upper bound on the number of cycles (>=3) for use with a sinusoidal weighting and an infinite end-point. See also dblquad double integral tplquad triple integral nquad n-dimensional integrals (uses quad recursively) fixed_quad fixed-order Gaussian quadrature quadrature adaptive Gaussian quadrature odeint ODE integrator ode ODE integrator simps integrator for sampled data romb integrator for sampled data scipy.special for coefficients and roots of orthogonal polynomials Notes Extra information for quad() inputs and outputs If full_output is non-zero, then the third output argument (infodict) is a dictionary with entries as tabulated below. For infinite limits, the range is transformed to (0,1) and the optional outputs are given with respect to this transformed range. Let M be the input argument limit and let K be infodict[‘last’]. The entries are: ‘neval’ The number of function evaluations. ‘last’ The number, K, of subintervals produced in the subdivision process. ‘alist’ A rank-1 array of length M, the first K elements of which are the left end points of the subintervals in the partition of the integration range. ‘blist’ A rank-1 array of length M, the first K elements of which are the right end points of the subintervals. ‘rlist’ A rank-1 array of length M, the first K elements of which are the integral approximations on the subintervals. ‘elist’ A rank-1 array of length M, the first K elements of which are the moduli of the absolute error estimates on the subintervals. ‘iord’ A rank-1 integer array of length M, the first L elements of which are pointers to the error estimates over the subintervals with L=K if K<=M/2+2 or L=M+1-K otherwise. Let I be the sequence infodict['iord'] and let E be the sequence infodict['elist']. Then E[I[1]], ..., E[I[L]] forms a decreasing sequence. If the input argument points is provided (i.e. it is not None), the following additional outputs are placed in the output dictionary. Assume the points sequence is of length P. ‘pts’ A rank-1 array of length P+2 containing the integration limits and the break points of the intervals in ascending order. This is an array giving the subintervals over which integration will occur. ‘level’ A rank-1 integer array of length M (=limit), containing the subdivision levels of the subintervals, i.e., if (aa,bb) is a subinterval of (pts[1], pts[2]) where pts[0] and pts[2] are adjacent elements of infodict['pts'], then (aa,bb) has level l if |bb-aa| = |pts[2]-pts[1]| * 2**(-l). ‘ndin’ A rank-1 integer array of length P+2. After the first integration over the intervals (pts[1], pts[2]), the error estimates over some of the intervals may have been increased artificially in order to put their subdivision forward. This array has ones in slots corresponding to the subintervals for which this happens. Weighting the integrand The input variables, weight and wvar, are used to weight the integrand by a select list of functions. Different integration methods are used to compute the integral with these weighting functions. The possible values of weight and the corresponding weighting functions are. weight Weight function used wvar ‘cos’ cos(w*x) wvar = w ‘sin’ sin(w*x) wvar = w ‘alg’ g(x) = ((x-a)**alpha)*((b-x)**beta) wvar = (alpha, beta) ‘alg-loga’ g(x)*log(x-a) wvar = (alpha, beta) ‘alg-logb’ g(x)*log(b-x) wvar = (alpha, beta) ‘alg-log’ g(x)*log(x-a)*log(b-x) wvar = (alpha, beta) ‘cauchy’ 1/(x-c) wvar = c wvar holds the parameter w, (alpha, beta), or c depending on the weight selected. In these expressions, a and b are the integration limits. For the ‘cos’ and ‘sin’ weighting, additional inputs and outputs are available. For finite integration limits, the integration is performed using a Clenshaw-Curtis method which uses Chebyshev moments. For repeated calculations, these moments are saved in the output dictionary: ‘momcom’ The maximum level of Chebyshev moments that have been computed, i.e., if M_c is infodict['momcom'] then the moments have been computed for intervals of length |b-a| * 2**(-l), l=0,1,...,M_c. ‘nnlog’ A rank-1 integer array of length M(=limit), containing the subdivision levels of the subintervals, i.e., an element of this array is equal to l if the corresponding subinterval is |b-a|* 2**(-l). ‘chebmo’ A rank-2 array of shape (25, maxp1) containing the computed Chebyshev moments. These can be passed on to an integration over the same interval by passing this array as the second element of the sequence wopts and passing infodict[‘momcom’] as the first element. If one of the integration limits is infinite, then a Fourier integral is computed (assuming w neq 0). If full_output is 1 and a numerical error is encountered, besides the error message attached to the output tuple, a dictionary is also appended to the output tuple which translates the error codes in the array info['ierlst'] to English messages. The output information dictionary contains the following entries instead of ‘last’, ‘alist’, ‘blist’, ‘rlist’, and ‘elist’: ‘lst’ The number of subintervals needed for the integration (call it K_f). ‘rslst’ A rank-1 array of length M_f=limlst, whose first K_f elements contain the integral contribution over the interval (a+(k-1)c, a+kc) where c = (2*floor(|w|) + 1) * pi / |w| and k=1,2,...,K_f. ‘erlst’ A rank-1 array of length M_f containing the error estimate corresponding to the interval in the same position in infodict['rslist']. ‘ierlst’ A rank-1 integer array of length M_f containing an error flag corresponding to the interval in the same position in infodict['rslist']. See the explanation dictionary (last entry in the output tuple) for the meaning of the codes. 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)  ### real_sph_harm dipy.reconst.csdeconv.real_sph_harm(m, n, theta, phi) Compute real spherical harmonics. 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. Parameters mint |m| <= n The order of the harmonic. nint >= 0 The degree of the harmonic. thetafloat [0, 2*pi] The azimuthal (longitudinal) coordinate. phifloat [0, pi] The polar (colatitudinal) coordinate. Returns y_mnreal float The real harmonic $$Y^m_n$$ sampled at theta and phi. See also scipy.special.sph_harm ### real_sym_sh_basis dipy.reconst.csdeconv.real_sym_sh_basis(sh_order, theta, phi) Samples a real symmetric spherical harmonic basis at point on the sphere 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. Parameters sh_orderint even int > 0, max spherical harmonic degree thetafloat [0, 2*pi] The azimuthal (longitudinal) coordinate. phifloat [0, pi] The polar (colatitudinal) coordinate. Returns y_mnreal float The real harmonic $$Y^m_n$$ sampled at theta and phi marray The order of the harmonics. narray The degree of the harmonics. References 1(1,2) Descoteaux, M., Angelino, E., Fitzgibbons, S. and Deriche, R. Regularized, Fast, and Robust Analytical Q-ball Imaging. Magn. Reson. Med. 2007;58:497-510. ### recursive_response dipy.reconst.csdeconv.recursive_response(gtab, data, mask=None, sh_order=8, peak_thr=0.01, init_fa=0.08, init_trace=0.0021, iter=8, convergence=0.001, parallel=True, nbr_processes=None, sphere=<dipy.core.sphere.HemiSphere object>) Recursive calibration of response function using peak threshold Parameters gtabGradientTable datandarray diffusion data maskndarray, optional mask for recursive calibration, for example a white matter mask. It has shape data.shape[0:3] and dtype=bool. Default: use the entire data array. sh_orderint, optional maximal spherical harmonics order. Default: 8 peak_thrfloat, optional peak threshold, how large the second peak can be relative to the first peak in order to call it a single fiber population [1]. Default: 0.01 init_fafloat, optional FA of the initial ‘fat’ response function (tensor). Default: 0.08 init_tracefloat, optional trace of the initial ‘fat’ response function (tensor). Default: 0.0021 iterint, optional maximum number of iterations for calibration. Default: 8. convergencefloat, optional convergence criterion, maximum relative change of SH coefficients. Default: 0.001. parallelbool, optional Whether to use parallelization in peak-finding during the calibration procedure. Default: True nbr_processes: int If parallel is True, the number of subprocesses to use (default multiprocessing.cpu_count()). sphereSphere, optional. The sphere used for peak finding. Default: default_sphere. Returns responsendarray response function in SH coefficients Notes In CSD there is an important pre-processing step: the estimation of the fiber response function. Using an FA threshold is not a very robust method. It is dependent on the dataset (non-informed used subjectivity), and still depends on the diffusion tensor (FA and first eigenvector), which has low accuracy at high b-value. This function recursively calibrates the response function, for more information see [1]. References 1 Tax, C.M.W., et al. NeuroImage 2014. Recursive calibration of the fiber response function for spherical deconvolution of diffusion MRI data. ### response_from_mask dipy.reconst.csdeconv.response_from_mask(gtab, data, mask) Estimate the response function from a given mask. Parameters gtabGradientTable datandarray Diffusion data maskndarray Mask to use for the estimation of the response function. For example a mask of the white matter voxels with FA values higher than 0.7 (see [1]). Returns responsetuple, (2,) (evals, S0) ratiofloat The ratio between smallest versus largest eigenvalue of the response. Notes See csdeconv.auto_response() or csdeconv.recursive_response() if you don’t have a computed mask for the response function estimation. References 1(1,2) Tournier, J.D., et al. NeuroImage 2004. Direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution ### sh_to_rh dipy.reconst.csdeconv.sh_to_rh(r_sh, m, n) Spherical harmonics (SH) to rotational harmonics (RH) Calculate the rotational harmonic decomposition up to harmonic order m, degree n for an axially and antipodally symmetric function. Note that all m != 0 coefficients will be ignored as axial symmetry is assumed. Hence, there will be (sh_order/2 + 1) non-zero coefficients. Parameters r_shndarray (N,) ndarray of SH coefficients for the single fiber response function. These coefficients must correspond to the real spherical harmonic functions produced by shm.real_sph_harm. mndarray (N,) The order of the spherical harmonic function associated with each coefficient. nndarray (N,) The degree of the spherical harmonic function associated with each coefficient. Returns r_rhndarray ((sh_order + 1)*(sh_order + 2)/2,) Rotational harmonics coefficients representing the input r_sh See also shm.real_sph_harm, shm.real_sym_sh_basis References 1 Tournier, J.D., et al. NeuroImage 2007. Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution ### single_tensor dipy.reconst.csdeconv.single_tensor(gtab, S0=1, evals=None, evecs=None, snr=None) Simulated Q-space signal with a single tensor. Parameters gtabGradientTable Measurement directions. S0double, Strength of signal in the presence of no diffusion gradient (also called the b=0 value). evals(3,) ndarray Eigenvalues of the diffusion tensor. By default, values typical for prolate white matter are used. evecs(3, 3) ndarray Eigenvectors of the tensor. You can also think of this as a rotation matrix that transforms the direction of the tensor. The eigenvectors need to be column wise. snrfloat Signal to noise ratio, assuming Rician noise. None implies no noise. Returns S(N,) ndarray Simulated signal: S(q, tau) = S_0 e^(-b g^T R D R.T g). References 1 M. Descoteaux, “High Angular Resolution Diffusion MRI: from Local Estimation to Segmentation and Tractography”, PhD thesis, University of Nice-Sophia Antipolis, p. 42, 2008. 2 E. Stejskal and J. Tanner, “Spin diffusion measurements: spin echos in the presence of a time-dependent field gradient”, Journal of Chemical Physics, nr. 42, pp. 288–292, 1965. ### sph_harm_ind_list dipy.reconst.csdeconv.sph_harm_ind_list(sh_order) Returns the degree (n) and order (m) 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_sph_harm(). Parameters sh_orderint even int > 0, max degree to return Returns m_listarray orders of even spherical harmonics n_listarray degrees of even spherical harmonics ### vec2vec_rotmat dipy.reconst.csdeconv.vec2vec_rotmat(u, v) rotation matrix from 2 unit vectors u, v being unit 3d vectors return a 3x3 rotation matrix R than aligns u to v. In general there are many rotations that will map u to v. If S is any rotation using v as an axis then R.S will also map u to v since (S.R)u = S(Ru) = Sv = v. The rotation R returned by vec2vec_rotmat leaves fixed the perpendicular to the plane spanned by u and v. The transpose of R will align v to u. Parameters uarray, shape(3,) varray, shape(3,) Returns Rarray, shape(3,3) Examples >>> import numpy as np >>> from dipy.core.geometry import vec2vec_rotmat >>> u=np.array([1,0,0]) >>> v=np.array([0,1,0]) >>> R=vec2vec_rotmat(u,v) >>> np.dot(R,u) array([ 0., 1., 0.]) >>> np.dot(R.T,v) array([ 1., 0., 0.])  ### DiffusionKurtosisFit class dipy.reconst.dki.DiffusionKurtosisFit(model, model_params) Class for fitting the Diffusion Kurtosis Model Attributes S0_hat directions For tracking - return the primary direction in each voxel evals Returns the eigenvalues of the tensor as an array evecs Returns the eigenvectors of the tensor as an array, columnwise kt Returns the 15 independent elements of the kurtosis tensor as an array quadratic_form Calculates the 3x3 diffusion tensor for each voxel shape Methods  ad() Axial diffusivity (AD) calculated from cached eigenvalues. adc(sphere) Calculate the apparent diffusion coefficient (ADC) in each direction on ak([min_kurtosis, max_kurtosis]) Axial Kurtosis (AK) of a diffusion kurtosis tensor [1]. akc(sphere) Calculates the apparent kurtosis coefficient (AKC) in each direction on the sphere for each voxel in the data color_fa() Color fractional anisotropy of diffusion tensor fa() Fractional anisotropy (FA) calculated from cached eigenvalues. ga() Geodesic anisotropy (GA) calculated from cached eigenvalues. kmax([sphere, gtol, mask]) Computes the maximum value of a single voxel kurtosis tensor linearity() Returns md() Mean diffusivity (MD) calculated from cached eigenvalues. mk([min_kurtosis, max_kurtosis]) Computes mean Kurtosis (MK) from the kurtosis tensor. mode() Tensor mode calculated from cached eigenvalues. odf(sphere) The diffusion orientation distribution function (dODF). planarity() Returns predict(gtab[, S0]) Given a DKI model fit, predict the signal on the vertices of a gradient table rd() Radial diffusivity (RD) calculated from cached eigenvalues. rk([min_kurtosis, max_kurtosis]) Radial Kurtosis (RK) of a diffusion kurtosis tensor [1]. sphericity() Returns trace() 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 Parameters modelDiffusionKurtosisModel Class instance Class instance containing the Diffusion Kurtosis Model for the fit model_paramsndarray (x, y, z, 27) or (n, 27) All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor ak(min_kurtosis=-0.42857142857142855, max_kurtosis=10) Axial Kurtosis (AK) of a diffusion kurtosis tensor [1]. Parameters min_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, axial kurtosis values that are smaller than min_kurtosis are replaced with -3./7 (theoretical kurtosis limit for regions that consist of water confined to spherical pores [2]) max_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, axial kurtosis values that are larger than max_kurtosis are replaced with max_kurtosis. Default = 10 Returns akarray Calculated AK. References 1(1,2,3) Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. Estimation of tensors and tensor-derived measures in diffusional kurtosis imaging. Magn Reson Med. 65(3), 823-836 2(1,2) Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: Robust estimation from DW-MRI using homogeneous polynomials. Proceedings of the 8th {IEEE} International Symposium on Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. doi: 10.1109/ISBI.2011.5872402 akc(sphere) Calculates the apparent kurtosis coefficient (AKC) in each direction on the sphere for each voxel in the data Parameters sphereSphere class instance Returns akcndarray The estimates of the apparent kurtosis coefficient in every direction on the input sphere Notes For each sphere direction with coordinates $$(n_{1}, n_{2}, n_{3})$$, the calculation of AKC is done using formula: $AKC(n)=\frac{MD^{2}}{ADC(n)^{2}}\sum_{i=1}^{3}\sum_{j=1}^{3} \sum_{k=1}^{3}\sum_{l=1}^{3}n_{i}n_{j}n_{k}n_{l}W_{ijkl}$ where $$W_{ijkl}$$ are the elements of the kurtosis tensor, MD the mean diffusivity and ADC the apparent diffusion coefficent computed as: $ADC(n)=\sum_{i=1}^{3}\sum_{j=1}^{3}n_{i}n_{j}D_{ij}$ where $$D_{ij}$$ are the elements of the diffusion tensor. kmax(sphere='repulsion100', gtol=1e-05, mask=None) Computes the maximum value of a single voxel kurtosis tensor Parameters sphereSphere class instance, optional The sphere providing sample directions for the initial search of the maximum value of kurtosis. gtolfloat, optional This input is to refine kurtosis maximum under the precision of the directions sampled on the sphere class instance. The gradient of the convergence procedure must be less than gtol before successful termination. If gtol is None, fiber direction is directly taken from the initial sampled directions of the given sphere object Returns max_valuefloat kurtosis tensor maximum value property kt Returns the 15 independent elements of the kurtosis tensor as an array mk(min_kurtosis=-0.42857142857142855, max_kurtosis=10) Computes mean Kurtosis (MK) from the kurtosis tensor. Parameters min_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, mean kurtosis values that are smaller than min_kurtosis are replaced with min_kurtosis. Default = -3./7 (theoretical kurtosis limit for regions that consist of water confined to spherical pores [2]) max_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, mean kurtosis values that are larger than max_kurtosis are replaced with max_kurtosis. Default = 10 Returns mkarray Calculated MK. Notes The MK analytical solution is calculated using the following equation [1]: $\begin{split}MK=F_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{1111}+ F_1(\lambda_2,\lambda_1,\lambda_3)\hat{W}_{2222}+ F_1(\lambda_3,\lambda_2,\lambda_1)\hat{W}_{3333}+ \\ F_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233}+ F_2(\lambda_2,\lambda_1,\lambda_3)\hat{W}_{1133}+ F_2(\lambda_3,\lambda_2,\lambda_1)\hat{W}_{1122}\end{split}$ where $$\hat{W}_{ijkl}$$ are the components of the $$W$$ tensor in the coordinates system defined by the eigenvectors of the diffusion tensor $$\mathbf{D}$$ and $\begin{split}F_1(\lambda_1,\lambda_2,\lambda_3)= \frac{(\lambda_1+\lambda_2+\lambda_3)^2} {18(\lambda_1-\lambda_2)(\lambda_1-\lambda_3)} [\frac{\sqrt{\lambda_2\lambda_3}}{\lambda_1} R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\ \frac{3\lambda_1^2-\lambda_1\lambda_2-\lambda_2\lambda_3- \lambda_1\lambda_3} {3\lambda_1 \sqrt{\lambda_2 \lambda_3}} R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-1 ]\end{split}$ and $\begin{split}F_2(\lambda_1,\lambda_2,\lambda_3)= \frac{(\lambda_1+\lambda_2+\lambda_3)^2} {3(\lambda_2-\lambda_3)^2} [\frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}} R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\ \frac{2\lambda_1-\lambda_2-\lambda_3}{3\sqrt{\lambda_2 \lambda_3}} R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-2]\end{split}$ where $$R_f$$ and $$R_d$$ are the Carlson’s elliptic integrals. References 1(1,2) Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. Estimation of tensors and tensor-derived measures in diffusional kurtosis imaging. Magn Reson Med. 65(3), 823-836 2(1,2) Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: Robust estimation from DW-MRI using homogeneous polynomials. Proceedings of the 8th {IEEE} International Symposium on Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. doi: 10.1109/ISBI.2011.5872402 predict(gtab, S0=1.0) Given a DKI model fit, predict the signal on the vertices of a gradient table Parameters gtaba GradientTable class instance The gradient table for this prediction S0float or ndarray (optional) The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1 Notes The predicted signal is given by: $S(n,b)=S_{0}e^{-bD(n)+\frac{1}{6}b^{2}D(n)^{2}K(n)}$ $$\mathbf{D(n)}$$ and $$\mathbf{K(n)}$$ can be computed from the DT and KT using the following equations: $D(n)=\sum_{i=1}^{3}\sum_{j=1}^{3}n_{i}n_{j}D_{ij}$ and $K(n)=\frac{MD^{2}}{D(n)^{2}}\sum_{i=1}^{3}\sum_{j=1}^{3} \sum_{k=1}^{3}\sum_{l=1}^{3}n_{i}n_{j}n_{k}n_{l}W_{ijkl}$ where $$D_{ij}$$ and $$W_{ijkl}$$ are the elements of the second-order DT and the fourth-order KT tensors, respectively, and $$MD$$ is the mean diffusivity. rk(min_kurtosis=-0.42857142857142855, max_kurtosis=10) Radial Kurtosis (RK) of a diffusion kurtosis tensor [1]. Parameters min_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, axial kurtosis values that are smaller than min_kurtosis are replaced with -3./7 (theoretical kurtosis limit for regions that consist of water confined to spherical pores [2]) max_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, axial kurtosis values that are larger than max_kurtosis are replaced with max_kurtosis. Default = 10 Returns rkarray Calculated RK. Notes RK is calculated with the following equation: $K_{\bot} = G_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2222} + G_1(\lambda_1,\lambda_3,\lambda_2)\hat{W}_{3333} + G_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233}$ where: $G_1(\lambda_1,\lambda_2,\lambda_3)= \frac{(\lambda_1+\lambda_2+\lambda_3)^2}{18\lambda_2(\lambda_2- \lambda_3)} \left (2\lambda_2 + \frac{\lambda_3^2-3\lambda_2\lambda_3}{\sqrt{\lambda_2\lambda_3}} \right)$ and $G_2(\lambda_1,\lambda_2,\lambda_3)= \frac{(\lambda_1+\lambda_2+\lambda_3)^2}{(\lambda_2-\lambda_3)^2} \left ( \frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}}-2 \right )$ References 1(1,2,3) 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 2(1,2) Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: Robust estimation from DW-MRI using homogeneous polynomials. Proceedings of the 8th {IEEE} International Symposium on Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. doi: 10.1109/ISBI.2011.5872402 ### DiffusionKurtosisModel class dipy.reconst.dki.DiffusionKurtosisModel(gtab, fit_method='WLS', *args, **kwargs) Class for the Diffusion Kurtosis Model Methods  fit(data[, mask]) Fit method of the DKI model class predict(dki_params[, S0]) Predict a signal for this DKI model class instance given parameters. __init__(gtab, fit_method='WLS', *args, **kwargs) Diffusion Kurtosis Tensor Model [1] Parameters gtabGradientTable class instance fit_methodstr or callable str can be one of the following: ‘OLS’ or ‘ULLS’ for ordinary least squares dki.ols_fit_dki ‘WLS’ or ‘UWLLS’ for weighted ordinary least squares dki.wls_fit_dki callable has to have the signature: fit_method(design_matrix, data, *args, **kwargs) args, kwargsarguments and key-word arguments passed to the fit_method. See dki.ols_fit_dki, dki.wls_fit_dki for details References 1 Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. Estimation of tensors and tensor-derived measures in diffusional kurtosis imaging. Magn Reson Med. 65(3), 823-836 fit(data, mask=None) Fit method of the DKI model class Parameters dataarray The measured signal from one voxel. maskarray A boolean array used to mark the coordinates in the data that should be analyzed that has the shape data.shape[-1] predict(dki_params, S0=1.0) Predict a signal for this DKI model class instance given parameters. Parameters dki_paramsndarray (x, y, z, 27) or (n, 27) All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor S0float or ndarray (optional) The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1 ### ReconstModel class dipy.reconst.dki.ReconstModel(gtab) Bases: object Abstract class for signal reconstruction models Methods  fit __init__(gtab) Initialization of the abstract class for signal reconstruction models Parameters gtabGradientTable class instance fit(data, mask=None, **kwargs) ### TensorFit class dipy.reconst.dki.TensorFit(model, model_params, model_S0=None) Bases: object Attributes S0_hat directions For tracking - return the primary direction in each voxel evals Returns the eigenvalues of the tensor as an array evecs Returns the eigenvectors of the tensor as an array, columnwise quadratic_form Calculates the 3x3 diffusion tensor for each voxel shape Methods  Axial diffusivity (AD) calculated from cached eigenvalues. adc(sphere) 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. Returns Mean diffusivity (MD) calculated from cached eigenvalues. Tensor mode calculated from cached eigenvalues. odf(sphere) The diffusion orientation distribution function (dODF). Returns predict(gtab[, S0, step]) Given a model fit, predict the signal on the vertices of a sphere Radial diffusivity (RD) calculated from cached eigenvalues. Returns Trace of the tensor calculated from cached eigenvalues.  lower_triangular __init__(model, model_params, model_S0=None) Initialize a TensorFit class instance. property S0_hat ad() Axial diffusivity (AD) calculated from cached eigenvalues. Returns adarray (V, 1) Calculated AD. Notes RD is calculated with the following equation: $AD = \lambda_1$ adc(sphere) Calculate the apparent diffusion coefficient (ADC) in each direction on the sphere for each voxel in the data Parameters sphereSphere class instance Returns adcndarray The estimates of the apparent diffusion coefficient in every direction on the input sphere ec{b} Q ec{b}^T Where Q is the quadratic form of the tensor. color_fa() Color fractional anisotropy of diffusion tensor property directions For tracking - return the primary direction in each voxel property evals Returns the eigenvalues of the tensor as an array property evecs Returns the eigenvectors of the tensor as an array, columnwise fa() Fractional anisotropy (FA) calculated from cached eigenvalues. ga() Geodesic anisotropy (GA) calculated from cached eigenvalues. linearity() Returns linearityarray Calculated linearity of the diffusion tensor 1. Notes Linearity is calculated with the following equation: $Linearity = \frac{\lambda_1-\lambda_2}{\lambda_1+\lambda_2+\lambda_3}$ References [1] Westin C.-F., Peled S., Gubjartsson H., Kikinis R., Jolesz F., “Geometrical diffusion measures for MRI from tensor basis analysis” in Proc. 5th Annual ISMRM, 1997. lower_triangular(b0=None) md() Mean diffusivity (MD) calculated from cached eigenvalues. Returns mdarray (V, 1) Calculated MD. Notes MD is calculated with the following equation: $MD = \frac{\lambda_1+\lambda_2+\lambda_3}{3}$ mode() Tensor mode calculated from cached eigenvalues. odf(sphere) The diffusion orientation distribution function (dODF). This is an estimate of the diffusion distance in each direction Parameters sphereSphere class instance. The dODF is calculated in the vertices of this input. Returns odfndarray The diffusion distance in every direction of the sphere in every voxel in the input data. Notes This is based on equation 3 in [Aganj2010]. To re-derive it from scratch, follow steps in [Descoteaux2008], Section 7.9 Equation 7.24 but with an $$r^2$$ term in the integral. References Aganj2010(1,2) 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 Descoteaux2008(1,2) 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() Returns sphericityarray Calculated sphericity of the diffusion tensor 1. Notes Sphericity is calculated with the following equation: $Sphericity = \frac{2 (\lambda_2 - \lambda_3)}{\lambda_1+\lambda_2+\lambda_3}$ References [1] Westin C.-F., Peled S., Gubjartsson H., Kikinis R., Jolesz F., “Geometrical diffusion measures for MRI from tensor basis analysis” in Proc. 5th Annual ISMRM, 1997. predict(gtab, S0=None, step=None) Given a model fit, predict the signal on the vertices of a sphere Parameters gtaba GradientTable class instance This encodes the directions for which a prediction is made S0float array The mean non-diffusion weighted signal in each voxel. Default: The fitted S0 value in all voxels if it was fitted. Otherwise 1 in all voxels. stepint The chunk size as a number of voxels. Optional parameter with default value 10,000. In order to increase speed of processing, tensor fitting is done simultaneously over many voxels. This parameter sets the number of voxels that will be fit at once in each iteration. A larger step value should speed things up, but it will also take up more memory. It is advisable to keep an eye on memory consumption as this value is increased. Notes The predicted signal is given by: $S( heta, b) = S_0 * e^{-b ADC}$ Where: .. math 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 property quadratic_form Calculates the 3x3 diffusion tensor for each voxel rd() Radial diffusivity (RD) calculated from cached eigenvalues. Returns rdarray (V, 1) Calculated RD. Notes RD is calculated with the following equation: $RD = \frac{\lambda_2 + \lambda_3}{2}$ property shape sphericity() Returns sphericityarray Calculated sphericity of the diffusion tensor 1. Notes Sphericity is calculated with the following equation: $Sphericity = \frac{3 \lambda_3}{\lambda_1+\lambda_2+\lambda_3}$ References [1] Westin C.-F., Peled S., Gubjartsson H., Kikinis R., Jolesz F., “Geometrical diffusion measures for MRI from tensor basis analysis” in Proc. 5th Annual ISMRM, 1997. trace() Trace of the tensor calculated from cached eigenvalues. Returns tracearray (V, 1) Calculated trace. Notes The trace is calculated with the following equation: $trace = \lambda_1 + \lambda_2 + \lambda_3$ ### Wcons dipy.reconst.dki.Wcons(k_elements) Construct the full 4D kurtosis tensors from its 15 independent elements Parameters k_elements(15,) elements of the kurtosis tensor in the following order: .. math:: begin{matrix} ( & W_{xxxx} & W_{yyyy} & W_{zzzz} & W_{xxxy} & W_{xxxz} & … \ & W_{xyyy} & W_{yyyz} & W_{xzzz} & W_{yzzz} & W_{xxyy} & … \ & W_{xxzz} & W_{yyzz} & W_{xxyz} & W_{xyyz} & W_{xyzz} & & )end{matrix} Returns Warray(3, 3, 3, 3) Full 4D kurtosis tensor ### Wrotate dipy.reconst.dki.Wrotate(kt, Basis) Rotate a kurtosis tensor from the standard Cartesian coordinate system to another coordinate system basis Parameters kt(15,) Vector with the 15 independent elements of the kurtosis tensor Basisarray (3, 3) Vectors of the basis column-wise oriented indsarray(m, 4) (optional) Array of vectors containing the four indexes of m specific elements of the rotated kurtosis tensor. If not specified all 15 elements of the rotated kurtosis tensor are computed. Returns Wrotarray (m,) or (15,) Vector with the m independent elements of the rotated kurtosis tensor. If ‘indices’ is not specified all 15 elements of the rotated kurtosis tensor are computed. Notes KT elements are assumed to be ordered as follows:  begin{matrix} ( & W_{xxxx} & W_{yyyy} & W_{zzzz} & W_{xxxy} & W_{xxxz} & … \ & W_{xyyy} & W_{yyyz} & W_{xzzz} & W_{yzzz} & W_{xxyy} & … \ & W_{xxzz} & W_{yyzz} & W_{xxyz} & W_{xyyz} & W_{xyzz} & & )end{matrix} References [1] Hui ES, Cheung MM, Qi L, Wu EX, 2008. Towards better MR characterization of neural tissues using directional diffusion kurtosis analysis. Neuroimage 42(1): 122-34 ### Wrotate_element dipy.reconst.dki.Wrotate_element(kt, indi, indj, indk, indl, B) Computes the the specified index element of a kurtosis tensor rotated to the coordinate system basis B. Parameters ktndarray (x, y, z, 15) or (n, 15) Array containing the 15 independent elements of the kurtosis tensor indiint Rotated kurtosis tensor element index i (0 for x, 1 for y, 2 for z) indjint Rotated kurtosis tensor element index j (0 for x, 1 for y, 2 for z) indkint Rotated kurtosis tensor element index k (0 for x, 1 for y, 2 for z) indl: int Rotated kurtosis tensor element index l (0 for x, 1 for y, 2 for z) B: array (x, y, z, 3, 3) or (n, 15) Vectors of the basis column-wise oriented Returns Wrefloat rotated kurtosis tensor element of index ind_i, ind_j, ind_k, ind_l Notes It is assumed that initial kurtosis tensor elementes are defined on the Cartesian coordinate system. References [1] Hui ES, Cheung MM, Qi L, Wu EX, 2008. Towards better MR characterization of neural tissues using directional diffusion kurtosis analysis. Neuroimage 42(1): 122-34 ### apparent_kurtosis_coef dipy.reconst.dki.apparent_kurtosis_coef(dki_params, sphere, min_diffusivity=0, min_kurtosis=-0.42857142857142855) Calculates the apparent kurtosis coefficient (AKC) in each direction of a sphere [1]. Parameters dki_paramsndarray (x, y, z, 27) or (n, 27) All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvectors respectively 3. Fifteen elements of the kurtosis tensor spherea Sphere class instance The AKC will be calculated for each of the vertices in the sphere min_diffusivityfloat (optional) Because negative eigenvalues are not physical and small eigenvalues cause quite a lot of noise in diffusion-based metrics, diffusivity values smaller than min_diffusivity are replaced with min_diffusivity. Default = 0 min_kurtosisfloat (optional) Because high-amplitude negative values of kurtosis are not physicaly and biologicaly pluasible, and these cause artefacts in kurtosis-based measures, directional kurtosis values smaller than min_kurtosis are replaced with min_kurtosis. Default = -3./7 (theoretical kurtosis limit for regions that consist of water confined to spherical pores [2]) Returns akcndarray (x, y, z, g) or (n, g) Apparent kurtosis coefficient (AKC) for all g directions of a sphere. Notes For each sphere direction with coordinates $$(n_{1}, n_{2}, n_{3})$$, the calculation of AKC is done using formula [1]: $AKC(n)=\frac{MD^{2}}{ADC(n)^{2}}\sum_{i=1}^{3}\sum_{j=1}^{3} \sum_{k=1}^{3}\sum_{l=1}^{3}n_{i}n_{j}n_{k}n_{l}W_{ijkl}$ where $$W_{ijkl}$$ are the elements of the kurtosis tensor, MD the mean diffusivity and ADC the apparent diffusion coefficent computed as: $ADC(n)=\sum_{i=1}^{3}\sum_{j=1}^{3}n_{i}n_{j}D_{ij}$ where $$D_{ij}$$ are the elements of the diffusion tensor. References 1(1,2,3,4) Neto Henriques R, Correia MM, Nunes RG, Ferreira HA (2015). Exploring the 3D geometry of the diffusion kurtosis tensor - Impact on the development of robust tractography procedures and novel biomarkers, NeuroImage 111: 85-99 2(1,2) Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: Robust estimation from DW-MRI using homogeneous polynomials. Proceedings of the 8th {IEEE} International Symposium on Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. doi: 10.1109/ISBI.2011.5872402 ### axial_kurtosis dipy.reconst.dki.axial_kurtosis(dki_params, min_kurtosis=-0.42857142857142855, max_kurtosis=10) Computes axial Kurtosis (AK) from the kurtosis tensor. Parameters dki_paramsndarray (x, y, z, 27) or (n, 27) All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor min_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, axial kurtosis values that are smaller than min_kurtosis are replaced with min_kurtosis. Default = -3./7 (theoretical kurtosis limit for regions that consist of water confined to spherical pores [1]) max_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, axial kurtosis values that are larger than max_kurtosis are replaced with max_kurtosis. Default = 10 Returns akarray Calculated AK. References 1(1,2) Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: Robust estimation from DW-MRI using homogeneous polynomials. Proceedings of the 8th {IEEE} International Symposium on Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. doi: 10.1109/ISBI.2011.5872402 ### carlson_rd dipy.reconst.dki.carlson_rd(x, y, z, errtol=0.0001) Computes the Carlson’s incomplete elliptic integral of the second kind defined as: $R_D = \frac{3}{2} \int_{0}^{\infty} (t+x)^{-\frac{1}{2}} (t+y)^{-\frac{1}{2}}(t+z) ^{-\frac{3}{2}}$ Parameters xndarray First independent variable of the integral. yndarray Second independent variable of the integral. zndarray Third independent variable of the integral. errtolfloat Error tolerance. Integral is computed with relative error less in magnitude than the defined value Returns RDndarray Value of the incomplete second order elliptic integral Notes x, y, and z have to be nonnegative and at most x or y is zero. ### carlson_rf dipy.reconst.dki.carlson_rf(x, y, z, errtol=0.0003) Computes the Carlson’s incomplete elliptic integral of the first kind defined as: $R_F = \frac{1}{2} \int_{0}^{\infty} \left [(t+x)(t+y)(t+z) \right ] ^{-\frac{1}{2}}dt$ Parameters xndarray First independent variable of the integral. yndarray Second independent variable of the integral. zndarray Third independent variable of the integral. errtolfloat Error tolerance. Integral is computed with relative error less in magnitude than the defined value Returns RFndarray Value of the incomplete first order elliptic integral Notes x, y, and z have to be nonnegative and at most one of them is zero. References 1 Carlson, B.C., 1994. Numerical computation of real or complex elliptic integrals. arXiv:math/9409227 [math.CA] ### cart2sphere dipy.reconst.dki.cart2sphere(x, y, z) Return angles for Cartesian 3D coordinates x, y, and z See doc for sphere2cart for angle conventions and derivation of the formulae. $$0\le\theta\mathrm{(theta)}\le\pi$$ and $$-\pi\le\phi\mathrm{(phi)}\le\pi$$ Parameters xarray_like x coordinate in Cartesian space yarray_like y coordinate in Cartesian space zarray_like z coordinate Returns rarray radius thetaarray inclination (polar) angle phiarray azimuth angle ### check_multi_b dipy.reconst.dki.check_multi_b(gtab, n_bvals, non_zero=True, bmag=None) Check if you have enough different b-values in your gradient table Parameters gtabGradientTable class instance. n_bvalsint The number of different b-values you are checking for. non_zerobool Whether to check only non-zero bvalues. In this case, we will require at least n_bvals non-zero b-values (where non-zero is defined depending on the gtab object’s b0_threshold attribute) bmagint The order of magnitude of the b-values used. The function will normalize the b-values relative $$10^{bmag}$$. Default: derive this value from the maximal b-value provided: $$bmag=log_{10}(max(bvals)) - 1$$. Returns boolWhether there are at least n_bvals different b-values in the gradient table used. ### decompose_tensor dipy.reconst.dki.decompose_tensor(tensor, min_diffusivity=0) Returns eigenvalues and eigenvectors given a diffusion tensor Computes tensor eigen decomposition to calculate eigenvalues and eigenvectors (Basser et al., 1994a). Parameters tensorarray (…, 3, 3) Hermitian matrix representing a diffusion tensor. min_diffusivityfloat Because negative eigenvalues are not physical and small eigenvalues, much smaller than the diffusion weighting, cause quite a lot of noise in metrics such as fa, diffusivity values smaller than min_diffusivity are replaced with min_diffusivity. Returns eigvalsarray (…, 3) Eigenvalues from eigen decomposition of the tensor. Negative eigenvalues are replaced by zero. Sorted from largest to smallest. eigvecsarray (…, 3, 3) Associated eigenvectors from eigen decomposition of the tensor. Eigenvectors are columnar (e.g. eigvecs[…, :, j] is associated with eigvals[…, j]) ### design_matrix dipy.reconst.dki.design_matrix(gtab) Construct B design matrix for DKI. Parameters gtabGradientTable Measurement directions. Returns Barray (N, 22) Design matrix or B matrix for the DKI model B[j, :] = (Bxx, Bxy, Bzz, Bxz, Byz, Bzz, Bxxxx, Byyyy, Bzzzz, Bxxxy, Bxxxz, Bxyyy, Byyyz, Bxzzz, Byzzz, Bxxyy, Bxxzz, Byyzz, Bxxyz, Bxyyz, Bxyzz, BlogS0) ### directional_diffusion dipy.reconst.dki.directional_diffusion(dt, V, min_diffusivity=0) Calculates the apparent diffusion coefficient (adc) in each direction of a sphere for a single voxel [1]. Parameters dtarray (6,) elements of the diffusion tensor of the voxel. Varray (g, 3) g directions of a Sphere in Cartesian coordinates min_diffusivityfloat (optional) Because negative eigenvalues are not physical and small eigenvalues cause quite a lot of noise in diffusion-based metrics, diffusivity values smaller than min_diffusivity are replaced with min_diffusivity. Default = 0 Returns adcndarray (g,) Apparent diffusion coefficient (adc) in all g directions of a sphere for a single voxel. References 1(1,2,3) Neto Henriques R, Correia MM, Nunes RG, Ferreira HA (2015). Exploring the 3D geometry of the diffusion kurtosis tensor - Impact on the development of robust tractography procedures and novel biomarkers, NeuroImage 111: 85-99 ### directional_diffusion_variance dipy.reconst.dki.directional_diffusion_variance(kt, V, min_kurtosis=-0.42857142857142855) Calculates the apparent diffusion variance (adv) in each direction of a sphere for a single voxel [1]. Parameters dtarray (6,) elements of the diffusion tensor of the voxel. ktarray (15,) elements of the kurtosis tensor of the voxel. Varray (g, 3) g directions of a Sphere in Cartesian coordinates min_kurtosisfloat (optional) Because high-amplitude negative values of kurtosis are not physicaly and biologicaly pluasible, and these cause artefacts in kurtosis-based measures, directional kurtosis values smaller than min_kurtosis are replaced with min_kurtosis. Default = -3./7 (theoretical kurtosis limit for regions that consist of water confined to spherical pores [2]_) adcndarray(g,) (optional) Apparent diffusion coefficient (adc) in all g directions of a sphere for a single voxel. advndarray(g,) (optional) Apparent diffusion variance coefficient (advc) in all g directions of a sphere for a single voxel. Returns advndarray (g,) Apparent diffusion variance (adv) in all g directions of a sphere for a single voxel. References 1(1,2,3) Neto Henriques R, Correia MM, Nunes RG, Ferreira HA (2015). Exploring the 3D geometry of the diffusion kurtosis tensor - Impact on the development of robust tractography procedures and novel biomarkers, NeuroImage 111: 85-99 ### directional_kurtosis dipy.reconst.dki.directional_kurtosis(dt, md, kt, V, min_diffusivity=0, min_kurtosis=-0.42857142857142855, adc=None, adv=None) Calculates the apparent kurtosis coefficient (akc) in each direction of a sphere for a single voxel [1]. Parameters dtarray (6,) elements of the diffusion tensor of the voxel. mdfloat mean diffusivity of the voxel ktarray (15,) elements of the kurtosis tensor of the voxel. Varray (g, 3) g directions of a Sphere in Cartesian coordinates min_diffusivityfloat (optional) Because negative eigenvalues are not physical and small eigenvalues cause quite a lot of noise in diffusion-based metrics, diffusivity values smaller than min_diffusivity are replaced with min_diffusivity. Default = 0 min_kurtosisfloat (optional) Because high-amplitude negative values of kurtosis are not physicaly and biologicaly pluasible, and these cause artefacts in kurtosis-based measures, directional kurtosis values smaller than min_kurtosis are replaced with min_kurtosis. Default = -3./7 (theoretical kurtosis limit for regions that consist of water confined to spherical pores [2]) adcndarray(g,) (optional) Apparent diffusion coefficient (adc) in all g directions of a sphere for a single voxel. advndarray(g,) (optional) Apparent diffusion variance (advc) in all g directions of a sphere for a single voxel. Returns akcndarray (g,) Apparent kurtosis coefficient (AKC) in all g directions of a sphere for a single voxel. References 1(1,2,3) Neto Henriques R, Correia MM, Nunes RG, Ferreira HA (2015). Exploring the 3D geometry of the diffusion kurtosis tensor - Impact on the development of robust tractography procedures and novel biomarkers, NeuroImage 111: 85-99 2(1,2) Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: Robust estimation from DW-MRI using homogeneous polynomials. Proceedings of the 8th {IEEE} International Symposium on Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. doi: 10.1109/ISBI.2011.5872402 ### dki_prediction dipy.reconst.dki.dki_prediction(dki_params, gtab, S0=1.0) Predict a signal given diffusion kurtosis imaging parameters. Parameters dki_paramsndarray (x, y, z, 27) or (n, 27) All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor gtaba GradientTable class instance The gradient table for this prediction S0float or ndarray (optional) The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1 Returns S(…, N) ndarray Simulated signal based on the DKI model: $S=S_{0}e^{-bD+$ rac{1}{6}b^{2}D^{2}K} ### from_lower_triangular dipy.reconst.dki.from_lower_triangular(D) Returns a tensor given the six unique tensor elements Given the six unique tensor elements (in the order: Dxx, Dxy, Dyy, Dxz, Dyz, Dzz) returns a 3 by 3 tensor. All elements after the sixth are ignored. Parameters Darray_like, (…, >6) Unique elements of the tensors Returns tensorndarray (…, 3, 3) 3 by 3 tensors ### get_sphere dipy.reconst.dki.get_sphere(name='symmetric362') provide triangulated spheres Parameters namestr which sphere - one of: * ‘symmetric362’ * ‘symmetric642’ * ‘symmetric724’ * ‘repulsion724’ * ‘repulsion100’ * ‘repulsion200’ Returns spherea dipy.core.sphere.Sphere class instance Examples >>> import numpy as np >>> from dipy.data import get_sphere >>> sphere = get_sphere('symmetric362') >>> verts, faces = sphere.vertices, sphere.faces >>> verts.shape == (362, 3) True >>> faces.shape == (720, 3) True >>> verts, faces = get_sphere('not a sphere name') #doctest: +IGNORE_EXCEPTION_DETAIL Traceback (most recent call last): ... DataError: No sphere called "not a sphere name"  ### kurtosis_maximum dipy.reconst.dki.kurtosis_maximum(dki_params, sphere='repulsion100', gtol=0.01, mask=None) Computes kurtosis maximum value Parameters dki_paramsndarray (x, y, z, 27) or (n, 27) All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows: 1. Three diffusion tensor’s eingenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor sphereSphere class instance, optional The sphere providing sample directions for the initial search of the maximal value of kurtosis. gtolfloat, optional This input is to refine kurtosis maximum under the precision of the directions sampled on the sphere class instance. The gradient of the convergence procedure must be less than gtol before successful termination. If gtol is None, fiber direction is directly taken from the initial sampled directions of the given sphere object maskndarray A boolean array used to mark the coordinates in the data that should be analyzed that has the shape dki_params.shape[:-1] Returns max_valuefloat kurtosis tensor maximum value max_dirarray (3,) Cartesian coordinates of the direction of the maximal kurtosis value ### local_maxima dipy.reconst.dki.local_maxima() Local maxima of a function evaluated on a discrete set of points. If a function is evaluated on some set of points where each pair of neighboring points is an edge in edges, find the local maxima. Parameters odfarray, 1d, dtype=double The function evaluated on a set of discrete points. edgesarray (N, 2) The set of neighbor relations between the points. Every edge, ie edges[i, :], is a pair of neighboring points. Returns peak_valuesndarray Value of odf at a maximum point. Peak values is sorted in descending order. peak_indicesndarray Indices of maximum points. Sorted in the same order as peak_values so odf[peak_indices[i]] == peak_values[i]. Notes A point is a local maximum if it is > at least one neighbor and >= all neighbors. If no points meet the above criteria, 1 maximum is returned such that odf[maximum] == max(odf). ### lower_triangular dipy.reconst.dki.lower_triangular(tensor, b0=None) Returns the six lower triangular values of the tensor and a dummy variable if b0 is not None Parameters tensorarray_like (…, 3, 3) a collection of 3, 3 diffusion tensors b0float if b0 is not none log(b0) is returned as the dummy variable Returns Dndarray If b0 is none, then the shape will be (…, 6) otherwise (…, 7) ### mean_diffusivity dipy.reconst.dki.mean_diffusivity(evals, axis=-1) Mean Diffusivity (MD) of a diffusion tensor. Parameters evalsarray-like Eigenvalues of a diffusion tensor. axisint Axis of evals which contains 3 eigenvalues. Returns mdarray Calculated MD. Notes MD is calculated with the following equation: $MD = \frac{\lambda_1 + \lambda_2 + \lambda_3}{3}$ ### mean_kurtosis dipy.reconst.dki.mean_kurtosis(dki_params, min_kurtosis=-0.42857142857142855, max_kurtosis=3) Computes mean Kurtosis (MK) from the kurtosis tensor [1]. Parameters dki_paramsndarray (x, y, z, 27) or (n, 27) All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor min_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, mean kurtosis values that are smaller than min_kurtosis are replaced with min_kurtosis. Default = -3./7 (theoretical kurtosis limit for regions that consist of water confined to spherical pores [2]) max_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, mean kurtosis values that are larger than max_kurtosis are replaced with max_kurtosis. Default = 10 Returns mkarray Calculated MK. Notes The MK analytical solution is calculated using the following equation [1]: $\begin{split}MK=F_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{1111}+ F_1(\lambda_2,\lambda_1,\lambda_3)\hat{W}_{2222}+ F_1(\lambda_3,\lambda_2,\lambda_1)\hat{W}_{3333}+ \\ F_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233}+ F_2(\lambda_2,\lambda_1,\lambda_3)\hat{W}_{1133}+ F_2(\lambda_3,\lambda_2,\lambda_1)\hat{W}_{1122}\end{split}$ where $$\hat{W}_{ijkl}$$ are the components of the $$W$$ tensor in the coordinates system defined by the eigenvectors of the diffusion tensor $$\mathbf{D}$$ and \begin{align}\begin{aligned}\begin{split}F_1(\lambda_1,\lambda_2,\lambda_3)= \frac{(\lambda_1+\lambda_2+\lambda_3)^2} {18(\lambda_1-\lambda_2)(\lambda_1-\lambda_3)} [\frac{\sqrt{\lambda_2\lambda_3}}{\lambda_1} R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\ \frac{3\lambda_1^2-\lambda_1\lambda_2-\lambda_2\lambda_3- \lambda_1\lambda_3} {3\lambda_1 \sqrt{\lambda_2 \lambda_3}} R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-1 ]\end{split}\\\begin{split}F_2(\lambda_1,\lambda_2,\lambda_3)= \frac{(\lambda_1+\lambda_2+\lambda_3)^2} {3(\lambda_2-\lambda_3)^2} [\frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}} R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\ \frac{2\lambda_1-\lambda_2-\lambda_3}{3\sqrt{\lambda_2 \lambda_3}} R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-2]\end{split}\end{aligned}\end{align} where $$R_f$$ and $$R_d$$ are the Carlson’s elliptic integrals. References 1(1,2,3,4) Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. Estimation of tensors and tensor-derived measures in diffusional kurtosis imaging. Magn Reson Med. 65(3), 823-836 2(1,2) Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: Robust estimation from DW-MRI using homogeneous polynomials. Proceedings of the 8th {IEEE} International Symposium on Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. doi: 10.1109/ISBI.2011.5872402 ### ndindex dipy.reconst.dki.ndindex(shape) An N-dimensional iterator object to index arrays. Given the shape of an array, an ndindex instance iterates over the N-dimensional index of the array. At each iteration a tuple of indices is returned; the last dimension is iterated over first. Parameters shapetuple of ints The dimensions of the array. Examples >>> from dipy.core.ndindex import ndindex >>> shape = (3, 2, 1) >>> for index in ndindex(shape): ... print(index) (0, 0, 0) (0, 1, 0) (1, 0, 0) (1, 1, 0) (2, 0, 0) (2, 1, 0)  ### ols_fit_dki dipy.reconst.dki.ols_fit_dki(design_matrix, data) Computes ordinary least squares (OLS) fit to calculate the diffusion tensor and kurtosis tensor using a linear regression diffusion kurtosis model 1. Parameters design_matrixarray (g, 22) Design matrix holding the covariants used to solve for the regression coefficients. dataarray (N, g) Data or response variables holding the data. Note that the last dimension should contain the data. It makes no copies of data. Returns dki_paramsarray (N, 27) All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor ### radial_kurtosis dipy.reconst.dki.radial_kurtosis(dki_params, min_kurtosis=-0.42857142857142855, max_kurtosis=10) Radial Kurtosis (RK) of a diffusion kurtosis tensor [1]. Parameters dki_paramsndarray (x, y, z, 27) or (n, 27) All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor min_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, radial kurtosis values that are smaller than min_kurtosis are replaced with min_kurtosis. Default = -3./7 (theoretical kurtosis limit for regions that consist of water confined to spherical pores [2]) max_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, radial kurtosis values that are larger than max_kurtosis are replaced with max_kurtosis. Default = 10 Returns rkarray Calculated RK. Notes RK is calculated with the following equation [1]: $K_{\bot} = G_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2222} + G_1(\lambda_1,\lambda_3,\lambda_2)\hat{W}_{3333} + G_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233}$ where: $G_1(\lambda_1,\lambda_2,\lambda_3)= \frac{(\lambda_1+\lambda_2+\lambda_3)^2}{18\lambda_2(\lambda_2- \lambda_3)} \left (2\lambda_2 + \frac{\lambda_3^2-3\lambda_2\lambda_3}{\sqrt{\lambda_2\lambda_3}} \right)$ and $G_2(\lambda_1,\lambda_2,\lambda_3)= \frac{(\lambda_1+\lambda_2+\lambda_3)^2}{(\lambda_2-\lambda_3)^2} \left ( \frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}}-2\right )$ References 1(1,2,3,4) Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. Estimation of tensors and tensor-derived measures in diffusional kurtosis imaging. Magn Reson Med. 65(3), 823-836 2(1,2) Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: Robust estimation from DW-MRI using homogeneous polynomials. Proceedings of the 8th {IEEE} International Symposium on Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. doi: 10.1109/ISBI.2011.5872402 ### sphere2cart dipy.reconst.dki.sphere2cart(r, theta, phi) Spherical to Cartesian coordinates This is the standard physics convention where theta is the inclination (polar) angle, and phi is the azimuth angle. Imagine a sphere with center (0,0,0). Orient it with the z axis running south-north, the y axis running west-east and the x axis from posterior to anterior. theta (the inclination angle) is the angle to rotate from the z-axis (the zenith) around the y-axis, towards the x axis. Thus the rotation is counter-clockwise from the point of view of positive y. phi (azimuth) gives the angle of rotation around the z-axis towards the y axis. The rotation is counter-clockwise from the point of view of positive z. Equivalently, given a point P on the sphere, with coordinates x, y, z, theta is the angle between P and the z-axis, and phi is the angle between the projection of P onto the XY plane, and the X axis. Geographical nomenclature designates theta as ‘co-latitude’, and phi as ‘longitude’ Parameters rarray_like radius thetaarray_like inclination or polar angle phiarray_like azimuth angle Returns xarray x coordinate(s) in Cartesion space yarray y coordinate(s) in Cartesian space zarray z coordinate Notes See these pages: for excellent discussion of the many different conventions possible. Here we use the physics conventions, used in the wikipedia page. Derivations of the formulae are simple. Consider a vector x, y, z of length r (norm of x, y, z). The inclination angle (theta) can be found from: cos(theta) == z / r -> z == r * cos(theta). This gives the hypotenuse of the projection onto the XY plane, which we will call Q. Q == r*sin(theta). Now x / Q == cos(phi) -> x == r * sin(theta) * cos(phi) and so on. We have deliberately named this function sphere2cart rather than sph2cart to distinguish it from the Matlab function of that name, because the Matlab function uses an unusual convention for the angles that we did not want to replicate. The Matlab function is trivial to implement with the formulae given in the Matlab help. ### split_dki_param dipy.reconst.dki.split_dki_param(dki_params) Extract the diffusion tensor eigenvalues, the diffusion tensor eigenvector matrix, and the 15 independent elements of the kurtosis tensor from the model parameters estimated from the DKI model Parameters dki_paramsndarray (x, y, z, 27) or (n, 27) All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor Returns eigvalsarray (x, y, z, 3) or (n, 3) Eigenvalues from eigen decomposition of the tensor. eigvecsarray (x, y, z, 3, 3) or (n, 3, 3) Associated eigenvectors from eigen decomposition of the tensor. Eigenvectors are columnar (e.g. eigvecs[:,j] is associated with eigvals[j]) ktarray (x, y, z, 15) or (n, 15) Fifteen elements of the kurtosis tensor ### vec_val_vect dipy.reconst.dki.vec_val_vect() Vectorize vecs.diag(vals).vecs.T for last 2 dimensions of vecs Parameters vecsshape (…, M, N) array containing tensor in last two dimensions; M, N usually equal to (3, 3) valsshape (…, N) array diagonal values carried in last dimension, ... shape above must match that for vecs Returns resshape (…, M, M) array For all the dimensions ellided by ..., loops to get (M, N) vec matrix, and (N,) vals vector, and calculates vec.dot(np.diag(val).dot(vec.T). Raises ValueErrornon-matching ... dimensions of vecs, vals ValueErrornon-matching N dimensions of vecs, vals Examples Make a 3D array where the first dimension is only 1 >>> vecs = np.arange(9).reshape((1, 3, 3)) >>> vals = np.arange(3).reshape((1, 3)) >>> vec_val_vect(vecs, vals) array([[[ 9., 24., 39.], [ 24., 66., 108.], [ 39., 108., 177.]]])  That’s the same as the 2D case (apart from the float casting): >>> vecs = np.arange(9).reshape((3, 3)) >>> vals = np.arange(3) >>> np.dot(vecs, np.dot(np.diag(vals), vecs.T)) array([[ 9, 24, 39], [ 24, 66, 108], [ 39, 108, 177]])  ### wls_fit_dki dipy.reconst.dki.wls_fit_dki(design_matrix, data) Computes weighted linear least squares (WLS) fit to calculate the diffusion tensor and kurtosis tensor using a weighted linear regression diffusion kurtosis model 1. Parameters design_matrixarray (g, 22) Design matrix holding the covariants used to solve for the regression coefficients. dataarray (N, g) Data or response variables holding the data. Note that the last dimension should contain the data. It makes no copies of data. min_signaldefault = 1 All values below min_signal are repalced with min_signal. This is done in order to avoid taking log(0) durring the tensor fitting. Returns dki_paramsarray (N, 27) All parameters estimated from the diffusion kurtosis model for all N voxels. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor ### DiffusionKurtosisFit class dipy.reconst.dki_micro.DiffusionKurtosisFit(model, model_params) Class for fitting the Diffusion Kurtosis Model Attributes S0_hat directions For tracking - return the primary direction in each voxel evals Returns the eigenvalues of the tensor as an array evecs Returns the eigenvectors of the tensor as an array, columnwise kt Returns the 15 independent elements of the kurtosis tensor as an array quadratic_form Calculates the 3x3 diffusion tensor for each voxel shape Methods  ad() Axial diffusivity (AD) calculated from cached eigenvalues. adc(sphere) Calculate the apparent diffusion coefficient (ADC) in each direction on ak([min_kurtosis, max_kurtosis]) Axial Kurtosis (AK) of a diffusion kurtosis tensor [1]. akc(sphere) Calculates the apparent kurtosis coefficient (AKC) in each direction on the sphere for each voxel in the data color_fa() Color fractional anisotropy of diffusion tensor fa() Fractional anisotropy (FA) calculated from cached eigenvalues. ga() Geodesic anisotropy (GA) calculated from cached eigenvalues. kmax([sphere, gtol, mask]) Computes the maximum value of a single voxel kurtosis tensor linearity() Returns md() Mean diffusivity (MD) calculated from cached eigenvalues. mk([min_kurtosis, max_kurtosis]) Computes mean Kurtosis (MK) from the kurtosis tensor. mode() Tensor mode calculated from cached eigenvalues. odf(sphere) The diffusion orientation distribution function (dODF). planarity() Returns predict(gtab[, S0]) Given a DKI model fit, predict the signal on the vertices of a gradient table rd() Radial diffusivity (RD) calculated from cached eigenvalues. rk([min_kurtosis, max_kurtosis]) Radial Kurtosis (RK) of a diffusion kurtosis tensor [1]. sphericity() Returns trace() 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 Parameters modelDiffusionKurtosisModel Class instance Class instance containing the Diffusion Kurtosis Model for the fit model_paramsndarray (x, y, z, 27) or (n, 27) All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor ak(min_kurtosis=-0.42857142857142855, max_kurtosis=10) Axial Kurtosis (AK) of a diffusion kurtosis tensor [1]. Parameters min_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, axial kurtosis values that are smaller than min_kurtosis are replaced with -3./7 (theoretical kurtosis limit for regions that consist of water confined to spherical pores [2]) max_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, axial kurtosis values that are larger than max_kurtosis are replaced with max_kurtosis. Default = 10 Returns akarray Calculated AK. References 1(1,2,3) Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. Estimation of tensors and tensor-derived measures in diffusional kurtosis imaging. Magn Reson Med. 65(3), 823-836 2(1,2) Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: Robust estimation from DW-MRI using homogeneous polynomials. Proceedings of the 8th {IEEE} International Symposium on Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. doi: 10.1109/ISBI.2011.5872402 akc(sphere) Calculates the apparent kurtosis coefficient (AKC) in each direction on the sphere for each voxel in the data Parameters sphereSphere class instance Returns akcndarray The estimates of the apparent kurtosis coefficient in every direction on the input sphere Notes For each sphere direction with coordinates $$(n_{1}, n_{2}, n_{3})$$, the calculation of AKC is done using formula: $AKC(n)=\frac{MD^{2}}{ADC(n)^{2}}\sum_{i=1}^{3}\sum_{j=1}^{3} \sum_{k=1}^{3}\sum_{l=1}^{3}n_{i}n_{j}n_{k}n_{l}W_{ijkl}$ where $$W_{ijkl}$$ are the elements of the kurtosis tensor, MD the mean diffusivity and ADC the apparent diffusion coefficent computed as: $ADC(n)=\sum_{i=1}^{3}\sum_{j=1}^{3}n_{i}n_{j}D_{ij}$ where $$D_{ij}$$ are the elements of the diffusion tensor. kmax(sphere='repulsion100', gtol=1e-05, mask=None) Computes the maximum value of a single voxel kurtosis tensor Parameters sphereSphere class instance, optional The sphere providing sample directions for the initial search of the maximum value of kurtosis. gtolfloat, optional This input is to refine kurtosis maximum under the precision of the directions sampled on the sphere class instance. The gradient of the convergence procedure must be less than gtol before successful termination. If gtol is None, fiber direction is directly taken from the initial sampled directions of the given sphere object Returns max_valuefloat kurtosis tensor maximum value property kt Returns the 15 independent elements of the kurtosis tensor as an array mk(min_kurtosis=-0.42857142857142855, max_kurtosis=10) Computes mean Kurtosis (MK) from the kurtosis tensor. Parameters min_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, mean kurtosis values that are smaller than min_kurtosis are replaced with min_kurtosis. Default = -3./7 (theoretical kurtosis limit for regions that consist of water confined to spherical pores [2]) max_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, mean kurtosis values that are larger than max_kurtosis are replaced with max_kurtosis. Default = 10 Returns mkarray Calculated MK. Notes The MK analytical solution is calculated using the following equation [1]: $\begin{split}MK=F_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{1111}+ F_1(\lambda_2,\lambda_1,\lambda_3)\hat{W}_{2222}+ F_1(\lambda_3,\lambda_2,\lambda_1)\hat{W}_{3333}+ \\ F_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233}+ F_2(\lambda_2,\lambda_1,\lambda_3)\hat{W}_{1133}+ F_2(\lambda_3,\lambda_2,\lambda_1)\hat{W}_{1122}\end{split}$ where $$\hat{W}_{ijkl}$$ are the components of the $$W$$ tensor in the coordinates system defined by the eigenvectors of the diffusion tensor $$\mathbf{D}$$ and $\begin{split}F_1(\lambda_1,\lambda_2,\lambda_3)= \frac{(\lambda_1+\lambda_2+\lambda_3)^2} {18(\lambda_1-\lambda_2)(\lambda_1-\lambda_3)} [\frac{\sqrt{\lambda_2\lambda_3}}{\lambda_1} R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\ \frac{3\lambda_1^2-\lambda_1\lambda_2-\lambda_2\lambda_3- \lambda_1\lambda_3} {3\lambda_1 \sqrt{\lambda_2 \lambda_3}} R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-1 ]\end{split}$ and $\begin{split}F_2(\lambda_1,\lambda_2,\lambda_3)= \frac{(\lambda_1+\lambda_2+\lambda_3)^2} {3(\lambda_2-\lambda_3)^2} [\frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}} R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\ \frac{2\lambda_1-\lambda_2-\lambda_3}{3\sqrt{\lambda_2 \lambda_3}} R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-2]\end{split}$ where $$R_f$$ and $$R_d$$ are the Carlson’s elliptic integrals. References 1(1,2) Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. Estimation of tensors and tensor-derived measures in diffusional kurtosis imaging. Magn Reson Med. 65(3), 823-836 2(1,2) Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: Robust estimation from DW-MRI using homogeneous polynomials. Proceedings of the 8th {IEEE} International Symposium on Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. doi: 10.1109/ISBI.2011.5872402 predict(gtab, S0=1.0) Given a DKI model fit, predict the signal on the vertices of a gradient table Parameters gtaba GradientTable class instance The gradient table for this prediction S0float or ndarray (optional) The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1 Notes The predicted signal is given by: $S(n,b)=S_{0}e^{-bD(n)+\frac{1}{6}b^{2}D(n)^{2}K(n)}$ $$\mathbf{D(n)}$$ and $$\mathbf{K(n)}$$ can be computed from the DT and KT using the following equations: $D(n)=\sum_{i=1}^{3}\sum_{j=1}^{3}n_{i}n_{j}D_{ij}$ and $K(n)=\frac{MD^{2}}{D(n)^{2}}\sum_{i=1}^{3}\sum_{j=1}^{3} \sum_{k=1}^{3}\sum_{l=1}^{3}n_{i}n_{j}n_{k}n_{l}W_{ijkl}$ where $$D_{ij}$$ and $$W_{ijkl}$$ are the elements of the second-order DT and the fourth-order KT tensors, respectively, and $$MD$$ is the mean diffusivity. rk(min_kurtosis=-0.42857142857142855, max_kurtosis=10) Radial Kurtosis (RK) of a diffusion kurtosis tensor [1]. Parameters min_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, axial kurtosis values that are smaller than min_kurtosis are replaced with -3./7 (theoretical kurtosis limit for regions that consist of water confined to spherical pores [2]) max_kurtosisfloat (optional) To keep kurtosis values within a plausible biophysical range, axial kurtosis values that are larger than max_kurtosis are replaced with max_kurtosis. Default = 10 Returns rkarray Calculated RK. Notes RK is calculated with the following equation: $K_{\bot} = G_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2222} + G_1(\lambda_1,\lambda_3,\lambda_2)\hat{W}_{3333} + G_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233}$ where: $G_1(\lambda_1,\lambda_2,\lambda_3)= \frac{(\lambda_1+\lambda_2+\lambda_3)^2}{18\lambda_2(\lambda_2- \lambda_3)} \left (2\lambda_2 + \frac{\lambda_3^2-3\lambda_2\lambda_3}{\sqrt{\lambda_2\lambda_3}} \right)$ and $G_2(\lambda_1,\lambda_2,\lambda_3)= \frac{(\lambda_1+\lambda_2+\lambda_3)^2}{(\lambda_2-\lambda_3)^2} \left ( \frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}}-2 \right )$ References 1(1,2,3) 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 2(1,2) Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: Robust estimation from DW-MRI using homogeneous polynomials. Proceedings of the 8th {IEEE} International Symposium on Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. doi: 10.1109/ISBI.2011.5872402 ### DiffusionKurtosisModel class dipy.reconst.dki_micro.DiffusionKurtosisModel(gtab, fit_method='WLS', *args, **kwargs) Class for the Diffusion Kurtosis Model Methods  fit(data[, mask]) Fit method of the DKI model class predict(dki_params[, S0]) Predict a signal for this DKI model class instance given parameters. __init__(gtab, fit_method='WLS', *args, **kwargs) Diffusion Kurtosis Tensor Model [1] Parameters gtabGradientTable class instance fit_methodstr or callable str can be one of the following: ‘OLS’ or ‘ULLS’ for ordinary least squares dki.ols_fit_dki ‘WLS’ or ‘UWLLS’ for weighted ordinary least squares dki.wls_fit_dki callable has to have the signature: fit_method(design_matrix, data, *args, **kwargs) args, kwargsarguments and key-word arguments passed to the fit_method. See dki.ols_fit_dki, dki.wls_fit_dki for details References 1 Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. Estimation of tensors and tensor-derived measures in diffusional kurtosis imaging. Magn Reson Med. 65(3), 823-836 fit(data, mask=None) Fit method of the DKI model class Parameters dataarray The measured signal from one voxel. maskarray A boolean array used to mark the coordinates in the data that should be analyzed that has the shape data.shape[-1] predict(dki_params, S0=1.0) Predict a signal for this DKI model class instance given parameters. Parameters dki_paramsndarray (x, y, z, 27) or (n, 27) All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor S0float or ndarray (optional) The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1 ### KurtosisMicrostructuralFit class dipy.reconst.dki_micro.KurtosisMicrostructuralFit(model, model_params) Class for fitting the Diffusion Kurtosis Microstructural Model Attributes S0_hat awf Returns the volume fraction of the restricted diffusion compartment also known as axonal water fraction. axonal_diffusivity Returns the axonal diffusivity defined as the restricted diffusion tensor trace 1. directions For tracking - return the primary direction in each voxel evals Returns the eigenvalues of the tensor as an array evecs Returns the eigenvectors of the tensor as an array, columnwise hindered_ad Returns the axial diffusivity of the hindered compartment. hindered_evals Returns the eigenvalues of the hindered diffusion compartment. hindered_rd Returns the radial diffusivity of the hindered compartment. kt Returns the 15 independent elements of the kurtosis tensor as an array quadratic_form Calculates the 3x3 diffusion tensor for each voxel restricted_evals Returns the eigenvalues of the restricted diffusion compartment. shape tortuosity Returns the tortuosity of the hindered diffusion which is defined by ADe / RDe, where ADe and RDe are the axial and radial diffusivities of the hindered compartment 1. Methods  ad() Axial diffusivity (AD) calculated from cached eigenvalues. adc(sphere) Calculate the apparent diffusion coefficient (ADC) in each direction on ak([min_kurtosis, max_kurtosis]) Axial Kurtosis (AK) of a diffusion kurtosis tensor [R0b1a747e81c9-1]. akc(sphere) Calculates the apparent kurtosis coefficient (AKC) in each direction on the sphere for each voxel in the data color_fa() Color fractional anisotropy of diffusion tensor fa() Fractional anisotropy (FA) calculated from cached eigenvalues. ga() Geodesic anisotropy (GA) calculated from cached eigenvalues. kmax([sphere, gtol, mask]) Computes the maximum value of a single voxel kurtosis tensor linearity() Returns md() Mean diffusivity (MD) calculated from cached eigenvalues. mk([min_kurtosis, max_kurtosis]) Computes mean Kurtosis (MK) from the kurtosis tensor. mode() Tensor mode calculated from cached eigenvalues. odf(sphere) The diffusion orientation distribution function (dODF). planarity() Returns predict(gtab[, S0]) Given a DKI microstructural model fit, predict the signal on the vertices of a gradient table rd() Radial diffusivity (RD) calculated from cached eigenvalues. rk([min_kurtosis, max_kurtosis]) Radial Kurtosis (RK) of a diffusion kurtosis tensor [Rc4101656d30e-1]. sphericity() Returns trace() Trace of the tensor calculated from cached eigenvalues.  lower_triangular __init__(model, model_params) Initialize a KurtosisMicrostructural Fit class instance. Parameters modelDiffusionKurtosisModel Class instance Class instance containing the Diffusion Kurtosis Model for the fit model_paramsndarray (x, y, z, 40) or (n, 40) All parameters estimated from the diffusion kurtosis microstructural model. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor 4. Six elements of the hindered diffusion tensor 5. Six elements of the restricted diffusion tensor 6. 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 1(1,2) 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 property awf Returns 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 1(1,2) 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 property axonal_diffusivity Returns the axonal diffusivity defined as the restricted diffusion tensor trace [1]. References 1(1,2) 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 property hindered_ad Returns 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 1(1,2) 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 property hindered_evals Returns 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 1(1,2) 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 property hindered_rd Returns 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 1(1,2) 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 gtaba GradientTable class instance The gradient table for this prediction S0float or ndarray (optional) The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1 Notes The predicted signal is given by: $$S(\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). property restricted_evals Returns 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 1(1,2) 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 property tortuosity Returns the tortuosity of the hindered diffusion which is defined by ADe / RDe, where ADe and RDe are the axial and radial diffusivities of the hindered compartment [1]. Notes The hindered diffusion tensor can be seem as the tissue’s extra-cellular diffusion compartment [1]. References 1(1,2,3) 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 ### KurtosisMicrostructureModel class dipy.reconst.dki_micro.KurtosisMicrostructureModel(gtab, fit_method='WLS', *args, **kwargs) Class for the Diffusion Kurtosis Microstructural Model Methods  fit(data[, mask, sphere, gtol, awf_only]) Fit method of the Diffusion Kurtosis Microstructural Model predict(params[, S0]) 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]. Parameters gtabGradientTable class instance fit_methodstr or callable 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 callable has to have the signature: fit_method(design_matrix, data, *args, **kwargs) args, kwargsarguments and key-word arguments passed to the fit_method. See dki.ols_fit_dki, dki.wls_fit_dki for details References 1(1,2) 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 Parameters dataarray An 4D matrix containing the diffusion-weighted data. maskarray A boolean array used to mark the coordinates in the data that should be analyzed that has the shape data.shape[-1] sphereSphere class instance, optional The sphere providing sample directions for the initial search of the maximal value of kurtosis. gtolfloat, optional This input is to refine kurtosis maxima under the precision of the directions sampled on the sphere class instance. The gradient of the convergence procedure must be less than gtol before successful termination. If gtol is None, fiber direction is directly taken from the initial sampled directions of the given sphere object awf_onlybool, optiomal If set to true only the axonal volume fraction is computed from the kurtosis tensor. Default = False predict(params, S0=1.0) Predict a signal for the DKI microstructural model class instance given parameters. Parameters paramsndarray (x, y, z, 40) or (n, 40) All parameters estimated from the diffusion kurtosis microstructural model. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor 4. Six elements of the hindered diffusion tensor 5. Six elements of the restricted diffusion tensor 6. Axonal water fraction S0float or ndarray (optional) 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 1(1,2) 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 ### axial_diffusivity dipy.reconst.dki_micro.axial_diffusivity(evals, axis=-1) Axial Diffusivity (AD) of a diffusion tensor. Also called parallel diffusivity. Parameters evalsarray-like Eigenvalues of a diffusion tensor, must be sorted in descending order along axis. axisint Axis of evals which contains 3 eigenvalues. Returns adarray Calculated AD. Notes AD is calculated with the following equation: $AD = \lambda_1$ ### axonal_water_fraction dipy.reconst.dki_micro.axonal_water_fraction(dki_params, sphere='repulsion100', gtol=0.01, mask=None) Computes the axonal water fraction from DKI [1]. Parameters dki_paramsndarray (x, y, z, 27) or (n, 27) All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor sphereSphere class instance, optional The sphere providing sample directions for the initial search of the maximal value of kurtosis. gtolfloat, optional This input is to refine kurtosis maxima under the precision of the directions sampled on the sphere class instance. The gradient of the convergence procedure must be less than gtol before successful termination. If gtol is None, fiber direction is directly taken from the initial sampled directions of the given sphere object maskndarray A boolean array used to mark the coordinates in the data that should be analyzed that has the shape dki_params.shape[:-1] Returns awfndarray (x, y, z) or (n) Axonal Water Fraction References 1(1,2,3) Fieremans E, Jensen JH, Helpern JA, 2011. White matter characterization with diffusional kurtosis imaging. Neuroimage 58(1):177-88. doi: 10.1016/j.neuroimage.2011.06.006 ### decompose_tensor 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). Parameters tensorarray (…, 3, 3) Hermitian matrix representing a diffusion tensor. min_diffusivityfloat Because negative eigenvalues are not physical and small eigenvalues, much smaller than the diffusion weighting, cause quite a lot of noise in metrics such as fa, diffusivity values smaller than min_diffusivity are replaced with min_diffusivity. Returns eigvalsarray (…, 3) Eigenvalues from eigen decomposition of the tensor. Negative eigenvalues are replaced by zero. Sorted from largest to smallest. eigvecsarray (…, 3, 3) Associated eigenvectors from eigen decomposition of the tensor. Eigenvectors are columnar (e.g. eigvecs[…, :, j] is associated with eigvals[…, j]) ### diffusion_components 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]. Parameters dki_paramsndarray (x, y, z, 27) or (n, 27) All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor sphereSphere class instance, optional The sphere providing sample directions to sample the restricted and hindered cellular diffusion tensors. For more details see Fieremans et al., 2011. awfndarray (optional) Array containing values of the axonal water fraction that has the shape dki_params.shape[:-1]. If not given this will be automatically computed using axonal_water_fraction()” with function’s default precision. maskndarray (optional) A boolean array used to mark the coordinates in the data that should be analyzed that has the shape dki_params.shape[:-1] Returns edtndarray (x, y, z, 6) or (n, 6) Parameters of the hindered diffusion tensor. idtndarray (x, y, z, 6) or (n, 6) Parameters of the restricted diffusion tensor. Notes In the original article of DKI microstructural model [1], the hindered and restricted tensors were definde as the intra-cellular and extra-cellular diffusion compartments respectively. References 1(1,2,3,4) Fieremans E, Jensen JH, Helpern JA, 2011. White matter characterization with diffusional kurtosis imaging. Neuroimage 58(1):177-88. doi: 10.1016/j.neuroimage.2011.06.006 ### directional_diffusion 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]. Parameters dtarray (6,) elements of the diffusion tensor of the voxel. Varray (g, 3) g directions of a Sphere in Cartesian coordinates min_diffusivityfloat (optional) Because negative eigenvalues are not physical and small eigenvalues cause quite a lot of noise in diffusion-based metrics, diffusivity values smaller than min_diffusivity are replaced with min_diffusivity. Default = 0 Returns adcndarray (g,) Apparent diffusion coefficient (adc) in all g directions of a sphere for a single voxel. References 1(1,2,3) Neto Henriques R, Correia MM, Nunes RG, Ferreira HA (2015). Exploring the 3D geometry of the diffusion kurtosis tensor - Impact on the development of robust tractography procedures and novel biomarkers, NeuroImage 111: 85-99 ### directional_kurtosis 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]. Parameters dtarray (6,) elements of the diffusion tensor of the voxel. mdfloat mean diffusivity of the voxel ktarray (15,) elements of the kurtosis tensor of the voxel. Varray (g, 3) g directions of a Sphere in Cartesian coordinates min_diffusivityfloat (optional) Because negative eigenvalues are not physical and small eigenvalues cause quite a lot of noise in diffusion-based metrics, diffusivity values smaller than min_diffusivity are replaced with min_diffusivity. Default = 0 min_kurtosisfloat (optional) Because high-amplitude negative values of kurtosis are not physicaly and biologicaly pluasible, and these cause artefacts in kurtosis-based measures, directional kurtosis values smaller than min_kurtosis are replaced with min_kurtosis. Default = -3./7 (theoretical kurtosis limit for regions that consist of water confined to spherical pores [2]) adcndarray(g,) (optional) Apparent diffusion coefficient (adc) in all g directions of a sphere for a single voxel. advndarray(g,) (optional) Apparent diffusion variance (advc) in all g directions of a sphere for a single voxel. Returns akcndarray (g,) Apparent kurtosis coefficient (AKC) in all g directions of a sphere for a single voxel. References 1(1,2,3) Neto Henriques R, Correia MM, Nunes RG, Ferreira HA (2015). Exploring the 3D geometry of the diffusion kurtosis tensor - Impact on the development of robust tractography procedures and novel biomarkers, NeuroImage 111: 85-99 2(1,2) Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: Robust estimation from DW-MRI using homogeneous polynomials. Proceedings of the 8th {IEEE} International Symposium on Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. doi: 10.1109/ISBI.2011.5872402 ### dkimicro_prediction dipy.reconst.dki_micro.dkimicro_prediction(params, gtab, S0=1) Signal prediction given the DKI microstructure model parameters. Parameters paramsndarray (x, y, z, 40) or (n, 40) All parameters estimated from the diffusion kurtosis microstructure model. Parameters are ordered as follows: 1. Three diffusion tensor’s eigenvalues 2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector 3. Fifteen elements of the kurtosis tensor 4. Six elements of the hindered diffusion tensor 5. Six elements of the restricted diffusion tensor 6. Axonal water fraction gtaba GradientTable class instance The gradient table for this prediction S0float or ndarray The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1 Returns S(…, N) ndarray 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.

### dti_design_matrix

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)

Parameters
gtabA GradientTable class instance
dtypestring

Parameter to control the dtype of returned designed matrix

Returns
design_matrixarray (g,7)

Design matrix or B matrix assuming Gaussian distributed tensor model design_matrix[j, :] = (Bxx, Byy, Bzz, Bxy, Bxz, Byz, dummy)

### from_lower_triangular

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.

Parameters
Darray_like, (…, >6)

Unique elements of the tensors

Returns
tensorndarray (…, 3, 3)

3 by 3 tensors

### get_sphere

dipy.reconst.dki_micro.get_sphere(name='symmetric362')

provide triangulated spheres

Parameters
namestr

which sphere - one of: * ‘symmetric362’ * ‘symmetric642’ * ‘symmetric724’ * ‘repulsion724’ * ‘repulsion100’ * ‘repulsion200’

Returns
spherea dipy.core.sphere.Sphere class instance

Examples

>>> import numpy as np
>>> from dipy.data import get_sphere
>>> sphere = get_sphere('symmetric362')
>>> verts, faces = sphere.vertices, sphere.faces
>>> verts.shape == (362, 3)
True
>>> faces.shape == (720, 3)
True
>>> verts, faces = get_sphere('not a sphere name') #doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
...
DataError: No sphere called "not a sphere name"


### kurtosis_maximum

dipy.reconst.dki_micro.kurtosis_maximum(dki_params, sphere='repulsion100', gtol=0.01, mask=None)

Computes kurtosis maximum value

Parameters
dki_paramsndarray (x, y, z, 27) or (n, 27)

All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows:

1. Three diffusion tensor’s eingenvalues

2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector

3. Fifteen elements of the kurtosis tensor

sphereSphere class instance, optional

The sphere providing sample directions for the initial search of the maximal value of kurtosis.

gtolfloat, optional

This input is to refine kurtosis maximum under the precision of the directions sampled on the sphere class instance. The gradient of the convergence procedure must be less than gtol before successful termination. If gtol is None, fiber direction is directly taken from the initial sampled directions of the given sphere object

A boolean array used to mark the coordinates in the data that should be analyzed that has the shape dki_params.shape[:-1]

Returns
max_valuefloat

kurtosis tensor maximum value

max_dirarray (3,)

Cartesian coordinates of the direction of the maximal kurtosis value

### lower_triangular

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

Parameters
tensorarray_like (…, 3, 3)

a collection of 3, 3 diffusion tensors

b0float

if b0 is not none log(b0) is returned as the dummy variable

Returns
Dndarray

If b0 is none, then the shape will be (…, 6) otherwise (…, 7)

### mean_diffusivity

dipy.reconst.dki_micro.mean_diffusivity(evals, axis=-1)

Mean Diffusivity (MD) of a diffusion tensor.

Parameters
evalsarray-like

Eigenvalues of a diffusion tensor.

axisint

Axis of evals which contains 3 eigenvalues.

Returns
mdarray

Calculated MD.

Notes

MD is calculated with the following equation:

$MD = \frac{\lambda_1 + \lambda_2 + \lambda_3}{3}$

### ndindex

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.

Parameters
shapetuple of ints

The dimensions of the array.

Examples

>>> from dipy.core.ndindex import ndindex
>>> shape = (3, 2, 1)
>>> for index in ndindex(shape):
...     print(index)
(0, 0, 0)
(0, 1, 0)
(1, 0, 0)
(1, 1, 0)
(2, 0, 0)
(2, 1, 0)


dipy.reconst.dki_micro.radial_diffusivity(evals, axis=-1)

Radial Diffusivity (RD) of a diffusion tensor. Also called perpendicular diffusivity.

Parameters
evalsarray-like

Eigenvalues of a diffusion tensor, must be sorted in descending order along axis.

axisint

Axis of evals which contains 3 eigenvalues.

Returns
rdarray

Calculated RD.

Notes

RD is calculated with the following equation:

$RD = \frac{\lambda_2 + \lambda_3}{2}$

### split_dki_param

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

Parameters
dki_paramsndarray (x, y, z, 27) or (n, 27)

All parameters estimated from the diffusion kurtosis model. Parameters are ordered as follows:

1. Three diffusion tensor’s eigenvalues

2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector

3. Fifteen elements of the kurtosis tensor

Returns
eigvalsarray (x, y, z, 3) or (n, 3)

Eigenvalues from eigen decomposition of the tensor.

eigvecsarray (x, y, z, 3, 3) or (n, 3, 3)

Associated eigenvectors from eigen decomposition of the tensor. Eigenvectors are columnar (e.g. eigvecs[:,j] is associated with eigvals[j])

ktarray (x, y, z, 15) or (n, 15)

Fifteen elements of the kurtosis tensor

### tortuosity

dipy.reconst.dki_micro.tortuosity(hindered_ad, hindered_rd)

Computes the tortuosity of the hindered diffusion compartment given its axial and radial diffusivities

Parameters

Array containing the values of the hindered axial diffusivity.

hindered_rd: ndarray

Array containing the values of the hindered radial diffusivity.

Returns
Tortuosity of the hindered diffusion compartment

### trace

dipy.reconst.dki_micro.trace(evals, axis=-1)

Trace of a diffusion tensor.

Parameters
evalsarray-like

Eigenvalues of a diffusion tensor.

axisint

Axis of evals which contains 3 eigenvalues.

Returns
tracearray

Calculated trace of the diffusion tensor.

Notes

Trace is calculated with the following equation:

$Trace = \lambda_1 + \lambda_2 + \lambda_3$

### vec_val_vect

dipy.reconst.dki_micro.vec_val_vect()

Vectorize vecs.diag(vals).vecs.T for last 2 dimensions of vecs

Parameters
vecsshape (…, M, N) array

containing tensor in last two dimensions; M, N usually equal to (3, 3)

valsshape (…, N) array

diagonal values carried in last dimension, ... shape above must match that for vecs

Returns
resshape (…, M, M) array

For all the dimensions ellided by ..., loops to get (M, N) vec matrix, and (N,) vals vector, and calculates vec.dot(np.diag(val).dot(vec.T).

Raises
ValueErrornon-matching ... dimensions of vecs, vals
ValueErrornon-matching N dimensions of vecs, vals

Examples

Make a 3D array where the first dimension is only 1

>>> vecs = np.arange(9).reshape((1, 3, 3))
>>> vals = np.arange(3).reshape((1, 3))
>>> vec_val_vect(vecs, vals)
array([[[   9.,   24.,   39.],
[  24.,   66.,  108.],
[  39.,  108.,  177.]]])


That’s the same as the 2D case (apart from the float casting):

>>> vecs = np.arange(9).reshape((3, 3))
>>> vals = np.arange(3)
>>> np.dot(vecs, np.dot(np.diag(vals), vecs.T))
array([[  9,  24,  39],
[ 24,  66, 108],
[ 39, 108, 177]])


### Cache

class dipy.reconst.dsi.Cache

Bases: object

Cache values based on a key object (such as a sphere or gradient table).

Notes

This class is meant to be used as a mix-in:

class MyModel(Model, Cache):
pass

class MyModelFit(Fit):
pass


Inside a method on the fit, typical usage would be:

def odf(sphere):
M = self.model.cache_get('odf_basis_matrix', key=sphere)

if M is None:
M = self._compute_basis_matrix(sphere)
self.model.cache_set('odf_basis_matrix', key=sphere, value=M)


Methods

 Clear the cache. cache_get(tag, key[, default]) Retrieve a value from the cache. cache_set(tag, key, value) Store a value in the cache.
__init__(self, /, *args, **kwargs) Initialize self. See help(type(self)) for accurate signature. cache_clear() Clear the cache. cache_get(tag, key, default=None) Retrieve a value from the cache. Parameters tagstr Description of the cached value. keyobject Key object used to look up the cached value. defaultobject Value to be returned if no cached entry is found. Returns vobject Value from the cache associated with (tag, key). Returns default if no cached entry is found. cache_set(tag, key, value) Store a value in the cache. Parameters tagstr Description of the cached value. keyobject Key object used to look up the cached value. valueobject Value stored in the cache for each unique combination of (tag, key). Examples >>> def compute_expensive_matrix(parameters): ... # Imagine the following computation is very expensive ... return (p**2 for p in parameters)  >>> 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  ### DiffusionSpectrumDeconvFit class dipy.reconst.dsi.DiffusionSpectrumDeconvFit(model, data) Methods  msd_discrete([normalized]) Calculates the mean squared displacement on the discrete propagator odf(sphere) 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 rtop_pdf([normalized]) Calculates the return to origin probability from the propagator, which is the propagator evaluated at zero (see Descoteaux et Al. rtop_signal([filtering]) 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 Parameters modelobject, DiffusionSpectrumModel data1d ndarray, signal values pdf() 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 ### DiffusionSpectrumDeconvModel class dipy.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) Methods  cache_clear() Clear the cache. cache_get(tag, key[, default]) Retrieve a value from the cache. cache_set(tag, key, value) Store a value in the cache. fit(data[, mask]) Fit method for every voxel in data __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. ..math:: nowrap begin{eqnarray*} 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]. Parameters gtabGradientTable, Gradient directions and bvalues container class qgrid_sizeint, 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). r_startfloat, ODF is sampled radially in the PDF. This parameters shows where the sampling should start. r_endfloat, Radial endpoint of ODF sampling r_stepfloat, Step size of the ODf sampling from r_start to r_end filter_widthfloat, Strength of the hanning filter References 1(1,2) Canales-Rodriguez E.J et al., “Deconvolution in Diffusion Spectrum Imaging”, Neuroimage, 2010. 2 Biggs David S.C. et al., “Acceleration of Iterative Image Restoration Algorithms”, Applied Optics, vol. 36, No. 8, p. 1766-1775, 1997. fit(data, mask=None) Fit method for every voxel in data ### DiffusionSpectrumFit class dipy.reconst.dsi.DiffusionSpectrumFit(model, data) Methods  msd_discrete([normalized]) Calculates the mean squared displacement on the discrete propagator odf(sphere) Calculates the real discrete odf for a given discrete sphere pdf([normalized]) Applies the 3D FFT in the q-space grid to generate the diffusion propagator rtop_pdf([normalized]) Calculates the return to origin probability from the propagator, which is the propagator evaluated at zero (see Descoteaux et Al. rtop_signal([filtering]) 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 Parameters modelobject, DiffusionSpectrumModel data1d ndarray, signal values msd_discrete(normalized=True) Calculates the mean squared displacement on the discrete propagator ..math:: nowrap begin{equation} 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]). Parameters normalizedboolean, optional Whether to normalize the propagator by its sum in order to obtain a pdf. Default: True Returns msdfloat the mean square displacement References 1(1,2) Wu Y. et al., “Hybrid diffusion imaging”, NeuroImage, vol 36, 1. 617-629, 2007. odf(sphere) Calculates the real discrete odf for a given discrete sphere ..math:: nowrap begin{equation} 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. pdf(normalized=True) Applies the 3D FFT in the q-space grid to generate the diffusion propagator 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) Parameters normalizedboolean, optional Whether to normalize the propagator by its sum in order to obtain a pdf. Default: True. Returns rtopfloat the return to origin probability References 1(1,2) Descoteaux M. et al., “Multiple q-shell diffusion propagator imaging”, Medical Image Analysis, vol 15, No. 4, p. 603-621, 2011. 2(1,2) Tuch D.S., “Diffusion MRI of Complex Tissue Structure”, PhD Thesis, 2002. 3(1,2) 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 rtop_signal(filtering=True) Calculates the return to origin probability (rtop) from the signal rtop equals to the sum of all signal values Parameters filteringboolean, optional Whether to perform Hanning filtering. Default: True Returns rtopfloat the return to origin probability ### DiffusionSpectrumModel class dipy.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) Methods  cache_clear() Clear the cache. cache_get(tag, key[, default]) Retrieve a value from the cache. cache_set(tag, key, value) Store a value in the cache. fit(data[, mask]) Fit method for every voxel in data __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})$$ ..math:: nowrap begin{eqnarray} 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]. Parameters gtabGradientTable, Gradient directions and bvalues container class qgrid_sizeint, 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). r_startfloat, ODF is sampled radially in the PDF. This parameters shows where the sampling should start. r_endfloat, Radial endpoint of ODF sampling r_stepfloat, Step size of the ODf sampling from r_start to r_end filter_widthfloat, Strength of the hanning filter See also dipy.reconst.gqi.GeneralizedQSampling Notes 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. 1. We assume that data only one b0 volume is provided. References 1(1,2) Wedeen V.J et al., “Mapping Complex Tissue Architecture With Diffusion Spectrum Magnetic Resonance Imaging”, MRM 2005. 2 Canales-Rodriguez E.J et al., “Deconvolution in Diffusion Spectrum Imaging”, Neuroimage, 2010. 3(1,2) 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  fit(data, mask=None) Fit method for every voxel in data ### OdfFit class dipy.reconst.dsi.OdfFit(model, data) Methods  odf(sphere) To be implemented but specific odf models __init__(model, data) Initialize self. See help(type(self)) for accurate signature. odf(sphere) To be implemented but specific odf models ### OdfModel class dipy.reconst.dsi.OdfModel(gtab) 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  fit(data) To be implemented by specific odf models __init__(gtab) Initialization of the abstract class for signal reconstruction models Parameters gtabGradientTable class instance fit(data) To be implemented by specific odf models ### LR_deconv dipy.reconst.dsi.LR_deconv(prop, psf, numit=5, acc_factor=1) Perform Lucy-Richardson deconvolution algorithm on a 3D array. Parameters prop3-D ndarray of dtype float The 3D volume to be deconvolve psf3-D ndarray of dtype float The filter that will be used for the deconvolution. numitint Number of Lucy-Richardson iteration to perform. acc_factorfloat Exponential acceleration factor as in [1]. References 1(1,2) Biggs David S.C. et al., “Acceleration of Iterative Image Restoration Algorithms”, Applied Optics, vol. 36, No. 8, p. 1766-1775, 1997. ### create_qspace dipy.reconst.dsi.create_qspace(gtab, origin) create the 3D grid which holds the signal values (q-space) Parameters gtabGradientTable origin(3,) ndarray center of qspace Returns qgridndarray qspace coordinates ### create_qtable dipy.reconst.dsi.create_qtable(gtab, origin) create a normalized version of gradients Parameters gtabGradientTable origin(3,) ndarray center of qspace Returns qtablendarray ### fftn 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. Parameters xarray_like The (n-dimensional) array to transform. shapetuple of ints, optional 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 i-th dimension is padded with zeros. If shape[i] < x.shape[i], the i-th dimension is truncated to length shape[i]. axesarray_like of ints, optional The axes of x (y if shape is not None) along which the transform is applied. overwrite_xbool, optional If True, the contents of x can be destroyed. Default is False. Returns ycomplex-valued n-dimensional numpy array The (n-dimensional) DFT of the input array. See also ifftn Notes 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  ### fftshift 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. Parameters xarray_like Input array. axesint or shape tuple, optional Axes over which to shift. Default is None, which shifts all axes. Returns yndarray The shifted array. See also ifftshift The inverse of fftshift. Examples >>> freqs = np.fft.fftfreq(10, 0.1) >>> freqs array([ 0., 1., 2., 3., 4., -5., -4., -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.]])  ### gen_PSF dipy.reconst.dsi.gen_PSF(qgrid_sampling, siz_x, siz_y, siz_z) Generate a PSF for DSI Deconvolution by taking the ifft of the binary q-space sampling mask and truncating it to keep only the center. ### half_to_full_qspace 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. Parameters dataarray, shape (X, Y, Z, W) where (X, Y, Z) volume size and W number of gradient directions gtabGradientTable container for b-values and b-vectors (gradient directions) Returns new_dataarray, shape (X, Y, Z, 2 * W -1) new_gtabGradientTable 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. ### hanning_filter 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. Parameters gtabGradientTable filter_widthint origin(3,) ndarray center of qspace Returns filter(N,) ndarray where N is the number of non-b0 gradient directions ### ifftshift 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. Parameters xarray_like Input array. axesint or shape tuple, optional Axes over which to calculate. Defaults to None, which shifts all axes. Returns yndarray The shifted array. See also fftshift Shift 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.]])  ### map_coordinates 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. Parameters inputarray_like The input array. coordinatesarray_like The coordinates at which input is evaluated. outputarray or dtype, optional 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. orderint, optional The order of the spline interpolation, default is 3. The order has to be in the range 0-5. mode{‘reflect’, ‘constant’, ‘nearest’, ‘mirror’, ‘wrap’}, optional The mode parameter determines how the input array is extended when the filter overlaps a border. Default is ‘reflect’. Behavior for each valid value is as follows: ‘reflect’ (d c b a | a b c d | d c b a) The input is extended by reflecting about the edge of the last pixel. ‘constant’ (k k k k | a b c d | k k k k) The input is extended by filling all values beyond the edge with the same constant value, defined by the cval parameter. ‘nearest’ (a a a a | a b c d | d d d d) The input is extended by replicating the last pixel. ‘mirror’ (d c b | a b c d | c b a) The input is extended by reflecting about the center of the last pixel. ‘wrap’ (a b c d | a b c d | a b c d) The input is extended by wrapping around to the opposite edge. cvalscalar, optional Value to fill past edges of input if mode is ‘constant’. Default is 0.0. prefilterbool, optional 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. Returns map_coordinatesndarray 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.interpolate Examples >>> 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)  ### multi_voxel_fit dipy.reconst.dsi.multi_voxel_fit(single_voxel_fit) Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition ### pdf_interp_coords dipy.reconst.dsi.pdf_interp_coords(sphere, rradius, origin) Precompute coordinates for ODF calculation from the PDF Parameters sphereobject, Sphere rradiusarray, shape (N,) line interpolation points originarray, shape (3,) center of the grid ### pdf_odf dipy.reconst.dsi.pdf_odf(Pr, rradius, interp_coords) Calculates the real ODF from the diffusion propagator(PDF) Pr Parameters Prarray, shape (X, X, X) probability density function rradiusarray, shape (N,) interpolation range on the radius interp_coordsarray, shape (3, M, N) coordinates in the pdf for interpolating the odf ### project_hemisph_bvecs dipy.reconst.dsi.project_hemisph_bvecs(gtab) Project any near identical bvecs to the other hemisphere Parameters gtabobject, GradientTable Notes Useful only when working with some types of dsi data. ### threshold_propagator dipy.reconst.dsi.threshold_propagator(P, estimated_snr=15.0) Applies hard threshold on the propagator to remove background noise for the deconvolution. ### ReconstModel class dipy.reconst.dti.ReconstModel(gtab) Bases: object Abstract class for signal reconstruction models Methods  fit __init__(gtab) Initialization of the abstract class for signal reconstruction models Parameters gtabGradientTable class instance fit(data, mask=None, **kwargs) ### TensorFit class dipy.reconst.dti.TensorFit(model, model_params, model_S0=None) Bases: object Attributes S0_hat directions For tracking - return the primary direction in each voxel evals Returns the eigenvalues of the tensor as an array evecs Returns the eigenvectors of the tensor as an array, columnwise quadratic_form Calculates the 3x3 diffusion tensor for each voxel shape Methods  Axial diffusivity (AD) calculated from cached eigenvalues. adc(sphere) 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. Returns Mean diffusivity (MD) calculated from cached eigenvalues. Tensor mode calculated from cached eigenvalues. odf(sphere) The diffusion orientation distribution function (dODF). Returns predict(gtab[, S0, step]) Given a model fit, predict the signal on the vertices of a sphere Radial diffusivity (RD) calculated from cached eigenvalues. Returns Trace of the tensor calculated from cached eigenvalues.  lower_triangular __init__(model, model_params, model_S0=None) Initialize a TensorFit class instance. property S0_hat ad() Axial diffusivity (AD) calculated from cached eigenvalues. Returns adarray (V, 1) Calculated AD. Notes RD is calculated with the following equation: $AD = \lambda_1$ adc(sphere) Calculate the apparent diffusion coefficient (ADC) in each direction on the sphere for each voxel in the data Parameters sphereSphere class instance Returns adcndarray The estimates of the apparent diffusion coefficient in every direction on the input sphere ec{b} Q ec{b}^T Where Q is the quadratic form of the tensor. color_fa() Color fractional anisotropy of diffusion tensor property directions For tracking - return the primary direction in each voxel property evals Returns the eigenvalues of the tensor as an array property evecs Returns the eigenvectors of the tensor as an array, columnwise fa() Fractional anisotropy (FA) calculated from cached eigenvalues. ga() Geodesic anisotropy (GA) calculated from cached eigenvalues. linearity() Returns linearityarray Calculated linearity of the diffusion tensor 1. Notes Linearity is calculated with the following equation: $Linearity = \frac{\lambda_1-\lambda_2}{\lambda_1+\lambda_2+\lambda_3}$ References [1] Westin C.-F., Peled S., Gubjartsson H., Kikinis R., Jolesz F., “Geometrical diffusion measures for MRI from tensor basis analysis” in Proc. 5th Annual ISMRM, 1997. lower_triangular(b0=None) md() Mean diffusivity (MD) calculated from cached eigenvalues. Returns mdarray (V, 1) Calculated MD. Notes MD is calculated with the following equation: $MD = \frac{\lambda_1+\lambda_2+\lambda_3}{3}$ mode() Tensor mode calculated from cached eigenvalues. odf(sphere) The diffusion orientation distribution function (dODF). This is an estimate of the diffusion distance in each direction Parameters sphereSphere class instance. The dODF is calculated in the vertices of this input. Returns odfndarray The diffusion distance in every direction of the sphere in every voxel in the input data. Notes This is based on equation 3 in [Aganj2010]. To re-derive it from scratch, follow steps in [Descoteaux2008], Section 7.9 Equation 7.24 but with an $$r^2$$ term in the integral. References Aganj2010(1,2) 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 Descoteaux2008(1,2) 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() Returns sphericityarray Calculated sphericity of the diffusion tensor 1. Notes Sphericity is calculated with the following equation: $Sphericity = \frac{2 (\lambda_2 - \lambda_3)}{\lambda_1+\lambda_2+\lambda_3}$ References [1] Westin C.-F., Peled S., Gubjartsson H., Kikinis R., Jolesz F., “Geometrical diffusion measures for MRI from tensor basis analysis” in Proc. 5th Annual ISMRM, 1997. predict(gtab, S0=None, step=None) Given a model fit, predict the signal on the vertices of a sphere Parameters gtaba GradientTable class instance This encodes the directions for which a prediction is made S0float array The mean non-diffusion weighted signal in each voxel. Default: The fitted S0 value in all voxels if it was fitted. Otherwise 1 in all voxels. stepint The chunk size as a number of voxels. Optional parameter with default value 10,000. In order to increase speed of processing, tensor fitting is done simultaneously over many voxels. This parameter sets the number of voxels that will be fit at once in each iteration. A larger step value should speed things up, but it will also take up more memory. It is advisable to keep an eye on memory consumption as this value is increased. Notes The predicted signal is given by: $S( heta, b) = S_0 * e^{-b ADC}$ Where: .. math 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 property quadratic_form Calculates the 3x3 diffusion tensor for each voxel rd() Radial diffusivity (RD) calculated from cached eigenvalues. Returns rdarray (V, 1) Calculated RD. Notes RD is calculated with the following equation: $RD = \frac{\lambda_2 + \lambda_3}{2}$ property shape sphericity() Returns sphericityarray Calculated sphericity of the diffusion tensor 1. Notes Sphericity is calculated with the following equation: $Sphericity = \frac{3 \lambda_3}{\lambda_1+\lambda_2+\lambda_3}$ References [1] Westin C.-F., Peled S., Gubjartsson H., Kikinis R., Jolesz F., “Geometrical diffusion measures for MRI from tensor basis analysis” in Proc. 5th Annual ISMRM, 1997. trace() Trace of the tensor calculated from cached eigenvalues. Returns tracearray (V, 1) Calculated trace. Notes The trace is calculated with the following equation: $trace = \lambda_1 + \lambda_2 + \lambda_3$ ### TensorModel class dipy.reconst.dti.TensorModel(gtab, fit_method='WLS', return_S0_hat=False, *args, **kwargs) Diffusion Tensor Methods  fit(data[, mask]) Fit method of the DTI model class predict(dti_params[, S0]) Predict a signal for this TensorModel class instance given parameters. __init__(gtab, fit_method='WLS', return_S0_hat=False, *args, **kwargs) A Diffusion Tensor Model [1], [2]. Parameters gtabGradientTable class instance fit_methodstr or callable str can be one of the following: ‘WLS’ for weighted least squares dti.wls_fit_tensor() ‘LS’ or ‘OLS’ for ordinary least squares dti.ols_fit_tensor() ‘NLLS’ for non-linear least-squares dti.nlls_fit_tensor() ‘RT’ or ‘restore’ or ‘RESTORE’ for RESTORE robust tensor fitting [3] dti.restore_fit_tensor() callable has to have the signature: fit_method(design_matrix, data, *args, **kwargs) return_S0_hatbool Boolean to return (True) or not (False) the S0 values for the fit. args, kwargsarguments and key-word arguments passed to the fit_method. See dti.wls_fit_tensor, dti.ols_fit_tensor for details min_signalfloat The minimum signal value. Needs to be a strictly positive number. Default: minimal signal in the data provided to fit. Notes In order to increase speed of processing, tensor fitting is done simultaneously over many voxels. Many fit_methods use the ‘step’ parameter to set the number of voxels that will be fit at once in each iteration. This is the chunk size as a number of voxels. A larger step value should speed things up, but it will also take up more memory. It is advisable to keep an eye on memory consumption as this value is increased. E.g., in iter_fit_tensor() we have a default step value of 1e4 References 1(1,2) Basser, P.J., Mattiello, J., LeBihan, D., 1994. Estimation of the effective self-diffusion tensor from the NMR spin echo. J Magn Reson B 103, 247-254. 2(1,2) Basser, P., Pierpaoli, C., 1996. Microstructural and physiological features of tissues elucidated by quantitative diffusion-tensor MRI. Journal of Magnetic Resonance 111, 209-219. 3(1,2) Lin-Ching C., Jones D.K., Pierpaoli, C. 2005. RESTORE: Robust estimation of tensors by outlier rejection. MRM 53: 1088-1095 fit(data, mask=None) Fit method of the DTI model class Parameters dataarray The measured signal from one voxel. maskarray A boolean array used to mark the coordinates in the data that should be analyzed that has the shape data.shape[:-1] predict(dti_params, S0=1.0) Predict a signal for this TensorModel class instance given parameters. Parameters dti_paramsndarray The last dimension should have 12 tensor parameters: 3 eigenvalues, followed by the 3 eigenvectors S0float or ndarray The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1 ### apparent_diffusion_coef dipy.reconst.dti.apparent_diffusion_coef(q_form, sphere) Calculate the apparent diffusion coefficient (ADC) in each direction of a sphere. Parameters q_formndarray The quadratic form of a tensor, or an array with quadratic forms of tensors. Should be of shape (…, 3, 3) spherea Sphere class instance The ADC will be calculated for each of the vertices in the sphere ec{b} Q ec{b}^T Where Q is the quadratic form of the tensor. ### auto_attr dipy.reconst.dti.auto_attr(func) Decorator to create OneTimeProperty attributes. Parameters funcmethod The method that will be called the first time to compute a value. Afterwards, the method’s name will be a standard attribute holding the value of this computation. Examples >>> class MagicProp(object): ... @auto_attr ... def a(self): ... return 99 ... >>> x = MagicProp() >>> 'a' in x.__dict__ False >>> x.a 99 >>> 'a' in x.__dict__ True  ### axial_diffusivity dipy.reconst.dti.axial_diffusivity(evals, axis=-1) Axial Diffusivity (AD) of a diffusion tensor. Also called parallel diffusivity. Parameters evalsarray-like Eigenvalues of a diffusion tensor, must be sorted in descending order along axis. axisint Axis of evals which contains 3 eigenvalues. Returns adarray Calculated AD. Notes AD is calculated with the following equation: $AD = \lambda_1$ ### color_fa dipy.reconst.dti.color_fa(fa, evecs) Color fractional anisotropy of diffusion tensor Parameters faarray-like Array of the fractional anisotropy (can be 1D, 2D or 3D) evecsarray-like eigen vectors from the tensor model Returns rgbArray with 3 channels for each color as the last dimension. Colormap of the FA with red for the x value, y for the green value and z for the blue value. ec{e})) imes fa ### decompose_tensor 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). Parameters tensorarray (…, 3, 3) Hermitian matrix representing a diffusion tensor. min_diffusivityfloat Because negative eigenvalues are not physical and small eigenvalues, much smaller than the diffusion weighting, cause quite a lot of noise in metrics such as fa, diffusivity values smaller than min_diffusivity are replaced with min_diffusivity. Returns eigvalsarray (…, 3) Eigenvalues from eigen decomposition of the tensor. Negative eigenvalues are replaced by zero. Sorted from largest to smallest. eigvecsarray (…, 3, 3) Associated eigenvectors from eigen decomposition of the tensor. Eigenvectors are columnar (e.g. eigvecs[…, :, j] is associated with eigvals[…, j]) ### design_matrix dipy.reconst.dti.design_matrix(gtab, dtype=None) Constructs design matrix for DTI weighted least squares or least squares fitting. (Basser et al., 1994a) Parameters gtabA GradientTable class instance dtypestring Parameter to control the dtype of returned designed matrix Returns design_matrixarray (g,7) Design matrix or B matrix assuming Gaussian distributed tensor model design_matrix[j, :] = (Bxx, Byy, Bzz, Bxy, Bxz, Byz, dummy) ### determinant dipy.reconst.dti.determinant(q_form) The determinant of a tensor, given in quadratic form Parameters q_formndarray The quadratic form of a tensor, or an array with quadratic forms of tensors. Should be of shape (x, y, z, 3, 3) or (n, 3, 3) or (3, 3). Returns detarray The determinant of the tensor in each spatial coordinate ### deviatoric dipy.reconst.dti.deviatoric(q_form) Calculate the deviatoric (anisotropic) part of the tensor [1]. Parameters q_formndarray The quadratic form of a tensor, or an array with quadratic forms of tensors. Should be of shape (x,y,z,3,3) or (n, 3, 3) or (3,3). Returns A_squigglendarray 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]): $\widetilde{A} = A - \bar{A}$ Where $$A$$ is the tensor quadratic form and $$\bar{A}$$ is the anisotropic part of the tensor. References 1(1,2,3,4) Daniel B. Ennis and G. Kindlmann, “Orthogonal Tensor Invariants and the Analysis of Diffusion Tensor Magnetic Resonance Images”, Magnetic Resonance in Medicine, vol. 55, no. 1, pp. 136-146, 2006. ### eig_from_lo_tri 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. Parameters dataarray_like (…, 6) diffusion tensors elements stored in lower triangular order min_diffusivityfloat See decompose_tensor() Returns dti_paramsarray (…, 12) Eigen-values and eigen-vectors of the same array. ### fractional_anisotropy dipy.reconst.dti.fractional_anisotropy(evals, axis=-1) Fractional anisotropy (FA) of a diffusion tensor. Parameters evalsarray-like Eigenvalues of a diffusion tensor. axisint Axis of evals which contains 3 eigenvalues. Returns faarray Calculated FA. Range is 0 <= FA <= 1. Notes FA is calculated using the following equation: $FA = \sqrt{\frac{1}{2}\frac{(\lambda_1-\lambda_2)^2+(\lambda_1- \lambda_3)^2+(\lambda_2-\lambda_3)^2}{\lambda_1^2+ \lambda_2^2+\lambda_3^2}}$ ### from_lower_triangular 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. Parameters Darray_like, (…, >6) Unique elements of the tensors Returns tensorndarray (…, 3, 3) 3 by 3 tensors ### geodesic_anisotropy dipy.reconst.dti.geodesic_anisotropy(evals, axis=-1) Geodesic anisotropy (GA) of a diffusion tensor. Parameters evalsarray-like Eigenvalues of a diffusion tensor. axisint Axis of evals which contains 3 eigenvalues. Returns gaarray Calculated GA. In the range 0 to +infinity Notes GA is calculated using the following equation given in [1]: $GA = \sqrt{\sum_{i=1}^3 \log^2{\left ( \lambda_i/<\mathbf{D}> \right )}}, \quad \textrm{where} \quad <\mathbf{D}> = (\lambda_1\lambda_2\lambda_3)^{1/3}$ 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 1(1,2,3,4,5) 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. 2(1,2,3) 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. 3(1,2,3) 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. 4 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. ### get_sphere dipy.reconst.dti.get_sphere(name='symmetric362') provide triangulated spheres Parameters namestr which sphere - one of: * ‘symmetric362’ * ‘symmetric642’ * ‘symmetric724’ * ‘repulsion724’ * ‘repulsion100’ * ‘repulsion200’ Returns spherea dipy.core.sphere.Sphere class instance Examples >>> import numpy as np >>> from dipy.data import get_sphere >>> sphere = get_sphere('symmetric362') >>> verts, faces = sphere.vertices, sphere.faces >>> verts.shape == (362, 3) True >>> faces.shape == (720, 3) True >>> verts, faces = get_sphere('not a sphere name') #doctest: +IGNORE_EXCEPTION_DETAIL Traceback (most recent call last): ... DataError: No sphere called "not a sphere name"  ### gradient_table dipy.reconst.dti.gradient_table(bvals, bvecs=None, big_delta=None, small_delta=None, b0_threshold=50, atol=0.01) 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. Parameters bvalscan be any of the four options 1. an array of shape (N,) or (1, N) or (N, 1) with the b-values. 2. a path for the file which contains an array like the above (1). 3. 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. 4. a path for the file which contains an array like the one at (3). bvecscan be any of two options 1. an array of shape (N, 3) or (3, N) with the b-vectors. 2. a path for the file which contains an array like the previous. big_deltafloat acquisition pulse separation time in seconds (default None) small_deltafloat acquisition pulse duration time in seconds (default None) b0_thresholdfloat All b-values with values less than or equal to bo_threshold are considered as b0s i.e. without diffusion weighting. atolfloat All b-vectors need to be unit vectors up to a tolerance. Returns gradientsGradientTable A GradientTable with all the gradient information. Notes 1. 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__. 2. We assume that the minimum number of b-values is 7. 3. 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  ### isotropic dipy.reconst.dti.isotropic(q_form) Calculate the isotropic part of the tensor [Rd0568a744381-1]. Parameters q_formndarray The quadratic form of a tensor, or an array with quadratic forms of tensors. Should be of shape (x,y,z,3,3) or (n, 3, 3) or (3,3). Returns A_hat: ndarray The isotropic part of the tensor in each spatial coordinate rac{1}{2} tr(A) I ### iter_fit_tensor 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. Parameters stepint The chunk size as a number of voxels. Optional parameter with default value 10,000. In order to increase speed of processing, tensor fitting is done simultaneously over many voxels. This parameter sets the number of voxels that will be fit at once in each iteration. A larger step value should speed things up, but it will also take up more memory. It is advisable to keep an eye on memory consumption as this value is increased. ### linearity dipy.reconst.dti.linearity(evals, axis=-1) The linearity of the tensor 1 Parameters evalsarray-like Eigenvalues of a diffusion tensor. axisint Axis of evals which contains 3 eigenvalues. Returns linearityarray Calculated linearity of the diffusion tensor. Notes Linearity is calculated with the following equation: $Linearity = \frac{\lambda_1-\lambda_2}{\lambda_1+\lambda_2+\lambda_3}$ References [1] Westin C.-F., Peled S., Gubjartsson H., Kikinis R., Jolesz F., “Geometrical diffusion measures for MRI from tensor basis analysis” in Proc. 5th Annual ISMRM, 1997. ### lower_triangular 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 Parameters tensorarray_like (…, 3, 3) a collection of 3, 3 diffusion tensors b0float if b0 is not none log(b0) is returned as the dummy variable Returns Dndarray If b0 is none, then the shape will be (…, 6) otherwise (…, 7) ### mean_diffusivity dipy.reconst.dti.mean_diffusivity(evals, axis=-1) Mean Diffusivity (MD) of a diffusion tensor. Parameters evalsarray-like Eigenvalues of a diffusion tensor. axisint Axis of evals which contains 3 eigenvalues. Returns mdarray Calculated MD. Notes MD is calculated with the following equation: $MD = \frac{\lambda_1 + \lambda_2 + \lambda_3}{3}$ ### mode dipy.reconst.dti.mode(q_form) Mode (MO) of a diffusion tensor [1]. Parameters q_formndarray The quadratic form of a tensor, or an array with quadratic forms of tensors. Should be of shape (x, y, z, 3, 3) or (n, 3, 3) or (3, 3). Returns modearray Calculated tensor mode in each spatial coordinate. Notes Mode ranges between -1 (planar anisotropy) and +1 (linear anisotropy) with 0 representing orthotropy. Mode is calculated with the following equation (equation 9 in [1]): $Mode = 3*\sqrt{6}*det(\widetilde{A}/norm(\widetilde{A}))$ Where $$\widetilde{A}$$ is the deviatoric part of the tensor quadratic form. References 1(1,2,3,4) Daniel B. Ennis and G. Kindlmann, “Orthogonal Tensor Invariants and the Analysis of Diffusion Tensor Magnetic Resonance Images”, Magnetic Resonance in Medicine, vol. 55, no. 1, pp. 136-146, 2006. ### nlls_fit_tensor dipy.reconst.dti.nlls_fit_tensor(design_matrix, data, weighting=None, sigma=None, jac=True, return_S0_hat=False) Fit the tensor params using non-linear least-squares. Parameters design_matrixarray (g, 7) Design matrix holding the covariants used to solve for the regression coefficients. dataarray ([X, Y, Z, …], g) Data or response variables holding the data. Note that the last dimension should contain the data. It makes no copies of data. weighting: str the weighting scheme to use in considering the squared-error. Default behavior is to use uniform weighting. Other options: ‘sigma’ ‘gmm’ sigma: float If the ‘sigma’ weighting scheme is used, a value of sigma needs to be provided here. According to [Chang2005], a good value to use is 1.5267 * std(background_noise), where background_noise is estimated from some part of the image known to contain no signal (only noise). jacbool Use the Jacobian? Default: True return_S0_hatbool Boolean to return (True) or not (False) the S0 values for the fit. Returns nlls_params: the eigen-values and eigen-vectors of the tensor in each voxel. ### norm dipy.reconst.dti.norm(q_form) Calculate the Frobenius norm of a tensor quadratic form Parameters q_form: ndarray The quadratic form of a tensor, or an array with quadratic forms of tensors. Should be of shape (x,y,z,3,3) or (n, 3, 3) or (3,3). Returns normndarray The Frobenius norm of the 3,3 tensor q_form in each spatial coordinate. See also np.linalg.norm Notes The Frobenius norm is defined as: Math ||A||_F = [sum_{i,j} abs(a_{i,j})^2]^{1/2} ### ols_fit_tensor 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 1. Parameters design_matrixarray (g, 7) Design matrix holding the covariants used to solve for the regression coefficients. dataarray ([X, Y, Z, …], g) Data or response variables holding the data. Note that the last dimension should contain the data. It makes no copies of data. return_S0_hatbool Boolean to return (True) or not (False) the S0 values for the fit. Returns eigvalsarray (…, 3) Eigenvalues from eigen decomposition of the tensor. eigvecsarray (…, 3, 3) Associated eigenvectors from eigen decomposition of the tensor. Eigenvectors are columnar (e.g. eigvecs[:,j] is associated with eigvals[j]) See also WLS_fit_tensor, decompose_tensor, design_matrix Notes \begin{align}\begin{aligned}\begin{split}y = \mathrm{data} \\ X = \mathrm{design matrix} \\\end{split}\\\hat{\beta}_\mathrm{OLS} = (X^T X)^{-1} X^T y\end{aligned}\end{align} References 1(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28) Chung, SW., Lu, Y., Henry, R.G., 2006. Comparison of bootstrap approaches for estimation of uncertainties of DTI parameters. NeuroImage 33, 531-541. ### pinv 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. Parameters aarray_like (…, M, N) Matrix to be pseudo-inverted. rcondfloat Cutoff for small singular values. Returns Bndarray (…, N, M) The pseudo-inverse of a. Raises LinAlgError If the SVD computation does not converge. See also np.linalg.pinv ### planarity dipy.reconst.dti.planarity(evals, axis=-1) The planarity of the tensor 1 Parameters evalsarray-like Eigenvalues of a diffusion tensor. axisint Axis of evals which contains 3 eigenvalues. Returns linearityarray Calculated linearity of the diffusion tensor. Notes Planarity is calculated with the following equation: $Planarity = \frac{2 (\lambda_2-\lambda_3)}{\lambda_1+\lambda_2+\lambda_3}$ References [1] Westin C.-F., Peled S., Gubjartsson H., Kikinis R., Jolesz F., “Geometrical diffusion measures for MRI from tensor basis analysis” in Proc. 5th Annual ISMRM, 1997. ### quantize_evecs dipy.reconst.dti.quantize_evecs(evecs, odf_vertices=None) Find the closest orientation of an evenly distributed sphere Parameters evecsndarray odf_verticesNone or ndarray If None, then set vertices from symmetric362 sphere. Otherwise use passed ndarray as vertices Returns INndarray ### radial_diffusivity dipy.reconst.dti.radial_diffusivity(evals, axis=-1) Radial Diffusivity (RD) of a diffusion tensor. Also called perpendicular diffusivity. Parameters evalsarray-like Eigenvalues of a diffusion tensor, must be sorted in descending order along axis. axisint Axis of evals which contains 3 eigenvalues. Returns rdarray Calculated RD. Notes RD is calculated with the following equation: $RD = \frac{\lambda_2 + \lambda_3}{2}$ ### restore_fit_tensor dipy.reconst.dti.restore_fit_tensor(design_matrix, data, sigma=None, jac=True, return_S0_hat=False) Use the RESTORE algorithm [Chang2005] to calculate a robust tensor fit Parameters design_matrixarray of shape (g, 7) Design matrix holding the covariants used to solve for the regression coefficients. dataarray of shape ([X, Y, Z, n_directions], g) Data or response variables holding the data. Note that the last dimension should contain the data. It makes no copies of data. sigmafloat An estimate of the variance. [Chang2005] 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). jacbool, optional Whether to use the Jacobian of the tensor to speed the non-linear optimization procedure used to fit the tensor parameters (see also nlls_fit_tensor()). Default: True return_S0_hatbool Boolean to return (True) or not (False) the S0 values for the fit. Returns restore_paramsan estimate of the tensor parameters in each voxel. References Chang, L-C, Jones, DK and Pierpaoli, C (2005). RESTORE: robust estimation of tensors by outlier rejection. MRM, 53: 1088-95. ### sphericity dipy.reconst.dti.sphericity(evals, axis=-1) The sphericity of the tensor 1 Parameters evalsarray-like Eigenvalues of a diffusion tensor. axisint Axis of evals which contains 3 eigenvalues. Returns sphericityarray Calculated sphericity of the diffusion tensor. Notes Sphericity is calculated with the following equation: $Sphericity = \frac{3 \lambda_3)}{\lambda_1+\lambda_2+\lambda_3}$ References [1] Westin C.-F., Peled S., Gubjartsson H., Kikinis R., Jolesz F., “Geometrical diffusion measures for MRI from tensor basis analysis” in Proc. 5th Annual ISMRM, 1997. ### tensor_prediction dipy.reconst.dti.tensor_prediction(dti_params, gtab, S0) Predict a signal given tensor parameters. Parameters dti_paramsndarray Tensor parameters. The last dimension should have 12 tensor parameters: 3 eigenvalues, followed by the 3 corresponding eigenvectors. gtaba GradientTable class instance The gradient table for this prediction S0float or ndarray The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1 Notes The predicted signal is given by: $$S( 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. ### trace dipy.reconst.dti.trace(evals, axis=-1) Trace of a diffusion tensor. Parameters evalsarray-like Eigenvalues of a diffusion tensor. axisint Axis of evals which contains 3 eigenvalues. Returns tracearray Calculated trace of the diffusion tensor. Notes Trace is calculated with the following equation: $Trace = \lambda_1 + \lambda_2 + \lambda_3$ ### vec_val_vect dipy.reconst.dti.vec_val_vect() Vectorize vecs.diag(vals).vecs.T for last 2 dimensions of vecs Parameters vecsshape (…, M, N) array containing tensor in last two dimensions; M, N usually equal to (3, 3) valsshape (…, N) array diagonal values carried in last dimension, ... shape above must match that for vecs Returns resshape (…, M, M) array For all the dimensions ellided by ..., loops to get (M, N) vec matrix, and (N,) vals vector, and calculates vec.dot(np.diag(val).dot(vec.T). Raises ValueErrornon-matching ... dimensions of vecs, vals ValueErrornon-matching N dimensions of vecs, vals Examples Make a 3D array where the first dimension is only 1 >>> vecs = np.arange(9).reshape((1, 3, 3)) >>> vals = np.arange(3).reshape((1, 3)) >>> vec_val_vect(vecs, vals) array([[[ 9., 24., 39.], [ 24., 66., 108.], [ 39., 108., 177.]]])  That’s the same as the 2D case (apart from the float casting): >>> vecs = np.arange(9).reshape((3, 3)) >>> vals = np.arange(3) >>> np.dot(vecs, np.dot(np.diag(vals), vecs.T)) array([[ 9, 24, 39], [ 24, 66, 108], [ 39, 108, 177]])  ### vector_norm dipy.reconst.dti.vector_norm(vec, axis=-1, keepdims=False) Return vector Euclidean (L2) norm Parameters vecarray_like Vectors to norm. axisint Axis over which to norm. By default norm over last axis. If axis is None, vec is flattened then normed. keepdimsbool If True, the output will have the same number of dimensions as vec, with shape 1 on axis. Returns normarray 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.])  ### wls_fit_tensor 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]. Parameters design_matrixarray (g, 7) Design matrix holding the covariants used to solve for the regression coefficients. dataarray ([X, Y, Z, …], g) Data or response variables holding the data. Note that the last dimension should contain the data. It makes no copies of data. return_S0_hatbool Boolean to return (True) or not (False) the S0 values for the fit. Returns eigvalsarray (…, 3) Eigenvalues from eigen decomposition of the tensor. eigvecsarray (…, 3, 3) Associated eigenvectors from eigen decomposition of the tensor. Eigenvectors are columnar (e.g. eigvecs[:,j] is associated with eigvals[j]) 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: 1. calculate OLS estimates of the data 2. apply the OLS estimates as weights to the WLS fit of the data This ensured heteroscedasticity could be properly modeled for various types of bootstrap resampling (namely residual bootstrap). $\begin{split}y = \mathrm{data} \\ X = \mathrm{design matrix} \\ \hat{\beta}_\mathrm{WLS} = \mathrm{desired regression coefficients (e.g. tensor)}\\ \\ \hat{\beta}_\mathrm{WLS} = (X^T W X)^{-1} X^T W y \\ \\ W = \mathrm{diag}((X \hat{\beta}_\mathrm{OLS})^2), \mathrm{where} \hat{\beta}_\mathrm{OLS} = (X^T X)^{-1} X^T y\end{split}$ References 1(1,2,3) Chung, SW., Lu, Y., Henry, R.G., 2006. Comparison of bootstrap approaches for estimation of uncertainties of DTI parameters. NeuroImage 33, 531-541. ### Cache class dipy.reconst.forecast.Cache Bases: object Cache values based on a key object (such as a sphere or gradient table). Notes This class is meant to be used as a mix-in: class MyModel(Model, Cache): pass class MyModelFit(Fit): pass  Inside a method on the fit, typical usage would be: def odf(sphere): M = self.model.cache_get('odf_basis_matrix', key=sphere) if M is None: M = self._compute_basis_matrix(sphere) self.model.cache_set('odf_basis_matrix', key=sphere, value=M)  Methods  Clear the cache. cache_get(tag, key[, default]) Retrieve a value from the cache. cache_set(tag, key, value) Store a value in the cache. __init__(self, /, *args, **kwargs)

Initialize self. See help(type(self)) for accurate signature.

cache_clear()

Clear the cache.

cache_get(tag, key, default=None)

Retrieve a value from the cache.

Parameters
tagstr

Description of the cached value.

keyobject

Key object used to look up the cached value.

defaultobject

Value to be returned if no cached entry is found.

Returns
vobject

Value from the cache associated with (tag, key). Returns default if no cached entry is found.

cache_set(tag, key, value)

Store a value in the cache.

Parameters
tagstr

Description of the cached value.

keyobject

Key object used to look up the cached value.

valueobject

Value stored in the cache for each unique combination of (tag, key).

Examples

>>> def compute_expensive_matrix(parameters):
...     # Imagine the following computation is very expensive
...     return (p**2 for p in parameters)

>>> 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


### ForecastFit

class dipy.reconst.forecast.ForecastFit(model, data, sh_coef, d_par, d_perp)
Attributes
dpar

The parallel diffusivity

dperp

The perpendicular diffusivity

sh_coeff

The FORECAST SH coefficients

Methods

 Calculates the fractional anisotropy. Calculates the mean diffusivity. odf(sphere[, clip_negative]) Calculates the fODF for a given discrete sphere. predict([gtab, S0]) Calculates the fODF for a given discrete sphere.
__init__(model, data, sh_coef, d_par, d_perp)

Calculates diffusion properties for a single voxel

Parameters
modelobject,

AnalyticalModel

data1d ndarray,

fitted data

sh_coef1d ndarray,

forecast sh coefficients

d_parfloat,

parallel diffusivity

d_perpfloat,

perpendicular diffusivity

property dpar

The parallel diffusivity

property dperp

The perpendicular diffusivity

fractional_anisotropy()

Calculates the fractional anisotropy.

mean_diffusivity()

Calculates the mean diffusivity.

odf(sphere, clip_negative=True)

Calculates the fODF for a given discrete sphere.

Parameters
sphereSphere,

the odf sphere

clip_negativeboolean, optional

if True clip the negative odf values to 0, default True

predict(gtab=None, S0=1.0)

Calculates the fODF for a given discrete sphere.

Parameters

gradient directions and bvalues container class.

S0float, optional

the signal at b-value=0

property sh_coeff

The FORECAST SH coefficients

### ForecastModel

class dipy.reconst.forecast.ForecastModel(gtab, sh_order=8, lambda_lb=0.001, dec_alg='CSD', sphere=None, lambda_csd=1.0)

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]

References

1

Anderson A. W., “Measurement of Fiber Orientation Distributions Using High Angular Resolution Diffusion Imaging”, Magnetic Resonance in Medicine, 2005.

2(1,2)

Kaden E. et al., “Quantitative Mapping of the Per-Axon Diffusion Coefficients in Brain White Matter”, Magnetic Resonance in Medicine, 2016.

3

Zucchelli E. et al., “A generalized SMT-based framework for Diffusion MRI microstructural model estimation”, MICCAI Workshop on Computational DIFFUSION MRI (CDMRI), 2017.

The implementation of FORECAST may require CVXPY (http://www.cvxpy.org/).

Methods

 cache_clear() Clear the cache. cache_get(tag, key[, default]) Retrieve a value from the cache. cache_set(tag, key, value) Store a value in the cache. fit(data[, mask]) Fit method for every voxel in data
__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})$$

..math::
nowrap
begin{equation}

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.

Parameters

gradient directions and bvalues container class.

sh_orderunsigned int,

an even integer that represent the SH order of the basis (max 12)

lambda_lb: float,

Laplace-Beltrami regularization weight.

dec_algstr,

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’.

spherearray, shape (N,3),

sphere points where to enforce positivity when ‘POS’ or ‘CSD’ dec_alg are selected.

lambda_csdfloat,

CSD regularization weight.

References

1(1,2)

Anderson A. W., “Measurement of Fiber Orientation Distributions Using High Angular Resolution Diffusion Imaging”, Magnetic Resonance in Medicine, 2005.

2

Kaden E. et al., “Quantitative Mapping of the Per-Axon Diffusion Coefficients in Brain White Matter”, Magnetic Resonance in Medicine, 2016.

3

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)

fit(data, mask=None)

Fit method for every voxel in data

### OdfFit

class dipy.reconst.forecast.OdfFit(model, data)

Methods

 odf(sphere) To be implemented but specific odf models
__init__(model, data)

Initialize self. See help(type(self)) for accurate signature.

odf(sphere)

To be implemented but specific odf models

### OdfModel

class dipy.reconst.forecast.OdfModel(gtab)

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

 fit(data) To be implemented by specific odf models
__init__(gtab)

Initialization of the abstract class for signal reconstruction models

Parameters
fit(data)

To be implemented by specific odf models

### cart2sphere

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$$

Parameters
xarray_like

x coordinate in Cartesian space

yarray_like

y coordinate in Cartesian space

zarray_like

z coordinate

Returns
rarray

thetaarray

inclination (polar) angle

phiarray

azimuth angle

### csdeconv

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.

Parameters
dwsignalarray

Diffusion weighted signals to be deconvolved.

Xarray

Prediction matrix which estimates diffusion weighted signals from FOD coefficients.

B_regarray (N, B)

SH basis matrix which maps FOD coefficients to FOD values on the surface of the sphere. B_reg should be scaled to account for lambda.

taufloat

Threshold controlling the amplitude below which the corresponding fODF is assumed to be zero. Ideally, tau should be set to zero. However, to improve the stability of the algorithm, tau is set to tau*100 % of the max fODF amplitude (here, 10% by default). This is similar to peak detection where peaks below 0.1 amplitude are usually considered noise peaks. Because SDT is based on a q-ball ODF deconvolution, and not signal deconvolution, using the max instead of mean (as in CSD), is more stable.

convergenceint

Maximum number of iterations to allow the deconvolution to converge.

Pndarray

This is an optimization to avoid computing dot(X.T, X) many times. If the same X is used many times, P can be precomputed and passed to this function.

Returns
fodf_shndarray ((sh_order + 1)*(sh_order + 2)/2,)

Spherical harmonics coefficients of the constrained-regularized fiber ODF.

num_itint

Number of iterations in the constrained-regularization used for convergence.

Notes

This section describes how the fitting of the SH coefficients is done. Problem is to minimise per iteration:

$$F(f_n) = ||Xf_n - S||^2 + \lambda^2 ||H_{n-1} f_n||^2$$

Where $$X$$ maps current FOD SH coefficients $$f_n$$ to DW signals $$s$$ and $$H_{n-1}$$ maps FOD SH coefficients $$f_n$$ to amplitudes along set of negative directions identified in previous iteration, i.e. the matrix formed by the rows of $$B_{reg}$$ for which $$Hf_{n-1}<0$$ where $$B_{reg}$$ maps $$f_n$$ to FOD amplitude on a sphere.

Solve by differentiating and setting to zero:

$$\Rightarrow \frac{\delta F}{\delta f_n} = 2X^T(Xf_n - S) + 2 \lambda^2 H_{n-1}^TH_{n-1}f_n=0$$

Or:

$$(X^TX + \lambda^2 H_{n-1}^TH_{n-1})f_n = X^Ts$$

Define $$Q = X^TX + \lambda^2 H_{n-1}^TH_{n-1}$$ , which by construction is a square positive definite symmetric matrix of size $$n_{SH} by n_{SH}$$. If needed, positive definiteness can be enforced with a small minimum norm regulariser (helps a lot with poorly conditioned direction sets and/or superresolution):

$$Q = X^TX + (\lambda H_{n-1}^T) (\lambda H_{n-1}) + \mu I$$

Solve $$Qf_n = X^Ts$$ using Cholesky decomposition:

$$Q = LL^T$$

where $$L$$ is lower triangular. Then problem can be solved by back-substitution:

$$L_y = X^Ts$$

$$L^Tf_n = y$$

To speeds things up further, form $$P = X^TX + \mu I$$, and update to form $$Q$$ by rankn update with $$H_{n-1}$$. The dipy implementation looks like:

form initially $$P = X^T X + \mu I$$ and $$\lambda B_{reg}$$

for each voxel: form $$z = X^Ts$$

estimate $$f_0$$ by solving $$Pf_0=z$$. We use a simplified $$l_{max}=4$$ solution here, but it might not make a big difference.

Then iterate until no change in rows of $$H$$ used in $$H_n$$

form $$H_{n}$$ given $$f_{n-1}$$

form $$Q = P + (\lambda H_{n-1}^T) (\lambda H_{n-1}$$) (this can be done by rankn update, but we currently do not use rankn update).

solve $$Qf_n = z$$ using Cholesky decomposition

We’d like to thanks Donald Tournier for his help with describing and implementing this algorithm.

References

1(1,2,3)

Tournier, J.D., et al. NeuroImage 2007. Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution.

### find_signal_means

dipy.reconst.forecast.find_signal_means(b_unique, data_norm, bvals, rho, lb_matrix, w=0.001)

Calculate the mean signal for each shell.

Parameters
b_unique1d ndarray,

unique b-values in a vector excluding zero

data_norm1d ndarray,

normalized diffusion signal

bvals1d ndarray,

the b-values

rho2d ndarray,

SH basis matrix for fitting the signal on each shell

lb_matrix2d ndarray,

Laplace-Beltrami regularization matrix

wfloat,

weight for the Laplace-Beltrami regularization

Returns
means1d ndarray

the average of the signal for each b-values

### forecast_error_func

dipy.reconst.forecast.forecast_error_func(x, b_unique, E)

Calculates the difference between the mean signal calculated using the parameter vector x and the average signal E using FORECAST and SMT

### forecast_matrix

dipy.reconst.forecast.forecast_matrix(sh_order, d_par, d_perp, bvals)

Compute the FORECAST radial matrix

### lb_forecast

dipy.reconst.forecast.lb_forecast(sh_order)

Returns the Laplace-Beltrami regularization matrix for FORECAST

### leastsq

dipy.reconst.forecast.leastsq(func, x0, args=(), Dfun=None, full_output=0, col_deriv=0, ftol=1.49012e-08, xtol=1.49012e-08, gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None)

Minimize the sum of squares of a set of equations.

x = arg min(sum(func(y)**2,axis=0))
y

Parameters
funccallable

should take at least one (possibly length N vector) argument and returns M floating point numbers. It must not return NaNs or fitting might fail.

x0ndarray

The starting estimate for the minimization.

argstuple, optional

Any extra arguments to func are placed in this tuple.

Dfuncallable, optional

A function or method to compute the Jacobian of func with derivatives across the rows. If this is None, the Jacobian will be estimated.

full_outputbool, optional

non-zero to return all optional outputs.

col_derivbool, optional

non-zero to specify that the Jacobian function computes derivatives down the columns (faster, because there is no transpose operation).

ftolfloat, optional

Relative error desired in the sum of squares.

xtolfloat, optional

Relative error desired in the approximate solution.

gtolfloat, optional

Orthogonality desired between the function vector and the columns of the Jacobian.

maxfevint, optional

The maximum number of calls to the function. If Dfun is provided then the default maxfev is 100*(N+1) where N is the number of elements in x0, otherwise the default maxfev is 200*(N+1).

epsfcnfloat, optional

A variable used in determining a suitable step length for the forward- difference approximation of the Jacobian (for Dfun=None). Normally the actual step length will be sqrt(epsfcn)*x If epsfcn is less than the machine precision, it is assumed that the relative errors are of the order of the machine precision.

factorfloat, optional

A parameter determining the initial step bound (factor * || diag * x||). Should be in interval (0.1, 100).

diagsequence, optional

N positive entries that serve as a scale factors for the variables.

Returns
xndarray

The solution (or the result of the last iteration for an unsuccessful call).

cov_xndarray

Uses the fjac and ipvt optional outputs to construct an estimate of the jacobian around the solution. None if a singular matrix encountered (indicates very flat curvature in some direction). This matrix must be multiplied by the residual variance to get the covariance of the parameter estimates – see curve_fit.

infodictdict

a dictionary of optional outputs with the key s:

nfev

The number of function calls

fvec

The function evaluated at the output

fjac

A permutation of the R matrix of a QR factorization of the final approximate Jacobian matrix, stored column wise. Together with ipvt, the covariance of the estimate can be approximated.

ipvt

An integer array of length N which defines a permutation matrix, p, such that fjac*p = q*r, where r is upper triangular with diagonal elements of nonincreasing magnitude. Column j of p is column ipvt(j) of the identity matrix.

qtf

The vector (transpose(q) * fvec).

mesgstr

A string message giving information about the cause of failure.

ierint

An integer flag. If it is equal to 1, 2, 3 or 4, the solution was found. Otherwise, the solution was not found. In either case, the optional output variable ‘mesg’ gives more information.

Notes

“leastsq” is a wrapper around MINPACK’s lmdif and lmder algorithms.

cov_x is a Jacobian approximation to the Hessian of the least squares objective function. This approximation assumes that the objective function is based on the difference between some observed target data (ydata) and a (non-linear) function of the parameters f(xdata, params)

func(params) = ydata - f(xdata, params)


so that the objective function is

  min   sum((ydata - f(xdata, params))**2, axis=0)
params


The solution, x, is always a 1D array, regardless of the shape of x0, or whether x0 is a scalar.

### multi_voxel_fit

dipy.reconst.forecast.multi_voxel_fit(single_voxel_fit)

Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition

### optional_package

dipy.reconst.forecast.optional_package(name, trip_msg=None)

Return package-like thing and module setup for package name

Parameters
namestr

package name

trip_msgNone or str

message to give when someone tries to use the return package, but we could not import it, and have returned a TripWire object instead. Default message if None.

Returns
pkg_likemodule or TripWire instance

If we can import the package, return it. Otherwise return an object raising an error when accessed

have_pkgbool

True if import for package was successful, false otherwise

module_setupfunction

callable usually set as setup_module in calling namespace, to allow skipping tests.

Examples

Typical use would be something like this at the top of a module using an optional package:

>>> from dipy.utils.optpkg import optional_package
>>> pkg, have_pkg, setup_module = optional_package('not_a_package')


Of course in this case the package doesn’t exist, and so, in the module:

>>> have_pkg
False


and

>>> pkg.some_function() #doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
...
TripWireError: We need package not_a_package for these functions, but
import not_a_package raised an ImportError


If the module does exist - we get the module

>>> pkg, _, _ = optional_package('os')
>>> hasattr(pkg, 'path')
True


Or a submodule if that’s what we asked for

>>> subpkg, _, _ = optional_package('os.path')
>>> hasattr(subpkg, 'dirname')
True


### psi_l

dipy.reconst.forecast.psi_l(l, b)

### real_sph_harm

dipy.reconst.forecast.real_sph_harm(m, n, theta, phi)

Compute real spherical harmonics.

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.

Parameters
mint |m| <= n

The order of the harmonic.

nint >= 0

The degree of the harmonic.

thetafloat [0, 2*pi]

The azimuthal (longitudinal) coordinate.

phifloat [0, pi]

The polar (colatitudinal) coordinate.

Returns
y_mnreal float

The real harmonic $$Y^m_n$$ sampled at theta and phi.

scipy.special.sph_harm

### rho_matrix

dipy.reconst.forecast.rho_matrix(sh_order, vecs)

Compute the SH matrix $$\rho$$

### warn

dipy.reconst.forecast.warn()

Issue a warning, or maybe ignore it or raise an exception.

### FreeWaterTensorFit

class dipy.reconst.fwdti.FreeWaterTensorFit(model, model_params)

Class for fitting the Free Water Tensor Model

Attributes
S0_hat
directions

For tracking - return the primary direction in each voxel

evals

Returns the eigenvalues of the tensor as an array

evecs

Returns the eigenvectors of the tensor as an array, columnwise

f

Returns the free water diffusion volume fraction f

quadratic_form

Calculates the 3x3 diffusion tensor for each voxel

shape

Methods

 ad() Axial diffusivity (AD) calculated from cached eigenvalues. adc(sphere) Calculate the apparent diffusion coefficient (ADC) in each direction on color_fa() Color fractional anisotropy of diffusion tensor fa() Fractional anisotropy (FA) calculated from cached eigenvalues. ga() Geodesic anisotropy (GA) calculated from cached eigenvalues. linearity() Returns md() Mean diffusivity (MD) calculated from cached eigenvalues. mode() Tensor mode calculated from cached eigenvalues. odf(sphere) The diffusion orientation distribution function (dODF). planarity() Returns predict(gtab[, S0]) Given a free water tensor model fit, predict the signal on the vertices of a gradient table rd() Radial diffusivity (RD) calculated from cached eigenvalues. sphericity() Returns trace() Trace of the tensor calculated from cached eigenvalues.
 lower_triangular
__init__(model, model_params)

Initialize a FreeWaterTensorFit class instance. Since the free water tensor model is an extension of DTI, class instance is defined as subclass of the TensorFit from dti.py

Parameters
modelFreeWaterTensorModel Class instance

Class instance containing the free water tensor model for the fit

model_paramsndarray (x, y, z, 13) or (n, 13)

All parameters estimated from the free water tensor model. Parameters are ordered as follows:

1. Three diffusion tensor’s eigenvalues

2. Three lines of the eigenvector matrix each containing the first, second and third coordinates of the eigenvector

3. The volume fraction of the free water compartment

References

1

Henriques, R.N., Rokem, A., Garyfallidis, E., St-Jean, S., Peterson E.T., Correia, M.M., 2017. [Re] Optimization of a free water elimination two-compartment model for diffusion tensor imaging. ReScience volume 3, issue 1, article number 2

property f

Returns the free water diffusion volume fraction f

predict(gtab, S0=1)

Given a free water tensor model fit, predict the signal on the vertices of a gradient table

Parameters
gtaba GradientTable class instance

The gradient table for this prediction

S0float array

The mean non-diffusion weighted signal in each voxel. Default: 1 in all voxels.

Returns
S(…, N) ndarray

Simulated signal based on the free water DTI model

### FreeWaterTensorModel

class dipy.reconst.fwdti.FreeWaterTensorModel(gtab, fit_method='NLS', *args, **kwargs)

Class for the Free Water Elimination Diffusion Tensor Model

Methods

 fit(data[, mask]) Fit method for every voxel in data predict(fwdti_params[, S0]) Predict a signal for this TensorModel class instance given parameters.
__init__(gtab, fit_method='NLS', *args, **kwargs)

Free Water Diffusion Tensor Model [1].

Parameters
fit_methodstr or callable

str can be one of the following:

‘WLS’ for weighted linear least square fit according to [1]

fwdti.wls_iter()

‘NLS’ for non-linear least square fit according to [1]

fwdti.nls_iter()

callable has to have the signature:

fit_method(design_matrix, data, *args, **kwargs)

args, kwargsarguments and key-word arguments passed to the

fit_method. See fwdti.wls_iter, fwdti.nls_iter for details

References

1(1,2,3,4)

Henriques, R.N., Rokem, A., Garyfallidis, E., St-Jean, S., Peterson E.T., Correia, M.M., 2017. [Re] Optimization of a free water elimination two-compartment model for diffusion tensor imaging. ReScience volume 3, issue 1, article number 2

fit(data, mask=None)

Fit method for every voxel in data

predict(fwdti_params, S0=1)

Predict a signal for this TensorModel class instance given parameters.

Parameters
fwdti_params(…, 13) ndarray

The last dimension should have 13 parameters: the 12 tensor parameters (3 eigenvalues, followed by the 3 corresponding eigenvectors) and the free water volume fraction.

S0float or ndarray

The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1

Returns
S(…, N) ndarray

Simulated signal based on the free water DTI model

### ReconstModel

class dipy.reconst.fwdti.ReconstModel(gtab)

Bases: object

Abstract class for signal reconstruction models

Methods

 fit
__init__(gtab)

Initialization of the abstract class for signal reconstruction models

Parameters
fit(data, mask=None, **kwargs)

### TensorFit

class dipy.reconst.fwdti.TensorFit(model, model_params, model_S0=None)

Bases: object

Attributes
S0_hat
directions

For tracking - return the primary direction in each voxel

evals

Returns the eigenvalues of the tensor as an array

evecs

Returns the eigenvectors of the tensor as an array, columnwise

quadratic_form

Calculates the 3x3 diffusion tensor for each voxel

shape

Methods

 Axial diffusivity (AD) calculated from cached eigenvalues. adc(sphere) 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. Returns Mean diffusivity (MD) calculated from cached eigenvalues. Tensor mode calculated from cached eigenvalues. odf(sphere) The diffusion orientation distribution function (dODF). Returns predict(gtab[, S0, step]) Given a model fit, predict the signal on the vertices of a sphere Radial diffusivity (RD) calculated from cached eigenvalues. Returns Trace of the tensor calculated from cached eigenvalues.
 lower_triangular
__init__(model, model_params, model_S0=None)

Initialize a TensorFit class instance.

property S0_hat
ad()

Axial diffusivity (AD) calculated from cached eigenvalues.

Returns

Notes

RD is calculated with the following equation:

$AD = \lambda_1$
adc(sphere)

Calculate the apparent diffusion coefficient (ADC) in each direction on the sphere for each voxel in the data

Parameters
sphereSphere class instance
Returns

The estimates of the apparent diffusion coefficient in every direction on the input sphere

ec{b} Q ec{b}^T

Where Q is the quadratic form of the tensor.

color_fa()

Color fractional anisotropy of diffusion tensor

property directions

For tracking - return the primary direction in each voxel

property evals

Returns the eigenvalues of the tensor as an array

property evecs

Returns the eigenvectors of the tensor as an array, columnwise

fa()

Fractional anisotropy (FA) calculated from cached eigenvalues.

ga()

Geodesic anisotropy (GA) calculated from cached eigenvalues.

linearity()
Returns
linearityarray

Calculated linearity of the diffusion tensor 1.

Notes

Linearity is calculated with the following equation:

$Linearity = \frac{\lambda_1-\lambda_2}{\lambda_1+\lambda_2+\lambda_3}$

References

[1] Westin C.-F., Peled S., Gubjartsson H., Kikinis R., Jolesz

F., “Geometrical diffusion measures for MRI from tensor basis analysis” in Proc. 5th Annual ISMRM, 1997.

lower_triangular(b0=None)
md()

Mean diffusivity (MD) calculated from cached eigenvalues.

Returns
mdarray (V, 1)

Calculated MD.

Notes

MD is calculated with the following equation:

$MD = \frac{\lambda_1+\lambda_2+\lambda_3}{3}$
mode()

Tensor mode calculated from cached eigenvalues.

odf(sphere)

The diffusion orientation distribution function (dODF). This is an estimate of the diffusion distance in each direction

Parameters
sphereSphere class instance.

The dODF is calculated in the vertices of this input.

Returns
odfndarray

The diffusion distance in every direction of the sphere in every voxel in the input data.

Notes

This is based on equation 3 in [Aganj2010]. To re-derive it from scratch, follow steps in [Descoteaux2008], Section 7.9 Equation 7.24 but with an $$r^2$$ term in the integral.

References

Aganj2010(1,2)

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

Descoteaux2008(1,2)

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()
Returns
sphericityarray

Calculated sphericity of the diffusion tensor 1.

Notes

Sphericity is calculated with the following equation:

$Sphericity = \frac{2 (\lambda_2 - \lambda_3)}{\lambda_1+\lambda_2+\lambda_3}$

References

[1] Westin C.-F., Peled S., Gubjartsson H., Kikinis R., Jolesz

F., “Geometrical diffusion measures for MRI from tensor basis analysis” in Proc. 5th Annual ISMRM, 1997.

predict(gtab, S0=None, step=None)

Given a model fit, predict the signal on the vertices of a sphere

Parameters
gtaba GradientTable class instance

This encodes the directions for which a prediction is made

S0float array

The mean non-diffusion weighted signal in each voxel. Default: The fitted S0 value in all voxels if it was fitted. Otherwise 1 in all voxels.

stepint

The chunk size as a number of voxels. Optional parameter with default value 10,000.

In order to increase speed of processing, tensor fitting is done simultaneously over many voxels. This parameter sets the number of voxels that will be fit at once in each iteration. A larger step value should speed things up, but it will also take up more memory. It is advisable to keep an eye on memory consumption as this value is increased.

Notes

The predicted signal is given by:

$S( heta, b) = S_0 * e^{-b ADC}$

Where: .. math

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

property quadratic_form

Calculates the 3x3 diffusion tensor for each voxel

rd()

Radial diffusivity (RD) calculated from cached eigenvalues.

Returns
rdarray (V, 1)

Calculated RD.

Notes

RD is calculated with the following equation:

$RD = \frac{\lambda_2 + \lambda_3}{2}$
property shape
sphericity()
Returns
sphericityarray

Calculated sphericity of the diffusion tensor 1.

Notes

Sphericity is calculated with the following equation:

$Sphericity = \frac{3 \lambda_3}{\lambda_1+\lambda_2+\lambda_3}$

References

[1] Westin C.-F., Peled S., Gubjartsson H., Kikinis R., Jolesz

F., “Geometrical diffusion measures for MRI from tensor basis analysis” in Proc. 5th Annual ISMRM, 1997.

trace()

Trace of the tensor calculated from cached eigenvalues.

Returns
tracearray (V, 1)

Calculated trace.

Notes

The trace is calculated with the following equation:

$trace = \lambda_1 + \lambda_2 + \lambda_3$

### check_multi_b

dipy.reconst.fwdti.check_multi_b(gtab, n_bvals, non_zero=True, bmag=None)

Check if you have enough different b-values in your gradient table

Parameters
n_bvalsint

The number of different b-values you are checking for.

non_zerobool

Whether to check only non-zero bvalues. In this case, we will require at least n_bvals non-zero b-values (where non-zero is defined depending on the gtab object’s b0_threshold attribute)

bmagint

The order of magnitude of the b-values used. The function will normalize the b-values relative $$10^{bmag}$$. Default: derive this value from the maximal b-value provided: $$bmag=log_{10}(max(bvals)) - 1$$.

Returns
boolWhether there are at least n_bvals different b-values in the

### cholesky_to_lower_triangular

dipy.reconst.fwdti.cholesky_to_lower_triangular(R)

Convert Cholesky decompostion elements to the diffusion tensor elements

Parameters
Rarray (6,)

Array containing the six Cholesky’s decomposition elements (R0, R1, R2, R3, R4, R5) [1].

Returns
tensor_elementsarray (6,)

Array containing the six elements of diffusion tensor’s lower triangular.

References

1(1,2)

Koay, C.G., Carew, J.D., Alexander, A.L., Basser, P.J., Meyerand, M.E., 2006. Investigation of anomalous estimates of tensor-derived quantities in diffusion tensor imaging. Magnetic Resonance in Medicine, 55(4), 930-936. doi:10.1002/mrm.20832

### decompose_tensor

dipy.reconst.fwdti.decompose_tensor(tensor, min_diffusivity=0)

Returns eigenvalues and eigenvectors given a diffusion tensor

Computes tensor eigen decomposition to calculate eigenvalues and eigenvectors (Basser et al., 1994a).

Parameters
tensorarray (…, 3, 3)

Hermitian matrix representing a diffusion tensor.

min_diffusivityfloat

Because negative eigenvalues are not physical and small eigenvalues, much smaller than the diffusion weighting, cause quite a lot of noise in metrics such as fa, diffusivity values smaller than min_diffusivity are replaced with min_diffusivity.

Returns
eigvalsarray (…, 3)

Eigenvalues from eigen decomposition of the tensor. Negative eigenvalues are replaced by zero. Sorted from largest to smallest.

eigvecsarray (…, 3, 3)

Associated eigenvectors from eigen decomposition of the tensor. Eigenvectors are columnar (e.g. eigvecs[…, :, j] is associated with eigvals[…, j])

### design_matrix

dipy.reconst.fwdti.design_matrix(gtab, dtype=None)

Constructs design matrix for DTI weighted least squares or least squares fitting. (Basser et al., 1994a)

Parameters
gtabA GradientTable class instance
dtypestring

Parameter to control the dtype of returned designed matrix

Returns
design_matrixarray (g,7)

Design matrix or B matrix assuming Gaussian distributed tensor model design_matrix[j, :] = (Bxx, Byy, Bzz, Bxy, Bxz, Byz, dummy)

### from_lower_triangular

dipy.reconst.fwdti.from_lower_triangular(D)

Returns a tensor given the six unique tensor elements

Given the six unique tensor elements (in the order: Dxx, Dxy, Dyy, Dxz, Dyz, Dzz) returns a 3 by 3 tensor. All elements after the sixth are ignored.

Parameters
Darray_like, (…, >6)

Unique elements of the tensors

Returns
tensorndarray (…, 3, 3)

3 by 3 tensors

### fwdti_prediction

dipy.reconst.fwdti.fwdti_prediction(params, gtab, S0=1, Diso=0.003)

Signal prediction given the free water DTI model parameters.

Parameters
params(…, 13) ndarray

Model parameters. The last dimension should have the 12 tensor parameters (3 eigenvalues, followed by the 3 corresponding eigenvectors) and the volume fraction of the free water compartment.

gtaba GradientTable class instance

The gradient table for this prediction

S0float or ndarray

The non diffusion-weighted signal in every voxel, or across all voxels. Default: 1

Disofloat, optional

Value of the free water isotropic diffusion. Default is set to 3e-3 $$mm^{2}.s^{-1}$$. Please adjust this value if you are assuming different units of diffusion.

Returns
S(…, N) ndarray

Simulated signal based on the free water DTI model

Notes

The predicted signal is given by: $$S(\theta, b) = S_0 * [(1-f) * e^{-b ADC} + f * e^{-b D_{iso}]$$, where $$ADC = \theta Q \theta^T$$, $$\theta$$ is a unit vector pointing at any direction on the sphere for which a signal is to be predicted, $$b$$ is the b value provided in the GradientTable input for that direction, $$Q$$ is the quadratic form of the tensor determined by the input parameters, $$f$$ is the free water diffusion compartment, $$D_{iso}$$ is the free water diffusivity which is equal to \$3 * 10^{-3} mm^{2}s^{-1} [1].

References

1(1,2)

Henriques, R.N., Rokem, A., Garyfallidis, E., St-Jean, S., Peterson E.T., Correia, M.M., 2017. [Re] Optimization of a free water elimination two-compartment model for diffusion tensor imaging. ReScience volume 3, issue 1, article number 2

### lower_triangular

dipy.reconst.fwdti.lower_triangular(tensor, b0=None)

Returns the six lower triangular values of the tensor and a dummy variable if b0 is not None

Parameters
tensorarray_like (…, 3, 3)

a collection of 3, 3 diffusion tensors

b0float

if b0 is not none log(b0) is returned as the dummy variable

Returns
Dndarray

If b0 is none, then the shape will be (…, 6) otherwise (…, 7)

### lower_triangular_to_cholesky

dipy.reconst.fwdti.lower_triangular_to_cholesky(tensor_elements)

Performs Cholesky decomposition of the diffusion tensor

Parameters
tensor_elementsarray (6,)

Array containing the six elements of diffusion tensor’s lower triangular.

Returns
cholesky_elementsarray (6,)

Array containing the six Cholesky’s decomposition elements (R0, R1, R2, R3, R4, R5) [1].

References

1(1,2)

Koay, C.G., Carew, J.D., Alexander, A.L., Basser, P.J., Meyerand, M.E., 2006. Investigation of anomalous estimates of tensor-derived quantities in diffusion tensor imaging. Magnetic Resonance in Medicine, 55(4), 930-936. doi:10.1002/mrm.20832

### multi_voxel_fit

dipy.reconst.fwdti.multi_voxel_fit(single_voxel_fit)

Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition

### ndindex

dipy.reconst.fwdti.ndindex(shape)

An N-dimensional iterator object to index arrays.

Given the shape of an array, an