# graft_hands.py — transplant a working 40-bone hand rig onto a body that was rigged without one.
#
#   blender --background --python tools/graft_hands.py -- \
#       --target <rigged_body.glb> --donor <Ariki_Female_QuatSkin.glb> --out <out.glb>
#   blender --background --python tools/graft_hands.py -- --selftest --donor <...QuatSkin.glb>
#
# WHY THIS EXISTS
# The rig-graft lane (.agents/plans/rig-graft-lane-2026-08-04.md) cuts the hands off the AccuRig
# bait, because close-packed fingers are where every historical hand-mangling came from. That
# leaves the body rigged and the hands unrigged — but the hands do not need solving at all:
#   * the game skeleton's hand chains are FIXED (40 of its 65 joints; verify_body_variant.py
#     gates on "65 joints in identical order"), so there is nothing to discover, only to place;
#   * the nude lane never touched the hands — measured 0.010 mm mean displacement from the Tripo
#     original, p50 exactly 0.000 — so a donor's hand weights fit this mesh as-is;
#   * the shipped clips articulate fingers up to 89 deg relative to each other, so a mitten
#     (fingers weighted as one mass) would visibly flatten four of the six dances.
# So: take the hands from a body that already ships with working ones, aligned at the wrist.
#
# WHY glTF JSON AND NOT BLENDER
# Blender's armature import/export re-derives bone rest orientation from edit-bone head/tail,
# which silently rotates rest poses. That is precisely the failure ariki-game/tools/rig_pose_gate.py
# was written to catch (Godot animation tracks store ABSOLUTE local transforms, so a rewritten
# rest orientation diverges under a clip while a rest-pose comparison still looks fine). Editing
# the node graph directly cannot introduce it. bpy is used only for its KD-tree.
#
# WHAT IT DOES
#   1. Reads the canonical joint ORDER and the hand subtree from the donor.
#   2. Builds a similarity transform M mapping donor space to target space such that the donor's
#      wrist frame lands exactly on the target's wrist frame, with bone lengths scaled by the
#      ratio of forearm lengths (the local scale that matters at the wrist, not global height).
#   3. Re-parents the transformed hand chain under the target's lowerarm, baking the scale into
#      translations so every joint keeps unit scale.
#   4. Copies hand skin weights donor -> target by nearest surface, remapped by joint NAME.
#   5. Rebuilds skin.joints in the donor's canonical order and recomputes inverse-bind matrices.
#
# The IBM convention is not assumed: it is recovered from the target's own body joints and
# asserted before anything is written (see check_ibm_convention).
import json, struct, sys, os, math, argparse
import numpy as np

try:
    from mathutils.kdtree import KDTree
except ImportError:
    KDTree = None

DT = {5120: np.int8, 5121: np.uint8, 5122: np.int16, 5123: np.uint16,
      5125: np.uint32, 5126: np.float32}
CT = {v: k for k, v in DT.items()}
NC = {"SCALAR": 1, "VEC2": 2, "VEC3": 3, "VEC4": 4, "MAT4": 16}
HAND_ROOTS = ("hand_l", "hand_r")


