refactor: use compas_cgal.isolines for contour extraction

Author: jf---

docs/concepts/slicing-algorithms.md

@@ -185,29 +185,21 @@ To build connected contours:
 
 ## Contour Assembly
 
-All slicers eventually produce contours via the `ScalarFieldContours` class:
-
-### From Crossings to Paths
+All slicers produce contours via `ScalarFieldContours`, which uses CGAL's isoline extraction:
 
 ```mermaid
 flowchart LR
-    A[Edge crossings] --> B[Face traversal]
-    B --> C[Connected polylines]
+    A[Scalar field on vertices] --> B[CGAL isolines]
+    B --> C[Sorted polylines]
     C --> D[Path objects]
 ```
 
-1. **Build crossing map**: Dictionary of edge → crossing point
-2. **Traverse faces**: Walk around faces connecting crossings
-3. **Handle branches**: Multiple paths per layer for complex geometry
-4. **Create Paths**: Wrap polylines in Path objects with metadata
-
-### Handling Complex Topology
-
-The algorithm handles:
+The CGAL backend (`compas_cgal.isolines`) handles:
 
-- **Multiple contours per layer**: Holes, disconnected regions
-- **Open contours**: When path hits mesh boundary
-- **Branching**: When contours merge or split
+- **Edge crossing detection**: Finding zero-crossings on mesh edges
+- **Polyline assembly**: Connecting crossings into coherent curves
+- **Multiple contours**: Holes, disconnected regions, branching
+- **Open/closed detection**: Identifying boundary-hitting paths
 
 ## Performance Considerations
 

docs/examples/05_scalar_field.md

@@ -193,9 +193,9 @@ Creates concentric circular layers (spiral vase mode).
 ### Geodesic Field
 
 ```python
-# Using igl for geodesic distance from boundary vertices
-import igl
-distances = igl.exact_geodesic(V, F, boundary_vertices)
+# Using CGAL for geodesic distance from boundary vertices
+from compas_cgal.geodesics import heat_geodesic_distances
+distances = heat_geodesic_distances((V, F), boundary_vertices)
 ```
 
 Creates layers that follow surface curvature.

examples/2_curved_slicing/ex2_curved_slicing.py

@@ -1,11 +1,12 @@
 import time
 from pathlib import Path
 
+import numpy as np
 from compas.datastructures import Mesh
 
 import compas_slicer.utilities as utils
 from compas_slicer.config import InterpolationConfig
-from compas_slicer.post_processing import seams_smooth, simplify_paths_rdp
+from compas_slicer.post_processing import seams_smooth
 from compas_slicer.pre_processing import InterpolationSlicingPreprocessor, create_mesh_boundary_attributes
 from compas_slicer.print_organization import (
     InterpolationPrintOrganizer,
@@ -28,9 +29,12 @@ def main(visualize: bool = False):
     # Load initial_mesh
     mesh = Mesh.from_obj(DATA_PATH / 'mesh.obj')
 
-    # Load targets (boundaries)
-    low_boundary_vs = utils.load_from_json(DATA_PATH, 'boundaryLOW.json')
-    high_boundary_vs = utils.load_from_json(DATA_PATH, 'boundaryHIGH.json')
+    # Identify boundaries from mesh topology
+    boundaries = [list(loop) for loop in mesh.vertices_on_boundaries()]
+    avg_zs = [np.mean([mesh.vertex_coordinates(v)[2] for v in loop]) for loop in boundaries]
+    low_idx = int(np.argmin(avg_zs))
+    low_boundary_vs = boundaries.pop(low_idx)
+    high_boundary_vs = [v for loop in boundaries for v in loop]  # flatten remaining
     create_mesh_boundary_attributes(mesh, low_boundary_vs, high_boundary_vs)
 
     avg_layer_height = 2.0
@@ -48,7 +52,6 @@ def main(visualize: bool = False):
     slicer = InterpolationSlicer(mesh, preprocessor, config)
     slicer.slice_model()
 
-    simplify_paths_rdp(slicer, threshold=0.25)
     seams_smooth(slicer, smooth_distance=3)
     slicer.printout_info()
     utils.save_to_json(slicer.to_data(), OUTPUT_PATH, 'curved_slicer.json')

src/compas_slicer/config.py

@@ -25,7 +25,6 @@
     "GcodeConfig",
     "PrintConfig",
     "OutputConfig",
-    "GeodesicsMethod",
     "UnionMethod",
     "load_defaults",
 ]
@@ -50,15 +49,6 @@ def load_defaults() -> dict[str, Any]:
 _DEFAULTS = load_defaults()
 
 
