"""Planar geometry helpers for the OSM → asset pipeline.

Pure Python: no `bpy`, so this runs and tests outside Blender.  All functions
work in projected metres (see `osmassets.osm.Projector`) unless the name says
otherwise; `geometry_rings` and `feature_in_bounds` take raw GeoJSON and are the
two exceptions, operating on lon/lat.

Rings are lists of (x, y) tuples.  A repeated closing point is tolerated
everywhere but never required.
"""

import math


def geometry_rings(geometry):
    """Exterior rings of a GeoJSON Polygon/MultiPolygon; holes are dropped."""
    if not geometry:
        return []
    kind = geometry.get("type")
    coordinates = geometry.get("coordinates", [])
    if kind == "Polygon":
        return coordinates[:1]
    if kind == "MultiPolygon":
        return [polygon[0] for polygon in coordinates if polygon]
    return []


def feature_in_bounds(feature, projector):
    """True when any coordinate of the feature falls inside the padded bounds."""
    def walk(value):
        if isinstance(value, list) and value and isinstance(value[0], (int, float)):
            return projector.inside(value)
        return any(walk(v) for v in value) if isinstance(value, list) else False
    return walk(feature.get("geometry", {}).get("coordinates", []))


def clip_polygon(ring, xmin, xmax, ymin, ymax):
    """Sutherland-Hodgman clip of a ring against an axis-aligned box."""
    if len(ring) < 3:
        return []

    def clip_edge(points, inside, intersection):
        if not points:
            return []
        result = []
        previous = points[-1]
        previous_inside = inside(previous)
        for current in points:
            current_inside = inside(current)
            if current_inside != previous_inside:
                result.append(intersection(previous, current))
            if current_inside:
                result.append(current)
            previous = current
            previous_inside = current_inside
        return result

    ring = clip_edge(
        ring, lambda p: p[0] >= xmin,
        lambda a, b: (xmin, a[1] + (b[1] - a[1]) * (xmin - a[0]) /
                      (b[0] - a[0]) if b[0] != a[0] else a[1]))
    ring = clip_edge(
        ring, lambda p: p[0] <= xmax,
        lambda a, b: (xmax, a[1] + (b[1] - a[1]) * (xmax - a[0]) /
                      (b[0] - a[0]) if b[0] != a[0] else a[1]))
    ring = clip_edge(
        ring, lambda p: p[1] >= ymin,
        lambda a, b: (a[0] + (b[0] - a[0]) * (ymin - a[1]) /
                      (b[1] - a[1]) if b[1] != a[1] else a[0], ymin))
    ring = clip_edge(
        ring, lambda p: p[1] <= ymax,
        lambda a, b: (a[0] + (b[0] - a[0]) * (ymax - a[1]) /
                      (b[1] - a[1]) if b[1] != a[1] else a[0], ymax))
    return ring


def sample_tree_row(points, spacing, height):
    """Evenly space (x, y, height) samples along a polyline.

    The trailing point is appended only when the last regular sample stops well
    short of it, so a row does not end in a double-planted tree.
    """
    if len(points) < 2:
        return []
    samples = [(points[0][0], points[0][1], height)]
    distance_until_next = spacing
    for start, end in zip(points, points[1:]):
        dx = end[0] - start[0]
        dy = end[1] - start[1]
        segment_length = math.hypot(dx, dy)
        if segment_length == 0:
            continue
        while distance_until_next <= segment_length:
            ratio = distance_until_next / segment_length
            samples.append((start[0] + dx * ratio, start[1] + dy * ratio, height))
            distance_until_next += spacing
        distance_until_next -= segment_length
    last = points[-1]
    if math.hypot(samples[-1][0] - last[0], samples[-1][1] - last[1]) > spacing * 0.45:
        samples.append((last[0], last[1], height))
    return samples


def polygon_area(ring):
    """Unsigned shoelace area; 0.0 for degenerate rings."""
    return abs(signed_polygon_area(ring))


def signed_polygon_area(ring):
    """Signed shoelace area; positive for counter-clockwise rings."""
    if len(ring) < 3:
        return 0.0
    area = 0.0
    for (x1, y1), (x2, y2) in zip(ring, ring[1:] + ring[:1]):
        area += x1 * y2 - x2 * y1
    return area * 0.5


