# 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 In :
import numpy as np

from dipy.viz import fvtk

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

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

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

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

In :
r=fvtk.ren()

vol.shape

(33, 33, 33, 65)