#!/usr/bin/env python3
"""Prepare topology-faithful, multi-resolution graphs from a triangle CAD mesh.

The requested resolutions are *vertex/node budgets*.  Open3D's quadric error
decimator takes triangle budgets, so this program deliberately records both the
requested node budget and the triangle budget passed to the decimator.  The
resulting number of graph nodes is always measured, never inferred.
"""
from __future__ import annotations

import argparse
import json
import os
import platform
import resource
import time
from pathlib import Path
from typing import Any

import numpy as np


DEFAULT_BUDGETS = (12000, 6000, 5000, 3000, 1500, 750, 300, 100)


def _root() -> Path:
    return Path(__file__).resolve().parents[1]


def _peak_memory_mb() -> float:
    # macOS reports ru_maxrss bytes; Linux reports KiB.
    rss = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss
    return float(rss / (1024 * 1024 if platform.system() == "Darwin" else 1024))


def _finite(a: np.ndarray) -> np.ndarray:
    return np.nan_to_num(a, nan=0.0, posinf=0.0, neginf=0.0).astype(np.float32)


def _safe_norm(v: np.ndarray) -> np.ndarray:
    return np.sqrt(np.maximum((v * v).sum(axis=-1), 0.0))


def _load_mesh(path: Path):
    import open3d as o3d

    mesh = o3d.io.read_triangle_mesh(str(path), enable_post_processing=False)
    if len(mesh.vertices) == 0 or len(mesh.triangles) == 0:
        raise ValueError(f"No triangles were read from {path}")
    # Binary/ASCII STL commonly carries one separate vertex record per triangle
    # corner.  Weld *exactly equal* coordinates to recover its encoded topology;
    # never use a tolerance/nearby-vertex merge, which could stitch CAD sheets.
    mesh.remove_duplicated_vertices()
    mesh.remove_duplicated_triangles()
    mesh.remove_degenerate_triangles()
    mesh.remove_unreferenced_vertices()
    mesh.compute_vertex_normals()
    return mesh


def _mesh_arrays(mesh) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    v = np.asarray(mesh.vertices, dtype=np.float64).copy()
    f = np.asarray(mesh.triangles, dtype=np.int64).copy()
    n = np.asarray(mesh.vertex_normals, dtype=np.float64).copy()
    nlen = _safe_norm(n)[:, None]
    n = np.divide(n, nlen, out=np.zeros_like(n), where=nlen > 1e-12)
    return v, f, n


def _topology(vertices: np.ndarray, faces: np.ndarray, normals: np.ndarray) -> dict[str, np.ndarray]:
    """Build edges solely from triangle topology (never spatial KNN)."""
    n = len(vertices)
    tri = faces
    directed = np.concatenate((tri[:, [0, 1]], tri[:, [1, 2]], tri[:, [2, 0]]), axis=0)
    undirected = np.sort(directed, axis=1)
    keys, inverse, counts = np.unique(undirected, axis=0, return_inverse=True, return_counts=True)
    boundary_u = counts == 1
    # A non-manifold edge is retained, but surfaced as a feature/uncertainty.
    face_normal = np.cross(vertices[tri[:, 1]] - vertices[tri[:, 0]], vertices[tri[:, 2]] - vertices[tri[:, 0]])
    flen = _safe_norm(face_normal)[:, None]
    face_normal = np.divide(face_normal, flen, out=np.zeros_like(face_normal), where=flen > 1e-12)
    # Map the three half edges back to triangle normals; this avoids an O(EF) walk.
    edge_face_normal = np.concatenate((face_normal, face_normal, face_normal), axis=0)
    sum_n = np.zeros((len(keys), 3), dtype=np.float64)
    np.add.at(sum_n, inverse, edge_face_normal)
    # For meaningful dihedral choose the first two incident face normals.
    first = np.full(len(keys), -1, dtype=np.int64)
    second = np.full(len(keys), -1, dtype=np.int64)
    for half_i, edge_i in enumerate(inverse):
        if first[edge_i] < 0:
            first[edge_i] = half_i
        elif second[edge_i] < 0:
            second[edge_i] = half_i
    dihedral = np.zeros(len(keys), dtype=np.float64)
    both = second >= 0
    dihedral[both] = np.arccos(np.clip((edge_face_normal[first[both]] * edge_face_normal[second[both]]).sum(1), -1, 1))
    # 30 degrees is a conventional, explicitly heuristic feature threshold.
    feature_u = boundary_u | (counts != 2) | (dihedral > np.deg2rad(30.0))
    src = np.concatenate((keys[:, 0], keys[:, 1]))
    dst = np.concatenate((keys[:, 1], keys[:, 0]))
    delta = vertices[dst] - vertices[src]
    length = _safe_norm(delta)
    direction = np.divide(delta, length[:, None], out=np.zeros_like(delta), where=length[:, None] > 1e-12)
    attrs = np.column_stack((length, direction, (normals[src] * normals[dst]).sum(1),
                             np.tile(dihedral, 2), np.tile(boundary_u, 2), np.tile(feature_u, 2)))
    degree = np.bincount(src, minlength=n).astype(np.int32)
    return {
        "edge_index": np.vstack((src, dst)).astype(np.int32), "edge_attr": _finite(attrs),
        "boundary_node": np.bincount(keys[boundary_u].ravel(), minlength=n).astype(bool),
        "feature_node": np.bincount(keys[feature_u].ravel(), minlength=n).astype(bool),
        "degree": degree, "unique_edges": keys.astype(np.int32), "edge_counts": counts,
    }


