"""HAND-SHAPE SOLVER — bakes hand poses as glTF morph targets (blend shapes).

My own method for the ariki-game hand-pose problem (2026-08-17), an alternative to the
bone-rotation + finger-weight-rebake lane that stalled at exp05:

  A hand pose is explicit SURFACE DEFORMATION, not skeleton rotation.

Why: the scan mesh has ~2.6k inter-digit bridge edges (web remnants). Weights only
choose which bone drags a shared vertex; when adjacent digits curl apart (fist/grip)
those bridges must tear — five weight-rebake iterations could not and cannot fix that.
Baking the pose offline turns the fight into ONE smoothable stretch; tearing is
impossible because a morph target IS the final vertex positions.

Pipeline (all offline, pure numpy on raw GLB bytes — no bpy, no Blender scene import):
  1. parse GLB, rest skeleton (node globals x IBM = skinning space)
  2. hand ROI: verts near the finger/hand bone segments, grown by edge rings
  2b. WELD the ROI graph: the scan is UV-chart soup, so hundreds of coincident
     duplicate verts sit on chart seams in separate graph components. Solve once per
     welded node and scatter the delta back to every duplicate (see weld_roi)
  3. per-bone geodesic fields: multi-source Dijkstra over the ROI subgraph
  4. weights: gaussian kernels on geodesic distance (top-4, renormalized, smoothed)
     — chain continuity is implicit: adjacent phalanges seed adjacent segments
  5. pose: parametric curl/spread/thumb-opposition per joint, world-axis rotations
     pivoted at each joint head, LBS with MY weights (the GLB's own weights are
     irrelevant — the shipped body is a rigid mitt and that is FINE)
  6. relax: stretch-gated Laplacian smoothing of the DELTA field until edge stretch
     bars pass — a converging diffusion, replacing the non-converging weight loop
  7. gates: numeric report (JSON) + OBJ dump for clay renders
  8. emit: splice morph-target accessors into the GLB (deltas added to base — exactly
     what Godot's gltf_document.cpp does: w = target + base) with extras.targetNames
     hand_<pose>_<l|r>; runtime driver = ariki-game src/Animation/HandMorphLayer.cs

Run under Blender's python for numpy (Blender is only the interpreter host — no bpy):
  "C:/Program Files/Blender Foundation/Blender 5.1/blender.exe" --background \
    --factory-startup --python tools/handshape_solve.py -- \
    --body <body.glb> --poses relaxed,fist,grip --out <out.glb> --workdir <dir>

  --selftest-bump   Spike A: splice one synthetic 3cm bump (no solve) to prove the
                    import+drive path end-to-end before trusting any solver output.
"""
import argparse
import heapq
import json
import math
import struct
import sys
from pathlib import Path

import numpy as np

DIGITS = ("thumb", "index", "middle", "ring", "pinky")
PHALANX = ("01", "02", "03")

# ───────────────────────── glTF container I/O ─────────────────────────

def read_glb(path):
    d = Path(path).read_bytes()
    ln = struct.unpack_from("<I", d, 8)[0]
    off = 12
    g = b = None
    while off < ln:
        clen, ct = struct.unpack_from("<II", d, off)
        off += 8
        if ct == 0x4E4F534A:
            g = json.loads(d[off:off + clen])
        elif ct == 0x004E4942:
            b = d[off:off + clen]
        off += clen
    return g, b


def load_accessor(g, bin_data, idx):
    """Numeric accessor -> float64 array (count,) or (count, ncomp)."""
    a = g["accessors"][idx]
    bv = g["bufferViews"][a["bufferView"]]
    nc = {"SCALAR": 1, "VEC2": 2, "VEC3": 3, "VEC4": 4, "MAT4": 16}[a["type"]]
    dtype = np.dtype({5120: "i1", 5121: "u1", 5122: "<i2", 5123: "<u2",
                      5125: "<u4", 5126: "<f4"}[a["componentType"]])
    size = dtype.itemsize * nc
    stride = bv.get("byteStride") or size
    off = bv.get("byteOffset", 0) + a.get("byteOffset", 0)
    raw = np.frombuffer(bin_data, dtype=np.uint8, count=a["count"] * stride,
                        offset=off).reshape(a["count"], stride)
    arr = raw[:, :size].copy().view(dtype).reshape(a["count"], nc).astype(np.float64)
    if a.get("normalized"):
        scale = {5120: 127.0, 5121: 255.0, 5122: 32767.0, 5123: 65535.0}
        arr /= scale.get(a["componentType"], 1.0)
    return arr


# ───────────────────────── skeleton ─────────────────────────

def quat_to_mat(q):
    x, y, z, w = q
    n = math.sqrt(x * x + y * y + z * z + w * w)
    x, y, z, w = x / n, y / n, z / n, w / n
    return np.array([
        [1 - 2 * (y * y + z * z), 2 * (x * y - z * w), 2 * (x * z + y * w)],
        [2 * (x * y + z * w), 1 - 2 * (x * x + z * z), 2 * (y * z - x * w)],
        [2 * (x * z - y * w), 2 * (y * z + x * w), 1 - 2 * (x * x + y * y)],
    ])


def node_local(nd):
    M = np.eye(4)
    M[:3, :3] = quat_to_mat(nd.get("rotation", [0, 0, 0, 1])) \
        * np.array(nd.get("scale", [1, 1, 1]))[None, :]
    M[:3, 3] = nd.get("translation", [0, 0, 0])
    return M


class Skeleton:
    def __init__(self, g, bin_data):
        skin = g["skins"][0]
        self.joints = skin["joints"]
        self.names = [g["nodes"][j].get("name", "") for j in self.joints]
        self.idx = {n: i for i, n in enumerate(self.names)}
        loc = [node_local(g["nodes"][j]) for j in self.joints]
        node_to_joint = {j: i for i, j in enumerate(self.joints)}
        parent_of = {}
        for i, j in enumerate(self.joints):
            for c in g["nodes"][j].get("children", []):
                if c in node_to_joint:
                    parent_of[node_to_joint[c]] = i
        self.parent = [parent_of.get(i, -1) for i in range(len(self.joints))]
        G = [None] * len(self.joints)

        def glob(i):
            if G[i] is None:
                G[i] = loc[i] if self.parent[i] < 0 else glob(self.parent[i]) @ loc[i]
            return G[i]

        self.rest = [glob(i) for i in range(len(self.joints))]
        ibm = load_accessor(g, bin_data, skin["inverseBindMatrices"]).reshape(-1, 4, 4)
        self.ibm = np.array([m.T for m in ibm])   # glTF is column-major

    def head(self, i):
        return self.rest[i][:3, 3].copy()

    def child_head(self, i, fallback_dir=None):
        for j, par in enumerate(self.parent):
            if par == i:
                return self.head(j), True
        h = self.head(i)
        if fallback_dir is not None:
            return h + fallback_dir, False
        return h + np.array([0.0, 0.0, 0.02]), False


