# Grouper: the hull's sections (exec'd by build_grouper.py after gr_lib.py and gr_measure.py; also by gr_collision.py on
# its own). Pure functions of the model's y (nose -Y), no scene changes.
#
# Each station's section is a PROFILE: the port half from the keel (lat 0, bottom) to the crest (lat 0, top), as a dense
# polyline in (lat, z), built from the measured tables of gr_measure.py section 1:
#   * the hull: two superellipse quadrants joined at the widest line (HULL_TOP / HULL_BOT / HULL_HALF / HULL_YW, the
#     exponents of SECTION_EXPONENTS; my addition: the upper exponent rises from 2.1 to 2.8 on the nose's front face so the
#     face is flat across, as FRONT draws the front window and the side windows' inner edges);
#   * the hull's own top under the dorsal ridge is the drawn top line less the ridge's height (gr_measure: "its top IS the
#     hull's top line from 548 to 630"); the ridge is a separate part (gr_fittings.py);
#   * the thruster pod (POD): a circle at the drawn axis (lat 1.017, z -1.382) united with the hull by two rolling-ball
#     fillets (concave), so it sits in the hull's lower corner as FRONT and BACK draw it (the concave fairings from the flank
#     down to the pod). Its radius follows POD["radius"] from the tail (x 733) to the full 19.5 px (x 800-836), then an
#     elliptic nose to the front face at x 843 (a station pair, gr_hull.py). Toward the tail the circle's centre is drawn
#     in toward the hull so that it always overlaps it (the drawn tail would hang free in the water 0.4 m from the hull):
#     the lobe fades into the lower flank at x 733.
# The features the grid's columns lock onto are not marked here (they are heights and lateral offsets read off SIDE, TOP
# and FRONT, found on the profile by gr_hull.py).
import math
from mathutils import Vector


def smootherstep(a, b, x):
    """(round 3) a C2 step (Perlin's): the lobe's growth along the hull without a curvature jump at its ends"""
    t = max(0.0, min(1.0, (x - a) / (b - a)))
    return t * t * t * (t * (6 * t - 15) + 10)


def _pc(tab, cx=Yi, cv=Zi):
    """a measured table [(image x, value)] as a C2 cubic spline of model y (round 3: was a monotone cubic, C1)"""
    return cspline(sorted((cx(x), cv(v)) for x, v in tab))


def _lin(tab, cx=Yi, cv=lambda v: v):
    pts = sorted((cx(x), cv(v)) for x, v in tab)
    def f(y):
        if y <= pts[0][0]:
            return pts[0][1]
        if y >= pts[-1][0]:
            return pts[-1][1]
        for (a, va), (b, vb) in zip(pts, pts[1:]):
            if a <= y <= b:
                return va + (vb - va) * (y - a) / (b - a)
        return pts[-1][1]
    return f


def ximg(y):
    """model y -> SIDE/TOP image x"""
    return X0 - y / S


# ------------------------------------------------------------------ the dorsal ridge's height (px) taken off the top line
_RH = [(548.0, 9.0)] + list(DORSAL_RIDGE["height_px_at"])
ridge_h_px = _lin([(x, h) for x, h in _RH], cx=lambda x: x)            # of image x (clamped: 9 px aft of 551, 0 at 630+)


def _ridge_off(x):
    return ridge_h_px(x) if 548.0 <= x <= 630.0 else 0.0


HULL_TOP_OWN = []
for x, yv in HULL_TOP:
    if x < 543.5:
        HULL_TOP_OWN.append((x, yv))
    elif x < 548.0:                                                  # the stern shoulder eased into the lowered deck
        continue
    else:
        HULL_TOP_OWN.append((x, yv + _ridge_off(x)))
HULL_TOP_OWN += [(544.5, 336.0), (546.0, 333.6)]
HULL_TOP_OWN.sort()


