class designer.fitting.dwipy.DWI(imPath, bvecPath=None, bvalPath=None, mask=None, nthreads=-1)

Bases: object

The DWI object handles tensor estimation and parameter extraction of dwiffusion weighted images

Variables

  • hdr (class) – Nibabel class object of input nifti file

  • img (ndarray) – 3D or 4D input image array

  • grad (ndarray) – [N x 4] gradient and bvalue scheme, where the first three columns are the X, Y, and Z gradient vectors respectively, and the fourth column is B-values

  • mask (ndarray(dtype=bool)) – 3D array of brain mask

  • MaskStatus (bool) – True if brain_mask.nii is present, False otherwise

  • workers (int) – Number of CPU workers to use in processing

akccorrect(akc_out, window=3, connectivity='face')

Applies AKC outlier map to DT to replace outliers with a moving median. Run this only after tensor fitting and akc outlier detection.

Parameters

  • akc_out (ndarray(dtype=bool)) – 3D map containing outliers from DWI.akcoutliers

  • window (int, optional) – Width of square matrix filter (Default: 5)

  • connectivity (str, {‘face’, ‘all’}, optional) – Specifies what kind of connected-component connectivity to use for median determination

Examples

dwi.akccorrect(akc_out), where dwi is the DWI class object

akcoutliers(iter=10)

Uses 100,000 direction in chunks of 10 to iteratively find outliers. Returns a mask of locations where said violations occur. Multiprocessing is disabled because this is a memory-intensive task. To be run only after tensor fitting.

Parameters

iter (int, optional) – number of iterations to perform out of 10. Reduce this number if your computer does not have sufficient RAM. (Default: 10)

Returns

akc_out – 3D map containing outliers where AKC falls fails the inequality test -2 < AKC < 10

Return type

ndarray(dtype=bool)

Examples

akc_out = dwi.akoutliers(), where dwi is the DWI class object

createConstraints(constraints=[0, 1, 0])

Generates constraint array for constrained minimization quadratic programming

Parameters

constraints (array_like(dtype=int)) – [1 X 3] logical vector indicating which constraints out of three to enable (Default: [0, 1, 0]) C1 is Dapp > 0 C1 is Kapp > 0 C3 is Kapp < 3/(b*Dapp)

Returns

C – Array containing constraints to consider during minimization, C is shaped [number of constraints enforced * number of directions, 22]

Return type

ndarray(dtype=float)

Examples

C = dwi.createConstraints([0, 1, 0])

createTensorOrder(order=None)

Creates tensor order array and indices

Parameters

order (2 or 4 (int or None)) – Tensor order number, 2 for diffusion and 4 for kurtosis. Default: None; auto-detect

Returns

  • cnt (ndarray(dtype=int)) – Count of number of times a tensor appears

  • ind (ndarray(dtype=int)) – Indices of count

Examples

(cnt, ind) = dwi.createTensorOrder(order)

Notes

The tensors for this pipeline are based on NYU’s designer layout as depicted in the table below. This will soon be depreciated and updated with MRTRIX3’s layout.

~~~~~~D~~~~~~
1  |    D11
2  |    D12
3  |    D13
4  |    D22
5  |    D23
6  |    D33
~~~~~~K~~~~~~
1  |   W1111
2  |   W1112
3  |   W1113
4  |   W1122
5  |   W1123
6  |   W1133
7  |   W1222
8  |   W1223
9  |   W1233
10 |   W1333
11 |   W2222
12 |   W2223
13 |   W2233
14 |   W2333
15 |   W3333
diffusionCoeff(dt, dir)

Computes apparent diffusion coefficient (ADC)

Parameters

  • dt (ndarray(dtype=float)) – [21 x nvoxel] array containing diffusion tensor

  • dir (ndarray(dtype=float)) – [n x 3] array containing gradient directions

Returns

adc – Array containing apparent diffusion coefficient

Return type

ndarray(dtype=float)

Examples

