"""Grouper, round 3 (2026-10-22): FAIR design curves from the sheet's own drawn strokes (plain Python: numpy + PIL).

The user sent the hull back: "fix the topology and lines needs to be smooth". The round-2 lines were the Designer's coarse
tables (a reading every 5-10 px, read to 0.5 px) lightly smoothed and joined by a monotone cubic: their curvature jumped at
every reading. Here every design line is
  1. TRACKED on the original sheet (reference/user_blueprint.jpg, red channel, blurred 0.6 px): at every 0.5 px along the
     predicted line the brightness profile across the line is sampled and the stroke's centre found to 0.05 px (the
     brightness-weighted centre of the peak nearest the prediction); points where the stroke is missing, crossed by another
     drawing (the seal, labels) or doubtful are dropped;
  2. FITTED as a smoothing cubic B-spline in arc length (x(t), y(t), uniform knots every KNOT px, Eilers-Marx difference
     penalty of order 3 on the coefficients ~ the curvature's variation), refitted twice with the outliers (> 1.2 px) out;
  3. CHECKED: the largest distance to the tracked stroke (limit 3 px), the curvature comb (no spikes) -> build/checks/lines_fit.json
     and previews/compare/lines_fit_*.png.
The fitted curves are stored in build/fair_curves.json (dense samples every 0.25 px of the sheet, with tangents and
curvature) for gr_shape.py / gr_hull.py, which read them through fair_y(name, x) and fair_curve(name).

    python build/gr_fair.py            (writes fair_curves.json, checks/lines_fit.json, the plots)
"""
import os, json, math
import numpy as np

HERE = os.path.dirname(os.path.abspath(__file__))
ROOT = os.path.dirname(HERE)
REF = os.path.join(ROOT, "reference")
OUT_JSON = os.path.join(HERE, "fair_curves.json")


# ------------------------------------------------------------------ B-spline machinery (uniform knots, cubic)
def _bspl_basis(t, t0, h, n_coef, deriv=0):
    """values (or derivatives) of the n_coef uniform cubic B-splines (knots t0 + h * (k - 3)) at the parameters t"""
    t = np.atleast_1d(np.asarray(t, float))
    u = (t - t0) / h                                   # in knot units; span s = floor(u)
    s = np.clip(np.floor(u).astype(int), 0, n_coef - 4)
    f = u - s
    if deriv == 0:
        b = np.stack([(1 - f) ** 3 / 6, (3 * f ** 3 - 6 * f ** 2 + 4) / 6, (-3 * f ** 3 + 3 * f ** 2 + 3 * f + 1) / 6, f ** 3 / 6], 1)
    elif deriv == 1:
        b = np.stack([-(1 - f) ** 2 / 2, (9 * f ** 2 - 12 * f) / 6, (-9 * f ** 2 + 6 * f + 3) / 6, f ** 2 / 2], 1) / h
    elif deriv == 2:
        b = np.stack([(1 - f), (18 * f - 12) / 6, (-18 * f + 6) / 6, f], 1) / h ** 2
    else:
        b = np.stack([-np.ones_like(f), 3 * np.ones_like(f), -3 * np.ones_like(f), np.ones_like(f)], 1) / h ** 3
    B = np.zeros((len(t), n_coef))
    for k in range(4):
        B[np.arange(len(t)), s + k] = b[:, k]
    return B


