"""Evaluate original CAD surface derivatives; no scan-derived curvature labels."""
from ac28_audit import *
from OCP.STEPControl import STEPControl_Reader
from OCP.BRepMesh import BRepMesh_IncrementalMesh
from OCP.TopExp import TopExp_Explorer
from OCP.TopAbs import TopAbs_FACE
from OCP.TopoDS import TopoDS
from OCP.BRep import BRep_Tool
from OCP.TopLoc import TopLoc_Location
from OCP.BRepAdaptor import BRepAdaptor_Surface
from OCP.BRepLProp import BRepLProp_SLProps
from scipy.spatial import cKDTree
import shutil

def main():
    r=STEPControl_Reader();r.ReadFile(str(next(DATA.glob('*.stp'))));r.TransferRoots();shape=r.OneShape()
    BRepMesh_IncrementalMesh(shape,.3,False,.2,True)
    ex=TopExp_Explorer(shape,TopAbs_FACE);rows=[];face_id=0;failed=0
    while ex.More():
        face=TopoDS.Face(ex.Current());tri=BRep_Tool.Triangulation_s(face,TopLoc_Location());surf=BRepAdaptor_Surface(face,True)
        if tri is not None and tri.HasUVNodes():
            for i in range(1,tri.NbTriangles()+1,max(1,tri.NbTriangles()//100)):
                ids=tri.Triangle(i).Get();uv=[tri.UVNode(j) for j in ids];u=sum(v.X() for v in uv)/3;v=sum(v.Y() for v in uv)/3
                try:
                    props=BRepLProp_SLProps(surf,u,v,2,1e-7)
                    if not props.IsCurvatureDefined():failed+=1;continue
                    p=props.Value();n=props.Normal();k1=props.MaxCurvature();k2=props.MinCurvature()
                    row=[p.X(),p.Y(),p.Z(),n.X(),n.Y(),n.Z(),k1,k2,face_id]
                    if np.all(np.isfinite(row)):rows.append(row)
                except Exception:failed+=1
        face_id+=1;ex.Next()
    a=np.asarray(rows);np.savez_compressed(OUT/'cad_surface_curvature.npz',points_mm=a[:,:3],parameter_normals=a[:,3:6],k1_per_mm=a[:,6],k2_per_mm=a[:,7],face_id=a[:,8].astype('int64'))
    g=json.loads((OUT/'design_graph.json').read_text());p=np.array(g['positions']).reshape(-1,3);n=np.array(g['normals']).reshape(-1,3)
    d,idx=cKDTree(a[:,:3]).query(p);valid=(d<3)&(np.abs(np.sum(n*a[idx,3:6],axis=1))>.8)
    g['curvature_max_abs_per_mm']=[float(max(abs(a[j,6]),abs(a[j,7]))) if ok else None for j,ok in zip(idx,valid)]
    g['curvature_mapping_distance_mm']=d.round(3).tolist()
    g['curvature_note']='Original STEP analytic surface derivatives sampled at tessellation UV triangle centroids; nearest sample mapped to simplified graph only within 3 mm and |normal dot| > 0.8. Sharp face boundaries do not have a unique smooth curvature.'
    (OUT/'design_graph.json').write_text(json.dumps(g,separators=(',',':')))
    graph=dict(np.load(OUT/'cad_graph.npz'));graph.update(k1_per_mm=a[idx,6].astype('float32'),k2_per_mm=a[idx,7].astype('float32'),curvature_valid=valid,curvature_mapping_distance_mm=d.astype('float32'))
    np.savez_compressed(OUT/'cad_graph.npz',**graph)
    for name in ['design_graph.json','cad_graph.npz','cad_surface_curvature.npz']:shutil.copy2(OUT/name,DATA/'viewer/explorer_data'/name)
    print('surface samples',len(a),'failed',failed,'mapped graph nodes',valid.sum(),'/',len(p),'k quantiles',np.quantile(np.max(np.abs(a[:,6:8]),axis=1),[.5,.9,.99]),flush=True)

if __name__=='__main__':main()