# --------------------------------------------------------------------------- glTF container
class Gltf:
    def __init__(self, path):
        data = open(path, "rb").read()
        total = struct.unpack("<I", data[8:12])[0]
        off, self.js, self.bin = 12, None, b""
        while off < total:
            ln, ty = struct.unpack("<II", data[off:off + 8]); off += 8
            ch = data[off:off + ln]; off += ln
            if ty == 0x4E4F534A: self.js = json.loads(ch)
            elif ty == 0x004E4942: self.bin = ch
        self.path = path
        self.extra = bytearray()          # appended payload for new accessors

    # ---- reading
    def read(self, i):
        a = self.js["accessors"][i]
        dt = np.dtype(DT[a["componentType"]]); nc = NC[a["type"]]; n = a["count"]
        if "bufferView" not in a:
            return np.zeros((n, nc), dtype=dt)
        bv = self.js["bufferViews"][a["bufferView"]]
        off = bv.get("byteOffset", 0) + a.get("byteOffset", 0)
        stride = bv.get("byteStride") or nc * dt.itemsize
        if stride == nc * dt.itemsize:
            return np.frombuffer(self.bin, dtype=dt, count=n * nc, offset=off).reshape(n, nc)
        raw = np.frombuffer(self.bin, dtype=np.uint8, count=stride * n, offset=off).reshape(n, stride)
        return raw[:, :nc * dt.itemsize].copy().view(dt).reshape(n, nc)

    # ---- writing (append-only: existing views are never disturbed)
    def add(self, arr, type_):
        arr = np.ascontiguousarray(arr)
        base = len(self.bin) + len(self.extra)
        pad = (-base) % 4
        self.extra += b"\x00" * pad
        off = base + pad
        raw = arr.tobytes()
        self.extra += raw
        self.js["bufferViews"].append({"buffer": 0, "byteOffset": off, "byteLength": len(raw)})
        acc = {"bufferView": len(self.js["bufferViews"]) - 1,
               "componentType": CT[arr.dtype.type], "count": len(arr), "type": type_}
        if type_ == "VEC3":
            acc["min"] = [float(x) for x in arr.min(axis=0)]
            acc["max"] = [float(x) for x in arr.max(axis=0)]
        self.js["accessors"].append(acc)
        return len(self.js["accessors"]) - 1

    def save(self, out):
        blob = bytes(self.bin) + bytes(self.extra)
        blob += b"\x00" * ((-len(blob)) % 4)
        self.js["buffers"] = [{"byteLength": len(blob)}]
        js = json.dumps(self.js, separators=(",", ":")).encode("utf-8")
        js += b" " * ((-len(js)) % 4)
        hdr = struct.pack("<III", 0x46546C67, 2, 12 + 8 + len(js) + 8 + len(blob))
        with open(out, "wb") as f:
            f.write(hdr)
            f.write(struct.pack("<II", len(js), 0x4E4F534A)); f.write(js)
            f.write(struct.pack("<II", len(blob), 0x004E4942)); f.write(blob)

    # ---- topology helpers
    def parents(self):
        p = {}
        for i, n in enumerate(self.js["nodes"]):
            for c in n.get("children", []): p[c] = i
        return p

    def by_name(self):
        return {n.get("name"): i for i, n in enumerate(self.js["nodes"]) if n.get("name")}

    def local(self, i):
        n = self.js["nodes"][i]
        if "matrix" in n:
            return np.array(n["matrix"], dtype=np.float64).reshape(4, 4).T
        T = np.eye(4); R = np.eye(4); S = np.eye(4)
        T[:3, 3] = n.get("translation", [0, 0, 0])
        R[:3, :3] = quat_mat(n.get("rotation", [0, 0, 0, 1]))
        S[:3, :3] = np.diag(n.get("scale", [1, 1, 1]))
        return T @ R @ S

    def world(self, i, par=None):
        par = par if par is not None else self.parents()
        M = np.eye(4); j = i
        chain = []
        while j is not None:
            chain.append(j); j = par.get(j)
        for j in reversed(chain): M = M @ self.local(j)
        return M

    def body_prim(self):
        best = None
        for mi, m in enumerate(self.js["meshes"]):
            for pi, pr in enumerate(m["primitives"]):
                n = self.js["accessors"][pr["attributes"]["POSITION"]]["count"]
                if best is None or n > best[0]: best = (n, mi, pi)
        return best[1], best[2]

    def skinned_node(self):
        for i, n in enumerate(self.js["nodes"]):
            if "skin" in n and "mesh" in n: return i
        return None


