"""Classify torn edges (posed stretch >5x) by the dominant joint of each endpoint.
usage: fin_bones.py posed.glb"""
import json, struct, sys, math
from pathlib import Path
from collections import Counter

def read_glb(path):
    d = Path(path).read_bytes()
    length = struct.unpack_from("<I", d, 8)[0]
    off = 12; g = None; 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+i*stride) for i 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)
    M = [[R[i][j]*s[j] for j in range(3)]+[t[i]] for i in range(3)]
    return M+[[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)]

g, b = read_glb(sys.argv[1])
names = [nd.get("name","") for nd in g["nodes"]]
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] if i not in parent else matmul(gm(parent[i]), loc[i])
    memo[i] = m; return m
G = [gm(i) for i in range(len(g["nodes"]))]
skin = g["skins"][0]; joints = skin["joints"]
ibm,_ = acc(g,b,skin["inverseBindMatrices"])
def m16(row): return [[row[c*4+r] for c in range(4)] for r in range(4)]
JM = [matmul(G[joints[j]], m16(ibm[j])) for j in range(len(joints))]
prim = g["meshes"][0]["primitives"][0]
P,_ = acc(g,b,prim["attributes"]["POSITION"])
J,_ = acc(g,b,prim["attributes"]["JOINTS_0"])
W,wt = acc(g,b,prim["attributes"]["WEIGHTS_0"])
wsc = 1.0 if wt==5126 else (1/255 if wt==5121 else 1/65535)
IDX,_ = acc(g,b,prim["indices"])
idx = [i[0] for i in IDX]

def skin_pos(vi):
    p = P[vi]; x=y=z=0.0
    for j,w in zip(J[vi],W[vi]):
        w*=wsc
        if w<=0: continue
        M=JM[j]
        x+=w*(M[0][0]*p[0]+M[0][1]*p[1]+M[0][2]*p[2]+M[0][3])
        y+=w*(M[1][0]*p[0]+M[1][1]*p[1]+M[1][2]*p[2]+M[1][3])
        z+=w*(M[2][0]*p[0]+M[2][1]*p[1]+M[2][2]*p[2]+M[2][3])
    return (x,y,z)

def dom(vi):
    best, bw = None, 0
    for j,w in zip(J[vi],W[vi]):
        w*=wsc
        if w>bw: bw, best = w, j
    return names[joints[best]] if best is not None else "?"

edges = set()
for t in range(0, len(idx), 3):
    a_,b_,c_ = idx[t], idx[t+1], idx[t+2]
    for u,v in ((a_,b_),(b_,c_),(c_,a_)):
        edges.add((u,v) if u<v else (v,u))

pos_cache = {}
def sp(vi):
    if vi not in pos_cache: pos_cache[vi] = skin_pos(vi)
    return pos_cache[vi]

pairs = Counter(); n_bad = 0; maxr = 0
for u,v in edges:
    rl = math.dist(P[u], P[v])
    if rl <= 1e-9: continue
    # cheap prefilter: only edges where an endpoint is finger/hand weighted
    dn_u, dn_v = dom(u), dom(v)
    lu, lv = dn_u.lower(), dn_v.lower()
    keys = ("thumb","index","middle","ring","pinky","hand","lower_arm","wrist")
    if not any(k in lu or k in lv for k in keys): continue
    r = math.dist(sp(u), sp(v)) / rl
    if r > 5:
        n_bad += 1; maxr = max(maxr, r)
        pairs[tuple(sorted((dn_u, dn_v)))] += 1
print(f"edges>5x: {n_bad}  max stretch {maxr:.0f}x")
for (a_,b_), n in pairs.most_common(20):
    print(f"  {n:5d}  {a_}  <->  {b_}")