def finger_bones(skel, side):
    sfx = "_l" if side == "l" else "_r"
    hand = "hand" + sfx
    bones = [f"{d}_{p}{sfx}" for d in DIGITS for p in PHALANX if f"{d}_{p}{sfx}" in skel.idx]
    return hand, bones


def refit_fingers(skel, P, side, lat=0.04, pct=99.5, verbose=True):
    """Shrink each digit's joint chain along its own axis until it spans the ACTUAL
    flesh, and rebuild that digit's inverse bind matrices to match.

    Measured on Ariki_Female_QuatSkin_LowPoly_40 (2026-08-18): Lena's mesh hand ends
    14.0cm from the wrist joint, but the Quaternius finger chain runs out to 19.8cm —
    the `_02` row sits 1.9-2.8cm outside the mesh and the `_03` row floats 4.7-5.9cm
    beyond the fingertips, in empty space. Curling about pivots outside the surface
    translates the whole mitt end sideways instead of bending anything, and it starves
    the weight solver: no vertex is within the 4mm seed radius of `_02`/`_03`, so the
    "200 euclidean-nearest" fallback fires for 10 of 15 bones per hand.

    The shipped skeleton is NOT touched — every clip pins all 65 bone positions and this
    body has to keep meeting it. This is solver-local scaffolding: a morph target is just
    final vertex positions, so the pose only has to be defined by pivots that lie inside
    the flesh they bend. Uniform scale about the wrist per digit (so the knuckle row
    moves inward too, not just the tips), and ibm := inv(rest) so the refit chain still
    reproduces the rest mesh exactly under LBS.
    """
    sfx = "_l" if side == "l" else "_r"
    hh = skel.head(skel.idx["hand" + sfx])
    rel_all = P - hh
    near = np.linalg.norm(rel_all, axis=1) < 0.30
    out = {}
    for d in DIGITS:
        chain = [f"{d}_{ph}{sfx}" for ph in PHALANX if f"{d}_{ph}{sfx}" in skel.idx]
        if not chain:
            continue
        tip = skel.head(skel.idx[chain[-1]])
        u = tip - hh
        L = float(np.linalg.norm(u))
        if L < 1e-6:
            continue
        u /= L
        t = rel_all @ u
        lat_d = np.linalg.norm(rel_all - t[:, None] * u[None, :], axis=1)
        sel = near & (t > 0.3 * L) & (lat_d < lat)
        if sel.sum() < 20:
            continue
        reach = float(np.percentile(t[sel], pct))
        k = float(np.clip(reach / L, 0.35, 1.0))
        # every joint at or below <digit>_01 scales about the wrist
        root = skel.idx[chain[0]]
        fam = [root]
        for i in range(len(skel.parent)):
            j, hops = i, 0
            while skel.parent[j] >= 0 and hops < 8:
                j = skel.parent[j]
                hops += 1
                if j == root:
                    fam.append(i)
                    break
        for i in fam:
            R = skel.rest[i].copy()
            R[:3, 3] = hh + k * (R[:3, 3] - hh)
            skel.rest[i] = R
            skel.ibm[i] = np.linalg.inv(R)
        out[d] = (k, reach, L)
    if verbose and out:
        print("[refit] " + "  ".join(
            f"{d}:{v[0]:.2f}({v[2] * 100:.1f}->{v[1] * 100:.1f}cm)" for d, v in out.items()))
    return out


def bone_segments(skel, hand, bones):
    """World (head, tail) per bone; leaf tails extrapolate one phalanx outward."""
    sfx = "_l" if hand.endswith("_l") else "_r"
    segs = {}
    hand_head = skel.head(skel.idx[hand])
    for b in bones:
        d, p, _ = b.rsplit("_", 2)
        i = skel.idx[b]
        h = skel.head(i)
        t, had_child = skel.child_head(i)
        if not had_child:
            pn = f"{d}_{int(p) - 1:02d}{sfx}"
            if pn in skel.idx:
                ph = skel.head(skel.idx[pn])
                t = h + (h - ph)
            else:
                t = h + (hand_head - h) * -0.3
        segs[b] = (h, t)
    first = bones[0] if bones else hand
    segs[hand] = (hand_head, skel.head(skel.idx[first]))
    return segs


def point_segment_dist(P, a, b):
    ab = b - a
    ab2 = ab @ ab
    if ab2 < 1e-12:
        return np.linalg.norm(P - a, axis=1)
    t = np.clip((P - a) @ ab / ab2, 0.0, 1.0)
    return np.linalg.norm(P - (a[None, :] + t[:, None] * ab[None, :]), axis=1)


# ───────────────────────── ROI + graph ─────────────────────────

def hand_joints(skel):
    """Wrist + every finger joint, both hands, in skinning space."""
    J = []
    for side in ("l", "r"):
        J.append(skel.head(skel.idx[f"hand_{side}"]))
        for d in DIGITS:
            for ph in PHALANX:
                nm = f"{d}_{ph}_{side}"
                if nm in skel.idx:
                    J.append(skel.head(skel.idx[nm]))
    return np.asarray(J)


def build_roi(P, F, skel, grow_rings=2, near=0.040):
    """ROI = flesh around the hand skeleton, as a union of spheres about the wrist and
    every (refit) finger joint — NOT a tube about the bone segments.

    A segment tube leaves this scan's ROI riddled with interior holes: any patch that
    happens to sit farther than `near` from a bone axis drops out, so the "ROI rim"
    (verts sharing a face with a non-ROI vert) becomes a fractal of interior chart holes
    rather than a wrist ring. Locking that as a boundary froze verts in the middle of the
    curling flesh — measured 48x edge stretch and 805 torn edges on fist_r. A union of
    joint spheres is blob-like, so its intersection with the surface is one clean band on
    the forearm: measured 766 rim verts, 95% of them at radius 4.3-6cm from the wrist,
    none among the fingers. That is a rim worth locking, and it is what makes the delta
    field fade to zero at the arm instead of cracking there.
    """
    mask = np.min(np.linalg.norm(P[:, None, :] - hand_joints(skel)[None, :, :], axis=2),
                  axis=1) < near
    for _ in range(grow_rings):
        fmask = mask[F].any(axis=1)
        cand = np.unique(F[fmask].reshape(-1))
        new = cand[~mask[cand]]
        if len(new) == 0:
            break
        mask[new] = True
    idx = np.where(mask)[0]
    remap = -np.ones(len(P), dtype=np.int64)
    remap[idx] = np.arange(len(idx))
    rf = remap[F]
    roi_faces = rf[(rf >= 0).all(axis=1)]
    # True ROI rim = an ROI vert sharing a face with a vert OUTSIDE the ROI. This has to
    # be measured against the FULL mesh: "low degree in the ROI subgraph" is a different
    # predicate and, on chart soup, mostly selects INTERIOR chart-corner dangles (verts
    # in a single triangle). Locking those froze their delta at zero while their
    # neighbours curled 4cm away — the entire remaining >5x needle population on the
    # right hand was exactly two such verts.
    part = F[(~mask[F]).any(axis=1) & mask[F].any(axis=1)]
    rim = np.zeros(len(P), dtype=bool)
    if len(part):
        rim[part.reshape(-1)] = True
    rim &= mask
    return idx, remap, roi_faces, rim[idx]


