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

"""
Isogrid triangulation — generates a proper isogrid rib pattern.

Uses Shewchuk's Triangle library (constrained Delaunay + quality refinement)
for boundary-conforming, clean triangulations.

Strategy:
1. Build PSLG (Planar Straight Line Graph):
   - Outer contour = inset boundary (w_frame offset inward)
   - Holes = keepout circles (d_keep offset outward from holes)
2. Triangulate with quality constraints:
   - Minimum angle ~25° (no slivers)
   - Maximum area varies by density heatmap (smaller near holes, larger away)
3. Result: clean mesh that fills 100% of valid area, with density-graded
   triangle sizes and edges parallel to boundary/hole contours.
"""

import numpy as np
import triangle as tr
from shapely.geometry import Polygon, Point, LinearRing, MultiPolygon
from shapely.geometry.base import BaseGeometry
from shapely.ops import unary_union

from src.shared.arc_utils import inset_arc, typed_segments_to_polyline
from .density_field import evaluate_density, density_to_spacing


# ---------------------------------------------------------------------------
# Helpers
# ---------------------------------------------------------------------------

def _geometry_to_list(geom: BaseGeometry) -> list:
    if geom.is_empty:
        return []
    if geom.geom_type == 'Polygon':
        return [geom]
    if geom.geom_type == 'MultiPolygon':
        return list(geom.geoms)
    return list(getattr(geom, 'geoms', []))


def _ring_to_segments(coords: np.ndarray, start_idx: int):
    """Convert a ring (Nx2 array, NOT closed) to vertex array + segment index pairs."""
    n = len(coords)
    segments = []
    for i in range(n):
        segments.append([start_idx + i, start_idx + (i + 1) % n])
    return segments


def _sample_ring(ring, spacing: float) -> np.ndarray:
    """Sample points along a Shapely ring at given spacing, returning Nx2 array (not closed)."""
    length = float(ring.length)
    if length < 1e-9:
        return np.empty((0, 2), dtype=np.float64)
    n = max(int(np.ceil(length / max(spacing, 1e-3))), 8)
    pts = []
    for i in range(n):
        p = ring.interpolate(i / n, normalized=True)
        pts.append([p.x, p.y])
    return np.array(pts, dtype=np.float64)


# ---------------------------------------------------------------------------
# Step 1 — Build the inset plate polygon (frame inner edge)
# ---------------------------------------------------------------------------

def _build_inner_plate(geometry, params) -> Polygon:
    """Offset sandbox boundary inward by w_frame."""
    w_frame = float(params.get('w_frame', 8.0))
    plate_poly = Polygon(geometry['outer_boundary'])

    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))
            else:
                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],
                })

        dense = typed_segments_to_polyline(inset_segments, arc_pts=32)
        if len(dense) >= 3:
            inner_plate = Polygon(dense)
            if not inner_plate.is_valid:
                inner_plate = inner_plate.buffer(0)
            if not inner_plate.is_empty:
                return inner_plate

    inner_plate = plate_poly.buffer(-w_frame)
    if inner_plate.is_empty or not inner_plate.is_valid:
        inner_plate = plate_poly
    return inner_plate


# ---------------------------------------------------------------------------
# Step 2 — Build hole keepout polygons
# ---------------------------------------------------------------------------

def _build_keepouts(geometry, params) -> list:
    """Build list of keepout Polygons (one per hole)."""
    d_keep = float(params['d_keep'])
    keepouts = []

    for hole in geometry.get('holes', []):
        diameter = float(hole.get('diameter', 10.0) or 10.0)
        if hole.get('is_circular', False) and 'center' in hole:
            cx, cy = hole['center']
            hole_radius = diameter / 2.0
            keepout_radius = hole_radius + d_keep * hole_radius
            keepout = Point(cx, cy).buffer(keepout_radius, resolution=32)
            keepouts.append(keepout)
        else:
            hb = np.asarray(hole.get('boundary', []), dtype=float)
            if len(hb) < 3:
                continue
            ring = LinearRing(hb)
            keepout = Polygon(ring).buffer(max(d_keep, 0.0))
            if not keepout.is_empty:
                keepouts.append(keepout)

    return keepouts


# ---------------------------------------------------------------------------
# Step 3 — Build PSLG for Triangle library
# ---------------------------------------------------------------------------

def _build_pslg(inner_plate: Polygon, keepouts: list, boundary_spacing: float):
    """
    Build a Planar Straight Line Graph (PSLG) for the Triangle library.

    Returns dict with:
      'vertices': Nx2 array
      'segments': Mx2 array of vertex index pairs
      'holes': Hx2 array of hole marker points
    """
    vertices = []
    segments = []
    hole_markers = []

    # Outer boundary of the inner plate
    # Handle MultiPolygon (take largest piece)
    polys = _geometry_to_list(inner_plate)
    if not polys:
        return None
    main_poly = max(polys, key=lambda p: p.area)

    # Sample the outer ring
    outer_pts = _sample_ring(main_poly.exterior, boundary_spacing)
    if len(outer_pts) < 3:
        return None

    start_idx = 0
    vertices.extend(outer_pts.tolist())
    segs = _ring_to_segments(outer_pts, start_idx)
    segments.extend(segs)
    start_idx += len(outer_pts)

    # Handle any interior rings of the inner plate (shouldn't happen normally)
    for interior_ring in main_poly.interiors:
        ring_pts = _sample_ring(interior_ring, boundary_spacing)
        if len(ring_pts) < 3:
            continue
        vertices.extend(ring_pts.tolist())
        segs = _ring_to_segments(ring_pts, start_idx)
        segments.extend(segs)
        # Hole marker inside this interior ring
        centroid = Polygon(interior_ring).centroid
        hole_markers.append([centroid.x, centroid.y])
        start_idx += len(ring_pts)

    # Keepout holes — clip to inner plate
    for keepout in keepouts:
        clipped = inner_plate.intersection(keepout)
        if clipped.is_empty:
            continue
        for kpoly in _geometry_to_list(clipped):
            if kpoly.area < 1.0:
                continue
            ring_pts = _sample_ring(kpoly.exterior, boundary_spacing)
            if len(ring_pts) < 3:
                continue
            vertices.extend(ring_pts.tolist())
            segs = _ring_to_segments(ring_pts, start_idx)
            segments.extend(segs)
            # Hole marker inside the keepout
            centroid = kpoly.centroid
            hole_markers.append([centroid.x, centroid.y])
            start_idx += len(ring_pts)

    return {
        'vertices': np.array(vertices, dtype=np.float64),
        'segments': np.array(segments, dtype=np.int32),
        'holes': np.array(hole_markers, dtype=np.float64) if hole_markers else np.empty((0, 2)),
    }


# ---------------------------------------------------------------------------
# Step 4 — Compute area constraint from density field
# ---------------------------------------------------------------------------

def _max_area_for_spacing(spacing: float) -> float:
    """Area of an equilateral triangle with given edge length."""
    return (np.sqrt(3) / 4.0) * spacing ** 2


def _refine_with_density(tri_result, geometry, params):
    """
    Iterative refinement: check each triangle's area against the density-based
    max area at its centroid. If too large, re-triangulate with tighter constraint.

    Uses Triangle's region-based area constraints via the 'a' switch.
    """
    s_min = float(params['s_min'])
    s_max = float(params['s_max'])

    verts = tri_result['vertices']
    tris = tri_result['triangles']

    # Compute max allowed area for each triangle based on density at centroid
    max_areas = np.empty(len(tris))
    for i, t in enumerate(tris):
        cx = np.mean(verts[t, 0])
        cy = np.mean(verts[t, 1])
        eta = evaluate_density(cx, cy, geometry, params)
        target_spacing = density_to_spacing(eta, params)
        max_areas[i] = _max_area_for_spacing(target_spacing)

    # Use the global maximum area as a single constraint
    # (Triangle library doesn't support per-triangle area easily in basic mode)
    # Instead, use a moderate constraint and let the quality refinement handle it
    global_max_area = _max_area_for_spacing(s_max)
    return global_max_area


# ---------------------------------------------------------------------------
# Main entry point
# ---------------------------------------------------------------------------