def prune_nodes(g, drop):
    """Delete nodes and remap every index that referred to them. Leaving them orphaned but
    present is not good enough: verify_body_variant.py compares the node-NAME SET against the
    canonical body, and stray nodes fail it (they would also ship as dead scene content)."""
    drop = set(drop)
    keep = [i for i in range(len(g.js["nodes"])) if i not in drop]
    remap = {old: new for new, old in enumerate(keep)}
    g.js["nodes"] = [g.js["nodes"][i] for i in keep]
    for n in g.js["nodes"]:
        if "children" in n:
            kids = [remap[c] for c in n["children"] if c in remap]
            if kids: n["children"] = kids
            else: n.pop("children")
    for sc in g.js.get("scenes", []):
        if "nodes" in sc:
            sc["nodes"] = [remap[i] for i in sc["nodes"] if i in remap]
    for sk in g.js.get("skins", []):
        sk["joints"] = [remap[i] for i in sk["joints"] if i in remap]
        if "skeleton" in sk:
            if sk["skeleton"] in remap: sk["skeleton"] = remap[sk["skeleton"]]
            else: sk.pop("skeleton")
    for an in g.js.get("animations", []):
        for ch in an.get("channels", []):
            t = ch.get("target", {})
            if "node" in t:
                if t["node"] in remap: t["node"] = remap[t["node"]]
                else: ch["_orphan"] = True
        an["channels"] = [c for c in an.get("channels", []) if not c.pop("_orphan", False)]
    return remap


def quat_mat(q):
    x, y, z, w = q
    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)]], dtype=np.float64)


def mat_quat(R):
    """Rotation matrix -> xyzw quaternion, via the numerically stable branch."""
    t = R[0, 0] + R[1, 1] + R[2, 2]
    if t > 0:
        s = math.sqrt(t + 1.0) * 2
        w = 0.25 * s
        x = (R[2, 1] - R[1, 2]) / s; y = (R[0, 2] - R[2, 0]) / s; z = (R[1, 0] - R[0, 1]) / s
    elif R[0, 0] > R[1, 1] and R[0, 0] > R[2, 2]:
        s = math.sqrt(1.0 + R[0, 0] - R[1, 1] - R[2, 2]) * 2
        w = (R[2, 1] - R[1, 2]) / s; x = 0.25 * s
        y = (R[0, 1] + R[1, 0]) / s; z = (R[0, 2] + R[2, 0]) / s
    elif R[1, 1] > R[2, 2]:
        s = math.sqrt(1.0 + R[1, 1] - R[0, 0] - R[2, 2]) * 2
        w = (R[0, 2] - R[2, 0]) / s; x = (R[0, 1] + R[1, 0]) / s
        y = 0.25 * s; z = (R[1, 2] + R[2, 1]) / s
    else:
        s = math.sqrt(1.0 + R[2, 2] - R[0, 0] - R[1, 1]) * 2
        w = (R[1, 0] - R[0, 1]) / s; x = (R[0, 2] + R[2, 0]) / s
        y = (R[1, 2] + R[2, 1]) / s; z = 0.25 * s
    q = np.array([x, y, z, w]); return q / np.linalg.norm(q)


def decompose_unit(M):
    """(translation, xyzw quaternion) with scale stripped — joints stay unit-scale so a scale
    factor never propagates down the finger chain."""
    R = M[:3, :3].copy()
    for k in range(3):
        n = np.linalg.norm(R[:, k])
        if n > 0: R[:, k] /= n
    return M[:3, 3].copy(), mat_quat(R)


def subtree(g, root, par=None):
    out, stack = [], [root]
    while stack:
        i = stack.pop(0); out.append(i)
        stack += g.js["nodes"][i].get("children", [])
    return out


def check_ibm_convention(g, tol=1e-4):
    """Recover, rather than assume, how this file relates inverse-bind matrices to rest poses.
    Returns the mesh-node world matrix that makes IBM == inv(world(joint)) @ Wmesh hold."""
    sk = g.js["skins"][0]
    if "inverseBindMatrices" not in sk: return np.eye(4), 0.0
    ibm = g.read(sk["inverseBindMatrices"]).reshape(-1, 4, 4).transpose(0, 2, 1)
    par = g.parents()
    node = g.skinned_node()
    Wmesh = g.world(node, par) if node is not None else np.eye(4)
    worst = 0.0
    for k, j in enumerate(sk["joints"]):
        pred = np.linalg.inv(g.world(j, par)) @ Wmesh
        worst = max(worst, float(np.abs(pred - ibm[k]).max()))
    return Wmesh, worst