def smooth_tab(tab, sig, step=0.5, out_step=0.5):
    """a measured table [(image x, value)] smoothed: resampled every `step` px, each sample the Gaussian-weighted mean over
    sig(x) px (the sheet's lines are read to 0.25-0.5 px, which shows as shading bands on a smooth hull), kept sharp where
    sig is small (the nose's face, the cheek); returned every out_step px.
    (round 3) the table is extended past its ends by an ODD reflection (v(x0 - d) = 2 v(x0) - v(x0 + d)), so the ends keep
    their value and slope (round 2's window was one-sided there and pulled the ends toward their neighbours); the samples are
    joined by a C2 natural spline (cspline), not the monotone cubic, whose curvature jumped at every 1 px sample"""
    xs = [a for a, b in tab]
    x0, x1 = xs[0], xs[-1]
    n = int(round((x1 - x0) / step))
    grid = [x0 + (x1 - x0) * k / n for k in range(n + 1)]

    def lin(x):
        for (a, va), (b, vb) in zip(tab, tab[1:]):
            if a <= x <= b:
                return va + (vb - va) * (x - a) / (b - a) if b > a else va
        return tab[0][1] if x < x0 else tab[-1][1]
    vals = [lin(x) for x in grid]

    def val(q):
        if q < 0:
            return 2.0 * vals[0] - vals[min(n, -q)]
        if q > n:
            return 2.0 * vals[n] - vals[max(0, 2 * n - q)]
        return vals[q]
    out = []
    m = max(1, int(round(out_step / step)))
    for k in range(0, n + 1, m):
        x = grid[k]
        sg = max(0.05, sig(x))
        w_sum = v_sum = 0.0
        r = int(math.ceil(3 * sg / step))
        for q in range(k - r, k + r + 1):
            w = math.exp(-0.5 * ((q - k) * step / sg) ** 2)
            w_sum += w; v_sum += w * val(q)
        out.append((x, v_sum / w_sum))
    if out[-1][0] < x1 - 1e-6:
        out.append((x1, vals[-1]))
    return out


def _blend(x, pts):
    """a smooth (C1) piecewise function through [(x, value)]: smoothstep between consecutive points, constant beyond"""
    if x <= pts[0][0]:
        return pts[0][1]
    for (xa, va), (xb, vb) in zip(pts, pts[1:]):
        if x <= xb:
            return va + (vb - va) * smoothstep(xa, xb, x)
    return pts[-1][1]


# (round 3) the smoothing widths vary smoothly along the hull (round 2's jumped at 548, 552, 872 and 883 px, which put kinks
# into the lines' curvature), with 2.0 px at the stern shoulder (was 0.3-0.4: the readings' kinks stayed in), 2.2 px over the
# nose's raked face and 2.8 px over the cheek and the step up to the nose's underside (x 860-878)
def _sig_top(x):
    return _blend(x, [(548.0, 2.0), (566.0, 7.0), (862.0, 7.0), (880.0, 2.2)])


def _sig_bot(x):
    return _blend(x, [(546.0, 2.0), (560.0, 5.0), (822.0, 5.0), (834.0, 2.8), (878.0, 2.8), (890.0, 2.0)])


def _sig_half(x):
    return _blend(x, [(548.0, 2.0), (570.0, 6.0)])


TOPf = _pc(smooth_tab(HULL_TOP_OWN, _sig_top))           # the hull's own crest z(y)
TOPDRAWNf = pchip(sorted((Yi(x), Zi(v)) for x, v in HULL_TOP))   # the drawn top line (with the ridge; raw: monotone, no overshoot)
BOTf = _pc(smooth_tab(HULL_BOT, _sig_bot))               # the keel z(y)
HALFf = _pc(smooth_tab(HULL_HALF, _sig_half), Yi, Lp)    # the half width (m) at the widest line
YWf = _pc(smooth_tab([(543.0, 379.1), (560.0, 378.0), (840.0, 378.0), (870.0, 372.0), (913.0, 368.0)], lambda x: 6.0))
NUf = _pc(smooth_tab([(543.0, 2.0), (562.0, 2.1), (862.0, 2.1), (890.0, 2.8), (913.0, 2.8)], lambda x: 5.0), Yi, lambda v: v)
NLf = _pc(smooth_tab([(x, b) for x, a, b in SECTION_EXPONENTS], lambda x: 5.0), Yi, lambda v: v)

