Made of little things.

"""Prepare BodyParts3D 4.0 named skeletal meshes for the bundled posing rig.
Run Blender in background/factory mode with --python this_file -- SOURCE_OBJ_DIR.
Source: official partof_BP3D_4.0_obj_99.zip, extracted. Only FJ3152..FJ3395
are used (bones, costal cartilage and intervertebral discs). No user scene touched.
"""
from pathlib import Path
import sys, json, struct, re, hashlib
import numpy as np
import bpy
from mathutils import Vector
sys.path.insert(0,str(Path(__file__).resolve().parent))
from anatomy_fitting import correct_skeleton_part, assert_stable_topology
ROOT = Path(__file__).resolve().parents[1]
SOURCE = Path(sys.argv[sys.argv.index('--') + 1])
raw = (ROOT/'assets/skin/body-skin.bin').read_bytes()
size = struct.unpack('<I',raw[8:12])[0]; skin = json.loads(raw[12:12+size]); start = 12+size
spec=skin['arrays']['heads4']; heads=np.frombuffer(raw,dtype='<f4',count=spec['length'],offset=start+spec['offset']).reshape(-1,3)
rig={b['name']:i for i,b in enumerate(skin['bones'])}
meshes={}
for path in sorted(SOURCE.glob('FJ*.obj')):
    if not path.stem[2:].isdigit() or not 3152<=int(path.stem[2:])<=3395: continue
    source=path.read_text(); vertices=[]; faces=[]
    for line in source.splitlines():
        if line.startswith('v '): vertices.append([float(v) for v in line.split()[1:4]])
        elif line.startswith('f '):
            f=[int(v.split('/')[0])-1 for v in line.split()[1:]]
            faces.extend([[f[0],f[i],f[i+1]] for i in range(1,len(f)-1)])
    name=re.search(r'# English name : (.+)',source).group(1)
    fma=re.search(r'# Concept ID : (.+)',source).group(1)
    p=np.array(vertices,dtype=float)/1000
    # Source millimeters: X patient left, Z up, -Y anterior. Feet at about -80 mm.
    p=p[:,[0,2,1]];p[:,1]+=.08;p[:,2]=-p[:,2]-.075
    meshes[name.lower()]={'id':fma,'label':name,'source':path.name,'p':p,'f':faces}
assert len(meshes)==238, f'Expected all 238 source meshes, got {len(meshes)}'

def endpoint(name,top=True):
    p=meshes[name]['p'];y=p[:,1]; lo,hi=np.quantile(y,[.025,.975]);band=min(.015,(hi-lo)*.18);cut=hi-band if top else lo+band
    points=p[y>=cut] if top else p[y<=cut]
    return points.mean(axis=0)

def segment(p,a,b,target_a,target_b):
    old=b-a;new=target_b-target_a;axis=old/np.linalg.norm(old);ratio=np.linalg.norm(new)/np.linalg.norm(old)
    stretched=(p-a)+(np.outer((p-a)@axis,axis)*(ratio-1))
    q=Vector(old).rotation_difference(Vector(new));rot=np.array(q.to_matrix())
    return stretched@rot.T+target_a

# Source hand frame uses independent knuckle landmarks to retain palm width/roll.
source_heads=heads.copy();parents=[b['parent'] for b in skin['bones']]
for side,suffix in [('left','l'),('right','r')]:
    source_heads[rig[f'lowerarm_{suffix}']]=endpoint(f'{side} ulna')
    source_heads[rig[f'hand_{suffix}']]=endpoint(f'{side} ulna',False)
    for digit,label in [('index','index finger'),('middle','middle finger'),('pinky','little finger')]:
        source_heads[rig[f'{digit}_01_{suffix}']]=endpoint(f'proximal phalanx of {side} {label}')

# One transform per anatomical segment keeps paired bones correctly related.
segments={}
for side,suffix in [('left','l'),('right','r')]:
    for name,driver,tail in [('humerus','upperarm','lowerarm'),('ulna','lowerarm','hand'),('femur','thigh','calf'),('tibia','calf','foot')]:
        a=endpoint(f'{side} {name}');b=endpoint(f'{side} {name}',False)
        segments[(side,name)]=(a,b,heads[rig[f'{driver}_{suffix}']],heads[rig[f'{tail}_{suffix}']])

# Fit each phalanx to the corresponding rig segment. Metacarpals connect the
# wrist to each finger root. Tiny carpals use a common wrist transform.
orders={'first':'thumb','second':'index','third':'middle','fourth':'ring','fifth':'pinky'}
foot_words=['toe','metatarsal','talus','calcaneus','cuneiform','cuboid','navicular bone']
carpal_words=['scaphoid','lunate','triquetral','pisiform','trapezium','trapezoid','capitate','hamate']
parts=[];blob=bytearray();catalog=[]
for name,m in sorted(meshes.items(),key=lambda item:int(item[1]['id'][3:])):
    p=m['p'].copy();side='left' if 'left' in name else 'right' if 'right' in name else '';suffix='l' if side=='left' else 'r'
    bone='spine_03';tail=None;part='torso'
    if side and any(w in name for w in ['humerus','radius','ulna','femur','tibia','fibula']):
        key='humerus' if 'humerus' in name else 'ulna' if ('radius' in name or 'ulna' in name) else 'femur' if 'femur' in name else 'tibia'
        driver,nextbone,region={'humerus':('upperarm','lowerarm','upperarm'),'ulna':('lowerarm','hand','forearm'),'femur':('thigh','calf','thigh'),'tibia':('calf','foot','calf')}[key]
        p=segment(p,*segments[(side,key)]);bone=f'{driver}_{suffix}';tail=f'{nextbone}_{suffix}';part=f'{side}-{region}'
    elif side and 'patella' in name:
        p+=heads[rig[f'calf_{suffix}']]+np.array([0,.007,.035])-p.mean(axis=0);bone=f'calf_{suffix}';part=f'{side}-knee'
    elif side and any(w in name for w in foot_words):
        bone=f'foot_{suffix}';part=f'{side}-foot'
        # Keep the original foot proportions; translate its talar center to ankle.
        p+=heads[rig[bone]]-meshes[f'{side} talus']['p'].mean(axis=0)
    elif side and ('finger' in name or 'thumb' in name or 'metacarpal' in name):
        # Keep this legacy pre-decimation fit stable. Correct the ambiguous
        # "middle phalanx" driver below, after reduction, so triangle IDs survive.
        digit=next((d for d in ['thumb','index','middle','ring','little'] if d in name),None)
        if not digit: digit=next(orders[o] for o in orders if o in name)
        if digit=='little': digit='pinky'
        if 'metacarpal' in name:
            bone=f'hand_{suffix}';tail=f'{digit}_01_{suffix}'
        else:
            segment_number=1 if 'proximal' in name else 2 if 'middle phalanx' in name else 3
            # The thumb's metacarpal starts at rig thumb_01; its two phalanges at 02/03.
            if digit=='thumb' and segment_number==1:segment_number=2
            bone=f'{digit}_{segment_number:02}_{suffix}';tail=f'{digit}_{segment_number+1:02}_{suffix}' if segment_number<3 else None
        a=endpoint(name);b=endpoint(name,False);target_a=heads[rig[bone]]
        if 'metacarpal' in name and digit=='thumb':bone=f'thumb_01_{suffix}';tail=f'thumb_02_{suffix}';target_a=heads[rig[bone]]
        target_b=heads[rig[tail]] if tail else target_a+(target_a-heads[skin['bones'][rig[bone]]['parent']])*.72
        p=segment(p,a,b,target_a,target_b);part=f'{side}-hand'
    elif side and any(w in name for w in carpal_words):
        bone=f'hand_{suffix}';part=f'{side}-hand';wrist=segments[(side,'ulna')][1]
        p=segment(p,wrist,wrist+np.array([0,-.08,0]),heads[rig[bone]],heads[rig[f'middle_01_{suffix}']])
    elif side and ('clavicle' in name or 'scapula' in name):
        bone=f'clavicle_{suffix}';part=f'{side}-shoulder'
        p+=heads[rig[f'upperarm_{suffix}']]-segments[(side,'humerus')][0]
    elif 'hip bone' in name or name in ['sacrum','coccyx']:
        bone='pelvis';part='pelvis'
    elif 'cervical' in name or name in ['atlas','axis']:
        bone='neck_01';part='neck'
    elif 'lumbar' in name or name=='intervertebral disk':
        bone='spine_01'
    elif 'thoracic' in name or 'rib' in name or 'costal' in name or 'sternum' in name or name in ['manubrium','xiphoid process']:
        bone='spine_03'
    else:
        bone='head';part='head'
    p-=heads[rig[bone]]
    assert np.isfinite(p).all(), f'Non-finite geometry: {name}'
    geometry=bpy.data.meshes.new(m['id']);geometry.from_pydata(p.tolist(),[],m['f']);geometry.update()
    obj=bpy.data.objects.new(m['id'],geometry);bpy.context.scene.collection.objects.link(obj)
    bpy.context.view_layer.objects.active=obj;obj.select_set(True)
    budget=2400 if part=='head' else 1600 if any(w in name for w in ['femur','humerus','tibia']) else 700
    if len(geometry.polygons)>budget:
        mod=obj.modifiers.new('Web reduction','DECIMATE');mod.ratio=budget/len(geometry.polygons)
        bpy.ops.object.modifier_apply(modifier=mod.name)
    geometry.calc_loop_triangles()
    points,corrected_tail=correct_skeleton_part(np.array([v.co[:] for v in geometry.vertices]),name,part,rig[bone],rig[tail] if tail else -1,heads,source_heads,rig,parents,endpoint)
    digit_match=re.search(r'middle phalanx of (left|right) (ring|little) finger',name)
    if digit_match:
        digit='ring' if digit_match[2]=='ring' else 'pinky'
        corrected_bone=f'{digit}_02_{suffix}';next_bone=f'{digit}_03_{suffix}'
        # Refit the already-reduced vertices to the intended segment without
        # rebuilding faces or changing barycentric annotation addresses.
        points=segment(points+heads[rig[bone]],heads[rig[bone]],heads[rig[tail]],heads[rig[corrected_bone]],heads[rig[next_bone]])-heads[rig[corrected_bone]]
        bone=corrected_bone;corrected_tail=rig[next_bone]
    positions=np.array(points,dtype='<f4').tobytes()
    faces=np.array([t.vertices[:] for t in geometry.loop_triangles],dtype='<u2').tobytes()
    assert len(geometry.vertices)<65536
    while len(blob)%4:blob.append(0)
    po=len(blob);blob+=positions
    io=len(blob);blob+=faces
    info={'id':m['id'],'label':m['label'],'part':part,'bone':rig[bone],'tail':corrected_tail,'vertices':len(geometry.vertices),'triangles':len(geometry.loop_triangles),'positions':po,'indices':io}
    parts.append(info);catalog.append({'id':m['id'],'label':m['label'],'triangles':info['triangles'],'part':part,'source':m['source']})
    bpy.data.objects.remove(obj,do_unlink=True);bpy.data.meshes.remove(geometry)
