Commit 09bd88ab authored by Sean Fitzgibbon's avatar Sean Fitzgibbon
Browse files

initial commit

parents
{
"header": {
"title": "TIRL MBP-to-NeuN image registration routine - adapted for histology-to-histology registration",
"description": "",
"author": "Istvan N Huszar - adapted by Istvan & Amy Howard Mar 2021 - Adapted by Sean Fitz 2021"
},
"general": {
"name": "slice-to-slice",
"system": "linux",
"loglevel": "debug",
"logfile": "./logfile.log",
"paramlogfile": "./paramlog.log",
"verbose": false,
"outputdir": "./reg",
"stages": ["rotation", "rigid", "affine", "nonlinear"],
"warnings": true,
"isotropic": true
},
"moving": {
"file": "",
"storage": "mem",
"dtype": "f4",
"resolution": 0.0064,
"mask": {
"file": null,
"normalise": true,
"function": null,
"automask": {
"thr": 0.00,
"uthr": 1.0
}
},
"preview": false,
"export": true,
"snapshot": true
},
"fixed": {
"file": "",
"storage": "mem",
"dtype": "f4",
"resolution": 0.0064,
"mask": {
"file": null,
"normalise": true,
"function": null,
"automask": {
"thr": 0.0,
"uthr": 1.0
}
},
"preview": false,
"export": false,
"snapshot": true
},
"preprocessing": {
"moving": ["moving_preprocessing", "match_fixed_resolution"],
"fixed": ["fixed_preprocessing"]
},
"regparams": {
"init": {
"scale": {
"x0": 1.0,
"lb": 0.95,
"ub": 1.05
},
"rotation": {
"x0": 0.0,
"mode": "deg",
"lb": -10,
"ub": 10
},
"translation": {
"x0": [0.0, 0.0],
"lb": [-5.0, -5.0],
"ub": [5.0, 5.0]
},
"affine": {
"x0": [1.0, 0.0, 0.0, 0.0, 1.0, 0.0],
"lb": [0.95, -0.2, -1.0, -0.2, 0.95, -1.0],
"ub": [1.05, 0.2, 1.0, 0.2, 1.05, 1.0]
},
"nonlinear": {
"x0": 0.0,
"lb": null,
"ub": null
}
},
"rotsearch": {
"coarse": 1,
"scale": 0.125,
"visualise": false,
"xtol_rel": 0.01,
"xtol_abs": [0.001, 0.001, 0.001, 0.001],
"opt_step": 0.5
},
"rigid": {
"scaling": [160, 80, 40, 20, 10, 5],
"smoothing": [0, 0, 0, 0, 0],
"xtol_rel": 0.01,
"xtol_abs": [0.001, 0.001, 0.001, 0.001],
"opt_step": 0.1,
"visualise": false
},
"affine": {
"scaling": [160, 80, 40, 20, 10, 5],
"smoothing": [0, 0, 0],
"xtol_rel": 0.01,
"xtol_abs": [0.001, 0.001, 0.001, 0.001, 0.001, 0.001],
"opt_step": 0.1,
"visualise": false
},
"nonlinear": {
"scaling1": [8, 4, 2, 1],
"scaling": [160, 80, 40, 20, 10, 5],
"smoothing": [0, 0, 0, 0, 0, 0],
"sigma": 1,
"truncate": 1.5,
"regweight": 0.4,
"regweight1": 0.6,
"maxiter": [20, 20, 20, 20, 10, 5],
"xtol_abs": 0.1,
"xtol_rel": 0.01,
"visualise": false
}
}
}
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# _______ _____ _____ _
# |__ __| |_ _| | __ \ | |
# | | | | | |__) | | |
# | | | | | _ / | |
# | | _| |_ | | \ \ | |____
# |_| |_____| |_| \_\ |______|
#
# Copyright (C) 2018-2020 University of Oxford
# Part of the FMRIB Software Library (FSL)
# Author: Istvan N. Huszar
# Date: 7 June 2020
# FMRIB Software Library, Release 6.0.4 (c) 2020, The University of Oxford
# (the "Software")
#
# The Software remains the property of the University of Oxford ("the
# University").
#
# The Software is distributed "AS IS" under this Licence solely for
# non-commercial use in the hope that it will be useful, but in order
# that the University as a charitable foundation protects its assets for
# the benefit of its educational and research purposes, the University
# makes clear that no condition is made or to be implied, nor is any
# warranty given or to be implied, as to the accuracy of the Software,
# or that it will be suitable for any particular purpose or for use
# under any specific conditions. Furthermore, the University disclaims
# all responsibility for the use which is made of the Software. It
# further disclaims any liability for the outcomes arising from using
# the Software.
#
# The Licensee agrees to indemnify the University and hold the
# University harmless from and against any and all claims, damages and
# liabilities asserted by third parties (including claims for
# negligence) which arise directly or indirectly from the use of the
# Software or the sale of any products based on the Software.
#
# No part of the Software may be reproduced, modified, transmitted or
# transferred in any form or by any means, electronic or mechanical,
# without the express permission of the University. The permission of
# the University is not required if the said reproduction, modification,
# transmission or transference is done without financial return, the
# conditions of this Licence are imposed upon the receiver of the
# product, and all original and amended source code is included in any
# transmitted product. You may be held legally responsible for any
# copyright infringement that is caused or encouraged by your failure to
# abide by these terms and conditions.
#
# You are not permitted under this Licence to use this Software
# commercially. Use for which any financial return is received shall be
# defined as commercial use, and includes (1) integration of all or part
# of the source code or the Software into a product for sale or license
# by or on behalf of Licensee to third parties or (2) use of the
# Software or any derivative of it for research with the final aim of
# developing software products for sale or license to a third party or
# (3) use of the Software or any derivative of it for research with the
# final aim of developing non-software products for sale or license to a
# third party, or (4) use of the Software to provide any service to an
# external organisation for which payment is received. If you are
# interested in using the Software commercially, please contact Oxford
# University Innovation ("OUI"), the technology transfer company of the
# University, to negotiate a licence. Contact details are:
# Innovation@innovation.ox.ac.uk quoting reference BS/9564.
#
"""
Oxford Post-Mortem MND Biomarkers Histology-MRI Registration Pipeline, Stage 1:
Registers a histology slide to the corresponding tissue block photograph.
After the input images are centralised, the registration employs sequential
rigid, affine and non-linear alignment. Depending on the chosen direction of
the registration (default: histology->block), the following TIRL chain is
optimised:
+-------------+
| Fixed image |
+-------------+
offset:
[centralisation]-->[resampling]
chain:
-->(2D anisotropic scale)-->(2D rigid rotation)-->(2D translation)
-->(2D full affine)
-->(2D full-resolution displacement field)
+--------------+
| Moving image |
+--------------+
offset:
[centralisation]
The above chain maps the full-resolution 2D histology coordinates to the 2D
space of the tissue block photograph.
"""
__tirlscript__ = True
# DEPENDENCIES
import argparse
import json
import logging
import os
import sys
import warnings
from math import radians, degrees
import numpy as np
import tirl
import tirl.scripts.mnd.image
import tirl.scripts.mnd.inout
import tirl.settings as ts
from attrdict import AttrMap
from skimage.measure import label, regionprops
from tirl.beta import beta_function
from tirl.chain import Chain
from tirl.constants import *
from tirl.costs.mind import CostMIND
from tirl.optimisation.gnoptdiff import GNOptimiserDiffusion
from tirl.optimisation.optgroup import OptimisationGroup
from tirl.optimisation.optnl import OptNL
from tirl.regularisers.diffusion import DiffusionRegulariser
from tirl.tfield import TField
from tirl.timage import TImage
from tirl.transformations.linear.affine import TxAffine
from tirl.transformations.linear.rotation import TxRotation2D
from tirl.transformations.linear.scale import TxScale, TxIsoScale
from tirl.transformations.linear.translation import TxTranslation
from tirl.transformations.nonlinear.displacement import TxDisplacementField
# TIRL IMPORTS
# DEFINITIONS
# Rotation search: number of best initialisations to test
N_BEST = 3
# File extension of the snapshot images
SNAPSHOT_EXT = "png"
# NumPy print formatting
np.set_printoptions(precision=4)
# IMPLEMENTATION
def run(cnf=None, **options):
"""
Runs TIRL histology-to-fixed registration (stage 1).
:param cnf:
Configuration file. If no file is specified (default), the default
parameters will be used that are embedded in the source code file
(see Definitions section above). Instead of a file, a dictionary with
suitable content may also be specified.
:type cnf: Union[str, dict, None]
:param options:
Overriding configuration parameters.
:type options: Any
"""
# Load script configuration
if cnf is not None:
if isinstance(cnf, dict):
cnf = dict(cnf)
elif isinstance(cnf, str):
with open(cnf, "r") as configfile:
cnf = dict(json.load(configfile))
else:
raise TypeError(
f"Unrecognised configuration format: {cnf.__.class__.__name__}")
cnf.update(options)
p, logger = tirl.scripts.mnd.general.initialise_script(**cnf)
p.logger = logger.name # avoid globals
# Load and configure input images
if p.moving.export is True:
ext = ts.EXTENSIONS["TImage"]
p.moving.export = \
os.path.join(p.general.outputdir, f"moving.{ext}")
moving = tirl.scripts.mnd.inout.load_image(scope=globals(), **p.moving)
if p.fixed.export is True:
ext = ts.EXTENSIONS["TImage"]
p.fixed.export = os.path.join(p.general.outputdir, f"fixed.{ext}")
fixed = tirl.scripts.mnd.inout.load_image(scope=globals(), **p.fixed)
# Having loaded both images, perform actions on the histology image prior
# to registration, unless it was loaded from a TImage.
# Actions are user-defined functions within the TIRL namespace. Actions
# can be chain-loaded to perform preparatory analysis steps on the images
# before registration begins.
isalternative = p.moving.file.lower().endswith(
(ts.EXTENSIONS["TImage"], ts.EXTENSIONS["TIRLObject"]))
if not isalternative:
moving = tirl.scripts.mnd.image.perform_image_operations(
moving, *p.preprocessing.moving, scope=globals(), other=fixed,
cnf=p)
# Initialise registration frame
moving.centralise(weighted=True)
# Perform actions on the fixed image prior to registration, unless it was
# loaded from a TImage file.
isalternative = p.fixed.file.lower().endswith(
(ts.EXTENSIONS["TImage"], ts.EXTENSIONS["TIRLObject"]))
if not isalternative:
fixed = tirl.scripts.mnd.image.perform_image_operations(
fixed, *p.preprocessing.fixed, scope=globals(), other=moving, cnf=p)
fixed.centralise(weighted=True)
# Run the registration routine
try:
register(fixed, moving, p)
except Exception as exc:
logger.error(exc.args)
logger.fatal(f"The registration terminated with an exception.")
raise exc
else:
logger.fatal("The registration was completed successfully.")
def labkmeans(histo, **kwargs):
"""
Segments the foreground (tissue) in a histological slide and sets it as the
TImage mask.
This method was tested on slides that were stained with DAB + haematoxylin,
and scanned with a bright white background. The method was also effectively
removing the shadow effect close to the slide edge.
"""
from skimage.color import rgb2lab
from sklearn.cluster import KMeans
from skimage.exposure import rescale_intensity
orig_order = histo.order
histo.order = VOXEL_MAJOR
imdata = np.asarray(
rescale_intensity(histo.data[..., :3], out_range=np.uint8),
dtype=np.uint8)
lab_a = rgb2lab(imdata)[..., 1]
X = lab_a.reshape(-1, 1)
km = KMeans(n_clusters=2, random_state=0).fit(X)
kc = km.cluster_centers_
mask = km.labels_.reshape(histo.vshape)
if kc[0] > kc[1]: # make sure that the lower intensity is labeled 0
mask = 1 - mask
histo.order = orig_order
return mask
def initialise_transformations(fixed, p):
"""
Create transformation chain that will be optimised.
"""
q = p.regparams
# Scale
lb = np.asarray(q.init.scale.lb)
ub = np.asarray(q.init.scale.ub)
if p.general.isotropic:
bounds = (float(lb), float(ub))
tx_scale = TxIsoScale(
float(q.init.scale.x0), dim=2, bounds=bounds, name="scale")
else:
tx_scale = TxScale(*q.init.scale.x0, bounds=(lb, ub), name="scale")
# Rotation
if str(q.init.rotation.mode).lower() == "deg":
lb = radians(float(q.init.rotation.lb))
ub = radians(float(q.init.rotation.ub))
else:
lb = float(q.init.rotation.lb)
ub = float(q.init.rotation.ub)
tx_rotation = TxRotation2D(
float(q.init.rotation.x0), mode=q.init.rotation.mode,
bounds=(lb, ub), name="rotation")
# Translation
lb = np.asarray(q.init.translation.lb)
ub = np.asarray(q.init.translation.ub)
tx_trans = TxTranslation(
*q.init.translation.x0, bounds=(lb, ub), name="translation")
# Affine
x0 = np.asarray(q.init.affine.x0).reshape((2, 3))
lb = np.asarray(q.init.affine.lb)
ub = np.asarray(q.init.affine.ub)
tx_affine = TxAffine(x0, bounds=(lb, ub), name="affine")
# Append linear transformations to the domain of the fixed image
linear_chain = Chain(tx_rotation, tx_scale, tx_trans, tx_affine)
domain = fixed.domain[:]
domain.chain.extend(linear_chain)
# Nonlinear
x0 = float(q.init.nonlinear.x0) * np.ones((2, *fixed.vshape))
if q.init.nonlinear.lb is None:
lb = None
else:
lb = float(q.init.nonlinear.lb) * np.ones_like(x0)
if q.init.nonlinear.ub is None:
ub = None
else:
ub = float(q.init.nonlinear.ub) * np.ones_like(x0)
field = TField.fromarray(
x0, tensor_axes=(0,), copy=False, domain=domain[:], order=TENSOR_MAJOR)
tx_nonlinear = TxDisplacementField(
field, bounds=(lb, ub), name="nonlinear", mode=NL_REL)
# Return the full transformation chain
return Chain(*linear_chain, tx_nonlinear)
def register(fixed, moving, cnf):
"""
Runs the four registration stages: rotation search, rigid, affine, and
non-linear. The function has no return value. The transformation chain is
attached to the Domain of the fixed TImage object and is optimised in situ.
:param fixed:
fixed image (to which the chain is attached)
:type fixed: TImage
:param moving:
moving image (that defines the coordinate space that the mapping
is into)
:type moving: TImage
:param cnf: all configuration options
:type cnf: dict or AttrMap
"""
p = AttrMap(cnf)
q = p.regparams
logger = logging.getLogger(p.logger)
# Create transformation chain (does not change the fixed image)
# rotation -> scale -> translation -> affine -> nonlinear
logger.info("Initialising transformation chain...")
chain = initialise_transformations(fixed, p)
logger.info("Transformation chain has been initialised.")
# Generate output: initial alignment
fixed.save(os.path.join(
p.general.outputdir, "fixed.timg"), overwrite=True)
fixed.snapshot(os.path.join(
p.general.outputdir, f"fixed.{SNAPSHOT_EXT}"), overwrite=True)
moving.save(os.path.join(
p.general.outputdir, "moving.timg"), overwrite=True)
moving.snapshot(os.path.join(
p.general.outputdir, f"moving.{SNAPSHOT_EXT}"), overwrite=True)
# Set the first part of the chain
fixed.domain.chain.extend(chain[:-2])
# Rotation search
if "rotation" in p.general.stages:
logger.info("Starting rotation search...")
rotation_search2d(fixed, moving, p)
logger.info("Completed rotation search.")
# Generate output
fixed.save(os.path.join(
p.general.outputdir, "fixed1_rotation.timg"), overwrite=True)
moving.evaluate(fixed.domain).snapshot(os.path.join(
p.general.outputdir, f"moving1_rotation.{SNAPSHOT_EXT}"),
overwrite=True)
else:
logger.info("Rotation search was skipped.")
# Rigid registration
if "rigid" in p.general.stages:
logger.info("Starting rigid registration...")
rigid2d(fixed, moving, p)
logger.info("Completed rigid registration.")
# Generate output
fixed.save(os.path.join(
p.general.outputdir, "fixed2_rigid.timg"), overwrite=True)
moving.evaluate(fixed.domain).snapshot(os.path.join(
p.general.outputdir, f"moving2_rigid.{SNAPSHOT_EXT}"),
overwrite=True)
else:
logger.info("Rigid registration was skipped.")
# Affine registration
fixed.domain.chain.append(chain[-2])
if "affine" in p.general.stages:
logger.info("Starting affine registration...")
affine2d(fixed, moving, p)
logger.info("Completed affine registration.")
# Generate output
fixed.save(os.path.join(
p.general.outputdir, "fixed3_affine.timg"), overwrite=True)
moving.evaluate(fixed.domain).snapshot(os.path.join(
p.general.outputdir, f"moving3_affine.{SNAPSHOT_EXT}"),
overwrite=True)
else:
logger.info("Affine registration was skipped.")
# Non-linear registration
tx_nonlinear = chain[-1]
tx_nonlinear.domain.chain = fixed.domain.chain[:]
fixed.domain.chain.append(tx_nonlinear)
if "nonlinear" in p.general.stages:
logger.info("Starting non-linear registration...")
diffreg2d(fixed, moving, p)
logger.info("Completed non-linear registration.")
# Generate output
fixed.save(os.path.join(
p.general.outputdir, "fixed4_nonlinear.timg"), overwrite=True)
moving.evaluate(fixed.domain).snapshot(os.path.join(
p.general.outputdir, f"moving4_nonlinear.{SNAPSHOT_EXT}"),
overwrite=True)
else:
logger.info("Non-linear registration was skipped.")
def rotation_search2d(fixed, moving, cnf):
"""
Finds the best relative orientation of the images.
:param fixed: Fixed image.
:type fixed: TImage
:param moving: Moving image.
:type moving: TImage
:param cnf: all configuration options
:type cnf: dict or AttrMap
"""
p = AttrMap(cnf)
q = p.regparams.rotsearch
logger = logging.getLogger(p.logger)
# Part 1: coarse search for N best orientations
# Coarse search at a predefined scale
fixed.resample(float(q.scale), copy=False)
moving.resample(float(q.scale), copy=False)
step = radians(q.coarse)
tx_rotation = fixed.domain.chain["rotation"]
lb, ub = tx_rotation.parameters.get_bounds()[0]
rotations = np.arange(lb, ub, step)
costvals = []
for i, angle in enumerate(rotations):
tx_rotation.parameters.parameters[0] = angle
tx_rotation.parameters.set_lower_bounds(angle - step / 2)
tx_rotation.parameters.set_upper_bounds(angle + step / 2)
cost = CostMIND(moving, fixed, normalise=True, kernel=MK_FULL)()
# cost = CostMI(moving, fixed, normalise=True, bins=32)()
logger.debug(f"{degrees(angle)} deg: {cost}")
costvals.append([cost, angle])
else:
costvals = np.asarray(costvals)
# Test the best N initial rotations
n_best = max(min(N_BEST, len(rotations)), 1)
best_angles = costvals[np.argsort(costvals[:, 0]), 1][:n_best].ravel()
logger.info(f"The {n_best} best initialisation angles: "
f"{np.rad2deg(best_angles)} deg")
# Part 2: fine-tune the rotation parameter by simultaneous additional
# optimisation of translation and scaling, starting from the best three
# orientations.
tx_scale = fixed.domain.chain["scale"]
tx_translation = fixed.domain.chain["translation"]
og = OptimisationGroup(tx_rotation, tx_scale, tx_translation)
scale_orig = tx_scale.parameters.parameters.copy()
translation_orig = tx_translation.parameters.parameters.copy()
costvals = []
for angle in best_angles:
tx_rotation.parameters.parameters[0] = angle
tx_rotation.parameters.set_lower_bounds(angle - step / 2)
tx_rotation.parameters.set_upper_bounds(angle + step / 2)
tx_rotation.parameters.unlock()
tx_scale.parameters.parameters[:] = scale_orig.copy()
tx_translation.parameters.parameters[:] = translation_orig.copy()
# Rescale to the same predefined resolution as above
fixed.resample(float(q.scale), copy=False)
moving.resample(float(q.scale), copy=False)
# Set cost function
cost = CostMIND(moving, fixed, maskmode="and", normalise=True)
# cost = CostMI(moving, fixed, maskmode="and", normalise=True)
# Start optimisation
logger.info(f"Co-optimising scale and translation at "
f"{degrees(angle)} deg...")