"""
Constrained Delaunay triangulation with density-adaptive refinement.

Uses Shewchuk's Triangle library to generate adaptive mesh that
respects plate boundary, hole keepouts, and density-driven spacing.
"""

import numpy as np
import triangle as tr
from shapely.geometry import Polygon

from .density_field import evaluate_density, density_to_spacing


def offset_polygon(coords, distance, inward=True):
    """
    Offset a polygon boundary by `distance`.
    inward=True shrinks, inward=False expands.
    
    Uses Shapely buffer (negative = inward for exterior ring).
    """
    poly = Polygon(coords)
    if inward:
        buffered = poly.buffer(-distance)
    else:
        buffered = poly.buffer(distance)
    
    if buffered.is_empty:
        return coords  # can't offset that much, return original
    
    if hasattr(buffered, 'exterior'):
        return list(buffered.exterior.coords)[:-1]  # remove closing duplicate
    return coords


def sample_circle(center, radius, num_points=32):
    """Sample a circle as a polygon with `num_points` vertices."""
    cx, cy = center
    angles = np.linspace(0.0, 2.0 * np.pi, num_points, endpoint=False)
    return [[cx + radius * np.cos(a), cy + radius * np.sin(a)] for a in angles]


def build_pslg(geometry, params):
    """
    Build Planar Straight Line Graph for Triangle library.
    
    Parameters
    ----------
    geometry : dict
        Plate geometry (outer_boundary, holes).
    params : dict
        Must contain d_keep (hole keepout multiplier).
    
    Returns
    -------
    dict : Triangle-compatible PSLG with vertices, segments, holes.
    """
    d_keep = params['d_keep']
    
    vertices = []
    segments = []
    hole_markers = []
    
    # Outer boundary (no offset needed — frame is handled in profile assembly)
    outer = geometry['outer_boundary']
    v_start = len(vertices)
    vertices.extend(outer)
    n = len(outer)
    for i in range(n):
        segments.append([v_start + i, v_start + (i + 1) % n])
    
    # Each hole with boss keepout reservation
    for hole in geometry['holes']:
        diameter = float(hole.get('diameter', 10.0) or 10.0)
        keepout_dist = d_keep * diameter / 2.0

        if hole.get('is_circular', False) and 'center' in hole:
            # Circular boss reservation around hole:
            # r_boss = r_hole + d_keep * hole_diameter / 2
            hole_radius = diameter / 2.0
            boss_radius = hole_radius + keepout_dist
            keepout_boundary = sample_circle(hole['center'], boss_radius, num_points=32)
        else:
            # Fallback for non-circular holes
            keepout_boundary = offset_polygon(hole['boundary'], keepout_dist, inward=False)

        v_start = len(vertices)
        vertices.extend(keepout_boundary)
        n_h = len(keepout_boundary)
        for i in range(n_h):
            segments.append([v_start + i, v_start + (i + 1) % n_h])

        # Marker inside hole tells Triangle to leave this keepout region empty
        hole_markers.append(hole['center'])
    
    result = {
        'vertices': np.array(vertices, dtype=np.float64),
        'segments': np.array(segments, dtype=np.int32),
    }
    if hole_markers:
        result['holes'] = np.array(hole_markers, dtype=np.float64)
    
    return result


def compute_triangle_areas(vertices, triangles):
    """Compute area of each triangle."""
    v0 = vertices[triangles[:, 0]]
    v1 = vertices[triangles[:, 1]]
    v2 = vertices[triangles[:, 2]]
    # Cross product / 2
    areas = 0.5 * np.abs(
        (v1[:, 0] - v0[:, 0]) * (v2[:, 1] - v0[:, 1]) -
        (v2[:, 0] - v0[:, 0]) * (v1[:, 1] - v0[:, 1])
    )
    return areas


def compute_centroids(vertices, triangles):
    """Compute centroid of each triangle."""
    v0 = vertices[triangles[:, 0]]
    v1 = vertices[triangles[:, 1]]
    v2 = vertices[triangles[:, 2]]
    return (v0 + v1 + v2) / 3.0


def filter_small_triangles(result, min_triangle_area):
    """Remove triangles smaller than the manufacturing threshold."""
    triangles = result.get('triangles')
    vertices = result.get('vertices')
    if triangles is None or vertices is None or len(triangles) == 0:
        return result

    areas = compute_triangle_areas(vertices, triangles)
    keep_mask = areas >= float(min_triangle_area)

    result['triangle_areas'] = areas
    result['small_triangle_mask'] = ~keep_mask
    result['triangles'] = triangles[keep_mask]
    return result


def generate_triangulation(geometry, params, max_refinement_passes=3):
    """
    Generate density-adaptive constrained Delaunay triangulation.
    
    Parameters
    ----------
    geometry : dict
        Plate geometry.
    params : dict
        Full parameter set (density field + spacing + manufacturing).
    max_refinement_passes : int
        Number of iterative refinement passes.
    
    Returns
    -------
    dict : Triangle result with 'vertices' and 'triangles'.
    """
    pslg = build_pslg(geometry, params)
    
    # Initial triangulation with global max area
    s_max = params['s_max']
    global_max_area = (np.sqrt(3) / 4.0) * s_max**2
    
    # Triangle options: p=PSLG, q30=min angle 30°, a=area constraint, D=conforming Delaunay
    result = tr.triangulate(pslg, f'pq30Da{global_max_area:.1f}')
    
    # Iterative refinement based on density field
    for iteration in range(max_refinement_passes):
        verts = result['vertices']
        tris = result['triangles']

        areas = compute_triangle_areas(verts, tris)
        centroids = compute_centroids(verts, tris)

        # Compute target area for each triangle based on density at centroid
        target_areas = np.array([
            (np.sqrt(3) / 4.0) * density_to_spacing(
                evaluate_density(cx, cy, geometry, params), params
            )**2
            for cx, cy in centroids
        ])

        # Check if all triangles satisfy constraints (20% tolerance)
        if np.all(areas <= target_areas * 1.2):
            break

        # Set per-triangle max area and refine
        result['triangle_max_area'] = target_areas
        result = tr.triangulate(result, 'rpq30D')

    min_triangle_area = params.get('min_triangle_area', 20.0)
    result = filter_small_triangles(result, min_triangle_area)
    return result