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

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

Strategy: lay a regular equilateral triangle grid, then adapt it:
1. Generate hex-packed vertex grid at base spacing
2. Scale spacing locally based on density field (growing away from features)
3. Delaunay triangulate the point set
4. Clip to plate boundary, exclude hole keepouts
5. Result: well-proportioned triangles everywhere, denser near holes/edges

This is NOT a mesh — it's an isogrid layout. Every triangle should be
a viable machining pocket.
"""

import numpy as np
from scipy.spatial import Delaunay
from shapely.geometry import Polygon, Point, LinearRing
from shapely.geometry.base import BaseGeometry
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


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


def _add_boundary_layer_seed_points(points, geometry, params, inner_plate, keepout_union):
    """Add a row of points just inside the inset boundary to align boundary triangles."""
    offset = float(params.get('boundary_layer_offset', 0.0) or 0.0)
    if offset <= 0.0:
        offset = max(params.get('s_min', 10.0) * 0.6, 0.5)

    layer_geom = inner_plate.buffer(-offset)
    if layer_geom.is_empty:
        return points

    boundary_pts = []
    for poly in _geometry_to_list(layer_geom):
        ring = poly.exterior
        length = float(ring.length)
        if length <= 1e-9:
            continue
        n_pts = max(int(np.ceil(length / max(params.get('s_min', 10.0), 1e-3))), 8)
        for i in range(n_pts):
            frac = i / n_pts
            p = ring.interpolate(frac, normalized=True)
            if not inner_plate.buffer(1e-6).covers(p):
                continue
            if not keepout_union.is_empty and keepout_union.buffer(1e-6).covers(p):
                continue
            boundary_pts.append([p.x, p.y])

    if boundary_pts:
        return np.vstack([points, np.asarray(boundary_pts, dtype=np.float64)])
    return points


def _np(pt):
    return np.asarray([float(pt[0]), float(pt[1])], dtype=float)


def _hole_keepout_seed_points(geometry, params):
    """Return sparse keepout seed points: 3 for circular holes, coarse ring for others."""
    d_keep = float(params['d_keep'])
    pts = []
    for hole in geometry.get('holes', []):
        if hole.get('is_circular', False) and 'center' in hole:
            cx, cy = hole['center']
            diameter = float(hole.get('diameter', 10.0) or 10.0)
            hole_radius = diameter / 2.0
            keepout_radius = hole_radius * (1.0 + d_keep)
            for k in range(3):
                a = (2.0 * np.pi * k) / 3.0
                pts.append([cx + keepout_radius * np.cos(a), cy + keepout_radius * np.sin(a)])
            continue

        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 keepout.is_empty:
            continue
        for geom in _geometry_to_list(keepout):
            ext = geom.exterior
            n = max(int(np.ceil(ext.length / max(params.get('s_min', 10.0), 1e-3))), 6)
            for i in range(n):
                p = ext.interpolate(i / n, normalized=True)
                pts.append([p.x, p.y])

    if not pts:
        return np.empty((0, 2), dtype=np.float64)
    return np.asarray(pts, dtype=np.float64)


def _generate_hex_grid(bbox, base_spacing):
    """
    Generate a regular hexagonal-packed point grid.
    
    This creates the classic isogrid vertex layout: rows of points
    offset by half-spacing, producing equilateral triangles when
    Delaunay triangulated.
    
    Returns ndarray of shape (N, 2).
    """
    x_min, y_min, x_max, y_max = bbox
    # Padding
    pad = base_spacing
    x_min -= pad
    y_min -= pad
    x_max += pad
    y_max += pad

    row_height = base_spacing * np.sqrt(3) / 2.0
    points = []
    row = 0
    y = y_min
    while y <= y_max:
        # Alternate rows offset by half spacing
        x_offset = (base_spacing / 2.0) if (row % 2) else 0.0
        x = x_min + x_offset
        while x <= x_max:
            points.append([x, y])
            x += base_spacing
        y += row_height
        row += 1

    return np.array(points, dtype=np.float64)


def _adaptive_hex_grid(geometry, params):
    """
    Generate density-adaptive hex grid.
    
    Start with a coarse grid at s_max spacing, then iteratively add
    points in high-density regions until local spacing matches the
    density-driven target.
    """
    outer = np.array(geometry['outer_boundary'])
    x_min, y_min = outer.min(axis=0)
    x_max, y_max = outer.max(axis=0)
    bbox = (x_min, y_min, x_max, y_max)

    s_min = params['s_min']
    s_max = params['s_max']
    plate_poly = Polygon(geometry['outer_boundary'])

    # Build hole keepout polygons
    d_keep = params['d_keep']
    keepout_polys = []
    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
            boss_radius = hole_radius + d_keep * hole_radius
            keepout = Point(cx, cy).buffer(boss_radius, resolution=16)
            keepout_polys.append(keepout)

    keepout_union = unary_union(keepout_polys) if keepout_polys else Polygon()

    # Valid region = plate minus keepouts
    valid_region = plate_poly.difference(keepout_union)

    # Inset valid region by frame width to keep grid points inside the frame
    w_frame = params.get('w_frame', 8.0)
    inner_region = plate_poly.buffer(-w_frame)
    if inner_region.is_empty:
        inner_region = plate_poly

    # Generate base grid at s_max (coarsest)
    base_pts = _generate_hex_grid(bbox, s_max)

    # Filter to plate interior (outside keepouts, inside frame)
    valid_pts = []
    for pt in base_pts:
        p = Point(pt[0], pt[1])
        if inner_region.contains(p) and not keepout_union.contains(p):
            valid_pts.append(pt)

    valid_pts = np.array(valid_pts, dtype=np.float64) if valid_pts else np.empty((0, 2))

    if s_min >= s_max * 0.95:
        # Uniform mode — no refinement needed
        return valid_pts, valid_region, keepout_union

    # Adaptive refinement: add denser points near high-density areas
    # Generate a fine grid at s_min
    fine_pts = _generate_hex_grid(bbox, s_min)

    # For each fine point, check if local density warrants adding it
    extra_pts = []
    for pt in fine_pts:
        p = Point(pt[0], pt[1])
        if not inner_region.contains(p) or keepout_union.contains(p):
            continue

        # What spacing does the density field want here?
        eta = evaluate_density(pt[0], pt[1], geometry, params)
        target_spacing = density_to_spacing(eta, params)

        # If target spacing < s_max * 0.8, this is a refinement zone
        if target_spacing < s_max * 0.8:
            # Check distance to nearest existing point
            if len(valid_pts) > 0:
                dists = np.sqrt(np.sum((valid_pts - pt)**2, axis=1))
                min_dist = np.min(dists)
                # Add point if nearest neighbor is farther than target spacing
                if min_dist > target_spacing * 0.85:
                    extra_pts.append(pt)
                    # Also add to valid_pts for future distance checks
                    valid_pts = np.vstack([valid_pts, pt])

    if extra_pts:
        all_pts = np.vstack([valid_pts, np.array(extra_pts)])
        # Deduplicate (shouldn't be needed but safe)
        all_pts = np.unique(np.round(all_pts, 6), axis=0)
    else:
        all_pts = valid_pts

    return all_pts, valid_region, keepout_union


def _add_boundary_vertices(points, geometry, params, keepout_union):
    """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'])

    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_plate = plate_poly.buffer(-w_frame)

        # Even spacing along inset boundary (as before)
        if not inner_plate.is_empty and inner_plate.is_valid:
            for inner_poly in _geometry_to_list(inner_plate):
                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:
        # v1 geometry fallback: inset from polygon and force all inset ring vertices.
        inner_plate = plate_poly.buffer(-w_frame)
        if not inner_plate.is_empty:
            for inner_poly in _geometry_to_list(inner_plate):
                ext = inner_poly.exterior
                # Always keep true inset corners/vertices from the offset polygon.
                coords = list(ext.coords)[:-1]
                for cx, cy in coords:
                    new_pts.append([cx, cy])
                # Plus evenly spaced boundary points for perimeter alignment.
                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 sparse points for hole keepouts (3 points per circular keepout).
    keepout_seeds = _hole_keepout_seed_points(geometry, params)
    if len(keepout_seeds) > 0:
        new_pts.extend(keepout_seeds.tolist())

    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):
    """
    Generate isogrid triangulation: hex-packed grid + Delaunay + clip.
    
    Returns dict with 'vertices' (ndarray Nx2) and 'triangles' (ndarray Mx3).
    """
    # Generate adaptive point grid
    grid_pts, valid_region, keepout_union = _adaptive_hex_grid(geometry, params)

    if len(grid_pts) < 3:
        return {'vertices': grid_pts, 'triangles': np.empty((0, 3), dtype=int)}

    # Add boundary-conforming vertices and get inset plate polygon for clipping
    all_pts, inner_plate = _add_boundary_vertices(grid_pts, geometry, params, keepout_union)

    # Add structured boundary-layer seed row along straight edges
    plate_poly = Polygon(geometry['outer_boundary'])
    all_pts = _add_boundary_layer_seed_points(all_pts, geometry, params, inner_plate, keepout_union)

    # Deduplicate close points
    all_pts = np.unique(np.round(all_pts, 4), axis=0)

    if len(all_pts) < 3:
        return {'vertices': all_pts, 'triangles': np.empty((0, 3), dtype=int)}

    # Delaunay triangulation of the point set
    tri = Delaunay(all_pts)
    triangles = tri.simplices

    # Filter: keep only triangles whose centroid is inside valid region
    plate_poly = Polygon(geometry['outer_boundary'])
    if inner_plate.is_empty:
        inner_plate = plate_poly

    inner_valid_region = inner_plate
    if not keepout_union.is_empty:
        inner_valid_region = inner_plate.difference(keepout_union)

    keep = []
    for i, t in enumerate(triangles):
        cx = np.mean(all_pts[t, 0])
        cy = np.mean(all_pts[t, 1])
        centroid = Point(cx, cy)
        # Keep if centroid is inside plate frame and outside keepouts,
        # AND all 3 vertices are inside the plate boundary (no crossovers)
        if not inner_plate.contains(centroid):
            continue
        if keepout_union.contains(centroid):
            continue
        all_inside = True
        for vi in t:
            if not plate_poly.contains(Point(all_pts[vi])):
                # Allow small tolerance (vertex on boundary is OK)
                if not plate_poly.buffer(0.5).contains(Point(all_pts[vi])):
                    all_inside = False
                    break
        if not all_inside:
            continue

        tri_poly = Polygon(all_pts[t])
        if tri_poly.is_empty or not tri_poly.is_valid:
            continue
        if not inner_valid_region.buffer(1e-6).covers(tri_poly):
            continue
        keep.append(i)

    triangles = triangles[keep] if keep else np.empty((0, 3), dtype=int)

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

    return {'vertices': all_pts, 'triangles': triangles}