def generate_triangulation(geometry, params, max_refinement_passes=3):
    """
    Generate isogrid triangulation using Triangle library.

    Returns dict with 'vertices' (Nx2) and 'triangles' (Mx3).
    """
    s_min = float(params['s_min'])
    s_max = float(params['s_max'])

    # Step 1: Build inner plate
    inner_plate = _build_inner_plate(geometry, params)
    if inner_plate.is_empty:
        return {'vertices': np.empty((0, 2)), 'triangles': np.empty((0, 3), dtype=int)}

    # Step 2: Build keepouts
    keepouts = _build_keepouts(geometry, params)
    keepout_union = unary_union(keepouts) if keepouts else Polygon()

    # Step 3: Build PSLG
    # Boundary sampling at intermediate spacing for clean boundary conformance
    boundary_spacing = max(s_min, min(s_max * 0.5, 30.0))
    pslg = _build_pslg(inner_plate, keepouts, boundary_spacing)

    if pslg is None or len(pslg['vertices']) < 3:
        return {'vertices': np.empty((0, 2)), 'triangles': np.empty((0, 3), dtype=int)}

    # Step 4: Initial triangulation with coarse area constraint
    # 'p' = use PSLG (segments + holes)
    # 'q25' = minimum angle 25° (no slivers)
    # 'a{area}' = max triangle area (coarse — s_max)
    max_area = _max_area_for_spacing(s_max)
    tri_switches = f'pq25a{max_area:.1f}'

    tri_input = {
        'vertices': pslg['vertices'],
        'segments': pslg['segments'],
    }
    if len(pslg['holes']) > 0:
        tri_input['holes'] = pslg['holes']

    try:
        result = tr.triangulate(tri_input, tri_switches)
    except Exception:
        return {'vertices': np.empty((0, 2)), 'triangles': np.empty((0, 3), dtype=int)}

    verts = result.get('vertices', np.empty((0, 2)))
    tris = result.get('triangles', np.empty((0, 3), dtype=int))

    if len(tris) == 0:
        return {'vertices': verts, 'triangles': tris}

    # Step 5: Density-based refinement using per-triangle area constraints
    # Triangle's 'r' switch refines an existing triangulation.
    # We set 'triangle_max_area' per triangle based on density at centroid.
    for pass_num in range(max_refinement_passes):
        # Compute per-triangle max area from density field
        per_tri_max = np.empty(len(tris))
        needs_refinement = False

        for i, t in enumerate(tris):
            cx = np.mean(verts[t, 0])
            cy = np.mean(verts[t, 1])
            eta = evaluate_density(cx, cy, geometry, params)
            target_spacing = density_to_spacing(eta, params)
            per_tri_max[i] = _max_area_for_spacing(target_spacing)

            # Check if this triangle actually needs refinement
            v0, v1, v2 = verts[t[0]], verts[t[1]], verts[t[2]]
            actual_area = 0.5 * abs(
                (v1[0] - v0[0]) * (v2[1] - v0[1]) -
                (v2[0] - v0[0]) * (v1[1] - v0[1])
            )
            if actual_area > per_tri_max[i] * 1.3:
                needs_refinement = True

        if not needs_refinement:
            break

        # Build refinement input with per-triangle area constraints
        refine_input = {
            'vertices': verts,
            'triangles': tris,
            'segments': pslg['segments'] if len(pslg['segments']) <= len(verts) else np.empty((0, 2), dtype=np.int32),
            'triangle_max_area': per_tri_max,
        }
        if len(pslg['holes']) > 0:
            refine_input['holes'] = pslg['holes']

        try:
            result = tr.triangulate(refine_input, 'rpq25')
            verts = result.get('vertices', verts)
            tris = result.get('triangles', tris)
        except Exception:
            break

    # Step 6: Filter degenerate / tiny triangles
    min_area_filter = float(params.get('min_triangle_area', 20.0))
    if len(tris) > 0:
        v0 = verts[tris[:, 0]]
        v1 = verts[tris[:, 1]]
        v2 = verts[tris[:, 2]]
        areas = 0.5 * np.abs(
            (v1[:, 0] - v0[:, 0]) * (v2[:, 1] - v0[:, 1]) -
            (v2[:, 0] - v0[:, 0]) * (v1[:, 1] - v0[:, 1])
        )
        tris = tris[areas >= min_area_filter]

    return {'vertices': verts, 'triangles': tris}