Added support of multiple boundary contours

slider/chart_reg.py

@@ -137,22 +137,22 @@ def register_chart_to_slide(chart, slide, slide_res, outdir, boundary_key=None,
         f"Rotation: {init.rotation}, Translation: {init.translation}, Scale: {init.scale}"
     )
 
-    # # calculate normal line to boundary points
-    # xfm_edge_coords = _apply_xfm(init, edge_crds)
-    # init_nrmls = normal(xfm_edge_coords)
-
-    # # refine edge_coords (to image)
-    # refined_edge_coords = _refine_edge_coord(
-    #     img, slide_res, xfm_edge_coords, init_nrmls
-    # )
-
-    # # estimate opimised affine transform
-    # opt = transform.SimilarityTransform()
-    # opt.estimate(edge_crds, refined_edge_coords)
-    # print(opt)
-    # # print(
-    # #     f"Rotation:\t{opt.rotation}\nTranslation:\t{opt.translation}\nScale:\t\t{list(opt.scale)}\nShear:\t\t{opt.shear}"
-    # # )
+    # calculate normal line to boundary points
+    xfm_edge_coords = apply_xfm(init, edge_crds)
+    init_nrmls = normal(xfm_edge_coords)
+
+    # refine edge_coords (to image)
+    refined_edge_coords, exclude_mask = refine_edge_coord(
+        img, slide_res, xfm_edge_coords, init_nrmls
+    )
+
+    # estimate opimised affine transform
+    opt = transform.SimilarityTransform()
+    opt.estimate(edge_crds[exclude_mask, :], refined_edge_coords[exclude_mask, :])
+    print(opt)
+    print(
+        f"Rotation: {opt.rotation}, Translation: {opt.translation}, Scale: {opt.scale}"
+    )
 
     # # save opt transform
     # np.savetxt(f"{out}/chart-to-image.xfm", opt.params)
@@ -160,11 +160,11 @@ def register_chart_to_slide(chart, slide, slide_res, outdir, boundary_key=None,
     # apply opt-xfm to contours and cells and save
 
     contour_xfm = [
-        (contour.name, apply_xfm(init, contour.points[:, :2] * [1, -1]).tolist(), contour.closed) for contour in chart.contours
+        (contour.name, apply_xfm(opt, contour.points[:, :2] * [1, -1]).tolist(), contour.closed) for contour in chart.contours
     ]
 
     if do_plots:
-        fig, ax = plt.subplots(figsize=(10, 10))
+        fig, ax = plt.subplots(figsize=(20, 30))
 
         extent = np.array([0, img.shape[1], img.shape[0], 0]) * slide_res
         ax.imshow(img, extent=extent, cmap='gray')
@@ -200,10 +200,11 @@ def register_chart_to_slide(chart, slide, slide_res, outdir, boundary_key=None,
     with open(f"{outdir}/contour.json", "w") as fp:
         json.dump(contour_xfm, fp, indent=4)
 
-    # if cells.size > 0:
-    #     cells_xfm = apply_xfm(init, cells[:, :2] * [1, -1])
-    #     with open(f"{outdir}/cells.json", "w") as fp:
-    #         json.dump(cells_xfm.tolist(), fp, indent=4)
+    if chart.n_cells > 0:
+        cells = np.concatenate([cell.point[:2][np.newaxis, :] for cell in chart.cells]) * [1, -1]
+        cells_xfm = apply_xfm(init, cells )
+        with open(f"{outdir}/cells.json", "w") as fp:
+            json.dump(cells_xfm.tolist(), fp, indent=4)
 
 
 def enhance(img0, kernel_size=None, lower_percentile=2, upper_percentile=98, sigma=5):
@@ -263,16 +264,18 @@ def refine_edge_coord(img, img_res, edge_coords, normals, do_plots=True):
 
     # TODO: move these line extents to the config file
     # line_smpl = np.linspace(-0.03, 0.15, 20)
-    line_smpl = np.linspace(-0.1, 0.15, 20)
+    line_smpl = np.linspace(-0.2, 0.2, 20)
 
     line_x = edge_x[:, np.newaxis] + nrml_x[:, np.newaxis] * line_smpl
     line_y = edge_y[:, np.newaxis] + nrml_y[:, np.newaxis] * line_smpl
 
+    brainmask = segment_foreground(img)
+
     if do_plots:
         fig, ax = plt.subplots(figsize=(20, 30))
 
         extent = np.array([0, img.shape[1], img.shape[0], 0]) * img_res
-        ax.imshow(img, extent=extent, cmap='gray')
+        ax.imshow(brainmask, extent=extent, cmap='gray')
 
         for line_x0, line_y0 in zip(line_x, line_y):
             ax.plot(line_x0, line_y0, "b.")
@@ -281,7 +284,7 @@ def refine_edge_coord(img, img_res, edge_coords, normals, do_plots=True):
 
         ax.set_xlabel("mm")
         ax.set_ylabel("mm")
-        ax.set_title("Largest connected component + normals")
+        ax.set_title("Brainmask + normals")
 
         plt.show()
         # fig.savefig(f"{OUTDIR}/normals.png", bbox_inches="tight", dpi=300)
@@ -291,10 +294,19 @@ def refine_edge_coord(img, img_res, edge_coords, normals, do_plots=True):
     line_x = np.round(line_x / img_res).astype(int)
     line_y = np.round(line_y / img_res).astype(int)
 
-    # line_int = brainmask[line_y, line_x]
-    line_int = img[line_y, line_x]
+    
+    line_int = brainmask[line_y, line_x]
+    # line_int = img[line_y, line_x]
+
+    
+
+    # exlude constant rows
+    constant_idx = np.all(line_int == line_int[:, 0][:, np.newaxis], axis=1)
+    # line_int = line_int[~constant_idx, :]
 
-    min_idx = np.argmax(np.diff(line_int, axis=-1), axis=-1)
+    # print(line_int)
+
+    min_idx = np.argmax(np.abs(np.diff(line_int, axis=-1)), axis=-1)
 
     refined_edge_coords = np.stack(
         [
@@ -304,15 +316,21 @@ def refine_edge_coord(img, img_res, edge_coords, normals, do_plots=True):
         axis=-1,
     )
 
+    print(constant_idx.shape)
+    print(refined_edge_coords.shape)
+    # refined_edge_coords = refined_edge_coords[~constant_idx, :]
+    print(refined_edge_coords.shape)
+
+
     # TODO: plot edge_coords + refined_edge_coords
     if do_plots:
         fig, ax = plt.subplots(figsize=(20, 30))
 
         extent = np.array([0, img.shape[1], img.shape[0], 0]) * img_res
-        ax.imshow(img, extent=extent, cmap='gray')
+        ax.imshow(brainmask, extent=extent, cmap='gray')
 
         ax.plot(edge_x, edge_y, "r.-")
-        ax.plot(refined_edge_coords[:, 0], refined_edge_coords[:, 1], "g.-")
+        ax.plot(refined_edge_coords[~constant_idx, 0], refined_edge_coords[~constant_idx, 1], "g.-")
 
         ax.set_xlabel("mm")
         ax.set_ylabel("mm")
@@ -323,7 +341,7 @@ def refine_edge_coord(img, img_res, edge_coords, normals, do_plots=True):
         plt.show()
         # fig.savefig(f"{OUTDIR}/refined_coords.png", bbox_inches="tight", dpi=150)
 
-    return refined_edge_coords
+    return refined_edge_coords, ~constant_idx
 
 
 def image_props(img, img_resolution, do_plots=False, plot_name=None):