# --------------------------------------------------------------------------- the graft
def unit_scale(M):
    """Same matrix with each basis vector normalised — scale removed, rotation kept."""
    out = M.copy()
    for k in range(3):
        n = np.linalg.norm(out[:3, k])
        if n > 0: out[:3, k] /= n
    return out


def wrist_transform(gt, gd, side, tn, dn, par_t, par_d):
    """Similarity transform mapping DONOR world space to TARGET world space so the donor wrist
    frame lands on the target wrist frame. Scale comes from forearm length — the length that
    governs how far the fingers reach; global body height would be wrong for a differently
    proportioned arm.

    Both wrist frames are stripped of their own scale before composing. Without that, a target
    whose nodes already carry a scale gets it applied TWICE (once inside Wt, once via the forearm
    ratio, which was measured in that same scaled space) — selftest B caught exactly this, as a
    0.055 u placement error that grew toward the finger tips."""
    hand, fore = f"hand_{side}", f"lowerarm_{side}"
    if hand not in tn:
        raise SystemExit(f"target has no '{hand}' node — cannot align. AccuRig must return at "
                         f"least a wrist joint, or use --wrist-from-forearm (not implemented).")
    Wt = gt.world(tn[hand], par_t)
    Wd = gd.world(dn[hand], par_d)
    s = 1.0
    if fore in tn and fore in dn:
        lt = np.linalg.norm(Wt[:3, 3] - gt.world(tn[fore], par_t)[:3, 3])
        ld = np.linalg.norm(Wd[:3, 3] - gd.world(dn[fore], par_d)[:3, 3])
        if ld > 1e-9: s = lt / ld
    Sc = np.eye(4); Sc[:3, :3] *= s
    return unit_scale(Wt) @ Sc @ np.linalg.inv(unit_scale(Wd)), s


