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, 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


def _boundary_layer_offset_for_segment(mid_pt, geometry, params):
    """Choose inward offset for boundary seed row."""
    explicit = params.get('boundary_layer_offset', None)
    if explicit is not None:
        return max(float(explicit), 0.0)
    eta = evaluate_density(mid_pt[0], mid_pt[1], geometry, params)
    return max(float(density_to_spacing(eta, params)), 1e-3)


def _add_boundary_layer_seed_points(points, geometry, params, plate_poly, keepout_union):
    """Add a structured point row offset inward from each straight outer edge."""
    boundary_pts = []
    ring = LinearRing(geometry['outer_boundary'])
    is_ccw = bool(ring.is_ccw)

    # Prefer typed segments to avoid treating discretized arcs as straight edges
    typed = geometry.get('outer_boundary_typed')
    if typed:
        segments = [seg for seg in typed if seg.get('type', 'line') == 'line']
        edge_pairs = [(_np(seg['start']), _np(seg['end'])) for seg in segments]
    else:
        coords = np.asarray(geometry['outer_boundary'], dtype=float)
        if len(coords) >= 2 and np.allclose(coords[0], coords[-1]):
            coords = coords[:-1]
        edge_pairs = []
        for i in range(len(coords)):
            edge_pairs.append((coords[i], coords[(i + 1) % len(coords)]))

    for a, b in edge_pairs:
        dx, dy = b[0] - a[0], b[1] - a[1]
        edge_len = float(np.hypot(dx, dy))
        if edge_len < 1e-9:
            continue

        mid = np.array([(a[0] + b[0]) * 0.5, (a[1] + b[1]) * 0.5], dtype=float)
        spacing = float(density_to_spacing(evaluate_density(mid[0], mid[1], geometry, params), params))
        spacing = max(spacing, 1e-3)
        offset = _boundary_layer_offset_for_segment(mid, geometry, params)

        nx_l, ny_l = (-dy / edge_len), (dx / edge_len)
        nx, ny = (nx_l, ny_l) if is_ccw else (-nx_l, -ny_l)

        n_pts = max(int(np.floor(edge_len / spacing)), 1)
        for k in range(1, n_pts + 1):
            t = k / (n_pts + 1)
            bx = a[0] + t * dx
            by = a[1] + t * dy
            px = bx + offset * nx
            py = by + offset * ny
            p = Point(px, py)
            if not plate_poly.buffer(1e-6).contains(p):
                continue
            if not keepout_union.is_empty and keepout_union.contains(p):
                continue
            boundary_pts.append([px, py])

    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 _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:
            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
            for cx, cy in list(ring.coords)[:-1]:
                new_pts.append([cx, cy])
            length = ring.length
            n_pts = max(int(length / (s_min * 0.7)), 6)
            for i in range(n_pts):
                frac = i / n_pts
                pt = ring.interpolate(frac, normalized=True)
                new_pts.append([pt.x, pt.y])

    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, plate_poly, 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

    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 all_inside:
            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}