def csr_from_edges(n, e):
    """Undirected adjacency in CSR form (nbr, start)."""
    deg = np.zeros(n, dtype=np.int64)
    np.add.at(deg, e[:, 0], 1)
    np.add.at(deg, e[:, 1], 1)
    start = np.zeros(n + 1, dtype=np.int64)
    np.cumsum(deg, out=start[1:])
    nbr = np.zeros(len(e) * 2, dtype=np.int64)
    fill = start[:-1].copy()
    nbr[fill[e[:, 0]]] = e[:, 1]
    fill[e[:, 0]] += 1
    nbr[fill[e[:, 1]]] = e[:, 0]
    fill[e[:, 1]] += 1
    return nbr, start


def roi_graph(P_roi, roi_faces, extra=None):
    e = np.concatenate([roi_faces[:, [0, 1]], roi_faces[:, [1, 2]], roi_faces[:, [2, 0]]])
    if extra is not None and len(extra):
        e = np.concatenate([e, extra])
    e = np.unique(np.sort(e, axis=1), axis=0)
    nbr, start = csr_from_edges(len(P_roi), e)
    return nbr, start, e


def components(n, e):
    par = np.arange(n)

    def find(a):
        while par[a] != a:
            par[a] = par[par[a]]
            a = par[a]
        return a

    for a, b in e:
        ra, rb = find(int(a)), find(int(b))
        if ra != rb:
            par[ra] = rb
    root = np.array([find(i) for i in range(n)])
    _, comp = np.unique(root, return_inverse=True)
    return comp.astype(np.int64).reshape(-1)


def stitch_components(P_w, comp, tol=0.012, verbose=True):
    """Bridge SEPARATE surface sheets that lie within `tol` of each other.

    Welding fixes coincident duplicates, but this scan is worse than duplicated: the hand
    is built from overlapping sheets with real 1-10mm gaps between them. Measured
    2026-08-18 on the refit hand ROI: 6 components, right hand 2111/1558/759 verts with
    4.0mm and 1.1mm gaps, left hand 2082/2028 with an 8.9mm gap. Dijkstra cannot cross a
    gap, so each bone's geodesic field covers whichever sheet its seeds landed on and
    goes blind on the others; diffusion likewise cannot move a delta between sheets. That
    is why one hand's fingertips curled 6.7cm while the other's moved 0.2cm from the same
    parameters — nothing to do with left/right handedness, just which sheet got seeded.

    Only CROSS-component pairs are joined, which is what makes this safe: a stitch can
    never short-circuit two points that already have a path through the surface, so it
    cannot fake a shortcut inside a sheet. It only restores coupling the scan lost.
    """
    n = len(P_w)
    out = []
    for lo in range(0, n, 512):
        hi = min(n, lo + 512)
        d = np.linalg.norm(P_w[lo:hi, None, :] - P_w[None, :, :], axis=2)
        d[comp[lo:hi, None] == comp[None, :]] = np.inf
        j = d.argmin(axis=1)
        dm = d[np.arange(hi - lo), j]
        keep = dm < tol
        if keep.any():
            out.append(np.stack([np.arange(lo, hi)[keep], j[keep]], axis=1))
    if not out:
        if verbose:
            print(f"[stitch] no cross-sheet pair within {tol * 1000:.0f}mm")
        return np.zeros((0, 2), dtype=np.int64)
    ex = np.unique(np.sort(np.concatenate(out), axis=1), axis=0)
    if verbose:
        L = np.linalg.norm(P_w[ex[:, 0]] - P_w[ex[:, 1]], axis=1)
        print(f"[stitch] {len(ex)} cross-sheet edges, len mm "
              f"med {np.median(L) * 1000:.2f} max {L.max() * 1000:.2f}")
    return ex


WELD_TOL = 1e-5     # 10um: 200x below the 1.96mm median ROI edge, so no real edge collapses


def weld_roi(P, tol=WELD_TOL):
    """Merge coincident ROI vertices into ONE graph node. Returns (wid, P_w, group_size).

    The Tripo scan is UV-chart soup. Measured on Ariki_Female_QuatSkin_LowPoly_40
    (2026-08-18): 356 groups / 732 of the 5,548 ROI verts are duplicate positions
    lying on chart seams, and the raw ROI graph therefore has 16 disconnected
    components instead of 2 hands. Three separate failures fall out of that:

      * duplicates are solved independently and get DIFFERENT deltas (measured up to
        2.55cm apart on fist_r) — the seam physically cracks open;
      * seam verts have a truncated one-ring, so `deg < 3` classes them as ROI
        BOUNDARY and locks their delta to zero — relaxation was forbidden from
        touching the exact verts that were tearing;
      * geodesic distance detours around every seam, fragmenting the weight field.

    Welding is done at the GRAPH level only — the mesh is never rewritten and morph
    deltas stay indexed by original vertex id (required, or the splice breaks). Since
    every member of a weld group then receives an identical delta, coincident pairs
    keep their length exactly and the whole failure class dies by construction instead
    of by tuning. Grid-bucket + 27-cell neighbour union-find rather than plain
    round-and-unique, so a pair straddling a cell boundary still merges.
    """
    cell = 2.0 * tol
    keys = np.floor(P / cell).astype(np.int64)
    buckets = {}
    for i, k in enumerate(map(tuple, keys)):
        buckets.setdefault(k, []).append(i)
    par = np.arange(len(P))

    def find(a):
        while par[a] != a:
            par[a] = par[par[a]]
            a = par[a]
        return a

    t2 = tol * tol
    offs = [(dx, dy, dz) for dx in (-1, 0, 1) for dy in (-1, 0, 1) for dz in (-1, 0, 1)]
    for k, ids in buckets.items():
        cand = []
        for dx, dy, dz in offs:
            b = buckets.get((k[0] + dx, k[1] + dy, k[2] + dz))
            if b:
                cand.extend(b)
        cand = np.asarray(cand)
        for i in ids:
            d2 = ((P[cand] - P[i]) ** 2).sum(axis=1)
            for j in cand[d2 <= t2]:
                ri, rj = find(i), find(int(j))
                if ri != rj:
                    par[ri] = rj
    root = np.array([find(i) for i in range(len(P))])
    _, wid = np.unique(root, return_inverse=True)
    wid = wid.astype(np.int64).reshape(-1)
    n_w = int(wid.max()) + 1
    grp = np.bincount(wid, minlength=n_w).astype(np.float64)
    P_w = np.stack([np.bincount(wid, weights=P[:, c], minlength=n_w) / grp
                    for c in range(3)], axis=1)
    return wid, P_w, grp


