Phantom Simulation¶
With a given fiber configuration, simulate the signal measured in S-space.
The volume, V, generated has shape (M, N, P, B), where B is the number of b-vectors.
Each measurement V[x, y, z, ...] is the integral of the diffusion probability density
function along the B different measurement directions. V[x, y, z, 0] is the reference
measurement, $S_0$, where no gradient is applied.
The signal is generated by tracing the fiber contour, determining its gradient at each position, and then adding the response of a single fiber along that direction. The response of a single fibre is given by
$$\Delta S = S_0 \exp(-b \mathbf{b}^\mathrm{T} R \Lambda R^\mathrm{T} \mathbf{b} ) $$where $S_0$ is the reference signal of the fibre, $\Lambda$ the diffusion tensor, and $R$ the rotation matrix to align the diffusion tensor with the fiber direction. For white matter, we typically take
$$\Lambda = \left(\begin{array}{ccc} \lambda_{\parallel} & 0 & 0\\\\ 0 & \lambda_{\perp} & 0\\\\ 0 & 0 & \lambda_{\perp} \end{array}\right) = 10^{-3}\left(\begin{array}{ccc} 1.4 & 0 & 0\\\\ 0 & 0.35 & 0\\\\ 0 & 0 & 0.35 \end{array}\right) \mathrm{mm}^2 \cdot \mathrm{s}^{-1}.$$Units:
- b: $\mathrm{s} / \mathrm{mm}^2$
- $\Lambda$: $\mathrm{mm}^2 / \mathrm{s}$
- Note that, with the above units, the exponent is dimensionless.
- $1\ \mathrm{T} = \frac{\mathrm{N} \cdot \mathrm{s}}{\mathrm{C} \cdot \mathrm{m}} = \frac{\mathrm{kg}}{\mathrm{C}\mathrm{s}}$
Also:
- $b = \gamma^2G \delta^2(\Delta - \delta / 3)$
- $\gamma = \frac{e}{2 m_p} g = 4257\ \mathrm{Hz}/\mathrm{gauss}$ : Gyromagnetic ratio of the proton
- $g \approx 5.586$: g-factor of the nucleus (dimensionless number)
- $m_p$: Mass of the proton
- $G$: field gradient amplitude
- $\Delta$, $\delta$: time constants related to the pulse sequence
import numpy as np
from dipy.sims.phantom import orbital_phantom, add_rician_noise
from dipy.viz import fvtk
def f1(t):
t = np.linspace(-1, 1, len(t)) * 0.95
x = t
y = 0 * t
z = t ** 2 + 0.05
return x, y, z
def f2(t):
t = np.linspace(-1, 1, len(t)) * 0.95
#x = (1 - np.cos(t)) ** 2
x = t ** 2 - 0.05
y = t
z = 0 * x
return x, y, z
def f3(t):
t = np.linspace(-1, 1, len(t)) * 0.95
x = -(t ** 2 + 0.05)
y = t
z = 0 * x
return x, y, z
t = np.linspace(0, 2 * np.pi, 100)
data_shape = (33, 33, 33, 65)
origin = np.array([16, 16, 16])
scale = 0.9 * origin
vol = orbital_phantom(func=f1, t=t, datashape=data_shape, origin=origin, scale=scale)
vol += orbital_phantom(func=f2, t=t, datashape=data_shape, origin=origin, scale=scale)
vol += orbital_phantom(func=f3, t=t, datashape=data_shape, origin=origin, scale=scale)
r=fvtk.ren()
fvtk.add(r, fvtk.volume(vol[..., 0]))
fvtk.show(r)
vol.shape