def graft(target, donor, out, hand_frac=0.756, verbose=True):
    gt, gd = Gltf(target), Gltf(donor)
    tn, dn = gt.by_name(), gd.by_name()
    par_t, par_d = gt.parents(), gd.parents()

    Wmesh_t, err_t = check_ibm_convention(gt)
    _, err_d = check_ibm_convention(gd)
    if verbose:
        print(f"[ibm] convention residual  target {err_t:.2e}  donor {err_d:.2e}")
    if err_d > 1e-3:
        raise SystemExit(f"donor IBMs do not follow inv(world(joint)) @ Wmesh (residual "
                         f"{err_d:.2e}); refusing to guess a different convention")

    canonical = [gd.js["nodes"][j].get("name") for j in gd.js["skins"][0]["joints"]]
    skin_t = gt.js["skins"][0]
    tj_names = [gt.js["nodes"][j].get("name") for j in skin_t["joints"]]

    # ---- 1. copy each donor hand subtree into the target, transformed to the target wrist
    new_nodes = {}
    dropped = []
    for side in ("l", "r"):
        M, s = wrist_transform(gt, gd, side, tn, dn, par_t, par_d)
        chain = subtree(gd, dn[f"hand_{side}"])
        if verbose:
            print(f"[wrist] {side}: forearm-length scale x{s:.5f}, {len(chain)} joints")
        # drop any hand chain the target already has, so this is a replacement not a duplicate
        if f"hand_{side}" in tn:
            old = subtree(gt, tn[f"hand_{side}"], par_t)
            p = par_t.get(old[0])
            if p is not None:
                gt.js["nodes"][p]["children"] = [c for c in gt.js["nodes"][p].get("children", [])
                                                 if c != old[0]]
            dropped.extend(old)                # removed for real at the end, once indices settle
        for j in chain:
            Wnew = M @ gd.world(j, par_d)
            t, q = decompose_unit(Wnew)
            gt.js["nodes"].append({"name": gd.js["nodes"][j].get("name"),
                                   "translation": [float(x) for x in t],
                                   "rotation": [float(x) for x in q]})
            new_nodes[gd.js["nodes"][j].get("name")] = len(gt.js["nodes"]) - 1
        # re-parent: the chain root goes under the target's forearm, children under their own
        for j in chain:
            nm = gd.js["nodes"][j].get("name")
            kids = [gd.js["nodes"][c].get("name") for c in gd.js["nodes"][j].get("children", [])]
            if kids:
                gt.js["nodes"][new_nodes[nm]]["children"] = [new_nodes[k] for k in kids]
        root_nm = gd.js["nodes"][dn[f'hand_{side}']].get("name")
        fore_i = tn.get(f"lowerarm_{side}")
        if fore_i is None:
            raise SystemExit(f"target has no lowerarm_{side} to parent the hand under")
        gt.js["nodes"][fore_i].setdefault("children", []).append(new_nodes[root_nm])
        # local transforms are currently WORLD; convert to parent-relative
        par_t = gt.parents()
        for j in chain:
            nm = gd.js["nodes"][j].get("name"); i = new_nodes[nm]
            Wnew = M @ gd.world(j, par_d)
            p = par_t.get(i)
            Lp = np.linalg.inv(gt.world(p, par_t)) @ Wnew if p is not None else Wnew
            t, q = decompose_unit(Lp)
            gt.js["nodes"][i]["translation"] = [float(x) for x in t]
            gt.js["nodes"][i]["rotation"] = [float(x) for x in q]
            par_t = gt.parents()

    # ---- 2. rebuild skin.joints in the donor's canonical order
    tn = gt.by_name(); par_t = gt.parents()
    joints_new, missing = [], []
    for nm in canonical:
        if nm in new_nodes: joints_new.append(new_nodes[nm])
        elif nm in tn: joints_new.append(tn[nm])
        else: missing.append(nm)
    if missing:
        raise SystemExit(f"target is missing non-hand joints the donor defines: {missing[:6]}")
    old_index = {nm: k for k, nm in enumerate(tj_names)}
    new_index = {nm: k for k, nm in enumerate(canonical)}

    # ---- 3. weights: donor hand -> target hand vertices, by nearest surface, remapped by NAME
    mi_t, pi_t = gt.body_prim(); prim_t = gt.js["meshes"][mi_t]["primitives"][pi_t]
    mi_d, pi_d = gd.body_prim(); prim_d = gd.js["meshes"][mi_d]["primitives"][pi_d]
    Pt = np.array(gt.read(prim_t["attributes"]["POSITION"]), dtype=np.float64)
    Pd = np.array(gd.read(prim_d["attributes"]["POSITION"]), dtype=np.float64)
    Jd = np.array(gd.read(prim_d["attributes"]["JOINTS_0"]), dtype=np.int64)
    Wd_ = np.array(gd.read(prim_d["attributes"]["WEIGHTS_0"]), dtype=np.float64)
    Jt = np.array(gt.read(prim_t["attributes"]["JOINTS_0"]), dtype=np.int64)
    Wt_ = np.array(gt.read(prim_t["attributes"]["WEIGHTS_0"]), dtype=np.float64)

    dj_names = [gd.js["nodes"][j].get("name") for j in gd.js["skins"][0]["joints"]]
    hand_joint_ids_d = {k for k, nm in enumerate(dj_names)
                        if nm and (nm.startswith(("index", "middle", "ring", "pinky", "thumb"))
                                   or nm in HAND_ROOTS)}
    # donor vertices that are actually skinned to the hand — the geometric definition of "hand"
    is_hand_d = np.array([any(Jd[i, k] in hand_joint_ids_d and Wd_[i, k] > 0 for k in range(4))
                          for i in range(len(Pd))])
    # target hand region by the same bbox-relative cut rigbait_decimate.py uses
    half = np.abs(Pt[:, 0]).max()
    is_hand_t = np.abs(Pt[:, 0]) > hand_frac * half

    # Transform donor hand verts into TARGET MESH-LOCAL space, which is the space POSITION data
    # lives in. M works in world space, so the round trip is:
    #     donor local -> donor world (Wmesh_d) -> target world (M) -> target local (inv Wmesh_t).
    # Skipping the mesh-node matrices only works when both are identity; selftest B caught that
    # as a 0.70 u mean nearest-donor distance where it should have been ~0.
    Wmesh_d, _ = check_ibm_convention(gd)
    inv_Wmesh_t = np.linalg.inv(Wmesh_t)
    Md = {}
    for side in ("l", "r"):
        Md[side], _ = wrist_transform(gt, gd, side, gt.by_name(), dn, gt.parents(), par_d)
    def donor_side(i):
        """Which hand a donor vertex belongs to, from the joint it is actually weighted to —
        not from the sign of x, which assumes a convention this file need not follow."""
        best, bw = None, -1.0
        for k in range(4):
            nm = dj_names[Jd[i, k]]
            if Wd_[i, k] > bw and nm and nm.endswith(("_l", "_r")):
                best, bw = nm[-1], Wd_[i, k]
        return best or "l"

    src_idx = np.where(is_hand_d)[0]
    src_pts = np.empty((len(src_idx), 3))
    for a, i in enumerate(src_idx):
        w = Wmesh_d @ np.append(Pd[i], 1.0)
        src_pts[a] = (inv_Wmesh_t @ (Md[donor_side(i)] @ w))[:3]

    if KDTree is None:
        raise SystemExit("mathutils unavailable — run this under blender --background --python")
    kd = KDTree(len(src_pts))
    for a, p in enumerate(src_pts.tolist()): kd.insert(p, a)
    kd.balance()

    Jt_new = np.zeros_like(Jt); Wt_new = np.zeros_like(Wt_)
    # body vertices keep their weights, remapped to the new joint ordering
    for i in range(len(Pt)):
        if is_hand_t[i]: continue
        for k in range(4):
            nm = tj_names[Jt[i, k]] if Jt[i, k] < len(tj_names) else None
            if nm and nm in new_index and Wt_[i, k] > 0:
                Jt_new[i, k] = new_index[nm]; Wt_new[i, k] = Wt_[i, k]
    moved = 0
    dists = []
    for i in np.where(is_hand_t)[0]:
        a = kd.find(tuple(Pt[i]))[1]
        dists.append(kd.find(tuple(Pt[i]))[2])
        s = src_idx[a]
        for k in range(4):
            nm = dj_names[Jd[s, k]]
            if Wd_[s, k] > 0 and nm in new_index:
                Jt_new[i, k] = new_index[nm]; Wt_new[i, k] = Wd_[s, k]
        moved += 1
    sums = Wt_new.sum(axis=1, keepdims=True)
    Wt_new = np.where(sums > 0, Wt_new / np.maximum(sums, 1e-12), Wt_new)
    if verbose and dists:
        d = np.array(dists)
        print(f"[weights] {moved:,} target hand verts sourced from {len(src_idx):,} donor hand "
              f"verts | nearest-donor distance mean {d.mean():.6f} p99 {np.percentile(d,99):.6f} "
              f"max {d.max():.6f}")

    # ---- 4. inverse-bind matrices for the whole (reordered) skin
    par_t = gt.parents()
    ibm = np.empty((len(joints_new), 4, 4))
    for k, j in enumerate(joints_new):
        ibm[k] = np.linalg.inv(gt.world(j, par_t)) @ Wmesh_t
    skin_t["joints"] = joints_new
    skin_t["inverseBindMatrices"] = gt.add(
        ibm.transpose(0, 2, 1).reshape(-1, 16).astype(np.float32), "MAT4")
    prim_t["attributes"]["JOINTS_0"] = gt.add(Jt_new.astype(np.uint16), "VEC4")
    prim_t["attributes"]["WEIGHTS_0"] = gt.add(Wt_new.astype(np.float32), "VEC4")

    # ---- 5. remove the hand chains we replaced. IBMs are keyed by position in skin.joints, and
    # prune preserves that order, so they stay valid across the reindex.
    if dropped:
        prune_nodes(gt, dropped)
        if verbose: print(f"[prune] removed {len(dropped)} replaced hand nodes")
    gt.save(out)
    if verbose:
        print(f"[done] {out}  ({os.path.getsize(out)/1e6:.2f} MB, {len(joints_new)} joints)")
    return out


