From 6842a923c1c4c388bd651afb0c1783eb0578b5f8 Mon Sep 17 00:00:00 2001 From: Yu-An Chen Date: Mon, 20 Jul 2020 12:44:33 -0400 Subject: [PATCH] Rescale edge-aligner to enable registering images of different scales Add an optional --scale argument in the CLI to allow user specify relative resolution scale for each cycle. Current implementation only downscale images and add fixed filename prefixes when there are multiple scales. --- ashlar/reg.py | 50 +++++++++++++++- ashlar/scripts/ashlar.py | 119 +++++++++++++++++++++++++++++++-------- ashlar/thumbnail.py | 3 +- ashlar/utils.py | 6 ++ 4 files changed, 150 insertions(+), 28 deletions(-) diff --git a/ashlar/reg.py b/ashlar/reg.py index d15573f1..bc2eb9b5 100644 --- a/ashlar/reg.py +++ b/ashlar/reg.py @@ -3,6 +3,7 @@ import re import xml.etree.ElementTree import io +import copy import uuid import struct import pathlib @@ -769,6 +770,43 @@ def debug(self, t1, t2, min_size=0): plt.tight_layout() +class ScaledEdgeAligner(EdgeAligner): + + def __init__(self, original_aligner, scale): + + super(ScaledEdgeAligner, self).__init__( + original_aligner.reader, channel=original_aligner.channel + ) + + self.original_aligner = original_aligner + self.scale = scale + + self.reader = copy.deepcopy(self.original_aligner.reader) + + # self.channel = self.original_aligner.channel + self.reader.thumbnail = skimage.transform.rescale( + self.original_aligner.reader.thumbnail, + self.scale, + multichannel=False, + anti_aliasing=False, + preserve_range=True, + ).astype(self.original_aligner.reader.thumbnail.dtype) + self.metadata._positions = ( + self.original_aligner.metadata.positions * self.scale + ) + self.positions = self.original_aligner.positions * self.scale + + self.lr = sklearn.linear_model.LinearRegression() + self.lr.fit(self.metadata.positions, self.positions) + self.reader.read = self._read + + def _read(self, series, c): + return skimage.transform.rescale( + self.original_aligner.reader.read(series, c), self.scale, + multichannel=False, anti_aliasing=True, preserve_range=True + ) + + class LayerAligner(object): def __init__(self, reader, reference_aligner, channel=None, max_shift=15, @@ -960,13 +998,14 @@ class Mosaic(object): def __init__( self, aligner, shape, filename_format, channels=None, - ffp_path=None, dfp_path=None, combined=False, tile_size=None, - first=False, verbose=False + ffp_path=None, dfp_path=None, scale=1., combined=False, + tile_size=None, first=False, verbose=False ): self.aligner = aligner self.shape = tuple(shape) self.filename_format = filename_format self.channels = self._sanitize_channels(channels) + self.scale = scale self.combined = combined self.tile_size = tile_size self.first = first @@ -1074,6 +1113,11 @@ def run(self, mode='write', debug=False): sys.stdout.flush() tile_image = self.aligner.reader.read(c=channel, series=tile) tile_image = self.correct_illumination(tile_image, channel) + if self.scale != 1: + tile_image = skimage.transform.rescale( + tile_image, self.scale, multichannel=False, + anti_aliasing=True, preserve_range=True + ).astype(tile_image.dtype) if debug: color_channel = node_colors[tile] rgb_image = np.zeros(tile_image.shape + (3,), @@ -1081,7 +1125,7 @@ def run(self, mode='write', debug=False): rgb_image[:,:,color_channel] = tile_image tile_image = rgb_image func = utils.pastefunc_blend if not debug else np.add - utils.paste(mosaic_image, tile_image, position, func=func) + utils.paste(mosaic_image, tile_image, position * self.scale, func=func) if debug: np.clip(mosaic_image, 0, 1, out=mosaic_image) w = int(1e6) diff --git a/ashlar/scripts/ashlar.py b/ashlar/scripts/ashlar.py index dabd192a..1d3e9af5 100644 --- a/ashlar/scripts/ashlar.py +++ b/ashlar/scripts/ashlar.py @@ -6,10 +6,11 @@ import blessed from .. import __version__ as VERSION from .. import reg -from ..reg import PlateReader, BioformatsReader +from ..reg import PlateReader, BioformatsReader, ScaledEdgeAligner from ..filepattern import FilePatternReader from ..fileseries import FileSeriesReader from ..zen import ZenReader +from ..utils import scale_shape def main(argv=sys.argv): @@ -84,6 +85,12 @@ def main(argv=sys.argv): help=('read dark field profile image from FILES; if specified must' ' be one common file for all cycles or one file for each cycle') ) + parser.add_argument( + '--scales', type=float, metavar='SCALE', nargs='*', + help=('float point numbers denote the desired output scales relative' + 'to the raw input images; must be one SCALE for all cycles or' + 'one SCALE for each cycle') + ) parser.add_argument( '--plates', default=False, action='store_true', help='enable plate mode for HTS data' @@ -145,6 +152,21 @@ def main(argv=sys.argv): if len(dfp_paths) == 1: dfp_paths = dfp_paths * len(filepaths) + scales = args.scales + if not scales: + scales = [1.] * len(filepaths) + if len(scales) == 1: + scales *= len(filepaths) + if len(scales) != len(filepaths): + print_error( + "Wrong number of scales. Must be 1, or {} (number of input files)" + " but get {}".format(len(filepaths), len(scales)) + ) + return 1 + if 0 in scales: + print_error("Any SCALE in scales ({}) must not be 0".format(scales)) + return 1 + aligner_args = {} aligner_args['channel'] = args.align_channel aligner_args['verbose'] = not args.quiet @@ -170,8 +192,8 @@ def main(argv=sys.argv): mosaic_path_format = str(output_path / args.filename_format) return process_single( filepaths, mosaic_path_format, args.flip_x, args.flip_y, - ffp_paths, dfp_paths, aligner_args, mosaic_args, args.pyramid, - args.quiet + ffp_paths, dfp_paths, scales, aligner_args, mosaic_args, + args.pyramid, args.quiet ) except ProcessingError as e: print_error(str(e)) @@ -180,21 +202,28 @@ def main(argv=sys.argv): def process_single( filepaths, mosaic_path_format, flip_x, flip_y, ffp_paths, dfp_paths, - aligner_args, mosaic_args, pyramid, quiet, plate_well=None + scales, aligner_args, mosaic_args, pyramid, quiet, plate_well=None ): - output_path_0 = format_cycle(mosaic_path_format, 0) if pyramid: - if output_path_0 != mosaic_path_format: + if format_cycle(mosaic_path_format, 0) != mosaic_path_format: raise ProcessingError( "For pyramid output, please use -f to specify an output" " filename without {cycle} or {channel} placeholders" ) + scale_set = set(scales) + if len(scale_set) > 1: + mosaic_name = pathlib.Path(mosaic_path_format).name + scale_format = 'scaled_{scale:4.2f}-' + mosaic_name + mosaic_path_format = mosaic_path_format.replace( + mosaic_name, scale_format + ) + mosaic_args = mosaic_args.copy() if pyramid: mosaic_args['combined'] = True - num_channels = 0 + num_channels = [0 for i in scale_set] if not quiet: print('Cycle 0:') @@ -208,16 +237,28 @@ def process_single( edge_aligner.run() mshape = edge_aligner.mosaic_shape mosaic_args_final = mosaic_args.copy() - mosaic_args_final['first'] = True if ffp_paths: mosaic_args_final['ffp_path'] = ffp_paths[0] if dfp_paths: mosaic_args_final['dfp_path'] = dfp_paths[0] - mosaic = reg.Mosaic( - edge_aligner, mshape, output_path_0, **mosaic_args_final - ) - mosaic.run() - num_channels += len(mosaic.channels) + + for idx, scale in enumerate(scale_set): + relative_scale = scale / scales[0] + if relative_scale <= 1: + output_path_0 = format_scale_and_cycle( + mosaic_path_format, scale=scale, cycle=0 + ) + is_first = num_channels[idx] == 0 + mosaic_args_final['first'] = is_first + mosaic_args_final['scale'] = relative_scale + mosaic = reg.Mosaic( + edge_aligner, + scale_shape(mshape, scale / scales[0]), + output_path_0, + **mosaic_args_final + ) + mosaic.run() + num_channels[idx] += len(mosaic.channels) for cycle, filepath in enumerate(filepaths[1:], 1): if not quiet: @@ -225,26 +266,52 @@ def process_single( print(' reading %s' % filepath) reader = build_reader(filepath, plate_well=plate_well) process_axis_flip(reader, flip_x, flip_y) - layer_aligner = reg.LayerAligner(reader, edge_aligner, **aligner_args) + + cycle_scale = scales[cycle] / scales[0] + if cycle_scale != 1: + scaled_edge_aligner = ScaledEdgeAligner(edge_aligner, cycle_scale) + layer_aligner = reg.LayerAligner( + reader, scaled_edge_aligner, **aligner_args + ) + else: + layer_aligner = reg.LayerAligner(reader, edge_aligner, **aligner_args) layer_aligner.run() mosaic_args_final = mosaic_args.copy() if ffp_paths: mosaic_args_final['ffp_path'] = ffp_paths[cycle] if dfp_paths: mosaic_args_final['dfp_path'] = dfp_paths[cycle] - mosaic = reg.Mosaic( - layer_aligner, mshape, format_cycle(mosaic_path_format, cycle), - **mosaic_args_final - ) - mosaic.run() - num_channels += len(mosaic.channels) + + for idx, scale in enumerate(scale_set): + relative_scale = scale / scales[cycle] + if relative_scale <= 1: + output_path_n = format_scale_and_cycle( + mosaic_path_format, scale=scale, cycle=cycle + ) + is_first = num_channels[idx] == 0 + mosaic_args_final['first'] = is_first + mosaic_args_final['scale'] = relative_scale + mosaic = reg.Mosaic( + layer_aligner, + scale_shape(mshape, scale / scales[0]), + output_path_n, + **mosaic_args_final + ) + mosaic.run() + num_channels[idx] += len(mosaic.channels) if pyramid: print("Building pyramid") - reg.build_pyramid( - output_path_0, num_channels, mshape, reader.metadata.pixel_dtype, - reader.metadata.pixel_size, mosaic_args['tile_size'], not quiet - ) + for idx, scale in enumerate(scale_set): + reg.build_pyramid( + format_scale_and_cycle(mosaic_path_format, scale, 0), + num_channels[idx], + scale_shape(mshape, scale / scales[0]), + reader.metadata.pixel_dtype, + reader.metadata.pixel_size * scales[-1] / scale, + mosaic_args['tile_size'], + not quiet + ) return 0 @@ -286,6 +353,10 @@ def format_cycle(f, cycle): return f.format(cycle=cycle, channel='{channel}') +def format_scale_and_cycle(f, scale, cycle): + return f.format(scale=scale, cycle=cycle, channel='{channel}') + + def process_axis_flip(reader, flip_x, flip_y): metadata = reader.metadata # Trigger lazy initialization. diff --git a/ashlar/thumbnail.py b/ashlar/thumbnail.py index de8c242f..ca08b3de 100644 --- a/ashlar/thumbnail.py +++ b/ashlar/thumbnail.py @@ -42,7 +42,8 @@ def calculate_cycle_offset(reader1, reader2, scale=0.05): img2 = reader2.thumbnail if img1.shape != img2.shape: padded_shape = np.array((img1.shape, img2.shape)).max(axis=0) - padded_img1, padded_img2 = np.zeros(padded_shape), np.zeros(padded_shape) + padded_img1 = np.zeros(padded_shape, dtype=img1.dtype) + padded_img2 = np.zeros(padded_shape, dtype=img2.dtype) utils.paste(padded_img1, img1, [0, 0]) utils.paste(padded_img2, img2, [0, 0]) img1 = padded_img1 diff --git a/ashlar/utils.py b/ashlar/utils.py index 4c029842..0d2003fa 100644 --- a/ashlar/utils.py +++ b/ashlar/utils.py @@ -212,3 +212,9 @@ def imsave(fname, arr, **kwargs): del kwargs["check_contrast"] import skimage.external.tifffile skimage.external.tifffile.imsave(fname, arr, **kwargs) + + +def scale_shape(shape, scale): + scaled_shape = np.array(shape).astype(np.float64) + scaled_shape *= scale + return tuple(scaled_shape.astype(np.int64)) \ No newline at end of file