Commit 6fc1186a authored by inhuszar's avatar inhuszar
Browse files

Pre-release v2.1

parent fa6d0484
......@@ -14,8 +14,8 @@ module_finv2d = \
"tirl/cmodules/finv3d_src.c",
"tirl/cmodules/finv.pyx"],
include_dirs=[np.get_include(),
os.path.join(os.environ.get("CONDA_PREFIX"),
"include")],
os.path.join(os.environ.get("CONDA_PREFIX"), "include"),
os.path.join(os.environ.get("CONDA_PREFIX"), "include", "gsl")],
extra_link_args=["-lcblas", "-lblas", "-llapacke"],
extra_compile_args=["-std=c99"],
language="c")
......@@ -80,6 +80,5 @@ setup(name="tirl",
"tirl/transformations/linear",
"tirl/transformations/nonlinear",
"tirl/usermodules"],
scripts=["tirl/tirl",
"tirl/tirlvision/sliceview.py"]
scripts=["tirl/tirl"]
)
#!/usr/bin/env python
import numpy as np
import nibabel as nib
from time import time
import scipy.sparse as sp
from scipy.signal.windows import gaussian
......@@ -12,7 +11,7 @@ import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
from tirl.cmodules.finv import invert2d, invert3d
from tirl.cmodules.finv import local_affine2d, local_affine3d
from tirl.cmodules.finv import local_affine
from tirl.cmodules.finv import jacobian2d, jacobian3d
from tirl.transformations.nonlinear.displacement import TxDisplacementField
......@@ -242,8 +241,9 @@ def jac2d():
h, w = shape = (256, 256)
f = generate_field(shape, type="sine")
locations = np.stack(np.mgrid[0:h, 0:w], axis=-1).reshape((-1, 2))
jac = local_affine2d(f, locations, output=None, pretx=None,
vectorder=(0, 1), maxiter=1000)
gridlocations = locations
jac = local_affine(f, locations, gridlocations, output=None,
vectorder=(0, 1), maxiter=1000)
jac = jac.reshape((-1, 2, 3), order="F")[..., :2]
#jac = np.pad(jac, pad_width=((0, 0), (0, 1), (0, 0)), constant_values=1)
#jacdet = np.linalg.det(jac[:, :2, :2])
......@@ -257,11 +257,13 @@ def jac2d():
def jac3d():
h, w, d = shape = (64, 64, 64)
f = generate_field(shape, type="gaussian")
f = generate_field(shape, type="sine")
locations = np.stack(np.mgrid[0:h, 0:w, 0:d], axis=-1).reshape((-1, 3))
jac = local_affine3d(f, locations, output=None, pretx=None,
vectorder=(0, 1, 2), maxiter=1000)
jac = jac.reshape((-1, 3, 4), order="F")[..., :3]
gridlocations = locations
jac = local_affine(f, locations, gridlocations, output=None,
vectorder=(0, 1, 2), maxiter=1000)
jac = jac.reshape(-1, 4, 3).transpose(0, 2, 1)[..., :3]
#jac = jac.reshape((-1, 3, 4), order="F")[..., :3]
jacdet = np.linalg.det(jac)
# jacdet[~np.isfinite(jacdet)] = 0
print(f"Jacobian range: {np.min(jacdet)} - {np.max(jacdet)}")
......
No preview for this file type
......@@ -5,14 +5,13 @@ from tirl.domain import Domain
from tirl import settings as ts
from tirl.transformations.basic.direct import TxDirect
from tirl.transformations.linear.scale import TxScale
from tirl.transformations.linear.rotation import TxRotation2D
from tirl.transformations.linear.translation import TxTranslation
from tirl.transformations.linear.shear import TxShear
# DEFINITIONS
ENABLE_HDD_TESTS = True
EPS = np.finfo(ts.DEFAULT_FLOAT_TYPE).eps
# IMPLEMENTATION
......@@ -38,7 +37,8 @@ class TestVoxelCoordinates(unittest.TestCase):
def setUp(self) -> None:
self.compact = Domain((200, 300))
self.noncompact = Domain(np.random.rand(16, 2))
self.rand = np.random.rand(16, 2)
self.noncompact = Domain(self.rand)
def test_get_all_voxel_coordinates(self):
# Compact
......@@ -51,9 +51,9 @@ class TestVoxelCoordinates(unittest.TestCase):
# Non-compact
voxels = self.noncompact.get_voxel_coordinates(dtype=np.int16)
self.assertTupleEqual(voxels.shape, (16, 2))
self.assertTrue(np.allclose(voxels[0] - [0, 0], 0))
self.assertTrue(np.allclose(voxels[-1] - [15, 0], 0))
self.assertEqual(voxels.dtype, np.int16)
self.assertTrue(np.allclose(voxels[0] - self.rand[0], 0, atol=EPS))
self.assertTrue(np.allclose(voxels[-1] - self.rand[15], 0, atol=EPS))
self.assertEqual(voxels.dtype, ts.DEFAULT_FLOAT_TYPE)
def test_get_linear_chunk_voxel_coordinates(self):
# Compact
......@@ -67,9 +67,9 @@ class TestVoxelCoordinates(unittest.TestCase):
# Non-compact
voxels = self.noncompact.get_voxel_coordinates(chunk=(8, 12))
self.assertTupleEqual(voxels.shape, (4, 2))
self.assertTrue(np.allclose(voxels[0] - [8, 0], 0))
self.assertTrue(np.allclose(voxels[-1] - [11, 0], 0))
self.assertEqual(voxels.dtype, np.int8)
self.assertTrue(np.allclose(voxels[0] - self.rand[8], 0, atol=EPS))
self.assertTrue(np.allclose(voxels[-1] - self.rand[11], 0, atol=EPS))
self.assertEqual(voxels.dtype, ts.DEFAULT_FLOAT_TYPE)
def test_get_multiindex_voxel_coordinates(self):
# Compact
......@@ -83,9 +83,9 @@ class TestVoxelCoordinates(unittest.TestCase):
# Non-compact
voxels = self.noncompact.get_voxel_coordinates(chunk=slice(8, 12))
self.assertTupleEqual(voxels.shape, (4, 2))
self.assertTrue(np.allclose(voxels[0] - [8, 0], 0))
self.assertTrue(np.allclose(voxels[-1] - [11, 0], 0))
self.assertEqual(voxels.dtype, np.int8)
self.assertTrue(np.allclose(voxels[0] - self.rand[8], 0, atol=EPS))
self.assertTrue(np.allclose(voxels[-1] - self.rand[11], 0, atol=EPS))
self.assertEqual(voxels.dtype, ts.DEFAULT_FLOAT_TYPE)
class TestSubchunkVoxelCoordinates(unittest.TestCase):
......@@ -103,7 +103,7 @@ class TestSubchunkVoxelCoordinates(unittest.TestCase):
self.assertEqual(vcoords.dtype, self.dtype)
def test_linear_chunk_coordinates_mem(self):
vcoords = self.domain.calculate_voxel_coordinates_mem(
vcoords = self.domain._calculate_voxel_coordinates_mem(
dtype="u2", chunk=(17000, 19000))
self.assertTupleEqual(vcoords.shape, (2000, 2))
self.assertTrue(np.allclose(vcoords[0] - [0, 17000], 0))
......@@ -113,7 +113,7 @@ class TestSubchunkVoxelCoordinates(unittest.TestCase):
@unittest.skipUnless(ENABLE_HDD_TESTS, "HDD tests disabled")
def test_linear_chunk_coordinates_hdd(self):
memlimit = 1 * 1024 ** 2 # 1 MB
vcoords = self.domain.calculate_voxel_coordinates_hdd(
vcoords = self.domain._calculate_voxel_coordinates_hdd(
dtype="u8", chunk=(17000, 90000), memlimit=memlimit)
self.assertTupleEqual(vcoords.shape, (73000, 2))
self.assertTrue(np.allclose(vcoords[0] - [0, 17000], 0))
......@@ -121,7 +121,7 @@ class TestSubchunkVoxelCoordinates(unittest.TestCase):
self.assertEqual(vcoords.dtype, "u8")
def test_multiindex_chunk_coordinates_mem(self):
vcoords = self.domain.calculate_voxel_coordinates_mem(
vcoords = self.domain._calculate_voxel_coordinates_mem(
dtype="u2", chunk=(slice(5000, 6000, 2), slice(7000, 10000, 3)))
self.assertTupleEqual(vcoords.shape, (5e5, 2))
self.assertTrue(np.allclose(vcoords[0] - [5000, 7000], 0))
......@@ -134,7 +134,7 @@ class TestSubchunkVoxelCoordinates(unittest.TestCase):
@unittest.skipUnless(ENABLE_HDD_TESTS, "HDD tests disabled")
def test_multiindex_chunk_coordinates_hdd(self):
memlimit = 1 * 1024 ** 2 # 1 MB
vcoords = self.domain.calculate_voxel_coordinates_hdd(
vcoords = self.domain._calculate_voxel_coordinates_hdd(
dtype="u2", chunk=(slice(5000, 6000, 2), slice(7000, 10000, 3)),
memlimit=memlimit)
self.assertTupleEqual(vcoords.shape, (5e5, 2))
......
......@@ -2,7 +2,7 @@ import os
import unittest
import numpy as np
import tirl
import tirl.fsl
from tirl.timage import TImage
from tirl.chain import Chain
from tirl.transformations.nonlinear.displacement import TxDisplacementField
......
......@@ -2,10 +2,11 @@ import unittest
import numpy as np
from pydoc import locate
import matplotlib.pyplot as plt
from scipy.interpolate import Rbf
# import matplotlib.pyplot as plt
# from scipy.interpolate import Rbf
from scipy.ndimage.interpolation import map_coordinates
from tirl.constants import *
import tirl.settings as ts
from tirl.domain import Domain
from tirl.tfield import TField
......@@ -15,6 +16,7 @@ from tirl.transformations.nonlinear.rbf_displacement import \
TxRbfDisplacementField
@unittest.skip
class TestConstruction2D(unittest.TestCase):
def setUp(self) -> None:
......@@ -84,6 +86,7 @@ def create_field():
return field, domain, points
@unittest.skip
class TestMapping2D(unittest.TestCase):
def setUp(self) -> None:
......@@ -110,6 +113,7 @@ class TestMapping2D(unittest.TestCase):
self.assertTrue(np.allclose(y, y0, 1e-3))
@unittest.skip
class TestInverse2D(unittest.TestCase):
def setUp(self) -> None:
......@@ -157,6 +161,7 @@ class TestInverse2D(unittest.TestCase):
self.assertTrue(np.isclose(np.median(x - xp), 0, atol=1e-1))
@unittest.skip
class TestUnidirectionalInverse2D3D(unittest.TestCase):
def setUp(self) -> None:
......@@ -204,5 +209,49 @@ class TestUnidirectionalInverse2D3D(unittest.TestCase):
print(x - xp)
class TestMapVector(unittest.TestCase):
def setUp(self):
import tirl
warpfile = tirl.home("..", "tests", "resources", "tirlfiles",
"rbfwarp.tx")
self.warp = tirl.load(warpfile)
def test_map_vector_rule_finite_strain(self):
# Define a horizontal voxel vector at the middle of the warp domain,
# which should coincide with the middle of the image being warped.
voxvect = [0, 1]
c = self.warp.domain.pcentre()
# To compute the mapped vector, the vector should be first mapped by the
# preceding transformations of the domain that the warp is defined on.
vect = self.warp.domain.map_voxel_vectors(voxvect, rule=RULE_FS)
self.assertEqual(vect.size, 3)
# Using RULE_FS, the size of the vectors should remain the same
self.assertAlmostEqual(np.linalg.norm(voxvect), np.linalg.norm(vect), 4)
warpedvect = self.warp.map_vector(vect, c, rule=RULE_FS)
self.assertAlmostEqual(
np.linalg.norm(vect), np.linalg.norm(warpedvect), 4)
# Using RULE_FS, the size of the vectors should remain the same
self.assertEqual(warpedvect.size, 3)
def test_map_vector_rule_ssr(self):
# Define a horizontal voxel vector at the middle of the warp domain,
# which should coincide with the middle of the image being warped.
voxvect = [0, 1]
c = self.warp.domain.pcentre()
# To compute the mapped vector, the vector should be first mapped by the
# preceding transformations of the domain that the warp is defined on.
vect = self.warp.domain.map_voxel_vectors(voxvect, rule=RULE_SSR)
self.assertEqual(vect.size, 3)
# Using RULE_SSR, the size of the vectors change will have changed
self.assertNotAlmostEqual(
np.linalg.norm(voxvect), np.linalg.norm(vect), 4)
warpedvect = self.warp.map_vector(vect, c, rule=RULE_SSR)
self.assertNotAlmostEqual(
np.linalg.norm(vect), np.linalg.norm(warpedvect), 4)
# Using RULE_SSR, the size of the vectors change will have changed
self.assertEqual(warpedvect.size, 3)
if __name__ == '__main__':
unittest.main()
......@@ -2,7 +2,7 @@ import unittest
import numpy as np
import tirl
from tirl import beta, BetaClassMeta
from tirl.beta import beta_function, BetaClassMeta
from tirl.timage import TImage
from tirl.tirlobject import TIRLObject
......@@ -17,7 +17,7 @@ class TestImage(unittest.TestCase):
"Test Image must not be all-zeros.")
@beta
@beta_function
def func_that_is_beta():
return True
......@@ -30,7 +30,7 @@ class BetaClass(TIRLObject, metaclass=BetaClassMeta):
class TestBetaMarkers(unittest.TestCase):
@beta
@beta_function
def method_that_is_beta(self):
return True
......
......@@ -214,18 +214,13 @@ own optimisation strategies.
# DEPENDENCIES
import os
import warnings
from functools import wraps
from tirl.tirlobject import InstanceCounterMeta
# PACKAGE IMPORTS
from tirl import fsl
from tirl import settings as ts
if ts.ENABLE_VISUALISATION:
from tirl import tirlvision
# from tirl import fsl
# from tirl import settings as ts
# if ts.ENABLE_VISUALISATION:
# from tirl import tirlvision
# DEFINITIONS
......@@ -251,6 +246,7 @@ def home(*args):
will be returned.
"""
import os
homedir = os.path.dirname(__file__)
if args:
return os.path.join(homedir, *args)
......@@ -261,7 +257,7 @@ def home(*args):
def testimg():
""" Returns a test TImage. """
from tirl.timage import TImage
impath = os.path.join(home(), "../tests/resources/testimage/testimg.png")
impath = home(*"../tests/resources/testimage/testimg.png".split("/"))
img = TImage.fromfile(impath, tensor_axes=(2,))
return img
......@@ -278,43 +274,3 @@ def load(fname):
"""
from tirl import tirlobject
return tirlobject.TIRLObject.load(fname)
def beta(func):
"""
Decorator that indicates when a certain function/method is untested. Warns
the user that results from the respective function/method might be
incorrect.
:param func: untested function or method
:type func: Callable
:returns: untested function or method with a warning message
:rtype: Callable
"""
@wraps(func)
def wrapped(*args, **kwargs):
warnings.warn(
f"{func.__name__} is a beta function in TIRL {__version__}. "
f"Use with caution: results might be incorrect.")
return func(*args, **kwargs)
return wrapped
class BetaClassMeta(InstanceCounterMeta):
"""
Metaclass that indicates when a certain class is untested. Warns the user
that results from the respective class might be incorrect.
"""
def __new__(cls, name, bases, attrs):
initfunc = attrs['__init__']
@wraps(initfunc)
def wrapped(*args, **kwargs):
self = args[0]
warnings.warn(f"{self.__class__.__name__} is a beta class in TIRL "
f"{__version__}. Use with caution: results might be incorrect.")
initfunc(*args, **kwargs)
attrs['__init__'] = wrapped
return super(BetaClassMeta, cls).__new__(cls, name, bases, attrs)
#!/usr/bin/env fslpython
# -*- coding: utf-8 -*-
# _______ _____ _____ _
# |__ __| |_ _| | __ \ | |
# | | | | | |__) | | |
# | | | | | _ / | |
# | | _| |_ | | \ \ | |____
# |_| |_____| |_| \_\ |______|
#
# Copyright (C) 2018-2020 University of Oxford
# Part of the FMRIB Software Library (FSL)
# Author: Istvan N. Huszar
# SHBASECOPYRIGHT
# DESCRIPTION
""" Beta function and beta class indicators. """
# DEPENDENCIES
import tirl
from tirl.tirlobject import InstanceCounterMeta
# DEFINITIONS
__version__ = tirl.__version__
# IMPLEMENTATION
def beta_function(func):
"""
Decorator that indicates when a certain function/method is untested. Warns
the user that results from the respective function/method might be
incorrect.
:param func: untested function or method
:type func: Callable
:returns: untested function or method with a warning message
:rtype: Callable
"""
import warnings
from functools import wraps
@wraps(func)
def wrapped(*args, **kwargs):
warnings.warn(
f"{func.__name__} is a beta function in TIRL {__version__}. "
f"Use with caution: results might be incorrect.")
return func(*args, **kwargs)
return wrapped
class BetaClassMeta(InstanceCounterMeta):
"""
Metaclass that indicates when a certain class is untested. Warns the user
that results from the respective class might be incorrect.
"""
def __new__(cls, name, bases, attrs):
import warnings
from functools import wraps
initfunc = attrs['__init__']
@wraps(initfunc)
def wrapped(*args, **kwargs):
self = args[0]
warnings.warn(f"{self.__class__.__name__} is a beta class in TIRL "
f"{__version__}. Use with caution: results might be incorrect.")
initfunc(*args, **kwargs)
attrs['__init__'] = wrapped
return super(BetaClassMeta, cls).__new__(cls, name, bases, attrs)
......@@ -107,7 +107,8 @@ class Cache(object):
# Identify caller function/method
if caller is None:
caller = str(inspect.stack()[1][0].f_code.co_name)
# caller = str(inspect.stack()[1][0].f_code.co_name)
caller = "default"
try:
fcache = self.cache[caller]
......
This diff is collapsed.
......@@ -172,90 +172,6 @@ def invert2d(field, output=None, vectorder=(0, 1), pretx=None, posttx=None,
invwarp.reshape((dim, *backward_shape)).ravel(), dtype=output.dtype)
def local_affine2d(field, locations, output=None, pretx=None,
vectorder=(0, 1), maxiter=1000):
# Obtain the properties of the input field (this must be absolute)
dim, *forward_shape = field.shape
cdef np.ndarray[int, ndim=1, mode="c"] fwd_shape = \
np.asarray(forward_shape, dtype=ctypes.c_int)
cdef int* fsh = <int *>malloc(2 * sizeof(int))
fsh = &fwd_shape[0]
# assert dim == 2, "Input is not a 2D vector field."
n_fwd_points = np.product(fwd_shape)
# Obtain locations
n_points, dim = locations.shape
cdef np.ndarray[double, ndim=2, mode="c"] c_locations = \
np.ascontiguousarray(locations, dtype=ctypes.c_double)
cdef double** loc = <double **>malloc(n_points * sizeof(double *))
cdef int point;
for point in range(n_points):
loc[point] = &c_locations[point, 0]
# Initialise output array for the local affines
return_output = False
if output is None:
output = np.zeros((n_points, 6), dtype=ctypes.c_double)
return_output = True
elif hasattr(output, "__array__"):
output = np.asanyarray(output)
assert output.shape == (n_points, 6)
else:
raise ValueError(f"Expected array for output, "
f"got {output.__class__.__name__} instead.")
cdef np.ndarray[double, ndim=2, mode="c"] c_jacobians = \
np.ascontiguousarray(output, dtype=ctypes.c_double)
cdef double** jac = <double **>malloc(n_points * sizeof(double *))
for point in range(n_points):
jac[point] = &c_jacobians[point, 0]
# Generate input coordinates of the grid points (physical space)
fwd_grid_coordinates = np.stack(
np.mgrid[tuple(slice(ax) for ax in forward_shape)],
axis=-1).reshape((-1, 2)).astype(ctypes.c_double)
if callable(pretx):
fwd_grid_coordinates = pretx(fwd_grid_coordinates)
cdef np.ndarray[double, ndim=2, mode="c"] c_input_coordinates = \
np.ascontiguousarray(fwd_grid_coordinates, dtype=ctypes.c_double)
cdef double** ic = <double **>malloc(n_points * sizeof(double *))
for point in range(n_fwd_points):
ic[point] = &c_input_coordinates[point, 0]
# Find nearest grid point to each location
tree = KDTree(fwd_grid_coordinates, leaf_size=2)
cdef np.ndarray[int, ndim=1, mode="c"] nearest_neighbour = \
np.ascontiguousarray(tree.query(locations, k=1)[1].ravel(),
dtype=ctypes.c_int)
cdef int* nn = <int *>malloc(n_points * sizeof(int))
nn = &nearest_neighbour[0]
# Generate mapped coordinates of the grid points (physical space)
mapped_coordinates = fwd_grid_coordinates.copy()
if vectorder != (0, 1):
for i, ax in enumerate(vectorder):
mapped_coordinates[:, ax] += field.reshape((dim, -1))[i]
else:
mapped_coordinates[...] += field.reshape((dim, -1)).T
cdef np.ndarray[double, ndim=2, mode="c"] c_mapped_coordinates = \
np.ascontiguousarray(mapped_coordinates, dtype=ctypes.c_double)
del fwd_grid_coordinates
del mapped_coordinates
cdef double** mc = <double **>malloc(n_fwd_points * sizeof(double *))
for point in range(n_fwd_points):
mc[point] = &c_mapped_coordinates[point, 0]
# Call C function to calculate the local affine for each location
CalculateLocalAffine2D(loc, nn, n_points, fsh, ic, mc, jac, maxiter)
# Generate output
if return_output:
return np.asanyarray(c_jacobians, dtype=field.dtype)
else:
output.flat[...] = \
np.asanyarray(c_jacobians, dtype=output.dtype).ravel()
def jacobian2d(field, locations, pretx=None, vectorder=(0, 1), maxiter=1000):
# Obtain the properties of the input field (this must be absolute)
......@@ -424,88 +340,109 @@ def invert3d(field, output=None, vectorder=(0, 1, 2), pretx=None, posttx=None,
invwarp.reshape((dim, *backward_shape)).ravel(), dtype=output.dtype)
def local_affine3d(field, locations, output=None, pretx=None,
vectorder=(0, 1, 2), maxiter=1000):
# Obtain the properties of the input field (this must be absolute)
dim, *forward_shape = field.shape
cdef np.ndarray[int, ndim=1, mode="c"] fwd_shape = \
np.asarray(forward_shape, dtype=ctypes.c_int)
cdef int* fsh = <int *>malloc(3 * sizeof(int))
fsh = &fwd_shape[0]
# assert dim == 3, "Input is not a 3D vector field."
n_fwd_points = np.product(fwd_shape)
# Obtain locations
n_points, dim = locations.shape
cdef np.ndarray[double, ndim=2, mode="c"] c_locations = \
np.ascontiguousarray(locations, dtype=ctypes.c_double)
cdef double** loc = <double **>malloc(n_points * sizeof(double *))
cdef int point;
for point in range(n_points):
loc[point] = &c_locations[point, 0]
# Initialise output array for the local affines
def local_affine(field, newlocations, gridlocations, output=None,
vectorder=None, maxiter=1000):
"""
Returns the local affine approximation of a field at the given
physical locations.
"""