Commit f29bab73 authored by inhuszar's avatar inhuszar
Browse files

TField/TImage resampling tested. Bugfixes, preparing TField for release.

parent 0a6990dc
#!/usr/bin/env python3
from tirl.timage import TImage
from tirl.transformations.linear import TxRotation2D
from itertools import product
import numpy as np
def main():
# Load a large image
hres = TImage("/mnt/nvme/resources/NP391-16-CD68-20x.svs", level=2,
name="histology", n_threads=8)
print("hres: {}".format(hres))
# Start simple: downscale the whole image
lowres = hres.resample(0.25)
print(lowres)
lowres.preview()
# Extract the bottom-right corner
hres_br = hres.voxels[4000:, 4000:]
print("hres_br: {}".format(hres_br))
# Resample the extracted corner
hres_br.resample(0.25, copy=False)
print("hres_br after resampling: {}".format(hres_br))
# Extract a rectangle from the inner area
rect = hres_br.voxels[400:1000, 600:1000]
print("sq: {}".format(rect))
# As a final trick, rotate the rectangle around its centre
rect.centralise()
tx_rot = TxRotation2D(180)
rect.domain.transformations.append(tx_rot)
# Generate high-res image for the rotated rectangle
sq_hr = rect.highres(interpolate=False)
print("sq_hr: {}".format(sq_hr))
# Compare the physical coordinates of the corners. These should be equal.
corners = np.stack(product(*zip((0, 0), rect.vshape)), axis=0)
print("Low-res coordinates:")
print(np.hstack((corners, rect.domain.map_voxel_coordinates(corners))))
corners_hr = np.stack(product(*zip((0, 0), sq_hr.vshape)), axis=0)
print("High-res coordinates:")
print(np.hstack(
(corners_hr, sq_hr.domain.map_voxel_coordinates(corners_hr))))
# Show the lowres and the highres tiles
sq.preview()
sq_hr.preview()
if __name__ == "__main__":
main()
......@@ -28,12 +28,13 @@ print("Opening time: {} s".format(time() - st))
img.preview()
st = time()
img.resample(0.25)
roi = img.voxels[2000:3000, 2000:3000]
print("Resampling time: {} s".format(time() - st))
print(img)
roi = img.voxels[1000:2000, 2000:3000]
roi.preview()
roi.highres.preview()
exit()
print("Resampling time: {} s".format(time() - st))
print(img)
print(img._data.fname)
img.preview()
# img.voxels[20000:22000, 30000:32000].preview()
......@@ -5,12 +5,17 @@ from tirl.timage import TImage
def main():
imfile = "/Users/inhuszar/nilsson/Histology.tif"
img = TImage(imfile, storage="hdd")
# imfile = "/Users/inhuszar/nilsson/Histology.tif"
# img = TImage(imfile, storage="hdd")
imfile = "../resources/testimg.png"
img = TImage(imfile)
print(img)
img.preview()
st = time()
img.resample(0.1)
print(img.shape)
print(img)
print("{} s.".format(time() - st))
img.preview()
if __name__ == "__main__":
......
import unittest
import numpy as np
from tirl.tfield import TField
class TFieldUfuncBehaviour(unittest.TestCase):
def setUp(self):
xx, yy, zz, ww = \
np.meshgrid(range(200), range(200), range(10), range(2),
indexing="ij")
arr = np.sin(xx / 200 * 4 * np.pi) + np.sin(yy / 200 * 2 * np.pi)
self.tf = TField.fromarray(arr, tensor_axes=(-2, -1))
self.tf.voxels[50, 50] = -10
self.tf.voxels[100, 100] = 10
def test_argmin_argmax(self):
minloc = np.argmin(self.tf.data)
maxloc = np.argmax(self.tf.data)
self.assertEqual(np.argmin(self.tf), minloc)
self.assertEqual(np.argmax(self.tf), maxloc)
self.assertEqual(np.argmin(self.tf.tensors),)
if __name__ == '__main__':
unittest.main()
......@@ -317,10 +317,12 @@ class Cost(TIRLObject):
return self.target.mask
# Return adapted source mask if the target mask is not defined
elif sm and not tm:
m = self.source.fromarray(self.source.mask)
return self.source.mask.evaluate(self.target.domain)
# Create composite mask (respecting maskmode) if both masks are defined
else:
source_mask = self.source.mask.evaluate(self.target.domain)
m = self.source.fromarray(self.source.mask)
source_mask = m.evaluate(self.target.domain)
if self.maskmode == "and":
cm = (source_mask * self.target.mask) ** 0.5
# TODO: This is for in-memory arrays only.
......
......@@ -147,7 +147,7 @@ class CostSSD(Cost):
n_tensor = reduce(mul, target.tshape)
# Regardless of the normalisation setting, the Euclidean distance
# is always normalised by the number of tensor elements.
result = calc_distance(result.asTField()) / n_tensor
result = calc_distance(result) / n_tensor
# Modulate distance (or difference) image by the composite mask
mask = self.combine_masks()
result = result * mask if mask is not None else result
......
......@@ -1892,6 +1892,38 @@ class Domain(TIRLObject):
else:
return tuple(txs)
def copy_transformations(self, old_domain, links=True):
"""
Copies the transfomation chain from an existing domain to the current
domain. This method is useful when a TField or a TImage is
resampled to a new grid.
:param old_domain: template domain
:type old_domain: Domain
:param links:
If True, Non-linear transformations that were dynamically linked to
the other domain, will be recreated so that they will be also
linked to the current domain. If False, these transformations will
be statically linked to the new domain.
:type links: bool
"""
self.transformations = old_domain.transformations
if links:
for i, tx in enumerate(old_domain.transformations):
if hasattr(tx, "domain"):
# If the transformation was dynamically linked to
# the current domain, re-grid that transformation.
if tx.domain == old_domain[:i]:
# Erase transformations on the tx domain
tx.metaparameters["domain"] = old_domain[:0]
# Resample among voxel domains
new_tx = tx.regrid(self[:0])
# Add back all new antecedent transformations
new_tx.metaparameters["domain"] = self[:i]
# Add the current transformation to the new domain
self.transformations[i] = new_tx
if __name__ == "__main__": # pragma: no cover
print("""This module is not intended for execution.""")
......@@ -436,7 +436,7 @@ class Interpolator(TIRLObject):
points = points - offset # note: dtype mismatch if -= is used
# Perform interpolation
n_threads = reduce(mul, tensor_shape)
n_threads = reduce(mul, tensor_shape) if tensor_shape else 1
p = mt.Pool(processes=n_threads)
args = (points, subfield, indices, out, tensor_shape)
job_func = partial(self._finish, *args)
......
......@@ -28,8 +28,12 @@ class ScipyNearestNeighbours(Interpolator):
coordinates = np.asarray(coordinates).reshape(1, -1)
try:
kwargs.pop("hold", None) # bugfix
vcoords = np.stack(
np.meshgrid(*tuple(range(dim) for dim in input_array.shape),
indexing="ij"), axis=-1)
ip = NearestNDInterpolator(
coordinates, input_array.ravel(), rescale=False)
vcoords.reshape(-1, input_array.ndim), input_array.ravel(),
rescale=False)
res = ip(coordinates)
except:
print("Interpolation failed with the following input: ",
......
......@@ -36,7 +36,7 @@ from tirl import settings as ts
from tirl import exceptions as te
from tirl.tirlobject import TIRLObject
from tirl.costs.cost import Cost
from tirl.transformations.basic import Transformation, TransformationGroup
from tirl.transformations import Transformation, TransformationGroup
from tirl.transformations.basic import TxIdentity
from tirl.regularisers.regulariser import Regulariser
......
......@@ -8,7 +8,7 @@
import os
# Temporary working directory (use absolute path!)
TWD = "/Users/inhuszar/temp"
TWD = "/mnt/nvme/temp"
if not os.path.isdir(TWD):
os.makedirs(TWD)
......@@ -87,6 +87,8 @@ DEFAULT_COMPACT_INTERPOLATOR = \
"tirl.interpolators.scipyinterpolator.ScipyInterpolator"
DEFAULT_NONCOMPACT_INTERPOLATOR = \
"tirl.interpolators.rbfinterpolator.RbfInterpolator"
TFIELD_PRESMOOTH = True
TFIELD_PRESMOOTH_KERNELSIZE_NSIGMA = 3
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
# TImage class #
......@@ -126,4 +128,4 @@ HALTON_SEED = 1
# Visualisation #
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
ENABLE_VISUALISATION = True
MPL_BACKEND = "macosx"
MPL_BACKEND = "tkagg"
This diff is collapsed.
This diff is collapsed.
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