/usr/share/doc/dipy/examples/simulate_dki.py is in python-dipy 0.10.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 | """
==========================
DKI MultiTensor Simulation
==========================
In this example we show how to simulate the diffusion kurtosis imaging (DKI)
data of a single voxel. DKI captures information about the non-Gaussian
properties of water diffusion which is a consequence of the existence of tissue
barriers and compartments. In these simulations compartmental heterogeneity is
taken into account by modeling different compartments for the intra- and
extra-cellular media of two populations of fibers. These simulations are
performed according to [RNH2015]_.
We first import all relevant modules.
"""
import numpy as np
import matplotlib.pyplot as plt
from dipy.sims.voxel import (multi_tensor_dki, single_tensor)
from dipy.data import get_data
from dipy.io.gradients import read_bvals_bvecs
from dipy.core.gradients import gradient_table
from dipy.reconst.dti import (decompose_tensor, from_lower_triangular)
"""
For the simulation we will need a GradientTable with the b-values and
b-vectors. Here we use the GradientTable of the sample Dipy dataset
'small_64D'.
"""
fimg, fbvals, fbvecs = get_data('small_64D')
bvals, bvecs = read_bvals_bvecs(fbvals, fbvecs)
"""
DKI requires data from more than one non-zero b-value. Since the dataset
'small_64D' was acquired with one non-zero bvalue we artificialy produce a
second non-zero b-value.
"""
bvals = np.concatenate((bvals, bvals * 2), axis=0)
bvecs = np.concatenate((bvecs, bvecs), axis=0)
"""
The b-values and gradient directions are then converted to Dipy's
GradientTable format.
"""
gtab = gradient_table(bvals, bvecs)
"""
In ``mevals`` we save the eigenvalues of each tensor. To simulate crossing
fibers with two different media (representing intra and extra-cellular media),
a total of four components have to be taken in to account (i.e. the first two
compartments correspond to the intra and extra cellular media for the first
fiber population while the others correspond to the media of the second fiber
population)
"""
mevals = np.array([[0.00099, 0, 0],
[0.00226, 0.00087, 0.00087],
[0.00099, 0, 0],
[0.00226, 0.00087, 0.00087]])
"""
In ``angles`` we save in polar coordinates (:math:`\theta, \phi`) the principal
axis of each compartment tensor. To simulate crossing fibers at 70 degrees
the compartments of the first fiber are aligned to the x-axis while the
compartments of the second fiber are aligned to the x-z plane with an angular
deviation of 70 degrees from the first one.
"""
angles = [(90, 0), (90, 0), (20, 0), (20, 0)]
"""
In ``fractions`` we save the percentage of the contribution of each
compartment, which is computed by multiplying the percentage of contribution
of each fiber population and the water fraction of each different medium
"""
fie = 0.49 # intra axonal water fraction
fractions = [fie*50, (1 - fie)*50, fie*50, (1 - fie)*50]
"""
Having defined the parameters for all tissue compartments, the elements of the
diffusion tensor (dt), the elements of the kurtosis tensor (kt) and the DW
signals simulated from the DKI model can be obtain using the function
``multi_tensor_dki``.
"""
signal_dki, dt, kt = multi_tensor_dki(gtab, mevals, S0=200, angles=angles,
fractions=fractions, snr=None)
"""
We can also add rician noise with a specific SNR.
"""
signal_noisy, dt, kt = multi_tensor_dki(gtab, mevals, S0=200,
angles=angles, fractions=fractions,
snr=10)
"""
For comparison purposes, we also compute the DW signal if only the diffusion
tensor components are taken into account. For this we use Dipy's function
single_tensor which requires that dt is decomposed into its eigenvalues and
eigenvectors.
"""
dt_evals, dt_evecs = decompose_tensor(from_lower_triangular(dt))
signal_dti = single_tensor(gtab, S0=200, evals=dt_evals, evecs=dt_evecs,
snr=None)
"""
Finally, we can visualize the values of the different version of simulated
signals for all assumed gradient directions and bvalues.
"""
plt.plot(signal_dti, label='noiseless dti')
plt.plot(signal_dki, label='noiseless dki')
plt.plot(signal_noisy, label='with noise')
plt.legend()
plt.show()
plt.savefig('simulated_dki_signal.png')
"""
.. figure:: simulated_dki_signal.png
:align: center
**Simulated signals obtain from the DTI and DKI models**.
Non-Gaussian diffusion properties in tissues are responsible to smaller signal
attenuations for larger bvalues when compared to signal attenuations from free
gaussian water diffusion. This feature can be shown from the figure above,
since signals simulated from the DKI models reveals larger DW signal
intensities than the signals obtained only from the diffusion tensor
components.
References:
[RNH2015] R. Neto Henriques et al., "Exploring the 3D geometry of the diffusion
kurtosis tensor - Impact on the development of robust tractography
procedures and novel biomarkers", NeuroImage (2015) 111, 85-99.
"""
|