Commit fa6d0484 authored by inhuszar's avatar inhuszar
Browse files

Updated histology_to_block with the library functions.

parent 17f40a8c
......@@ -80,5 +80,6 @@ setup(name="tirl",
"tirl/transformations/linear",
"tirl/transformations/nonlinear",
"tirl/usermodules"],
scripts=["tirl/tirl"]
scripts=["tirl/tirl",
"tirl/tirlvision/sliceview.py"]
)
No preview for this file type
......@@ -29,6 +29,6 @@ The package includes the following modules:
"""
from tirl.scripts.mnd import io
from tirl.scripts.mnd import inout
from tirl.scripts.mnd import image
from tirl.scripts.mnd import general
\ No newline at end of file
......@@ -5,47 +5,54 @@
"author": "Istvan N Huszar"
},
"general": {
"name": "histology_to_block",
"direction": "h2b",
"system": "linux",
"verbosity": "debug",
"logfile": null,
"outdir": "/Users/inhuszar/temp/histo_to_block/4La/reg_svs_automask",
"outputdir": "/mnt/rio/temp/histreg/histo_to_block/4La/reg",
"stages": ["rotation", "rigid", "affine", "nonlinear"],
"warnings": false
},
"histology": {
"file": "/Users/inhuszar/temp/histo_to_block/4La/NP140-14-4La-2-PLP.svs",
"alternative": null,
"file": "/mnt/rio/temp/histreg/histo_to_block/4La/NP140-14-4La-2-PLP.svs",
"storage": "mem",
"dtype": "<f4",
"resolution": 0.008,
"mask": {
"file": null,
"function": null,
"normalise": false
},
"automask": {
"thr": 0,
"uthr": 1
"normalise": true,
"function": "labkmeans",
"automask": {
"thr": 0,
"uthr": 1
}
},
"resolution": 0.008,
"preview": false,
"actions": ["resample_histology", "labkmeans", "histology_preproc"]
"export": true,
"snapshot": true
},
"block": {
"file": "/Users/inhuszar/temp/histo_to_block/4La/block_4La_sup_seg.tif",
"file": "/mnt/rio/temp/histreg/histo_to_block/4La/block_4La_sup_seg.tif",
"storage": "mem",
"dtype": "<f4",
"resolution": 0.050,
"mask": {
"file": null,
"normalise": true,
"function": null,
"normalise": false
},
"automask": {
"thr": 0,
"uthr": 1
"automask": {
"thr": 0,
"uthr": 1
}
},
"resolution": 0.050,
"preview": false,
"actions": ["block_preproc"]
"export": false,
"snapshot": false
},
"preprocessing": {
"histology": ["match_block_resolution", "histology_preprocessing"],
"block": ["block_preprocessing"]
},
"regparams": {
"init": {
......@@ -96,7 +103,7 @@
"scaling": [20, 10, 5, 2],
"smoothing": [0, 0, 0, 0],
"xtol_rel": 0.01,
"xtol_abs": [0.001, 0.001, 0.1, 0.001, 0.001, 0.1],
"xtol_abs": [0.001, 0.1, 0.1, 0.1, 0.001, 0.1],
"opt_step": [0.01, 0.01, 1, 0.01, 0.01, 1],
"visualise": false
},
......
......@@ -159,14 +159,14 @@ def run_stage(target, source, stage_no, counter, scope, config):
# Export stage result as TImage
ext = ts.EXTENSIONS["TImage"]
scr.io.export(target, q.export.timage, default=os.path.join(
scr.inout.export(target, q.export.timage, default=os.path.join(
p.general.outputdir, f"{counter}_stage{stage_no}.{ext}"))
# Create snapshot showing end-stage alignment
snpfile = os.path.join(
p.general.outputdir, f"{counter}_stage{stage_no}.{SNAPSHOT_EXT}")
scr.io.snapshot(source.evaluate(target.domain), q.export.snapshot,
default=snpfile)
scr.inout.snapshot(source.evaluate(target.domain), q.export.snapshot,
default=snpfile)
# Save mask if requested
if target.tmask() is not None:
......@@ -175,14 +175,14 @@ def run_stage(target, source, stage_no, counter, scope, config):
targetmaskfile = os.path.join(
p.general.outputdir,
f"{counter}_stage{stage_no}_targetmask.{SNAPSHOT_EXT}")
scr.io.snapshot(mask, q.export.target_mask, default=targetmaskfile)
scr.inout.snapshot(mask, q.export.target_mask, default=targetmaskfile)
if source.tmask() is not None:
from tirl.timage import TImage
mask = TImage.fromTField(target.tmask())
sourcemaskfile = os.path.join(
p.general.outputdir,
f"{counter}_stage{stage_no}_sourcemask.{SNAPSHOT_EXT}")
scr.io.snapshot(mask, q.export.source_mask, default=sourcemaskfile)
scr.inout.snapshot(mask, q.export.source_mask, default=sourcemaskfile)
# Visualise end-stage alignment
if q.visualise:
......
......@@ -59,17 +59,17 @@ import logging
import argparse
import warnings
import numpy as np
from pydoc import locate
from attrdict import AttrMap
from math import radians, degrees
from skimage.morphology import dilation
from skimage.color import rgb2yiq
from skimage.measure import label, regionprops
# TIRL IMPORTS
import tirl
import tirl.settings as ts
import tirl.scripts.mnd as scr
from tirl.chain import Chain
from tirl.tfield import TField
from tirl.timage import TImage
......@@ -92,12 +92,14 @@ from tirl.transformations.nonlinear.displacement import TxDisplacementField
from tirl.constants import *
# Resolution of the histological image by magnification
HISTO_RESOLUTIONS = {20: 5.019e-4}
# 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)
logger = logging.getLogger("tirl.scripts.mnd.histology_to_block")
p = AttrMap({})
# IMPLEMENTATION
......@@ -117,66 +119,58 @@ def run(cnf=None, **options):
:type options: Any
"""
# Prepare for running stage 1, obtain script attribute map (p).
p = housekeeping(cnf, **options)
# Load histology slide or alternative image (with pre-optimised chain)
# and perform initial actions prior to registration.
# Also centralise histo image, unless it is from an alternative source.
if p.histology.alternative is None:
alternative = False
histo = load_histology_slide(p.histology)
logger.info(f"Loaded histology image: {histo}")
else:
alternative = True
histo = load_alternative_input(p.histology.alternative)
logger.info(f"Loaded histology image from alternative source: {histo}")
# Preview histology image if requested
if p.histology.preview:
histo.preview()
# Load block photograph
block = load_block_photo(p.block)
# Report the end of initialisation of the block photograph.
logger.info(f"Loaded block photograph: {block}")
# Preview block photograph if requested
if p.block.preview:
block.preview()
global p
global logger
global verbosity
# Having loaded both images, perform actions on the histology slide prior
# to registration, unless it was loaded from an alternative source.
# 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 = scr.general.initialise_script(**cnf)
verbosity = scr.general.verbosity(p.general.verbosity)
# Load and configure input images
if p.histology.export is True:
ext = ts.EXTENSIONS["TImage"]
p.histology.export = \
os.path.join(p.general.outputdir, f"histology.{ext}")
histo = scr.inout.load_histology(scope=globals(), **p.histology)
if p.block.export is True:
ext = ts.EXTENSIONS["TImage"]
p.block.export = os.path.join(p.general.outputdir, f"block.{ext}")
block = scr.inout.load_image(scope=globals(), **p.block)
# 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.
if not alternative:
histo = perform_actions(histo, *p.histology.actions, other=block)
isalternative = p.histology.file.lower().endswith(
(ts.EXTENSIONS["TImage"], ts.EXTENSIONS["TIRLObject"]))
if not isalternative:
histo = scr.image.perform_image_operations(
histo, *p.preprocessing.histology, scope=globals(), other=block)
# Initialise registration frame
histo.centralise()
# Check validity and centralise histology slide.
assert isinstance(histo, TImage)
# Perform actions on the tissue block photograph prior to registration.
block = perform_actions(block, *p.block.actions, other=histo)
# Check validity and also centralise the block photograph.
assert isinstance(block, TImage), "Block must be a TImage."
block.centralise()
# Set masks
set_histology_mask(histo, p.histology)
set_block_mask(block, p.block)
# Export mask images
hm = histo.tmask()
if hm is not None:
TImage.fromTField(hm).snapshot(os.path.join(
p.general.outdir, f"histology_mask.{SNAPSHOT_EXT}"), overwrite=True)
bm = block.tmask()
if bm is not None:
TImage.fromTField(bm).snapshot(os.path.join(
p.general.outdir, f"block_mask.{SNAPSHOT_EXT}"), overwrite=True)
# Run the selected registration stages (rot, rigid, affine, non-linear)
# Perform actions on the block image prior to registration, unless it was
# loaded from a TImage file.
isalternative = p.block.file.lower().endswith(
(ts.EXTENSIONS["TImage"], ts.EXTENSIONS["TIRLObject"]))
if not isalternative:
block = scr.image.perform_image_operations(
block, *p.preprocessing.block, scope=globals(), other=histo)
block.centralise()
# Run the registration routine
try:
if str(p.general.direction).lower() == "h2b":
register(histo, block, p.regparams)
......@@ -193,114 +187,6 @@ def run(cnf=None, **options):
logger.fatal("The registration was completed successfully.")
def housekeeping(cnf, **options):
"""
Prepares the run of stage 1. Does the following:
1. Unifies command and configuration file options, sets global params
2. Disables warnings (if requested)
3. Creates the log file and initialises the logger
4. Creates the output directory if it does not yet exist
:returns: script-wide attribute map
:rtype: AttrMap
"""
# Load pipeline configuration and make parameters global at the module level
if cnf is not None:
if isinstance(cnf, dict):
options_from_file = dict(cnf)
elif isinstance(cnf, str):
with open(cnf, "r") as cnf:
options_from_file = json.load(cnf)
else:
raise TypeError(
f"Unrecognised configuration format: {cnf.__.class__.__name__}")
options_from_file.update(options)
options = options_from_file
global p
p = AttrMap(options)
# Warnings
if not p.general.warnings:
import warnings
warnings.filterwarnings("ignore")
# Create output directory if it does not yet exist
if not os.path.isdir(p.general.outdir):
os.makedirs(p.general.outdir)
# Create logfile (defaults to the output directory)
if p.general.logfile is None:
p.general.logfile = os.path.join(
p.general.outdir, "histology_to_block.log")
# Create logfile
global logger
logging.basicConfig(format='[%(asctime)s] %(message)s',
datefmt='%d-%m-%Y %H:%M:%S',
filename=p.general.logfile,
filemode="w")
logger = logging.getLogger("tirl.protocols.histology_to_block")
logger.setLevel(verbosity(p.general.verbosity))
logger.critical(f"The program started with the command: "
f"{' '.join(sys.argv)}")
# Create a local copy of the configuration file in the output directory
with open(os.path.join(p.general.outdir, "configuration.json"), "w") as fp:
json.dump(cnf, fp, sort_keys=True, indent=4)
return p
def load_histology_slide(q):
"""
Loads histology slide from file as TImage.
:param q: configurations related to the histology input
:type q: AttrMap
"""
logger.info(f"Loading histology slide from file {q.file}...")
if not os.path.isfile(q.file):
raise FileNotFoundError(
f"Histology slide could not be loaded because the following input "
f"file does not exist: {q.file}.")
# Load slice
if q.file.endswith(".svs"):
import tirl.loader
histo = TField.fromfile(q.file, storage=MEM, dtype=q.dtype,
loader=tirl.loader.OpenSlideLoader,
loader_kwargs={"level": 2})
histo = TImage.fromTField(histo, copy=False)
else:
histo = TImage.fromfile(q.file, storage=q.storage, dtype=q.dtype)
# Set resolution based on user input
histo.resolution = q.resolution
logger.info(f"Pixel size of the histology slide was set per user "
f"specification to {q.resolution} [a.u./px].")
return histo
def set_histology_mask(histo, q):
if q.mask.file is not None:
mask = TField.fromfile(q.mask.file)
mask.reduce_tensors()
histo.mask = mask / np.max(mask.data)
# histo.apply_mask(mask.data, normalise=q.mask.normalise)
elif q.mask.function is not None:
func = globals().get(q.mask.function, locate(q.mask.function))
histo.mask = func(histo)
else:
mask = automask(histo, thr=q.automask.thr, uthr=q.automask.uthr)
histo.mask = mask
# histo.apply_mask(mask.data, normalise=False)
def labkmeans(histo, **kwargs):
"""
Segments the foreground (tissue) in a histological slide and sets it as the
......@@ -326,109 +212,13 @@ def labkmeans(histo, **kwargs):
mask = km.labels_.reshape(histo.vshape)
if kc[0] > kc[1]: # make sure that the lower intensity is labeled 0
mask = 1 - mask
histo.mask = mask
histo.order = orig_order
return histo
def set_block_mask(block, q):
if q.mask.file is not None:
mask = TField.fromfile(q.mask.file)
mask.reduce_tensors()
# block.apply_mask(mask.data, normalise=q.mask.normalise)
block.mask = mask / np.max(mask.data)
elif q.mask.function is not None:
func = globals().get(q.mask.function, locate(q.mask.function))
block.mask = func(block)
else:
mask = automask(block, thr=q.automask.thr, uthr=q.automask.uthr)
# block.apply_mask(mask.data, normalise=False)
block.mask = mask
def load_alternative_input(fpath):
"""
Loads histology slide from alternative source. Use this functionality to
load a histology slide from an existing TImage file. The following
assumptions are made about the input:
1. it is a histology slide
2. all pre-reg actions have been carried out on the image
3. the image has a suitable transformation chain to optimise
"""
return tirl.load(fpath)
def load_block_photo(q):
"""
Loads tissue block photograph from file as TImage.
:param q: configurations related to the histology input
:type q: AttrMap
"""
logger.info(f"Loading tissue block photograph from file {q.file}...")
if not os.path.isfile(q.file):
raise FileNotFoundError(
f"Tissue block photograph could not be loaded because the "
f"following input file does not exist: {q.file}.")
# Load slice
block = TImage.fromfile(q.file, storage=q.storage, dtype=q.dtype)
# Set resolution based on user input
block.resolution = q.resolution
logger.info(f"Pixel size of the tissue block photograph was set per user "
f"specification to {q.resolution} [a.u./px].")
return block
def perform_actions(img, *actions, other=None):
"""
Calls the sequence of action functions on the specified TImage.
The input image should not be modified by the action functions. The result
of each action is used as the input for the next action. The action
functions should have the following signature:
actionfunc(img: TImage, other=None, history=()) -> result: TImage
The result of the last action is returned by this function.
:param img: TImage on which the actions are performed
:type img: TImage
:param actions:
Sequence of functions from the TIRL namespace defined by their full
module path. For example, to call the 'filter' function from the user
module 'histology' in package 'pkg', enter
'tirl.usermodules.pkg.histology.filter'.
:type actions: str
:param other: other TImage that is used in the registration (optional)
:type other: TImage
:returns: pre-processed image
:rtype: TImage
"""
if not actions:
return img
action_history = []
lastimg = img
for action in actions:
logger.info(f"Performing action \"{action}\"...")
func = globals().get(action, locate(action))
res = func(lastimg, other=other, history=action_history)
lastimg = res or lastimg
action_history.append((action, lastimg))
else:
assert isinstance(lastimg, TImage)
return lastimg
# histo.mask = mask
# return histo
return mask
def initialise_chain(fixed, q):
def initialise_transformations(fixed, q):
"""
Create transformation chain that will be optimised.
......@@ -505,18 +295,18 @@ def register(fixed, moving, q):
# Create transformation chain (does not change the fixed image)
# rotation -> scale -> translation -> affine -> nonlinear
logger.info("Initialising transformation chain...")
chain = initialise_chain(fixed, q)
chain = initialise_transformations(fixed, q)
logger.info("Transformation chain has been initialised.")
# Generate output: initial alignment
fixed.save(os.path.join(
p.general.outdir, "fixed.timg"), overwrite=True)
p.general.outputdir, "fixed.timg"), overwrite=True)
fixed.snapshot(os.path.join(
p.general.outdir, f"fixed.{SNAPSHOT_EXT}"), overwrite=True)
p.general.outputdir, f"fixed.{SNAPSHOT_EXT}"), overwrite=True)
moving.save(os.path.join(
p.general.outdir, "moving.timg"), overwrite=True)
p.general.outputdir, "moving.timg"), overwrite=True)
moving.snapshot(os.path.join(
p.general.outdir, f"moving.{SNAPSHOT_EXT}"), overwrite=True)
p.general.outputdir, f"moving.{SNAPSHOT_EXT}"), overwrite=True)
# Set the first part of the chain
fixed.domain.chain.extend(chain[:-2])
......@@ -528,9 +318,9 @@ def register(fixed, moving, q):
logger.info("Completed rotation search.")
# Generate output
fixed.save(os.path.join(
p.general.outdir, "fixed1_rotation.timg"), overwrite=True)
p.general.outputdir, "fixed1_rotation.timg"), overwrite=True)
moving.evaluate(fixed.domain).snapshot(os.path.join(
p.general.outdir, f"moving1_rotation.{SNAPSHOT_EXT}"),
p.general.outputdir, f"moving1_rotation.{SNAPSHOT_EXT}"),
overwrite=True)
else:
logger.info("Rotation search was skipped.")
......@@ -542,9 +332,9 @@ def register(fixed, moving, q):
logger.info("Completed rigid registration.")
# Generate output
fixed.save(os.path.join(
p.general.outdir, "fixed2_rigid.timg"), overwrite=True)
p.general.outputdir, "fixed2_rigid.timg"), overwrite=True)
moving.evaluate(fixed.domain).snapshot(os.path.join(
p.general.outdir, f"moving2_rigid.{SNAPSHOT_EXT}"), overwrite=True)
p.general.outputdir, f"moving2_rigid.{SNAPSHOT_EXT}"), overwrite=True)
else:
logger.info("Rigid registration was skipped.")
......@@ -556,9 +346,9 @@ def register(fixed, moving, q):
logger.info("Completed affine registration.")
# Generate output
fixed.save(os.path.join(
p.general.outdir, "fixed3_affine.timg"), overwrite=True)
p.general.outputdir, "fixed3_affine.timg"), overwrite=True)
moving.evaluate(fixed.domain).snapshot(os.path.join(
p.general.outdir, f"moving3_affine.{SNAPSHOT_EXT}"), overwrite=True)
p.general.outputdir, f"moving3_affine.{SNAPSHOT_EXT}"), overwrite=True)
else:
logger.info("Affine registration was skipped.")
......@@ -572,9 +362,9 @@ def register(fixed, moving, q):
logger.info("Completed non-linear registration.")
# Generate output
fixed.save(os.path.join(
p.general.outdir, "fixed4_nonlinear.timg"), overwrite=True)
p.general.outputdir, "fixed4_nonlinear.timg"), overwrite=True)
moving.evaluate(fixed.domain).snapshot(os.path.join(
p.general.outdir, f"moving4_nonlinear.{SNAPSHOT_EXT}"),
p.general.outputdir, f"moving4_nonlinear.{SNAPSHOT_EXT}"),
overwrite=True)
else:
logger.info("Non-linear registration was skipped.")
......@@ -750,9 +540,9 @@ def diffreg2d(fixed, moving, q):
moving.resample(1, copy=False)