diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 3cf1ddb..94b3bbe 100644
--- a/.pre-commit-config.yaml
+++ b/.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"]
diff --git a/pyproject.toml b/pyproject.toml
index 93c9438..1e33656 100644
--- a/pyproject.toml
+++ b/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 = []
diff --git a/settings_template.yaml b/settings_template.yaml
index fff9417..bba193c 100644
--- a/settings_template.yaml
+++ b/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
diff --git a/src/indi/extensions/crop_fov.py b/src/indi/extensions/crop_fov.py
index ed78634..034326b 100644
--- a/src/indi/extensions/crop_fov.py
+++ b/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(
diff --git a/src/indi/extensions/export_vectors_tensors_vtk.py b/src/indi/extensions/export_vectors_tensors_vtk.py
index 6ff988f..94019dc 100644
--- a/src/indi/extensions/export_vectors_tensors_vtk.py
+++ b/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"))
 
diff --git a/src/indi/extensions/extensions.py b/src/indi/extensions/extensions.py
index 3199931..f7ed10b 100644
--- a/src/indi/extensions/extensions.py
+++ b/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")
diff --git a/src/indi/extensions/heart_segmentation.py b/src/indi/extensions/heart_segmentation.py
index 318b51b..041a174 100644
--- a/src/indi/extensions/heart_segmentation.py
+++ b/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
diff --git a/src/indi/extensions/image_registration.py b/src/indi/extensions/image_registration.py
index beee8cf..34a561c 100755
--- a/src/indi/extensions/image_registration.py
+++ b/src/indi/extensions/image_registration.py
@@ -57,7 +57,7 @@ def denoise_img_nlm(c_img: NDArray) -> NDArray:
     # estimate the noise standard deviation from the noisy image
     sigma_est = np.mean(estimate_sigma(c_img, channel_axis=None))
     # fast algorithm, sigma provided
-    denoised_img = denoise_nl_means(c_img, h=10 * sigma_est, sigma=sigma_est, fast_mode=True, **patch_kw)
+    denoised_img = denoise_nl_means(c_img, h=5 * sigma_est, sigma=sigma_est, fast_mode=True, **patch_kw)
 
     return denoised_img
 
@@ -147,6 +147,26 @@ def registration_loop(
         #             "current_parameters{0}.txt".format(index),
         #         ),
         #     )
+    if settings["registration"] == "elastix_non_rigid_fb":
+        parameter_object = itk.ParameterObject.New()
+        parameter_object.AddParameterFile(os.path.join(script_path, "image_registration_recipes", "Elastix_rigid.txt"))
+        if settings["registration_speed"] == "slow":
+            parameter_object.SetParameter("MaximumNumberOfIterations", "2000")
+            parameter_object.SetParameter("NumberOfResolutions", "4")
+
+        parameter_object.AddParameterFile(
+            os.path.join(script_path, "image_registration_recipes", "Elastix_affine.txt")
+        )
+        if settings["registration_speed"] == "slow":
+            parameter_object.SetParameter("MaximumNumberOfIterations", "2000")
+            parameter_object.SetParameter("NumberOfResolutions", "4")
+
+        parameter_object.AddParameterFile(
+            os.path.join(script_path, "image_registration_recipes", "Elastix_bspline_fb.txt")
+        )
+        if settings["registration_speed"] == "slow":
+            parameter_object.SetParameter("MaximumNumberOfIterations", "2000")
+            parameter_object.SetParameter("NumberOfResolutions", "4")
 
     if settings["registration"] == "elastix_groupwise":
         parameter_object = itk.ParameterObject.New()
@@ -175,6 +195,9 @@ def registration_loop(
         #         ),
         #     )
 
+    if settings["registration"] == "ants_non_rigid":
+        import ants
+
     # dicts to store information about the registration
     registration_image_data = {}
 
@@ -202,16 +225,27 @@ def registration_loop(
         settings["registration"] == "elastix_rigid"
         or settings["registration"] == "elastix_affine"
         or settings["registration"] == "elastix_non_rigid"
+        or settings["registration"] == "elastix_non_rigid_fb"
     ):
         ref = itk.GetImageFromArray(ref)
+    if settings["registration"] == "ants_non_rigid":
+        ref = ants.from_numpy(ref)
 
     # images to be registered
     mov_all = np.transpose(np.dstack(data["image"].values), (2, 0, 1))
     if settings["complex_data"]:
         mov_all_phase = np.transpose(np.dstack(data["image_phase"].values), (2, 0, 1))
 
-    # mask only regions of interest in the middle of the FOV.
-    mask = itk.GetImageFromArray(mask)
+    if (
+        settings["registration"] == "elastix_rigid"
+        or settings["registration"] == "elastix_affine"
+        or settings["registration"] == "elastix_non_rigid"
+        or settings["registration"] == "elastix_non_rigid_fb"
+    ):
+        # mask only regions of interest in the middle of the FOV.
+        mask = itk.GetImageFromArray(mask)
+
+        # masking makes things worse for ANTs registration
 
     # Groupwise registration
     if settings["registration"] == "elastix_groupwise":
@@ -262,10 +296,16 @@ def registration_loop(
                 settings["registration"] == "elastix_rigid"
                 or settings["registration"] == "elastix_affine"
                 or settings["registration"] == "elastix_non_rigid"
+                or settings["registration"] == "elastix_non_rigid_fb"
             ):
                 # apply the registration to a denoised version (helps with registration of low SNR images)
-                mov_norm = (mov - np.min(mov)) / (np.max(mov) - np.min(mov))
-                denoised_mov = denoise_img_nlm(mov_norm)
+                # denoisng used only for STEAM images, for SE the SNR is usually good enough without denoising
+                if settings["sequence_type"] == "steam":
+                    mov_norm = (mov - np.min(mov)) / (np.max(mov) - np.min(mov))
+                    denoised_mov = denoise_img_nlm(mov_norm)
+                else:
+                    mov_norm = (mov - np.min(mov)) / (np.max(mov) - np.min(mov))
+                    denoised_mov = mov_norm
 
                 denoised_mov = itk.GetImageFromArray(denoised_mov)
 
@@ -310,6 +350,63 @@ def registration_loop(
                     img_reg = itk.transformix_filter(mov, result_transform_parameters)
                     img_reg = itk.GetArrayFromImage(img_reg)
 
+            # ANTs registration
+            elif settings["registration"] == "ants_non_rigid":
+                if settings["complex_data"]:
+                    # complex data registration
+                    c_real = np.multiply(mov, np.cos(mov_phase))
+                    c_real = ants.from_numpy(c_real)
+                    c_imag = np.multiply(mov, np.sin(mov_phase))
+                    c_imag = ants.from_numpy(c_imag)
+                    mov = ants.from_numpy(mov)
+                    mytx = ants.registration(
+                        fixed=ref,
+                        moving=mov,
+                        type_of_transform="antsRegistrationSyN[b]",
+                        mask=ants.from_numpy(mask),
+                        moving_mask=ants.from_numpy(mask),
+                    )
+                    img_reg_real = ants.apply_transforms(
+                        fixed=ref,
+                        moving=c_real,
+                        transformlist=mytx["fwdtransforms"],
+                        whichtoinvert=[False, False],
+                    )
+                    img_reg_real = img_reg_real.numpy()
+                    img_reg_imag = ants.apply_transforms(
+                        fixed=ref,
+                        moving=c_imag,
+                        transformlist=mytx["fwdtransforms"],
+                        whichtoinvert=[False, False],
+                    )
+                    img_reg_imag = img_reg_imag.numpy()
+                    img_reg = np.sqrt(np.square(img_reg_real) + np.square(img_reg_imag))
+                    img_phase_reg = np.arctan2(img_reg_imag, img_reg_real)
+
+                else:
+                    # get moving image and register it to the reference
+                    mov = ants.from_numpy(mov)
+                    mytx = ants.registration(
+                        fixed=ref,
+                        moving=mov,
+                        type_of_transform="antsRegistrationSyN[b]",
+                        mask=ants.from_numpy(mask),
+                        moving_mask=ants.from_numpy(mask),
+                    )
+                    img_reg = mytx["warpedmovout"].numpy()
+
+                # apply the deformation field to the grid image
+                grid_img = get_grid_image(info["img_size"], 6)
+                grid_img = ants.from_numpy(grid_img)
+                grid_img = ants.apply_transforms(
+                    fixed=ref,
+                    moving=grid_img,
+                    transformlist=mytx["fwdtransforms"],
+                    whichtoinvert=[False, False],
+                )
+                grid_img = grid_img.numpy()
+                registration_image_data["deformation_field"]["grid"][i] = grid_img
+
             # basic quick rigid
             elif settings["registration"] == "quick_rigid":
                 shift, error, diffphase = phase_cross_correlation(
@@ -439,8 +536,11 @@ def get_ref_image(current_entries: pd.DataFrame, slice_idx: int, settings: dict,
             # normalise 0 to 1
             c_img = (c_img - np.min(c_img)) / (np.max(c_img) - np.min(c_img))
 
-            # denoise image
-            denoised_img = denoise_img_nlm(c_img)
+            # denoise image (only for STEAM, for SE the SNR is usually good enough without denoising)
+            if settings["sequence_type"] == "steam":
+                denoised_img = denoise_img_nlm(c_img)
+            else:
+                denoised_img = c_img
 
             ref_images["image"] = denoised_img
             ref_images["index"] = index_pos[np.argmax(image_stack_sum)]
@@ -460,27 +560,27 @@ def get_ref_image(current_entries: pd.DataFrame, slice_idx: int, settings: dict,
             image_stack = np.stack(current_entries["image"][index_pos].values)
 
             # groupwise registration recipe
-            parameter_object = itk.ParameterObject.New()
-            parameter_object.AddParameterFile(
-                os.path.join(
-                    os.path.dirname(os.path.dirname(__file__)),
-                    "image_registration_recipes",
-                    "Elastix_groupwise.txt",
-                )
-            )
-
             # parameter_object = itk.ParameterObject.New()
-            # groupwise_parameter_map = parameter_object.GetDefaultParameterMap("groupwise")
-            # parameter_object.AddParameterMap(groupwise_parameter_map)
+            # parameter_object.AddParameterFile(
+            #     os.path.join(
+            #         os.path.dirname(os.path.dirname(__file__)),
+            #         "image_registration_recipes",
+            #         "Elastix_groupwise.txt",
+            #     )
+            # )
+            parameter_object = itk.ParameterObject.New()
+            groupwise_parameter_map = parameter_object.GetDefaultParameterMap("groupwise")
+            parameter_object.AddParameterMap(groupwise_parameter_map)
 
             # modify type for itk
             image_stack = np.ascontiguousarray(np.array(image_stack, dtype=np.float32))
             # store images before registration
             img_pre = image_stack
 
-            # denoise stack before masking
-            for i in range(image_stack.shape[0]):
-                image_stack[i] = denoise_img_nlm(image_stack[i])
+            # denoise stack before masking (for STEAM images, for SE the SNR is usually good enough without denoising)
+            if settings["sequence_type"] == "steam":
+                for i in range(image_stack.shape[0]):
+                    image_stack[i] = denoise_img_nlm(image_stack[i])
 
             # create mask stack of the FOV central region
             mask = np.zeros([image_stack.shape[1], image_stack.shape[2]])
@@ -606,7 +706,7 @@ def plot_ref_images(
                 img_reg = ref_images[slice_idx]["groupwise_reg_info"]["post"]
                 c_ref = ref_images[slice_idx]["image"]
 
-                plt.figure()
+                plt.figure(figsize=(n_images, 5))
                 for i in range(n_images):
                     plt.subplot(4, n_images, i + 1)
                     plt.imshow(img_pre[i], vmin=np.min(img_pre[i]), vmax=np.max(img_pre[i]) * 0.3, cmap="Greys_r")
diff --git a/src/indi/extensions/image_registration_recipes/Elastix_bspline_fb.txt b/src/indi/extensions/image_registration_recipes/Elastix_bspline_fb.txt
new file mode 100644
index 0000000..2bbda16
--- /dev/null
+++ b/src/indi/extensions/image_registration_recipes/Elastix_bspline_fb.txt
@@ -0,0 +1,182 @@
+// Example parameter file for rotation registration
+// C-style comments: //
+
+// The internal pixel type, used for internal computations
+// Leave to float in general.
+// NB: this is not the type of the input images! The pixel
+// type of the input images is automatically read from the
+// images themselves.
+// This setting can be changed to "short" to save some memory
+// in case of very large 3D images.
+(FixedInternalImagePixelType "float")
+(MovingInternalImagePixelType "float")
+
+// The dimensions of the fixed and moving image
+// NB: This has to be specified by the user. The dimension of
+// the images is currently NOT read from the images.
+// Also note that some other settings may have to specified
+// for each dimension separately.
+(FixedImageDimension 2)
+(MovingImageDimension 2)
+
+// Specify whether you want to take into account the so-called
+// direction cosines of the images. Recommended: true.
+// In some cases, the direction cosines of the image are corrupt,
+// due to image format conversions for example. In that case, you
+// may want to set this option to "false".
+(UseDirectionCosines "true")
+
+// **************** Main Components **************************
+
+// The following components should usually be left as they are:
+// (Registration "MultiResolutionRegistration")
+(Registration "MultiMetricMultiResolutionRegistration")
+(Interpolator "BSplineInterpolator")
+(ResampleInterpolator "FinalBSplineInterpolator")
+(Resampler "DefaultResampler")
+
+// These may be changed to Fixed/MovingSmoothingImagePyramid.
+// See the manual.
+(FixedImagePyramid "FixedRecursiveImagePyramid")
+(MovingImagePyramid "MovingRecursiveImagePyramid")
+//(FixedImagePyramid "FixedShrinkingImagePyramid")
+//(MovingImagePyramid "MovingShrinkingImagePyramid")
+//(WriteResultImageAfterEachResolution "true" "true" "true" "true")
+
+// The following components are most important:
+// The optimizer AdaptiveStochasticGradientDescent (ASGD) works
+// quite ok in general. The Transform and Metric are important
+// and need to be chosen careful for each application. See manual.
+(Optimizer "AdaptiveStochasticGradientDescent")
+(Transform "BSplineTransform")
+//(Metric "AdvancedMattesMutualInformation")
+//(Metric "AdvancedNormalizedCorrelation")
+
+//pf_rbh
+//(Metric "AdvancedNormalizedCorrelation" "TransformBendingEnergyPenalty")
+(Metric "AdvancedMattesMutualInformation" "TransformBendingEnergyPenalty")
+
+// ***************** Transformation **************************
+
+// Scales the rotations compared to the translations, to make
+// sure they are in the same range. In general, it's best to
+// use automatic scales estimation:
+//(AutomaticScalesEstimation "true")
+
+// Automatically guess an initial translation by aligning the
+// geometric centers of the fixed and moving.
+//(AutomaticTransformInitialization "true")
+
+//pf_rbh
+(FinalGridSpacingInVoxels  16)
+
+// Whether transforms are combined by composition or by addition.
+// In generally, Compose is the best option in most cases.
+// It does not influence the results very much.
+(HowToCombineTransforms "Compose")
+
+// ******************* Similarity measure *********************
+
+// Number of grey level bins in each resolution level,
+// for the mutual information. 16 or 32 usually works fine.
+// You could also employ a hierarchical strategy:
+//(NumberOfHistogramBins 16 32 64)
+(NumberOfHistogramBins 32)
+
+// If you use a mask, this option is important.
+// If the mask serves as region of interest, set it to false.
+// If the mask indicates which pixels are valid, then set it to true.
+// If you do not use a mask, the option doesn't matter.
+//pf_rbh
+(ErodeMask "true")
+
+//pf_rbh
+(Metric0Weight 1.0)
+(Metric1Weight 10.0)
+// or define (experimental):
+//(Metric0RelativeWeight 1.0)
+//(Metric1RelativeWeight 0.1)
+
+
+
+// ******************** Multiresolution **********************
+
+// The number of resolutions. 1 Is only enough if the expected
+// deformations are small. 3 or 4 mostly works fine. For large
+// images and large deformations, 5 or 6 may even be useful.
+//pf_rbh
+(NumberOfResolutions 2)
+
+// The downsampling/blurring factors for the image pyramids.
+// By default, the images are downsampled by a factor of 2
+// compared to the next resolution.
+// So, in 2D, with 4 resolutions, the following schedule is used:
+//(ImagePyramidSchedule 8 8  4 4  2 2  1 1 )
+// And in 3D:
+//(ImagePyramidSchedule 1 1 1  1 1 1)
+// You can specify any schedule, for example:
+//(ImagePyramidSchedule 4 4  4 3  2 1  1 1 )
+// Make sure that the number of elements equals the number
+// of resolutions times the image dimension.
+
+// ******************* Optimizer ****************************
+
+// Maximum number of iterations in each resolution level:
+// 200-500 works usually fine for rigid registration.
+// For more robustness, you may increase this to 1000-2000.
+(MaximumNumberOfIterations 256)
+
+// The step size of the optimizer, in mm. By default the voxel size is used.
+// which usually works well. In case of unusual high-resolution images
+// (eg histology) it is necessary to increase this value a bit, to the size
+// of the "smallest visible structure" in the image:
+(MaximumStepLength 0.117188)
+
+// **************** Image sampling **********************
+
+// Number of spatial samples used to compute the mutual
+// information (and its derivative) in each iteration.
+// With an AdaptiveStochasticGradientDescent optimizer,
+// in combination with the two options below, around 2000
+// samples may already suffice.
+//(RandomSeed 30101983) Doesn't seem to work
+(NumberOfSpatialSamples 2048)
+
+// Refresh these spatial samples in every iteration, and select
+// them randomly. See the manual for information on other sampling
+// strategies.
+(NewSamplesEveryIteration "true")
+
+//(ImageSampler "Random")
+(ImageSampler "Grid")
+(SampleGridSpacing 2 2 1 1)
+//(ImageSampler "Full")
+
+// ************* Interpolation and Resampling ****************
+
+// Order of B-Spline interpolation used during registration/optimisation.
+// It may improve accuracy if you set this to 3. Never use 0.
+// An order of 1 gives linear interpolation. This is in most
+// applications a good choice.
+(BSplineInterpolationOrder 1)
+
+// Order of B-Spline interpolation used for applying the final
+// deformation.
+// 3 gives good accuracy; recommended in most cases.
+// 1 gives worse accuracy (linear interpolation)
+// 0 gives worst accuracy, but is appropriate for binary images
+// (masks, segmentations); equivalent to nearest neighbor interpolation.
+(FinalBSplineInterpolationOrder 3)
+
+//Default pixel value for pixels that come from outside the picture:
+(DefaultPixelValue 0)
+
+// Choose whether to generate the deformed moving image.
+// You can save some time by setting this to false, if you are
+// only interested in the final (nonrigidly) deformed moving image
+// for example.
+(WriteResultImage "true")
+
+// The pixel type and format of the resulting deformed moving image
+(ResultImagePixelType "float")
+(ResultImageFormat "mhd")
\ No newline at end of file
diff --git a/src/indi/extensions/manual_lv_segmentation.py b/src/indi/extensions/manual_lv_segmentation.py
index 21c43be..6c64e1e 100644
--- a/src/indi/extensions/manual_lv_segmentation.py
+++ b/src/indi/extensions/manual_lv_segmentation.py
@@ -216,6 +216,65 @@ def plot_manual_lv_segmentation(
         save_path (str): Directory where the output PNG files are saved.
     """
 
+    for slice_idx in slices:
+        # alpha mask
+        alphas_myocardium = np.copy(mask_3c[slice_idx])
+        alphas_myocardium[alphas_myocardium == 2] = 0
+        alphas_myocardium[alphas_myocardium > 0.1] = 0.3
+
+        # plot average images and respective masks
+        plt.figure(figsize=(5, 5))
+        plt.imshow(average_maps[slice_idx], cmap="Greys_r")
+        plt.imshow(alphas_myocardium, alpha=alphas_myocardium, vmin=0, vmax=1, cmap="hot")
+        plt.axis("off")
+        if segmentation[slice_idx]["epicardium"].size != 0:
+            plt.scatter(
+                segmentation[slice_idx]["epicardium_true_border"][:, 0],
+                segmentation[slice_idx]["epicardium_true_border"][:, 1],
+                marker=".",
+                s=2,
+                color="tab:blue",
+                alpha=0.5,
+            )
+        if segmentation[slice_idx]["endocardium"].size != 0:
+            plt.scatter(
+                segmentation[slice_idx]["endocardium_true_border"][:, 0],
+                segmentation[slice_idx]["endocardium_true_border"][:, 1],
+                marker=".",
+                s=2,
+                color="tab:red",
+                alpha=0.5,
+            )
+        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=10,
+                alpha=0.5,
+            )
+        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=10,
+                alpha=0.5,
+            )
+        plt.savefig(
+            os.path.join(
+                save_path,
+                filename + "_true_border_slice_" + str(slice_idx).zfill(2) + ".png",
+            ),
+            dpi=100,
+            bbox_inches="tight",
+            pad_inches=0,
+            transparent=False,
+        )
+        plt.close()
+
     for slice_idx in slices:
         # alpha mask
         alphas_myocardium = np.copy(mask_3c[slice_idx])
@@ -260,7 +319,7 @@ def plot_manual_lv_segmentation(
                 segmentation[slice_idx]["inferior_ip"][1],
                 "1",
                 color="tab:orange",
-                markersize=5,
+                markersize=10,
                 alpha=0.5,
             )
         plt.savefig(
@@ -630,7 +689,7 @@ def manual_lv_segmentation(
             2,
             dpi=my_dpi,
             figsize=(settings["screen_size"][0] / my_dpi, (settings["screen_size"][1] - 52) / my_dpi),
-            num="Slice " + str(slice_idx) + " of " + str(len(slices) - 1),
+            num="Slice " + str(slice_idx + 1) + " of " + str(len(slices)),
         )
     else:
         fig, ax = plt.subplots(
@@ -638,7 +697,7 @@ def manual_lv_segmentation(
             1,
             dpi=my_dpi,
             figsize=(settings["screen_size"][0] / my_dpi, (settings["screen_size"][1] - 52) / my_dpi),
-            num="Slice " + str(slice_idx) + " of " + str(len(slices) - 1),
+            num="Slice " + str(slice_idx + 1) + " of " + str(len(slices)),
         )
 
     # leave some space for the buttons
diff --git a/src/indi/extensions/read_data/read_and_pre_process_data.py b/src/indi/extensions/read_data/read_and_pre_process_data.py
index 4dd7277..fc398b7 100644
--- a/src/indi/extensions/read_data/read_and_pre_process_data.py
+++ b/src/indi/extensions/read_data/read_and_pre_process_data.py
@@ -683,7 +683,7 @@ def adjust_b_val_and_dir(
         and (data_type == "dicom" or data_type == "pandas")
         and info["manufacturer"] == "siemens"
     ):
-        if info["image_comments"]:
+        if "image_comments" in info.keys():
             logger.debug("Dicom header comment found: " + info["image_comments"])
             # get all numbers from comment field
             m = re.findall(r"[-+]?(?:\d*\.*\d+)", info["image_comments"])
@@ -981,10 +981,10 @@ def create_2d_montage_from_database(
         if segmentation:
             # create montage with segmentation
             seg_img = np.zeros((info["img_size"][0], info["img_size"][1], 3))
-            pts = np.array(segmentation[slice_int]["epicardium"], dtype=int)
+            pts = np.array(segmentation[slice_int]["epicardium_true_border"], dtype=int)
             seg_img[pts[:, 1], pts[:, 0]] = [1.0, 1.0, 0.33]
             if segmentation[slice_int]["endocardium"].size != 0:
-                pts = np.array(segmentation[slice_int]["endocardium"], dtype=int)
+                pts = np.array(segmentation[slice_int]["endocardium_true_border"], dtype=int)
                 seg_img[pts[:, 1], pts[:, 0]] = [1.0, 1.0, 0.33]
             # repeat image for the entire stack
             seg_img = np.tile(seg_img, (len(c_img_stack), max_number_of_images, 1))
diff --git a/src/indi/extensions/select_outliers.py b/src/indi/extensions/select_outliers.py
index 40b9ff1..9dc8b04 100644
--- a/src/indi/extensions/select_outliers.py
+++ b/src/indi/extensions/select_outliers.py
@@ -157,16 +157,16 @@ def manual_image_removal(
                     line_colour = c_colour
 
                     axs[idx, idx2].plot(
-                        segmentation[slice_idx]["epicardium"][:, 0],
-                        segmentation[slice_idx]["epicardium"][:, 1],
+                        segmentation[slice_idx]["epicardium_true_border"][:, 0],
+                        segmentation[slice_idx]["epicardium_true_border"][:, 1],
                         lw=0.5,
                         color=line_colour,
                         alpha=1.0,
                     )
                     if segmentation[slice_idx]["endocardium"].size != 0:
                         axs[idx, idx2].plot(
-                            segmentation[slice_idx]["endocardium"][:, 0],
-                            segmentation[slice_idx]["endocardium"][:, 1],
+                            segmentation[slice_idx]["endocardium_true_border"][:, 0],
+                            segmentation[slice_idx]["endocardium_true_border"][:, 1],
                             lw=0.5,
                             color=line_colour,
                             alpha=1.0,
diff --git a/src/indi/scripts/main.py b/src/indi/scripts/main.py
index 0e0117d..c558528 100644
--- a/src/indi/scripts/main.py
+++ b/src/indi/scripts/main.py
@@ -76,12 +76,6 @@ def main() -> None:
     # tf_keras package.
     os.environ["TF_USE_LEGACY_KERAS"] = "1"
 
-    # ITK
-    # import itk
-
-    # limit the amount of parallel threads during registration
-    # itk.MultiThreaderBase.SetGlobalMaximumNumberOfThreads(1)
-
     # DTCMR tailored colormaps
     colormaps = get_colourmaps()