Atomizer/tools/adaptive-isogrid/archive/deprecated-triangle-mesher/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.

    Uses Shapely simplify() to reduce vertex count on curved buffer segments,
    then adds vertices from the simplified ring plus interpolated points on
    long edges. This preserves corners/notches while avoiding vertex clusters.
    """
    # Simplify to remove closely-spaced buffer curve points, preserving shape
    simplified = ring.simplify(spacing * 0.15, preserve_topology=True)
    coords = np.array(simplified.coords)
    if len(coords) > 1 and np.allclose(coords[0], coords[-1]):
        coords = coords[:-1]
    if len(coords) < 3:
        # Fallback to uniform interpolation
        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)

    pts = []
    n = len(coords)
    for i in range(n):
        p1 = coords[i]
        p2 = coords[(i + 1) % n]
        pts.append(p1.tolist())

        # Add interpolated points on long edges
        edge_len = np.linalg.norm(p2 - p1)
        if edge_len > spacing * 1.5:
            n_sub = int(np.ceil(edge_len / spacing))
            for j in range(1, n_sub):
                t = j / n_sub
                pts.append((p1 + t * (p2 - p1)).tolist())

    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.

    Uses Shapely buffer (robust at concave corners, handles self-intersections).
    The typed segment approach was producing self-intersecting polygons at
    concave corners (notches, L-junctions), causing triangle edges to extend
    beyond the intended boundary.
    """
    w_frame = float(params.get('w_frame', 8.0))
    plate_poly = Polygon(geometry['outer_boundary'])

    inner_plate = plate_poly.buffer(-w_frame, resolution=16)
    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
    # _sample_ring now uses actual polygon vertices (preserving tight features)
    # and only adds interpolated points on long straight edges.
    boundary_spacing = max(s_min, min(s_max * 0.4, 25.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]

    # Step 7: Snap out-of-bounds vertices to nearest boundary point
    # Only snap vertices that are clearly outside (> 0.1mm), not boundary vertices
    snap_tol = 0.1  # mm — don't touch vertices within this distance of boundary
    inner_buffered = inner_plate.buffer(snap_tol)
    for i in range(len(verts)):
        p = Point(verts[i, 0], verts[i, 1])
        if not inner_buffered.contains(p):
            nearest = inner_plate.exterior.interpolate(
                inner_plate.exterior.project(p)
            )
            verts[i, 0] = nearest.x
            verts[i, 1] = nearest.y

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