class PSpline:
    """a parametric smoothing cubic B-spline (x(t), y(t)) over t in [0, L] (arc length, px); closed curves are periodic"""

    def __init__(self, t, P, knot=4.0, lam=1.0, w=None, closed=False, L=None):
        t = np.asarray(t, float)
        P = np.asarray(P, float)
        self.closed = closed
        self.L = float(L if L is not None else t.max())
        if closed:
            n_int = max(4, int(round(self.L / knot)))
            self.h = self.L / n_int
            self.n_int = n_int
            self.t0 = 0.0
            self.n = n_int + 3
        else:
            n_int = max(1, int(math.ceil(self.L / knot)))
            self.h = self.L / n_int
            self.t0 = 0.0
            self.n = n_int + 3
        B = self._B(t)
        W = np.ones(len(t)) if w is None else np.asarray(w, float)
        D = np.diff(np.eye(self.n if not closed else n_int), 3, axis=0)
        if closed:                                     # periodic: coefficients wrap (n_int free ones)
            Pm = np.zeros((self.n, n_int))
            for k in range(self.n):
                Pm[k, k % n_int] = 1.0
            Bf = B @ Pm
            Dp = np.zeros((n_int, n_int))
            for k in range(n_int):                     # periodic third differences
                for j, c in enumerate((-1, 3, -3, 1)):
                    Dp[k, (k + j) % n_int] += c
            A = Bf.T @ (W[:, None] * Bf) + lam * Dp.T @ Dp
            cf = np.linalg.solve(A, Bf.T @ (W[:, None] * P))
            self.c = Pm @ cf
        else:
            A = B.T @ (W[:, None] * B) + lam * D.T @ D
            self.c = np.linalg.solve(A, B.T @ (W[:, None] * P))

    def _B(self, t, d=0):
        t = np.atleast_1d(np.asarray(t, float))
        if self.closed:
            t = np.mod(t, self.L)
        else:
            t = np.clip(t, 0.0, self.L)
        return _bspl_basis(t, self.t0, self.h, self.n, d)

    def __call__(self, t, d=0):
        return self._B(t, d) @ self.c

    def curvature(self, t):
        d1, d2 = self(t, 1), self(t, 2)
        return (d1[:, 0] * d2[:, 1] - d1[:, 1] * d2[:, 0]) / np.maximum(1e-12, np.hypot(d1[:, 0], d1[:, 1]) ** 3)


# ------------------------------------------------------------------ the sheet
def load_sheet():
    from PIL import Image, ImageFilter
    im = Image.open(os.path.join(REF, "user_blueprint.jpg")).convert("RGB")
    r = np.asarray(im.split()[0].filter(ImageFilter.GaussianBlur(0.6)), float)
    return r


def bilinear(img, x, y):
    x = np.asarray(x, float); y = np.asarray(y, float)
    x0 = np.clip(np.floor(x).astype(int), 0, img.shape[1] - 2)
    y0 = np.clip(np.floor(y).astype(int), 0, img.shape[0] - 2)
    fx, fy = x - x0, y - y0
    return ((1 - fx) * (1 - fy) * img[y0, x0] + fx * (1 - fy) * img[y0, x0 + 1] + (1 - fx) * fy * img[y0 + 1, x0]
            + fx * fy * img[y0 + 1, x0 + 1])


def resample(poly, step, closed=False):
    P = np.asarray(poly, float)
    if closed:
        P = np.vstack([P, P[:1]])
    seg = np.hypot(*np.diff(P, axis=0).T)
    s = np.concatenate([[0], np.cumsum(seg)])
    n = max(2, int(round(s[-1] / step)))
    ts = np.linspace(0, s[-1], n + (0 if closed else 1), endpoint=not closed)
    x = np.interp(ts, s, P[:, 0]); y = np.interp(ts, s, P[:, 1])
    return np.stack([x, y], 1), ts, s[-1]