def weld_faces(roi_faces, wid):
    """Re-index faces onto welded nodes, dropping the ones that collapse. Winding is
    preserved (never sort face index triples)."""
    fw = wid[roi_faces]
    keep = ((fw[:, 0] != fw[:, 1]) & (fw[:, 1] != fw[:, 2]) & (fw[:, 2] != fw[:, 0]))
    return fw[keep]


def dijkstra_multi(P_roi, nbr, start, seeds):
    n = len(P_roi)
    INF = float("inf")
    dist = [INF] * n
    heap = [(0.0, int(s)) for s in seeds]
    for _, s in heap:
        dist[s] = 0.0
    heapq.heapify(heap)
    while heap:
        d, v = heapq.heappop(heap)
        if d > dist[v]:
            continue
        pv = P_roi[v]
        for k in range(start[v], start[v + 1]):
            u = int(nbr[k])
            du = d + float(np.linalg.norm(pv - P_roi[u]))
            if du < dist[u]:
                dist[u] = du
                heapq.heappush(heap, (du, u))
    return np.asarray(dist)


# ───────────────────────── weights ─────────────────────────

def digit_chain(skel, side, digit):
    """Polyline [wrist, ph01, ph02, ph03, tip] for one digit, plus per-segment arc."""
    sfx = "_l" if side == "l" else "_r"
    pts = [skel.head(skel.idx[f"hand{sfx}"])]
    last = None
    for ph in PHALANX:
        nm = f"{digit}_{ph}{sfx}"
        if nm in skel.idx:
            pts.append(skel.head(skel.idx[nm]))
            last = skel.idx[nm]
    if last is not None:
        tail, had = skel.child_head(last)
        if not had and len(pts) >= 3:
            tail = pts[-1] + (pts[-1] - pts[-2])
        pts.append(tail)
    return np.asarray(pts)


def chain_project(P, pts):
    """Per-vertex (arc coordinate along the polyline, lateral distance to it)."""
    seg_len = np.linalg.norm(np.diff(pts, axis=0), axis=1)
    cum = np.concatenate([[0.0], np.cumsum(seg_len)])
    best_d = np.full(len(P), np.inf)
    best_s = np.zeros(len(P))
    for k in range(len(pts) - 1):
        a, b = pts[k], pts[k + 1]
        ab = b - a
        ab2 = float(ab @ ab)
        if ab2 < 1e-12:
            continue
        t = np.clip((P - a) @ ab / ab2, 0.0, 1.0)
        proj = a[None, :] + t[:, None] * ab[None, :]
        d = np.linalg.norm(P - proj, axis=1)
        upd = d < best_d
        best_d[upd] = d[upd]
        best_s[upd] = cum[k] + t[upd] * seg_len[k]
    return best_s, best_d, cum


def solve_weights(skel, P_roi, nbr, start, edges, side, sig_lat=0.004, support=0.05):
    """Weights as a partition of unity along each digit's CHAIN ARC, not as isotropic
    kernels around bone segments.

    The kernel version could not produce an owned vertex. Its width (sig_f = 2.5 x median
    edge = 8mm) is wider than the surface's distance to the NEIGHBOURING phalanx (~1-1.5cm
    on a hand this size), so every vertex came out a 3-way blend: measured 2026-08-18,
    ZERO of 17,480 ROI verts had more than 0.85 weight on any phalanx and only 12 had
    more than 0.70. Blending three phalanges that rotate by cumulatively different
    amounts averages the curl away — which is exactly what the numbers showed: the raw
    pose reached only 2.6cm of tip motion on one hand, and relaxation then flattened it
    to 0.3cm while happily reporting every stretch bar as passed.

    Along-chain arc position is the natural coordinate for a finger: tent functions
    centred on each segment's midpoint give 1.0 mid-phalanx, a clean 50/50 at each joint,
    and a smooth handover to the rigid hand bone behind the knuckle. Lateral distance to
    the chain then picks WHICH digit owns the vertex. Two coordinates, no tuning race.
    """
    hand, bones = finger_bones(skel, side)
    all_bones = [hand] + bones
    col = {b: i for i, b in enumerate(all_bones)}
    sfx = "_l" if side == "l" else "_r"
    W = np.zeros((len(P_roi), len(all_bones)))
    lat_min = np.full(len(P_roi), np.inf)
    chains, arcs = {}, {}
    for d in DIGITS:
        pts = digit_chain(skel, side, d)
        if len(pts) < 3:
            continue
        s_arc, lat, cum = chain_project(P_roi, pts)
        chains[d] = (pts, s_arc, lat, cum)
        lat_min = np.minimum(lat_min, lat)
    for d in DIGITS:
        if d not in chains:
            continue
        pts, s_arc, lat, cum = chains[d]
        # segment midpoints in arc space: index 0 is the palm (wrist->knuckle) segment,
        # then one per phalanx present
        mid = 0.5 * (cum[:-1] + cum[1:])
        names = [hand] + [f"{d}_{ph}{sfx}" for ph in PHALANX if f"{d}_{ph}{sfx}" in skel.idx]
        names = names[:len(mid)]
        # 1D partition of unity along the arc: a narrow handover ramp centred on each
        # JOINT, not a tent spanning midpoint-to-midpoint. A midpoint tent is 1.0 only at
        # the exact midpoint and decays linearly across the whole phalanx, so the typical
        # vertex still ends up a 2-way blend (measured p50 ownership 0.55). Confining the
        # blend to +/-25% of a phalanx around each joint gives 1.0 through the middle of
        # each segment and a clean 0.5 exactly at the joint, which is what a hand rig
        # looks like and what lets a curl actually accumulate.
        seg = np.diff(cum)
        u = [np.ones(len(P_roi))]
        for k in range(1, len(mid)):
            w = 0.25 * min(seg[k - 1], seg[k])
            u.append(np.clip((s_arc - (cum[k] - w)) / max(2 * w, 1e-9), 0.0, 1.0))
        u.append(np.zeros(len(P_roi)))
        T = np.stack([u[k] * (1.0 - u[k + 1]) for k in range(len(mid))], axis=1)
        T /= np.maximum(T.sum(axis=1, keepdims=True), 1e-12)
        # Digit ownership from lateral distance RELATIVE to the nearest chain, not
        # absolute. The whole surface sits ~1cm off its own chain (that is just the
        # finger's radius), so an absolute kernel suppresses every digit equally and
        # renormalization hands the vertex back as a 5-way blend — measured p50 = 0.47
        # ownership, which averages the curl away exactly as the isotropic version did.
        # Against the nearest chain, a vertex 4mm farther from digit B than from digit A
        # scores 0.37 on B, 8mm farther scores 0.02, and A itself gets 1.0. Only genuine
        # near-ties (the inter-digit webs) blend, which is what should blend.
        aff = np.exp(-((lat - lat_min) / sig_lat) ** 2)
        for k, nm in enumerate(names):
            if nm in col:
                W[:, col[nm]] += aff * T[:, k]
    # a floor on the hand bone so verts no digit claims stay rigid instead of
    # normalizing a zero row into noise
    W[:, 0] += 1e-3
    # A sheet Dijkstra never reaches from this hand is detached geometry (see
    # stitch_components): give it to the hand bone rather than letting euclidean lateral
    # distance hand it a phalanx it has no surface path to.
    seeds = np.where(lat_min < 0.012)[0]
    if len(seeds) < 20:
        seeds = np.argsort(lat_min)[:200]
    geo = dijkstra_multi(P_roi, nbr, start, seeds.tolist())
    far = (~np.isfinite(geo)) | (lat_min > support)
    if far.any():
        W[far] = 0.0
        W[far, 0] = 1.0
    W /= np.maximum(W.sum(axis=1, keepdims=True), 1e-12)
    # one light smoothing pass: kills per-facet noise without un-owning anything
    S = np.zeros_like(W)
    cnt = (start[1:] - start[:-1]).clip(1)
    for c in range(W.shape[1]):
        S[:, c] = np.add.reduceat(W[nbr, c], start[:-1]) / cnt
    W = 0.85 * W + 0.15 * S
    W /= np.maximum(W.sum(axis=1, keepdims=True), 1e-12)
    return all_bones, W