def sample_ring_boundary(ring, spacing, inset=0.0, max_samples=None):
    """Evenly sample a closed ring's boundary.

    Returns (x, y, angle, index) samples.  ``angle`` follows the local edge
    direction, and ``inset`` moves the sample toward the polygon interior.
    """
    if len(ring) > 1 and ring[0] == ring[-1]:
        ring = ring[:-1]
    if len(ring) < 3 or spacing <= 0.0:
        return []

    edges = []
    perimeter = 0.0
    winding = signed_polygon_area(ring)
    for index, (start, end) in enumerate(zip(ring, ring[1:] + ring[:1])):
        dx = end[0] - start[0]
        dy = end[1] - start[1]
        length = math.hypot(dx, dy)
        if length <= 1e-9:
            continue
        ux = dx / length
        uy = dy / length
        # Counter-clockwise rings have their interior on the left side of each
        # edge; clockwise rings have it on the right.
        inward = (-uy, ux) if winding >= 0.0 else (uy, -ux)
        edges.append((perimeter, start, ux, uy, length, inward, index))
        perimeter += length
    if not edges:
        return []

    count = max(1, int(perimeter / spacing))
    if max_samples:
        count = min(count, max_samples)
    step = perimeter / count
    samples = []
    edge_cursor = 0
    for sample_index in range(count):
        target = (sample_index + 0.5) * step
        while edge_cursor + 1 < len(edges) and (
            edges[edge_cursor][0] + edges[edge_cursor][4] < target
        ):
            edge_cursor += 1
        edge_start, start, ux, uy, length, inward, _ = edges[edge_cursor]
        along = max(0.0, min(length, target - edge_start))
        x = start[0] + ux * along
        y = start[1] + uy * along
        sx = x + inward[0] * inset
        sy = y + inward[1] * inset
        if inset > 0.0 and not point_in_polygon((sx, sy), ring):
            sx, sy = x, y
        samples.append((sx, sy, math.atan2(uy, ux), sample_index))
    return samples


def sample_polygon_interior(ring, spacing, edge_clearance=0.0, max_samples=None,
                            seed=0):
    """Jittered interior samples for sparse planting inside a polygon."""
    if len(ring) > 1 and ring[0] == ring[-1]:
        ring = ring[:-1]
    if len(ring) < 3 or spacing <= 0.0 or polygon_area(ring) <= 1e-9:
        return []

    xmin = min(x for x, _ in ring)
    xmax = max(x for x, _ in ring)
    ymin = min(y for _, y in ring)
    ymax = max(y for _, y in ring)
    cols = max(1, int(math.ceil((xmax - xmin) / spacing)))
    rows = max(1, int(math.ceil((ymax - ymin) / spacing)))
    samples = []
    for col in range(cols):
        for row in range(rows):
            sample_seed = ((col + 1) * 73856093) ^ ((row + 1) * 19349663) ^ seed
            jx = ((sample_seed * 0.61803398875) % 1.0 - 0.5) * spacing * 0.7
            jy = ((sample_seed * 0.41421356237) % 1.0 - 0.5) * spacing * 0.7
            x = xmin + (col + 0.5) * spacing + jx
            y = ymin + (row + 0.5) * spacing + jy
            if not point_in_polygon((x, y), ring):
                continue
            if edge_clearance > 0.0 and distance_to_ring((x, y), ring) < edge_clearance:
                continue
            samples.append((x, y, sample_seed))

    if max_samples and len(samples) > max_samples:
        samples.sort(key=lambda item: (item[2] * 0.754877666) % 1.0)
        samples = samples[:max_samples]
    return samples


def point_in_polygon(point, ring):
    x, y = point
    inside = False
    j = len(ring) - 1
    for i, (xi, yi) in enumerate(ring):
        xj, yj = ring[j]
        crosses = ((yi > y) != (yj > y))
        if crosses:
            x_at_y = (xj - xi) * (y - yi) / (yj - yi) + xi
            if x < x_at_y:
                inside = not inside
        j = i
    return inside


def distance_to_ring(point, ring):
    """Shortest distance from a point to the ring's edges (not its interior)."""
    px, py = point
    best = float("inf")
    count = len(ring)
    for index in range(count):
        ax, ay = ring[index]
        bx, by = ring[(index + 1) % count]
        dx = bx - ax
        dy = by - ay
        length_sq = dx * dx + dy * dy
        if length_sq <= 1e-9:
            distance = math.hypot(px - ax, py - ay)
        else:
            t = ((px - ax) * dx + (py - ay) * dy) / length_sq
            t = max(0.0, min(1.0, t))
            distance = math.hypot(px - (ax + t * dx), py - (ay + t * dy))
        if distance < best:
            best = distance
    return best