"""Export auditable comparison fields and a topology-based graph scaffold."""
from ac28_audit import *
from register_ac28 import transform,fit
from scipy.spatial import cKDTree
import shutil, struct
from PIL import Image

WEB=DATA/'viewer'

def closest(scene,q):
    r=scene.compute_closest_points(o3d.core.Tensor(q.astype('float32')))
    cp=r['points'].numpy();normal=r['primitive_normals'].numpy();delta=q-cp
    return cp,np.linalg.norm(delta,axis=1),np.einsum('ij,ij->i',delta,normal)

def main():
    web=WEB/'explorer_data';web.mkdir(exist_ok=True)
    mesh=o3d.io.read_triangle_mesh(str(OUT/'design_fine.stl'));mesh.remove_duplicated_vertices();mesh.compute_triangle_normals()
    scene=o3d.t.geometry.RaycastingScene();scene.add_triangles(o3d.t.geometry.TriangleMesh.from_legacy(mesh))
    comparison=json.loads((OUT/'registration_comparison.json').read_text())
    graph=np.load(OUT/'cad_graph.npz');gp=graph['pos_mm']
    baseline=json.loads((ROOT/'analysis/baseline_transforms.json').read_text())
    variants={
        '16h_mirror':('16h',np.array(comparison['reflection_diagnostic']['transform']),True),
        '16h_rigid':('16h',np.array(comparison['rigid']['transform']),True),
        '16h_previous':('16h',np.array(baseline['16h_previous']),False),
        '18h01':('18h01',np.array(baseline['18h01']),False)}
    metadata={'variants':{},'design':json.loads((OUT/'design_metadata.json').read_text()),'metric':'Euclidean distance to nearest CAD triangle; full 2 mm voxel sample; no distance exclusions','seed':42}
    scans={}
    base=next(p for p in DATA.glob('1. PLY*') if p.is_dir())
    for name in ['16h','18h01']:
        folder=next(base.glob('*16*' if name=='16h' else '*18*'))
        c=o3d.io.read_point_cloud(str(folder/'all_merge_pointcloud_full.ply'))
        # Preserve actual measured coordinates in the display sample (no voxel averaging).
        raw=np.asarray(c.points)[::12];rgb=np.asarray(c.colors)[::12] if c.has_colors() else np.full_like(raw,.7)
        _,idx=np.unique(np.floor(raw/.8).astype(np.int32),axis=0,return_index=True)
        scans[name]=(raw[idx].copy(),rgb[idx].copy(),len(c.points))
        del c
        im=Image.open(folder/'Merged_Image.png');dims=im.size
        im.thumbnail((800,1800));im.save(web/f'image_{name}.png')
        metadata.setdefault('images',{})[name]={'size':dims,'dent_pixel':[1223,8110] if name=='16h' else [1655,1754],'label':'찍힘 (dent)','coordinate_origin':'CSV convention unconfirmed'}
    for key,(name,t,refine) in variants.items():
        p=np.load(OUT/f'{name}_all_merge_pointcloud_full.npz')['points']
        if refine:
            for _ in range(35):
                q=transform(p,t);cp,d,_=closest(scene,q);mask=d<12
                delta=fit(q[mask],cp[mask]);t=delta@t
                if np.linalg.norm(delta[:3,3])<.001:break
        q=transform(p,t);cp,d,signed=closest(scene,q)
        raw,rgb,n=scans[name];v=transform(raw,t);vcp,vd,vs=closest(scene,v)
        # Local plane RMS at an explicit physical scale; not standardized surface Ra.
        tree=cKDTree(q);dd,ii=tree.query(q,k=20,workers=2)
        near=q[ii];centered=near-near.mean(1,keepdims=True)
        cov=np.einsum('nki,nkj->nij',centered,centered)/20
        rms=np.sqrt(np.maximum(np.linalg.eigvalsh(cov)[:,0],0));rms[dd[:,-1]>8]=np.nan
        nearest=tree.query(v,workers=2)[1];vrms=rms[nearest]
        # Rows = x y z Euclidean_mm local_normal_projection_mm plane_RMS_mm r g b.
        data=np.column_stack([v,vd,vs,vrms,rgb]).astype('<f4')
        with (web/f'{key}.bin').open('wb') as f:f.write(struct.pack('<I',len(v)));f.write(data.tobytes())
        gdist=tree.query(gp,workers=2)[0].astype('float32');gdist.tofile(web/f'{key}_coverage.bin')
        region={}
        for region_name,mask in [('long_arm',q[:,2]>700),('bend',(q[:,2]>600)&(q[:,2]<=700)),('leg',q[:,2]<=600)]:
            region[region_name]={'n':int(mask.sum()),'median_mm':float(np.median(d[mask])) if mask.any() else None,'p90_mm':float(np.quantile(d[mask],.9)) if mask.any() else None}
        sensors=[]
        for path in sorted(OUT.glob(f'{name}_merged_pointcloud_*.npz')):
            sp=transform(np.load(path)['points'],t);_,sd,_=closest(scene,sp)
            sensors.append({'sensor':path.stem.split('merged_pointcloud_')[-1],'median_mm':float(np.median(sd)),'p90_mm':float(np.quantile(sd,.9)),'within5':float(np.mean(sd<5))})
        row={'batch':name,'n_raw':n,'n_display':len(v),'n_evaluation':len(q),'transform':t.tolist(),'determinant':float(np.linalg.det(t[:3,:3])),'median_mm':float(np.median(d)),'p90_mm':float(np.quantile(d,.9)),'p95_mm':float(np.quantile(d,.95)),'within2':float(np.mean(d<2)),'within5':float(np.mean(d<5)),'within10':float(np.mean(d<10)),'outside20':float(np.mean(d>20)),'regions':region,'sensors':sensors,'independently_validated':False}
        metadata['variants'][key]=row
        # Correspondence is a distance target only, not tracked physical displacement.
        np.savez_compressed(OUT/f'{key}_graph_observations.npz',cad_to_scan_distance_mm=gdist,coverage_within5=gdist<5,training_valid=np.zeros(len(gp),dtype=bool),transform_raw_to_cad=t)
        shutil.copy2(OUT/f'{key}_graph_observations.npz',web/f'{key}_graph_observations.npz')
        print(key, {k:row[k] for k in ['median_mm','p90_mm','within5','outside20']},flush=True)
    for name in ['process.json','design_graph.json','cad_graph.npz','registration_comparison.json','raw_geometry.png','rotation_vs_reflection.png']:
        shutil.copy2(OUT/name,web/name)
    (web/'metadata.json').write_text(json.dumps(metadata,indent=2,ensure_ascii=False))
    (OUT/'explorer_metadata.json').write_text(json.dumps(metadata,indent=2,ensure_ascii=False))

if __name__=='__main__':main()