import struct, math, numpy as np
from PIL import Image
data=open("content/VIDEO/GEO/BLX_COP.BGF","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))
V=np.array(verts)
# color each face by its normal: R=|nx| G=|ny| B=|nz|.  Render from centroid looking -Z.
def render(eye,look,W=400,H=300,fov=90):
    eye=np.array(eye,float);fwd=np.array(look,float)-eye;fwd/=np.linalg.norm(fwd)
    up=np.array([0,1,0.]);right=np.cross(fwd,up);right/=np.linalg.norm(right);upv=np.cross(right,fwd)
    f=(H/2)/math.tan(math.radians(fov)/2);img=np.zeros((H,W,3));zb=np.full((H,W),1e18)
    hist={'+z':0,'-z':0,'+x':0,'-x':0,'+y':0,'-y':0}
    for face in faces:
        tri=V[list(face)]
        nrm=np.cross(tri[1]-tri[0],tri[2]-tri[0]);nl=np.linalg.norm(nrm)
        if nl<1e-9:continue
        nrm/=nl
        col=[abs(nrm[0]),abs(nrm[1]),abs(nrm[2])]  # R=x G=y B=z
        ax=np.argmax(np.abs(nrm)); s='+' if nrm[ax]>0 else '-'; hist[s+'xyz'[ax]]+=1
        cam=[(np.dot(pt-eye,right),np.dot(pt-eye,upv),np.dot(pt-eye,fwd)) for pt in tri]
        if any(cz<=0.03 for _,_,cz in cam):continue
        pts=[(W/2+f*cx/cz,H/2-f*cy/cz,cz) for cx,cy,cz in cam]
        xs=[p[0] for p in pts];ys=[p[1] for p in pts]
        mnx=max(int(min(xs)),0);mxx=min(int(max(xs))+1,W-1);mny=max(int(min(ys)),0);mxy=min(int(max(ys))+1,H-1)
        if mnx>=mxx or mny>=mxy:continue
        (x0,y0,z0),(x1,y1,z1),(x2,y2,z2)=pts
        den=((y1-y2)*(x0-x2)+(x2-x1)*(y0-y2))
        if abs(den)<1e-9:continue
        ys2,xs2=np.mgrid[mny:mxy+1,mnx:mxx+1];xf=xs2+.5;yf=ys2+.5
        l0=((y1-y2)*(xf-x2)+(x2-x1)*(yf-y2))/den;l1=((y2-y0)*(xf-x2)+(x0-x2)*(yf-y2))/den;l2=1-l0-l1
        m=(l0>=0)&(l1>=0)&(l2>=0)
        z=l0*z0+l1*z1+l2*z2
        yi=ys2[m];xi=xs2[m];zz=z[m];cl=zz<zb[yi,xi]
        zb[yi[cl],xi[cl]]=zz[cl];img[yi[cl],xi[cl]]=col
    return (np.clip(img,0,1)*255).astype(np.uint8),hist
cen=V.mean(0)
img,hist=render(cen.tolist(),(cen+[0,0,-30]).tolist())
Image.fromarray(img).save("scratchpad/cop_normals.png")
print("face-normal dominant-axis histogram:",hist)
print("saved cop_normals.png [R=|nx| G=|ny| B=|nz|], view from centroid look -Z")