Commit 7a44410b authored by ihuszar's avatar ihuszar
Browse files

Slice-to-volume works with affine re-estimation.

parent 1845fb4c
from protocols import tirl_affine2d
from protocols import tirl_affine3d
from protocols import tirl_rigid2d
from protocols import tirl_rigid3d
from protocols import tirl_rigid2d3d
from protocols import tirl_affine2d3d
from protocols import tirl_nonlin2d
# from protocols import tirl_affine2d
# from protocols import tirl_affine3d
# from protocols import tirl_rigid2d
# from protocols import tirl_rigid3d
# from protocols import tirl_rigid2d3d
# from protocols import tirl_affine2d3d
# from protocols import tirl_nonlin2d
{
"header": {
"title": "Slice-to-volume image registration routine (BigMac dataset)",
"description": "Pipeline to register histology whole slides to volumetric MRI data",
"author": "Istvan N Huszar, Amy FD Howard"
},
"slice": {
"file": "/Users/inhuszar/bigmac/histo/H092x/mosaic_colour_mrires.tif",
"alternative": "/Users/inhuszar/bigmac/reg/stage4.timg",
"storage": "mem",
"mask": {
"file": null,
"normalise": false
},
"resolution": 0.3,
"preview": false,
"actions": ["resample_image", "lab_b"]
},
"volume": {
"file": "/Users/inhuszar/bigmac/MRI_correct/struct.nii.gz",
"storage": "mem",
"mask": {
"file": null,
"normalise": false
},
"resolution": 0.3,
"preview": false,
"actions": [],
"export": true
},
"general": {
"verbosity": "debug",
"system": "macos",
"logfile": "/Users/inhuszar/bigmac/reg/tirl_slice_to_volume.log",
"outdir": "/Users/inhuszar/bigmac/reg",
"stages": [5, 2, 3, 4],
"stage_index": 8,
"snapshot_ext": "png",
"randomseed": 1,
"warnings": false
},
"regparams": {
"init": {
"scale2d": [1.0, 1.0],
"rot2d": -90,
"trans2d": [0, 0],
"rot3d": [0, 0, 90],
"affine3d": [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0],
"trans3d": [0, 0, 0]
},
"stage_1": {
"export": {
"image": true,
"snapshot": true
},
"visualise": true,
"verbose": 4,
"slab": {
"centre": [0, -4.3, 0],
"normal": [0, 1, 0],
"thickness": 1.5,
"inits": 5
},
"iterations": 1,
"scaling": [8, 4, 4, 2, 2, 1, 1],
"smoothing": [0, 1, 0, 1, 0, 1, 0],
"constrained": true,
"opt_step": 0.1,
"x0": {
"scale2d": [1.0, 1.0],
"rot2d": -90,
"trans2d": [0, 0],
"rot3d": [0, 0, 90],
"affine3d": [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0],
"trans3d": [0, 0, 0]
},
"dx0": {
"scale2d_lower_delta": [0.2, 0.2],
"scale2d_upper_delta": [0.2, 0.2],
"rot2d_lower_delta": 30.0,
"rot2d_upper_delta": 30.0,
"trans2d_lower_delta": [10.0, 10.0],
"trans2d_upper_delta": [10.0, 10.0],
"warp_lower_dxy": 2.5,
"warp_lower_dz": 2.5,
"warp_upper_dxy": 2.5,
"warp_upper_dz": 2.5,
"rot3d_lower_delta": [20.0, 20.0, 20.0],
"rot3d_upper_delta": [20.0, 20.0, 20.0],
"affine3d_lower_delta": [0.1, 0.1, 0.1, 2.5, 0.1, 0.1, 0.1, 2.5, 0.1, 0.1, 0.1, 2.5],
"affine3d_upper_delta": [0.1, 0.1, 0.1, 2.5, 0.1, 0.1, 0.1, 2.5, 0.1, 0.1, 0.1, 2.5],
"trans3d_lower_delta": [0.3, 0.3, 0.3],
"trans3d_upper_delta": [0.3, 0.3, 0.3]
},
"dx": {
"scale2d_lower_delta" : [0.05, 0.05],
"scale2d_upper_delta" : [0.05, 0.05],
"rot2d_lower_delta" : 5,
"rot2d_upper_delta" : 5,
"trans2d_lower_delta" : [1.0, 1.0],
"trans2d_upper_delta" : [1.0, 1.0],
"warp_lower_dxy": 2.0,
"warp_lower_dz": 1.0,
"warp_upper_dxy": 2.0,
"warp_upper_dz": 1.0,
"rot3d_lower_delta" : [5.0, 5.0, 5.0],
"rot3d_upper_delta" : [5.0, 5.0, 5.0],
"affine3d_lower_delta": [0.1, 0.1, 0.1, 1.0, 0.1, 0.1, 0.1, 1.0, 0.1, 0.1, 0.1, 1.0],
"affine3d_upper_delta": [0.1, 0.1, 0.1, 1.0, 0.1, 0.1, 0.1, 1.0, 0.1, 0.1, 0.1, 1.0],
"trans3d_lower_delta" : [0.3, 0.3, 0.3],
"trans3d_upper_delta" : [0.3, 0.3, 0.3]
}
},
"stage_2": {
"export": {
"image": true,
"snapshot": true
},
"visualise": true,
"verbose": 4,
"scaling": [4, 2, 1, 1],
"smoothing": [0, 0, 1, 0],
"lower_delta": 0.1,
"upper_delta": 0.1,
"xtol_rel": 0.01,
"opt_step": 0.5
},
"stage_3": {
"export": {
"image": true,
"snapshot": true,
"halton": true
},
"n_points": 16,
"visualise": true,
"verbose": 4,
"smoothing": [2, 1, 0],
"lower_dz": 3,
"upper_dz": 3,
"reg_weight": 0,
"opt_step": 0.2,
"xtol_abs": 0.1,
"randomseed": 1
},
"stage_4": {
"export": {
"image": true,
"snapshot": true,
"halton": true
},
"n_points": 16,
"visualise": false,
"verbose": 4,
"smoothing": [2],
"lower_dxy": 3,
"lower_dz": 3,
"upper_dxy": 3,
"upper_dz": 3,
"reg_weight": 0,
"opt_step": 0.2,
"xtol_abs": 0.1
},
"stage_5": {
"export": {
"image": true,
"snapshot": true,
"mask": true
},
"visualise": false,
"mask": {
"lthr": 0.1,
"uthr": 1
}
}
}
}
This diff is collapsed.
......@@ -9,11 +9,12 @@ import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
import tirl
import tirl.transformations.basic as ttb
import tirl.transformations.linear as ttl
from tirl.domain import Domain
from tirl.timage import TImage
from tirl.transformations.basic.identity import TxIdentity
from tirl.transformations.linear.affine import TxAffine
from tirl.transformations.linear.translation import TxTranslation
from tirl.transformations.linear.axisangle import TxRotationAxisAngle
def set_args(parser):
......@@ -102,8 +103,8 @@ def load_mri(mrifile, with_tx=True):
mri = TImage(nifti.get_data(), name="mri")
if with_tx:
affine = nifti.header.get_best_affine()
tx_identity = ttb.TxIdentity(1)
tx_affine = ttl.TxDirectAffine(affine, name="sform")
tx_identity = TxIdentity(1)
tx_affine = TxAffine(affine, name="sform")
mri.domain.transformations += [tx_identity, tx_affine]
# Recognise invalid input format
......@@ -142,7 +143,7 @@ def centralise(timg):
TxTranslation object. """
assert isinstance(timg, tirl.timage.TImage)
centre = tuple(-pos for pos in timg.centre(weighted=False))
tx_centre = ttl.TxTranslation(*centre, name="centralise")
tx_centre = TxTranslation(*centre, name="centralise")
timg.domain.transformations.append(tx_centre)
......@@ -155,7 +156,6 @@ def _calc_rotation(a, b):
s = np.linalg.norm(axis)
c = np.sum(a * b)
angle = np.arctan2(s, c)
from tirl.transformations.linear import TxRotationAxisAngle
return TxRotationAxisAngle(angle, axis, mode="rad")
......@@ -166,7 +166,7 @@ def main(args):
mri = load_mri(args.mri, with_tx=(not args.notx))
# Set the resolution of the MRI
resolution = float(args.voxel[0])
resolution = float(args.voxel)
if resolution != 1:
from tirl.transformations.linear import TxScale
try:
......
......@@ -51,6 +51,7 @@ setup(name="tirl",
"protocols/tirl_rigid2d",
"protocols/tirl_rigid2d3d",
"protocols/tirl_rigid3d",
"protocols/tirl_slice_to_volume",
"scripts/tirlinfo",
"scripts/tirlconfig",
"scripts/slice_explorer"])
#!/usr/bin/env python3
import tempfile
import openslide
import numpy as np
from time import time
from functools import partial
import multiprocessing as mp
mp.set_start_method("fork", force=True)
import multiprocessing.dummy as mt
import tirl.settings as ts
from tirl.timage import TImage
def main():
""" Main program code. """
# Open MRI image
mri = TImage("../resources/McLaren.nii.gz")
print(mri)
del mri
# Open MRI warpfield
warp = TImage("../resources/F99_2_struct_fnirt_field.nii.gz")
print(warp)
del warp
# Open TIFF
tif = TImage("../resources/cropped_colour.tif")
print(tif)
del tif
# Open SVS
svs = TImage.fromfile("/Volumes/SSD_T5/temp/H092a.svs",
loader_kwargs={"series": 0}, storage="hdd")
print(svs)
del svs
def main2():
memlimit = 2048 # MiB
n_cpu = 4
objname = "/Volumes/SSD_T5/temp/H092a.svs"
obj = openslide.OpenSlide(objname)
dtype = np.uint8
x, y = obj.level_dimensions[0]
shape = (y, x, 3)
pixelsize = np.dtype(dtype).itemsize * 3
chunkpixels = memlimit * 1024 ** 2 // n_cpu // pixelsize
xsize = x
ysize = chunkpixels // x
if ysize == 0:
xsize = ysize = np.sqrt(chunkpixels).astype(np.int)
jobs = create_parallel_jobs(obj, xsize, ysize)
# fileno, fname = tempfile.mkstemp(prefix="svs_", dir=ts.TWD)
fname = "/Volumes/SSD_T5/temp/svs__tu6kzm5"
# print("Creating storage...")
# st = time()
# arr = np.memmap(fname, dtype=dtype, mode="r+",
# offset=0, shape=shape, order="C")
# print("Storage created: {} sec".format(time() - st))
print("Reading SVS file...")
st = time()
with mp.Pool(processes=2) as p:
worker_func = partial(worker, objname, fname, dtype, shape)
p.map(worker_func, jobs)
print("Reading complete: {} sec".format(time() - st))
def worker(objname, fname, dtype, shape, job):
x, y, xsize, ysize = job
print(f"Worker {mp.current_process().pid}: "
f"reading region (x={x}:{x+xsize}, y={y}:{y+ysize})...")
obj = openslide.OpenSlide(objname)
region = obj.read_region((x, y), level=0, size=(xsize, ysize))
arr = np.memmap(fname, dtype=dtype, mode="r+",
offset=0, shape=shape, order="C")
arr[y:y + ysize, x:x + xsize, :] = np.asarray(region)[..., :3]
arr.flush()
return 0
def create_parallel_jobs(obj, xsize, ysize):
width, height = obj.level_dimensions[0]
for y in range(0, height, ysize):
for x in range(0, width, xsize):
yield x, y, xsize, ysize
def visualise():
fname = "/Volumes/SSD_T5/temp/svs__tu6kzm5"
objname = "/Volumes/SSD_T5/temp/H092a.svs"
obj = openslide.OpenSlide(objname)
dtype = np.uint8
x, y = obj.level_dimensions[0]
shape = (y, x, 3)
arr = np.memmap(fname, dtype=dtype, mode="r+",
offset=0, shape=shape, order="C")
TImage(arr).voxels[50000:52000, 63000:65000].preview()
if __name__ == "__main__":
main2()
# visualise()
import os
from pydoc import locate
from tirl import fsl
from tirl import settings as ts
from tirl import tirlfile
from tirl import tirlobject
if ts.ENABLE_VISUALISATION:
from tirl import tirlvision
......@@ -11,6 +9,7 @@ if ts.ENABLE_VISUALISATION:
__version__ = 1.0
__author__ = "Istvan N. Huszar"
def home():
"""
Returns the home directory of TIRL.
......@@ -21,7 +20,8 @@ def home():
def testimg():
from tirl.timage import TImage
img = TImage.fromfile(os.path.join(home(), "../tests/resources/testimg.png"))
impath = os.path.join(home(), "../tests/resources/testimg.png")
img = TImage.fromfile(impath, tensor_axes=(2,))
return img
......@@ -35,50 +35,26 @@ def load(fname):
:rtype: TIRLObject
"""
from tirl import tirlobject
return tirlobject.TIRLObject.load(fname)
# if not os.path.isfile(fname):
# raise FileNotFoundError("File at {} does not exist.".format(fname))
# dump = tirlfile.load(fname)
# try:
# txtype = locate(dump["type"])
# except Exception as exc:
# raise TypeError("Unsupported type: {}".format(dump["type"])) from exc
# try:
# return txtype.load(fname)
# except (NotImplementedError, AttributeError):
# raise NotImplementedError("Class-specific load method is not "
# "implemented by the requested "
# "object ({}).".format(dump["type"]))
def expose_package_contents(baseclass, pkg, path, globals=None):
from importlib import import_module
from inspect import isclass
from glob import glob
module_name_pattern = "{}/*.py".format(path[0])
modules = [os.path.split(m)[-1] for m in glob(module_name_pattern)]
for m in modules:
if m.startswith("__"):
continue
m = "." + os.path.splitext(m)[0]
imported_module = import_module(m, package=pkg)
for class_name in dir(imported_module):
if not class_name.startswith("__"):
c = getattr(imported_module, class_name)
if isclass(c) and issubclass(c, baseclass):
if globals:
globals.update({class_name: c})
else:
globals().update({class_name: c})
# from tirl import transformations
# from tirl.domain import Domain
# from tirl import tfield
# from tirl.tfield import TField
# from tirl import timage
# from tirl.timage import TImage
# from tirl import operations
# from tirl import costs
# from tirl import regularisers
# from tirl import optimisers
# def expose_package_contents(baseclass, pkg, path, globals=None):
# from importlib import import_module
# from inspect import isclass
# from glob import glob
# module_name_pattern = "{}/*.py".format(path[0])
# modules = [os.path.split(m)[-1] for m in glob(module_name_pattern)]
# for m in modules:
# if m.startswith("__"):
# continue
# m = "." + os.path.splitext(m)[0]
# imported_module = import_module(m, package=pkg)
# for class_name in dir(imported_module):
# if not class_name.startswith("__"):
# c = getattr(imported_module, class_name)
# if isclass(c) and issubclass(c, baseclass):
# if globals:
# globals.update({class_name: c})
# else:
# globals().update({class_name: c})
......@@ -17,16 +17,13 @@
# A uniform interface object for ndarrays and memory-mapped arrays with
# a copy method, a destructor, and a signature.
import re
import os
import shutil
import joblib
import warnings
import tempfile
import numpy as np
from tirl import settings as ts
from tirl import utils as tu
class Buffer(object):
......
This diff is collapsed.
from tirl import expose_package_contents
# from tirl import expose_package_contents
from tirl.costs.cost import Cost
# Expose all cost function objects at the module level
# Note: this routine imports the Cost base class, and every
# subclass thereof from all modules within the "costs" package.
expose_package_contents(
baseclass=Cost, pkg="tirl.costs", path=__path__, globals=globals())
#
# # Expose all cost function objects at the module level
# # Note: this routine imports the Cost base class, and every
# # subclass thereof from all modules within the "costs" package.
#
# expose_package_contents(
# baseclass=Cost, pkg="tirl.costs", path=__path__, globals=globals())
......@@ -19,19 +19,14 @@
# DEPENDENCIES
from time import time
import warnings
import numpy as np
from functools import partial
import multiprocessing.dummy as mt
from itertools import product, permutations
from scipy.ndimage.interpolation import shift
from scipy.ndimage.filters import gaussian_filter
from tirl.tfield import TField
from tirl import settings as ts
from tirl.operations.robust import Robust
from tirl.costs.cost import DEFAULTS
from tirl.operations.tensor import TensorOperator
from tirl.operations.spatial import SpatialOperator
......@@ -270,7 +265,7 @@ class CostMIND(CostSSD):
mind = np.exp(-mind / np.mean(mind, ax, keepdims=True))
mind = mind / np.max(mind, ax, keepdims=True)
mind[~np.isfinite(mind)] = 0
mind[~np.isfinite(mind)] = 1
# Return MIND as a TImage
if img.order == ts.TENSOR_MAJOR:
......
......@@ -95,10 +95,12 @@ class CostSSD(Cost):
:rtype: float
"""
if self.metaparameters.get("normalise", DEFAULT_NORMALISE):
ip = self.source.interpolator.copy()
ip._values = None
ip.tensor_axes = ()
# ip = self.source.interpolator.copy()
ip = self.source.interpolator.__class__
# ip._values = None
# ip.tensor_axes = ()
ones = TField(self.source.domain, order=self.source.order,
dtype=self.source.dtype, interpolator=ip)
ones[...] = 1
......
......@@ -29,7 +29,6 @@ import warnings
import tempfile
import numpy as np
from math import ceil
from pydoc import locate
from attrdict import AttrDict
from operator import mul, add
import multiprocessing.dummy as mt
......@@ -43,9 +42,11 @@ from tirl.tirlobject import TIRLObject
from tirl.parallelism import map_chunk_coordinates
from tirl.parallelism import copy_chunk_coordinates
from tirl.transformations import Transformation, TransformationGroup
from tirl.transformations.basic import TxIdentity, TxEmbed, TxReduce
from tirl.transformations.linear import TxTranslation, TxScale
from tirl.transformations.auxiliary import TxRotationField
from tirl.transformations.basic.identity import TxIdentity
from tirl.transformations.basic.embedding import TxEmbed, TxReduce
from tirl.transformations.linear.scale import TxScale
from tirl.transformations.linear.translation import TxTranslation
from tirl.transformations.auxiliary.rotfield import TxRotationField
from tirl.parallelism import calculate_chunk_voxel_coordinates
......
#!/usr/bin/env python
# Tensor Image Registration Library
# FSL Compatibility layer
# Author: Istvan N. Huszar <istvan.huszar@dtc.ox.ac.uk>
# DEPENDENCIES
import os
# IMPLEMENTATION
def load(src):
"""
Universal load method for files/folders created by FSL tools.
:param src: source file/folder
:type src: str
:returns: equivalent TIRL object
:rtype: TIRLObject
"""
obj = None
# Load a single file
if os.path.isfile(src):
if src.endswith(".mat"):
obj = load_flirt_matrix(src)
elif src.endswith(("_field.nii.gz", "_field.nii")):
obj = load_fnirt_field(src)
elif src.endswith(("_coef.nii.gz", "_coef.nii")):
obj = load_fnirt_coef(src