-class GeodesicsMethod(str, Enum):
-    """Method for computing geodesic distances."""
-
-    EXACT_IGL = "exact_igl"
-    HEAT_IGL = "heat_igl"
-    HEAT_CGAL = "heat_cgal"
-    HEAT = "heat"
-
-
 class UnionMethod(str, Enum):
     """Method for combining target boundaries."""
 
@@ -162,10 +152,6 @@ class InterpolationConfig(Data):
         Maximum layer height.
     vertical_layers_max_centroid_dist : float
         Maximum distance for grouping paths into vertical layers.
-    target_low_geodesics_method : GeodesicsMethod
-        Method for computing geodesics to low boundary.
-    target_high_geodesics_method : GeodesicsMethod
-        Method for computing geodesics to high boundary.
     target_high_union_method : UnionMethod
         Method for combining high target boundaries.
     target_high_union_params : list[float]
@@ -181,12 +167,6 @@ class InterpolationConfig(Data):
     vertical_layers_max_centroid_dist: float = field(
         default_factory=lambda: _interpolation_defaults().get("vertical_layers_max_centroid_dist", 25.0)
     )
-    target_low_geodesics_method: GeodesicsMethod = field(
-        default_factory=lambda: GeodesicsMethod(_interpolation_defaults().get("target_low_geodesics_method", "heat_igl"))
-    )
-    target_high_geodesics_method: GeodesicsMethod = field(
-        default_factory=lambda: GeodesicsMethod(_interpolation_defaults().get("target_high_geodesics_method", "heat_igl"))
-    )
     target_high_union_method: UnionMethod = field(
         default_factory=lambda: UnionMethod(_interpolation_defaults().get("target_high_union_method", "min"))
     )
@@ -199,11 +179,6 @@ class InterpolationConfig(Data):
 
     def __post_init__(self) -> None:
         super().__init__()
-        # Convert string enums if needed
-        if isinstance(self.target_low_geodesics_method, str):
-            self.target_low_geodesics_method = GeodesicsMethod(self.target_low_geodesics_method)
-        if isinstance(self.target_high_geodesics_method, str):
-            self.target_high_geodesics_method = GeodesicsMethod(self.target_high_geodesics_method)
         if isinstance(self.target_high_union_method, str):
             self.target_high_union_method = UnionMethod(self.target_high_union_method)
 
@@ -214,8 +189,6 @@ def __data__(self) -> dict[str, Any]:
             "min_layer_height": self.min_layer_height,
             "max_layer_height": self.max_layer_height,
             "vertical_layers_max_centroid_dist": self.vertical_layers_max_centroid_dist,
-            "target_low_geodesics_method": self.target_low_geodesics_method.value,
-            "target_high_geodesics_method": self.target_high_geodesics_method.value,
             "target_high_union_method": self.target_high_union_method.value,
             "target_high_union_params": self.target_high_union_params,
             "uneven_upper_targets_offset": self.uneven_upper_targets_offset,
@@ -231,12 +204,6 @@ def __from_data__(cls, data: dict[str, Any]) -> InterpolationConfig:
             vertical_layers_max_centroid_dist=data.get(
                 "vertical_layers_max_centroid_dist", d.get("vertical_layers_max_centroid_dist", 25.0)
             ),
