# Polygon helpers shared by the model builders:  boundary tracing of a bitmap mask, simplification,
# offsetting, orientation and triangulation.  Rings are Nx2 numpy arrays of (x, y) points.
import numpy as np


def traceBoundary(mask):
    # Moore-neighbour tracing, clockwise on screen, starting from the top-left foreground pixel.
    padded = np.pad(mask, 1)
    ys, xs = np.nonzero(padded)
    start = (ys[0], xs[0])
    steps = [(0, 1), (1, 1), (1, 0), (1, -1), (0, -1), (-1, -1), (-1, 0), (-1, 1)]
    boundary = [start]
    current = start
    came = 6
    while True:
        found = False
        for k in range(8):
            d = (came + 1 + k) % 8
            ny = current[0] + steps[d][0]
            nx = current[1] + steps[d][1]
            if padded[ny, nx]:
                came = (d + 4) % 8
                # Next search starts just after the backtrack direction.
                came = (d + 5) % 8
                current = (ny, nx)
                found = True
                break
        if not found or (current == start and len(boundary) > 1):
            break
        boundary.append(current)
        if len(boundary) > 100000:
            break
    return np.array([(x - 1, y - 1) for y, x in boundary], dtype=float)


def simplify(points, tolerance):
    # Douglas-Peucker on a closed ring, split at the two points farthest apart.
    def dp(pts):
        if len(pts) < 3:
            return pts
        a = pts[0]
        b = pts[-1]
        d = b - a
        n = np.hypot(*d)
        if n < 1e-9:
            dist = np.hypot(*(pts - a).T)
        else:
            dist = np.abs(d[0] * (pts[:, 1] - a[1]) - d[1] * (pts[:, 0] - a[0])) / n
        i = int(np.argmax(dist))
        if dist[i] > tolerance:
            return np.vstack([dp(pts[:i + 1])[:-1], dp(pts[i:])])
        return np.array([a, b])
    far = int(np.argmax(np.hypot(*(points - points[0]).T)))
    first = dp(points[:far + 1])
    second = dp(np.vstack([points[far:], points[:1]]))
    ring = np.vstack([first[:-1], second[:-1]])
    return ring


def signedArea(ring):
    n = len(ring)
    return 0.5 * sum(ring[i][0] * ring[(i + 1) % n][1] - ring[(i + 1) % n][0] * ring[i][1] for i in range(n))


def offsetRing(ring, amount):
    # Grow (positive) or shrink a ring by a distance; mitres are clamped so acute tips do not spike.
    n = len(ring)
    sign = 1.0 if signedArea(ring) > 0 else -1.0
    result = []
    for i in range(n):
        p = ring[i]
        d0 = p - ring[i - 1]
        d1 = ring[(i + 1) % n] - p
        d0 = d0 / max(1e-9, np.hypot(*d0))
        d1 = d1 / max(1e-9, np.hypot(*d1))
        n0 = np.array([d0[1], -d0[0]]) * sign
        n1 = np.array([d1[1], -d1[0]]) * sign
        bis = n0 + n1
        length = np.hypot(*bis)
        if length < 1e-6:
            result.append(p + n0 * amount)
            continue
        bis = bis / length
        mitre = amount / max(0.35, bis @ n0)
        result.append(p + bis * mitre)
    return np.array(result)


def earClip(ring):
    # Triangulates a simple polygon; returns index triples into ring.
    n = len(ring)
    indices = list(range(n))
    if signedArea(ring) < 0:
        indices.reverse()
    triangles = []

    def convex(a, b, c):
        return (ring[b][0] - ring[a][0]) * (ring[c][1] - ring[a][1]) - (ring[b][1] - ring[a][1]) * (ring[c][0] - ring[a][0]) > 1e-9

    def inside(p, a, b, c):
        def side(u, v):
            return (v[0] - u[0]) * (p[1] - u[1]) - (v[1] - u[1]) * (p[0] - u[0])
        return side(ring[a], ring[b]) >= 0 and side(ring[b], ring[c]) >= 0 and side(ring[c], ring[a]) >= 0

    guard = 0
    while len(indices) > 3 and guard < 10 * n:
        guard += 1
        clipped = False
        for k in range(len(indices)):
            a = indices[k - 1]
            b = indices[k]
            c = indices[(k + 1) % len(indices)]
            if not convex(a, b, c):
                continue
            if any(inside(ring[m], a, b, c) for m in indices if m not in (a, b, c)):
                continue
            triangles.append((a, b, c))
            del indices[k]
            clipped = True
            break
        if not clipped:
            break
    if len(indices) == 3:
        triangles.append(tuple(indices))
    return triangles


def isConvex(ring):
    n = len(ring)
    sign = 0.0
    for i in range(n):
        a = ring[i]
        b = ring[(i + 1) % n]
        c = ring[(i + 2) % n]
        cross = (b[0] - a[0]) * (c[1] - b[1]) - (b[1] - a[1]) * (c[0] - b[0])
        if abs(cross) < 1e-9:
            continue
        if sign == 0.0:
            sign = np.sign(cross)
        elif np.sign(cross) != sign:
            return False
    return sign != 0.0


def selfIntersects(ring):
    # True when any two non-adjacent edges of the ring cross.
    n = len(ring)
    for i in range(n):
        a0 = ring[i]
        a1 = ring[(i + 1) % n]
        for j in range(i + 2, n):
            if i == 0 and j == n - 1:
                continue
            b0 = ring[j]
            b1 = ring[(j + 1) % n]
            d1 = (a1[0] - a0[0]) * (b0[1] - a0[1]) - (a1[1] - a0[1]) * (b0[0] - a0[0])
            d2 = (a1[0] - a0[0]) * (b1[1] - a0[1]) - (a1[1] - a0[1]) * (b1[0] - a0[0])
            d3 = (b1[0] - b0[0]) * (a0[1] - b0[1]) - (b1[1] - b0[1]) * (a0[0] - b0[0])
            d4 = (b1[0] - b0[0]) * (a1[1] - b0[1]) - (b1[1] - b0[1]) * (a1[0] - b0[0])
            if ((d1 > 0) != (d2 > 0)) and ((d3 > 0) != (d4 > 0)):
                return True
    return False


def removeSpikes(ring, degrees=8.0):
    # Drops vertices where the outline turns back on itself, the hair-thin fingers a raster boundary
    # leaves at one-pixel spurs.
    limit = -np.cos(np.radians(degrees))
    changed = True
    while changed and len(ring) > 3:
        changed = False
        n = len(ring)
        for i in range(n):
            d0 = ring[i] - ring[i - 1]
            d1 = ring[(i + 1) % n] - ring[i]
            l0 = np.hypot(*d0)
            l1 = np.hypot(*d1)
            if l0 < 1e-9 or l1 < 1e-9 or (d0 @ d1) / (l0 * l1) < limit:
                ring = np.delete(ring, i, axis=0)
                changed = True
                break
    return ring


def simplifyClean(points, tolerance):
    # Simplify, drop spikes, and back off the tolerance until the ring does not cross itself.
    while True:
        ring = removeSpikes(simplify(points, tolerance))
        if not selfIntersects(ring) or tolerance < 0.2:
            return ring
        tolerance = tolerance / 2