Avoid uneven contour points error

CorpusCallosum/shape/thickness.py

@@ -113,7 +113,7 @@ def insert_point_with_thickness(
     point: np.ndarray,
     thickness_value: float,
     return_index: Literal[True],
-) -> tuple[np.ndarray, np.ndarray, int, bool]:
+) -> tuple[np.ndarray, np.ndarray, int]:
     ...
 
 
@@ -123,7 +123,7 @@ def insert_point_with_thickness(
     point: np.ndarray,
     thickness_value: float,
     return_index: bool = False
-) -> tuple[np.ndarray, np.ndarray, int, bool] | tuple[np.ndarray, np.ndarray]:
+) -> tuple[np.ndarray, np.ndarray, int] | tuple[np.ndarray, np.ndarray]:
     """Inserts a point and its thickness value into the contour.
 
     Parameters
@@ -142,23 +142,52 @@ def insert_point_with_thickness(
     Returns
     -------
     contour_in_as_space : np.ndarray
-        Updated contour. Shape is unchanged if the point already exists, otherwise (N+1, 2).
+        Updated contour of shape (N+1, 2).
     contour_thickness : np.ndarray
-        Updated thickness values. Shape is unchanged if the point already exists, otherwise (N+1,).
+        Updated thickness values of shape (N+1,).
     insertion_index : int
         The index, where the point was inserted (only if return_index is True).
-    inserted : bool
-        Whether a new contour point was inserted (only if return_index is True).
     """
     existing_distances = np.linalg.norm(contour_in_as_space[:, :2] - point[:2], axis=1)
     existing_matches = np.where(existing_distances <= 1e-10)[0]
     if len(existing_matches) > 0:
-        point_idx = existing_matches[0]
+        point_idx = int(existing_matches[0])
         if np.isnan(contour_thickness[point_idx]):
             contour_thickness[point_idx] = thickness_value
 
+        # Keep the contour cardinality stable across slices, but avoid inserting
+        # an exact duplicate vertex. Exact duplicates create zero-length edges in
+        # visualization/mesh code, while skipping the point breaks the equal
+        # point-count invariant expected by CCMesh.from_contours().
+        before_idx = (point_idx - 1) % len(contour_in_as_space)
+        after_idx = (point_idx + 1) % len(contour_in_as_space)
+        before_vector = contour_in_as_space[before_idx, :2] - contour_in_as_space[point_idx, :2]
+        after_vector = contour_in_as_space[after_idx, :2] - contour_in_as_space[point_idx, :2]
+        before_length = np.linalg.norm(before_vector)
+        after_length = np.linalg.norm(after_vector)
+
+        if after_length >= before_length and after_length > 0:
+            insert_idx = point_idx + 1
+            direction = after_vector / after_length
+            edge_length = after_length
+        elif before_length > 0:
+            insert_idx = point_idx
+            direction = before_vector / before_length
+            edge_length = before_length
+        else:
+            insert_idx = point_idx + 1
+            direction = np.zeros_like(point[:2])
+            edge_length = 0.0
+
+        if edge_length > 0:
+            step = min(max(edge_length * 1e-6, 1e-8), edge_length * 0.5)
+            point = contour_in_as_space[point_idx, :2] + direction * step
+
+        contour_in_as_space = np.insert(contour_in_as_space, insert_idx, point, axis=0)
+        contour_thickness = np.insert(contour_thickness, insert_idx, thickness_value)
+
         if return_index:
-            return contour_in_as_space, contour_thickness, point_idx, False
+            return contour_in_as_space, contour_thickness, insert_idx
         else:
             return contour_in_as_space, contour_thickness
 
@@ -170,7 +199,7 @@ def insert_point_with_thickness(
     contour_thickness = np.insert(contour_thickness, edge_idx + 1, thickness_value)
 
     if return_index:
-        return contour_in_as_space, contour_thickness, edge_idx + 1, True
+        return contour_in_as_space, contour_thickness, edge_idx + 1
     else:
         return contour_in_as_space, contour_thickness
 
@@ -335,25 +364,23 @@ def cc_thickness(
         levelpath_start = levelpath_asz[0, :2]
         levelpath_end = levelpath_asz[-1, :2]
 
-        contour_2d, contour_thickness, inserted_idx_start, inserted_start = insert_point_with_thickness(
+        contour_2d, contour_thickness, inserted_idx_start = insert_point_with_thickness(
             contour_2d, contour_thickness, levelpath_start, lvlpath_length, return_index=True,
         )
         # keep track of start index
-        if inserted_start:
-            if inserted_idx_start <= anterior_endpoint_idx:
-                anterior_endpoint_idx += 1
-            if inserted_idx_start <= posterior_endpoint_idx:
-                posterior_endpoint_idx += 1
+        if inserted_idx_start <= anterior_endpoint_idx:
+            anterior_endpoint_idx += 1
+        if inserted_idx_start <= posterior_endpoint_idx:
+            posterior_endpoint_idx += 1
 
-        contour_2d, contour_thickness, inserted_idx_end, inserted_end = insert_point_with_thickness(
+        contour_2d, contour_thickness, inserted_idx_end = insert_point_with_thickness(
             contour_2d, contour_thickness, levelpath_end, lvlpath_length, return_index=True,
         )
         # keep track of end index
-        if inserted_end:
-            if inserted_idx_end <= anterior_endpoint_idx:
-                anterior_endpoint_idx += 1
-            if inserted_idx_end <= posterior_endpoint_idx:
-                posterior_endpoint_idx += 1
+        if inserted_idx_end <= anterior_endpoint_idx:
+            anterior_endpoint_idx += 1
+        if inserted_idx_end <= posterior_endpoint_idx:
+            posterior_endpoint_idx += 1
 
     contour_2d_with_thickness = np.concatenate([contour_2d, contour_thickness[:, None]], axis=1)