Added apply_xfm support, changed coordinate system to mm (instead of pixels),...

Added apply_xfm support, changed coordinate system to mm (instead of pixels), added infrastructure for changing resolution

slider/chart_reg.py

@@ -14,19 +14,20 @@ from matplotlib import patches
 import numpy as np
 import matplotlib.pyplot as plt
 
-
-def register_chart_to_slide(chart, slide, out, rlevel=0, boundary_key='outline'):
+def register_chart_to_slide(chart, slide, slide_res, out, rlevel=4, boundary_key='outline'):
 
     # load chart
     contour, cells = neurolucida.read(chart)
 
+    slide_res = slide_res * (2**rlevel)
+
     # load slide
     jp2 = glymur.Jp2k(slide)
     img = jp2.read(rlevel=rlevel)
 
     # initial scaling based on boundng boxes
     edge_crds = contour[boundary_key][:, :2] * [1, -1]
-    init = _init_scale(img, edge_crds, verbose=True, plots=False)
+    init = _init_scale(img, slide_res, edge_crds, verbose=True, plots=False)
     print(init)
     print(
         f'Rotation: {init.rotation}, Translation: {init.translation}, Scale: {init.scale}, Shear: {init.shear}')
@@ -36,7 +37,7 @@ def register_chart_to_slide(chart, slide, out, rlevel=0, boundary_key='outline')
     init_nrmls = normal(xfm_edge_coords, plot=False, distance=1, verbose=True)
 
     # refine edge_coords (to image)
-    refined_edge_coords = _refine_edge_coord(img, xfm_edge_coords, init_nrmls)
+    refined_edge_coords = _refine_edge_coord(img, slide_res, xfm_edge_coords, init_nrmls)
 
     # estimate opimised affine transform
     opt = transform.AffineTransform()
@@ -65,7 +66,7 @@ def register_chart_to_slide(chart, slide, out, rlevel=0, boundary_key='outline')
             json.dump(cells_xfm, fp, indent=4)
 
 
-def _refine_edge_coord(img, edge_coords, normals):
+def _refine_edge_coord(img, img_res, edge_coords, normals):
     """
     Refine edge_coord by sampling binarised image along normal (to edge) and looking for big step change.
     """
@@ -74,15 +75,14 @@ def _refine_edge_coord(img, edge_coords, normals):
     edge_x, edge_y = np.round(edge_coords.T).astype(int)
     nrml_x, nrml_y = normals.T
 
-    line_smpl = np.arange(-20, 20)
+    line_smpl = np.arange(-20, 20) * img_res
+
+    line_x = edge_x[:, np.newaxis] + nrml_x[:, np.newaxis] * line_smpl
+    line_y = edge_y[:, np.newaxis] + nrml_y[:, np.newaxis] * line_smpl
 
-    line_x = np.round(
-        edge_x[:, np.newaxis] + nrml_x[:, np.newaxis] * line_smpl
-    ).astype(int)
+    line_x = np.round(line_x / img_res).astype(int)
+    line_y = np.round(line_y / img_res).astype(int)
     
-    line_y = np.round(
-        edge_y[:, np.newaxis] + nrml_y[:, np.newaxis] * line_smpl
-    ).astype(int)
 
     # find threshold for background
     image = rgb2gray(img)
@@ -114,11 +114,11 @@ def _refine_edge_coord(img, edge_coords, normals):
         line_y[np.arange(min_idx.size), min_idx],
     ], axis=-1)
 
-    return refined_edge_coords
+    return refined_edge_coords * img_res
 
 
 
-def _img_props(img, verbose=False, plots=False):
+def _img_props(img, img_resolution, verbose=False, plots=False):
 
     # convert to grayscale
     image = rgb2gray(img)
@@ -142,7 +142,7 @@ def _img_props(img, verbose=False, plots=False):
     ridx = np.argmax([p0.area for p0 in p])
     p = p[ridx]
 
-    bbox = p.bbox
+    bbox = np.array(p.bbox) * img_resolution
 
     centroid = (
         (bbox[0] + bbox[2])/2,
@@ -150,25 +150,29 @@ def _img_props(img, verbose=False, plots=False):
     )
 
     # TODO: make plots work again