Y_TAIL = Yi(543.0)                               # the hull's last station (the duct's front)
Y_NOSE = Yi(913.0)                               # the nose's front frame

# ------------------------------------------------------------------ the pod
# Critique round 1 (M2, M3): the pod is a STRAIGHT nacelle, its axis at (+-1.017, z -1.382) from the tail cone to the nose,
# with the measured radii (SIDE's bottom line). Its middle is part of the hull grid (the section's lobe, joined to the hull by
# the two rolling-ball fillets); its two ends are revolved pieces of their own (gr_fittings.py): the TAIL CONE (x 735-768,
# BACK's rings r 15/11/8 px and an end cap r 5 px, concentric on the axis, hanging free aft of x 747 where the pod leaves the
# hull) and the NOSE (x 831-846, SIDE's rounded bottom outline, the lamp's bezel). The grid's lobe is full on the axis from
# x 762 to 834 and sinks into the hull only INSIDE those pieces (x 745-762 and 834-846), where its faces are hidden and the
# finish pass prunes them. Before: the lobe's centre drifted in toward the hull aft of x 790 and closed elliptically at the
# nose (the bent tails BACK did not show, the folds at the pod's nose).
POD_L, POD_Z = Lp(POD["lateral_px"]), Zi(POD["y_c"])
X_POD_TAIL, X_POD_FULL0, X_POD_FULL1, X_POD_FACE = 745.0, 768.0, 834.0, 846.0   # lobe: fades in, full, full, faded out
# (round 3) the fade-in runs to the tail cone's front lip (x 768, was 762) and the fillets grow with it (below): the upper
# fillet swept up past the lower groove within 0.29 m and kinked the groove's line in 3D
X_TAILPIECE = (735.0, 768.0)                      # the tail cone piece (its end cap .. its front lip)
X_NOSEPIECE = (831.0, 846.0)                      # the nose piece (its aft lip .. the bezel)
_POD_R = [(733.0, 3.5), (736.0, 5.6), (740.0, 8.0), (745.0, 10.2)] + [(x, r) for x, r in POD["radius"] if 750.0 <= x <= 836.0]
_podr = pchip(_POD_R)
POD_NOSE_L = 10.0                                 # px: SIDE's outline rounds up from (836, 475.5) through (841, 468), (845, 459)
                                                  # to the axis height at 845.5: an elliptic nose 10 px long


def pod_r_true(x):
    """the drawn pod's radius (m) at image x (the nacelle the pieces and the lobe follow); 0 outside"""
    if x < 733.0 or x > 836.0 + POD_NOSE_L:
        return 0.0
    if x <= 836.0:
        return Lp(_podr(x))
    t = (x - 836.0) / POD_NOSE_L
    return Lp(19.5 * math.sqrt(max(0.0, 1.0 - t * t)))


def pod_r(x):
    """the grid lobe's radius (m) at image x: the drawn radius, 0 outside the lobe's span"""
    if x < X_POD_TAIL or x > X_POD_FACE + 1e-6:
        return 0.0
    return pod_r_true(min(x, 836.0))


# ------------------------------------------------------------------ the section
NSEC = 160                                        # points per quadrant of the hull curve
NL_KEEL_ADD = 0.6                                 # the lower exponent's growth from the widest line to the keel


def hull_params(y):
    zt, zb = TOPf(y), BOTf(y)
    zw = YWf(y)
    h = zt - zb
    zw = min(max(zw, zb + 0.30 * h), zt - 0.30 * h)
    return zt, zb, zw, HALFf(y), NUf(y), NLf(y)


