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
from shapely.ops import unary_union

from .density_field import evaluate_density, density_to_spacing


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 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:
        # Handle MultiPolygon from buffer on complex shapes
        if inner_frame.geom_type == 'MultiPolygon':
            inner_polys = list(inner_frame.geoms)
        else:
            inner_polys = [inner_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]:
                    new_pts.append([cx, cy])

    # 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
            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])

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


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
    all_pts = _add_boundary_vertices(grid_pts, geometry, params, 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'])
    w_frame = params.get('w_frame', 8.0)
    inner_plate = plate_poly.buffer(-w_frame)
    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}