-    if plots:
-
-        plt.imshow(cc == p.label)
-        plt.colorbar()
-
-        rect = patches.Rectangle(
-            (bbox[1], bbox[0]),
-            bbox[3]-bbox[1],
-            bbox[2]-bbox[0],
-            linewidth=1,
-            edgecolor='r',
-            facecolor='none'
-        )
-        plt.gca().add_patch(rect)
-
-        plt.plot(centroid[1], centroid[0], 'ro')
-        plt.show()
-
-    return {'bbox': bbox, 'bbox_centroid': centroid, 'shape': (bbox[2]-bbox[0], bbox[3]-bbox[1])}
+    # if plots:
+
+    #     plt.imshow(cc == p.label)
+    #     plt.colorbar()
+
+    #     rect = patches.Rectangle(
+    #         (bbox[1], bbox[0]),
+    #         bbox[3]-bbox[1],
+    #         bbox[2]-bbox[0],
+    #         linewidth=1,
+    #         edgecolor='r',
+    #         facecolor='none'
+    #     )
+    #     plt.gca().add_patch(rect)
+
+    #     plt.plot(centroid[1], centroid[0], 'ro')
+    #     plt.show()
+
+    return {
+        'bbox': bbox, 
+        'bbox_centroid': centroid, 
+        'shape': (bbox[2]-bbox[0], bbox[3]-bbox[1])
+        }
 
 
 def _pnt_props(pnts, verbose=False, plots=False):
@@ -213,9 +217,9 @@ def _pnt_props(pnts, verbose=False, plots=False):
     return {'bbox': bbox, 'bbox_centroid': centroid, 'shape': (bbox[2]-bbox[0], bbox[3]-bbox[1])}
 
 
-def _init_scale(img, crd, verbose=False, plots=False):
+def _init_scale(img, img_resolution, crd, verbose=False, plots=False):
 
-    img_p = _img_props(img)
+    img_p = _img_props(img, img_resolution)
     crd_p = _pnt_props(crd)
 
     if verbose:

slider/slide_reg.py

@@ -141,25 +141,44 @@ SNAPSHOT_EXT = "png"
 # NumPy print formatting
 np.set_printoptions(precision=4)
 
-def _load_image(p):
-
-    if p.file.lower().endswith('.jp2'): # if jpg2k
-
+def _load_jpg2k(file, resolution, resolution_level, dtype=None):
     # load jp2
     import glymur
-        from tirl.scripts.mnd.image import set_mask
-        jp2= glymur.Jp2k(p.file)
-        img = jp2.read(rlevel=p.resolution_level)
+    jp2= glymur.Jp2k(file)
+    img = jp2.read(rlevel=resolution_level)
 
     # adjust resolution by resolution_level
-        p.resolution = p.resolution * (2**p.resolution_level)
+    resolution = resolution * (2**resolution_level)
 
     # create timg
-        timg = TImage.fromarray(img, dtype=p.dtype, tensor_axes=(2,))
-        timg.resolution = p.resolution
+    timg = TImage.fromarray(img, tensor_axes=(2,), dtype=dtype)
+    timg.resolution = resolution
+    
+    return timg
+
+
+
+def _load_image(p):
+
+    if p.file.lower().endswith('.jp2'): # if jpg2k
+
+        # # load jp2
+        # import glymur
+        # from tirl.scripts.mnd.image import set_mask
+        # jp2= glymur.Jp2k(p.file)
+        # img = jp2.read(rlevel=p.resolution_level)
+
+        # # adjust resolution by resolution_level
+        # p.resolution = p.resolution * (2**p.resolution_level)
+
+        # # create timg
+        # timg = TImage.fromarray(img, dtype=p.dtype, tensor_axes=(2,))
+        # timg.resolution = p.resolution
+
+        timg = _load_jpg2k(p.file, p.resolution, p.resolution_level, p.dtype)
 
         # add mask
-        set_mask(timg, scope=globals(), **p.mask)
+        tirl.scripts.mnd.image.set_mask(timg, scope=globals(), **p.mask)
 
         # Export
         tirl.scripts.mnd.inout.export(timg, p.export, default=None)
@@ -888,3 +907,50 @@ def register_slide_to_slide(moving, moving_res, fixed, fixed_res, out, config):
 
     # Run registration script
     run(config)
