/usr/lib/python2.7/dist-packages/dipy/sims/tests/test_phantom.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 | from __future__ import division
import numpy as np
import nose
import nibabel as nib
import numpy.testing.decorators as dec
from numpy.testing import (assert_, assert_equal, assert_array_equal,
assert_array_almost_equal, assert_almost_equal,
run_module_suite)
from dipy.core.geometry import vec2vec_rotmat
from dipy.data import get_data
from dipy.reconst.dti import TensorModel
from dipy.sims.phantom import orbital_phantom, add_noise
from dipy.sims.voxel import single_tensor
from dipy.core.gradients import gradient_table
fimg,fbvals,fbvecs=get_data('small_64D')
bvals=np.load(fbvals)
bvecs=np.load(fbvecs)
bvecs[np.isnan(bvecs)]=0
gtab = gradient_table(bvals, bvecs)
def f(t):
"""
Helper function used to define a mapping time => xyz
"""
x = np.linspace(-1,1,len(t))
y = np.linspace(-1,1,len(t))
z = np.linspace(-1,1,len(t))
return x,y,z
def test_phantom():
N = 50
vol = orbital_phantom(gtab,
func=f,
t=np.linspace(0, 2 * np.pi, N),
datashape=(10, 10, 10, len(bvals)),
origin=(5, 5, 5),
scale=(3, 3, 3),
angles=np.linspace(0, 2 * np.pi, 16),
radii=np.linspace(0.2, 2, 6),
S0=100)
m = TensorModel(gtab)
t = m.fit(vol)
FA = t.fa
# print vol
FA[np.isnan(FA)] = 0
# 686 -> expected FA given diffusivities of [1500, 400, 400]
l1, l2, l3 = 1500e-6, 400e-6, 400e-6
expected_fa = (np.sqrt(0.5) * np.sqrt((l1 - l2)**2 + (l2-l3)**2 + (l3-l1)**2 )/np.sqrt(l1**2 + l2**2 + l3**2))
assert_array_almost_equal(FA.max(), expected_fa, decimal=2)
def test_add_noise():
np.random.seed(1980)
N = 50
S0 = 100
options = dict(func=f,
t=np.linspace(0, 2 * np.pi, N),
datashape=(10, 10, 10, len(bvals)),
origin=(5, 5, 5),
scale=(3, 3, 3),
angles=np.linspace(0, 2 * np.pi, 16),
radii=np.linspace(0.2, 2, 6),
S0=S0)
vol = orbital_phantom(gtab, **options)
for snr in [10, 20, 30, 50]:
vol_noise = orbital_phantom(gtab, snr=snr, **options)
sigma = S0 / snr
assert_(np.abs(np.var(vol_noise - vol) - sigma ** 2) < 1)
if __name__ == "__main__":
run_module_suite()
|