Update spline_volume.py for contour: vertical alignment based on Eucl…

…idean distance

src/pyhexspline/spline_volume.py

@@ -1023,19 +1023,28 @@ def process_slice(self, mask, slice_idx):
     def get_line_sets(self, contour):
         num_points = contour[0].shape[0]
         line_sets = []
+
+        # Align slices based on the first point
+        first_point = contour[0][0]
+        for i in range(1, len(contour)):
+            distances = np.linalg.norm(contour[i] - first_point, axis=1)
+            min_idx = np.argmin(distances)
+            contour[i] = np.roll(contour[i], -min_idx, axis=0)
+
         for pt_idx in range(num_points):
             vertical_ll = []
             for _slice in contour:
                 # Collecting the pt_idx-th point from each slice
                 vertical_ll.append(_slice[pt_idx])
             line_sets.append(vertical_ll)
 
-        # plot every line_set with a different color
         if self.show_plots:
             for line_set in line_sets:
                 line_set = np.array(line_set)
                 plt.plot(line_set[:, 0], line_set[:, 1])
+                plt.title("Matched vertical lines")
             plt.show()
+
         return line_sets
 
     def plot_slices_with_slider(self, imnp, total_slices, results):
@@ -1093,13 +1102,9 @@ def volume_splines_optimized(self, imsitk_pad) -> Tuple[ndarray, ndarray]:
         height = imnp.shape[2]
         self.logger.debug(f"Height:\t{height}")
 
-        NUM_SLICES = height // 10  # TODO: parametrize this 10
+        NUM_SLICES = height // 10  # TODO: parametrize this value
         total_slices = np.linspace(0, height - 1, NUM_SLICES, dtype=int)
 
-        # results = Parallel(n_jobs=-1)(
-        #     delayed(self.process_slice)(imnp, slice_idx) for slice_idx in total_slices
-        # )
-
         processed_slices = []
         for slice_idx in total_slices:
             print(f"Processing slice {slice_idx}")

standalone.py

@@ -46,7 +46,7 @@ def main():
         "outside_val": int(1),  # threshold value for the outside of the mask
         "lower_thresh": float(0),  # lower threshold for the mask
         "upper_thresh": float(0.9),  # upper threshold for the mask
-        "s": 100,  # smoothing factor of the spline
+        "s": 50,  # smoothing factor of the spline
         "k": 3,  # degree of the spline
         "interp_points": 500,  # number of points to interpolate the spline
         "dp_simplification": 2,  # Ramer-Douglas-Peucker simplification factor
@@ -59,7 +59,7 @@ def main():
         "n_elms_transverse_cort": 3,  # number of elements in the transverse direction for the cortical compartment
         "n_elms_radial": 20,  # number of elements in the radial direction # ! Should be 10 if trab_refinement is True
         "ellipsoid_fitting": True,  # True: perform ellipsoid fitting
-        "show_plots": False,  # show plots during construction
+        "show_plots": True,  # show plots during construction
         "show_gmsh": True,  # show gmsh GUI
         "write_mesh": True,  # write mesh to file
         "trab_refinement": False,  # True: refine trabecular mesh at the center
@@ -74,7 +74,7 @@ def main():
     # print(sitk_image_s.GetSize())
 
     sitk_image_s = sitk.ReadImage(
-        "/home/simoneponcioni/Desktop/repro-to-remesh/C0002211_CORTMASK.mhd"
+        "/home/simoneponcioni/Desktop/repro-to-remesh/C0002215_CORTMASK.mhd"
     )
     print(sitk_image_s.GetSize())
     # sitk_image_s = sitk_image_s[:, :, :-40]