# ───────────────────────── posing ─────────────────────────

def palm_normal(skel, side):
    sfx = "_l" if side == "l" else "_r"
    h = skel.head(skel.idx["hand" + sfx])
    mid = skel.head(skel.idx["middle_01" + sfx])
    ix = skel.head(skel.idx["index_01" + sfx])
    pk = skel.head(skel.idx["pinky_01" + sfx])
    n = np.cross(mid - h, pk - ix)
    n /= np.linalg.norm(n)
    # NO artificial sign flip: the cross product already mirrors correctly between
    # hands (left/right palms face opposite ways in the bind pose). A z-flip heuristic
    # here breaks exactly one side and curls its fingers backward — measured the hard
    # way (left raw stretch 2000x vs right 5x). Per-joint curl DIRECTION is instead
    # verified against this normal by curl_axis's self-check.
    return n


def axis_angle(axis, theta):
    x, y, z = axis / max(np.linalg.norm(axis), 1e-12)
    c, s = math.cos(theta), math.sin(theta)
    C = 1 - c
    return np.array([
        [c + x * x * C, x * y * C - z * s, x * z * C + y * s, 0],
        [y * x * C + z * s, c + y * y * C, y * z * C - x * s, 0],
        [z * x * C - y * s, z * y * C + x * s, c + z * z * C, 0],
        [0, 0, 0, 1],
    ])


def curl_axis(skel, side, n, digit, ph, sfx):
    """Axis = normalize(u x n) at the joint; sign-checked so +theta curls INTO the palm."""
    name = f"{digit}_{ph}{sfx}"
    i = skel.idx[name]
    h = skel.head(i)
    child, _ = skel.child_head(i, fallback_dir=None)
    u = child - h
    u /= max(np.linalg.norm(u), 1e-9)
    axis = np.cross(u, n)
    an = np.linalg.norm(axis)
    if an < 1e-6:
        return None, None
    axis /= an
    # self-check: rotating a probe point on the finger by +0.5 rad must move it
    # toward the palm (component along -n grows)
    moved = axis_angle(axis, 0.5)[:3, :3] @ u
    if (moved @ (-n)) <= 0:
        axis = -axis
    return h, axis


def pose_globals(skel, side, params):
    """{bone: 4x4 posed global}. Each finger joint rotates ONLY its own curl (plus
    spread at the MCP); ancestors' rotations enter through the parent recursion:
    G_pose(j) = G_pose(p) . G_rest(p)^-1 . Rot_j . G_rest(j)."""
    sfx = "_l" if side == "l" else "_r"
    n = palm_normal(skel, side)
    hand, bones = finger_bones(skel, side)
    own = {}      # bone -> list[(axis, theta, pivot)]
    for d in DIGITS:
        p = params["fingers"].get(d)
        if p is None:
            continue
        for ph in PHALANX:
            name = f"{d}_{ph}{sfx}"
            if name not in skel.idx:
                continue
            h, axis = curl_axis(skel, side, n, d, ph, sfx)
            if axis is None:
                continue
            own.setdefault(name, [])
            if ph == "01" and p.get("spread", 0.0):
                own[name].append((n, p["spread"], h))
            if axis is not None:
                own[name].append((axis, p[f"curl_{ph}"], h))
    tp = params.get("thumb", {})
    for ph in PHALANX:
        name = f"thumb_{ph}{sfx}"
        if name not in skel.idx:
            continue
        h, axis = curl_axis(skel, side, n, "thumb", ph, sfx)
        if axis is not None:
            own.setdefault(name, []).append((axis, tp.get(f"curl_{ph}", 0.0), h))
    if f"thumb_01{sfx}" in own:
        palm_c = skel.head(skel.idx[f"hand{sfx}"])
        idx_mcp = skel.head(skel.idx[f"index_01{sfx}"])
        op = idx_mcp - palm_c
        op /= max(np.linalg.norm(op), 1e-9)
        own[f"thumb_01{sfx}"].insert(0, (op, tp.get("opposition", 0.0), palm_c))

    Gp = {}

    def glob(i):
        name = skel.names[i]
        if name in Gp:
            return Gp[name]
        G = skel.rest[i].copy()
        for axis, theta, pivot in own.get(name, []):
            T1 = np.eye(4); T1[:3, 3] = pivot
            T2 = np.eye(4); T2[:3, 3] = -pivot
            G = T1 @ axis_angle(axis, theta) @ T2 @ G
        par = skel.parent[i]
        if par >= 0:
            G = glob(par) @ np.linalg.inv(skel.rest[par]) @ G
        Gp[name] = G
        return G

    for b in bones:
        glob(skel.idx[b])
    Gp[hand] = skel.rest[skel.idx[hand]].copy()
    return Gp


def lbs(P_roi, skel, all_bones, W, Gp):
    out = np.zeros_like(P_roi)
    Ph = np.concatenate([P_roi, np.ones((len(P_roi), 1))], axis=1)
    for bi, bn in enumerate(all_bones):
        w = W[:, bi]
        sel = np.where(w > 1e-9)[0]
        if len(sel) == 0:
            continue
        G = Gp.get(bn, skel.rest[skel.idx[bn]])
        D = G @ skel.ibm[skel.idx[bn]]
        out[sel] += w[sel, None] * (Ph[sel] @ D.T)[:, :3]
    # verts with no weight on this hand's bones (other hand, wrist edge of the ROI)
    # must KEEP their rest position — a zero fallback flings them to the origin,
    # which reads as 1500x edge stretch and poisons the relaxation. A weight row
    # that doesn't sum to ~1 (degenerate isolated verts) is garbage too — same rule.
    unw = W.sum(axis=1) <= 0.5
    out[unw] = P_roi[unw]
    return out


# ───────────────────────── relaxation ─────────────────────────

