/usr/lib/python2.7/dist-packages/dipy/viz/tests/test_fvtk_utils.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 | import numpy as np
import numpy.testing as npt
from dipy.viz.utils import map_coordinates_3d_4d
def trilinear_interp_numpy(input_array, indices):
""" Evaluate the input_array data at the given indices
"""
if input_array.ndim <= 2 or input_array.ndim >= 5:
raise ValueError("Input array can only be 3d or 4d")
x_indices = indices[:, 0]
y_indices = indices[:, 1]
z_indices = indices[:, 2]
x0 = x_indices.astype(np.integer)
y0 = y_indices.astype(np.integer)
z0 = z_indices.astype(np.integer)
x1 = x0 + 1
y1 = y0 + 1
z1 = z0 + 1
# Check if xyz1 is beyond array boundary:
x1[np.where(x1 == input_array.shape[0])] = x0.max()
y1[np.where(y1 == input_array.shape[1])] = y0.max()
z1[np.where(z1 == input_array.shape[2])] = z0.max()
if input_array.ndim == 3:
x = x_indices - x0
y = y_indices - y0
z = z_indices - z0
elif input_array.ndim == 4:
x = np.expand_dims(x_indices - x0, axis=1)
y = np.expand_dims(y_indices - y0, axis=1)
z = np.expand_dims(z_indices - z0, axis=1)
output = (input_array[x0, y0, z0] * (1 - x) * (1 - y) * (1 - z) +
input_array[x1, y0, z0] * x * (1 - y) * (1 - z) +
input_array[x0, y1, z0] * (1 - x) * y * (1-z) +
input_array[x0, y0, z1] * (1 - x) * (1 - y) * z +
input_array[x1, y0, z1] * x * (1 - y) * z +
input_array[x0, y1, z1] * (1 - x) * y * z +
input_array[x1, y1, z0] * x * y * (1 - z) +
input_array[x1, y1, z1] * x * y * z)
return output
def test_trilinear_interp():
A = np.zeros((5, 5, 5))
A[2, 2, 2] = 1
indices = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2], [1.5, 1.5, 1.5]])
values = trilinear_interp_numpy(A, indices)
values2 = map_coordinates_3d_4d(A, indices)
npt.assert_almost_equal(values, values2)
B = np.zeros((5, 5, 5, 3))
B[2, 2, 2] = np.array([1, 1, 1])
values = trilinear_interp_numpy(B, indices)
values_4d = map_coordinates_3d_4d(B, indices)
npt.assert_almost_equal(values, values_4d)
if __name__ == '__main__':
npt.run_module_suite()
|