#!/usr/bin/env python3
"""Measure a 600 dpi scan of the MeerK40t calibration burn (mk_cal_S*_L*.svg on black paper) and
tell exactly what to change in MeerK40t (lens_size, flip/swap) and build a residual-corrected .cor.

  python3 mk_cal_analyze.py SCAN.png --lens 131.5 --cor-in /home/hank/laser/meerk40t/jcz110.cor \
        [--flip-x 0 --flip-y 1 --swap-xy 1 --rotate 0] [--prev sheet1_result.json] --out OUTDIR

Scan: face down on the glass, colour or grey, 600 dpi, the paper's snipped corner = burn's BACK-LEFT corner.
Nothing here talks to the laser. Grok 2026-10-31 for Hank.
"""
import argparse, json, math, os, struct, sys, time
import numpy as np
import cv2

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import mk_cal_design as D
from mk_view_table import VIEW_MATRIX

PPM_X_DEFAULT = 600 / 25.4          # scanner X (carriage) px/mm, nominal
PPM_Y_DEFAULT = 599.70 / 25.4       # scanner Y (motor) px/mm, measured on Hank's HP 2600 (2026-10)
D4 = [np.array(m, float) for m in ([[1, 0], [0, 1]], [[0, -1], [1, 0]], [[-1, 0], [0, -1]], [[0, 1], [-1, 0]],
                                   [[-1, 0], [0, 1]], [[1, 0], [0, -1]], [[0, 1], [1, 0]], [[0, -1], [-1, 0]])]


def log(*a):
    print(*a, flush=True)


# ----------------------------------------------------------------------------------------- image prep
def load_gray(path):
    img = cv2.imread(path, cv2.IMREAD_UNCHANGED)
    if img is None:
        raise SystemExit(f'cannot read {path}')
    if img.dtype == np.uint16:
        img = (img / 257).astype(np.uint8)
    if img.ndim == 3:
        img = cv2.cvtColor(img[:, :, :3], cv2.COLOR_BGR2GRAY)
    return img