def _geometry_features(v: np.ndarray, f: np.ndarray, topo: dict[str, np.ndarray]) -> tuple[np.ndarray, np.ndarray, list[str]]:
    """Finite, intrinsic-ish local features; curvature is unsigned due to STL orientation."""
    n = len(v)
    e01, e12, e20 = v[f[:, 1]] - v[f[:, 0]], v[f[:, 2]] - v[f[:, 1]], v[f[:, 0]] - v[f[:, 2]]
    dbl_area = _safe_norm(np.cross(e01, -e20))
    area = dbl_area * 0.5
    vertex_area = np.zeros(n, dtype=np.float64)
    np.add.at(vertex_area, f.ravel(), np.repeat(area / 3.0, 3))
    # Triangle angles and Gaussian angle deficit.  Boundary vertices use pi.
    angles = []
    for a, b in ((1, 2), (2, 0), (0, 1)):
        x, y = v[f[:, a]] - v[f[:, 0]], v[f[:, b]] - v[f[:, 0]]
        # Recompute each corner below in the same compact pattern.
    for corner in range(3):
        a, b = f[:, (corner + 1) % 3], f[:, (corner + 2) % 3]
        x, y = v[a] - v[f[:, corner]], v[b] - v[f[:, corner]]
        denom = _safe_norm(x) * _safe_norm(y)
        angles.append(np.arccos(np.clip(np.divide((x * y).sum(1), denom, out=np.ones_like(denom), where=denom > 1e-15), -1, 1)))
    angle_sum = np.zeros(n, dtype=np.float64)
    for corner, ang in enumerate(angles): np.add.at(angle_sum, f[:, corner], ang)
    target = np.where(topo["boundary_node"], np.pi, 2 * np.pi)
    gaussian = np.divide(target - angle_sum, vertex_area, out=np.zeros(n), where=vertex_area > 1e-12)
    # Cotangent Laplace mean-curvature vector.  Degenerate triangles contribute 0.
    lap = np.zeros((n, 3), dtype=np.float64)
    for corner, ang in enumerate(angles):
        cot = np.divide(1.0, np.tan(ang), out=np.zeros_like(ang), where=np.abs(np.tan(ang)) > 1e-12)
        i, j = f[:, (corner + 1) % 3], f[:, (corner + 2) % 3]
        d = v[j] - v[i]
        np.add.at(lap, i, cot[:, None] * d)
        np.add.at(lap, j, -cot[:, None] * d)
    mean_mag = np.divide(_safe_norm(lap), 4.0 * vertex_area, out=np.zeros(n), where=vertex_area > 1e-12)
    curvature_valid = (vertex_area > 1e-10) & (topo["degree"] >= 2) & np.isfinite(mean_mag) & np.isfinite(gaussian)
    x = np.column_stack((vertex_area, angle_sum / np.maximum(topo["degree"], 1), mean_mag,
                         np.abs(gaussian), topo["boundary_node"], topo["feature_node"], topo["degree"]))
    valid = np.isfinite(x)
    valid[:, 2:4] &= curvature_valid[:, None]
    return _finite(x), valid.astype(bool), ["vertex_area_mm2", "mean_incident_angle_rad", "mean_curvature_abs_1_per_mm", "gaussian_curvature_abs_1_per_mm2", "boundary_node", "feature_node", "topological_degree"]


