#!/usr/bin/env python3 """Measure three additional sizes per CAD representation for both parts. Large image NPZs contain the exact six-view z-buffer/visibility diagnostic, not a materialized 17-channel CNN tensor. Its byte cost is calculated exactly. Large sparse voxel NPZs retain model features and masks, omitting the millions of per-triangle sampling correspondences used to build occupied cells. """ from __future__ import annotations import csv import json import sys import time from pathlib import Path import numpy as np from scipy.spatial import cKDTree HERE = Path(__file__).resolve().parent ROOT = HERE.parents[1] sys.path.insert(0, str(ROOT / "전처리")) from compare_graph_representations import prepared_cad from representations import cad_points, surface_voxels, _view_basis PARTS = {"AC28": ROOT, "ME1a_DT": ROOT / "experiments/ME1a_DT"} IMAGE_SIZES = (512, 1024, 2048) VOXEL_SIZES = {256: (3.75, (256, 32, 160)), 512: (1.875, (512, 64, 320)), 1024: (.9375, (1024, 128, 640))} POINT_SIZES = (16384, 32768, 65536) GRAPH_SIZES = (24000, 48000, 96000) FOV_MM = 980. DIRECTIONS = np.asarray([[1,0,0],[-1,0,0],[0,1,0],[0,-1,0],[0,0,1],[0,0,-1]]) def raster_geometry(cad: dict, view: np.ndarray, resolution: int, field_of_view_mm: float = FOV_MM) -> tuple[np.ndarray, np.ndarray]: """Same triangle-center z-buffer convention as representations.rasterize_view.""" p = np.asarray(cad["pos_mm"], np.float64) f = np.asarray(cad["face"], np.int64).T v = np.asarray(view, np.float64); v /= np.linalg.norm(v) u, w = _view_basis(v) pu, pw, depth = p @ u, p @ w, p @ v uc, wc = (pu.min()+pu.max())/2, (pw.min()+pw.max())/2 pitch = field_of_view_mm/(resolution-1) px = (pu-(uc-field_of_view_mm/2))/pitch py = ((wc+field_of_view_mm/2)-pw)/pitch zbuf = np.full((resolution,resolution), -np.inf, np.float32) for tri in f: qx, qy, qz = px[tri], py[tri], depth[tri] xmin, xmax = max(0, int(np.ceil(qx.min()-.5))), min(resolution-1, int(np.floor(qx.max()-.5))) ymin, ymax = max(0, int(np.ceil(qy.min()-.5))), min(resolution-1, int(np.floor(qy.max()-.5))) if xmin > xmax or ymin > ymax: continue denom = (qy[1]-qy[2])*(qx[0]-qx[2])+(qx[2]-qx[1])*(qy[0]-qy[2]) if abs(denom) < 1e-12: continue yy, xx = np.mgrid[ymin:ymax+1,xmin:xmax+1] xx = xx+.5; yy = yy+.5 b0 = ((qy[1]-qy[2])*(xx-qx[2])+(qx[2]-qx[1])*(yy-qy[2]))/denom b1 = ((qy[2]-qy[0])*(xx-qx[2])+(qx[0]-qx[2])*(yy-qy[2]))/denom b2 = 1-b0-b1 current = zbuf[ymin:ymax+1,xmin:xmax+1] take = (b0>=-1e-7)&(b1>=-1e-7)&(b2>=-1e-7)&(b0*qz[0]+b1*qz[1]+b2*qz[2]>current) if np.any(take): current[take] = (b0*qz[0]+b1*qz[1]+b2*qz[2])[take] vis = np.isfinite(zbuf) return np.where(vis,zbuf,0).astype(np.float32), vis def save_rows(path: Path, rows: list[dict]) -> None: with path.open("w",newline="",encoding="utf8") as f: w=csv.DictWriter(f,fieldnames=list(dict.fromkeys(k for r in rows for k in r))) w.writeheader();w.writerows(rows) def process_part(part: str, root: Path) -> list[dict]: folder=HERE/"candidates"/part cad=prepared_cad(root/"data_processed/graphs/graph_5000.npz") rows=[] with np.load(folder/"image_256.npz") as z: ref_vis=z["visibility"].copy();ref_depth=z["depth_mm"].copy() for size in IMAGE_SIZES: start=time.perf_counter() frames=[raster_geometry(cad,v,size) for v in DIRECTIONS] depth=np.stack([a for a,b in frames]);vis=np.stack([b for a,b in frames]) path=folder/f"image_{size}.npz" np.savez_compressed(path,depth_mm=depth,visibility=vis, field_of_view_mm=np.full(6,FOV_MM,np.float32), pixel_pitch_mm=np.full(6,FOV_MM/(size-1),np.float32), channel_count=np.asarray(17), diagnostic_only=np.asarray(True), source_graph=np.asarray("graph_5000")) idx=np.rint(np.linspace(0,size-1,256)).astype(int) low=vis[:,idx][:,:,idx]; dlow=depth[:,idx][:,:,idx] common=low&ref_vis rows.append(dict(part=part,modality="image",candidate=size, fixed_tensor_shape=f"6x{size}x{size}",physical_unit_mm=FOV_MM/(size-1), fidelity_metric="silhouette_iou_to_256", fidelity=float((low&ref_vis).sum()/max(1,(low|ref_vis).sum())), secondary_metric="depth_mae_mm_on_common_pixels", secondary=float(np.mean(np.abs(dlow[common]-ref_depth[common]))), native_elements=int(6*size*size), raw_feature_mb=round(6*17*size*size*4/1e6,3), measured_occupancy=int(vis.sum()), generation_seconds=round(time.perf_counter()-start,3), candidate_storage="six-view depth and visibility; full 17-channel tensor cost estimated", candidate_file=str(path.relative_to(HERE)))) print(part,"image",size,"seconds",rows[-1]["generation_seconds"],flush=True) del depth,vis,frames with np.load(folder/"voxel_128.npz") as z: f=z["surface_sample_source_triangle"] b=z["surface_sample_barycentric"] fine=np.einsum("ni,nij->nj",b,cad["pos_mm"][cad["face"].T[f]]).astype(np.float32) for size,(pitch,shape) in VOXEL_SIZES.items(): start=time.perf_counter() vox=surface_voxels(cad,size,voxel_pitch_mm=pitch,grid_shape=shape) path=folder/f"voxel_{size}.npz" retained={k:v for k,v in vox.items() if not k.startswith("surface_sample_")} retained["correspondence_storage_omitted"] = np.asarray(True) np.savez_compressed(path,**retained) centers=vox["grid_origin_mm"]+(vox["voxel_coords"]+.5)*float(vox["voxel_pitch_mm"]) distance=cKDTree(centers).query(fine,workers=-1)[0] rows.append(dict(part=part,modality="voxel",candidate=size, fixed_tensor_shape="x".join(map(str,shape)),physical_unit_mm=pitch, fidelity_metric="surface_to_occupied_center_p95_mm", fidelity=float(np.quantile(distance,.95)), secondary_metric="occupied_surface_cells",secondary=int(len(centers)), native_elements=int(np.prod(shape)), raw_feature_mb=round(12*np.prod(shape)*4/1e6,3), generation_seconds=round(time.perf_counter()-start,3), candidate_storage="sparse 12-channel features and masks; correspondence omitted", candidate_file=str(path.relative_to(HERE)))) print(part,"voxel",size,"occupied",len(centers),"seconds",rows[-1]["generation_seconds"],flush=True) del vox,retained,centers,distance for size in POINT_SIZES: start=time.perf_counter();points=cad_points(cad,size) path=folder/f"points_{size}.npz";np.savez_compressed(path,**points) dist=cKDTree(points["pos_mm"]).query(cad["pos_mm"],workers=-1)[0] rows.append(dict(part=part,modality="point",candidate=size, fixed_tensor_shape=f"{size}x{points['x'].shape[1]}",physical_unit_mm="", fidelity_metric="cad_node_to_sample_p95_mm",fidelity=float(np.quantile(dist,.95)), secondary_metric="sampled_surface_points",secondary=size,native_elements=size, raw_feature_mb=round(points["x"].nbytes/1e6,3), generation_seconds=round(time.perf_counter()-start,3), candidate_file=str(path.relative_to(HERE)))) print(part,"point",size,flush=True) index=json.loads((folder/"high_graphs/index.json").read_text()) for rec in index["resolutions"]: size=rec["requested_node_budget"] if size not in GRAPH_SIZES: continue path=folder/"high_graphs"/f"graph_{size}.npz" rows.append(dict(part=part,modality="graph",candidate=size, fixed_tensor_shape=f"variable_actual_{rec['actual_nodes']}",physical_unit_mm="", fidelity_metric="source_to_graph_p95_mm", fidelity=rec["surface_distance_mm"]["source_to_simplified"]["p95_mm"], secondary_metric="six_view_silhouette_iou_to_12k",secondary="", native_elements=rec["actual_nodes"], raw_feature_mb=round(rec["actual_nodes"]*16*4/1e6,3), generation_seconds=round(rec["elapsed_s"],3), candidate_file=str(path.relative_to(HERE)))) return rows def main() -> None: rows=[] for part,root in PARTS.items():rows+=process_part(part,root) save_rows(HERE/"extended_measurements.csv",rows) print("Saved",len(rows),"new measured candidates",flush=True) if __name__=="__main__":main()