def paper_mask_and_rotation(gray, ppm):
    """dark paper on light lid. returns (mask, k_rot90_ccw so that snipped corner goes to top-left, info)"""
    f = 8
    small = cv2.resize(gray, (gray.shape[1] // f, gray.shape[0] // f), interpolation=cv2.INTER_AREA)
    small = cv2.GaussianBlur(small, (0, 0), 3)
    _, th = cv2.threshold(small, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)
    n, lab, st, _ = cv2.connectedComponentsWithStats(th)
    if n < 2:
        return None, 0, {'paper': 'not found'}
    i = 1 + int(np.argmax(st[1:, cv2.CC_STAT_AREA]))
    m = (lab == i).astype(np.uint8) * 255
    m = cv2.morphologyEx(m, cv2.MORPH_CLOSE, np.ones((9, 9), np.uint8))
    cnts, _ = cv2.findContours(m, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
    c = max(cnts, key=cv2.contourArea)
    hull = cv2.convexHull(c).reshape(-1, 2).astype(float)
    box = cv2.boxPoints(cv2.minAreaRect(c))
    dists = [np.min(np.hypot(*(hull - b).T)) for b in box]
    k = int(np.argmax(dists))
    clip_mm = dists[k] * f / ppm
    info = {'paper_area_mm2': float(cv2.contourArea(c) * f * f / ppm / ppm), 'snip_size_mm': round(float(clip_mm), 1)}
    mask = cv2.resize(m, (gray.shape[1], gray.shape[0]), interpolation=cv2.INTER_NEAREST)
    mask = cv2.erode(mask, np.ones((25, 25), np.uint8))
    if clip_mm < 3.0:
        info['snip'] = 'no snipped corner found - assuming the scan is already back-edge-up'
        return mask, 0, info
    cx, cy = c.reshape(-1, 2).mean(0)
    bx, by = box[k]
    quad = ('T' if by < cy else 'B') + ('L' if bx < cx else 'R')
    krot = {'TL': 0, 'TR': 1, 'BR': 2, 'BL': 3}[quad]
    info['snip'] = f'snipped corner found at scan {quad}; rotating scan {90*krot} deg CCW so it is top-left'
    return mask, krot, info


def strength_image(gray, mask, ppm):
    k = int(round(1.3 * ppm)) | 1
    bg = cv2.medianBlur(gray, k)
    g = gray.astype(np.float32)
    b = bg.astype(np.float32)
    inner = mask > 0
    d1, d2 = g - b, b - g
    p1, p2 = np.percentile(d1[inner], 99.7), np.percentile(d2[inner], 99.7)
    s = d1 if p1 >= p2 else d2
    polarity = 'marks LIGHTER than paper' if p1 >= p2 else 'marks DARKER than paper'
    s[~inner] = 0
    np.maximum(s, 0, out=s)
    vals = s[inner]
    med = float(np.median(vals))
    sigma = float(1.4826 * np.median(np.abs(vals - med))) + 0.5
    return s, polarity, sigma


# ----------------------------------------------------------------------------------------- coarse registration
def find_outer_square(s, sigma, ppm):
    thr = max(5 * sigma, 0.2 * float(np.percentile(s[s > 0], 99.9)) if np.any(s > 0) else 10)
    b = (s > thr).astype(np.uint8) * 255
    b = cv2.morphologyEx(b, cv2.MORPH_CLOSE, np.ones((7, 7), np.uint8))
    cnts, _ = cv2.findContours(b, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    best = None
    for c in cnts:
        (cx, cy), (w, h), a = cv2.minAreaRect(c)
        w, h = w / ppm, h / ppm
        if 60 < w < 145 and 60 < h < 145 and max(w, h) / min(w, h) < 1.2:
            if best is None or cv2.contourArea(c) > cv2.contourArea(best):
                best = c
    if best is None:
        raise SystemExit('outer 100 mm square not found in the scan (check polarity / that the burn is visible)')
    hull = cv2.convexHull(best)
    peri = cv2.arcLength(hull, True)
    quad = None
    for e in np.linspace(0.005, 0.08, 40):
        ap = cv2.approxPolyDP(hull, e * peri, True)
        if len(ap) == 4:
            quad = ap.reshape(4, 2).astype(float)
            break
    if quad is None:
        quad = cv2.boxPoints(cv2.minAreaRect(best))
    ctr = quad.mean(0)
    ang = np.arctan2(quad[:, 1] - ctr[1], quad[:, 0] - ctr[0])
    quad = quad[np.argsort(ang)]          # clockwise on screen (y down), starting near -180deg (left/top-left)
    return quad, thr


def homography_from(design_pts, img_pts):
    H, _ = cv2.findHomography(np.float64(design_pts), np.float64(img_pts), 0)
    return H


def apply_h(H, p):
    p = np.asarray(p, float)
    q = np.c_[p, np.ones(len(p))] @ H.T
    return q[:, :2] / q[:, 2:3]


def sample(s, pts):
    """bilinear sample s at float pixel coords pts (...,2)"""
    sh = pts.shape[:-1]
    mx = pts[..., 0].ravel().astype(np.float32)
    my = pts[..., 1].ravel().astype(np.float32)
    n = len(mx)
    w = 4096
    pad = (-n) % w
    mx = np.pad(mx, (0, pad), constant_values=-10).reshape(-1, w)
    my = np.pad(my, (0, pad), constant_values=-10).reshape(-1, w)
    v = cv2.remap(s, mx, my, cv2.INTER_LINEAR, borderMode=cv2.BORDER_CONSTANT, borderValue=0)
    return v.ravel()[:n].reshape(sh)


def disk_max(s, f, centre, r=1.1):
    a = np.linspace(0, 2 * np.pi, 24, endpoint=False)
    pts = [centre]
    for rr in (0.4, 0.8, r):
        pts += [(centre[0] + rr * math.cos(t), centre[1] + rr * math.sin(t)) for t in a]
    return float(np.max(sample(s, f(np.array(pts)))))


def choose_orientation(s, quad, thr):
    dc = np.array([[-50, -50], [50, -50], [50, 50], [-50, 50]], float)   # TL TR BR BL clockwise
    best = None
    for rev in (False, True):
        for sh in range(4):
            q = np.roll(quad, -sh, axis=0)
            if rev:
                q = q[[0, 3, 2, 1]]
            H = homography_from(dc, q)
            f = lambda p, H=H: apply_h(H, p)
            score = 0
            for key, pos, present in D.pip_positions():
                m = disk_max(s, f, pos)
                score += (1 if (m > thr) == present else -1)
            if best is None or score > best[0]:
                best = (score, H, sh, rev)
    return best


# ----------------------------------------------------------------------------------------- feature sampling
def features():
    """list of (name, design point (2,), design unit normal (2,)) for every profile to measure"""
    F = []
    h = D.OUTER / 2
    for t in np.arange(-47, 47.01, 1.0):
        F += [('outer_top', (t, -h), (0, 1)), ('outer_bottom', (t, h), (0, 1)),
              ('outer_left', (-h, t), (1, 0)), ('outer_right', (h, t), (1, 0))]
    hi = D.INNER / 2
    for t in np.arange(-36, 36.01, 1.0):
        if abs(t) < 2.5:
            continue
        F += [('inner_top', (t, -hi), (0, 1)), ('inner_bottom', (t, hi), (0, 1)),
              ('inner_left', (-hi, t), (1, 0)), ('inner_right', (hi, t), (1, 0))]
    for t in np.arange(-36, 36.01, 1.0):           # along diagonals, t = x
        r = abs(t) * math.sqrt(2)
        if r < 2.5 or abs(r - D.R_CEN) < 2 or abs(r - D.R_BIG) < 2:
            continue
        n = 1 / math.sqrt(2)
        F += [('diag_main', (t, t), (n, -n)), ('diag_anti', (t, -t), (n, n))]
    avoid = [0] + list(D.TICKS) + [D.R_CEN, -D.R_CEN, D.R_BIG, -D.R_BIG, hi, -hi]
    for t in np.arange(-45, 45.01, 0.5):
        if min(abs(t - a) for a in avoid) < (2.5 if abs(t) < 1 else 1.4):
            continue
        F += [('axis_x', (t, 0), (0, 1)), ('axis_y', (0, t), (1, 0))]
    for tk in D.TICKS:
        for o in (-1.3, 1.3):
            F += [(f'tick_x{tk:+d}', (tk, o), (1, 0)), (f'tick_y{tk:+d}', (o, tk), (0, 1))]
    for name, r, cx, cy, step, skip in [('circle_big', D.R_BIG, 0, 0, 2, True), ('circle_centre', D.R_CEN, 0, 0, 4, True)] + \
            [(f'circle_corner_{"T" if sy < 0 else "B"}{"L" if sx < 0 else "R"}', D.R_CORNER, sx * D.CORNER_C, sy * D.CORNER_C, 6, False)
             for sx in (-1, 1) for sy in (-1, 1)]:
        for a in np.arange(0, 360, step):
            if skip and min(abs(((a - b) + 180) % 360 - 180) for b in range(0, 360, 45)) < 5:
                continue
            ar = math.radians(a)
            F += [(name, (cx + r * math.cos(ar), cy + r * math.sin(ar)), (math.cos(ar), math.sin(ar)))]
    names = [f[0] for f in F]
    P = np.array([f[1] for f in F], float)
    N = np.array([f[2] for f in F], float)
    return names, P, N


def measure(s, f, P, N, half, sigma):
    """profile across each feature through mapping f (design mm -> px). returns image pts (n,2), ok mask, peak"""
    st = np.arange(-half, half + 1e-9, 0.025)
    pts = P[:, None, :] + st[None, :, None] * N[:, None, :]
    px = f(pts.reshape(-1, 2)).reshape(len(P), len(st), 2)
    prof = sample(s, px)
    base = np.median(prof, axis=1, keepdims=True)
    v = prof - base
    peak = v.max(1)
    w = np.clip(v - 0.45 * peak[:, None], 0, None)
    sw = w.sum(1) + 1e-9
    sstar = (w * st[None, :]).sum(1) / sw
    imax = v.argmax(1)
    ok = (peak > max(6 * sigma, 6)) & (imax > 2) & (imax < len(st) - 3)
    # reject wide / double responses
    width = (w > 0).sum(1) * 0.025
    ok &= width < 0.8
    q = f(P + sstar[:, None] * N)
    return q, ok, peak


# ----------------------------------------------------------------------------------------- models
def poly_terms(p, deg=3):
    x, y = p[:, 0], p[:, 1]
    cols = [np.ones_like(x)]
    for d in range(1, deg + 1):
        for i in range(d, -1, -1):
            cols.append(x ** i * y ** (d - i))
    return np.stack(cols, 1)


def fit_normal_poly(P, Nq, Q, deg=3, scale=50.0, w=None):
    """fit F(p)=[Px(p),Py(p)] (polynomial in p/scale) minimising sum (nq.(F(p)-q))^2. returns coef (2,nt)"""
    T = poly_terms(P / scale, deg)
    A = np.concatenate([T * Nq[:, :1], T * Nq[:, 1:2]], 1)
    b = (Nq * Q).sum(1)
    if w is not None:
        A, b = A * w[:, None], b * w
    c, *_ = np.linalg.lstsq(A, b, rcond=None)
    nt = T.shape[1]
    return np.stack([c[:nt], c[nt:]])


def eval_poly(coef, P, scale=50.0):
    T = poly_terms(np.asarray(P, float) / scale, int(round((math.sqrt(8 * coef.shape[1] + 1) - 3) / 2)))
    return T @ coef.T


def robust_poly(P, Nq, Q, deg=3, iters=3):
    keep = np.ones(len(P), bool)
    for _ in range(iters):
        coef = fit_normal_poly(P[keep], Nq[keep], Q[keep], deg)
        r = ((eval_poly(coef, P) - Q) * Nq).sum(1)
        sd = 1.4826 * np.median(np.abs(r[keep])) + 1e-6
        keep = np.abs(r) < max(4 * sd, 1e-3)
    return coef, keep, r


def fit_rigid(P, Nq, Q, allow_scale=False, iters=8):
    """q ~ s R(th) p + T using normal residuals (Gauss-Newton)"""
    th, s, T = 0.0, 1.0, Q.mean(0) - P.mean(0)
    for _ in range(iters):
        c, si = math.cos(th), math.sin(th)
        R = np.array([[c, -si], [si, c]])
        Rp = P @ R.T
        F = s * Rp + T
        r = ((F - Q) * Nq).sum(1)
        dR = np.array([[-si, -c], [c, -si]])
        J = [s * ((P @ dR.T) * Nq).sum(1), Nq[:, 0], Nq[:, 1]]
        if allow_scale:
            J.append((Rp * Nq).sum(1))
        J = np.stack(J, 1)
        d, *_ = np.linalg.lstsq(J, -r, rcond=None)
        th += d[0]
        T = T + d[1:3]
        if allow_scale:
            s += d[3]
    c, si = math.cos(th), math.sin(th)
    R = np.array([[c, -si], [si, c]])
    r = ((s * P @ R.T + T - Q) * Nq).sum(1)
    return th, s, T, r


# ----------------------------------------------------------------------------------------- .cor IO
def read_cor(path):
    b = open(path, 'rb').read()
    label = b[:0x16]
    try:
        lab = label.decode('utf-16')
    except UnicodeDecodeError:
        lab = label.decode('utf-8', errors='ignore')
    if lab == 'LMC1COR_1.0':
        raise SystemExit('float LMC1COR_1.0 .cor not supported by this tool (JCZ_COR_2_1 int format expected)')
    head = b[:0x24]
    scale = struct.unpack('d', b[0x1c:0x24])[0]
    tab = np.frombuffer(b[0x24:0x24 + 65 * 65 * 8], '<i4').reshape(65, 65, 2).astype(float)
    tail = b[0x24 + 65 * 65 * 8:]
    return head, tab, tail, scale


def write_cor(path, head, tab, tail):
    t = np.round(tab).astype(np.int64)
    if np.abs(t).max() >= 0x8000:
        raise SystemExit('correction too large for .cor (|value| >= 32768) - refusing to write')
    open(path, 'wb').write(head + t.astype('<i4').tobytes() + tail)


# ----------------------------------------------------------------------------------------- main analysis
def describe_S(S):
    names = {
        ((1, 0), (0, 1)): 'correct (not mirrored, not rotated)',
        ((0, -1), (1, 0)): 'rotated 90 deg clockwise', ((-1, 0), (0, -1)): 'rotated 180 deg',
        ((0, 1), (-1, 0)): 'rotated 90 deg counter-clockwise',
        ((-1, 0), (0, 1)): 'MIRRORED left-right', ((1, 0), (0, -1)): 'MIRRORED top-bottom',
        ((0, 1), (1, 0)): 'MIRRORED across the TL-BR diagonal (x/y swapped)',
        ((0, -1), (-1, 0)): 'MIRRORED across the TR-BL diagonal'}
    return names[tuple(tuple(int(v) for v in row) for row in S)]


def flags_for(M):
    out = []
    for k, v in VIEW_MATRIX.items():
        if np.array_equal(np.array(v), M):
            out.append(k)
    out.sort(key=lambda k: (k[3] != 0, k))
    return out


def eps_at_ticks(res, tx, ty):
    """residual distortion epsilon (mm, physical paper frame, rotation removed) at galvo ticks (input side)"""
    L = res['config']['lens_mm']
    k = 65535.0 / L
    R = np.array(res['config']['view_matrix'], float)
    t = np.stack([tx - 32768.0, ty - 32768.0], -1) / k
    p = t @ R            # R^T applied: p = R^T t  (row vectors)
    p = np.clip(p, -52, 52)
    coef = np.array(res['eps_coef'])
    return eval_poly(coef, p.reshape(-1, 2)).reshape(p.shape)


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('scan')
    ap.add_argument('--lens', type=float, required=True, help='MeerK40t lens_size used for this burn (mm)')
    ap.add_argument('--cor-in', help='the .cor that was ACTIVE for this burn (JCZ_COR_2_1)')
    ap.add_argument('--flip-x', type=int, default=0)
    ap.add_argument('--flip-y', type=int, default=1)
    ap.add_argument('--swap-xy', type=int, default=1)
    ap.add_argument('--rotate', type=int, default=0)
    ap.add_argument('--ppm-x', type=float, default=PPM_X_DEFAULT, help='scanner px/mm across (carriage)')
    ap.add_argument('--ppm-y', type=float, default=PPM_Y_DEFAULT, help='scanner px/mm along (motor)')
    ap.add_argument('--scan-rotate', type=int, default=None, help='force: rotate scan CCW by this (0/90/180/270)')
    ap.add_argument('--gain', type=float, default=1.0, help='fraction of the measured residual to put in the new .cor')
    ap.add_argument('--cor-sign', type=float, default=1.0, help='table += sign*... (+1: board out = t - table)')
    ap.add_argument('--prev', help='result JSON of the previous sheet (whose .cor was burned now) -> gain/sign check')
    ap.add_argument('--tag', default=None)
    ap.add_argument('--out', default='.')
    a = ap.parse_args()
    os.makedirs(a.out, exist_ok=True)
    tag = a.tag or os.path.splitext(os.path.basename(a.scan))[0]
    t0 = time.time()
    gray = load_gray(a.scan)
    ppm0 = (a.ppm_x + a.ppm_y) / 2
    mask, krot, pinfo = paper_mask_and_rotation(gray, ppm0)
    if a.scan_rotate is not None:
        krot = (a.scan_rotate // 90) % 4
        pinfo['snip'] = f'forced scan rotation {a.scan_rotate}'
    if mask is None:
        mask = np.full(gray.shape, 255, np.uint8)
    gray = np.ascontiguousarray(np.rot90(gray, krot))
    mask = np.ascontiguousarray(np.rot90(mask, krot))
    ppmx, ppmy = (a.ppm_x, a.ppm_y) if krot % 2 == 0 else (a.ppm_y, a.ppm_x)
    log(f'[scan] {gray.shape[1]}x{gray.shape[0]} px, {pinfo.get("snip")}')
    s, polarity, sigma = strength_image(gray, mask, ppm0)
    log(f'[scan] {polarity}, noise sigma {sigma:.1f}')
    quad, thr = find_outer_square(s, sigma, ppm0)
    score, H, sh, rev = choose_orientation(s, quad, thr)
    log(f'[orient] pip score {score}/12')
    names, P, N = features()
    names = np.array(names)
    # pass 1 (homography), pass 2+3 (cubic image model)
    f = lambda p: apply_h(H, p)
    half = 1.6
    for it in range(3):
        Q, ok, peak = measure(s, f, P, N, half, sigma)
        # image-space normals of each feature
        dn = f(P + 0.05 * N) - f(P)
        Nq = dn / np.linalg.norm(dn, axis=1, keepdims=True)
        coef_img, keep, _ = robust_poly(P[ok], Nq[ok], Q[ok], 3)
        f = lambda p, c=coef_img: eval_poly(c, p)
        half = 0.8
    ok_idx = np.where(ok)[0][keep]
    log(f'[measure] {len(ok_idx)} / {len(P)} profile points used')
    # physical mm (paper frame, y down)
    Pd = P[ok_idx]
    Nd = N[ok_idx]
    Qmm = Q[ok_idx] / np.array([ppmx, ppmy])
    dn = (f(Pd + 0.05 * Nd) - f(Pd)) / np.array([ppmx, ppmy])
    Nm = dn / np.linalg.norm(dn, axis=1, keepdims=True)
    # orientation S from the jacobian of the image model at the centre
    J = (f(np.array([[1.0, 0], [0, 1.0]])) - f(np.array([[0.0, 0]]))) / np.array([ppmx, ppmy])
    J = J.T   # columns = images of design e1, e2
    S = max(D4, key=lambda M: float(np.trace(M.T @ J)))
    Pp = Pd @ S.T            # design expressed in the physical orientation of the burn
    # fits
    th_r, _, T_r, r_rigid = fit_rigid(Pp, Nm, Qmm)
    th_s, s_sim, T_s, r_sim = fit_rigid(Pp, Nm, Qmm, allow_scale=True)
    A_aff = fit_normal_poly(Pp, Nm, Qmm, 1)
    r_aff = ((eval_poly(A_aff, Pp) - Qmm) * Nm).sum(1)
    coef3, keep3, r_cub = robust_poly(Pp, Nm, Qmm, 3)
    A = A_aff[:, 1:3] / 50.0     # linear part (per mm)
    U, sv, Vt = np.linalg.svd(A)
    Rpol = U @ Vt
    Psym = Rpol.T @ A
    rot_deg = math.degrees(math.atan2(Rpol[1, 0], Rpol[0, 0]))
    sx, sy, shear = Psym[0, 0], Psym[1, 1], Psym[0, 1]
    ex, ey = A[:, 0], A[:, 1]
    skew_deg = math.degrees(math.acos(np.dot(ex, ey) / np.linalg.norm(ex) / np.linalg.norm(ey))) - 90
    s_mean = math.sqrt(abs(np.linalg.det(Psym)))
    # model-based corner / side quantities in the rotation-removed physical frame (origin = centre cross)
    Fm = lambda p: (eval_poly(coef3, np.asarray(p, float)) - eval_poly(coef3, np.zeros((1, 2)))) @ Rpol
    C = Fm([[-50, -50], [50, -50], [50, 50], [-50, 50]])
    top, bot = np.linalg.norm(C[1] - C[0]), np.linalg.norm(C[2] - C[3])
    lef, rig = np.linalg.norm(C[3] - C[0]), np.linalg.norm(C[2] - C[1])
    trap_x = (top - bot) / ((top + bot) / 2) * 100
    trap_y = (lef - rig) / ((lef + rig) / 2) * 100
    sq_ctr = C.mean(0)
    centre_offset = -sq_ctr          # centre cross relative to the outer-square centre
    mids = Fm([[0, -50], [50, 0], [0, 50], [-50, 0]])
    sag = {}
    for nm, m, c1, c2, outn in [('top', mids[0], C[0], C[1], (0, -1)), ('right', mids[1], C[1], C[2], (1, 0)),
                                ('bottom', mids[2], C[2], C[3], (0, 1)), ('left', mids[3], C[3], C[0], (-1, 0))]:
        sag[nm] = float(np.dot(m - (c1 + c2) / 2, outn))
    # scale on centre line (ruler) vs at the edges
    ax = Fm([[-44, 0], [44, 0], [0, -44], [0, 44]])
    ruler_x = np.linalg.norm(ax[1] - ax[0]) / 88
    ruler_y = np.linalg.norm(ax[3] - ax[2]) / 88
    edge_x = (top + bot) / 200
    edge_y = (lef + rig) / 200
    # residual field epsilon (rotation & mean scale & translation removed) on a grid -> worst points
    g = np.linspace(-50, 50, 21)
    GX, GY = np.meshgrid(g, g)
    G = np.c_[GX.ravel(), GY.ravel()]
    rigid_err = Fm(G) - G
    rigid_err -= rigid_err.mean(0)
    worst_rigid = float(np.max(np.linalg.norm(rigid_err, axis=1)))
    eps = Fm(G) - s_mean * G
    eps -= eps.mean(0)
    worst_scaled = float(np.max(np.linalg.norm(eps, axis=1)))
    # circles (direct ellipse fit on measured points, rotation removed)
    circles = {}
    for cname in sorted(set(n for n in names[ok_idx] if n.startswith('circle'))):
        sel = names[ok_idx] == cname
        pts = ((Qmm[sel] - T_r) @ np.array([[math.cos(th_r), -math.sin(th_r)], [math.sin(th_r), math.cos(th_r)]]))
        if sel.sum() >= 12:
            (ecx, ecy), (d1, d2), ang = cv2.fitEllipse((pts * 100).astype(np.float32))
            nom = {'circle_big': 2 * D.R_BIG, 'circle_centre': 2 * D.R_CEN}.get(cname, 2 * D.R_CORNER)
            circles[cname] = {'nominal_d': nom, 'd_major': round(max(d1, d2) / 100, 3), 'd_minor': round(min(d1, d2) / 100, 3),
                              'ovality_mm': round(abs(d1 - d2) / 100, 3), 'angle_deg': round(ang, 1), 'n': int(sel.sum())}
    # tick positions on the rulers (along-axis error after similarity fit)
    ticks = {}
    for tk in D.TICKS:
        for axn, sel_name in (('x', f'tick_x{tk:+d}'), ('y', f'tick_y{tk:+d}')):
            sel = names[ok_idx] == sel_name
            if sel.sum() == 0:
                continue
            pred = s_sim * Pp[sel] @ np.array([[math.cos(th_s), -math.sin(th_s)], [math.sin(th_s), math.cos(th_s)]]).T + T_s
            r = ((Qmm[sel] - pred) * Nm[sel]).sum(1).mean()
            ticks[f'{axn}{tk:+d}'] = round(float(r), 3)
    # weak-corner check: median peak strength per region
    pk = peak[ok_idx]
    region = {}
    for nm, (cx, cy) in {'TL': (-1, -1), 'TR': (1, -1), 'BL': (-1, 1), 'BR': (1, 1)}.items():
        sel = (Pp[:, 0] * cx > 25) & (Pp[:, 1] * cy > 25)
        region[nm] = round(float(np.median(pk[sel])), 1) if sel.any() else None
    region['centre'] = round(float(np.median(pk[(np.abs(Pp[:, 0]) < 20) & (np.abs(Pp[:, 1]) < 20)])), 1)
    # --------------------------------------------------------------------------- MeerK40t recommendations
    key = (bool(a.flip_x), bool(a.flip_y), bool(a.swap_xy), int(a.rotate))
    R1 = np.array(VIEW_MATRIX[key], float)
    R2 = R1 @ S.T                     # config that makes the burn come out un-mirrored / un-rotated
    new_flags = flags_for(R2.round().astype(int))
    new_lens = a.lens * s_mean
    # epsilon as polynomial in sheet-1 *design* coordinates p (not Pp): eps_phys(p) = Fm(p) - s_mean*S p (then remove mean)
    Pg = G @ S                        # design points p = S^T g whose physical-orientation coords are the grid G
    eps_g = Fm(G) - s_mean * G
    eps_g -= eps_g.mean(0)
    T_eps = poly_terms(Pg / 50.0, 3)
    eps_coef = np.linalg.lstsq(T_eps, eps_g, rcond=None)[0].T
    res = {
        'tag': tag, 'scan': os.path.abspath(a.scan), 'time': time.strftime('%Y-%m-%d %H:%M:%S %Z'),
        'config': {'lens_mm': a.lens, 'flip_x': key[0], 'flip_y': key[1], 'swap_xy': key[2], 'rotate': key[3],
                   'view_matrix': R1.tolist(), 'cor_in': a.cor_in},
        'scan_info': dict(pinfo, polarity=polarity, noise_sigma=round(sigma, 2), ppm_x=ppmx, ppm_y=ppmy,
                          scan_rot_ccw_deg=90 * krot, pip_score=f'{score}/12', points=int(len(ok_idx))),
        'orientation': {'S': S.astype(int).tolist(), 'burn_appears': describe_S(S),
                        'recommended_flags(flip_x,flip_y,swap_xy,rotate)': new_flags[:3]},
        'scale': {'x_pct': round((sx - 1) * 100, 3), 'y_pct': round((sy - 1) * 100, 3), 'mean_pct': round((s_mean - 1) * 100, 3),
                  'ruler_x_pct': round((ruler_x - 1) * 100, 3), 'ruler_y_pct': round((ruler_y - 1) * 100, 3),
                  'edge_x_pct': round((edge_x - 1) * 100, 3), 'edge_y_pct': round((edge_y - 1) * 100, 3)},
        'skew_deg': round(skew_deg, 4), 'shear_sym': round(float(shear), 5), 'paper_rotation_deg': round(rot_deg, 3),
        'trapezoid_pct': {'x_top_minus_bottom': round(trap_x, 3), 'y_left_minus_right': round(trap_y, 3)},
        'outer_square_mm': {'top': round(top, 3), 'bottom': round(bot, 3), 'left': round(lef, 3), 'right': round(rig, 3)},
        'side_bow_mm(+ = bows outward/barrel, - = inward/pincushion)': {k: round(v, 3) for k, v in sag.items()},
        'centre_cross_offset_from_square_centre_mm': [round(float(v), 3) for v in centre_offset],
        'circles': circles, 'tick_error_mm': ticks, 'mark_strength_by_region': region,
        'fit_residual_rms_mm': {'rigid': round(float(np.sqrt(np.mean(r_rigid ** 2))), 4), 'similarity': round(float(np.sqrt(np.mean(r_sim ** 2))), 4),
                                'affine': round(float(np.sqrt(np.mean(r_aff ** 2))), 4), 'cubic': round(float(np.sqrt(np.mean(r_cub[keep3] ** 2))), 4)},
        'worst_point_mm': {'rigid_only(true geometry)': round(worst_rigid, 3), 'after_lens_fix_only': round(worst_scaled, 3)},
        'recommend': {'lens_size_mm': round(new_lens, 3)},
        'eps_coef': eps_coef.tolist(),
    }
    # --------------------------------------------------------------------------- residual .cor
    if a.cor_in:
        head, tab, tail, hscale = read_cor(a.cor_in)
        kk, jj = np.meshgrid(np.arange(65), np.arange(65))      # tab[j,k] : j = galvo y index, k = galvo x index
        tx, ty = kk * 1024.0, jj * 1024.0
        gain = a.gain
        note = ''
        if a.prev:
            prev = json.load(open(a.prev))
            e_prev = eps_at_ticks(prev, tx, ty)
            e_now = eps_at_ticks(res, tx, ty)
            g_prev = prev.get('cor_built', {}).get('gain', 1.0)
            gamma = float(np.sum((e_prev - e_now) * e_prev) / (g_prev * np.sum(e_prev * e_prev) + 1e-12))
            res['prev_check'] = {'prev': a.prev, 'prev_gain': g_prev, 'response_gamma': round(gamma, 3),
                                 'meaning': '1.0 = correction acted exactly as predicted; ~0 = ignored/transposed; <0 = SIGN INVERTED'}
            log(f'[prev] correction response gamma = {gamma:.3f} (1 = as predicted)')
            if gamma < 0:
                note = 'previous correction made things WORSE: sign inverted -> this build uses the opposite sign'
                a.cor_sign = -a.cor_sign * (prev.get('cor_built', {}).get('sign', 1.0))
                gamma = -gamma
            if gamma < 0.3:
                note += ' | response too weak - not building an automatic .cor; check transposition'
                gain = 0.0
            else:
                gain = a.gain / min(max(gamma, 0.5), 1.5)
        k2 = 65535.0 / new_lens
        Rinv = R1 @ S.T                                     # galvo_delta = R1 S^T phys_delta
        e = eps_at_ticks(res, tx, ty)                       # (65,65,2) mm, physical
        dtick = k2 * (e @ Rinv.T)                           # (65,65,2) galvo ticks
        new_tab = tab + a.cor_sign * gain * dtick
        base = os.path.join(a.out, f'{tag}_corrected_L{new_lens:.3f}.cor')
        if gain > 0:
            write_cor(base, head, new_tab, tail)
            half_path = base.replace('.cor', '_half.cor')
            write_cor(half_path, head, tab + a.cor_sign * 0.5 * gain * dtick, tail)
            res['cor_built'] = {'path': os.path.abspath(base), 'half_strength_path': os.path.abspath(half_path), 'gain': gain,
                                'sign': a.cor_sign, 'from': os.path.abspath(a.cor_in), 'max_change_ticks': round(float(np.abs(gain * dtick).max()), 1),
                                'header_scale': hscale, 'note': note,
                                'IMPORTANT': f'after selecting this .cor in MeerK40t the GUI resets Lens Size to {65536/hscale:.3f} - type {new_lens:.3f}mm back in'}
        else:
            res['cor_built'] = {'path': None, 'note': note}
    # --------------------------------------------------------------------------- outputs
    jp = os.path.join(a.out, f'{tag}_result.json')
    json.dump(res, open(jp, 'w'), indent=1)
    write_overlay(os.path.join(a.out, f'{tag}_overlay.jpg'), gray, Q[ok_idx], r_rigid, Nm, ppmx, ppmy, res)
    write_report(os.path.join(a.out, f'{tag}_REPORT.txt'), res)
    log(open(os.path.join(a.out, f'{tag}_REPORT.txt')).read())
    log(f'[done] {time.time()-t0:.1f}s -> {jp}')


def write_overlay(path, gray, Qpx, r, Nm, ppmx, ppmy, res):
    f = 0.25
    img = cv2.cvtColor(cv2.resize(gray, None, fx=f, fy=f, interpolation=cv2.INTER_AREA), cv2.COLOR_GRAY2BGR)
    img = cv2.normalize(img, None, 0, 255, cv2.NORM_MINMAX)
    for q, rr, n in zip(Qpx, r, Nm):
        col = (0, 255, 0) if abs(rr) < 0.1 else (0, 200, 255) if abs(rr) < 0.25 else (0, 0, 255)
        p0 = (int(q[0] * f), int(q[1] * f))
        p1 = (int((q[0] - n[0] * rr * 5 * ppmx) * f), int((q[1] - n[1] * rr * 5 * ppmy) * f))
        cv2.line(img, p0, p1, col, 1, cv2.LINE_AA)
        cv2.circle(img, p0, 1, col, -1)
    y = 30
    for line in [f"{res['tag']}  rigid worst {res['worst_point_mm']['rigid_only(true geometry)']} mm  (arrows = rigid residual x5)",
                 f"scale X {res['scale']['x_pct']}%  Y {res['scale']['y_pct']}%  skew {res['skew_deg']} deg  -> lens {res['recommend']['lens_size_mm']}",
                 f"orientation: {res['orientation']['burn_appears']}"]:
        cv2.putText(img, line, (20, y), cv2.FONT_HERSHEY_SIMPLEX, 0.6, (255, 255, 0), 2, cv2.LINE_AA)
        y += 28
    cv2.imwrite(path, img, [cv2.IMWRITE_JPEG_QUALITY, 85])


def write_report(path, r):
    c = r['config']
    L = []
    L.append(f"MeerK40t calibration scan report - {r['tag']} - {r['time']}")
    L.append(f"burned with: lens_size {c['lens_mm']} mm, flip_x {c['flip_x']}, flip_y {c['flip_y']}, swap_xy {c['swap_xy']}, rotate {c['rotate']}, cor {c['cor_in']}")
    si = r['scan_info']
    L.append(f"scan: {si.get('snip')}; {si['polarity']}; {si['points']} points; pips {si['pip_score']}")
    L.append('')
    L.append(f"ORIENTATION: burn appears {r['orientation']['burn_appears']}")
    fl = r['orientation']['recommended_flags(flip_x,flip_y,swap_xy,rotate)']
    if 'no snipped corner' in str(si.get('snip')):
        L.append('  WARNING: no snipped paper corner found, so a ROTATION of the burn cannot be told apart from how the paper lay on the glass.')
        L.append('  Mirror detection is still valid (det %+d). Do NOT change flip/swap from this sheet unless it says MIRRORED; re-scan with the snip visible.' % int(round(np.linalg.det(np.array(r['orientation']['S'])))))
        L.append('  (The .cor and lens numbers do not depend on this and are still valid.)')
    if r['orientation']['S'] != [[1, 0], [0, 1]]:
        L.append(f"  -> change MeerK40t Device config to flip_x={fl[0][0]}, flip_y={fl[0][1]}, swap_xy={fl[0][2]}, rotate={fl[0][3]}")
    else:
        L.append('  -> no flip/swap change needed')
    s = r['scale']
    L.append(f"SCALE: X {s['x_pct']:+.3f} %   Y {s['y_pct']:+.3f} %   mean {s['mean_pct']:+.3f} %  (+ = burn too big)")
    L.append(f"  centre-line rulers X {s['ruler_x_pct']:+.3f} %  Y {s['ruler_y_pct']:+.3f} % ;  square edges X {s['edge_x_pct']:+.3f} %  Y {s['edge_y_pct']:+.3f} %")
    L.append(f"  -> MeerK40t Lens Size: {c['lens_mm']} -> {r['recommend']['lens_size_mm']} mm   (only global-scale knob MeerK40t has)")
    L.append(f"SKEW: {r['skew_deg']:+.4f} deg (angle between burned X and Y axes minus 90)")
    t = r['trapezoid_pct']
    L.append(f"TRAPEZOID: top-bottom {t['x_top_minus_bottom']:+.3f} %   left-right {t['y_left_minus_right']:+.3f} %")
    o = r['outer_square_mm']
    L.append(f"  outer square sides (nominal 100, rotation removed): top {o['top']} bottom {o['bottom']} left {o['left']} right {o['right']}")
    L.append(f"BOW of outer sides (mm, + outward = barrel, - inward = pincushion): {r['side_bow_mm(+ = bows outward/barrel, - = inward/pincushion)']}")
    L.append(f"CENTRE: cross is {r['centre_cross_offset_from_square_centre_mm']} mm (x,y) from the outer-square centre")
    for k, v in r['circles'].items():
        L.append(f"CIRCLE {k}: nominal {v['nominal_d']} -> {v['d_major']} x {v['d_minor']} mm (ovality {v['ovality_mm']} at {v['angle_deg']} deg)")
    L.append(f"RULER tick errors after best similarity fit (mm): {r['tick_error_mm']}")
    L.append(f"MARK STRENGTH by region (weak-corner check): {r['mark_strength_by_region']}")
    L.append(f"FIT RMS (mm): {r['fit_residual_rms_mm']}")
    L.append(f"WORST POINT (mm, inside the 100 mm square): rigid-only {r['worst_point_mm']['rigid_only(true geometry)']} ; after lens fix only {r['worst_point_mm']['after_lens_fix_only']}")
    if 'prev_check' in r:
        L.append(f"PREVIOUS CORRECTION RESPONSE gamma = {r['prev_check']['response_gamma']} ({r['prev_check']['meaning']})")
    cb = r.get('cor_built')
    if cb and cb.get('path'):
        L.append(f"NEW .cor: {cb['path']}  (max table change {cb['max_change_ticks']} ticks, gain {cb['gain']:.2f}, sign {cb['sign']:+.0f})")
        L.append(f"  half-strength alternative: {cb['half_strength_path']}")
        L.append(f"  {cb['IMPORTANT']}")
        if cb.get('note'):
            L.append('  NOTE: ' + cb['note'])
    elif cb:
        L.append('NO .cor built: ' + cb.get('note', ''))
    open(path, 'w').write('\n'.join(L) + '\n')


if __name__ == '__main__':
    main()