def _components(n: int, unique_edges: np.ndarray) -> int:
    parent = np.arange(n)
    def find(a):
        while parent[a] != a:
            parent[a] = parent[parent[a]]; a = parent[a]
        return a
    for a, b in unique_edges:
        a, b = find(int(a)), find(int(b))
        if a != b: parent[a] = b
    return len({find(i) for i in range(n)})


def _area_uniform_points(mesh, count: int, rng: np.random.Generator) -> np.ndarray:
    """Deterministic area-uniform triangle samples, rather than vertex-biased samples."""
    v, f = np.asarray(mesh.vertices), np.asarray(mesh.triangles)
    a = _safe_norm(np.cross(v[f[:, 1]] - v[f[:, 0]], v[f[:, 2]] - v[f[:, 0]])) * .5
    ids = rng.choice(len(f), size=count, replace=True, p=a / a.sum())
    r1, r2 = rng.random(count), rng.random(count)
    flip = r1 + r2 > 1; r1[flip] = 1-r1[flip]; r2[flip] = 1-r2[flip]
    return v[f[ids, 0]] + r1[:, None] * (v[f[ids, 1]] - v[f[ids, 0]]) + r2[:, None] * (v[f[ids, 2]] - v[f[ids, 0]])


def _sample_surface_distance(source_mesh, simplified_mesh, sample_count: int = 20000) -> tuple[dict[str, dict[str, float]], str]:
    """Bidirectional, area-uniform surface sampling queried against mesh triangles."""
    import open3d as o3d
    rng = np.random.default_rng(20260914)
    src_points = _area_uniform_points(source_mesh, sample_count, rng)
    simp_points = _area_uniform_points(simplified_mesh, sample_count, rng)
    try:
        target_s = o3d.t.geometry.RaycastingScene(); target_s.add_triangles(o3d.t.geometry.TriangleMesh.from_legacy(simplified_mesh))
        target_o = o3d.t.geometry.RaycastingScene(); target_o.add_triangles(o3d.t.geometry.TriangleMesh.from_legacy(source_mesh))
        d_forward = _safe_norm(src_points - target_s.compute_closest_points(o3d.core.Tensor(src_points.astype(np.float32)))["points"].numpy())
        d_reverse = _safe_norm(simp_points - target_o.compute_closest_points(o3d.core.Tensor(simp_points.astype(np.float32)))["points"].numpy())
        method = "Open3D RaycastingScene closest point on triangles; deterministic area-uniform samples"
    except Exception as exc:  # still truthful when an old Open3D lacks raycasting
        from scipy.spatial import cKDTree
        d_forward, _ = cKDTree(np.asarray(simplified_mesh.vertices)).query(src_points, workers=-1)
        d_reverse, _ = cKDTree(np.asarray(source_mesh.vertices)).query(simp_points, workers=-1)
        method = f"fallback nearest vertices (NOT surface distance): {type(exc).__name__}"
    def stats(d): return {"mean_mm": float(d.mean()), "p95_mm": float(np.percentile(d, 95)), "max_mm": float(d.max()), "samples": int(len(d))}
    return {"source_to_simplified": stats(d_forward), "simplified_to_source": stats(d_reverse)}, method