def track(img, pred, half=1.8, thr=60.0, closed=False, step=0.5, dmax=None):
    """the stroke's centre near each point of the predicted curve pred (dense polyline): brightness profile along the normal,
    the peak nearest the prediction, its brightness-weighted centre. Returns (t, points, ok-mask, offsets)"""
    P = np.asarray(pred, float)
    n = len(P)
    if closed:
        T = np.roll(P, -1, 0) - np.roll(P, 1, 0)
    else:
        T = np.gradient(P, axis=0)
    T /= np.maximum(1e-9, np.hypot(T[:, 0], T[:, 1]))[:, None]
    N = np.stack([-T[:, 1], T[:, 0]], 1)
    offs = np.arange(-half, half + 1e-9, 0.05)
    X = P[:, 0][:, None] + N[:, 0][:, None] * offs[None, :]
    Y = P[:, 1][:, None] + N[:, 1][:, None] * offs[None, :]
    V = bilinear(img, X, Y)
    out = np.zeros(n); ok = np.zeros(n, bool); peak = np.zeros(n)
    for i in range(n):
        v = V[i]
        # local maxima above thr; the one nearest the prediction
        mx = [k for k in range(1, len(v) - 1) if v[k] >= v[k - 1] and v[k] >= v[k + 1] and v[k] > thr]
        if not mx:
            continue
        k = min(mx, key=lambda k: abs(offs[k]))
        half_max = 0.5 * (v[k] + max(0.0, np.percentile(v, 10)))
        a = k
        while a > 0 and v[a - 1] > half_max:
            a -= 1
        b = k
        while b < len(v) - 1 and v[b + 1] > half_max:
            b += 1
        w = v[a:b + 1] - half_max
        if w.sum() <= 0:
            continue
        c = float((offs[a:b + 1] * w).sum() / w.sum())
        if (b - a) * 0.05 > 3.2:                       # a wide blob, not a line (a crossing, a filled shape)
            continue
        if dmax is not None and abs(c) > dmax:
            continue
        out[i] = c; ok[i] = True; peak[i] = v[k]
    pts = P + N * out[:, None]
    return pts, ok, out, peak


def fit_line(img, pred, knot=4.0, lam=2.0, half=1.8, iters=3, closed=False, thr=60.0, exclude=(), keep_ends=True,
             outlier=1.2, step=0.5):
    """track and fit; exclude: (x0, y0, x1, y1) boxes where the stroke is crossed by other drawings (dropped)"""
    P0, ts, L = resample(pred, step, closed)
    cur = P0
    info = {}
    for it in range(iters):
        pts, ok, off, peak = track(img, cur, half=half if it == 0 else min(half, 1.2), thr=thr, closed=closed)
        for (a, b, c, d) in exclude:
            m = (pts[:, 0] >= a) & (pts[:, 0] <= c) & (pts[:, 1] >= b) & (pts[:, 1] <= d)
            ok &= ~m
        # arc-length parameter of the tracked points: along the current curve
        t = ts.copy()
        w = ok.astype(float)
        if keep_ends and not closed:                   # anchor the ends lightly to the prediction (lines that end on others)
            for e in (0, -1):
                if not ok[e]:
                    pts[e] = cur[e]; w[e] = 0.3
        sp = PSpline(t[w > 0], pts[w > 0], knot=knot, lam=lam, w=w[w > 0], closed=closed, L=L)
        fit = sp(t)
        res = np.hypot(*(fit - pts).T)
        if it < iters - 1:
            bad = ok & (res > outlier)
            ok &= ~bad
            # re-track along the fitted curve
            cur = fit
            pts2, _, _, _ = track(img, cur, half=1.0, thr=thr, closed=closed)
        info = dict(n_track=int(ok.sum()), n_all=len(ok), rms=float(np.sqrt(np.mean(res[ok] ** 2))) if ok.any() else None,
                    max=float(res[ok].max()) if ok.any() else None)
    sp.tracked = pts[ok]
    sp.info = info
    return sp


