CMPlanning/scripts/verify_geom/final_slice.py

import sys, json, pickle, collections
sys.path += ['/root/.local/lib/python3.13/site-packages','/usr/local/lib/python3.13/dist-packages','/usr/lib/python3/dist-packages']
import numpy as np
from build123d import import_step, Shape
from OCP.BRepAlgoAPI import BRepAlgoAPI_Section
from OCP.gp import gp_Pln, gp_Pnt, gp_Dir
from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeFace
from OCP.TopExp import TopExp_Explorer
from OCP.TopAbs import TopAbs_EDGE
solids = import_step("/mnt/user-data/uploads/97.과제/Code/3d-drawing/model/tunnel_solids.step").solids()
mp = pickle.load(open('/tmp/verify_daesan2_run/mapping.pkl','rb'))['mapping']
caps = json.load(open('/tmp/verify_daesan2_run/caps.json'))
ch = json.load(open('/tmp/verify_daesan2_run/chain.json'))
meta = json.load(open('/mnt/user-data/uploads/97.과제/Code/3d-drawing/model/tunnel_solids_meta.json'))
comps = {c['id']: c for c in meta['components'] if c['group']!='시설한계(참조선)'}
solid_comp = {v:k for k,v in mp.items()}
# 1) boundary points from lining chain
pts = []
for k in ch['lin']:
    pts.append(np.array(ch['starts'][k])); pts.append(np.array(ch['ends'][k]))
P = np.array(pts)[:, :2]
uniq = []
for p in P:
    if not any(np.linalg.norm(p-q) < 1e-3 for q in uniq): uniq.append(p)
uniq = np.array(uniq)
print("unique boundary points", len(uniq), "(expect 41)")
D = P[-1]-P[0] if False else (np.array(ch['ends'][ch['lin'][-1]])[:2]-np.array(ch['starts'][ch['lin'][0]])[:2])
proj = uniq@D/np.linalg.norm(D)
order = np.argsort(proj); B = uniq[order]
cum = np.concatenate([[0], np.cumsum(np.linalg.norm(np.diff(B,axis=0),axis=1))])
Sb = 4745 + cum
M = sorted(set([c['from_sta'] for c in comps.values()] + [5480]))
print("boundary count", len(Sb), "meta boundary count", len(M))
print("max |station_chain - meta_boundary| (m):", round(float(np.max(np.abs(Sb-np.array(M)))),3) if len(Sb)==len(M) else 'count mismatch')
# 2) station for each solid by centroid projection
def station_of(p):
    best=None
    for k in range(len(B)-1):
        a,b=B[k],B[k+1]; t=b-a; L2=t@t
        f=np.clip(((p-a)@t)/L2,0,1); q=a+f*t; d=np.linalg.norm(p-q)
        if best is None or d<best[0]: best=(d, Sb[k]+f*(Sb[k+1]-Sb[k]), k, f)
    return best
sol_info={}
for i,s in enumerate(solids):
    c=s.center(); d,st,k,f=station_of(np.array([c.X,c.Y]))
    cid=solid_comp[i]; sol_info[i]=dict(cid=cid, group=comps[cid]['group'], st=st, lat=d)
# 3) pair within group by sorted mid
assign={}
groups=collections.defaultdict(list)
for i,v in sol_info.items(): groups[v['group']].append(i)
maxres=0
for g,idxs in groups.items():
    cs=[c for c in comps.values() if c['group']==g]
    idxs_s=sorted(idxs, key=lambda i: sol_info[i]['st'])
    cs_s=sorted(cs, key=lambda c: (c['from_sta']+c['to_sta'])/2)
    if len(idxs_s)!=len(cs_s): print("COUNT MISMATCH", g, len(idxs_s), len(cs_s))
    for i,c in zip(idxs_s, cs_s):
        mid=(c['from_sta']+c['to_sta'])/2
        res=abs(sol_info[i]['st']-mid); maxres=max(maxres,res)
        assign[i]=dict(id=c['id'], group=g, frm=c['from_sta'], to=c['to_sta'], mid_chain=sol_info[i]['st'], lat=sol_info[i]['lat'])
print("max residual solid-centre station vs comp mid (m):", round(maxres,3))
# lateral offset sanity
print("max lateral offset of solid centroid from chain (m):", round(max(v['lat'] for v in assign.values()),3))
# 4) slice
def slice_solid(idx, Q, T):
    nrm=np.array([T[0],T[1],0.0])
    pl=gp_Pln(gp_Pnt(float(Q[0]),float(Q[1]),0.0), gp_Dir(float(nrm[0]),float(nrm[1]),0.0))
    face=BRepBuilderAPI_MakeFace(pl,-80,80,-80,80).Face()
    sec=BRepAlgoAPI_Section(solids[idx].wrapped, face); sec.Build()
    polys=[]
    exp=TopExp_Explorer(sec.Shape(), TopAbs_EDGE)
    while exp.More():
        e=Shape.cast(exp.Current()); L=e.length; N=max(2,int(L/0.02))
        polys.append([[p.X,p.Y,p.Z] for p in (e.position_at(j/N) for j in range(N+1))])
        exp.Next()
    return polys
def plane_at(s):
    k=int(np.searchsorted(Sb, s, side='right')-1); k=min(max(k,0),len(B)-2)
    f=(s-Sb[k])/(Sb[k+1]-Sb[k])
    Q=(1-f)*B[k]+f*B[k+1]; T=(B[k+1]-B[k]); T=T/np.linalg.norm(T)
    return Q, T
out={}
for s in [4745.0005, 5100.0, 5469.9995]:
    Q,T=plane_at(s); eh=np.array([-T[1],T[0]])
    print(f"== station {s}: Q=({Q[0]:.3f},{Q[1]:.3f}) T=({T[0]:.5f},{T[1]:.5f})")
    out[str(s)]={'Q':Q.tolist(),'T':T.tolist(),'eh':eh.tolist(),'solids':{}}
    cnt=collections.Counter()
    for i,a in assign.items():
        if not (a['frm'] <= s <= a['to']): continue
        polys=slice_solid(i,Q,T)
        uv=[]
        for pl_ in polys:
            arr=np.array(pl_); u=(arr[:,:2]-Q)@eh; uv.append(np.stack([u,arr[:,2]],1).tolist())
        out[str(s)]['solids'][str(i)]={'id':a['id'],'group':a['group'],'polys':uv}
        cnt[a['group']]+=1
        allp=np.vstack([np.array(x) for x in uv]) if uv else np.zeros((0,2))
        if len(allp): print(f"   {a['group']:<12} {a['id']:<16} nloops={len(uv)} u[{allp[:,0].min():.3f},{allp[:,0].max():.3f}] v[{allp[:,1].min():.3f},{allp[:,1].max():.3f}]")
        else: print("   EMPTY", a['group'], a['id'])
    print("   counts", dict(cnt))
json.dump(out, open('/tmp/verify_daesan2_run/sections_final.json','w'))
json.dump({str(i):a for i,a in assign.items()}, open('/tmp/verify_daesan2_run/assign.json','w'), ensure_ascii=False)
print("saved")