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 [10]:
import numpy as np

from dipy.sims.phantom import orbital_phantom, add_rician_noise
from dipy.viz import fvtk
In [11]:
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 [12]:
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 [13]:
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 [14]:
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 [15]:
r=fvtk.ren()
fvtk.add(r, fvtk.volume(vol[..., 0]))
fvtk.show(r)
In [16]:
vol.shape
Out[16]:
(33, 33, 33, 65)
In [ ]: