"""LmGirl fixed Mixamo-41 preparation, weight variants and motion review data.

All website weights are discarded. Fitting preserves bone axes and hierarchy;
the source animation can keep its rotations with target bind translations.
"""
import json
import importlib.util
import hashlib
import math
import sys
from pathlib import Path

import bpy
import numpy as np
from mathutils import Matrix, Vector

ROOT=Path(__file__).resolve().parents[2]
OUT=ROOT/'Tools/output/lmgirl_tpose_studio/rig_mixamo41'
STATIC=ROOT/'Tools/output/lmgirl_tpose_studio/static/LmGirlTpose.fbx'
REFERENCE=ROOT/'AIClient/Assets/Bundles/Characters/Human/test/Boy_mixamo/CampGuideStudio.fbx'

def write(name,data):
    OUT.mkdir(parents=True,exist_ok=True)
    (OUT/name).write_text(json.dumps(data,indent=2),encoding='utf-8')

def mesh_points(mesh):return np.array([mesh.matrix_world@v.co for v in mesh.data.vertices])

def load():
    bpy.ops.wm.read_factory_settings(use_empty=True)
    bpy.ops.import_scene.fbx(filepath=str(STATIC))
    mesh=next(o for o in bpy.context.scene.objects if o.type=='MESH')
    assert not any(o.type=='ARMATURE' for o in bpy.context.scene.objects)
    # Project static FBX contains Y-up mesh coordinates; glTF needs Blender Z-up.
    if np.argmax(np.ptp(mesh_points(mesh),axis=0))==1:
        mesh.matrix_world=Matrix.Rotation(math.pi/2,4,'X')@mesh.matrix_world
    bpy.ops.object.select_all(action='DESELECT');mesh.select_set(True);bpy.context.view_layer.objects.active=mesh
    bpy.ops.object.transform_apply(location=True,rotation=True,scale=True)
    pts=mesh_points(mesh)
    lo,hi=pts.min(0),pts.max(0)
    mesh.location-=Vector(((lo[0]+hi[0])/2,(lo[1]+hi[1])/2,lo[2]))
    bpy.ops.object.transform_apply(location=True,rotation=False,scale=False)
    bpy.ops.import_scene.fbx(filepath=str(REFERENCE))
    arm=next(o for o in bpy.context.scene.objects if o.type=='ARMATURE')
    for obj in list(bpy.context.scene.objects):
        if obj not in (mesh,arm):bpy.data.objects.remove(obj,do_unlink=True)
    arm.animation_data_clear()
    arm.matrix_world=Matrix.Rotation(math.pi/2,4,'X')@arm.matrix_world
    bpy.context.view_layer.update()
    return mesh,arm

def inspect():
    mesh,arm=load();pts=mesh_points(mesh);lo,hi=pts.min(0),pts.max(0);h=hi[2]-lo[2]
    bands=[]
    for z in np.linspace(.03,.97,48):
        a=pts[np.abs(pts[:,2]-z*h)<h*.012]
        if len(a):bands.append({'z':round(float(z),3),'n':len(a),'min':a.min(0).tolist(),'max':a.max(0).tolist(),'mean':a.mean(0).tolist()})
    data={'bounds':[lo.tolist(),hi.tolist()],'vertices':len(pts),'triangles':sum(len(p.vertices)-2 for p in mesh.data.polygons),'bands':bands,
          'reference_bones':{b.name:list(arm.matrix_world@b.head_local) for b in arm.data.bones}}
    write('inspection.json',data)
    print('INSPECT',json.dumps(data))
    # Plain clay inspection; source UV remains untouched.
    scene=bpy.context.scene;scene.render.engine='BLENDER_WORKBENCH'
    scene.display.shading.light='STUDIO';scene.display.shading.color_type='SINGLE';scene.display.shading.single_color=(.65,.68,.72)
    scene.display.shading.show_shadows=True;scene.display.shading.show_cavity=True
    camera=bpy.data.objects.new('InspectionCamera',bpy.data.cameras.new('InspectionCamera'));scene.collection.objects.link(camera)
    camera.location=(0,-h*3,h*.5);camera.rotation_euler=(Vector((0,0,h*.5))-camera.location).to_track_quat('-Z','Y').to_euler()
    camera.data.type='ORTHO';camera.data.ortho_scale=h*1.45;scene.camera=camera
    scene.render.resolution_x=1100;scene.render.resolution_y=1000;scene.render.resolution_percentage=100
    scene.render.filepath=str(OUT/'static_front.png');bpy.ops.render.render(write_still=True)

def prepare():
    mesh,arm=load()
    old={b.name:arm.matrix_world@b.matrix_local for b in arm.data.bones}
    positions={
        'Hips':(0,-.005,.610),'Spine':(0,-.005,.660),'Spine1':(0,-.005,.730),
        'Spine2':(0,.005,.825),'Neck':(0,.010,.925),'Head':(0,.010,.965),'HeadTop_End':(0,.010,1.150)}
    for side,sign in [('Left',1),('Right',-1)]:
        for suffix,p in {'Shoulder':(.043,.018,.900),'Arm':(.120,.018,.877),
                         'ForeArm':(.248,.018,.862),'Hand':(.352,.018,.866),
                         'UpLeg':(.057,-.005,.571),'Leg':(.060,.006,.312),
                         'Foot':(.056,.018,.060),'ToeBase':(.056,-.039,.016),
                         'Toe_End':(.056,-.058,.016)}.items():
            positions[side+suffix]=(p[0]*sign,p[1],p[2])
        anchor=old['mixamorig:'+side+'Hand'].translation
        for name,m in old.items():
            if name.startswith('mixamorig:'+side+'Hand') and not name.endswith('Hand'):
                d=m.translation-anchor
                positions[name.split(':')[-1]]=Vector(positions[side+'Hand'])+Vector((d.x*.95,d.y*.35,d.z*.6))
    bpy.context.view_layer.objects.active=arm
    bpy.ops.object.mode_set(mode='EDIT')
    inv=arm.matrix_world.inverted()
    shapes={b.name:(b.tail-b.head,b.roll) for b in arm.data.edit_bones}
    for b in arm.data.edit_bones:b.use_connect=False
    for bone in arm.data.edit_bones:
        target=inv@Vector(positions[bone.name.split(':')[-1]])
        bone.head=target;bone.tail=target+shapes[bone.name][0];bone.roll=shapes[bone.name][1]
    bpy.ops.object.mode_set(mode='OBJECT');bpy.context.view_layer.update()
    drift=max(old[b.name].to_quaternion().rotation_difference((arm.matrix_world@b.matrix_local).to_quaternion()).angle for b in arm.data.bones)
    assert drift<.001,drift
    for modifier in list(mesh.modifiers):
        if modifier.type=='ARMATURE':mesh.modifiers.remove(modifier)
    mesh.vertex_groups.clear()
    for b in arm.data.bones:mesh.vertex_groups.new(name=b.name)
    mesh.vertex_groups['mixamorig:Hips'].add(list(range(len(mesh.data.vertices))),1,'REPLACE')
    for i,b in enumerate(arm.data.bones):
        if b.name!='mixamorig:Hips':
            mesh.vertex_groups['mixamorig:Hips'].remove([i]);mesh.vertex_groups[b.name].add([i],1,'REPLACE')
    world=mesh.matrix_world.copy();mesh.parent=arm;mesh.matrix_world=world
    modifier=mesh.modifiers.new('TokenRigInputSkeleton','ARMATURE');modifier.object=arm
    arm.data.pose_position='REST'
    bpy.context.view_layer.update()
    pts=mesh_points(mesh)
    # Keep material payload small; final textured assets reuse the exact source
    # UV and local baked maps. Transfer is still exercised on this original mesh.
    write('input_manifest.json',{'static_mesh_source':str(STATIC),'skeleton_source':str(REFERENCE),
        'bones':{b.name:{'parent':b.parent.name if b.parent else None,'world_matrix':[list(row) for row in arm.matrix_world@b.matrix_local],
                        'head':list(arm.matrix_world@b.head_local)} for b in arm.data.bones},
        'vertices':len(pts),'triangles':sum(len(p.vertices)-2 for p in mesh.data.polygons),
        'height':float(np.ptp(pts[:,2])),'max_preserved_axis_rotation_error_radians':float(drift),
        'input_skin':'Only non-anatomical keepalive placeholders; no website weights',
        'fit':'Manually specified anatomical landmark positions for this mesh; original bone axes and hierarchy retained',
        'landmarks':{n:list(p) for n,p in positions.items()},'blender_up':'Z','blender_front':'-Y'})
    bpy.ops.wm.save_as_mainfile(filepath=str(OUT/'prepared.blend'),compress=True)
    bpy.ops.export_scene.gltf(filepath=str(OUT/'source.glb'),export_format='GLB',export_yup=True,export_skins=True,
                              export_animations=False,export_def_bones=False)
    print('PREPARED',len(arm.data.bones),len(pts),'axis drift',drift)

def process():
    from mathutils.kdtree import KDTree
    from collections import deque
    bpy.ops.wm.open_mainfile(filepath=str(OUT/'prepared.blend'))
    arm=next(o for o in bpy.context.scene.objects if o.type=='ARMATURE')
    mesh=next(o for o in bpy.context.scene.objects if o.type=='MESH')
    arm.data.pose_position='POSE'
    names=[b.name for b in arm.data.bones];pts=mesh_points(mesh)
    parents=[names.index(b.parent.name) if b.parent else -1 for b in arm.data.bones]
    heads=np.array([arm.matrix_world@b.head_local for b in arm.data.bones])
    matrices={b.name:arm.matrix_world@b.matrix_local for b in arm.data.bones}
    variants={}
    for tag in ['raw','post']:
        decoded=json.loads((OUT/(tag+'_decoded_weights.json')).read_text())
        tree=KDTree(len(decoded))
        for i,d in enumerate(decoded):
            x,y,z=d['position_y_up'];tree.insert(Vector((x,-z,y)),i)
        tree.balance();w=np.zeros((len(pts),41));dist=[]
        for i,p in enumerate(pts):
            _,j,distance=tree.find(Vector(p));dist.append(distance)
            for name,value in decoded[j]['weights'].items():w[i,names.index(name)]+=value
        assert max(dist)<1e-5
        w/=w.sum(1,keepdims=True);variants[tag]=w
    # Run the upstream pure-array repair on the RESTORED original graph and
    # world coordinates. TokenRig's temporary graph/quantized axes are excluded.
    spec=importlib.util.spec_from_file_location('weight_repair_upstream',OUT/'weight_repair_upstream.py')
    repair=importlib.util.module_from_spec(spec);sys.modules[spec.name]=repair;spec.loader.exec_module(repair)
    hop=np.full((41,41),99,dtype=int);adj=[[] for _ in names]
    for i,p in enumerate(parents):
        if p>=0:adj[i].append(p);adj[p].append(i)
    for i in range(41):
        hop[i,i]=0;q=deque([i])
        while q:
            a=q.popleft()
            for b in adj[a]:
                if hop[i,b]==99:hop[i,b]=hop[i,a]+1;q.append(b)
    ids=np.argsort(-variants['post'],axis=1)[:,:4]
    values=np.take_along_axis(variants['post'],ids,axis=1)
    clean_ids,clean_weights,stats=repair.repair_weights_array(pts,ids,values,heads,hop,joint_names=names)
    clean=np.zeros_like(variants['post'])
    np.put_along_axis(clean,clean_ids,clean_weights,axis=1);clean/=clean.sum(1,keepdims=True)
    variants['clean']=clean
    eps=clean.copy();eps[eps<.01]=0;eps/=eps.sum(1,keepdims=True);variants['clean_eps']=eps
    # Conservative character-specific alternative: the low-poly hands have no
    # separated finger loops, and the braid must follow the head. Preserve 41
    # joints, but use rigid hand/head influences in those identifiable regions.
    local=variants['raw'].copy()
    local_changes=[]
    for i,p in enumerate(pts):
        before=local[i].copy()
        if p[2]>.945:
            local[i]=0;local[i,names.index('mixamorig:Head')]=1
        elif abs(p[0])>.335:
            side='Left' if p[0]>0 else 'Right'
            local[i]=0;local[i,names.index('mixamorig:'+side+'Hand')]=1
        elif abs(p[0])>.15 and .77<p[2]<.94:
            side='Left' if p[0]>0 else 'Right'
            allowed=np.array([n in ['mixamorig:'+side+'Arm','mixamorig:'+side+'ForeArm','mixamorig:'+side+'Hand'] for n in names])
            # Merge finger channels into wrist before applying the sleeve mask.
            amount=sum(local[i,j] for j,n in enumerate(names) if n.startswith('mixamorig:'+side+'Hand'))
            local[i,~allowed]=0;local[i,names.index('mixamorig:'+side+'Hand')]=amount
            if local[i].sum()<1e-8:local[i,names.index('mixamorig:'+side+'ForeArm')]=1
            local[i]/=local[i].sum()
        elif p[2]<.48:
            side='Left' if p[0]>0 else 'Right'
            allowed=np.array([n=='mixamorig:Hips' or n in ['mixamorig:'+side+x for x in ['UpLeg','Leg','Foot','ToeBase','Toe_End']] for n in names])
            local[i,~allowed]=0
            if local[i].sum()<1e-8:local[i,names.index('mixamorig:'+side+'Leg')]=1
            local[i]/=local[i].sum()
        if np.abs(before-local[i]).sum()>.001:local_changes.append(i)
    local[local<.01]=0;local/=local.sum(1,keepdims=True);variants['local']=local
    write('cleanup.json',{'source':'https://github.com/stevelittlefish/ComfyUI-SkinTokens-NoBlender/blob/main/skintokens/weight_repair.py',
        'algorithm':'Upstream repair_weights_array, original restored skeleton graph and world-space coordinates',
        'min_hops':5,'repaired_vertices':stats.n_vertices_repaired,'removed_influences':stats.n_influences_removed,
        'removed_by_bone':stats.removed_by_bone,'pairs':{' | '.join(k):v for k,v in stats.comingled_pairs.items()},
        'epsilon':.01,'epsilon_changed_vertices':int(np.any(np.abs(eps-clean)>1e-7,axis=1).sum())})
    write('local_constraints.json',{'base':'raw predicted weights; no voxel sharpening',
        'rules':['head/braid above 0.945m rigid to Head','hands outside |x|=0.335m rigid to same-side Hand',
                 'sleeves outside |x|=0.15m and 0.77<z<0.94m allow same-side Arm/ForeArm/Hand only',
                 'below z=0.48m allow Hips and same-side UpLeg/Leg/Foot/Toe only; crotch blending retained','epsilon=0.01; renormalize'],
        'changed_vertices_over_0_001':len(local_changes),'vertices':local_changes,
        'limitation':'41 bones retained; separated finger articulation is intentionally not used for this mitten-like low-poly hand mesh; no mesh/topology changes'})
    C=Matrix(((1,0,0,0),(0,0,1,0),(0,1,0,0),(0,0,0,1)))
    def trs(m):
        p,q,s=m.decompose();return list(p)+[q.x,q.y,q.z,q.w]+list(s)
    world={n:C@m@C for n,m in matrices.items()}
    bones=[]
    for i,n in enumerate(names):
        local=world[names[parents[i]]].inverted()@world[n] if parents[i]>=0 else world[n]
        bones.append({'name':n,'parent':parents[i],'rest':trs(local),'bindpose':[list(r) for r in world[n].inverted()]})
    source=json.loads((ROOT/'Tools/output/campguide_studio/rig_mixamo41/unity_review/review_data.json').read_text())
    assert names==[b['name'] for b in source['bones']]
    # Foot/hip span ratio determines root travel and vertical motion scaling.
    ratio=(.610-.060)/(.5381249189-.0935430303)
    offsets=[(np.array(b['rest'][:3])-np.array(s['rest'][:3])*ratio).tolist() for b,s in zip(bones,source['bones'])]
    from mathutils import Quaternion
    clip=source['clips'][0]
    poses=[]
    for sample in clip['samples']:
        m=[]
        for i,values in enumerate(sample):
            p=Vector(np.array(values[:3])*ratio+offsets[i]);q=Quaternion((values[6],*values[3:6]))
            local=Matrix.LocRotScale(p,q,Vector(values[7:10]));m.append(m[parents[i]]@local if parents[i]>=0 else local)
        poses.append(m)
    rest=np.array([C@Vector(p) for p in pts]);h=float(np.ptp(pts[:,2]))
    edges=np.array([list(e.vertices) for e in mesh.data.edges]);lengths=np.linalg.norm(rest[edges[:,0]]-rest[edges[:,1]],axis=1)
    masks={'all':lengths>h*.002,'hands':(np.abs(pts[:,0])>.32),'neck':(pts[:,2]>.91)&(pts[:,2]<.99)&(np.abs(pts[:,0])<.1)}
    masks={k:(v if k=='all' else v[edges[:,0]]&v[edges[:,1]]&(lengths>h*.002)) for k,v in masks.items()}
    report={};result={'bones':bones,'translationScale':ratio,'translationOffsets':offsets,'variants':[],
                      'height':h,'walkingFrames':len(poses),'reuseClips':'CampGuide original native AnimationClips, angular curves unchanged'}
    for tag,w in variants.items():
        mesh.vertex_groups.clear()
        for n in names:mesh.vertex_groups.new(name=n)
        for i,row in enumerate(w):
            for j in np.flatnonzero(row>1e-8):mesh.vertex_groups[names[j]].add([i],float(row[j]),'REPLACE')
        arm.animation_data_clear()
        for b in arm.pose.bones:b.matrix_basis=Matrix.Identity(4)
        bpy.context.view_layer.update()
        bpy.ops.wm.save_as_mainfile(filepath=str(OUT/('lmgirl_'+tag+'.blend')),compress=True)
        bpy.ops.export_scene.gltf(filepath=str(OUT/('lmgirl_'+tag+'.glb')),export_format='GLB',export_yup=True,export_skins=True,export_animations=False,export_def_bones=False)
        maxs={k:0 for k in masks};over={k:0 for k in masks};samples=[];floors=[];worst={};p99=0
        p4=np.column_stack((rest,np.ones(len(rest))))
        for frame,pose in enumerate(poses):
            deformed=np.zeros_like(rest)
            for j,n in enumerate(names):
                transform=np.array(pose[j]@world[n].inverted())
                deformed+=w[:,j,None]*(p4@transform.T)[:,:3]
            ratios=np.linalg.norm(deformed[edges[:,0]]-deformed[edges[:,1]],axis=1)/np.maximum(lengths,1e-10)
            floors.append(float(deformed[:,1].min()))
            ei=int(np.argmax(np.where(masks['all'],ratios,0)))
            if float(ratios[ei])>worst.get('ratio',0):
                worst={'frame':frame+1,'ratio':float(ratios[ei]),'vertices':edges[ei].tolist(),
                       'rest':pts[edges[ei]].tolist(),'weights':[{names[j]:float(w[i,j]) for j in np.flatnonzero(w[i]>1e-8)} for i in edges[ei]]}
            p99=max(p99,float(np.percentile(ratios[masks['all']],99)))
            for k,mask in masks.items():
                maxs[k]=max(maxs[k],float(ratios[mask].max()));over[k]=max(over[k],int((ratios[mask]>2).sum()))
            samples.append({'frame':frame+1,'maxStretch':float(ratios[masks['all']].max()),'floorMin':floors[-1]})
        # Mesh exporter: exact source topology/UV/corner normals, reflected Y-up.
        transform=C@mesh.matrix_world;normal_matrix=transform.to_3x3().inverted().transposed()
        mesh.data.calc_loop_triangles();uv_layer=mesh.data.uv_layers.active
        vertex=[];normal=[];uv=[];weights=[];tri=[];lookup={};first={}
        for triangle in mesh.data.loop_triangles:
            ids=[]
            for li in triangle.loops:
                vi=mesh.data.loops[li].vertex_index;n=normal_matrix@mesh.data.corner_normals[li].vector;n.normalize()
                tex=list(uv_layer.data[li].uv);key=(vi,*[round(v,6) for v in tex],*[round(v,6) for v in n])
                if key not in lookup:
                    index=len(vertex);lookup[key]=index;first.setdefault(vi,index)
                    vertex.append(list(transform@mesh.data.vertices[vi].co));normal.append(list(n));uv.append(tex)
                    weights.append([[int(j),float(w[vi,j])] for j in np.argsort(-w[vi]) if w[vi,j]>1e-8])
                ids.append(lookup[key])
            tri.extend([ids[0],ids[2],ids[1]])
        audit_edges=[[first[a],first[b],float(lengths[i]),1|(2 if masks['hands'][i] else 0)|(4 if masks['neck'][i] else 0)] for i,(a,b) in enumerate(edges) if masks['all'][i]]
        result['variants'].append({'name':tag,'vertices':vertex,'normals':normal,'uv':uv,'triangles':tri,'weights':weights,
                                   'auditEdges':audit_edges,'originalVertices':len(pts),'trianglesCount':len(tri)//3})
        foot_ids=[i for i,n in enumerate(names) if any(x in n for x in ['Leg','Foot','Toe'])]
        report[tag]={'maxStretch':maxs,'maxEdgesOver2x':over,'floorMinRange':[min(floors),max(floors)],'samples':samples,'worstEdge':worst,'maxP99':p99,
            'handMaxLowerBodyInfluence':float(w[np.abs(pts[:,0])>.32][:,foot_ids].sum(1).max()),
            'maxInfluences':int((w>1e-8).sum(1).max()),'minimumNonzeroWeight':float(w[w>1e-8].min())}
        # Bake Walking for external review; native Unity uses the shared clip.
        for b in arm.pose.bones:b.rotation_mode='QUATERNION'
        for f,pose in enumerate(poses,1):
            bpy.context.scene.frame_set(f)
            for i,b in enumerate(arm.pose.bones):
                b.matrix=arm.matrix_world.inverted()@C@pose[i]@C;bpy.context.view_layer.update()
                for attr in ['location','rotation_quaternion','scale']:b.keyframe_insert(attr,frame=f)
        bpy.context.scene.frame_start=1;bpy.context.scene.frame_end=len(poses);bpy.context.scene.render.fps=30
        bpy.ops.export_scene.gltf(filepath=str(OUT/('Walking_'+tag+'.glb')),export_format='GLB',export_skins=True,export_animations=True,
                                  export_animation_mode='ACTIVE_ACTIONS',export_force_sampling=True,export_def_bones=False)
        print('VARIANT',tag,json.dumps({k:v for k,v in report[tag].items() if k!='samples'}))
    write('walking_comparison.json',report)
    (OUT/'unity_review_data.json').write_text(json.dumps(result,separators=(',',':')),encoding='utf-8')
    print('PROCESSED',len(variants),'variants')

def topology():
    from collections import deque
    bpy.ops.wm.open_mainfile(filepath=str(OUT/'prepared.blend'))
    mesh=next(o for o in bpy.context.scene.objects if o.type=='MESH');pts=mesh_points(mesh)
    # Restrict the graph BELOW the crotch: otherwise two separate pant legs
    # would legitimately connect through the pelvis.
    allowed=pts[:,2]<.50;adj=[[] for _ in pts]
    for e in mesh.data.edges:
        a,b=e.vertices
        if allowed[a] and allowed[b]:adj[a].append(b);adj[b].append(a)
    left=int(np.argmin(np.linalg.norm(pts-np.array((.06,.01,.26)),axis=1)))
    right=int(np.argmin(np.linalg.norm(pts-np.array((-.06,.01,.26)),axis=1)))
    prev={left:None};queue=deque([left])
    while queue and right not in prev:
        a=queue.popleft()
        for b in adj[a]:
            if b not in prev:prev[b]=a;queue.append(b)
    path=[];a=right if right in prev else None
    while a is not None:path.append(a);a=prev[a]
    path.reverse()
    faces=[p.index for p in mesh.data.polygons if all(pts[i,2]<.5 for i in p.vertices)
           and min(pts[i,0] for i in p.vertices)<=1e-6 and max(pts[i,0] for i in p.vertices)>=-1e-6]
    write('topology_audit.json',{'source_sha256':hashlib.sha256(STATIC.read_bytes()).hexdigest(),
        'vertices':len(pts),'triangles':sum(len(p.vertices)-2 for p in mesh.data.polygons),
        'test':'Mesh-edge connectivity restricted to z<0.50m, below the 0.56m crotch',
        'left_right_lower_legs_connected_below_crotch':right in prev,'path_vertices':path,'path_positions':pts[path].tolist(),
        'faces_touching_or_crossing_center_below_0_5m':faces,
        'source_uv_sha256':hashlib.sha256(np.array([tuple(v.uv) for v in mesh.data.uv_layers.active.data],dtype='<f4').tobytes()).hexdigest(),
        'conclusion':'Lower pant legs are topologically joined; separate the inner-leg mesh before final weight painting. Original mesh was not changed.' if right in prev else 'No below-crotch connection detected.'})
    print('TOPOLOGY',right in prev,'path',len(path),'central faces',len(faces))

def verify_exports():
    from mathutils.kdtree import KDTree
    bpy.ops.wm.open_mainfile(filepath=str(OUT/'prepared.blend'))
    original=next(o for o in bpy.context.scene.objects if o.type=='MESH');rest=mesh_points(original)
    edges=np.array([list(e.vertices) for e in original.data.edges]);length=np.linalg.norm(rest[edges[:,0]]-rest[edges[:,1]],axis=1)
    valid=length>np.ptp(rest[:,2])*.002
    expected=json.loads((OUT/'walking_comparison.json').read_text())
    results=[]
    for tag in ['raw','post','clean','clean_eps','local']:
        bpy.ops.wm.read_factory_settings(use_empty=True);bpy.context.scene.render.fps=30
        bpy.ops.import_scene.gltf(filepath=str(OUT/('Walking_'+tag+'.glb')))
        arm=next(o for o in bpy.context.scene.objects if o.type=='ARMATURE')
        mesh=next(o for o in bpy.context.scene.objects if o.type=='MESH');points=mesh_points(mesh)
        tree=KDTree(len(points))
        for i,p in enumerate(points):tree.insert(Vector(p),i)
        tree.balance();mapping=[];dist=[]
        for p in rest:
            _,index,distance=tree.find(Vector(p));mapping.append(index);dist.append(distance)
        assert max(dist)<1e-5,(tag,'rest mismatch',max(dist))
        actual=0;first,last=map(int,arm.animation_data.action.frame_range)
        for frame in range(first,last+1):
            bpy.context.scene.frame_set(frame);bpy.context.view_layer.update()
            obj=mesh.evaluated_get(bpy.context.evaluated_depsgraph_get());geometry=obj.to_mesh()
            deformed=np.array([obj.matrix_world@v.co for v in geometry.vertices])[mapping];obj.to_mesh_clear()
            ratios=np.linalg.norm(deformed[edges[:,0]]-deformed[edges[:,1]],axis=1)/np.maximum(length,1e-10)
            actual=max(actual,float(ratios[valid].max()))
        diff=abs(actual-expected[tag]['maxStretch']['all'])
        assert diff<.003,(tag,actual,diff)
        assert [first,last]==[1,32]
        results.append({'variant':tag,'frames':[1,32],'maxStretch':actual,'difference':diff})
    write('export_verification.json',{'passed':True,'results':results})

if __name__=='__main__':
    args=sys.argv[sys.argv.index('--')+1:]
    if args[0]=='inspect':inspect()
    elif args[0]=='prepare':prepare()
    elif args[0]=='process':process()
    elif args[0]=='topology':topology()
    elif args[0]=='verify_exports':verify_exports()
    else:raise SystemExit('Unknown stage')