def vertex_normals(P, F):
    """Area-weighted vertex normals, vectorized, over the FULL mesh (so ROI-rim verts
    get their complete one-ring — a ROI-local pass would shade the rim differently
    than the imported rest normals and leave a discontinuity ring at the wrist)."""
    v0, v1, v2 = P[F[:, 0]], P[F[:, 1]], P[F[:, 2]]
    fn = np.cross(v1 - v0, v2 - v0)          # area-weighted (unnormalized cross)
    n = np.zeros_like(P)
    for col in range(3):
        np.add.at(n, F[:, col], fn)
    lens = np.linalg.norm(n, axis=1, keepdims=True)
    lens[lens < 1e-12] = 1.0
    return n / lens


NEEDLE_MM = 1.0     # absolute growth above which an over-stretched edge is a real spike


def stretch_stats(P, Q, e):
    """Edge stretch with NO rest-length floor, plus the absolute over-stretch a bare
    ratio hides.

    The old 1mm floor was justified as "a 0.1mm edge doubling is invisible", which is
    true — but it also silently excluded a real population: sub-mm edges on chart seams
    that grew to 3-4cm in fist/grip, i.e. hairline needles sticking out of the hand,
    sub-pixel in the screenshots that "passed" and entirely absent from the numbers. So
    ratios are reported over EVERY edge, and a tear is judged by ratio AND absolute
    growth together, which is scale-honest for both a 2mm web bridge at 5x and a 0.1mm
    decimation sliver at 5x. `max`/`p999`/`n_gt*` stay scoped to >=1mm edges so the
    table remains comparable to the pre-weld baseline."""
    lr = np.linalg.norm(P[e[:, 0]] - P[e[:, 1]], axis=1)
    lq = np.linalg.norm(Q[e[:, 0]] - Q[e[:, 1]], axis=1)
    r = lq / np.maximum(lr, 1e-9)
    grow_mm = (lq - lr) * 1000.0
    slv = lr < 0.001
    big = ~slv
    return {"max": float(r[big].max()) if big.any() else 1.0,
            "p999": float(np.percentile(r[big], 99.9)) if big.any() else 1.0,
            "n_gt2": int((big & (r > 2)).sum()), "n_gt5": int((big & (r > 5)).sum()),
            # the sliver population the floor used to hide
            "slv_n": int(slv.sum()),
            "slv_gt5x": int((slv & (r > 5)).sum()),
            "slv_grow_gt1mm": int((slv & (grow_mm > NEEDLE_MM)).sum()),
            "slv_max_grow_mm": round(float(grow_mm[slv].max()), 2) if slv.any() else 0.0,
            # scale-honest tear count over ALL edges: >5x AND >1mm of real growth
            "needles": int(((r > 5) & (grow_mm > NEEDLE_MM)).sum()),
            "max_grow_mm": round(float(grow_mm.max()), 2)}


def relax(P_roi, delta, e, lr, pin, seam=None, tau=1.35, iters=1500, omega=0.6):
    """Strain-only edge projection: move ONLY the endpoints of edges that exceed `tau`.

    This replaces the gated-Laplacian-diffusion relaxation, which was structurally unable
    to do the job. Diffusing the DELTA field has a null space — constants — and the
    edge-stretch gate is blind to every member of it, because a rigid translation stretches
    no edge and neither does a collapse to zero. So the diffusion always found one of those
    two exits, and reported perfect bars on the way out. Measured 2026-08-18 on
    Ariki_Female_QuatSkin_LowPoly_40: the left hand's delta decayed to 0.3cm of fingertip
    travel (max 1.93x, p99.9 1.60x, zero torn edges — a flawless report for a morph that
    does nothing), while the right hand converged to a near-constant 6.5cm delta at EVERY
    arc position from wrist to fingertip, i.e. the whole hand translated 6.5cm sideways
    with its shape intact (max 2.43x, p99.9 1.46x, also "passing"). That is the real
    reason five weight-rebake iterations and every gate before this one signed off on
    hands that do not make a fist.

    Projection has no such exit: a conforming edge contributes no correction, so regions
    that are not over-stretched are left exactly as the pose put them, and the pose can
    only be modified where it actually tears. Jacobi-style (accumulate, average by
    incidence count, under-relax by `omega`) so dense web clusters cannot oscillate.

    `seam` = (inside_vert, fixed_outside_position, rest_len) constrains ROI-border edges
    against the un-morphed body. Those edges are NOT in `e` — they leave the ROI — so
    without them nothing measures or limits the crack at the wrist, and the un-anchored
    hand is free to walk away from the arm. Same 70-140mm drift as above, from the other
    side of the same blind spot.
    """
    d = delta.copy()
    n = len(P_roi)
    free = (~pin).astype(np.float64)[:, None]
    hist = []
    for _ in range(iters):
        Q = P_roi + d
        lq = np.linalg.norm(Q[e[:, 0]] - Q[e[:, 1]], axis=1)
        viol = lq > tau * lr
        corr = np.zeros_like(d)
        cnt = np.zeros(n)
        if viol.any():
            a, b = e[viol, 0], e[viol, 1]
            dv = Q[b] - Q[a]
            L = np.linalg.norm(dv, axis=1)
            ex = (L - tau * lr[viol])[:, None] * (dv / np.maximum(L, 1e-12)[:, None]) * 0.5
            np.add.at(corr, a, ex)
            np.add.at(cnt, a, 1)
            np.add.at(corr, b, -ex)
            np.add.at(cnt, b, 1)
        n_seam = 0
        if seam is not None:
            si, sp, sl = seam
            dvs = Q[si] - sp
            Ls = np.linalg.norm(dvs, axis=1)
            vs = Ls > tau * sl
            n_seam = int(vs.sum())
            if n_seam:
                exs = (Ls[vs] - tau * sl[vs])[:, None] * (dvs[vs] / np.maximum(Ls[vs], 1e-12)[:, None])
                np.add.at(corr, si[vs], -exs)
                np.add.at(cnt, si[vs], 1)
        hist.append((int(viol.sum()), n_seam))
        if not viol.any() and n_seam == 0:
            break
        d += omega * free * corr / np.maximum(cnt, 1)[:, None]
    return d, hist


# ───────────────────────── pose parameters ─────────────────────────

DEFAULT_PARAMS = {
    "flat":    {"fingers": {d: {"curl_01": 0.02, "curl_02": 0.02, "curl_03": 0.01,
                                "spread": 0.0} for d in DIGITS},
                "thumb": {"curl_01": 0.05, "curl_02": 0.05, "curl_03": 0.02,
                          "opposition": 0.0}},
    "relaxed": {"fingers": {d: {"curl_01": 0.15, "curl_02": 0.25, "curl_03": 0.15,
                                "spread": 0.05} for d in DIGITS},
                "thumb": {"curl_01": 0.15, "curl_02": 0.15, "curl_03": 0.1,
                          "opposition": 0.25}},
    "fist":    {"fingers": {d: {"curl_01": 1.15, "curl_02": 1.05, "curl_03": 0.75,
                                "spread": -0.05} for d in DIGITS},
                "thumb": {"curl_01": 0.9, "curl_02": 0.9, "curl_03": 0.5,
                          "opposition": 0.8}},
    "grip":    {"fingers": {d: {"curl_01": 0.75, "curl_02": 0.85, "curl_03": 0.7,
                                "spread": 0.1} for d in DIGITS},
                "thumb": {"curl_01": 1.0, "curl_02": 0.9, "curl_03": 0.6,
                          "opposition": 1.2}},
}


