SE FB registration

.pre-commit-config.yaml

@@ -1,6 +1,6 @@
 repos:
   - repo: https://github.com/timothycrosley/isort
-    rev: "6.0.1"
+    rev: "8.0.1"
     hooks:
       - id: isort
         args: ["--profile", "black"]

pyproject.toml

@@ -52,6 +52,7 @@ dependencies = [
   "PyWavelets",
   "pynrrd",
   "python-gdcm",
+  "antspyx",
 ]
 description = "INDI is a command line tool to process in-vivo cardiac diffusion tensor imaging."
 keywords = []

settings_template.yaml

@@ -23,8 +23,8 @@ registration_extra_debug: False
 # ========================================================================================
 # IMAGE REGISTRATION
 # ========================================================================================
-# none, quick_rigid, elastix_rigid, elastix_affine, elastix_non_rigid
-# elastix_groupwise (not working well at the moment)
+# none, quick_rigid, elastix_rigid, elastix_affine, elastix_non_rigid, elastix_non_rigid_fb
+# elastix_groupwise (not working well at the moment), ants_non_rigid
 registration: elastix_non_rigid
 # fast, slow
 registration_speed: fast

src/indi/extensions/crop_fov.py

@@ -88,10 +88,18 @@ def crop_images(
             segmentation[slice_name]["epicardium"] = np.array(segmentation[slice_name]["epicardium"]) - np.flip(
                 first_corner
             )
+        if segmentation[slice_name]["epicardium_true_border"].size != 0:
+            segmentation[slice_name]["epicardium_true_border"] = np.array(
+                segmentation[slice_name]["epicardium_true_border"]
+            ) - np.flip(first_corner)
         if segmentation[slice_name]["endocardium"].size != 0:
             segmentation[slice_name]["endocardium"] = np.array(segmentation[slice_name]["endocardium"]) - np.flip(
                 first_corner
             )
+        if segmentation[slice_name]["endocardium_true_border"].size != 0:
+            segmentation[slice_name]["endocardium_true_border"] = np.array(
+                segmentation[slice_name]["endocardium_true_border"]
+            ) - np.flip(first_corner)
 
     # add crop info to info dictionary
     temp_val = list(first_corner)
@@ -185,7 +193,7 @@ def record_image_registration(
         ]
         store_v_lp_post = np.mean(store_v_lp_post, axis=2)
 
-        plt.figure(figsize=(5, 5))
+        plt.figure(figsize=(10, 10))
         plt.subplot(2, 2, 1)
         plt.imshow(store_h_lp_pre, cmap="inferno")
         plt.axis("off")
@@ -246,6 +254,14 @@ def record_image_registration(
         # step of the vector field for the displacement transform
         step = 3
 
+        def fake_color_reg(a, b):
+            factor = 2.0
+            fake_colour_pre = np.zeros((b.shape[0], b.shape[1], 3), dtype=np.uint8)
+            fake_colour_pre[..., 0] = np.clip(a / np.max(a) * factor * 255, 0, 255)
+            fake_colour_pre[..., 2] = np.clip(b / np.max(b) * factor * 255, 0, 255)
+            fake_colour_pre[..., 1] = 0
+            return fake_colour_pre
+
         for slice_idx in slices:
             X, Y = np.meshgrid(
                 np.arange(0, registration_image_data[slice_idx]["deformation_field"]["grid"][1].shape[1], step),
@@ -271,8 +287,8 @@ def record_image_registration(
                     np.max(registration_image_data[slice_idx]["img_post_reg"][img_idx]),
                 )
                 comp_1 = compare_images(c_ref, c_img_post, method="checkerboard")
-                comp_2 = compare_images(c_ref, c_img_post, method="diff")
-                comp_3 = compare_images(c_ref, c_img_post, method="blend")
+                comp_2 = fake_color_reg(c_ref, c_img_post)
+                comp_3 = fake_color_reg(c_ref, c_img_pre)
 
                 plt.figure(figsize=(5, 5))
 
@@ -299,12 +315,12 @@ def record_image_registration(
                 plt.subplot(3, 3, 5)
                 plt.imshow(comp_2)
                 plt.axis("off")
-                plt.title("diff", fontsize=7)
+                plt.title("fake colour blend", fontsize=7)
 
                 plt.subplot(3, 3, 6)
                 plt.imshow(comp_3)
                 plt.axis("off")
-                plt.title("blend", fontsize=7)
+                plt.title("fake colour blend", fontsize=7)
 
                 plt.subplot(3, 3, 7)
                 plt.imshow(

src/indi/extensions/export_vectors_tensors_vtk.py

@@ -66,7 +66,7 @@ def save_vtk_file(
     # save surface mesh as png
     if save_mesh:
         vol = mesh.threshold(value=1, scalars="mask", invert=False)
-        surf = vol.extract_geometry()
+        surf = vol.extract_surface(algorithm=None)
         surf_smooth = surf.smooth_taubin(n_iter=50, pass_band=0.05)
 
         surf_smooth.save(os.path.join(folder_path, "surface_mesh" + ".ply"))
@@ -75,7 +75,7 @@ def save_vtk_file(
         # save primary eigenvector field
         # it will only save the vectors in the surface, so if multislice we won't see all the vectors
         vol = mesh.threshold(value=1, scalars="mask", invert=False)
-        surf = vol.extract_geometry()
+        surf = vol.extract_surface(algorithm=None)
         surf.set_active_vectors("primary_evs")
         surf.arrows.save(os.path.join(folder_path, "primary_eigenvectors" + ".ply"))
 

src/indi/extensions/extensions.py

@@ -1007,25 +1007,43 @@ def plot_results_maps(
         # plt.imshow(average_images[slice_idx], cmap="Blues_r", vmin=0, vmax=1)
         vmin, vmax = get_window(dti["s0"][slice_idx], mask_3c[slice_idx])
         plt.imshow(dti["s0"][slice_idx], cmap="Greys_r", vmin=vmin, vmax=vmax)
-        if segmentation[slice_idx]["epicardium"].size != 0:
-            plt.plot(
-                segmentation[slice_idx]["epicardium"][:, 0],
-                segmentation[slice_idx]["epicardium"][:, 1],
-                ".",
+        plt.colorbar(fraction=0.046, pad=0.04)
+        if segmentation[slice_idx]["epicardium_true_border"].size != 0:
+            plt.scatter(
+                segmentation[slice_idx]["epicardium_true_border"][:, 0],
+                segmentation[slice_idx]["epicardium_true_border"][:, 1],
+                marker=".",
+                s=10,
                 color="tab:blue",
-                markersize=0.5,
-                alpha=1.0,
+                alpha=0.8,
             )
-        if segmentation[slice_idx]["endocardium"].size != 0:
-            plt.plot(
-                segmentation[slice_idx]["endocardium"][:, 0],
-                segmentation[slice_idx]["endocardium"][:, 1],
-                ".",
+        if segmentation[slice_idx]["endocardium_true_border"].size != 0:
+            plt.scatter(
+                segmentation[slice_idx]["endocardium_true_border"][:, 0],
+                segmentation[slice_idx]["endocardium_true_border"][:, 1],
+                marker=".",
+                s=10,
                 color="tab:red",
-                markersize=0.5,
-                alpha=1.0,
+                alpha=0.8,
+            )
+        if segmentation[slice_idx]["anterior_ip"].size != 0:
+            plt.plot(
+                segmentation[slice_idx]["anterior_ip"][0],
+                segmentation[slice_idx]["anterior_ip"][1],
+                "2",
+                color="tab:orange",
+                markersize=20,
+                alpha=0.8,
+            )
+        if segmentation[slice_idx]["inferior_ip"].size != 0:
+            plt.plot(
+                segmentation[slice_idx]["inferior_ip"][0],
+                segmentation[slice_idx]["inferior_ip"][1],
+                "1",
+                color="tab:orange",
+                markersize=20,
+                alpha=0.8,
             )
-        plt.colorbar(fraction=0.046, pad=0.04)
         plt.tight_layout(pad=1.0)
         plt.axis("off")
         plt.title("S0")

src/indi/extensions/heart_segmentation.py

@@ -218,8 +218,22 @@ def heart_segmentation(
                 allow_pickle=True,
             )
             mask_3c[slice_idx] = npzfile["mask_3c"]
-            segmentation[slice_idx] = npzfile["segmentation"]
-            segmentation[slice_idx] = segmentation[slice_idx].item()
+            segmentation[slice_idx] = npzfile["segmentation"].item()
+
+            # Backwards-compatibility: ensure *_true_border keys exist for downstream plotting/cropping.
+            if segmentation[slice_idx].get("epicardium", np.array([])).size != 0 and (
+                "epicardium_true_border" not in segmentation[slice_idx]
+                or segmentation[slice_idx]["epicardium_true_border"].size == 0
+            ):
+                mask_lv = mask_3c[slice_idx].copy()
+                mask_lv[mask_lv != 1] = 0
+                epi_contour, endo_contour = get_sa_contours(mask_lv)
+                segmentation[slice_idx]["epicardium_true_border"] = epi_contour
+                if segmentation[slice_idx].get("endocardium", np.array([])).size != 0:
+                    segmentation[slice_idx]["endocardium_true_border"] = endo_contour
+
+            segmentation[slice_idx].setdefault("epicardium_true_border", np.empty((0, 2)))
+            segmentation[slice_idx].setdefault("endocardium_true_border", np.empty((0, 2)))
 
             # if there is no epicardial border defined, mark this slice to be removed in the dataframe
             if segmentation[slice_idx]["epicardium"].size == 0:
@@ -277,13 +291,16 @@ def heart_segmentation(
                     epi_contour = get_epi_contour(mask_lv)
                     endo_contour = np.array([])
 
+                segmentation[slice_idx]["epicardium_true_border"] = epi_contour
+                segmentation[slice_idx]["endocardium_true_border"] = endo_contour
+
                 epi_len = len(epi_contour)
                 endo_len = len(endo_contour)
-                epi_contour = spline_interpolate_contour(epi_contour, 20, join_ends=False)
+                epi_contour = spline_interpolate_contour(epi_contour, 10, join_ends=False)
                 epi_contour = spline_interpolate_contour(epi_contour, epi_len, join_ends=False)
 
                 if segmentation[slice_idx]["endocardium"].size != 0:
-                    endo_contour = spline_interpolate_contour(endo_contour, 20, join_ends=False)
+                    endo_contour = spline_interpolate_contour(endo_contour, 10, join_ends=False)
                     endo_contour = spline_interpolate_contour(endo_contour, endo_len, join_ends=False)
 
                 segmentation[slice_idx]["epicardium"] = epi_contour