qsiprep.utils.shm module

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.

class qsiprep.utils.shm.SphHarmFit(model, shm_coef, mask)

Bases: OdfFit

Diffusion data fit to a spherical harmonic model

gfa()

odf(sphere)

Samples the odf function on the points of a sphere

Parameters:

sphere (Sphere) – The points on which to sample the odf.

Returns:

values – The value of the odf on each point of sphere.

Return type:

ndarray

predict(gtab=None, S0=1.0)

Predict the diffusion signal from the model coefficients.

Parameters:

• gtab (a GradientTable class instance) – The directions and bvalues on which prediction is desired

• S0 (float 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

class qsiprep.utils.shm.SphHarmModel(gtab)

Bases: OdfModel, Cache

To be subclassed by all models that return a SphHarmFit when fit.

sampling_matrix(sphere)

The matrix needed to sample ODFs from coefficients of the model.

Parameters:

sphere (Sphere) – Points used to sample ODF.

Returns:

sampling_matrix – 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.

Return type:

array

qsiprep.utils.shm.anisotropic_power(sh_coeffs, norm_factor=1e-05, power=2, non_negative=True)

Calculates anisotropic power map with a given SH coefficient matrix

Parameters:

• sh_coeffs (ndarray) – A ndarray where the last dimension is the SH coefficients estimates for that voxel.

• norm_factor (float, optional) – The value to normalize the ap values. Default is 10^-5.

• power (int, optional) – The degree to which power maps are calculated. Default: 2.

• non_negative (bool, optional) – Whether to rectify the resulting map to be non-negative. Default: True.

Returns:

log_ap – The log of the resulting power image.

Return type:

ndarray

Notes

Calculate AP image based on a IxJxKxC SH coefficient matrix based on the following equation.

\[AP = \sum_{l=2,4,6,...}{\frac{1}{2l+1} \sum_{m=-l}^l{|a_{l,m}|^n}}\]

Where the last dimension, C, is made of a flattened array of $l$x$m$ coefficients, where $l$ are the SH orders, and $m = 2l+1$, So l=1 has 1 coeffecient, l=2 has 5, … l=8 has 17 and so on. A l=2 SH coefficient matrix will then be composed of a IxJxKx6 volume. The power, $n$ is usually set to $n=2$.

The final AP image is then shifted by -log(norm_factor), to be strictly non-negative. Remaining values < 0 are discarded (set to 0), per default, and this option is controlled through the non_negative keyword argument.

References

[1]

Dell’Acqua, F., Lacerda, L., Catani, M., Simmons, A., 2014. Anisotropic Power Maps: A diffusion contrast to reveal low anisotropy tissues from HARDI data, in: Proceedings of International Society for Magnetic Resonance in Medicine. Milan, Italy.

qsiprep.utils.shm.bootstrap_data_array(data, H, R, permute=None)

Applies the Residual Bootstraps to the data given H and R

data must be normalized, ie 0 < data <= 1

This function, and the bootstrap_data_voxel function, calculate residual-bootsrap samples given a Hat matrix and a Residual matrix. These samples can be used for non-parametric statistics or for bootstrap probabilistic tractography:

References

[1]

J. I. Berman, et al., “Probabilistic streamline q-ball tractography using the residual bootstrap” 2008.

[2]

HA Haroon, et al., “Using the model-based residual bootstrap to quantify uncertainty in fiber orientations from Q-ball analysis” 2009.

[3]

B. Jeurissen, et al., “Probabilistic Fiber Tracking Using the Residual Bootstrap with Constrained Spherical Deconvolution” 2011.

qsiprep.utils.shm.bootstrap_data_voxel(data, H, R, permute=None)

Like bootstrap_data_array but faster when for a single voxel

data must be 1d and normalized

qsiprep.utils.shm.calculate_max_order(n_coeffs)

Calculate the maximal harmonic order, given that you know the

number of parameters that were estimated.

n_coeffsint

The number of SH coefficients

Lint

The maximal SH order, given the number of coefficients

The calculation in this function proceeds according to the following logic.

\[n =\]

rac{1}{2} (L+1) (L+2)

arrow 2n = L^2 + 3L + 2

arrow L^2 + 3L + 2 - 2n = 0

arrow L^2 + 3L + 2(1-n) = 0

arrow L_{1,2} = rac{-3 pm sqrt{9 - 8 (1-n)}}{2}

arrow L{1,2} = rac{-3 pm sqrt{1 + 8n}}{2}

Finally, the positive value is chosen between the two options.

qsiprep.utils.shm.forward_sdeconv_mat(r_rh, n)

Build forward spherical deconvolution matrix

Parameters:

• r_rh (ndarray) – Rotational harmonics coefficients for the single fiber response function. Each element rh[i] is associated with spherical harmonics of degree 2*i.

• n (ndarray) – The degree of spherical harmonic function associated with each row of the deconvolution matrix. Only even degrees are allowed

Returns:

R – Deconvolution matrix with shape (N, N)

Return type:

ndarray (N, N)

qsiprep.utils.shm.gen_dirac(m, n, theta, phi)

Generate Dirac delta function orientated in (theta, phi) on the sphere

The spherical harmonics (SH) representation of this Dirac is returned as coefficients to spherical harmonic functions produced by shm.real_sph_harm.

Parameters:

• m (ndarray (N,)) – The order of the spherical harmonic function associated with each coefficient.

• n (ndarray (N,)) – The degree of the spherical harmonic function associated with each coefficient.

• theta (float [0, 2*pi]) – The azimuthal (longitudinal) coordinate.

• phi (float [0, pi]) – The polar (colatitudinal) coordinate.

See also

shm.real_sph_harm, shm.real_sym_sh_basis

Returns:

dirac – SH coefficients representing the Dirac function. The shape of this is (m + 2) * (m + 1) / 2.

Return type:

ndarray

qsiprep.utils.shm.hat(B)

Returns the hat matrix for the design matrix B

qsiprep.utils.shm.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.

qsiprep.utils.shm.lcr_matrix(H)

Returns a matrix for computing leveraged, centered residuals from data

if r = (d-Hd), the leveraged centered residuals are lcr = (r/l)-mean(r/l) ruturns the matrix R, such lcr = Rd

qsiprep.utils.shm.normalize_data(data, where_b0, min_signal=1.0, out=None)

Normalizes the data with respect to the mean b0

qsiprep.utils.shm.order_from_ncoef(ncoef)

Given a number n of coefficients, calculate back the sh_order

qsiprep.utils.shm.randint(low, high=None, size=None, dtype=int)

Return random integers from low (inclusive) to high (exclusive).

Return random integers from the “discrete uniform” distribution of the specified dtype in the “half-open” interval [low, high). If high is None (the default), then results are from [0, low).

Note

New code should use the ~numpy.random.Generator.integers method of a ~numpy.random.Generator instance instead; please see the Quick start.

Parameters:

• low (int or array-like of ints) – Lowest (signed) integers to be drawn from the distribution (unless high=None, in which case this parameter is one above the highest such integer).

• high (int or array-like of ints, optional) – If provided, one above the largest (signed) integer to be drawn from the distribution (see above for behavior if high=None). If array-like, must contain integer values

• size (int or tuple of ints, optional) – Output shape. If the given shape is, e.g., (m, n, k), then m * n * k samples are drawn. Default is None, in which case a single value is returned.

• dtype (dtype, optional) – Desired dtype of the result. Byteorder must be native. The default value is int.

  Added in version 1.11.0.

Returns:

out – size-shaped array of random integers from the appropriate distribution, or a single such random int if size not provided.

Return type:

int or ndarray of ints

See also

random_integers

similar to randint, only for the closed interval [low, high], and 1 is the lowest value if high is omitted.

random.Generator.integers

which should be used for new code.

Examples

>>> np.random.randint(2, size=10)
array([1, 0, 0, 0, 1, 1, 0, 0, 1, 0]) # random
>>> np.random.randint(1, size=10)
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

Generate a 2 x 4 array of ints between 0 and 4, inclusive:

>>> np.random.randint(5, size=(2, 4))
array([[4, 0, 2, 1], # random
       [3, 2, 2, 0]])

Generate a 1 x 3 array with 3 different upper bounds

>>> np.random.randint(1, [3, 5, 10])
array([2, 2, 9]) # random

Generate a 1 by 3 array with 3 different lower bounds

>>> np.random.randint([1, 5, 7], 10)
array([9, 8, 7]) # random

Generate a 2 by 4 array using broadcasting with dtype of uint8

>>> np.random.randint([1, 3, 5, 7], [[10], [20]], dtype=np.uint8)
array([[ 8,  6,  9,  7], # random
       [ 1, 16,  9, 12]], dtype=uint8)

qsiprep.utils.shm.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:

• m (int |m| <= n) – The order of the harmonic.

• n (int >= 0) – The degree of the harmonic.

• theta (float [0, 2*pi]) – The azimuthal (longitudinal) coordinate.

• phi (float [0, pi]) – The polar (colatitudinal) coordinate.

Returns:

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

Return type:

real float

See also

scipy.special.sph_harm

qsiprep.utils.shm.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 fibernavigator [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_order (int) – even int > 0, max spherical harmonic degree

• theta (float [0, 2*pi]) – The azimuthal (longitudinal) coordinate.

• phi (float [0, pi]) – The polar (colatitudinal) coordinate.

Returns:

• y_mn (real float) – The real harmonic $Y^m_n$ sampled at theta and phi

• m (array) – The order of the harmonics.

• n (array) – The degree of the harmonics.

References

[1]

https://github.com/scilus/fibernavigator

qsiprep.utils.shm.real_sym_sh_brainsuite(sh_order, theta, phi)

Compute the real spherical harmonics used in BrainSuite [1]

Parameters:

• sh_order (int) – The maximum degree of the spherical harmonic basis.

• theta (float [0, 2*pi]) – The azimuthal (longitudinal) coordinate.

• phi (float [0, pi]) – The polar (colatitudinal) coordinate.

Returns:

• y_mn (real float) – The real harmonic $Y^m_n$ sampled at theta and phi as implemented in BrainSuite.

• m (array) – The order of the harmonics.

• n (array) – The degree of the harmonics.

References

..[1] Justin P. Haldar, Richard M. Leahy, “Linear transforms for Fourier

data on the sphere: Application to high angular resolution diffusion MRI of the brain”, NeuroImage, 2013.

qsiprep.utils.shm.real_sym_sh_mrtrix(sh_order, theta, phi)

Compute real spherical harmonics as in mrtrix, where the real harmonic $Y^m_n$ is defined to be:

Real($Y^m_n$)       if m > 0
$Y^0_n$             if m = 0
Imag($Y^|m|_n$)     if m < 0

This may take scalar or array arguments. The inputs will be broadcasted against each other.

Parameters:

• sh_order (int) – The maximum degree or the spherical harmonic basis.

• theta (float [0, pi]) – The polar (colatitudinal) coordinate.

• phi (float [0, 2*pi]) – The azimuthal (longitudinal) coordinate.

Returns:

• y_mn (real float) – The real harmonic $Y^m_n$ sampled at theta and phi as implemented in mrtrix. Warning: the basis is Tournier et al 2004 and 2007 is slightly different.

• m (array) – The order of the harmonics.

• n (array) – The degree of the harmonics.

qsiprep.utils.shm.sf_to_sh(sf, sphere, sh_order=4, basis_type=None, smooth=0.0)

Spherical function to spherical harmonics (SH).

Parameters:

• sf (ndarray) – Values of a function on the given sphere.

• sphere (Sphere) – The points on which the sf is defined.

• sh_order (int, optional) – Maximum SH order in the SH fit. For sh_order, there will be (sh_order + 1) * (sh_order_2) / 2 SH coefficients (default 4).

• basis_type ({None, ‘mrtrix’, ‘fibernav’}) – None for the default dipy basis, mrtrix for the MRtrix basis, and fibernav for the FiberNavigator basis (default None).

• smooth (float, optional) – Lambda-regularization in the SH fit (default 0.0).

Returns:

sh – SH coefficients representing the input function.

Return type:

ndarray

qsiprep.utils.shm.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_sh (ndarray (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.

• m (ndarray (N,)) – The order of the spherical harmonic function associated with each coefficient.

• n (ndarray (N,)) – The degree of the spherical harmonic function associated with each coefficient.

Returns:

r_rh – Rotational harmonics coefficients representing the input r_sh

Return type:

ndarray ((sh_order + 1)*(sh_order + 2)/2,)

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

qsiprep.utils.shm.sh_to_sf(sh, sphere, sh_order, basis_type=None)

Spherical harmonics (SH) to spherical function (SF).

Parameters:

• sh (ndarray) – SH coefficients representing a spherical function.

• sphere (Sphere) – The points on which to sample the spherical function.

• sh_order (int, optional) – Maximum SH order in the SH fit. For sh_order, there will be (sh_order + 1) * (sh_order_2) / 2 SH coefficients (default 4).

• basis_type ({None, ‘mrtrix’, ‘fibernav’}) – None for the default dipy basis, mrtrix for the MRtrix basis, and fibernav for the FiberNavigator basis (default None).

Returns:

sf – Spherical function values on the sphere.

Return type:

ndarray

qsiprep.utils.shm.sh_to_sf_matrix(sphere, sh_order, basis_type=None, return_inv=True, smooth=0)

Matrix that transforms Spherical harmonics (SH) to spherical function (SF).

Parameters:

• sphere (Sphere) – The points on which to sample the spherical function.

• sh_order (int, optional) – Maximum SH order in the SH fit. For sh_order, there will be (sh_order + 1) * (sh_order_2) / 2 SH coefficients (default 4).

• basis_type ({None, ‘mrtrix’, ‘fibernav’}) – None for the default dipy basis, mrtrix for the MRtrix basis, and fibernav for the FiberNavigator basis (default None).

• return_inv (bool) – If True then the inverse of the matrix is also returned

• smooth (float, optional) – Lambda-regularization in the SH fit (default 0.0).

Returns:

• B (ndarray) – Matrix that transforms spherical harmonics to spherical function sf = np.dot(sh, B).

• invB (ndarray) – Inverse of B.

qsiprep.utils.shm.smooth_pinv(B, L)

Regularized pseudo-inverse

Computes a regularized least square inverse of B

Parameters:

• B (array_like (n, m)) – Matrix to be inverted

• L (array_like (n,))

Returns:

inv – regularized least square inverse of B

Return type:

ndarray (m, n)

Notes

In the literature this inverse is often written $(B^{T}B+L^{2})^{-1}B^{T}$. However here this inverse is implemented using the pseudo-inverse because it is more numerically stable than the direct implementation of the matrix product.

qsiprep.utils.shm.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_order (int) – even int > 0, max degree to return

Returns:

• m_list (array) – orders of even spherical harmonics

• n_list (array) – degrees of even spherical harmonics

See also

real_sph_harm

qsiprep.utils.shm.spherical_harmonics(m, n, theta, phi)

Compute spherical harmonics

This may take scalar or array arguments. The inputs will be broadcasted against each other.

Parameters:

• m (int |m| <= n) – The order of the harmonic.

• n (int >= 0) – The degree of the harmonic.

• theta (float [0, 2*pi]) – The azimuthal (longitudinal) coordinate.

• phi (float [0, pi]) – The polar (colatitudinal) coordinate.

Returns:

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

Return type:

complex float

Notes

This is a faster implementation of scipy.special.sph_harm for scipy version < 0.15.0. For scipy 0.15 and onwards, we use the scipy implementation of the function