adc = dwi.diffusionCoeff(dt, dir)

dkiTensorParams(v1, dt)

Uses average directional statistics to approximate axial kurtosis(AK) and radial kurtosis (RK)

Parameters

  • v1 (ndarray(dtype=float)) – Array of first eigenvectors from DWI.dtiTensorParams()

  • dt (ndarray(dtype=float)) – Array of diffusion tensor

Returns

  • rk (ndarray(dtype=float)) – Radial Kurtosis

  • ak (ndarray(dtype=float)) – Axial Kurtosis

  • kfa (ndarray(dtype=float)) – Kurtosis Fractional Anisotropy

  • mkt (ndarray(dtype=float)) – Mean Kurtosis Tensor

Examples

(rk, ak, kfa, mkt) = dwi.dkiTensorParams(v1, dt)

dtiTensorParams(nn)

Computes sorted DTI tensor eigenvalues and eigenvectors

Parameters

DT (ndarray(dtype=float)) – Diffusion tensor array

Returns

  • values (ndarray(dtype=float)) – Array of sorted eigenvalues

  • vectors (ndarray(dtype=float)) – Array pf sorted eigenvectors

Examples

(values, vectors) = dwi.dtiTensorParams(DT)

extractDKI()

Extract all DKI parameters from DT tensor. Warning, this can only be run after tensor fitting dwi.fit()

Returns

  • mk (ndarray(dtype=float)) – Mean Diffusivity

  • rk (ndarray(dtype=float)) – Radial Diffusivity

  • ak (ndarray(dtype=float)) – Axial Diffusivity

  • kfa (ndarray(dtype=float)) – Kurtosis Fractional Anisotropy

  • mkt (ndarray(dtype=float)) – Mean Kurtosis Tensor

  • trace (ndarray(dtype=float)) – Sum of first eigenvalues

Examples

(mk, rk, ak, fe, trace) = dwi.extractDTI(), where dwi is the DWI class object

extractDTI()

Extract all DTI parameters from DT tensor. Warning, this can only be run after tensor fitting dwi.fit()

Returns

  • md (ndarray(dtype=float)) – Mean Diffusivity

  • rd (ndarray(dtype=float)) – Radial Diffusivity

  • ad (ndarray(dtype=float)) – Axial Diffusivity

  • fa (ndarray(dtype=float)) – Fractional Anisotropy

  • fe (ndarray(dtype=float)) – First Eigenvectors

  • trace (ndarray(dtype=float)) – Sum of first eigenvalues

Examples

(md, rd, ad, fa) = dwi.extractDTI(), where dwi is the DWI class object

extractWMTI()

Returns white matter tract integrity (WMTI) parameters. Warning: this can only be run after fitting and DWI.extractDTI().

Returns

  • awf (ndarray(dtype=float)) – Axonal Water Fraction

  • eas_ad (ndarray(dtype=float)) – Extra-axonal space Axial Diffusivity

  • eas_rd (ndarray(dtype=float)) – Extra-axonal Space Radial Diffusivity

  • eas_md (ndarray(dtype=float)) – Extra-axonal Space Mean Diffusivity

  • eas_tort (ndarray(dtype=float)) – Extra-axonal Space Tortuosity

  • ias_ad (ndarray(dtype=float)) – Intra-axonal Space Axial Diffusivity

  • ias_rd (ndarray(dtype=float)) – Intra-axonal Space Radial Diffusivity

  • ias_da (ndarray(dtype=float)) – Intra-axonal Space Intrinsic Diffusivity

  • ias_tort (ndarray(dtype=float)) – Intra-axonal Space Tortuosity

fbi(l_max=6, fbwm=True, rectify=True, res='med')

Perform fiber ball imaging (FBI) and FBI white matter model (FBWM) analyses

