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)
# boundary-edge count => holes?
from collections import Counter
ec=Counter()
for a,b,c in faces:
    for e in ((a,b),(b,c),(c,a)):
        ec[tuple(sorted(e))]+=1
bound=sum(1 for k,v in ec.items() if v==1)
print(f"faces={len(faces)} edges={len(ec)} boundary_edges(used once)={bound}  (many boundary edges => open holes/lattice)")

def render(eye,look,cull,fov=95,W=380,H=290):
    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)
    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; fc=tri.mean(0)
        facing=np.dot(nrm,fc-eye)   # >0 => normal points away from eye (backface)
        if cull=='front' and facing>0: continue   # keep faces facing eye
        if cull=='back' and facing<0: continue
        col=[0.85,0.85,0.85]
        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)
cen=V.mean(0);eye=cen.tolist();look=(cen+[0,0,-30]).tolist()
row=[render(eye,look,'none'),render(eye,look,'front'),render(eye,look,'back')]
Image.fromarray(np.concatenate(row,axis=1)).save("scratchpad/cop_cull.png")
print("saved cop_cull.png [no-cull | keep front-facing | keep back-facing]  black=open")