# --------------------------------------------------------------------------- self-tests
def selftest(donor, tmp):
    """Two synthetic tests, because the real input (an AccuRig FBX with no hands) does not exist
    yet. Both use the donor as its own target, so the correct answer is known exactly.

      A. IDENTITY  — strip the hand chains, graft them back, expect the original rest poses and
                     weights to return. Validates ordering, re-parenting, IBMs and weight remap.
      B. SIMILARITY — same, but the target is first scaled and rotated by a known amount. The
                     graft must land the hands on the transformed wrist, which is what the real
                     cross-body case needs (AccuRig output differs in scale and orientation).
    """
    ok = True
    ref = Gltf(donor)
    ref_names = [ref.js["nodes"][j].get("name") for j in ref.js["skins"][0]["joints"]]
    par = ref.parents()
    ref_world = {nm: ref.world(ref.js["skins"][0]["joints"][k], par)
                 for k, nm in enumerate(ref_names)}

    for label, scale, deg in (("A identity", 1.0, 0.0), ("B similarity", 0.55, 7.0)):
        tgt = os.path.join(tmp, f"selftest_{label.split()[0]}_target.glb")
        g = Gltf(donor)
        # transform the whole target by a known similarity, applied at the scene roots
        if scale != 1.0 or deg != 0.0:
            c, s_ = math.cos(math.radians(deg)), math.sin(math.radians(deg))
            R = np.array([[c, 0, s_, 0], [0, 1, 0, 0], [-s_, 0, c, 0], [0, 0, 0, 1]])
            S = np.eye(4); S[:3, :3] *= scale
            X = R @ S
            roots = set(range(len(g.js["nodes"]))) - set(g.parents().keys())
            for r in roots:
                L = X @ g.local(r)
                t, q = decompose_unit(L)
                sc = np.linalg.norm(L[:3, 0])
                g.js["nodes"][r].pop("matrix", None)
                g.js["nodes"][r]["translation"] = [float(v) for v in t]
                g.js["nodes"][r]["rotation"] = [float(v) for v in q]
                g.js["nodes"][r]["scale"] = [float(sc)] * 3
        # strip the hand chains from the target's skin (simulating the hands-off bait)
        keep = [j for j, nm in zip(g.js["skins"][0]["joints"],
                                   [g.js["nodes"][x].get("name") for x in g.js["skins"][0]["joints"]])
                if not (nm.startswith(("index", "middle", "ring", "pinky", "thumb")))]
        names_keep = [g.js["nodes"][j].get("name") for j in keep]
        mi, pi = g.body_prim(); prim = g.js["meshes"][mi]["primitives"][pi]
        J = np.array(g.read(prim["attributes"]["JOINTS_0"]), dtype=np.int64)
        W = np.array(g.read(prim["attributes"]["WEIGHTS_0"]), dtype=np.float64)
        old_names = [g.js["nodes"][j].get("name") for j in g.js["skins"][0]["joints"]]
        ni = {nm: k for k, nm in enumerate(names_keep)}
        J2 = np.zeros_like(J); W2 = np.zeros_like(W)
        for i in range(len(J)):
            for k in range(4):
                nm = old_names[J[i, k]]
                # finger weights collapse onto the wrist, as a hands-off rig would have them
                nm = nm if nm in ni else ("hand_l" if nm.endswith("_l") else "hand_r")
                J2[i, k] = ni[nm]; W2[i, k] = W[i, k]
        ibm_old = g.read(g.js["skins"][0]["inverseBindMatrices"]).reshape(-1, 4, 4)
        keepidx = [old_names.index(nm) for nm in names_keep]
        g.js["skins"][0]["joints"] = keep
        g.js["skins"][0]["inverseBindMatrices"] = g.add(
            ibm_old[keepidx].reshape(-1, 16).astype(np.float32), "MAT4")
        prim["attributes"]["JOINTS_0"] = g.add(J2.astype(np.uint16), "VEC4")
        prim["attributes"]["WEIGHTS_0"] = g.add(W2.astype(np.float32), "VEC4")
        g.save(tgt)

        out = os.path.join(tmp, f"selftest_{label.split()[0]}_out.glb")
        print(f"\n=== selftest {label}  (target scaled x{scale}, rotated {deg} deg)")
        graft(tgt, donor, out)

        r = Gltf(out)
        names = [r.js["nodes"][j].get("name") for j in r.js["skins"][0]["joints"]]
        if names != ref_names:
            print(f"  FAIL joint order differs ({len(names)} vs {len(ref_names)})"); ok = False
        else:
            print(f"  PASS joint order — {len(names)} joints, canonical")
        # hand rest poses, compared in the target's own frame (undo the known transform)
        parr = r.parents()
        worst, worstn = 0.0, ""
        for k, nm in enumerate(names):
            if not (nm.startswith(("index", "middle", "ring", "pinky", "thumb")) or nm in HAND_ROOTS):
                continue
            Wg = r.world(r.js["skins"][0]["joints"][k], parr)
            # expected: reference world transformed by the same similarity, scale stripped
            c, s_ = math.cos(math.radians(deg)), math.sin(math.radians(deg))
            R = np.array([[c, 0, s_, 0], [0, 1, 0, 0], [-s_, 0, c, 0], [0, 0, 0, 1]])
            S = np.eye(4); S[:3, :3] *= scale
            exp = R @ S @ ref_world[nm]
            d = np.linalg.norm(Wg[:3, 3] - exp[:3, 3])
            if d > worst: worst, worstn = d, nm
        tol = 1e-5 * max(scale, 1e-3)
        print(f"  {'PASS' if worst < 1e-4 else 'FAIL'} hand joint placement — worst origin error "
              f"{worst:.3e} u at {worstn}")
        ok &= worst < 1e-4
        # weights: every hand vertex should recover the donor's own weights
        mi, pi = r.body_prim(); pr = r.js["meshes"][mi]["primitives"][pi]
        Jn = np.array(r.read(pr["attributes"]["JOINTS_0"]), dtype=np.int64)
        Wn = np.array(r.read(pr["attributes"]["WEIGHTS_0"]), dtype=np.float64)
        mi0, pi0 = ref.body_prim(); pr0 = ref.js["meshes"][mi0]["primitives"][pi0]
        J0 = np.array(ref.read(pr0["attributes"]["JOINTS_0"]), dtype=np.int64)
        W0 = np.array(ref.read(pr0["attributes"]["WEIGHTS_0"]), dtype=np.float64)
        P0 = np.array(ref.read(pr0["attributes"]["POSITION"]), dtype=np.float64)
        half = np.abs(P0[:, 0]).max()
        hand = np.abs(P0[:, 0]) > 0.756 * half
        def as_dict(J, W, i):
            return {J[i, k]: round(float(W[i, k]), 4) for k in range(4) if W[i, k] > 1e-6}
        same = sum(1 for i in np.where(hand)[0] if as_dict(Jn, Wn, i) == as_dict(J0, W0, i))
        tot = int(hand.sum())
        print(f"  {'PASS' if same == tot else 'WARN'} hand weights recovered exactly on "
              f"{same:,}/{tot:,} hand verts ({100*same/max(tot,1):.2f}%)")
        wsum = Wn.sum(axis=1)
        print(f"  {'PASS' if abs(wsum-1).max() < 1e-3 else 'FAIL'} weights normalised "
              f"(max deviation {abs(wsum-1).max():.2e})")
        ok &= abs(wsum - 1).max() < 1e-3
    print(f"\nSELFTEST {'PASS' if ok else 'FAIL'}")
    return 0 if ok else 2


def main():
    argv = sys.argv[sys.argv.index("--") + 1:] if "--" in sys.argv else sys.argv[1:]
    ap = argparse.ArgumentParser()
    ap.add_argument("--target"); ap.add_argument("--donor", required=True)
    ap.add_argument("--out"); ap.add_argument("--selftest", action="store_true")
    ap.add_argument("--tmp", default=".")
    a = ap.parse_args(argv)
    if a.selftest:
        raise SystemExit(selftest(a.donor, a.tmp))
    if not a.target or not a.out:
        raise SystemExit("--target and --out are required unless --selftest")
    graft(a.target, a.donor, a.out)


main()