def hull_curve(y, n=NSEC):
    """the hull's port half (no pod): points from the keel to the crest, uneven in angle (finer at the widest line)"""
    zt, zb, zw, W, nu, nl = hull_params(y)
    pts = []
    for k in range(n + 1):                          # lower quadrant, keel -> widest; its exponent grows toward the keel
        t = math.pi / 2 * (1 - k / n)               # (FRONT: the belly between the pods 10 cm lower than one exponent gives)
        c, s = math.cos(t), math.sin(t)
        e = nl + NL_KEEL_ADD * (t / (math.pi / 2)) ** 2
        pts.append((W * c ** (2 / e) if c > 0 else 0.0, zw - (zw - zb) * s ** (2 / e)))
    for k in range(1, n + 1):                       # upper quadrant, widest -> crest
        t = math.pi / 2 * k / n
        c, s = math.cos(t), math.sin(t)
        pts.append((W * c ** (2 / nu) if c > 1e-12 else 0.0, zw + (zt - zw) * s ** (2 / nu)))
    pts[0] = (0.0, zb)
    pts[-1] = (0.0, zt)
    return pts


def _normals(pts):
    """outward normals of a port-half profile (keel -> crest, counter-clockwise seen from aft... lat right, z up)"""
    out = []
    n = len(pts)
    for i in range(n):
        a = pts[max(0, i - 1)]; b = pts[min(n - 1, i + 1)]
        tx, tz = b[0] - a[0], b[1] - a[1]
        L = math.hypot(tx, tz) or 1.0
        out.append((tz / L, -tx / L))               # (lat, z): the tangent turned clockwise points outward
    out[0] = (0.0, -1.0)
    out[-1] = (0.0, 1.0)
    return out


def _seg_circle(p, q, c, R):
    """parameters t in [0, 1) where segment p->q crosses the circle (c, R)"""
    dx, dz = q[0] - p[0], q[1] - p[1]
    fx, fz = p[0] - c[0], p[1] - c[1]
    a = dx * dx + dz * dz
    b = 2 * (fx * dx + fz * dz)
    cc = fx * fx + fz * fz - R * R
    disc = b * b - 4 * a * cc
    if a < 1e-14 or disc < 0:
        return []
    s = math.sqrt(disc)
    return sorted(t for t in ((-b - s) / (2 * a), (-b + s) / (2 * a)) if 0.0 <= t < 1.0)


def _arc(c, r, a0, a1, step=0.01):
    """points on a circle from angle a0 to a1 (signed sweep), excluding the start"""
    n = max(2, int(math.ceil(abs(a1 - a0) * r / step)))
    return [(c[0] + r * math.cos(a0 + (a1 - a0) * k / n), c[1] + r * math.sin(a0 + (a1 - a0) * k / n)) for k in range(1, n + 1)]


def _nearest_on(pts, p):
    best = (1e9, 0, 0.0)
    for i in range(len(pts) - 1):
        a, b = pts[i], pts[i + 1]
        dx, dz = b[0] - a[0], b[1] - a[1]
        L2 = dx * dx + dz * dz or 1e-12
        t = max(0.0, min(1.0, ((p[0] - a[0]) * dx + (p[1] - a[1]) * dz) / L2))
        q = (a[0] + dx * t, a[1] + dz * t)
        d = math.hypot(p[0] - q[0], p[1] - q[1])
        if d < best[0]:
            best = (d, i, t)
    return best


def _inside_poly(pts, p):
    """inside the closed port half (the profile closed along the centre line)"""
    poly = list(pts) + [(0.0, pts[-1][1]), (0.0, pts[0][1])]
    x, y = p
    ins = False
    n = len(poly)
    for i in range(n):
        x1, y1 = poly[i]; x2, y2 = poly[(i + 1) % n]
        if (y1 > y) != (y2 > y) and x < (x2 - x1) * (y - y1) / (y2 - y1) + x1:
            ins = not ins
    return ins


POD_R_MIN = 0.12                                  # m: at its ends the lobe keeps this radius and sinks into the hull instead


