import struct, glob, os, numpy as np
from collections import Counter
def load(path):
    data=open(path,"rb").read()
    verts=[];faces=[];cb=[0]
    def addpoly(idx,b):
        for k in range(1,len(idx)-1): faces.append((b+idx[0],b+idx[k],b+idx[k+1]))
    def parse(p,end):
        while p+3<=end:
            tw=struct.unpack_from("<H",data,p)[0];p+=2
            cid=tw&0x2fff;lw=4 if(tw&0x8000)else(2 if(tw&0x4000)else 1)
            if p+lw>end:return
            ln=data[p] if lw==1 else(struct.unpack_from("<H",data,p)[0] if lw==2 else struct.unpack_from("<I",data,p)[0])
            p+=lw
            if p+ln>end:return
            if cid==0x88:
                cb[0]=len(verts)
                for i in range(ln//20):
                    o=p+i*20;x,y,z=struct.unpack_from("<fff",data,o);verts.append((x,y,z))
            elif cid==0x47:
                n=ln//4;idx=[struct.unpack_from("<I",data,p+i*4)[0] for i in range(n)]
                for f in range(0,(n//3)*3,3): faces.append((cb[0]+idx[f],cb[0]+idx[f+1],cb[0]+idx[f+2]))
            elif cid==0x4d:
                if ln>=1:
                    ppf=data[p];n=(ln-1)//4
                    idx=[struct.unpack_from("<I",data,p+1+i*4)[0] for i in range(n)]
                    for f in range(0,(n//ppf)*ppf,ppf): addpoly([idx[f+j] for j in range(ppf)],cb[0])
            if cid in {0x3,0x70,0x40,0x41,0x42,0x46,0x48,0x10,0x20,0x30}: parse(p,p+ln)
            p+=ln
            if cid==0x5:break
    parse(8,len(data))
    return np.array(verts),faces
for path in sorted(glob.glob("content/VIDEO/GEO/*_COP.BGF")):
    V,faces=load(path)
    if len(V)==0: print(os.path.basename(path),'-- no verts parsed'); continue
    nv=len(V); faces=[(a,b,c) for (a,b,c) in faces if a<nv and b<nv and c<nv]
    # boundary edges
    ec=Counter()
    for a,b,c in faces:
        for edg in ((a,b),(b,c),(c,a)): ec[tuple(sorted(edg))]+=1
    bound=sum(1 for v in ec.values() if v==1)
    # normal dominant-axis area (weighted by face area)
    axarea={'+x':0.,'-x':0.,'+y':0.,'-y':0.,'+z':0.,'-z':0.}
    for a,b,c in faces:
        n=np.cross(V[b]-V[a],V[c]-V[a]); area=np.linalg.norm(n)
        if area<1e-9: continue
        nn=n/area; ax=int(np.argmax(np.abs(nn))); s='+' if nn[ax]>0 else '-'
        axarea[s+'xyz'[ax]]+=area
    tot=sum(axarea.values()) or 1
    dom=max(axarea,key=axarea.get)
    name=os.path.basename(path)
    frac={k:round(v/tot,2) for k,v in axarea.items()}
    print(f"{name:14s} verts={len(V):4d} faces={len(faces):4d} boundary={bound:4d} | area-by-axis {frac} | DOMINANT={dom}")