meta=json.dumps({'id':'bp3d-skeleton-1','parts':parts},separators=(',',':')).encode()
meta+=b' '*((-len(meta))%4)
output=b'BJBONE01'+struct.pack('<I',len(meta))+meta+blob
assert_stable_topology(ROOT/'assets/skeleton/skeleton.bin',parts,blob)
(ROOT/'assets/skeleton/skeleton.bin').write_bytes(output)
(ROOT/'src/skeleton-catalog.ts').write_text('// Generated by scripts/prepare-skeleton.py. Stable topology: bump model ID when rebuilding.\nexport const skeletonCatalog = '+json.dumps(catalog,separators=(',',':'))+';\n')
provenance={'id':'bp3d-skeleton-1','source':'BodyParts3D 4.0, PART-OF tree, 99% reduction archive','url':'https://dbarchive.biosciencedbc.jp/data/bodyparts3d/LATEST/partof_BP3D_4.0_obj_99.zip','attribution':'BodyParts3D, © The Database Center for Life Science licensed under CC Attribution 4.0 International','license':'CC-BY-4.0','licenseSource':'https://dbarchive.biosciencedbc.jp/en/bodyparts3d/lic.html','licenseUpdated':'2025-02-27','preparationBlenderVersion':bpy.app.version_string,'archiveSha256':'9fbc713fffeee924a5a657d9813d84d7eb957bded63adb854931dd5e3eb61c97','fittingRevision':4,'modifications':'Coordinate conversion, geometry decimation, topology-preserving regional fit; shared palm frame and distal forearm roll; corrected ring/little middle-phalanx drivers; binary packaging. Bones, discs and costal cartilage remain separately named.','structures':len(parts),'triangles':sum(p['triangles'] for p in parts),'bytes':len(output),'sha256':hashlib.sha256(output).hexdigest(),'sourceFiles':{m['source']:hashlib.sha256((SOURCE/m['source']).read_bytes()).hexdigest() for m in meshes.values()}}
(ROOT/'assets/skeleton/provenance.json').write_text(json.dumps(provenance,indent=2)+'\n')
print('SKELETON', {k:v for k,v in provenance.items() if k!='sourceFiles'})