def pod_lobe(y, H):
    """the pod's circle at y: (centre, radius) or None. Amidships (x 790-836) the drawn axis and radius. Toward the tail and at
    the nose the lobe's PROTRUSION beyond the hull (along the line from the hull's nearest point to the drawn axis) fades to
    zero (the tail eased in from x 722, the nose elliptic to x 847, as SIDE's outline rounds up), while its radius stays at
    least 0.12 m: the circle sinks into the hull rather than shrinking to a point (a vanishing lobe pinches the grid)"""
    x = ximg(y)
    r = pod_r(x)
    if x < X_POD_TAIL or x > X_POD_FACE:
        return None
    A = (POD_L, POD_Z)
    d, i, t = _nearest_on(H, A)
    a_, b_ = H[i], H[i + 1]
    P = (a_[0] + (b_[0] - a_[0]) * t, a_[1] + (b_[1] - a_[1]) * t)
    if _inside_poly(H, A):
        d = -d
    if abs(d) < 1e-6:
        nx, nz = 0.6, -0.8
    else:
        nx, nz = (A[0] - P[0]) / d, (A[1] - P[1]) / d
    if X_POD_FULL0 <= x <= X_POD_FULL1:
        return A, r                                             # the straight nacelle: its own axis and radius
    if x < X_POD_FULL0:                                         # the tail: sinks into the hull inside the tail cone piece
        e = max(0.0, d + r) * smootherstep(X_POD_TAIL, X_POD_FULL0, x)
    else:                                                        # the nose: sinks into the hull inside the nose piece
        e = max(0.0, d + r) * (1.0 - smootherstep(X_POD_FULL1, X_POD_FACE, x))
    re = max(r, POD_R_MIN)
    g = e - re
    return (P[0] + nx * g, P[1] + nz * g), re


def union_pod(H, c, r, R1, R2):
    """the outline of the hull curve H (keel -> crest) united with the circle (c, r), the two concave corners rounded by
    rolling balls: R1 below (between the hull's bottom and the pod), R2 above (the fairing from the flank down to the pod).
    Returns (points, (i0, i1), (s1, s2)): i0/i1 the outline indices where it leaves and rejoins H, s1/s2 those places as
    fractional indices of H."""
    N = _normals(H)

    def first_hit(R, last):
        HR = [(p[0] + n[0] * R, p[1] + n[1] * R) for p, n in zip(H, N)]
        hits = []
        for i in range(len(HR) - 1):
            for t in _seg_circle(HR[i], HR[i + 1], c, r + R):
                hits.append((i + t, i, t))
        if not hits:
            return None, HR
        return (hits[-1] if last else hits[0]), HR
    s1, HR1 = first_hit(R1, False)
    s2, HR2 = first_hit(R2, True)
    if s1 is None or s2 is None or s2[0] <= s1[0]:
        return list(H), None, None

    def at(seq, s):
        i, t = s[1], s[2]
        a, b = seq[i], seq[min(i + 1, len(seq) - 1)]
        return (a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t)
    F1, F2 = at(HR1, s1), at(HR2, s2)
    T1, T2 = at(H, s1), at(H, s2)
    C1 = (c[0] + r * (F1[0] - c[0]) / (r + R1), c[1] + r * (F1[1] - c[1]) / (r + R1))
    C2 = (c[0] + r * (F2[0] - c[0]) / (r + R2), c[1] + r * (F2[1] - c[1]) / (r + R2))

    def ang(p, o):
        return math.atan2(p[1] - o[1], p[0] - o[0])

    def short(a0, a1):
        d = (a1 - a0 + math.pi) % (2 * math.pi) - math.pi
        return a0, a0 + d
    out = list(H[:s1[1] + 1]) + [T1]
    i_leave = len(out) - 1
    a0, a1 = short(ang(T1, F1), ang(C1, F1))
    out += _arc(F1, R1, a0, a1, 0.008)
    # the pod's arc from C1 to C2 the outer way (not through the direction of the hull's inside)
    b0, b1 = ang(C1, c), ang(C2, c)
    d_ccw = (b1 - b0) % (2 * math.pi)
    # of the two arcs from C1 to C2 take the one whose middle lies outside the hull (the lobe's outer side)
    mid_ccw = (c[0] + r * math.cos(b0 + 0.5 * d_ccw), c[1] + r * math.sin(b0 + 0.5 * d_ccw))
    sweep = d_ccw if not _inside_poly(H, mid_ccw) else d_ccw - 2 * math.pi
    out += _arc(c, r, b0, b0 + sweep, 0.008)
    a0, a1 = short(ang(C2, F2), ang(T2, F2))
    out += _arc(F2, R2, a0, a1, 0.008)
    i_rejoin = len(out) - 1
    out += list(H[s2[1] + 1:])
    return out, (i_leave, i_rejoin), (s1[0], s2[0])


