import numpy as np
from numpy import dot
from scipy.special import jn
from dipy.core.geometry import sphere2cart, vec2vec_rotmat
from dipy.reconst.utils import dki_design_matrix
from dipy.testing.decorators import warning_for_keywords
# Diffusion coefficients for white matter tracts, in mm^2/s
#
# Based roughly on values from:
#
# Pierpaoli, Basser, "Towards a Quantitative Assessment of Diffusion
# Anisotropy", Magnetic Resonance in Medicine, 1996; 36(6):893-906.
#
diffusion_evals = np.array([1500e-6, 400e-6, 400e-6])
def _check_directions(angles):
"""
Helper function to check if direction ground truth have the right format
and are in cartesian coordinates
Parameters
----------
angles : array (K,2) or (K, 3)
List of K polar angles (in degrees) for the sticks or array of K
sticks as unit vectors.
Returns
-------
sticks : (K,3)
Sticks in cartesian coordinates.
"""
angles = np.array(angles)
if angles.shape[-1] == 3:
sticks = angles
else:
sticks = [
sphere2cart(1, np.deg2rad(pair[0]), np.deg2rad(pair[1])) for pair in angles
]
sticks = np.array(sticks)
return sticks
def _add_gaussian(sig, noise1, noise2):
"""
Helper function to add_noise
This one simply adds one of the Gaussians to the sig and ignores the other
one.
"""
return sig + noise1
def _add_rician(sig, noise1, noise2):
"""
Helper function to add_noise.
This does the same as abs(sig + complex(noise1, noise2))
"""
return np.sqrt((sig + noise1) ** 2 + noise2**2)
def _add_rayleigh(sig, noise1, noise2):
r"""Helper function to add_noise.
The Rayleigh distribution is $\sqrt\{Gauss_1^2 + Gauss_2^2}$.
"""
return sig + np.sqrt(noise1**2 + noise2**2)
[docs]
@warning_for_keywords()
def add_noise(signal, snr, S0, *, noise_type="rician", rng=None):
r"""Add noise of specified distribution to the signal from a single voxel.
Parameters
----------
signal : 1-d ndarray
The signal in the voxel.
snr : float
The desired signal-to-noise ratio. (See notes below.)
If `snr` is None, return the signal as-is.
S0 : float
Reference signal for specifying `snr`.
noise_type : string, optional
The distribution of noise added. Can be either 'gaussian' for Gaussian
distributed noise, 'rician' for Rice-distributed noise (default) or
'rayleigh' for a Rayleigh distribution.
rng : numpy.random.Generator, optional
Random number generator for the noise. If ``None``, uses NumPy's default
random generator.
Returns
-------
signal : array, same shape as the input
Signal with added noise.
Notes
-----
SNR is defined here, following :footcite:p:`Descoteaux2007`, as
``S0 / sigma``, where ``sigma`` is the standard deviation of the two
Gaussian distributions forming the real and imaginary components of the
Rician noise distribution (see :footcite:p:`Gudbjartsson1995`).
References
----------
.. footbibliography::
Examples
--------
>>> signal = np.arange(800).reshape(2, 2, 2, 100)
>>> signal_w_noise = add_noise(signal, 10., 100., noise_type='rician')
"""
if snr is None:
return signal
if rng is None:
rng = np.random.default_rng()
sigma = S0 / snr
noise_adder = {
"gaussian": _add_gaussian,
"rician": _add_rician,
"rayleigh": _add_rayleigh,
}
noise1 = rng.normal(0, sigma, size=signal.shape)
if noise_type == "gaussian":
noise2 = None
else:
noise2 = rng.normal(0, sigma, size=signal.shape)
return noise_adder[noise_type](signal, noise1, noise2)
[docs]
@warning_for_keywords()
def sticks_and_ball(
gtab, *, d=0.0015, S0=1.0, angles=((0, 0), (90, 0)), fractions=(35, 35), snr=20
):
"""Simulate the signal for a Sticks & Ball model.
See :footcite:p:`Behrens2007` for a definition of the Sticks & Ball model.
Parameters
----------
gtab : GradientTable
Signal measurement directions.
d : float, optional
Diffusivity value.
S0 : float, optional
Unweighted signal value.
angles : array (K, 2) or (K, 3), optional
List of K polar angles (in degrees) for the sticks or array of K
sticks as unit vectors.
fractions : array-like, optional
Percentage of each stick. Remainder to 100 specifies isotropic
component.
snr : float, optional
Signal to noise ratio, assuming Rician noise. If set to None, no
noise is added.
Returns
-------
S : (N,) ndarray
Simulated signal.
sticks : (M,3)
Sticks in cartesian coordinates.
References
----------
.. footbibliography::
"""
fractions = [f / 100.0 for f in fractions]
f0 = 1 - np.sum(fractions)
S = np.zeros(len(gtab.bvals))
sticks = _check_directions(angles)
for i, g in enumerate(gtab.bvecs[1:]):
S[i + 1] = f0 * np.exp(-gtab.bvals[i + 1] * d) + np.sum(
[
fractions[j] * np.exp(-gtab.bvals[i + 1] * d * np.dot(s, g) ** 2)
for (j, s) in enumerate(sticks)
]
)
S[i + 1] = S0 * S[i + 1]
S[gtab.b0s_mask] = S0
S = add_noise(S, snr, S0)
return S, sticks
[docs]
def callaghan_perpendicular(q, radius):
"""Calculates the perpendicular diffusion signal E(q) in a cylinder of
radius R using the Soderman model.
Assumes that the pulse length is infinitely short and the diffusion time is
infinitely long.
See :footcite:p:`Soderman1995` for details about the Soderman model.
Parameters
----------
q : array, shape (N,)
q-space value in 1/mm
radius : float
cylinder radius in mm
Returns
-------
E : array, shape (N,)
signal attenuation
References
----------
.. footbibliography::
"""
# Eq. [6] in the paper
numerator = (2 * jn(1, 2 * np.pi * q * radius)) ** 2
denom = (2 * np.pi * q * radius) ** 2
E = np.divide(numerator, denom, out=np.zeros_like(q), where=denom != 0)
return E
[docs]
@warning_for_keywords()
def gaussian_parallel(q, tau, *, D=0.7e-3):
r"""Calculates the parallel Gaussian diffusion signal.
Parameters
----------
q : array, shape (N,)
q-space value in 1/mm
tau : float
diffusion time in s
D : float, optional
diffusion constant
Returns
-------
E : array, shape (N,)
signal attenuation
"""
return np.exp(-((2 * np.pi * q) ** 2) * tau * D)
[docs]
@warning_for_keywords()
def cylinders_and_ball_soderman(
gtab,
tau,
*,
radii=(5e-3, 5e-3),
D=0.7e-3,
S0=1.0,
angles=((0, 0), (90, 0)),
fractions=(35, 35),
snr=20,
):
"""Calculates the three-dimensional signal attenuation E(q) originating
from within a cylinder of radius R using the Soderman approximation.
The diffusion signal is assumed to be separable perpendicular and parallel
to the cylinder axis :footcite:p:`Assaf2004`.
This function is basically an extension of the ball and stick model.
Setting the radius to zero makes them equivalent.
See :footcite:p:`Soderman1995` for details about the Soderman model.
Parameters
----------
gtab : GradientTable
Signal measurement directions.
tau : float
diffusion time in s
radii : array-like, optional
cylinder radius in mm
D : float, optional
diffusion constant
S0 : float, optional
Unweighted signal value.
angles : array (K, 2) or (K, 3), optional
List of K polar angles (in degrees) for the sticks or array of K
sticks as unit vectors.
fractions : array-like, optional
Percentage of each stick. Remainder to 100 specifies isotropic
component.
snr : float, optional
Signal to noise ratio, assuming Rician noise. If set to None, no
noise is added.
Returns
-------
E : array, shape (N,)
signal attenuation
References
----------
.. footbibliography::
"""
qvals = np.sqrt(gtab.bvals / tau) / (2 * np.pi)
qvecs = qvals[:, None] * gtab.bvecs
q_norm = np.sqrt(np.einsum("ij,ij->i", qvecs, qvecs))
fractions = [f / 100.0 for f in fractions]
f0 = 1 - np.sum(fractions)
S = np.zeros(len(gtab.bvals))
sticks = _check_directions(angles)
for i, f in enumerate(fractions):
q_par = abs(np.dot(qvecs, sticks[i]))
q_perp = np.sqrt(q_norm**2 - q_par**2)
S_cylinder = callaghan_perpendicular(q_perp, radii[i]) * gaussian_parallel(
q_par, tau, D=D
)
S += f * S_cylinder
S += f0 * np.exp(-gtab.bvals * D)
S *= S0
S[gtab.b0s_mask] = S0
S = add_noise(S, snr, S0)
return S, sticks
[docs]
@warning_for_keywords()
def single_tensor(gtab, S0=1, *, evals=None, evecs=None, snr=None, rng=None):
"""Simulate diffusion-weighted signals with a single tensor.
See :footcite:p:`Stejskal1965`, :footcite:p:`Descoteaux2008b`.
Parameters
----------
gtab : GradientTable
Table with information of b-values and gradient directions g.
Note that if gtab has a btens attribute, simulations will be performed
according to the given b-tensor B information.
S0 : double, optional
Strength of signal in the presence of no diffusion gradient (also
called the ``b=0`` value).
evals : (3,) ndarray, optional
Eigenvalues of the diffusion tensor. By default, values typical for
prolate white matter are used.
evecs : (3, 3) ndarray, optional
Eigenvectors of the tensor. You can also think of this as a rotation
matrix that transforms the direction of the tensor. The eigenvectors
need to be column wise.
snr : float, optional
Signal to noise ratio, assuming Rician noise. None implies no noise.
rng : numpy.random.Generator, optional
Random number generator for the noise. If ``None``, uses NumPy's default
random generator.
Returns
-------
S : (N,) ndarray
Simulated signal:
``S(b, g) = S_0 e^(-b g^T R D R.T g)``, if gtab.tens=None
``S(B) = S_0 e^(-B:D)``, if gtab.tens information is given
References
----------
.. footbibliography::
"""
if rng is None:
rng = np.random.default_rng()
if evals is None:
evals = diffusion_evals
if evecs is None:
evecs = np.eye(3)
out_shape = gtab.bvecs.shape[: gtab.bvecs.ndim - 1]
gradients = gtab.bvecs.reshape(-1, 3)
R = np.asarray(evecs)
S = np.zeros(len(gradients))
D = dot(dot(R, np.diag(evals)), R.T)
if gtab.btens is None:
for i, g in enumerate(gradients):
S[i] = S0 * np.exp(-gtab.bvals[i] * dot(dot(g.T, D), g))
else:
for i, b in enumerate(gtab.btens):
S[i] = S0 * np.exp(-np.sum(b * D))
S = add_noise(S, snr, S0, rng=rng)
return S.reshape(out_shape)
[docs]
@warning_for_keywords()
def multi_tensor(
gtab,
mevals,
*,
S0=1.0,
angles=((0, 0), (90, 0)),
fractions=(50, 50),
snr=20,
rng=None,
):
r"""Simulate a Multi-Tensor signal.
Parameters
----------
gtab : GradientTable
Table with information of b-values and gradient directions.
Note that if gtab has a btens attribute, simulations will be performed
according to the given b-tensor information.
mevals : array (K, 3)
each tensor's eigenvalues in each row
S0 : float, optional
Unweighted signal value (b0 signal).
angles : array (K, 2) or (K, 3), optional
List of K tensor directions in polar angles (in degrees) or unit
vectors
fractions : array-like, optional
Percentage of the contribution of each tensor. The sum of fractions
should be equal to 100%.
snr : float, optional
Signal to noise ratio, assuming Rician noise. If set to None, no
noise is added.
rng : numpy.random.Generator, optional
Random number generator for the noise. If ``None``, uses NumPy's default
random generator.
Returns
-------
S : (N,) ndarray
Simulated signal.
sticks : (M,3)
Sticks in cartesian coordinates.
Examples
--------
>>> import numpy as np
>>> from dipy.sims.voxel import multi_tensor
>>> from dipy.data import get_fnames
>>> from dipy.core.gradients import gradient_table
>>> from dipy.io.gradients import read_bvals_bvecs
>>> fimg, fbvals, fbvecs = get_fnames(name='small_101D')
>>> bvals, bvecs = read_bvals_bvecs(fbvals, fbvecs)
>>> gtab = gradient_table(bvals, bvecs=bvecs)
>>> mevals=np.array(([0.0015, 0.0003, 0.0003],[0.0015, 0.0003, 0.0003]))
>>> e0 = np.array([1, 0, 0.])
>>> e1 = np.array([0., 1, 0])
>>> S = multi_tensor(gtab, mevals)
"""
if rng is None:
rng = np.random.default_rng()
if np.round(np.sum(fractions), 2) != 100.0:
raise ValueError("Fractions should sum to 100")
fractions = [f / 100.0 for f in fractions]
S = np.zeros(len(gtab.bvals))
sticks = _check_directions(angles)
for i in range(len(fractions)):
S = S + fractions[i] * single_tensor(
gtab, S0=S0, evals=mevals[i], evecs=all_tensor_evecs(sticks[i]), snr=None
)
return add_noise(S, snr, S0, rng=rng), sticks
[docs]
@warning_for_keywords()
def multi_tensor_dki(
gtab,
mevals,
*,
S0=1.0,
angles=((90.0, 0.0), (90.0, 0.0)),
fractions=(50, 50),
snr=20,
):
r"""Simulate the diffusion-weight signal, diffusion and kurtosis tensors
based on the DKI model
See :footcite:p:`NetoHenriques2015` for further details.
Parameters
----------
gtab : GradientTable
Gradient table.
mevals : array (K, 3)
eigenvalues of the diffusion tensor for each individual compartment
S0 : float, optional
Unweighted signal value (b0 signal).
angles : array (K,2) or (K,3), optional
List of K tensor directions of the diffusion tensor of each compartment
in polar angles (in degrees) or unit vectors
fractions : float (K,), optional
Percentage of the contribution of each tensor. The sum of fractions
should be equal to 100%.
snr : float, optional
Signal to noise ratio, assuming Rician noise. If set to None, no
noise is added.
Returns
-------
S : (N,) ndarray
Simulated signal based on the DKI model.
dt : (6,)
elements of the diffusion tensor.
kt : (15,)
elements of the kurtosis tensor.
Notes
-----
Simulations are based on multicompartmental models which assumes that
tissue is well described by impermeable diffusion compartments
characterized by their only diffusion tensor. Since simulations are based
on the DKI model, coefficients larger than the fourth order of the signal's
taylor expansion approximation are neglected.
Examples
--------
>>> import numpy as np
>>> from dipy.sims.voxel import multi_tensor_dki
>>> from dipy.data import get_fnames
>>> from dipy.core.gradients import gradient_table
>>> from dipy.io.gradients import read_bvals_bvecs
>>> fimg, fbvals, fbvecs = get_fnames(name='small_64D')
>>> bvals, bvecs = read_bvals_bvecs(fbvals, fbvecs)
>>> bvals_2s = np.concatenate((bvals, bvals * 2), axis=0)
>>> bvecs_2s = np.concatenate((bvecs, bvecs), axis=0)
>>> gtab = gradient_table(bvals_2s, bvecs=bvecs_2s)
>>> mevals = np.array([[0.00099, 0, 0],[0.00226, 0.00087, 0.00087]])
>>> S, dt, kt = multi_tensor_dki(gtab, mevals)
References
----------
.. footbibliography::
"""
if np.round(np.sum(fractions), 2) != 100.0:
raise ValueError("Fractions should sum to 100")
fractions = [f / 100.0 for f in fractions]
# S = np.zeros(len(gtab.bvals))
sticks = _check_directions(angles)
# computing a 3D matrix containing the individual DT components
D_comps = np.zeros((len(fractions), 3, 3))
for i in range(len(fractions)):
R = all_tensor_evecs(sticks[i])
D_comps[i] = dot(dot(R, np.diag(mevals[i])), R.T)
# compute voxel's DT
DT = np.zeros((3, 3))
for i in range(len(fractions)):
DT = DT + fractions[i] * D_comps[i]
dt = np.array([DT[0][0], DT[0][1], DT[1][1], DT[0][2], DT[1][2], DT[2][2]])
# compute voxel's MD
MD = (DT[0][0] + DT[1][1] + DT[2][2]) / 3
# compute voxel's KT
kt = np.zeros(15)
kt[0] = kurtosis_element(D_comps, fractions, 0, 0, 0, 0, DT=DT, MD=MD)
kt[1] = kurtosis_element(D_comps, fractions, 1, 1, 1, 1, DT=DT, MD=MD)
kt[2] = kurtosis_element(D_comps, fractions, 2, 2, 2, 2, DT=DT, MD=MD)
kt[3] = kurtosis_element(D_comps, fractions, 0, 0, 0, 1, DT=DT, MD=MD)
kt[4] = kurtosis_element(D_comps, fractions, 0, 0, 0, 2, DT=DT, MD=MD)
kt[5] = kurtosis_element(D_comps, fractions, 0, 1, 1, 1, DT=DT, MD=MD)
kt[6] = kurtosis_element(D_comps, fractions, 1, 1, 1, 2, DT=DT, MD=MD)
kt[7] = kurtosis_element(D_comps, fractions, 0, 2, 2, 2, DT=DT, MD=MD)
kt[8] = kurtosis_element(D_comps, fractions, 1, 2, 2, 2, DT=DT, MD=MD)
kt[9] = kurtosis_element(D_comps, fractions, 0, 0, 1, 1, DT=DT, MD=MD)
kt[10] = kurtosis_element(D_comps, fractions, 0, 0, 2, 2, DT=DT, MD=MD)
kt[11] = kurtosis_element(D_comps, fractions, 1, 1, 2, 2, DT=DT, MD=MD)
kt[12] = kurtosis_element(D_comps, fractions, 0, 0, 1, 2, DT=DT, MD=MD)
kt[13] = kurtosis_element(D_comps, fractions, 0, 1, 1, 2, DT=DT, MD=MD)
kt[14] = kurtosis_element(D_comps, fractions, 0, 1, 2, 2, DT=DT, MD=MD)
# compute S based on the DT and KT
S = dki_signal(gtab, dt, kt, S0=S0, snr=snr)
return S, dt, kt
[docs]
@warning_for_keywords()
def kurtosis_element(D_comps, frac, ind_i, ind_j, ind_k, ind_l, *, DT=None, MD=None):
r"""Computes the diffusion kurtosis tensor element (with indexes i, j, k
and l) based on the individual diffusion tensor components of a
multicompartmental model.
Parameters
----------
D_comps : (K,3,3) ndarray
Diffusion tensors for all K individual compartment of the
multicompartmental model.
frac : [float]
Percentage of the contribution of each tensor. The sum of fractions
should be equal to 100%.
ind_i : int
Element's index i (0 for x, 1 for y, 2 for z)
ind_j : int
Element's index j (0 for x, 1 for y, 2 for z)
ind_k : int
Element's index k (0 for x, 1 for y, 2 for z)
ind_l: int
Elements index l (0 for x, 1 for y, 2 for z)
DT : (3,3) ndarray, optional
Voxel's global diffusion tensor.
MD : float, optional
Voxel's global mean diffusivity.
Returns
-------
wijkl : float
kurtosis tensor element of index i, j, k, l
Notes
-----
wijkl is calculated using equation 8 given in
:footcite:p:`NetoHenriques2015`.
References
----------
.. footbibliography::
"""
if DT is None:
DT = np.zeros((3, 3))
for i in range(len(frac)):
DT = DT + frac[i] * D_comps[i]
if MD is None:
MD = (DT[0][0] + DT[1][1] + DT[2][2]) / 3
wijkl = 0
for f in range(len(frac)):
wijkl = wijkl + frac[f] * (
D_comps[f][ind_i][ind_j] * D_comps[f][ind_k][ind_l]
+ D_comps[f][ind_i][ind_k] * D_comps[f][ind_j][ind_l]
+ D_comps[f][ind_i][ind_l] * D_comps[f][ind_j][ind_k]
)
wijkl = (
wijkl
- DT[ind_i][ind_j] * DT[ind_k][ind_l]
- DT[ind_i][ind_k] * DT[ind_j][ind_l]
- DT[ind_i][ind_l] * DT[ind_j][ind_k]
) / (MD**2)
return wijkl
[docs]
@warning_for_keywords()
def dki_signal(gtab, dt, kt, *, S0=150, snr=None):
r"""Simulated signal based on the diffusion and diffusion kurtosis
tensors of a single voxel. Simulations are performed assuming the DKI
model.
See :footcite:p:`NetoHenriques2015` for further details.
Parameters
----------
gtab : GradientTable
Measurement directions.
dt : (6,) ndarray
Elements of the diffusion tensor.
kt : (15, ) ndarray
Elements of the diffusion kurtosis tensor.
S0 : float, optional
Strength of signal in the presence of no diffusion gradient.
snr : float, optional
Signal to noise ratio, assuming Rician noise. None implies no noise.
Returns
-------
S : (N,) ndarray
Simulated signal based on the DKI model:
.. math::
S=S_{0}e^{-bD+\frac{1}{6}b^{2}D^{2}K}
References
----------
.. footbibliography::
"""
dt = np.array(dt)
kt = np.array(kt)
A = dki_design_matrix(gtab)
# define vector of DKI parameters
MD = (dt[0] + dt[2] + dt[5]) / 3
X = np.concatenate((dt, kt * MD * MD, np.array([-np.log(S0)])), axis=0)
# Compute signals based on the DKI model
S = np.exp(dot(A, X))
S = add_noise(S, snr, S0)
return S
[docs]
@warning_for_keywords()
def single_tensor_odf(r, *, evals=None, evecs=None):
"""Simulated ODF with a single tensor.
See :footcite:p:`Aganj2010` for further details.
Parameters
----------
r : (N,3) or (M,N,3) ndarray
Measurement positions in (x, y, z), either as a list or on a grid.
evals : (3,)
Eigenvalues of diffusion tensor. By default, use values typical for
prolate white matter.
evecs : (3, 3) ndarray
Eigenvectors of the tensor, written column-wise. You can also think
of these as the rotation matrix that determines the orientation of
the diffusion tensor.
Returns
-------
ODF : (N,) ndarray
The diffusion probability at ``r`` after time ``tau``.
References
----------
.. footbibliography::
"""
if evals is None:
evals = diffusion_evals
if evecs is None:
evecs = np.eye(3)
out_shape = r.shape[: r.ndim - 1]
R = np.asarray(evecs)
D = dot(dot(R, np.diag(evals)), R.T)
Di = np.linalg.inv(D)
r = r.reshape(-1, 3)
P = np.zeros(len(r))
for i, u in enumerate(r):
P[i] = (dot(dot(u.T, Di), u)) ** (3 / 2)
return (1 / (4 * np.pi * np.prod(evals) ** (1 / 2) * P)).reshape(out_shape)
[docs]
def all_tensor_evecs(e0):
"""Given the principle tensor axis, return the array of all
eigenvectors column-wise (or, the rotation matrix that orientates the
tensor).
Parameters
----------
e0 : (3,) ndarray
Principle tensor axis.
Returns
-------
evecs : (3,3) ndarray
Tensor eigenvectors, arranged column-wise.
"""
axes = np.eye(3)
mat = vec2vec_rotmat(axes[0], e0)
e1 = np.dot(mat, axes[1])
e2 = np.dot(mat, axes[2])
# Return the eigenvectors column-wise:
return np.array([e0, e1, e2]).T
[docs]
def multi_tensor_odf(odf_verts, mevals, angles, fractions):
"""Simulate a Multi-Tensor ODF.
Parameters
----------
odf_verts : (N,3) ndarray
Vertices of the reconstruction sphere.
mevals : sequence of 1D arrays,
Eigen-values for each tensor.
angles : sequence of 2d tuples,
Sequence of principal directions for each tensor in polar angles
or cartesian unit coordinates.
fractions : sequence of floats,
Percentages of the fractions for each tensor.
Returns
-------
ODF : (N,) ndarray
Orientation distribution function.
Examples
--------
Simulate a MultiTensor ODF with two peaks and calculate its exact ODF.
>>> import numpy as np
>>> from dipy.sims.voxel import multi_tensor_odf, all_tensor_evecs
>>> from dipy.data import default_sphere
>>> vertices, faces = default_sphere.vertices, default_sphere.faces
>>> mevals = np.array(([0.0015, 0.0003, 0.0003],[0.0015, 0.0003, 0.0003]))
>>> angles = [(0, 0), (90, 0)]
>>> odf = multi_tensor_odf(vertices, mevals, angles, [50, 50])
"""
mf = [f / 100.0 for f in fractions]
sticks = _check_directions(angles)
odf = np.zeros(len(odf_verts))
mevecs = []
for s in sticks:
mevecs += [all_tensor_evecs(s)]
for j, f in enumerate(mf):
odf += f * single_tensor_odf(odf_verts, evals=mevals[j], evecs=mevecs[j])
return odf
[docs]
@warning_for_keywords()
def single_tensor_rtop(*, evals=None, tau=1.0 / (4 * np.pi**2)):
"""Simulate a Single-Tensor rtop.
See :footcite:p:`Cheng2012` for further details.
Parameters
----------
evals : 1D arrays, optional
Eigen-values for the tensor. By default, values typical for prolate
white matter are used.
tau : float, optional
diffusion time. By default the value that makes q=sqrt(b).
Returns
-------
rtop : float,
Return to origin probability.
References
----------
.. footbibliography::
"""
if evals is None:
evals = diffusion_evals
rtop = 1.0 / np.sqrt((4 * np.pi * tau) ** 3 * np.prod(evals))
return rtop
[docs]
@warning_for_keywords()
def multi_tensor_rtop(mf, *, mevals=None, tau=1 / (4 * np.pi**2)):
"""Simulate a Multi-Tensor rtop.
See :footcite:p:`Cheng2012` for further details.
Parameters
----------
mf : sequence of floats, bounded [0,1]
Percentages of the fractions for each tensor.
mevals : sequence of 1D arrays, optional
Eigen-values for each tensor. By default, values typical for prolate
white matter are used.
tau : float, optional
diffusion time. By default the value that makes q=sqrt(b).
Returns
-------
rtop : float,
Return to origin probability.
References
----------
.. footbibliography::
"""
rtop = 0
if mevals is None:
mevals = [
None,
] * len(mf)
for j, f in enumerate(mf):
rtop += f * single_tensor_rtop(evals=mevals[j], tau=tau)
return rtop
[docs]
@warning_for_keywords()
def single_tensor_pdf(r, *, evals=None, evecs=None, tau=1 / (4 * np.pi**2)):
"""Simulated ODF with a single tensor.
See :footcite:p:`Cheng2012` for further details.
Parameters
----------
r : (N,3) or (M,N,3) ndarray
Measurement positions in (x, y, z), either as a list or on a grid.
evals : (3,), optional
Eigenvalues of diffusion tensor. By default, use values typical for
prolate white matter.
evecs : (3, 3) ndarray, optional
Eigenvectors of the tensor. You can also think of these as the
rotation matrix that determines the orientation of the diffusion
tensor.
tau : float, optional
diffusion time. By default the value that makes q=sqrt(b).
Returns
-------
pdf : (N,) ndarray
The diffusion probability at ``r`` after time ``tau``.
References
----------
.. footbibliography::
"""
if evals is None:
evals = diffusion_evals
if evecs is None:
evecs = np.eye(3)
out_shape = r.shape[: r.ndim - 1]
R = np.asarray(evecs)
D = dot(dot(R, np.diag(evals)), R.T)
Di = np.linalg.inv(D)
r = r.reshape(-1, 3)
P = np.zeros(len(r))
for i, u in enumerate(r):
P[i] = (-dot(dot(u.T, Di), u)) / (4 * tau)
pdf = (1 / np.sqrt((4 * np.pi * tau) ** 3 * np.prod(evals))) * np.exp(P)
return pdf.reshape(out_shape)
[docs]
@warning_for_keywords()
def multi_tensor_pdf(pdf_points, mevals, angles, fractions, *, tau=1 / (4 * np.pi**2)):
"""Simulate a Multi-Tensor ODF.
See :footcite:p:`Cheng2012` for further details.
Parameters
----------
pdf_points : (N, 3) ndarray
Points to evaluate the PDF.
mevals : sequence of 1D arrays,
Eigen-values for each tensor. By default, values typical for prolate
white matter are used.
angles : sequence,
Sequence of principal directions for each tensor in polar angles
or cartesian unit coordinates.
fractions : sequence of floats,
Percentages of the fractions for each tensor.
tau : float, optional
diffusion time. By default the value that makes q=sqrt(b).
Returns
-------
pdf : (N,) ndarray,
Probability density function of the water displacement.
References
----------
.. footbibliography::
"""
mf = [f / 100.0 for f in fractions]
sticks = _check_directions(angles)
pdf = np.zeros(len(pdf_points))
mevecs = []
for s in sticks:
mevecs += [all_tensor_evecs(s)]
for j, f in enumerate(mf):
pdf += f * single_tensor_pdf(
pdf_points, evals=mevals[j], evecs=mevecs[j], tau=tau
)
return pdf
[docs]
@warning_for_keywords()
def single_tensor_msd(*, evals=None, tau=1 / (4 * np.pi**2)):
"""Simulate a Multi-Tensor rtop.
See :footcite:p:`Cheng2012` for further details.
Parameters
----------
evals : 1D arrays, optional
Eigen-values for the tensor. By default, values typical for prolate
white matter are used.
tau : float, optional
diffusion time. By default the value that makes q=sqrt(b).
Returns
-------
msd : float,
Mean square displacement.
References
----------
.. footbibliography::
"""
if evals is None:
evals = diffusion_evals
msd = 2 * tau * np.sum(evals)
return msd
[docs]
@warning_for_keywords()
def multi_tensor_msd(mf, *, mevals=None, tau=1 / (4 * np.pi**2)):
"""Simulate a Multi-Tensor rtop.
See :footcite:p:`Cheng2012` for further details.
Parameters
----------
mf : sequence of floats, bounded [0,1]
Percentages of the fractions for each tensor.
mevals : sequence of 1D arrays, optional
Eigen-values for each tensor. By default, values typical for prolate
white matter are used.
tau : float, optional
diffusion time. By default the value that makes q=sqrt(b).
Returns
-------
msd : float,
Mean square displacement.
References
----------
.. footbibliography::
"""
msd = 0
if mevals is None:
mevals = [
None,
] * len(mf)
for j, f in enumerate(mf):
msd += f * single_tensor_msd(evals=mevals[j], tau=tau)
return msd