-            target_low_geodesics_method=data.get(
-                "target_low_geodesics_method", d.get("target_low_geodesics_method", "heat_igl")
-            ),
-            target_high_geodesics_method=data.get(
-                "target_high_geodesics_method", d.get("target_high_geodesics_method", "heat_igl")
-            ),
             target_high_union_method=data.get("target_high_union_method", d.get("target_high_union_method", "min")),
             target_high_union_params=data.get("target_high_union_params", d.get("target_high_union_params", [])),
             uneven_upper_targets_offset=data.get(

src/compas_slicer/data/defaults.toml

@@ -11,8 +11,6 @@ avg_layer_height = 5.0
 min_layer_height = 0.5
 max_layer_height = 10.0
 vertical_layers_max_centroid_dist = 25.0
-target_low_geodesics_method = "heat_cgal"
-target_high_geodesics_method = "heat_cgal"
 target_high_union_method = "min"
 target_high_union_params = []
 uneven_upper_targets_offset = 0.0

src/compas_slicer/pre_processing/interpolation_slicing_preprocessor.py

@@ -63,23 +63,15 @@ def create_compound_targets(self) -> None:
         """Create target_LOW and target_HIGH and compute geodesic distances."""
 
         # --- low target
-        geodesics_method = self.config.target_low_geodesics_method.value
-        method = 'min'  # no other union methods currently supported for lower target
-        params: list[float] = []
-        self.target_LOW = CompoundTarget(self.mesh, 'boundary', 1, self.DATA_PATH,
-                                         union_method=method,
-                                         union_params=params,
-                                         geodesics_method=geodesics_method)
+        self.target_LOW = CompoundTarget(self.mesh, 'boundary', 1, self.DATA_PATH)
 
         # --- high target
-        geodesics_method = self.config.target_high_geodesics_method.value
         method = self.config.target_high_union_method.value
         params = self.config.target_high_union_params
         logger.info(f"Creating target with union type: {method} and params: {params}")
-        self.target_HIGH = CompoundTarget(self.mesh, 'boundary', 2, self.DATA_PATH,
-                                          union_method=method,
-                                          union_params=params,
-                                          geodesics_method=geodesics_method)
+        self.target_HIGH = CompoundTarget(
+            self.mesh, 'boundary', 2, self.DATA_PATH, union_method=method, union_params=params
+        )
 
         # --- uneven boundaries of high target
         self.target_HIGH.offset = self.config.uneven_upper_targets_offset

src/compas_slicer/pre_processing/preprocessing_utils/compound_target.py

@@ -11,14 +11,8 @@
 from numpy.typing import NDArray
 
 import compas_slicer.utilities as utils
-from compas_slicer.pre_processing.preprocessing_utils.geodesics import (
-    get_cgal_HEAT_geodesic_distances,
-    get_custom_HEAT_geodesic_distances,
-    get_igl_EXACT_geodesic_distances,
-    get_igl_HEAT_geodesic_distances,
-)
-
-GeodesicsMethod = Literal['exact_igl', 'heat_igl', 'heat_cgal', 'heat']
+from compas_slicer.pre_processing.preprocessing_utils.geodesics import get_heat_geodesic_distances
+
 UnionMethod = Literal['min', 'smooth', 'chamfer', 'stairs']
 
 
@@ -57,11 +51,6 @@ class CompoundTarget:
     DATA_PATH: str
     has_blend_union: bool
     blend_radius : float
-    geodesics_method: str
-        'heat_cgal'  CGAL heat geodesic distances (recommended)
-        'heat'       custom heat geodesic distances
-    anisotropic_scaling: bool
-        This is not yet implemented
     """
 
     def __init__(
@@ -72,8 +61,6 @@ def __init__(
         DATA_PATH: str,
         union_method: UnionMethod = 'min',
         union_params: list[Any] | None = None,
-        geodesics_method: GeodesicsMethod = 'heat_cgal',
-        anisotropic_scaling: bool = False,
     ) -> None:
 
         if union_params is None:
@@ -89,9 +76,6 @@ def __init__(
         self.union_method = union_method
         self.union_params = union_params
 
-        self.geodesics_method = geodesics_method
-        self.anisotropic_scaling = anisotropic_scaling  # Anisotropic scaling not yet implemented
-
         self.offset = 0
         self.VN = len(list(self.mesh.vertices()))
 
@@ -145,20 +129,9 @@ def compute_geodesic_distances(self) -> None:
         Computes the geodesic distances from each of the target's neighborhoods  to all the mesh vertices.
         Fills in the distances attributes.
         """
-        if self.geodesics_method == 'exact_igl':
-            distances_lists = [get_igl_EXACT_geodesic_distances(self.mesh, vstarts) for vstarts in
-                               self.clustered_vkeys]
-        elif self.geodesics_method == 'heat_igl':
-            distances_lists = [get_igl_HEAT_geodesic_distances(self.mesh, vstarts) for vstarts in
-                               self.clustered_vkeys]
-        elif self.geodesics_method == 'heat_cgal':
-            distances_lists = [get_cgal_HEAT_geodesic_distances(self.mesh, vstarts) for vstarts in
-                               self.clustered_vkeys]
-        elif self.geodesics_method == 'heat':
-            distances_lists = [get_custom_HEAT_geodesic_distances(self.mesh, vstarts, str(self.OUTPUT_PATH)) for vstarts in
-                               self.clustered_vkeys]
-        else:
-            raise ValueError('Unknown geodesics method : ' + self.geodesics_method)
+        distances_lists = [
+            get_heat_geodesic_distances(self.mesh, vstarts) for vstarts in self.clustered_vkeys
+        ]
 
         distances_lists = [list(dl) for dl in distances_lists]  # number_of_boundaries x #V
         self.update_distances_lists(distances_lists)