# ------------------------------------------------------------------ the lines
def build():
    import importlib.util
    G = {}
    exec(open(os.path.join(HERE, "gr_measure.py"), encoding="utf-8").read(), G)
    img = load_sheet()
    lines = {}
    URIM, LRIM = G["WAIST_UPPER_RIM"], G["WAIST_LOWER_RIM"]
    FT = [(839.5, 345.0), (845.0, 341.0), (870.0, 339.0), (887.0, 339.0), (893.0, 341.0)]
    # the upper rim, both lines, from the duct to the nose's upper corner (forward of the band the rim runs along the window
    # frame's top: its lower line is the frame's top edge, the upper line ~4 px above it)
    ur_u = [(x, a) for x, a, b in URIM] + [(845.0, 336.5), (870.0, 335.0), (887.0, 335.5)]
    ur_l = [(x, b) for x, a, b in URIM] + [(845.0, 341.0), (870.0, 339.0), (887.0, 339.0)]
    SEAL = (684.0, 334.0, 711.0, 361.0)                # the seal's circle crosses nothing of the rims; its text box does
    lines["ur_u"] = fit_line(img, ur_u, knot=5.0, lam=8.0, exclude=[SEAL])
    lines["ur_l"] = fit_line(img, ur_l, knot=5.0, lam=8.0, exclude=[SEAL])
    lr_u = [(x, a) for x, a, b in LRIM]
    lr_l = [(x, b) for x, a, b in LRIM]
    lines["lr_u"] = fit_line(img, lr_u, knot=4.0, lam=6.0)
    lines["lr_l"] = fit_line(img, lr_l, knot=4.0, lam=6.0)
    # the window frame's outer outline (SIDE_WINDOW_FRAME_OUT) and the side glass (closed)
    lines["frame_out"] = fit_line(img, G["SIDE_WINDOW_FRAME_OUT"], knot=4.0, lam=4.0, closed=True, half=1.5)
    lines["side_glass"] = fit_line(img, G["SIDE_WINDOW_GLASS"], knot=4.0, lam=4.0, closed=True, half=1.5)
    return lines, img, G


def save(lines):
    out = {}
    for k, sp in lines.items():
        n = int(math.ceil(sp.L / 0.25))
        t = np.linspace(0, sp.L, n + (0 if sp.closed else 1), endpoint=not sp.closed)
        P = sp(t); d1 = sp(t, 1); kap = sp.curvature(t)
        out[k] = dict(closed=sp.closed, L=sp.L, knot=sp.h, t=t.round(4).tolist(), x=P[:, 0].round(4).tolist(),
                      y=P[:, 1].round(4).tolist(), tx=d1[:, 0].round(5).tolist(), ty=d1[:, 1].round(5).tolist(),
                      k=kap.round(6).tolist(), fit=sp.info)
    json.dump(out, open(OUT_JSON, "w"))
    return out


if __name__ == "__main__":
    lines, img, G = build()
    for k, sp in lines.items():
        print(k, sp.info, "L", round(sp.L, 1))
    save(lines)


# ------------------------------------------------------------------ plots (the sheet with the fits and their curvature combs)
def plot(lines, img, path, box=(520, 315, 920, 485), K=4, comb=60.0, names=None):
    from PIL import Image, ImageDraw
    x0, y0, x1, y1 = box
    sheet = Image.open(os.path.join(REF, "user_blueprint.jpg")).convert("RGB").crop(box).resize(((x1 - x0) * K, (y1 - y0) * K), Image.LANCZOS)
    sheet = Image.blend(sheet, Image.new("RGB", sheet.size, (0, 0, 0)), 0.45)
    d = ImageDraw.Draw(sheet)
    f = lambda p: ((p[0] - x0) * K, (p[1] - y0) * K)
    for k, sp in lines.items():
        if names and k not in names:
            continue
        n = int(sp.L / 0.5)
        t = np.linspace(0, sp.L, n + 1)
        P = sp(t); d1 = sp(t, 1); kap = sp.curvature(t)
        Tn = d1 / np.maximum(1e-9, np.hypot(d1[:, 0], d1[:, 1]))[:, None]
        N = np.stack([-Tn[:, 1], Tn[:, 0]], 1)
        for i in range(0, n + 1, 2):                    # comb teeth (green) and their tips' envelope
            q = P[i] - N[i] * kap[i] * comb
            d.line([f(P[i]), f(q)], fill=(40, 200, 80), width=1)
        env = [f(P[i] - N[i] * kap[i] * comb) for i in range(n + 1)]
        d.line(env, fill=(120, 255, 140), width=1)
        for p in getattr(sp, "tracked", []):
            u, v = f(p)
            d.point((u, v), fill=(255, 230, 0))
        d.line([f(p) for p in P], fill=(255, 40, 40), width=2)
    sheet.save(path)


if __name__ == "__main__":
    os.makedirs(os.path.join(ROOT, "previews", "compare"), exist_ok=True)
    plot(lines, img, os.path.join(ROOT, "previews", "compare", "lines_fit_side.png"))