def _source_correspondence(source_vertices: np.ndarray, vertices: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    from scipy.spatial import cKDTree
    d, ix = cKDTree(source_vertices).query(vertices, workers=-1)
    return ix.astype(np.int32), d.astype(np.float32)


def _external_curvature(path: Path | None, vertices: np.ndarray, normals: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray, str | None]:
    """Map supplied analytic samples only when explicitly provided for this mesh.

    Mapping is nearest sample correspondence, with a conservative distance mask;
    it never silently substitutes samples from the reference CAD for an alternate
    ``--mesh``.  Signed principal values are retained as supplied, while the
    discrete unsigned channels remain available independently.
    """
    values = np.zeros((len(vertices), 3), np.float32); valid = np.zeros((len(vertices), 3), bool)
    distance = np.full(len(vertices), np.nan, np.float32)
    if path is None: return values, valid, distance, None
    from scipy.spatial import cKDTree
    a = np.load(path, allow_pickle=False)
    required = {"points_mm", "k1_per_mm", "k2_per_mm"}
    if not required.issubset(a.files):
        return values, valid, distance, f"ignored: missing {sorted(required-set(a.files))}"
    d, ix = cKDTree(a["points_mm"]).query(vertices, workers=-1)
    k1, k2 = a["k1_per_mm"][ix], a["k2_per_mm"][ix]
    sampled_normal = a["parameter_normals"][ix] if "parameter_normals" in a.files else np.full_like(vertices, np.nan)
    alignment = np.abs((normals * sampled_normal).sum(1))
    # Thin opposite sheets can be close in Euclidean space.  Both a 3 mm position
    # bound and orientation-free normal agreement are required for this mapping.
    close = (d < 3.0) & (alignment > .8)
    # Principal signs depend on face orientation; expose magnitudes.  Their product
    # is Gaussian curvature and retains its sign (saddle vs elliptic surface).
    raw = np.column_stack((np.abs(k1), np.abs(k2), k1*k2))
    values[:] = np.nan_to_num(raw, nan=0., posinf=0., neginf=0.)
    valid[:] = close[:, None] & np.isfinite(raw)
    distance[:] = d
    return values, valid, distance, "nearest supplied analytic samples; valid iff distance < 3 mm and abs(normal dot parameter_normal) > 0.8; principal values stored as magnitudes, Gaussian retains sign"


def _save_figure(graphs: list[dict[str, Any]], out: Path) -> None:
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    rows = (len(graphs) + 4) // 4
    fig = plt.figure(figsize=(13, 3.5 * rows))
    ax = fig.add_subplot(rows, 4, 1)
    ax.plot([g["actual_nodes"] for g in graphs], [g["distance"]["source_to_simplified"]["p95_mm"] for g in graphs], "o-", label="source → simplified")
    ax.plot([g["actual_nodes"] for g in graphs], [g["distance"]["simplified_to_source"]["p95_mm"] for g in graphs], "s--", label="simplified → source")
    ax.set_xscale("log"); ax.set_xlabel("actual graph nodes"); ax.set_ylabel("surface p95 (mm)"); ax.grid(True, alpha=.3); ax.legend(fontsize=7)
    for i, g in enumerate(graphs, 2):
        ax = fig.add_subplot(rows, 4, i, projection="3d")
        p, c, f = g["pos_norm"], g["curvature"], g["faces"]
        # A bounded, deterministic triangle subset keeps the all-budget figure
        # legible while preserving actual mesh edges rather than point-only views.
        use = f if len(f) <= 8000 else f[np.linspace(0, len(f)-1, 8000, dtype=int)]
        ax.plot_trisurf(p[:, 0], p[:, 1], p[:, 2], triangles=use, cmap="magma", shade=False, linewidth=.08, edgecolor="#3b3b3b", alpha=.96, antialiased=False)
        ax.set_title(f"request {g['requested_nodes']} nodes\nactual {g['actual_nodes']}", fontsize=8)
        ax.view_init(elev=22, azim=-55); ax.set_box_aspect((1, 1, 1)); ax.set_xlim(-.52,.52); ax.set_ylim(-.52,.52); ax.set_zlim(-.52,.52)
        ax.set_axis_off()
    fig.tight_layout(); fig.savefig(out, dpi=160); plt.close(fig)


def _write_viewer(graphs: list[dict[str, Any]], out_dir: Path) -> None:
    # The small data file is local and deliberately separate from NPZ training data.
    payload = {"graphs": [{"budget": g["requested_nodes"], "actual": g["actual_nodes"],
               # ``faces`` is already [F, 3].  Do not transpose: flattening a
               # [3, F] transpose in Three.js would connect unrelated vertices.
               "pos": g["pos_norm"].round(6).tolist(), "face": g["faces"].tolist(),
               "edge": g["edge_index"].T.tolist(), "curvature": g["curvature"].round(6).tolist()} for g in graphs]}
    (out_dir / "viewer_data.js").write_text("window.CAD_GRAPH_DATA=" + json.dumps(payload, separators=(",", ":")) + ";", encoding="utf8")
    html = '''<!doctype html><meta charset="utf-8"><title>CAD graph resolutions</title>
<style>body{margin:0;background:#e8edf2;color:#15212b;font:14px sans-serif}#bar{position:fixed;z-index:2;padding:10px;background:#ffffffdd;border-radius:0 0 6px 0}canvas{display:block}label{margin-right:9px}</style>
<div id="bar"><label>Budget <select id="budget"></select></label><label><input id="wire" type="checkbox" checked> wireframe</label><label><input id="nodes" type="checkbox"> nodes</label><label><input id="curv" type="checkbox" checked> curvature</label><span id="info"></span></div>
<script src="../../datas_2/viewer/vendor/three.min.js"></script><script src="../../datas_2/viewer/vendor/OrbitControls.js"></script><script src="viewer_data.js"></script><script>
const D=window.CAD_GRAPH_DATA,sel=document.querySelector('#budget'), scene=new THREE.Scene(),cam=new THREE.PerspectiveCamera(45,innerWidth/innerHeight,.01,100),r=new THREE.WebGLRenderer({antialias:true});r.setClearColor(0xe8edf2);r.setSize(innerWidth,innerHeight);document.body.append(r.domElement);cam.up.set(0,0,1);cam.position.set(0,-1.4,.2);cam.lookAt(0,0,0);const controls=new THREE.OrbitControls(cam,r.domElement);controls.target.set(0,0,0);controls.update();scene.add(new THREE.HemisphereLight(0xffffff,.6,2));let group;D.graphs.forEach((g,i)=>sel.add(new Option(g.budget+' requested / '+g.actual+' actual',i)));function dispose(o){o.traverse(x=>{if(x.geometry)x.geometry.dispose();if(x.material){let a=Array.isArray(x.material)?x.material:[x.material];a.forEach(m=>m.dispose())}})}function draw(){if(group){scene.remove(group);dispose(group)}let g=D.graphs[sel.value], geo=new THREE.BufferGeometry();geo.setAttribute('position',new THREE.Float32BufferAttribute(g.pos.flat(),3));geo.setIndex(g.face.flat());geo.computeVertexNormals();let col=[],mx=Math.max(...g.curvature,1e-9);for(let x of g.curvature){let t=Math.min(1,x/mx);col.push(t,.1,.9-t)}geo.setAttribute('color',new THREE.Float32BufferAttribute(col,3));group=new THREE.Group();let mat=new THREE.MeshStandardMaterial({vertexColors:curv.checked,wireframe:wire.checked,side:THREE.DoubleSide,roughness:.72,color:curv.checked?0xffffff:0x7da9bd});group.add(new THREE.Mesh(geo,mat));if(nodes.checked){let pg=new THREE.BufferGeometry();pg.setAttribute('position',geo.getAttribute('position'));group.add(new THREE.Points(pg,new THREE.PointsMaterial({size:.009,color:0x15212b})));}scene.add(group);info.textContent=' '+g.actual+' actual nodes; requested budget is a node budget. Legend: blue→magenta = discrete |mean curvature|; analytic values are masked and not rendered.'}for(let e of [sel,wire,nodes,curv])e.onchange=draw;addEventListener('resize',()=>{cam.aspect=innerWidth/innerHeight;cam.updateProjectionMatrix();r.setSize(innerWidth,innerHeight)});(function loop(){requestAnimationFrame(loop);controls.update();r.render(scene,cam)})();draw();</script>'''
    vendor_path = os.path.relpath(_root() / "datas_2/viewer/vendor", out_dir)
    html = html.replace("../../datas_2/viewer/vendor", vendor_path)
    (out_dir / "viewer.html").write_text(html, encoding="utf8")


def prepare(mesh_path: Path, curvature_path: Path | None, out_dir: Path, budgets: tuple[int, ...] = DEFAULT_BUDGETS, results_dir: Path | None = None) -> dict[str, Any]:
    import open3d as o3d
    started = time.perf_counter(); out_dir.mkdir(parents=True, exist_ok=True)
    source = _load_mesh(mesh_path); source_v, source_f, _ = _mesh_arrays(source)
    center = (source_v.min(0) + source_v.max(0)) / 2; scale = float(np.max(source_v.max(0) - source_v.min(0)))
    if scale <= 0: raise ValueError("mesh has zero bounding-box scale")
    source_area = float(source.get_surface_area())
    graphs: list[dict[str, Any]] = []; index: dict[str, Any] = {"source": str(mesh_path), "normalization": {"center_mm": center.tolist(), "scale_mm": scale}, "resolutions": []}
    for requested in budgets:
        t0 = time.perf_counter()
        # QEM's input is a triangle count. ~2V is an Euler-characteristic heuristic,
        # not a claim about output nodes; actual output is measured below.
        target_tri = min(len(source_f), max(4, int(requested * 2)))
        simple = source.simplify_quadric_decimation(target_tri)
        simple.remove_degenerate_triangles(); simple.remove_unreferenced_vertices(); simple.compute_vertex_normals()
        v, f, normal = _mesh_arrays(simple); topo = _topology(v, f, normal)
        base_x, base_valid, base_names = _geometry_features(v, f, topo)
        # Curvature and normal channels have separate explicit masks.  A singular
        #/zero vertex normal must not make a curvature estimate appear reliable.
        normal_good = np.isfinite(normal).all(1) & (_safe_norm(normal) > .5)
        base_valid[:, 2:4] &= normal_good[:, None]
        external_x, external_valid, external_distance, external_note = _external_curvature(curvature_path, v, normal)
        pos_norm = (v-center)/scale
        x = np.column_stack((pos_norm, normal, base_x, external_x)).astype(np.float32)
        feature_valid = np.column_stack((np.isfinite(pos_norm), np.repeat(normal_good[:, None], 3, axis=1), base_valid, external_valid))
        names = ["pos_norm_x", "pos_norm_y", "pos_norm_z", "normal_x", "normal_y", "normal_z"] + base_names + ["analytic_k1_abs_1_per_mm", "analytic_k2_abs_1_per_mm", "analytic_gaussian_signed_1_per_mm2"]
        original_index, correspondence_mm = _source_correspondence(source_v, v)
        distance, distance_method = _sample_surface_distance(source, simple)
        area = float(simple.get_surface_area()); actual = len(v)
        exceeded = actual > requested; floor = exceeded
        floor_nodes = None
        if floor:
            # Only call a budget a floor after asking QEM for its minimum legal
            # triangle target as a separate check; this mesh retains 329 nodes.
            minimum = source.simplify_quadric_decimation(4); minimum.remove_unreferenced_vertices()
            floor_nodes = int(len(minimum.vertices))
            floor = floor_nodes > requested
        record = {"requested_node_budget": requested, "decimator_target_triangles": target_tri, "actual_nodes": actual,
                  "actual_triangles": int(len(f)), "node_budget_exceeded": exceeded,
                  "node_budget_note": "measured output count may differ because QEM accepts a triangle target, not a node target",
                  "confirmed_qem_node_floor": floor_nodes if floor else None,
                  "floor_note": "QEM minimum triangle request still retained this many nodes" if floor else None,
                  "surface_distance_mm": distance, "surface_distance_method": distance_method,
                  "components": _components(actual, topo["unique_edges"]), "boundary_edges": int((topo["edge_counts"] == 1).sum()),
                  "nonmanifold_edges": int((topo["edge_counts"] > 2).sum()), "curvature_valid_fraction": float(base_valid[:, 2].mean()),
                  "analytic_curvature_valid_fraction": float(external_valid[:, 0].mean()), "analytic_curvature_mapping": external_note,
                  "surface_area_mm2": area, "area_change_fraction": (area - source_area) / source_area,
                  "elapsed_s": time.perf_counter() - t0, "peak_process_memory_mb": _peak_memory_mb()}
        np.savez_compressed(out_dir / f"graph_{requested}.npz", pos_mm=v.astype(np.float32), pos_norm=pos_norm.astype(np.float32), normal=normal.astype(np.float32),
                            x=x, feature_valid=feature_valid, feature_names=np.asarray(names, dtype="<U48"), edge_index=topo["edge_index"], edge_attr=topo["edge_attr"],
                            edge_feature_names=np.asarray(["length_mm","direction_x","direction_y","direction_z","normal_dot","dihedral_rad","boundary_edge","feature_edge"], dtype="<U32"),
                            face=f.T.astype(np.int32), original_vertex_index=original_index, original_vertex_distance_mm=correspondence_mm,
                            analytic_curvature_mapping_distance_mm=external_distance,
                            training_valid=np.zeros(actual, dtype=bool))
        index["resolutions"].append(record)
        graphs.append({"requested_nodes": requested, "actual_nodes": actual, "pos_norm": pos_norm, "faces": f, "edge_index": topo["edge_index"], "curvature": base_x[:, 2], "distance": distance})
    index["source_vertices"] = int(len(source_v)); index["source_triangles"] = int(len(source_f)); index["source_area_mm2"] = source_area
    index["curvature_provenance"] = "x includes distinct discrete mesh mean-curvature magnitude/absolute Gaussian and (when explicitly supplied) analytic principal/gaussian samples. STL normal orientation is not assumed globally outward."
    index["external_curvature_input"] = str(curvature_path) if curvature_path and curvature_path.exists() else None
    index["total_elapsed_s"] = time.perf_counter() - started
    index["peak_process_memory_mb"] = _peak_memory_mb()
    (out_dir / "index.json").write_text(json.dumps(index, indent=2), encoding="utf8")
    results_dir = results_dir or _root()/"results/graphs"; results_dir.mkdir(parents=True, exist_ok=True)
    _save_figure(graphs, results_dir / "comparison.png"); _write_viewer(graphs, results_dir)
    return index


def main() -> None:
    p = argparse.ArgumentParser(description=__doc__)
    p.add_argument("--mesh", type=Path, default=_root()/"analysis/results/design_fine.stl", help="STL/mesh input (default reference CAD mesh)")
    p.add_argument("--curvature", type=Path, default=None, help="Explicit analytic curvature NPZ corresponding to --mesh (not implicitly reused for changed CAD)")
    p.add_argument("--out", type=Path, default=_root()/"data_processed/graphs"); p.add_argument("--results", type=Path, default=_root()/"results/graphs")
    p.add_argument("--budgets", type=int, nargs="+", default=list(DEFAULT_BUDGETS), help="Requested NODE budgets; output counts may differ")
    a = p.parse_args()
    # Reference analytic samples are only selected for the reference mesh.  An
    # alternate --mesh must opt into its own --curvature input.
    curvature = a.curvature or (_root()/"analysis/results/cad_surface_curvature.npz" if a.mesh.resolve() == (_root()/"analysis/results/design_fine.stl").resolve() else None)
    summary = prepare(a.mesh, curvature, a.out, tuple(a.budgets), a.results); print(json.dumps(summary, indent=2))

if __name__ == "__main__": main()