import sys; sys.path.insert(0, "./pylib")
sys.stdout.reconfigure(encoding='utf-8')
import rhino3dm, numpy as np
m = rhino3dm.File3dm.Read('m.3dm')
byidx = {l.Index: l.FullPath for l in m.Layers}
def meshpts(g, n=1):
    V=[];F=[];off=0
    for fi in range(len(g.Faces)):
        k=g.Faces[fi].GetMesh(rhino3dm.MeshType.Any)
        if k is None: continue
        v=np.array([[p.X,p.Y,p.Z] for p in k.Vertices]); V.append(v)
        for i in range(len(k.Faces)):
            f=k.Faces[i]; F.append((f[0]+off,f[1]+off,f[2]+off))
            if f[2]!=f[3]: F.append((f[0]+off,f[2]+off,f[3]+off))
        off+=len(v)
    return np.vstack(V), np.array(F)
cop=None; bars=[]
for o in m.Objects:
    p=byidx.get(o.Attributes.LayerIndex)
    if p=='P1::코핑': cop=o.Geometry
    if p=='P1::철근::코핑::C1': bars.append((o.Attributes.Name,o.Geometry))
V,F=meshpts(cop)
print('coping verts',len(V),'faces',len(F))
# sample points on coping triangles
A=V[F[:,0]];B=V[F[:,1]];C=V[F[:,2]]
N=np.cross(B-A,C-A); area=np.linalg.norm(N,axis=1)/2; N=N/np.maximum(np.linalg.norm(N,axis=1)[:,None],1e-12)
rng=np.random.default_rng(1)
pts=[];nrm=[]
for k in range(len(F)):
    n=max(1,int(area[k]/400))   # ~20mm grid
    r1=rng.random(n);r2=rng.random(n); s=np.sqrt(r1)
    P=A[k]*(1-s)[:,None]+B[k]*(s*(1-r2))[:,None]+C[k]*(s*r2)[:,None]
    pts.append(P);nrm.append(np.repeat(N[k][None],n,0))
pts=np.vstack(pts);nrm=np.vstack(nrm)
pts=np.vstack([pts,V]);nrm=np.vstack([nrm,np.zeros_like(V)])
print('sample pts',len(pts))
print('name | min clear cover | nearest-face normal | pt z of bar | #verts')
for name,g in bars:
    bv,_=meshpts(g)
    # subsample
    bv=bv[::2]
    best=1e9;bi=None
    for s in range(0,len(bv),400):
        c=bv[s:s+400]
        d=np.linalg.norm(c[:,None,:]-pts[None,:,:],axis=2)
        j=d.argmin(axis=1); dm=d[np.arange(len(c)),j]
        t=dm.argmin()
        if dm[t]<best: best=dm[t];bi=(c[t],pts[j[t]],nrm[j[t]])
    print(name,'|',round(best,1),'|',np.round(bi[2],2),'| bar',np.round(bi[0],0),'cop',np.round(bi[1],0))