Parameters

  • l_max (int) – Maximum spherical harmonic degree specified as an even integer (Default: 6)

  • fbwm (bool) – Perform FBWM parameterization if True (Default: True)

  • rectify (bool) – Perform fODF rectification if True (Default: True)

  • res (str) – Resolution of spherical sampling distribution (Default: ‘med’) ‘low’ defines the spherical grid defined by 3 fold quadrisection of the isocahedron, or 8 fold tesselation of icosahedron. ‘med’ defines the spherical grid defined by 4 fold quadrisection of the isocahedron, or 16 fold tesselation of icosahedron. ‘high’ defines the spherical grid defined by 5 fold quadrisection of the isocahedron, or 24 fold tesselation of icosahedron. Default: “med”

Returns

  • zeta (array_like(dtype=float)) – Zeta parameter

  • faa (array_like(dtype=float)) – Intra-axonal fractional anisotropy

  • fodf (array_like(dtype=float)) – fodf from spherical harmonic expansion

  • min_awf (array_like(dtype=float)) – Axonal water fraction

  • Da (array_like(dtype=float)) – Intrinsic intra-axonal diffusivity

  • De_mean (array_like(dtype=float)) – Mean extra-axonal diffusion

  • De_ax (array_like(dtype=float)) – Axial extra-axonal diffusion

  • De_rad (array_like(dtype=float)) – Radial extra-axonal diffusion

  • De_fa (array_like(dtype=float)) – Extra-axonal FA

  • min_cost (array_like(dtype=float)) – Minimum cost of the cost function (first index of min_cost_fn)

  • min_cost_fn (array_like(dtype=float)) – Cost function

fibonacciSphere(samples=1, randomize=True)

Returns evenly spaced points on a sphere

Parameters

  • samples (int) – Number of points to compute from sphere, must be a positive and real integer (Default: 1)

  • randomize (bool) – True if sampling is randomized; False otherwise (Default: True)

Returns

points – [3 x samples] array containing evenly spaced points from a sphere

Return type

ndarray(dtype=float)

Examples

dirs = dwi.fibonacciSphere(256, True)

findViols(c=[0, 1, 0])

Returns a 3D violation map of voxels that violate constraints.

Parameters

  • img (ndarray(dtype=float)) – 3D metric array such as mk or fa

  • c (array_like(dtype=int)) – [3 x 1] vector that toggles which constraints to check c[0]: Check D < 0 constraint c[1]: Check K < 0 constraint (Default) c[2]: Check K > 3/(b*D) constraint

Returns

map – 3D array containing locations of voxels that incur directional violations. Voxels with values contain violaions and voxel values represent proportion of directional violations.

Return type

ndarray(dtype=bool)

Examples

