brain: add arc-aware inset boundary handling

tools/adaptive-isogrid/src/brain/__main__.py

@@ -14,6 +14,7 @@ import numpy as np
 from shapely.geometry import Polygon
 
 from src.atomizer_study import DEFAULT_PARAMS
+from src.shared.arc_utils import typed_segments_to_polyline
 from .density_field import evaluate_density_grid
 from .triangulation import generate_triangulation
 from .pocket_profiles import generate_pockets
@@ -40,6 +41,15 @@ def _merge_params(geometry: Dict[str, Any], params_file: Path | None) -> Dict[st
     return params
 
 
+def _plot_boundary_polyline(geometry: Dict[str, Any], arc_pts: int = 64) -> np.ndarray:
+    typed = geometry.get("outer_boundary_typed")
+    if typed:
+        pts = typed_segments_to_polyline(typed, arc_pts=arc_pts)
+        if len(pts) >= 3:
+            return np.asarray(pts, dtype=float)
+    return np.asarray(geometry["outer_boundary"], dtype=float)
+
+
 def _plot_density(geometry: Dict[str, Any], params: Dict[str, Any], out_path: Path, resolution: float = 3.0) -> None:
     X, Y, eta = evaluate_density_grid(geometry, params, resolution=resolution)
 
@@ -47,7 +57,7 @@ def _plot_density(geometry: Dict[str, Any], params: Dict[str, Any], out_path: Pa
     m = ax.pcolormesh(X, Y, eta, shading="auto", cmap="viridis", vmin=0.0, vmax=1.0)
     fig.colorbar(m, ax=ax, label="Density η")
 
-    outer = np.asarray(geometry["outer_boundary"])
+    outer = _plot_boundary_polyline(geometry)
     ax.plot(np.r_[outer[:, 0], outer[0, 0]], np.r_[outer[:, 1], outer[0, 1]], "w-", lw=1.5)
 
     for hole in geometry.get("holes", []):
@@ -73,7 +83,7 @@ def _plot_triangulation(geometry: Dict[str, Any], triangulation: Dict[str, Any],
     verts = triangulation["vertices"]
     tris = triangulation["triangles"]
 
-    outer = np.asarray(geometry["outer_boundary"])
+    outer = _plot_boundary_polyline(geometry)
     plate_poly = ShapelyPolygon(outer)
     w_frame = (params or {}).get("w_frame", 2.0)
     inner_plate = plate_poly.buffer(-w_frame)
@@ -133,7 +143,7 @@ def _plot_final_profile(geometry, pockets, ribbed_plate, out_path: Path, params:
     """
     from shapely.geometry import Polygon as ShapelyPolygon
 
-    outer = np.asarray(geometry["outer_boundary"])
+    outer = _plot_boundary_polyline(geometry)
     plate_poly = ShapelyPolygon(outer)
     w_frame = (params or {}).get("w_frame", 2.0)
     inner_plate = plate_poly.buffer(-w_frame)

tools/adaptive-isogrid/src/brain/triangulation.py

@@ -14,9 +14,10 @@ a viable machining pocket.
 
 import numpy as np
 from scipy.spatial import Delaunay
-from shapely.geometry import Polygon, Point, MultiPoint
+from shapely.geometry import Polygon, Point, MultiPoint, LinearRing
 from shapely.ops import unary_union
 
+from src.shared.arc_utils import inset_arc, typed_segments_to_polyline, typed_segments_to_sparse
 from .density_field import evaluate_density, density_to_spacing
 
 
@@ -149,58 +150,93 @@ def _adaptive_hex_grid(geometry, params):
 
 
 def _add_boundary_vertices(points, geometry, params, keepout_union):
-    """
+    """Add inset boundary vertices (arc-aware) and keepout boundary points."""
-    Add vertices along the outer boundary and hole keepout boundaries.
-    
-    This ensures triangles conform to boundaries rather than just being
-    clipped. Points are spaced at approximately s_min along boundaries.
-    
-    KEY: Enforce explicit vertices at every corner of the inset boundary.
-    This guarantees no triangle can cross a corner — the Delaunay triangulation
-    is forced to use these corner points as vertices.
-    """
     s_min = params['s_min']
     w_frame = params.get('w_frame', 8.0)
 
     new_pts = list(points)
-
     plate_poly = Polygon(geometry['outer_boundary'])
-    inner_frame = plate_poly.buffer(-w_frame)
+
-    if not inner_frame.is_empty:
+    typed_segments = geometry.get('outer_boundary_typed')
-        # Handle MultiPolygon from buffer on complex shapes
+    if typed_segments:
-        if inner_frame.geom_type == 'MultiPolygon':
+        ring = LinearRing(geometry['outer_boundary'])
-            inner_polys = list(inner_frame.geoms)
+        is_ccw = bool(ring.is_ccw)
+        inset_segments = []
+
+        for seg in typed_segments:
+            stype = seg.get('type', 'line')
+            if stype == 'arc':
+                center_inside = plate_poly.contains(Point(seg['center']))
+                inset_segments.append(inset_arc({**seg, 'center_inside': center_inside}, w_frame))
+                continue
+
+            x1, y1 = seg['start']
+            x2, y2 = seg['end']
+            dx, dy = (x2 - x1), (y2 - y1)
+            ln = np.hypot(dx, dy)
+            if ln < 1e-12:
+                continue
+            nx_l, ny_l = (-dy / ln), (dx / ln)
+            nx, ny = (nx_l, ny_l) if is_ccw else (-nx_l, -ny_l)
+            inset_segments.append({
+                'type': 'line',
+                'start': [x1 + w_frame * nx, y1 + w_frame * ny],
+                'end': [x2 + w_frame * nx, y2 + w_frame * ny],
+            })
+
+        # Sparse points forced into Delaunay (3/arc, 2/line)
+        sparse_pts = typed_segments_to_sparse(inset_segments)
+        new_pts.extend(sparse_pts)
+
+        # Dense inset polygon for containment checks
+        dense = typed_segments_to_polyline(inset_segments, arc_pts=16)
+        if len(dense) >= 3:
+            inner_plate = Polygon(dense)
+            if not inner_plate.is_valid:
+                inner_plate = inner_plate.buffer(0)
+            if inner_plate.is_empty:
+                inner_plate = plate_poly.buffer(-w_frame)
         else:
-            inner_polys = [inner_frame]
+            inner_plate = plate_poly.buffer(-w_frame)
 
-        for inner_poly in inner_polys:
+        # Even spacing along inset boundary (as before)
-            ring = inner_poly.exterior
+        if not inner_plate.is_empty and inner_plate.is_valid:
-
+            if inner_plate.geom_type == 'MultiPolygon':
-            # 1) ENFORCE corner vertices: add every vertex of the inset boundary
+                inner_polys = list(inner_plate.geoms)
-            #    These are the actual corner points — critical for preventing crossovers
+            else:
-            coords = list(ring.coords)[:-1]  # skip closing duplicate
+                inner_polys = [inner_plate]
-            for cx, cy in coords:
+            for inner_poly in inner_polys:
-                new_pts.append([cx, cy])
+                ext = inner_poly.exterior
-
+                length = ext.length
-            # 2) Add evenly spaced points along edges for density
+                n_pts = max(int(length / s_min), 4)
-            length = ring.length
+                for i in range(n_pts):
-            n_pts = max(int(length / s_min), 4)
+                    frac = i / n_pts
-            for i in range(n_pts):
+                    pt = ext.interpolate(frac, normalized=True)
-                frac = i / n_pts
+                    new_pts.append([pt.x, pt.y])
-                pt = ring.interpolate(frac, normalized=True)
+    else:
-                new_pts.append([pt.x, pt.y])
+        inner_plate = plate_poly.buffer(-w_frame)
-
+        if not inner_plate.is_empty:
-            # Also add inner ring vertices (for any holes in the inset boundary)
+            if inner_plate.geom_type == 'MultiPolygon':
-            for interior in inner_poly.interiors:
+                inner_polys = list(inner_plate.geoms)
-                for cx, cy in list(interior.coords)[:-1]:
+            else:
+                inner_polys = [inner_plate]
+            for inner_poly in inner_polys:
+                ext = inner_poly.exterior
+                coords = list(ext.coords)[:-1]
+                for cx, cy in coords:
                     new_pts.append([cx, cy])
+                length = ext.length
+                n_pts = max(int(length / s_min), 4)
+                for i in range(n_pts):
+                    frac = i / n_pts
+                    pt = ext.interpolate(frac, normalized=True)
+                    new_pts.append([pt.x, pt.y])
 
     # Add points along hole keepout boundaries
     if not keepout_union.is_empty:
         geoms = [keepout_union] if keepout_union.geom_type == 'Polygon' else list(keepout_union.geoms)
         for geom in geoms:
             ring = geom.exterior
-            # Enforce corner vertices on keepout boundaries too
             for cx, cy in list(ring.coords)[:-1]:
                 new_pts.append([cx, cy])
             length = ring.length
@@ -210,7 +246,10 @@ def _add_boundary_vertices(points, geometry, params, keepout_union):
                 pt = ring.interpolate(frac, normalized=True)
                 new_pts.append([pt.x, pt.y])
 
-    return np.array(new_pts, dtype=np.float64)
+    if inner_plate.is_empty or not inner_plate.is_valid:
+        inner_plate = plate_poly
+
+    return np.array(new_pts, dtype=np.float64), inner_plate
 
 
 def generate_triangulation(geometry, params, max_refinement_passes=3):
@@ -225,8 +264,8 @@ def generate_triangulation(geometry, params, max_refinement_passes=3):
     if len(grid_pts) < 3:
         return {'vertices': grid_pts, 'triangles': np.empty((0, 3), dtype=int)}
 
-    # Add boundary-conforming vertices
+    # Add boundary-conforming vertices and get inset plate polygon for clipping
-    all_pts = _add_boundary_vertices(grid_pts, geometry, params, keepout_union)
+    all_pts, inner_plate = _add_boundary_vertices(grid_pts, geometry, params, keepout_union)
 
     # Deduplicate close points
     all_pts = np.unique(np.round(all_pts, 4), axis=0)
@@ -240,8 +279,6 @@ def generate_triangulation(geometry, params, max_refinement_passes=3):
 
     # Filter: keep only triangles whose centroid is inside valid region
     plate_poly = Polygon(geometry['outer_boundary'])
-    w_frame = params.get('w_frame', 8.0)
-    inner_plate = plate_poly.buffer(-w_frame)
     if inner_plate.is_empty:
         inner_plate = plate_poly
 

tools/adaptive-isogrid/src/shared/__init__.py

@@ -0,0 +1 @@
+"""Shared geometry utilities."""

tools/adaptive-isogrid/src/shared/arc_utils.py

@@ -0,0 +1,126 @@
+"""Utilities for working with typed boundary segments (line/arc)."""
+
+from __future__ import annotations
+
+import math
+from typing import Any, Dict, List
+
+
+def _as_xy(pt):
+    return [float(pt[0]), float(pt[1])]
+
+
+def _arc_angles(seg: Dict[str, Any]):
+    cx, cy = seg["center"]
+    sx, sy = seg["start"]
+    ex, ey = seg["end"]
+    a0 = math.atan2(sy - cy, sx - cx)
+    a1 = math.atan2(ey - cy, ex - cx)
+    cw = bool(seg.get("clockwise", False))
+    if cw:
+        if a1 > a0:
+            a1 -= 2.0 * math.pi
+    else:
+        if a1 < a0:
+            a1 += 2.0 * math.pi
+    return a0, a1
+
+
+def arc_to_3pt(seg: Dict[str, Any]) -> List[List[float]]:
+    """Sparse arc representation for triangulation constraints."""
+    start = _as_xy(seg["start"])
+    end = _as_xy(seg["end"])
+    mid = _as_xy(seg.get("mid") or arc_to_polyline(seg, n_pts=3)[1])
+    return [start, mid, end]
+
+
+def arc_to_polyline(seg: Dict[str, Any], n_pts: int = 16) -> List[List[float]]:
+    """Discretize an arc segment into polyline points (inclusive endpoints)."""
+    cx, cy = _as_xy(seg["center"])
+    radius = float(seg["radius"])
+    if radius <= 0.0:
+        return [_as_xy(seg["start"]), _as_xy(seg["end"])]
+
+    n_pts = max(int(n_pts), 2)
+    a0, a1 = _arc_angles(seg)
+
+    # Handle full-circle arc representation (start == end): use one revolution.
+    s = _as_xy(seg["start"])
+    e = _as_xy(seg["end"])
+    if abs(s[0] - e[0]) < 1e-9 and abs(s[1] - e[1]) < 1e-9:
+        a1 = a0 - 2.0 * math.pi if bool(seg.get("clockwise", False)) else a0 + 2.0 * math.pi
+
+    out = []
+    for i in range(n_pts):
+        t = i / (n_pts - 1)
+        a = a0 + t * (a1 - a0)
+        out.append([cx + radius * math.cos(a), cy + radius * math.sin(a)])
+    return out
+
+
+def inset_arc(seg: Dict[str, Any], offset: float) -> Dict[str, Any]:
+    """Return a parallel offset arc.
+
+    Uses seg['center_inside'] when provided to choose convex vs concave behavior.
+    """
+    cx, cy = _as_xy(seg["center"])
+    radius = float(seg["radius"])
+    center_inside = bool(seg.get("center_inside", True))
+    new_radius = radius - offset if center_inside else radius + offset
+    new_radius = max(new_radius, 1e-6)
+
+    def scale_point(pt):
+        px, py = _as_xy(pt)
+        vx, vy = px - cx, py - cy
+        vn = math.hypot(vx, vy)
+        if vn < 1e-12:
+            return [cx + new_radius, cy]
+        k = new_radius / vn
+        return [cx + vx * k, cy + vy * k]
+
+    new_seg = dict(seg)
+    new_seg["radius"] = new_radius
+    new_seg["start"] = scale_point(seg["start"])
+    new_seg["end"] = scale_point(seg["end"])
+    if "mid" in seg:
+        new_seg["mid"] = scale_point(seg["mid"])
+    else:
+        new_seg["mid"] = arc_to_polyline(new_seg, n_pts=3)[1]
+    return new_seg
+
+
+def typed_segments_to_polyline(segments, arc_pts: int = 16) -> List[List[float]]:
+    """Convert typed segments to a dense boundary polyline."""
+    out: List[List[float]] = []
+    for seg in segments or []:
+        stype = seg.get("type", "line")
+        if stype == "arc":
+            pts = arc_to_polyline(seg, n_pts=arc_pts)
+        else:
+            pts = [_as_xy(seg["start"]), _as_xy(seg["end"])]
+        if out and pts:
+            if abs(out[-1][0] - pts[0][0]) < 1e-9 and abs(out[-1][1] - pts[0][1]) < 1e-9:
+                out.extend(pts[1:])
+            else:
+                out.extend(pts)
+        else:
+            out.extend(pts)
+    return out
+
+
+def typed_segments_to_sparse(segments) -> List[List[float]]:
+    """Convert typed segments to sparse points (3 for arcs, 2 for lines)."""
+    out: List[List[float]] = []
+    for seg in segments or []:
+        if seg.get("type") == "arc":
+            pts = arc_to_3pt(seg)
+        else:
+            pts = [_as_xy(seg["start"]), _as_xy(seg["end"])]
+        if out and pts:
+            if abs(out[-1][0] - pts[0][0]) < 1e-9 and abs(out[-1][1] - pts[0][1]) < 1e-9:
+                out.extend(pts[1:])
+            else:
+                out.extend(pts)
+        else:
+            out.extend(pts)
+    return out