From 7d5bd33bb55957c61f9a6b68e6c845e97e077c2e Mon Sep 17 00:00:00 2001 From: Antoine Date: Tue, 17 Feb 2026 14:05:28 +0000 Subject: [PATCH] brain: add arc-aware inset boundary handling --- tools/adaptive-isogrid/src/brain/__main__.py | 16 ++- .../src/brain/triangulation.py | 125 +++++++++++------ tools/adaptive-isogrid/src/shared/__init__.py | 1 + .../adaptive-isogrid/src/shared/arc_utils.py | 126 ++++++++++++++++++ 4 files changed, 221 insertions(+), 47 deletions(-) create mode 100644 tools/adaptive-isogrid/src/shared/__init__.py create mode 100644 tools/adaptive-isogrid/src/shared/arc_utils.py diff --git a/tools/adaptive-isogrid/src/brain/__main__.py b/tools/adaptive-isogrid/src/brain/__main__.py index e9ca75be..267bd877 100644 --- a/tools/adaptive-isogrid/src/brain/__main__.py +++ b/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) diff --git a/tools/adaptive-isogrid/src/brain/triangulation.py b/tools/adaptive-isogrid/src/brain/triangulation.py index 532015c7..0a607a7c 100644 --- a/tools/adaptive-isogrid/src/brain/triangulation.py +++ b/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 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. - """ + """Add inset boundary vertices (arc-aware) and keepout boundary points.""" 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: - # Handle MultiPolygon from buffer on complex shapes - if inner_frame.geom_type == 'MultiPolygon': - inner_polys = list(inner_frame.geoms) + + typed_segments = geometry.get('outer_boundary_typed') + if typed_segments: + ring = LinearRing(geometry['outer_boundary']) + 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: - ring = inner_poly.exterior - - # 1) ENFORCE corner vertices: add every vertex of the inset boundary - # These are the actual corner points — critical for preventing crossovers - coords = list(ring.coords)[:-1] # skip closing duplicate - for cx, cy in coords: - new_pts.append([cx, cy]) - - # 2) Add evenly spaced points along edges for density - length = ring.length - n_pts = max(int(length / s_min), 4) - for i in range(n_pts): - frac = i / n_pts - pt = ring.interpolate(frac, normalized=True) - new_pts.append([pt.x, pt.y]) - - # Also add inner ring vertices (for any holes in the inset boundary) - for interior in inner_poly.interiors: - for cx, cy in list(interior.coords)[:-1]: + # Even spacing along inset boundary (as before) + if not inner_plate.is_empty and inner_plate.is_valid: + if inner_plate.geom_type == 'MultiPolygon': + inner_polys = list(inner_plate.geoms) + else: + inner_polys = [inner_plate] + for inner_poly in inner_polys: + ext = inner_poly.exterior + 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]) + else: + inner_plate = plate_poly.buffer(-w_frame) + if not inner_plate.is_empty: + if inner_plate.geom_type == 'MultiPolygon': + inner_polys = list(inner_plate.geoms) + 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 - all_pts = _add_boundary_vertices(grid_pts, geometry, params, keepout_union) + # Add boundary-conforming vertices and get inset plate polygon for clipping + 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 diff --git a/tools/adaptive-isogrid/src/shared/__init__.py b/tools/adaptive-isogrid/src/shared/__init__.py new file mode 100644 index 00000000..790e69dd --- /dev/null +++ b/tools/adaptive-isogrid/src/shared/__init__.py @@ -0,0 +1 @@ +"""Shared geometry utilities.""" diff --git a/tools/adaptive-isogrid/src/shared/arc_utils.py b/tools/adaptive-isogrid/src/shared/arc_utils.py new file mode 100644 index 00000000..c49f0495 --- /dev/null +++ b/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