class Prof:
    """a section profile: points (lat, z) keel -> crest; with a pod: its span on the outline (indices where it leaves and
    rejoins the hull curve), the same places on the hull curve (fractional indices of H), the circle"""
    def __init__(self, pts, pod=None, c=None, r=0.0, H=None, hs=None):
        self.p = pts
        self.pod = pod
        self.c, self.r = c, r
        self.H, self.hs = H, hs


def profile(y, pod=True):
    """the section at y (the pod united unless pod=False)"""
    H = hull_curve(y)
    lobe = pod_lobe(y, H) if pod else None
    if lobe is None:
        return Prof(H, H=H)
    c, r = lobe
    # the fillets scale with the lobe (no floor: a ball larger than a vanishing lobe makes the outline loop at the pod's ends)
    R1 = min(0.12, 0.36 * r)                                  # below: the notch FRONT and BACK draw at the belly, softened
    R2 = min(0.20, 0.60 * r)                                  # above: the concave fairing from the flank down to the pod
    x_ = ximg(y)
    g_ = smootherstep(X_POD_TAIL, X_POD_FULL0, x_) * (1.0 - smootherstep(X_POD_FULL1, X_POD_FACE, x_))
    R1, R2 = max(1e-4, R1 * g_), max(1e-4, R2 * g_)            # (round 3) the fillets grow with the lobe, from nothing
    d = _nearest_on(H, c)[0]
    inside = _inside_poly(H, c)
    if inside and d >= r - 1e-4:                              # sunk wholly into the hull (the nose and tail fades): no lobe
        return Prof(H, H=H)
    if not inside:
        R1 = max(R1, (0.6 * (d - r) + 0.1 * r) * max(g_, 0.15))
        R2 = max(R2, (0.6 * (d - r) + 0.1 * r) * max(g_, 0.15))
    pts, span, hs = union_pod(H, c, r, R1, R2)
    return Prof(pts, span, c, r, H, hs)


def arc_of(P):
    arc = [0.0]
    for i in range(1, len(P.p)):
        arc.append(arc[-1] + math.hypot(P.p[i][0] - P.p[i - 1][0], P.p[i][1] - P.p[i - 1][1]))
    return arc


def at_arc(P, arc, s):
    if s <= 0:
        return P.p[0]
    if s >= arc[-1]:
        return P.p[-1]
    lo, hi = 0, len(arc) - 1
    while hi - lo > 1:
        m = (lo + hi) // 2
        if arc[m] <= s:
            lo = m
        else:
            hi = m
    t = (s - arc[lo]) / max(1e-12, arc[hi] - arc[lo])
    (a0, b0), (a1, b1) = P.p[lo], P.p[hi]
    return (a0 + (a1 - a0) * t, b0 + (b1 - b0) * t)


def arc_at_z_outer(P, arc, z):
    """arc position of the outermost (largest lat) point at height z: what SIDE sees at that height"""
    best = None
    for i in range(1, len(P.p)):
        (l0, z0), (l1, z1) = P.p[i - 1], P.p[i]
        if (z0 - z) * (z1 - z) <= 0 and z0 != z1:
            t = (z - z0) / (z1 - z0)
            lat = l0 + (l1 - l0) * t
            if best is None or lat > best[0]:
                best = (lat, arc[i - 1] + (arc[i] - arc[i - 1]) * t)
    if best is None:
        return arc[-1] if z > max(p[1] for p in P.p) else 0.0
    return best[1]