map = findViols(img, [0 1 0]

findVoxelViol(adcVox, akcVox, maxB, c)

Returns the proportions of violations occurring at a voxel.

Parameters

  • img (ndarray(dtype=float)) – 3D metric array such as mk or fa

  • c (array_like(dtype=int)) – [3 x 1] vector that toggles which constraints to check c[0]: Check D < 0 constraint c[1]: Check K < 0 constraint (Default) c[2]: Check K > 3/(b*D) constraint

Returns

n – percentaghhe ranging from 0 to 1 that indicates proportion of violations occuring at voxel.

Return type

ndarray(dtype=float)

Examples

map = findViols(voxel, [0 1 0]

fit(constraints=None, reject=None)

Returns fitted diffusion or kurtosis tensor

Parameters

  • constraints (array_like(dtype=int)) – [1 x 3] vector that specifies which constraints to use (Default: None)

  • reject (ndarray(dtype=bool)) – 4D array containing information on voxels to exclude from DT estimation (Default: None)

Examples

dwi.fit() dwi.fit(constraints=[0,1,0], reject=irlls_output)

getBvals()

Returns a vector of b-values, requires no input arguments

Returns

Vector array of b-values

Return type

ndarray(dtype=float)

Examples

bvals = dwi.getBvals(), where dwi is the DWI class object

getBvecs()

Returns an array of gradient vectors, requires no input parameters

Returns

[N x 3] array of gradient vectors

Return type

ndarray(dtype=float)

Examples

bvecs = dwi.getBvecs(), where dwi is the DWI class object

getndirs()

Returns the number of gradient directions acquired from the scanner

Returns

number of gradient directions

Return type

ndarray

Examples

n = dwi.getndirs(), where dwi is the DWI class object

goodDirections(outliers)

Creates a 3D maps of good directions from IRLLS outlier map For any given voxel, a violation is computed using logical or operant for all b-values. Whether an outlier occurs at b1000 or b1000 and b2000, that voxel is still a violation unless none of the b-values have outliers.

Parameters

outliers (ndarray(dtype=bool)) – 4D maps of outliers from IRLLS

Returns

map – 3D map of number of good directions

Return type

ndarray(dtype=int)

Examples

map = dwi.goodDirections(outliers)

idxb0()

Returns the index of all B-zeros according to bvals in record

Returns

idx – Index of DTI/DKI b-values

Return type

bool

idxdki()

Returns the index of all DTI/DKI B-values according to bvals in record

Returns

idx – Index of DTI/DKI b-values

Return type

bool

idxdti()

Returns the index of all DTI/DKI B-values according to bvals in record

Returns

idx – Index of DTI/DKI b-values

Return type

bool

idxfbi()

Returns the index of all FBI B-values according to bvals in record

Returns

Index of DTI/DKI b-values.

Return type

bool

irlls(excludeb0=True, maxiter=25, convcrit=0.001, mode='DKI', leverage=0.85, bounds=3)

This functions performs outlier detection and robust parameter estimation for diffusion MRI using the iterative reweigthed linear least squares (IRLLS) approach.

Parameters

  • exludeb0 (bool, optional) – Exlude the b0 images when removing outliers (Default: True)

  • maxiter (int, optional Integer; default: 25) – Maximum number of iterations in the iterative reweighting loop (Default: 25)

  • convcrit (float, optional) – Fraction of L2-norm of estimated diffusion parameter vector that the L2-norm of different vector should get in order to reach convergence in the iterative reweighted loop (Default: 1e-3)

  • mode (str, {‘DKI’, ‘DTI’}, optional) – Specifies whether to use DTI or DKI model (Default: ‘DKI’)

  • leverage (float, optional) – Measurement ranging from 0 to 1 where a leverage above this threshold will not be excluded in estimation of DT after outlier (Default: 0.85)

  • Bounds (int, optional) – Set the threshold of the number of standard deviation that are needed to exclude a measurement (Default: 3)

Returns

  • outliers (ndarray(dtype=bool)) – 4D image same size as input DWI marking voxels that are outliers

  • dt (ndarray(dtype=float)) – IRLLS method of DT estimation

Examples

outliers = dwi.irlls()

irllsviolmask(reject)

Computes 3D violation mask of outliers detected from IRLLS method

Parameters

reject (ndarray(dtype=bool)) – 4D input outlier map from IRLLS

Returns

propviol – 3D mask where voxel value is the percentage of directional violations

Return type

ndarray(dtype=float)

Examples

mask = dwi.irllsviolmask(outliers)

isdki()

Returns logical value to answer the mystical question whether the input image is DKI

Returns

ans – True if DKI; false otherwise

Return type

bool

Examples

ans = dwi.isdki(), where dwi is the DWI class object

isdti()

Returns logical value to answer the mystical question whether the input image is DTI

Returns

ans – True if DTI; false otherwise

Return type

bool

Examples

ans = dwi.isdki(), where dwi is the DWI class object

isfbi()

Returns bool value to specify whether image input image is FBI

Returns

and – True if FBI; false otherwise

Return type

bool

Examples

ans = dwi.isfbi(), where dwi is the DWI class object

isfbwm()

Returns bool value to specify whether image input image is FBWM

Returns

and – True if FBWM; false otherwise

Return type

bool

Examples

ans = dwi.isfbi(), where dwi is the DWI class object

kurtosisCoeff(dt, dir)

Computes apparent kurtosis coefficient (AKC)

Parameters

  • dt (ndarray(dtype=float)) – [21 x nvoxel] array containing diffusion tensor

  • dir (ndarray(dtype=float)) – [n x 3] array containing gradient directions

Returns

adc – Array containing apparent kurtosis coefficient

Return type

ndarray(dtype=float)

Examples

adc = dwi.kurtosisCoeff(dt, dir)

maxBval()

Returns the maximum b-value in a dataset to determine between DTI and DKI, requires no input parameters

Returns

maximum B-value in DWI

Return type

float

Examples

a = dwi.maxBval(), where dwi is the DWI class object

maxDKIBval()

Returns the maximum DKI b-value in a dataset

Returns

maximum DKI B-value in DWI

Return type

float

Examples

a = dwi.maxDKIBval(), where dwi is the DWI class object

maxDTIBval()

Returns the maximum DTI b-value in a dataset

Returns

maximum DTI B-value in DWI

Return type

float

Examples

a = dwi.maxDKIBval(), where dwi is the DWI class object

maxFBIBval()

Returns the maximum FBI b-value in a dataset

Returns

maximum DKI B-value in DWI

Return type

float

Examples

a = dwi.maxDKIBval(), where dwi is the DWI class object

multiplyMask(img)

Multiplies a 3D image by the brain mask

Parameters

img (ndarray(dtype=float)) – 3D image to be multiplied

Returns

multiplied image

Return type

ndarray(dtype=float)

optimal_lmax()

Computes the highest harmonic order (l_max) for spherical harmonic expansion. This is adapted from the information posted at https://mrtrix.readthedocs.io/en/dev/concepts/sh_basis_lmax.html

This function runs successfully only if input DWI is an FBI or HARDI acquisition.

Returns

l_max suitable for DWI

Return type

int

parfindViols(c=[0, 0, 0])
radialSampling(dir, n)

Get the radial component of a metric from a set of directions

Parameters

  • dir (ndarray(dtype=float)) – [n x 3] input array of directions

  • n (int) – number of rows, n

Returns

dirs

Return type

Matrix containing radial components

Examples

grad = dwi.radialSampling(dir, number_of_dirs)

tensorReorder(dwiType)

Reorders tensors in DT to those of MRTRIX in accordance to the table below

Parameters

dwiType (str, {‘dti’, ‘dki’}) – Indicates whether image is DTI or DKI

Returns

  • DT (ndarray(dtype=float)) – 4D image containing DT tensor

  • KT (ndarray(dtype=float)) – 4D image containing KT tensor

Examples

dt = dwi.tensorReorder()

Notes

MRTRIX3 and Designer tensors are described below.

MRTRIX3 Tensors                     DESIGNER Tensors
~~~~~~~~~~~~~~~                     ~~~~~~~~~~~~~~~~

0   D0      1   1                       1   1
1   D1      2   2                       1   2
2   D2      3   3                       1   3
3   D3      1   2                       2   2
4   D4      1   3                       2   3
5   D5      2   3                       3   3

6   K0      1   1   1   1               1   1   1   1
7   K1      2   2   2   2               1   1   1   2
8   K2      3   3   3   3               1   1   1   3
9   K3      1   1   1   2               1   1   2   2
10  K4      1   1   1   3               1   1   2   3
11  K5      1   2   2   2               1   1   3   3
12  K6      1   3   3   3               1   2   2   2
13  K7      2   2   2   3               1   2   2   3
14  K8      2   3   3   3               1   2   3   3
15  K9      1   1   2   2               1   3   3   3
16  K10     1   1   3   3               2   2   2   2
17  K11     2   2   3   3               2   2   2   3
18  K12     1   1   2   3               2   2   3   3
19  K13     1   2   2   3               2   3   3   3
20  K14     1   2   3   3               3   3   3   3

Value Assignment
~~~~~~~~~~~~~~~~

MRTRIX3         DESIGNER
~~~~~~~         ~~~~~~~~
    0               0
    1               3
    2               5
    3               1
    4               2
    5               4

    6               6
    7               16
    8               20
    9               7
    10              8
    11              12
    12              15
    13              17
    14              19
    15              9
    16              11
    17              18
    18              10
    19              13
    20              14
tensorType()

Returns whether input image is DTI or DKI compatible, requires no input parameters

Returns

contains list of string ‘dti’, ‘dki’, or ‘fbi’ based on the protocols the input DWI represents

Return type

list of str

Examples

a = dwi.tensorType(), where dwi is the DWI class object

wlls(shat, dwi, b, cons=None, warmup=None)

Estimates diffusion and kurtosis tenor at voxel with unconstrained Moore-Penrose pseudoinverse or constrained quadratic convex optimization. This is a helper function for dwi.fit() so a multiprocessing parallel loop can be iterated over voxels

Parameters

  • shat (ndarray(dtype=float)) – [ndir x 1] array of S_hat, approximated signal intensity at voxel

  • dwi (ndarray(dtype=float)) – [ndir x 1] array of diffusion weighted intensity values at voxel, for all b-values

  • b (ndarray(dtype=float)) – [ndir x 1] array of b-values vector

  • cons (ndarray(dtype=float)) – [(n * dir) x 22) array containing inequality constraints for fitting (Default: None)

  • warmup (ndarray(dtype=float)) – Estimated dt vector (22, 0) at each voxel for warm starting constrianed tensor fitting (Default: None)

Returns

dt – Diffusion tensor

Return type

ndarray(dtype=float)

Examples

dt = dwi.wlls(shat, dwi, b, constraints)

Notes

For Unconstrained Fitting: In the absence of constraints, an exact formulation in the form Cx = b is produced. This is further simplified to x_hat = C^+ * b. One can use the Moore-Penrose method to compute the pseudoinverse to approximate diffusion tensors.

For Constrained Fitting: .. code-block:: none

The equation |Cx -b|^2 expands to 0.5*x.T(C.T*A)*x -(C.T*b).T
~~~~~ ~~~~~

P q

where A is denoted by multiplier matrix (w * b) Multiplying by a positive constant (0.5) does not change the value of optimum x*. Similarly, the constant offset b.T*b does not affect x*, therefore we can leave these out.

Minimize: || C*x -b ||_2^2

subject to A*x <= b No lower or upper bounds

designer.fitting.dwipy.fit_regime(input, output, prefix=None, suffix=None, ext=None, irlls=True, akc=True, qcpath=None, fit_constraints=[0, 1, 0], l_max=None, rectify=True, res='med', n_fibers=5, mask=None, nthreads=None)

Performs the entire tensor fitting regime and writes out maps. Uses auto-detections methods to determine the types of protocols encoded by DWI and extract their metrics.

Parameters

  • input (str) – Path to DWI

  • output (str) – Output directory to write all outputs

  • prefix (str) – Prefix to append to output file names. (Default: None)

  • suffix (str) – Suffix to append to output file names. (Default: None)

  • ext (str) – Specify output image extension type. (Default: None)

  • irlls (bool) – Specify whether to perform IRLLS outlier detection. (Default: True)

  • akc (bool) – Specify whether to perform brute-forced AKC correction. (Default: True)

  • qcpath (str) – Specify output directory to write QC metrics. (Default: None)

  • fit_constraints (list) – List of 3 bool elements specifying which fit contraints to use. See DWI.createConstraints() on usage. (Default: [0, 1, 0])

  • l_max (int) – Maximum spherical harminic degree for FBI/FBWM fit. (Default: None)

  • rectify (bool) – Specify whether to rectify FBI fODF (Default: True)

  • mask (str) – Path to brain mask (Default: None)

  • nthreads (int) – Number of workers to use in processing. Default value uses all available workers. (Default: None)