# ───────────────────────── morph splice emit ─────────────────────────

def emit_morph_glb(g, bin_data, shapes, out_path):
    """Splice {name: {"POSITION": (N,3) delta, "NORMAL": (N,3) delta|None}} into a NEW
    GLB. Godot's gltf importer computes w = target + base for BOTH attributes
    (gltf_document.cpp), so deltas are written verbatim; names ride in mesh
    extras.targetNames. Existing bufferViews are untouched — new data is appended to
    the BIN chunk and buffers[0].byteLength grows. Never writes the source body."""
    mesh = g["meshes"][0]
    prim = mesh["primitives"][0]
    vert_count = g["accessors"][prim["attributes"]["POSITION"]]["count"]
    new_bin = bytearray()
    targets, names = [], []
    for name, attrs in shapes.items():
        target_entry = {}
        for attr in ("POSITION", "NORMAL"):
            delta = attrs.get(attr)
            if delta is None:
                continue
            assert delta.shape == (vert_count, 3), f"{name}.{attr}: {delta.shape} vs {vert_count}"
            pad = (4 - len(new_bin) % 4) % 4
            new_bin += b"\x00" * pad
            f32 = delta.astype("<f4")
            bv_idx = len(g["bufferViews"])
            g["bufferViews"].append({"buffer": 0, "byteOffset": len(bin_data) + len(new_bin),
                                     "byteLength": f32.nbytes})
            acc_idx = len(g["accessors"])
            g["accessors"].append({"bufferView": bv_idx, "componentType": 5126,
                                   "count": vert_count, "type": "VEC3"})
            new_bin += f32.tobytes()
            target_entry[attr] = acc_idx
        targets.append(target_entry)
        names.append(name)
    new_bin += b"\x00" * ((4 - len(new_bin) % 4) % 4)
    g["buffers"][0]["byteLength"] = len(bin_data) + len(new_bin)
    # G7 bounds self-check: every bufferView must fit the declared buffer
    total = g["buffers"][0]["byteLength"]
    for bv in g["bufferViews"]:
        assert bv["byteOffset"] + bv["byteLength"] <= total, \
            f"bufferView out of bounds: {bv} vs {total}"

    prim["targets"] = targets
    mesh["weights"] = [0.0] * len(targets)
    mesh.setdefault("extras", {})["targetNames"] = names

    js = json.dumps(g, separators=(",", ":")).encode("utf-8")
    js += b" " * ((4 - len(js) % 4) % 4)
    bin_total = bytes(bin_data) + bytes(new_bin)
    body = struct.pack("<II", len(js), 0x4E4F534A) + js
    body += struct.pack("<II", len(bin_total), 0x004E4942) + bin_total
    hdr = struct.pack("<III", 0x46546C67, 2, 12 + len(body))
    Path(out_path).write_bytes(hdr + body)
    return names


def dump_obj(path, verts, faces):
    with open(path, "w") as f:
        for v in verts:
            f.write(f"v {v[0]:.6f} {v[1]:.6f} {v[2]:.6f}\n")
        for a, b, c in faces:
            f.write(f"f {a + 1} {b + 1} {c + 1}\n")


# ───────────────────────── main ─────────────────────────