+
+
+def apply_slide_xfm(moving, moving_res, fixed, fixed_res, moving_reg, fixed_reg, out, rlevel=5, inverse=False):
+
+    moving_reg = tirl.load(moving_reg)
+    fixed_reg = tirl.load(fixed_reg)
+
+    # get transformation chain from registration
+
+    if inverse:
+        raise NotImplementedError('Inverse not yet implemented')
+    else:
+        chain = fixed_reg.domain.chain
+        chain += moving_reg.domain.chain.inverse()
+
+    # update input/output resolution in transformation chain
+
+    new_res_moving = moving_res * (2**rlevel)
+    new_res_fixed = fixed_res * (2**rlevel)
+
+    chain[0].parameters[:] = [new_res_fixed, new_res_fixed]
+    chain[-1].parameters[:] = [1 / new_res_moving, 1 / new_res_moving]
+
+    # load native resolution images
+
+    moving_native = _load_jpg2k(moving, moving_res, rlevel)
+    fixed_native = _load_jpg2k(fixed, fixed_res, rlevel)
+
+    # fixed pixel coordinates in native res
+
+    v = fixed_native.domain.get_voxel_coordinates()
+
+    # Map through transformation chain
+
+    fixed_coords_in_moving_native = chain.map(v)
+
+    # Interpolate
+
+    moving_value = moving_native.interpolator(fixed_coords_in_moving_native)
+    moving_value_reshaped = moving_value.reshape((fixed_native.shape[0],fixed_native.shape[1],3))
+
+    import glymur
+    glymur.Jp2k(out, data=np.round(moving_value_reshaped).astype(moving_native.dtype))
+
+    # import imageio
+    # imageio.imsave('test_moving.png', np.round(moving_value_reshaped).astype(moving_native.dtype))
+    # imageio.imsave('test_fixed.png', fixed_native.data)
\ No newline at end of file

slider_app.py

@@ -2,7 +2,7 @@
 import argparse
 import sys
 
-from slider.slide_reg import register_slide_to_slide
+from slider.slide_reg import register_slide_to_slide, apply_slide_xfm
 from slider.chart_reg import register_chart_to_slide
 
 def add_slide_cli(subparsers):
@@ -26,7 +26,7 @@ def add_slide_cli(subparsers):
                         help="Fixed image resolution (mm)", type=float)
 
     parser.add_argument("--out", metavar="<dir>",
-                        help="Output directory", default='./slider.reg', type=str,
+                        help="Output directory", default='./slide-to-slide.reg', type=str,
                         required=False)
     parser.add_argument("--config", metavar="<config.json>",
                         help="configuration file", default=None, type=str,
@@ -47,12 +47,44 @@ def add_chart_cli(subparsers):
                         help="Chart file", type=str)
     parser.add_argument("slide", metavar="<slide>",
                         help="Fixed slide", type=str)
+    parser.add_argument("slide_res", metavar="<slide-resolution>",
+                        help="Slide image resolution (mm)", type=float)
     parser.add_argument("--out", metavar="<dir>",
-                        help="Output directory", default='./slider.reg', type=str,
+                        help="Output directory", default='./chart-to-slide.reg', type=str,
                         required=False)
     parser.set_defaults(method='chart')
 
 
+def add_applyxfm_cli(subparsers):
+    """
+    Set up applyxfm subparser instance.
+    """
+    parser = subparsers.add_parser(
+        'APPLYXFM',
+        description='Apply transform to image.',
+        formatter_class=lambda prog: argparse.HelpFormatter(prog, max_help_position=55, width=100)
+    )
+
+    parser.add_argument("moving", metavar="<moving>",
+                        help="Moving slide image", type=str)
+    parser.add_argument("moving_res", metavar="<moving-resolution>",
+                        help="Moving image resolution (mm)", type=float)
+
+    parser.add_argument("fixed", metavar="<fixed>",
+                        help="Fixed slide image", type=str)
+    parser.add_argument("fixed_res", metavar="<fixed-resolution>",
+                        help="Fixed image resolution (mm)", type=float)
+
+    parser.add_argument("moving_reg", metavar="<moving-reg>",
+                        help="Moving timg from registration", type=str)     
+    parser.add_argument("fixed_reg", metavar="<fixed-reg>",
+                        help="Fixed timg from registration", type=str)                
+
+    parser.add_argument("out", metavar="<dir>",
+                        help="Name for resampled output image", type=str)
+    parser.set_defaults(method='applyxfm')
+
+
 if __name__ == "__main__":
     """ Main program code. """
 
@@ -61,6 +93,7 @@ if __name__ == "__main__":
     subparsers = parser.add_subparsers()
     add_slide_cli(subparsers)
     add_chart_cli(subparsers)
+    add_applyxfm_cli(subparsers)
 
     # ---
 
@@ -75,5 +108,7 @@ if __name__ == "__main__":
         register_slide_to_slide(**args)
     elif method == 'chart':
         register_chart_to_slide(**args)
+    elif method == 'applyxfm':
+        apply_slide_xfm(**args)
     else:
         raise RuntimeError(f'Unknown method: {method}')