"""Find corrupt skin bindings: verts bound to a joint that is implausibly far away in
REST pose (long bind lever arm), and verts bound across the body midline to the opposite
hand's bones. These are invisible at rest and only explode once the joint rotates.

usage: skin_lever_audit.py body.glb [--lever CM] [--dump N]
"""
import json, struct, sys, math
from pathlib import Path
from collections import Counter, defaultdict

FING = ("thumb", "index", "middle", "ring", "pinky")
LEVER_CM = 20.0
DUMP = 0

args = [a for a in sys.argv[1:]]
path = args[0]
if "--lever" in args:
    LEVER_CM = float(args[args.index("--lever") + 1])
if "--dump" in args:
    DUMP = int(args[args.index("--dump") + 1])


def read_glb(p):
    d = Path(p).read_bytes()
    length = struct.unpack_from("<I", d, 8)[0]
    off = 12
    g = b_ = None
    while off < length:
        clen, ct = struct.unpack_from("<II", d, off)
        off += 8
        if ct == 0x4E4F534A:
            g = json.loads(d[off:off + clen])
        else:
            b_ = d[off:off + clen]
        off += clen
    return g, b_


def acc(g, b, i):
    a = g["accessors"][i]
    bv = g["bufferViews"][a["bufferView"]]
    nc = {"SCALAR": 1, "VEC2": 2, "VEC3": 3, "VEC4": 4, "MAT4": 16}[a["type"]]
    fmt = {5121: "B", 5123: "H", 5125: "I", 5126: "f"}[a["componentType"]]
    size = struct.calcsize(fmt) * nc
    stride = bv.get("byteStride") or size
    off = bv.get("byteOffset", 0) + a.get("byteOffset", 0)
    return [struct.unpack_from("<%d%s" % (nc, fmt), b, off + k * stride)
            for k in range(a["count"])], a["componentType"]


def quat_mat(q):
    x, y, z, w = q
    return [[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):
    t = nd.get("translation", [0, 0, 0])
    r = nd.get("rotation", [0, 0, 0, 1])
    s = nd.get("scale", [1, 1, 1])
    R = quat_mat(r)
    return [[R[i][j] * s[j] for j in range(3)] + [t[i]] for i in range(3)] + [[0, 0, 0, 1]]


def matmul(A, B):
    return [[sum(A[i][k] * B[k][j] for k in range(4)) for j in range(4)] for i in range(4)]


def global_mats(g):
    loc = [node_local(nd) for nd in g["nodes"]]
    parent = {}
    for i, nd in enumerate(g["nodes"]):
        for c in nd.get("children", []):
            parent[c] = i
    memo = {}

    def gm(i):
        if i in memo:
            return memo[i]
        m = loc[i]
        p = parent.get(i)
        if p is not None:
            m = matmul(gm(p), m)
        memo[i] = m
        return m

    return [gm(i) for i in range(len(g["nodes"]))]


g, b = read_glb(path)
names = [nd.get("name", "") for nd in g["nodes"]]
GM = global_mats(g)

print(f"== {Path(path).name} ==")
print(f"   lever threshold {LEVER_CM:.0f} cm\n")

for mi, mesh in enumerate(g.get("meshes", [])):
    for pi, prim in enumerate(mesh.get("primitives", [])):
        att = prim["attributes"]
        if "JOINTS_0" not in att:
            continue
        skin_idx = next((nd.get("skin") for nd in g["nodes"]
                         if nd.get("mesh") == mi and "skin" in nd), None)
        if skin_idx is None:
            continue
        joints = g["skins"][skin_idx]["joints"]
        jname = [names[j] for j in joints]
        # joint rest world positions
        jpos = [(GM[j][0][3], GM[j][1][3], GM[j][2][3]) for j in joints]

        P, _ = acc(g, b, att["POSITION"])
        J, _ = acc(g, b, att["JOINTS_0"])
        W, wt = acc(g, b, att["WEIGHTS_0"])
        wsc = 1.0 if wt == 5126 else (1 / 255 if wt == 5121 else 1 / 65535)

        long_lever = []
        cross = []
        by_joint = Counter()
        cross_by_joint = Counter()

        for vi, (p, jrow, wrow) in enumerate(zip(P, J, W)):
            for j, w in zip(jrow, wrow):
                w *= wsc
                if w <= 0.001:
                    continue
                n = jname[j]
                nl = n.lower()
                if not any(t in nl for t in FING):
                    continue
                jp = jpos[j]
                d = math.dist(p, jp) * 100.0  # cm (glTF metres)
                if d > LEVER_CM:
                    long_lever.append((d, vi, n, w, p))
                    by_joint[n] += 1
                # cross-hand: vert on opposite side of midline from the joint
                if nl.endswith("_r") and p[0] > 0.02:
                    cross.append((d, vi, n, w, p)); cross_by_joint[n] += 1
                elif nl.endswith("_l") and p[0] < -0.02:
                    cross.append((d, vi, n, w, p)); cross_by_joint[n] += 1

        print(f"  mesh[{mi}] '{mesh.get('name','')}' prim{pi}: {len(P)} verts")
        print(f"    finger refs with lever > {LEVER_CM:.0f} cm : {len(long_lever)}")
        if long_lever:
            mx = max(long_lever)
            print(f"      worst {mx[0]:.1f} cm  vert {mx[1]}  joint {mx[2]}  w={mx[3]:.3f}")
            for n, c in by_joint.most_common(10):
                print(f"        {n:22s} {c}")
        print(f"    cross-midline finger refs        : {len(cross)}")
        if cross:
            mx = max(cross)
            print(f"      worst {mx[0]:.1f} cm  vert {mx[1]}  joint {mx[2]}  w={mx[3]:.3f}")
            for n, c in cross_by_joint.most_common(10):
                print(f"        {n:22s} {c}")
        for d, vi, n, w, p in sorted(long_lever, reverse=True)[:DUMP]:
            print(f"      v{vi} {n} w={w:.3f} lever={d:.1f}cm pos=({p[0]*100:.1f},{p[1]*100:.1f},{p[2]*100:.1f})cm")