def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--body", required=True)
    ap.add_argument("--out", required=True)
    ap.add_argument("--poses", default="flat,relaxed,fist,grip")
    ap.add_argument("--workdir", default=None)
    ap.add_argument("--selftest-bump", action="store_true")
    ap.add_argument("--no-refit", action="store_true",
                    help="skip the finger-chain refit (see refit_fingers)")
    ap.add_argument("--no-stitch", action="store_true",
                    help="skip cross-sheet stitching (see stitch_components)")
    ap.add_argument("--no-weld", action="store_true",
                    help="solve on the raw (chart-split) ROI graph — reproduces the "
                         "pre-weld baseline for before/after comparison")
    argv = sys.argv
    argv = argv[argv.index("--") + 1:] if "--" in argv else argv[1:]
    args = ap.parse_args(argv)

    g, bin_data = read_glb(args.body)
    prim = g["meshes"][0]["primitives"][0]
    P = load_accessor(g, bin_data, prim["attributes"]["POSITION"])
    N_rest = load_accessor(g, bin_data, prim["attributes"]["NORMAL"]) \
        if "NORMAL" in prim["attributes"] else None
    F = load_accessor(g, bin_data, prim["indices"]).astype(np.int64).reshape(-1, 3)
    skel = Skeleton(g, bin_data)

    if args.selftest_bump:
        c = skel.head(skel.idx["hand_r"])
        delta = np.zeros_like(P)
        d = np.linalg.norm(P - c, axis=1)
        delta[:, 2] += 0.03 * np.exp(-(d / 0.05) ** 2)
        names = emit_morph_glb(g, bin_data,
                               {"hand_bump_r": {"POSITION": delta}}, args.out)
        print(f"[selftest] spliced {names} -> {args.out}; run under BODY_OVERRIDE, the")
        print(f"[selftest] layer should log \"found 1 hand pose(s): bump\" and K raises")
        print(f"[selftest] a 3cm bump on the right palm — proves import+drive end to end")
        return

    poses = [p.strip() for p in args.poses.split(",") if p.strip()]
    if not args.no_refit:
        for side in ("l", "r"):
            refit_fingers(skel, P, side)
    segs_by_side = {}
    for side in ("l", "r"):
        hand, bones = finger_bones(skel, side)
        segs_by_side[side] = bone_segments(skel, hand, bones)
    idx, remap, roi_faces_raw, rim_dup = build_roi(P, F, skel)
    P_dup = P[idx]
    print(f"[roi] {len(idx)} verts / {len(roi_faces_raw)} faces "
          f"({100.0 * len(idx) / len(P):.1f}% of {len(P)})")

    # Weld coincident verts into single graph nodes (see weld_roi). Everything from here
    # on — graph, geodesics, weights, pose, relax — runs on welded nodes; deltas scatter
    # back to every duplicate at emit time.
    if args.no_weld:
        wid = np.arange(len(P_dup))
        P_roi, wgrp, roi_faces = P_dup, np.ones(len(P_dup)), roi_faces_raw
    else:
        wid, P_roi, wgrp = weld_roi(P_dup)
        roi_faces = weld_faces(roi_faces_raw, wid)
    nbr, start, e = roi_graph(P_roi, roi_faces)
    comp = components(len(P_roi), e)
    stitch = (np.zeros((0, 2), dtype=np.int64) if args.no_stitch
              else stitch_components(P_roi, comp))
    if len(stitch):
        nbr, start, e = roi_graph(P_roi, roi_faces, extra=stitch)
        comp2 = components(len(P_roi), e)
        print(f"[stitch] surface sheets {comp.max() + 1} -> {comp2.max() + 1}")
    rim = np.zeros(len(P_roi), dtype=bool)
    np.logical_or.at(rim, wid, rim_dup)
    print(f"[weld] {len(P_dup)} -> {len(P_roi)} nodes "
          f"({int((wgrp > 1).sum())} merged groups, max {int(wgrp.max())}); "
          f"rim {int(rim.sum())}; faces {len(roi_faces_raw)} -> {len(roi_faces)}"
          + ("  [--no-weld]" if args.no_weld else ""))
    lr_roi = np.linalg.norm(P_roi[e[:, 0]] - P_roi[e[:, 1]], axis=1)
    # ROI-border constraints. These edges leave the ROI, so they are absent from `e`;
    # they are the ONLY thing tying the morph to the un-morphed arm (see relax).
    in_roi = np.zeros(len(P), dtype=bool)
    in_roi[idx] = True
    all_e = np.unique(np.sort(np.concatenate(
        [F[:, [0, 1]], F[:, [1, 2]], F[:, [2, 0]]]), axis=1), axis=0)
    xe = all_e[in_roi[all_e[:, 0]] != in_roi[all_e[:, 1]]]
    x_in = np.where(in_roi[xe[:, 0]], xe[:, 0], xe[:, 1])
    x_out = np.where(in_roi[xe[:, 0]], xe[:, 1], xe[:, 0])
    w_of = -np.ones(len(P), dtype=np.int64)
    w_of[idx] = wid
    seam = (w_of[x_in], P[x_out], np.linalg.norm(P[x_in] - P[x_out], axis=1))
    print(f"[seam] {len(xe)} ROI-border edges constrained against the un-morphed body")

    workdir = Path(args.workdir) if args.workdir else None
    if workdir:
        workdir.mkdir(parents=True, exist_ok=True)
        dump_obj(workdir / "rest_hands.obj", P_roi, roi_faces)

    shapes = {}
    report = {"body": args.body, "roi_verts": int(len(idx)),
              "weld": {"enabled": not args.no_weld, "tol_m": WELD_TOL,
                       "nodes": int(len(P_roi)), "merged_groups": int((wgrp > 1).sum()),
                       "max_group": int(wgrp.max()), "rim": int(rim.sum())},
              "stitch_edges": int(len(stitch)), "seam_edges": int(len(xe)),
              "poses": {}}
    # NORMAL deltas: position-only morphs leave lighting on the REST shape — curled
    # fingers shade flat/stale and read as "torn texture". Deltas are measured in the
    # solver's own recomputed-rest frame (exporter normal conventions cancel), over
    # the FULL mesh so ROI-rim verts keep a complete one-ring.
    N_mine_rest = vertex_normals(P, F) if N_rest is not None else None
    for side in ("l", "r"):
        all_bones, W = solve_weights(skel, P_roi, nbr, start, e, side)
        hand_i = skel.idx[f"hand_{side}"]
        hh = skel.head(hand_i)
        side_mask = np.linalg.norm(P_roi - hh, axis=1) < 0.30
        W[~side_mask] = 0.0
        lock_other = ~side_mask
        for pose_name in poses:
            Gp = pose_globals(skel, side, DEFAULT_PARAMS[pose_name])
            posed = lbs(P_roi, skel, all_bones, W, Gp)
            delta0 = posed - P_roi
            delta0[lock_other] = 0.0
            raw = stretch_stats(P_roi, P_roi + delta0, e)
            relaxed, hist = relax(P_roi, delta0, e, lr_roi, lock_other, seam=seam)
            fin = stretch_stats(P_roi, P_roi + relaxed, e)
            seam_mm = float(np.abs(np.linalg.norm(
                (P_roi + relaxed)[seam[0]] - seam[1], axis=1) - seam[2]).max()) * 1000
            sfx = "_l" if side == "l" else "_r"
            travel = {}
            mesh_travel = {}
            for d in DIGITS:
                tn = f"{d}_03{sfx}"
                if tn in Gp:
                    tip = skel.head(skel.idx[tn])
                    pv = (Gp[tn] @ skel.ibm[skel.idx[tn]]) @ np.append(tip, 1.0)
                    travel[d] = float(np.linalg.norm(pv[:3] - tip)) * 100
                    # MESH travel: what the surface near that joint actually does. The
                    # bone number above is scaffolding — it says nothing about whether
                    # any flesh moved, and on this body it read 10cm/finger while the
                    # `_03` joints floated ~5cm outside the mesh entirely.
                    sel = np.linalg.norm(P_roi - tip, axis=1) < 0.020
                    if sel.any():
                        mesh_travel[d] = float(
                            np.linalg.norm(relaxed[sel], axis=1).mean()) * 100
            key = f"{pose_name}_{side}"
            report["poses"][key] = {
                "raw": raw, "relaxed": fin,
                "tip_bone_cm": {k: round(v, 2) for k, v in travel.items()},
                "tip_mesh_cm": {k: round(v, 2) for k, v in mesh_travel.items()},
                "relax_iters": len(hist),
                "viol_edges_left": hist[-1][0], "viol_seam_left": hist[-1][1],
                "seam_max_mm": round(seam_mm, 2),
            }
            full = np.zeros_like(P)
            # scatter: every duplicate of a welded node gets the SAME delta, so
            # coincident pairs keep their rest length exactly and chart seams cannot
            # crack open. Deltas stay indexed by ORIGINAL vertex id — required, or the
            # accessor splice desyncs from the GLB's attribute order.
            full[idx] = relaxed[wid]
            entry = {"POSITION": full.astype(np.float32)}
            if N_mine_rest is not None:
                P_full = P.copy()
                P_full[idx] = P_dup + relaxed[wid]
                n_delta = vertex_normals(P_full, F) - N_mine_rest
                entry["NORMAL"] = n_delta.astype(np.float32)
            shapes[f"hand_{key}"] = entry
            if workdir:
                dump_obj(workdir / f"{key}.obj", P_roi + relaxed, roi_faces)
            print(f"[{key}] raw max {raw['max']:.2f}x p999 {raw['p999']:.2f}x n>5x {raw['n_gt5']:4d}"
                  f" -> relaxed max {fin['max']:.2f}x p999 {fin['p999']:.2f}x n>5x {fin['n_gt5']:4d}"
                  f" needles {fin['needles']:3d}"
                  f" seam {seam_mm:5.1f}mm"
                  f" ({len(hist):2d} it) meshtip cm: "
                  + " ".join(f"{d[:2]}={mesh_travel[d]:.1f}" for d in mesh_travel))

    names = emit_morph_glb(g, bin_data, shapes, args.out)
    if workdir:
        (workdir / "handmorph_report.json").write_text(json.dumps(report, indent=1))
    print(f"[emit] {len(names)} shapes -> {args.out}: {', '.join(names)}")


if __name__ == "__main__":
    main()