Rescale edge-aligner to enable registering images of different scales

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,14 +1113,19 @@ 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,),
                                          tile_image.dtype)
                     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)

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,43 +237,81 @@ 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:
             print('Cycle %d:' % cycle)
             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.

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

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))