Commit 0cad5e38 authored by inhuszar's avatar inhuszar
Browse files

Fixed 3D field inversion bug in C code.

parent a530cf6d
......@@ -2,3 +2,4 @@
/dist/
/examples/
/tirl.egg-info/
/.idea/
#!/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
import matplotlib
matplotlib.use("tkagg")
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 jacobian2d, jacobian3d
from tirl.transformations.nonlinear.displacement import TxDisplacementField
def sine_field(xx, yy, zz=None):
h, w, *d = xx.shape
xx = h / 10. * np.sin(xx / (h - 1) * 2 * np.pi)
yy = w / 10. * np.sin(yy / (w - 1) * 2 * np.pi)
if zz is not None:
zz = d[0] / 10. * np.sin(zz / (d[0] - 1) * 2 * np.pi)
return xx, yy, zz
else:
return xx, yy
def cosine_field(xx, yy, zz=None):
h, w, *d = xx.shape
xx = h / 10. * np.cos(xx / (h - 1) * 2 * np.pi)
yy = w / 10. * np.cos(yy / (w - 1) * 2 * np.pi)
if zz is not None:
zz = d[0] / 10. * np.cos(zz / (d[0] - 1) * 2 * np.pi)
return xx, yy, zz
else:
return xx, yy
def gaussian_field(xx, yy, zz=None):
h, w, *d = xx.shape
d = d[0] if d else d
print(d)
g1 = gaussian(h, std=h/10.)
g2 = gaussian(w, std=w/10.)
if d:
g3 = gaussian(d, std=d/10.)
if zz is not None:
xx = h / 10. * g1[:, np.newaxis, np.newaxis] \
* g1[np.newaxis, :, np.newaxis] \
* g1[np.newaxis, np.newaxis, :]
yy = w / 10. * g2[:, np.newaxis, np.newaxis] \
* g2[np.newaxis, :, np.newaxis] \
* g2[np.newaxis, np.newaxis, :]
zz = d / 10. * g3[:, np.newaxis, np.newaxis] \
* g3[np.newaxis, :, np.newaxis] \
* g3[np.newaxis, np.newaxis, :]
return xx, yy, zz
else:
xx = h / 10. * g1[:, np.newaxis] * g1[np.newaxis, :]
yy = w / 10. * g2[:, np.newaxis] * g2[np.newaxis, :]
return xx, yy
def jacmap(xx, yy, zz=None):
gradx = np.stack(np.gradient(xx), axis=-1)
grady = np.stack(np.gradient(yy), axis=-1)
if zz is not None:
gradz = np.stack(np.gradient(zz), axis=-1)
J = np.eye(3)[np.newaxis] + \
np.stack((gradx.reshape((-1, 1, 3)),
grady.reshape((-1, 1, 3)),
gradz.reshape((-1, 1, 3))), axis=1)
else:
J = np.eye(2)[np.newaxis] + \
np.stack((gradx.reshape((-1, 1, 2)),
grady.reshape((-1, 1, 2))), axis=1)
return np.linalg.det(J)
def generate_field(shape, type="sine", display=True):
"""
Returns (?, height, width) field.
"""
# Create grid
h, w, *d = shape
d = d[0] if d else d
if d:
xx, yy, zz = np.mgrid[0:h, 0:w, 0:d]
else:
xx, yy = np.mgrid[0:h, 0:w]
zz = None
# Sine (Diffeomorphic)
if type == "sine":
xx, yy, *zz = sine_field(xx, yy, zz)
# Cosine (Diffeomorphic, but tears edges)
elif type == "cosine":
xx, yy, *zz = cosine_field(xx, yy, zz)
# Gaussian (non-diffeomorphic)
elif type == "gaussian":
xx, yy, *zz = gaussian_field(xx, yy, zz)
else:
raise ValueError("Available field prototypes: sine, cosine, gaussian.")
zz = zz[0] if zz else None
# Calculate Jacobian determinant map
jac = jacmap(xx, yy, zz)
jac = np.linalg.norm(jac.reshape(*shape, -1), axis=-1)
# Save Jacobian determinant map as NIfTI
# nifti = nib.Nifti1Image(jac, np.eye(4), nib.Nifti1Header())
# nib.save(nifti, "jac.nii.gz")
if display:
# Display deformation map
fig, axes = plt.subplots(2, 2, figsize=(12, 12))
axes = axes.ravel()
if d:
centre = (slice(None), slice(None), d//2)
else:
centre = Ellipsis
# Vertical warp
xvec = axes[0].imshow(xx[centre], cmap="gray")
axes[0].set_title("Vertical")
# Horizontal warp
yvec = axes[1].imshow(yy[centre], cmap="gray")
axes[1].set_title("Horizontal")
# Depth warp
if zz is not None:
zvec = axes[2].imshow(zz[centre], cmap="gray")
axes[2].set_title("Depth")
# Jacobian
jacplot = axes[3].imshow(jac[centre], cmap="gray")
axes[3].set_title("Jacobian")
# Vertical warp colorbar
divider = make_axes_locatable(axes[0])
cax = divider.append_axes("right", size="5%", pad=0.05)
fig.colorbar(xvec, cax=cax)
# Horizontal warp colorbar
divider = make_axes_locatable(axes[1])
cax = divider.append_axes("right", size="5%", pad=0.05)
fig.colorbar(yvec, cax=cax)
# Depth warp colorbar
if zz is not None:
divider = make_axes_locatable(axes[2])
cax = divider.append_axes("right", size="5%", pad=0.05)
fig.colorbar(zvec, cax=cax)
# Jacobian
divider = make_axes_locatable(axes[3])
cax = divider.append_axes("right", size="5%", pad=0.05)
fig.colorbar(jacplot, cax=cax)
plt.tight_layout()
plt.show()
if zz is not None:
return np.stack((xx, yy, zz), axis=0)
else:
return np.stack((xx, yy), axis=0)
def inv2d():
""" Main program code. """
h, w = shape = (128, 128)
# h, w = shape = (1000, 1000)
field = generate_field(shape, type="gaussian")
print("Inverting field...")
invwarp = np.zeros((2, *shape), dtype=np.float64)
st = time()
invert2d(field, invwarp, fill_holes=True)
print(f"Inversion time: {time() - st} seconds")
# Test
# x = np.random.rand(200, 2) * np.asarray(shape)[np.newaxis, :]
xx, yy = np.mgrid[0:h, 0:w]
x = np.stack((xx.ravel(), yy.ravel()), axis=1)
fwd = TxDisplacementField(field)
bwd = TxDisplacementField(invwarp)
err = bwd.map(fwd.map(x)) - x
print(err)
print(np.median(np.abs(err), axis=0))
print(np.max(np.abs(err), axis=0))
fig, axes = plt.subplots(1, 3)
axes[0].imshow(invwarp[0], cmap="gray")
axes[1].imshow(invwarp[1], cmap="gray")
axes[2].imshow(np.linalg.norm(err.reshape(*shape, 2), axis=-1), cmap="gray")
plt.show()
def inv3d():
""" Main program code. """
# h, w, d = shape = (32, 32, 32)
h, w, d = shape = (64, 64, 64)
field = generate_field(shape, type="gaussian")
print("Inverting field...")
invwarp = np.zeros((3, *shape), dtype=np.float64)
st = time()
invert3d(field, invwarp, fill_holes=True)
print(f"Inversion time: {time() - st} seconds")
# Test
xx, yy, zz = np.mgrid[0:h, 0:w, 0:d]
x = np.stack((xx.ravel(), yy.ravel(), zz.ravel()), axis=1)
fwd = TxDisplacementField(field, mode="abs")
print(f"Invalid: {np.count_nonzero(~np.isfinite(invwarp))}")
# invwarp[~np.isfinite(invwarp)] = 0
bwd = TxDisplacementField(invwarp, mode="abs")
err = bwd.map(fwd.map(x)) - x
print(err)
print(np.median(np.abs(err), axis=0))
print(np.max(np.abs(err), axis=0))
fig, axes = plt.subplots(2, 2, figsize=(12, 12))
axes = axes.ravel()
centre = (slice(None), slice(None), d//2)
axes[0].imshow(invwarp[0][centre], cmap="gray")
axes[1].imshow(invwarp[1][centre], cmap="gray")
axes[2].imshow(invwarp[2][centre], cmap="gray")
err = np.linalg.norm(err.reshape(*shape, 3), axis=-1)
axes[3].imshow(err[centre], cmap="gray")
plt.show()
# Save inverse consistency map as NIfTI
# nifti = nib.Nifti1Image(err, np.eye(4), nib.Nifti1Header())
# nib.save(nifti, "err.nii.gz")
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)
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])
jacdet = np.linalg.det(jac)
print(f"Jacobian range: {np.min(jacdet)} - {np.max(jacdet)}")
plt.imshow(jacdet.reshape(shape), cmap="gray")
plt.colorbar()
plt.show()
def jac3d():
h, w, d = shape = (64, 64, 64)
f = generate_field(shape, type="gaussian")
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]
jacdet = np.linalg.det(jac)
# jacdet[~np.isfinite(jacdet)] = 0
print(f"Jacobian range: {np.min(jacdet)} - {np.max(jacdet)}")
plt.imshow(jacdet.reshape(shape)[:, :, d // 2], cmap="gray")
plt.colorbar()
plt.show()
def txjac2d():
h, w = shape = (256, 256)
f = generate_field(shape, type="sine", display=False)
locations = np.stack(np.mgrid[0:h, 0:w], axis=-1).reshape((-1, 2))
locations = locations + 0.1 * np.ones((1, 2), dtype=np.float64)
jac, vi = jacobian2d(f, locations, pretx=None, vectorder=(0, 1),
maxiter=1000)
print(jac[:20])
print(vi[:20])
# Create stacked sparse parameter-Jacobian matrix
mask = np.all(vi >= 0, axis=1)
vi = vi[mask, ...]
u = jac[mask, 0, :]
v = jac[mask, 1, :]
print(f"Negatives: {np.count_nonzero(vi < 0)}")
n_points = np.count_nonzero(mask)
n_nodes = h * w
n_parameters = 2 * n_nodes
combined = np.concatenate((u, v), axis=0).ravel()
col = np.tile(np.concatenate((vi, n_nodes + vi), axis=1), (2, 1)).ravel()
row = np.repeat(np.arange(2 * n_points).reshape(-1, 1), 6, axis=1).ravel()
txjac = sp.csr_matrix((combined, (row, col)), shape=(n_points * 2, n_parameters))
print(txjac)
def txjac3d():
h, w, d = shape = (64, 64, 64)
n_nodes = h * w * d
f = generate_field(shape, type="sine", display=False)
locations = np.stack(np.mgrid[0:h, 0:w, 0:d], axis=-1).reshape((-1, 3))
locations = locations + 0.1 * np.ones((1, 3), dtype=np.float64)
jac, vi = jacobian3d(f, locations, pretx=None, vectorder=(0, 1, 2),
maxiter=1000)
print(jac[:-20])
print(vi[:-20])
# Create stacked sparse parameter-Jacobian matrix
print(f"Negatives: {np.count_nonzero(np.any(vi < 0, axis=-1))}")
mask = np.all(vi >= 0, axis=1)
vi = vi[mask, ...]
u = jac[mask, 0, :]
v = jac[mask, 1, :]
w = jac[mask, 2, :]
print(f"Negatives: {np.count_nonzero(vi < 0)}")
n_points = np.count_nonzero(mask)
n_parameters = 3 * n_nodes
combined = np.concatenate((u, v, w), axis=0).ravel()
print(combined[:-36])
print(vi[:-36])
col = np.tile(np.concatenate((vi, n_nodes + vi, 2 * n_nodes + vi), axis=1), (3, 1)).ravel()
row = np.repeat(np.arange(3 * n_points).reshape(-1, 1), 12, axis=1).ravel()
txjac = sp.csr_matrix((combined, (row, col)), shape=(n_points * 3, n_parameters))
print(txjac)
if __name__ == "__main__":
inv2d()
inv3d()
jac2d()
jac3d()
txjac2d()
txjac3d()
......@@ -36,7 +36,7 @@ def main():
cov = coverage.coverage()
cov.start()
loader = unittest.TestLoader()
unit_tests = unittest.TestLoader().discover("unit_testing", "test_*.py")
unit_tests = unittest.TestLoader().discover("unittests", "test_*.py")
result = unittest.TextTestRunner(verbosity=2).run(unit_tests)
if result.wasSuccessful():
cov.stop()
......
......@@ -261,7 +261,8 @@ class TestRepresentation(unittest.TestCase):
tx_embed, tx_scale3d, tx_trans3d)
def test_repr(self):
expected = f"Chain{str(self.transformations)}"
listing = "\n".join(tuple(str(tx) for tx in self.transformations))
expected = f"Chain[\n{listing}\n]"
self.assertEqual(str(self.chain), expected)
......
......@@ -8,6 +8,9 @@ from tirl.costs.mi import CostMI
from tirl.transformations.linear.rotation import TxRotation2D
resourcedir = tirl.home("..", "tests", "resources")
class TestConstruction(unittest.TestCase):
def setUp(self) -> None:
......@@ -40,7 +43,7 @@ class TestScalar(unittest.TestCase):
cost = CostMI(self.rot, self.orig)
fullcost = cost()
self.orig.mask = TImage(os.path.join(
tirl.home(), "../tests/resources/testimg_edge_mask.tif"))
resourcedir, "testimage", "testimg_edge_mask.tif"))
cost = CostMI(self.rot, self.orig)
self.assertNotEqual(fullcost, cost())
......@@ -66,7 +69,7 @@ class TestGradient(unittest.TestCase):
fulljac = cost.dx()
self.assertTrue(np.all(np.isfinite(fulljac)))
self.orig.mask = TImage(os.path.join(
tirl.home(), "../tests/resources/testimg_edge_mask.tif"))
resourcedir, "testimage", "testimg_edge_mask.tif"))
cost = CostMI(self.rot, self.orig)
jac = cost.dx()
self.assertFalse(np.allclose(jac, fulljac))
......
......@@ -8,6 +8,9 @@ from tirl.costs.mind import CostMIND
from tirl.transformations.linear.rotation import TxRotation2D
resourcedir = tirl.home("..", "tests", "resources")
class TestConstruction(unittest.TestCase):
def setUp(self) -> None:
......@@ -40,7 +43,7 @@ class TestScalar(unittest.TestCase):
cost = CostMIND(self.rot, self.orig)
fullcost = cost()
self.orig.mask = TImage(os.path.join(
tirl.home(), "../tests/resources/testimg_edge_mask.tif"))
resourcedir, "testimage", "testimg_edge_mask.tif"))
cost = CostMIND(self.rot, self.orig)
self.assertNotEqual(fullcost, cost())
......@@ -69,7 +72,7 @@ class TestGradient(unittest.TestCase):
fulljac = cost.dx()
self.assertTrue(np.all(np.isfinite(fulljac)))
self.orig.mask = TImage(os.path.join(
tirl.home(), "../tests/resources/testimg_edge_mask.tif"))
resourcedir, "testimage", "testimg_edge_mask.tif"))
cost = CostMIND(self.rot, self.orig)
jac = cost.dx()
self.assertFalse(np.allclose(jac, fulljac))
......
import unittest
import os
import warnings
import numpy as np
from tirl.parameters import ParameterVector
from tirl.parameters import ParameterGroup
warnings.filterwarnings("default")
class TestParameterVector(unittest.TestCase):
def test_construction_from_scalararray(self):
......@@ -99,7 +103,7 @@ class TestParameterVector(unittest.TestCase):
self.assertTrue(np.all(p.get_lower_bounds() == q.get_lower_bounds()))
self.assertTrue(np.all(p.get_upper_bounds() == q.get_upper_bounds()))
self.assertTrue(np.all(q.get_bounds(7) == [-10, 10]))
self.assertEqual(q.name, name + "_copy")
self.assertEqual(q.name, name)
def test_saveload(self):
p = ParameterVector(range(10))
......@@ -421,7 +425,7 @@ class TestParameterGroup(unittest.TestCase):
self.assertTrue(np.all(gc.get_lower_bounds() == g.get_lower_bounds()))
self.assertTrue(np.all(gc.get_upper_bounds() == g.get_upper_bounds()))
self.assertTrue(np.all(gc.get_bounds(23) == [-100, 100]))
self.assertEqual(gc.name, name + "_copy")
self.assertEqual(gc.name, name)
def test_saveload(self):
p = ParameterVector(range(10))
......
......@@ -261,7 +261,7 @@ def home(*args):
def testimg():
""" Returns a test TImage. """
from tirl.timage import TImage
impath = os.path.join(home(), "../tests/resources/testimg.png")
impath = os.path.join(home(), "../tests/resources/testimage/testimg.png")
img = TImage.fromfile(impath, tensor_axes=(2,))
return img
......
......@@ -244,8 +244,8 @@ int CalculateInverseWarp2D(double ** SrcLocations, double ** SrcMapped,
cblas_dgemv(CblasColMajor, CblasNoTrans, 2, 3, 1.0, ReducedAffine, 2, HomoTargetPoint, 1, 0.0, BwdMappedTargetPoint, 1);
InvWarp[0][Point] = BwdMappedTargetPoint[0] - HomoTargetPoint[0];
InvWarp[1][Point] = BwdMappedTargetPoint[1] - HomoTargetPoint[1];
free(Affine);
}
free(Affine);
// If the matching simplex could not be found,
// leave these values undefined and fill them later by local averaging.
......
......@@ -321,8 +321,8 @@ int CalculateInverseWarp3D(double ** SrcLocations, double ** SrcMapped,
for (int Vertex = 0; Vertex < 4; Vertex++)
GridSimplex[Vertex] = calloc(sizeof(int), 3);
// Declare simplex to operate on the space of forward-mapped grid points
double** SrcMappedSimplex = calloc(sizeof(double *), 3);
double** SrcLocationSimplex = calloc(sizeof(double *), 3);
double** SrcMappedSimplex = calloc(sizeof(double *), 4);
double** SrcLocationSimplex = calloc(sizeof(double *), 4);
double * TargetPoint;
// Loop over all grid points of the output domain
......@@ -418,8 +418,8 @@ int CalculateInverseWarp3D(double ** SrcLocations, double ** SrcMapped,
InvWarp[0][Point] = BwdMappedTargetPoint[0] - HomoTargetPoint[0];
InvWarp[1][Point] = BwdMappedTargetPoint[1] - HomoTargetPoint[1];
InvWarp[2][Point] = BwdMappedTargetPoint[2] - HomoTargetPoint[2];
free(Affine);
}
free(Affine);
// If the matching simplex could not be found,
// leave these values undefined and fill them later by local averaging.
......
......@@ -844,7 +844,7 @@ class ParameterVector(TIRLObject, np.lib.mixins.NDArrayOperatorsMixin):
# If the parameter bounds close up on certain parameters, lock them.
if (lower_bounds is not None) and (upper_bounds is not None):
equality = (lower_bounds - upper_bounds) >= EPS
equality = (lower_bounds - upper_bounds) >= 0
if np.any(equality):
import warnings
newlock = np.flatnonzero(equality)
......@@ -928,7 +928,7 @@ class ParameterVector(TIRLObject, np.lib.mixins.NDArrayOperatorsMixin):
"the current lower bounds specification.")
# If the parameter bounds close up on certain parameters, lock them.
equality = (lower_bounds - self.get_upper_bounds()) >= EPS
equality = (lower_bounds - self.get_upper_bounds()) >= 0
if np.any(equality):
import warnings
newlock = np.flatnonzero(equality)
......@@ -1014,7 +1014,7 @@ class ParameterVector(TIRLObject, np.lib.mixins.NDArrayOperatorsMixin):
"the current upper bounds specification.")
# If the parameter bounds close up on certain parameters, lock them.
equality = (self.get_lower_bounds() - upper_bounds) >= EPS
equality = (self.get_lower_bounds() - upper_bounds) >= 0
if np.any(equality):
import warnings
newlock = np.flatnonzero(equality)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment