refactor: use compas_cgal.isolines for contour extraction

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.

docs/index.md

@@ -19,7 +19,7 @@ COMPAS Slicer was developed at [ETH Zurich](https://ethz.ch/) by:
 
 - **[Andrei Jipa](https://github.com/stratocaster)** - [Gramazio Kohler Research](https://gramaziokohler.arch.ethz.ch/)
 
-- **[Jelle Feringa](https://github.com/jf---)** - [Gramazio Kohler Research](https://terrestrial.construction)
+- **[Jelle Feringa](https://github.com/jf---)** - [Terrestrial](https://terrestrial.construction)
 
 The package emerged from research on non-planar 3D printing and robotic fabrication at the Institute of Technology in Architecture.
 

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/__main__.py

@@ -3,7 +3,7 @@
 
 import compas_slicer
 
-if __name__ == '__main__':
-    logger.info(f'COMPAS: {compas.__version__}')
-    logger.info(f'COMPAS Slicer: {compas_slicer.__version__}')
-    logger.info('Awesome! Your installation worked! :)')
+if __name__ == "__main__":
+    logger.info(f"COMPAS: {compas.__version__}")
+    logger.info(f"COMPAS Slicer: {compas_slicer.__version__}")
+    logger.info("Awesome! Your installation worked! :)")

src/compas_slicer/_numpy_ops.py

@@ -143,9 +143,7 @@ def face_gradient_from_scalar_field(
     cross2 = np.cross(v1 - v0, face_normals)  # (F, 3)
 
     # Compute gradient
-    grad = (
-        (u1 - u0)[:, np.newaxis] * cross1 + (u2 - u0)[:, np.newaxis] * cross2
-    ) / (2 * face_areas[:, np.newaxis])
+    grad = ((u1 - u0)[:, np.newaxis] * cross1 + (u2 - u0)[:, np.newaxis] * cross2) / (2 * face_areas[:, np.newaxis])
 
     return grad
 
@@ -187,9 +185,9 @@ def per_vertex_divergence(
     e2 = v0 - v1  # edge opposite to v2
 
     # Compute dot products with gradient
-    dot0 = np.einsum('ij,ij->i', X, e0)  # (F,)
-    dot1 = np.einsum('ij,ij->i', X, e1)  # (F,)
-    dot2 = np.einsum('ij,ij->i', X, e2)  # (F,)
+    dot0 = np.einsum("ij,ij->i", X, e0)  # (F,)
+    dot1 = np.einsum("ij,ij->i", X, e1)  # (F,)
+    dot2 = np.einsum("ij,ij->i", X, e2)  # (F,)
 
     # Cotangent contributions (cotans[f, i] is cotan of angle at vertex i)
     # For vertex i: contrib = cotan[k] * dot(X, e_i) + cotan[j] * dot(X, -e_k)

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/geometry/__init__.py

@@ -5,4 +5,4 @@
 from .print_point import *  # noqa: F401 E402 F403
 from .printpoints_collection import *  # noqa: F401 E402 F403
 
-__all__ = [name for name in dir() if not name.startswith('_')]
+__all__ = [name for name in dir() if not name.startswith("_")]

src/compas_slicer/geometry/path.py

@@ -49,10 +49,7 @@ def __from_data__(cls, data: dict[str, Any]) -> Path:
         points_data = data["points"]
         # Handle both list format and legacy dict format
         if isinstance(points_data, dict):
-            pts = [
-                Point.__from_data__(points_data[key])
-                for key in sorted(points_data.keys(), key=lambda x: int(x))
-            ]
+            pts = [Point.__from_data__(points_data[key]) for key in sorted(points_data.keys(), key=lambda x: int(x))]
         else:
             pts = [Point.__from_data__(p) for p in points_data]
         return cls(points=pts, is_closed=data["is_closed"])

src/compas_slicer/post_processing/__init__.py

@@ -21,4 +21,4 @@
 from .unify_paths_orientation import *  # noqa: F401 E402 F403
 from .zig_zag_open_paths import *  # noqa: F401 E402 F403
 
-__all__ = [name for name in dir() if not name.startswith('_')]
+__all__ = [name for name in dir() if not name.startswith("_")]

src/compas_slicer/post_processing/generate_brim.py

@@ -14,6 +14,7 @@
 try:
     from compas_cgal.straight_skeleton_2 import offset_polygon as _cgal_offset
     from compas_cgal.straight_skeleton_2 import offset_polygon_with_holes as _cgal_offset_with_holes
+
     _USE_CGAL = True
 except ImportError:
     _cgal_offset = None
@@ -23,7 +24,7 @@
     from compas_slicer.slicers import BaseSlicer
 
 
-__all__ = ['generate_brim', 'offset_polygon', 'offset_polygon_with_holes']
+__all__ = ["generate_brim", "offset_polygon", "offset_polygon_with_holes"]
 
 
 def _offset_polygon_cgal(points: list[Point], offset: float, z: float) -> list[Point]:
@@ -83,16 +84,12 @@ def _offset_polygon_pyclipper(points: list[Point], offset: float, z: float) -> l
     import pyclipper
     from pyclipper import scale_from_clipper, scale_to_clipper
 
-    SCALING_FACTOR = 2 ** 32
+    SCALING_FACTOR = 2**32
 
     xy_coords = [[p[0], p[1]] for p in points]
 
     pco = pyclipper.PyclipperOffset()
-    pco.AddPath(
-        scale_to_clipper(xy_coords, SCALING_FACTOR),
-        pyclipper.JT_MITER,
-        pyclipper.ET_CLOSEDPOLYGON
-    )
+    pco.AddPath(scale_to_clipper(xy_coords, SCALING_FACTOR), pyclipper.JT_MITER, pyclipper.ET_CLOSEDPOLYGON)
 
     result = scale_from_clipper(pco.Execute(offset * SCALING_FACTOR), SCALING_FACTOR)
 
@@ -132,10 +129,7 @@ def offset_polygon(points: list[Point], offset: float, z: float) -> list[Point]:
 
 
 def offset_polygon_with_holes(
-    outer: list[Point],
-    holes: list[list[Point]],
-    offset: float,
-    z: float
+    outer: list[Point], holes: list[list[Point]], offset: float, z: float
 ) -> list[tuple[list[Point], list[list[Point]]]]:
     """Offset a polygon with holes using CGAL straight skeleton.
 
@@ -224,7 +218,7 @@ def generate_brim(slicer: BaseSlicer, layer_width: float, number_of_brim_offsets
         has_vertical_layers = False
 
     if len(paths_to_offset) == 0:
-        raise ValueError('Brim generator did not find any path on the base. Please check the paths of your slicer.')
+        raise ValueError("Brim generator did not find any path on the base. Please check the paths of your slicer.")
 
     # (2) --- create new empty brim_layer
     brim_layer = Layer(paths=[])