def arc_at_lat_top(P, arc, lat):
    """arc position on the upper side (counted from the crest backward) where the profile reaches lateral offset lat"""
    for i in range(len(P.p) - 1, 0, -1):
        l0, l1 = P.p[i][0], P.p[i - 1][0]
        if (l0 - lat) * (l1 - lat) <= 0 and l0 != l1:
            return arc[i] + (arc[i - 1] - arc[i]) * (lat - l0) / (l1 - l0)
    return arc[-1]


def widest_arc(P, arc):
    """arc position of the widest point on the upper half (above the pod)"""
    i = max(range(len(P.p)), key=lambda k: (P.p[k][0] if P.p[k][1] > POD_Z + 0.35 else -1.0))
    return arc[i]


def section_pts(y, n=16):
    """n points round the whole section at y with outward 2D normals (for collision hulls when no grid is built)"""
    P = profile(y)
    arc = arc_of(P)
    L = arc[-1]
    out = []
    half = n // 2
    for k in range(half + 1):
        p = at_arc(P, arc, L * k / half)
        out.append(Vector((p[0], y, p[1])))
    return out + [Vector((-p.x, p.y, p.z)) for p in reversed(out[1:-1])]


def side_poly(img):
    return [(Yi(x), Zi(v)) for x, v in img]


def pip(pt, poly):
    x, y = pt
    inside = False
    n = len(poly)
    for i in range(n):
        x1, y1 = poly[i]; x2, y2 = poly[(i + 1) % n]
        if (y1 > y) != (y2 > y) and x < (x2 - x1) * (y - y1) / (y2 - y1) + x1:
            inside = not inside
    return inside


def span_at(poly, y):
    """the (min, max) of a closed 2D polygon's crossings at abscissa y (None if it misses)"""
    zs = []
    n = len(poly)
    for i in range(n):
        (y1, z1), (y2, z2) = poly[i], poly[(i + 1) % n]
        if (y1 - y) * (y2 - y) <= 0 and y1 != y2:
            zs.append(z1 + (z2 - z1) * (y - y1) / (y2 - y1))
    return (min(zs), max(zs)) if len(zs) >= 2 else None


def smooth_closed(pts, r_px, step=1.0):
    """round the corners of a closed polyline (image px) by moving along it: a polyline with each corner replaced by a
    quadratic arc reaching r_px along both edges"""
    n = len(pts)
    out = []
    for i in range(n):
        a, b, c = Vector(pts[i - 1]), Vector(pts[i]), Vector(pts[(i + 1) % n])
        la, lc = (a - b).length, (c - b).length
        ra, rc = min(r_px, 0.45 * la), min(r_px, 0.45 * lc)
        p0 = b + (a - b).normalized() * ra
        p2 = b + (c - b).normalized() * rc
        m = 6
        for k in range(m + 1):
            t = k / m
            q = (1 - t) ** 2 * p0 + 2 * (1 - t) * t * b + t * t * p2
            out.append((q.x, q.y))
    return out


def resample_closed(pts, step):
    """a closed polyline resampled every `step` (same units)"""
    P_ = [Vector(p) for p in pts] + [Vector(pts[0])]
    L = [0.0]
    for i in range(1, len(P_)):
        L.append(L[-1] + (P_[i] - P_[i - 1]).length)
    n = max(8, int(round(L[-1] / step)))
    out = []
    j = 0
    for k in range(n):
        s = L[-1] * k / n
        while j < len(L) - 2 and L[j + 1] < s:
            j += 1
        t = (s - L[j]) / max(1e-12, L[j + 1] - L[j])
        q = P_[j].lerp(P_[j + 1], t)
        out.append((q.x, q.y))
    return out