#!/usr/bin/env python3
"""Create deterministic, local HTML/PNG galleries for processed POC artifacts.

This is deliberately a reader: it never changes NPZ inputs, their schemas, or
their family index.json files.  Run from the repository root with:

    .venv/bin/python 전처리/visualize_processed.py
"""
from __future__ import annotations

import argparse
import html
import json
from pathlib import Path
from typing import Any

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
import numpy as np

# Process channel names are Korean in the supplied source. Apple SD Gothic Neo is
# available on the target macOS environment; DejaVu remains a portable fallback.
plt.rcParams["font.family"] = ["Apple SD Gothic Neo", "DejaVu Sans"]
plt.rcParams["axes.unicode_minus"] = False

FAMILIES = ("images", "voxels", "points", "graphs", "brep", "tabular", "paired")
CSS = """
body{font:15px/1.45 system-ui,sans-serif;margin:0;background:#f7f8fa;color:#18212f}
main{max-width:1280px;margin:auto;padding:28px}.card{background:#fff;border:1px solid #dfe5ec;border-radius:8px;padding:18px;margin:16px 0}
img{max-width:100%;height:auto;border:1px solid #dfe5ec}.thumb{max-height:250px;object-fit:contain}a{color:#0758aa}table{border-collapse:collapse;width:100%;font-size:13px}th,td{border:1px solid #dfe5ec;padding:5px;text-align:left}th{background:#edf3fa}code{font-size:12px} .muted{color:#526170}
"""


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


def rel(target: Path, start: Path) -> str:
    return str(Path(__import__("os").path.relpath(target, start))).replace("\\", "/")


def scalar(value: np.ndarray | Any) -> Any:
    a = np.asarray(value)
    return a.item() if a.ndim == 0 else a


def strings(value: np.ndarray) -> list[str]:
    return [str(x) for x in np.asarray(value).tolist()]


def safe_limits(values: np.ndarray, symmetric: bool = False) -> dict[str, float]:
    a = np.asarray(values, dtype=float); a = a[np.isfinite(a)]
    if not len(a): return {"vmin": 0., "vmax": 1.}
    if symmetric:
        q = max(float(np.percentile(np.abs(a), 99)), 1e-12); return {"vmin": -q, "vmax": q}
    lo, hi = np.percentile(a, [1, 99])
    if hi <= lo: hi = lo + max(abs(lo)*.01, 1e-12)
    return {"vmin": float(lo), "vmax": float(hi)}


def save(fig: plt.Figure, path: Path) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    # Do not call tight_layout: it conflicts with shared colorbars and 3D axes.
    # Fixed generous figure sizes plus the tight bounding box retain labels safely.
    fig.savefig(path, dpi=150, bbox_inches="tight"); plt.close(fig)


def plot_images(d: dict[str, np.ndarray], out: Path, title: str) -> str:
    depth, curv, ok = d["depth_mm"], d["curvature_magnitude_per_mm"], d["curvature_valid"]
    visible = d.get("visibility", depth != 0)
    names = ["+X", "−X", "+Y", "−Y", "+Z", "−Z"]
    fig = plt.figure(figsize=(21, 6.7), layout="constrained")
    grid = fig.add_gridspec(2, 7, width_ratios=[1,1,1,1,1,1,.055])
    axes = np.asarray([[fig.add_subplot(grid[row, col]) for col in range(6)] for row in range(2)])
    depth_lim, lim = safe_limits(depth[visible]), safe_limits(curv[ok]); extent = None
    if "field_of_view_mm" in d:
        side = float(d["field_of_view_mm"][0]); extent = [-side/2, side/2, -side/2, side/2]
    for i, name in enumerate(names):
        dm = np.ma.masked_where(~visible[i], depth[i]); cm = np.ma.masked_where(~ok[i], curv[i])
        im0 = axes[0,i].imshow(dm, cmap="viridis", extent=extent, origin="upper", **depth_lim)
        im1 = axes[1,i].imshow(cm, cmap="magma", extent=extent, origin="upper", **lim)
        for a in (axes[0,i], axes[1,i]): a.set_aspect("equal"); a.set_title(name); a.set_xticks([]); a.set_yticks([])
        if i == 0: axes[0,i].set_ylabel("depth\nmm"); axes[1,i].set_ylabel("|curvature|\n1/mm")
    fig.colorbar(im0, cax=fig.add_subplot(grid[0,6]), label="front-most depth (mm), shared scale")
    fig.colorbar(im1, cax=fig.add_subplot(grid[1,6]), label="curvature magnitude (1/mm); white = unknown")
    fig.suptitle(f"{title}: six orthographic CAD views; square physical FOV, no axis stretching")
    save(fig, out)
    return "Six front-most z-buffer CAD views. Curvature white means source curvature mapping is invalid/unknown; no thickness, thermal, or stress field is shown."


def _set_3d_aspect(ax: Any, p: np.ndarray) -> None:
    span = np.maximum(np.ptp(p, axis=0), 1e-6); ax.set_box_aspect(span)


def plot_points(d: dict[str, np.ndarray], out: Path, title: str) -> str:
    pkey = "pos_mm" if "pos_mm" in d else "pos_native"; p = d[pkey]
    names = strings(d.get("feature_names", np.asarray([])))
    x = d.get("x")
    if "curvature_magnitude_per_mm" in names:
        j = names.index("curvature_magnitude_per_mm"); value = x[:,j]; label = "curvature magnitude (1/mm)"
        valid = d.get("feature_valid", np.ones_like(x, bool))[:,j]
    elif "local_surface_variation" in names:
        j = names.index("local_surface_variation"); value = x[:,j]; label = "local surface variation (unitless)"
        valid = d.get("feature_valid", np.ones_like(x, bool))[:,j]
    else: value, valid, label = np.zeros(len(p)), np.zeros(len(p),bool), "no scalar descriptor"
    fig = plt.figure(figsize=(14, 4.5)); a0 = fig.add_subplot(1,3,1, projection="3d"); a1=fig.add_subplot(1,3,2); a2=fig.add_subplot(1,3,3)
    color = np.where(valid, value, np.nan); kw = safe_limits(value[valid])
    a0.scatter(p[:,0],p[:,1],p[:,2], c=color, s=.8, cmap="viridis", **kw); _set_3d_aspect(a0,p); a0.set_title("3D sample; physical aspect")
    im = a1.scatter(p[:,0],p[:,2], c=color, s=.8, cmap="viridis", **kw); a1.set_aspect("equal"); a1.set_xlabel(f"X ({'mm' if pkey=='pos_mm' else 'native'})"); a1.set_ylabel(f"Z ({'mm' if pkey=='pos_mm' else 'native'})"); a1.set_title("XZ broadside")
    a2.scatter(p[:,0],p[:,1], c=color, s=.8, cmap="viridis", **kw); a2.set_aspect("equal"); a2.set_xlabel(f"X ({'mm' if pkey=='pos_mm' else 'native'})"); a2.set_ylabel(f"Y ({'mm' if pkey=='pos_mm' else 'native'})"); a2.set_title("XY plan")
    fig.colorbar(im, ax=[a0,a1,a2], shrink=.72, label=label)
    frame = str(scalar(d["frame"])) if "frame" in d else "CAD millimetre frame"
    fig.suptitle(f"{title}: {len(p):,} points — {frame}", y=1.02); save(fig,out)
    diagnostic = bool(scalar(d.get("diagnostic_post_process", False)))
    return ("Diagnostic post-process scan in its native unregistered frame; it is not a CAD-aligned predictor." if diagnostic else
            "Deterministic CAD surface samples with triangle/barycentric correspondence. Scalar is geometric curvature, not stress or a defect label.")


def plot_sparse_voxels(d: dict[str, np.ndarray], out: Path, title: str) -> str:
    c, x = d["voxel_coords"], d["x"]; names = strings(d["feature_names"])
    j = names.index("curvature_magnitude_per_mm") if "curvature_magnitude_per_mm" in names else 0
    valid = d.get("feature_valid", np.ones_like(x,bool))[:,j]; val=np.where(valid,x[:,j],np.nan)
    fig=plt.figure(figsize=(11,5)); a0=fig.add_subplot(1,2,1,projection="3d"); a1=fig.add_subplot(1,2,2)
    kw=safe_limits(x[valid,j]); im=a0.scatter(c[:,0],c[:,1],c[:,2],s=4,c=val,cmap="magma",**kw); _set_3d_aspect(a0,np.vstack([np.zeros(3),np.asarray(d["grid_dims"])])); a0.set_title("occupied surface voxels")
    a1.scatter(c[:,0],c[:,2],s=4,c=val,cmap="magma",**kw); a1.set_aspect("equal"); a1.set_xlabel("voxel X"); a1.set_ylabel("voxel Z"); a1.set_title("XZ occupied surface projection")
    fig.colorbar(im,ax=[a0,a1],label="curvature magnitude (1/mm); white = unknown")
    pitch=float(scalar(d.get("voxel_pitch_mm", np.nan))); dims=np.asarray(d["grid_dims"]).tolist()
    fig.suptitle(f"{title}: {len(c):,} occupied surface voxels; grid {dims}, pitch {pitch:.4g} mm",y=1.02); save(fig,out)
    return "Surface-only triangle-sampled occupancy, not filled-solid occupancy, SDF, thickness, or an interior field."


def plot_dense_voxels(d: dict[str, np.ndarray], out: Path, title: str) -> str:
    x=d["x"]; names=strings(d["feature_names"]); occ=names.index("surface_occupancy") if "surface_occupancy" in names else 0
    cur=names.index("curvature_magnitude_per_mm") if "curvature_magnitude_per_mm" in names else occ
    occupied=x[occ]>0; coords=np.argwhere(occupied); values=x[cur][occupied]
    fig=plt.figure(figsize=(12,4.5)); a0=fig.add_subplot(1,3,1,projection="3d"); a1=fig.add_subplot(1,3,2); a2=fig.add_subplot(1,3,3)
    im=a0.scatter(coords[:,0],coords[:,1],coords[:,2],s=5,c=values,cmap="magma",**safe_limits(values)); _set_3d_aspect(a0,np.vstack([np.zeros(3),np.asarray(d["grid_dims"])])); a0.set_title("occupied surface voxels")
    y=x.shape[2]//2; a1.imshow(np.ma.masked_where(~occupied[:,y,:],x[cur,:,y,:]),cmap="magma",**safe_limits(values),origin="lower",aspect="equal"); a1.set_title(f"curvature XZ slice y={y}"); a1.set_xlabel("Z voxel"); a1.set_ylabel("X voxel")
    a2.imshow(occupied.max(axis=1),cmap="Greys",origin="lower",aspect="equal"); a2.set_title("surface occupancy XZ projection"); a2.set_xlabel("Z voxel"); a2.set_ylabel("X voxel")
    fig.colorbar(im,ax=[a0,a1],label="curvature magnitude (1/mm)"); fig.suptitle(f"{title}: dense tensor {list(x.shape)}; surface-only",y=1.02); save(fig,out)
    return "Bounded dense tensor with surface occupancy only; it is not a filled CAD solid or a signed-distance field."


def graph_requested(family_dir: Path, stem: str) -> int | None:
    idx=family_dir/"index.json"
    if not idx.exists(): return None
    try:
        for row in json.loads(idx.read_text(encoding="utf-8")).get("resolutions",[]):
            if stem == f"graph_{row.get('requested_node_budget')}": return int(row["requested_node_budget"])
    except (json.JSONDecodeError, KeyError, TypeError): pass
    try: return int(stem.rsplit("_",1)[1])
    except ValueError: return None


def plot_graph(d: dict[str, np.ndarray], out: Path, title: str, requested: int | None) -> str:
    p=d["pos_mm"]; e=d["edge_index"]; x=d["x"]; names=strings(d.get("feature_names",np.asarray([]))); valid=d.get("feature_valid",np.ones_like(x,bool))
    # Use analytic curvature only when that mapped source field is valid.
    candidates=[i for i,n in enumerate(names) if "analytic_k1" in n]
    if not candidates: candidates=[i for i,n in enumerate(names) if "mean_curvature" in n]
    j=candidates[0] if candidates else 0
    vv=valid[:,j]; val=np.where(vv,x[:,j],np.nan)
    fig,ax=plt.subplots(1,2,figsize=(13,5))
    edges=np.stack([p[e[0]][:,[0,2]], p[e[1]][:,[0,2]]],axis=1)
    ax[0].add_collection(LineCollection(edges,colors="#386a8c",linewidths=.15,alpha=.48)); ax[0].scatter(p[:,0],p[:,2],s=1,c="#13283a"); ax[0].autoscale(); ax[0].set_aspect("equal"); ax[0].set_xlabel("X (mm)"); ax[0].set_ylabel("Z (mm)"); ax[0].set_title("topology wireframe, XZ broadside")
    im=ax[1].scatter(p[:,0],p[:,2],s=max(.6, 3000/len(p)),c=val,cmap="magma",**safe_limits(x[vv,j])); ax[1].set_aspect("equal"); ax[1].set_xlabel("X (mm)"); ax[1].set_ylabel("Z (mm)"); ax[1].set_title("source-mapped curvature; gray/white = invalid")
    fig.colorbar(im,ax=ax[1],label=f"{names[j] if names else 'curvature'}")
    fig.suptitle(f"{title}: requested {requested if requested is not None else 'unknown'} nodes; actual {len(p):,}; {e.shape[1]:,} directed edges",y=1.02); save(fig,out)
    return f"Topology uses the exported edge_index. Curvature panel displays {names[j] if names else 'a scalar'} only where its source-valid mask is true; no physical label is implied."


def plot_brep(d: dict[str, np.ndarray], out: Path, title: str) -> str:
    pos, trim=d["pos_mm"],d["trim_valid"]; normal_ok=d["normal_valid"]; curv=d["curvature"]; curv_ok=d["curvature_valid"]
    cent=np.zeros((len(pos),3));
    centroid_valid=np.any(trim,axis=(1,2))
    for i in range(len(pos)):
        q=pos[i][trim[i]]; cent[i]=q.mean(0) if len(q) else 0
    fig=plt.figure(figsize=(15,8)); a0=fig.add_subplot(2,3,1); edge=d.get("edge_index")
    omitted = int((~centroid_valid).sum())
    if edge is not None and edge.shape[1]:
        keep=centroid_valid[edge[0]] & centroid_valid[edge[1]]
        seg=np.stack([cent[edge[0,keep]][:,[0,2]],cent[edge[1,keep]][:,[0,2]]],axis=1); a0.add_collection(LineCollection(seg,colors="#497b96",linewidths=.12,alpha=.42))
        a0.scatter(cent[centroid_valid,0],cent[centroid_valid,2],s=.5,c="#172b3a"); note=f"real directed face adjacency: {edge.shape[1]:,} edges ({int(keep.sum()):,} drawable); {omitted} faces without UV samples omitted from preview"
    else:
        a0.scatter(cent[:,0],cent[:,2],s=.5,c="#172b3a"); note="no positive adjacency exported; points only"
    a0.autoscale(); a0.set_aspect("equal"); a0.set_title("face-centroid graph, XZ (mm)"); a0.set_xlabel("X (mm)"); a0.set_ylabel("Z (mm)")
    kinds=strings(d.get("surface_kind",np.asarray([""]*len(pos)))); area=d.get("face_area_mm2",np.ones(len(pos)))
    chosen=[]
    for kind in dict.fromkeys(kinds):
        ids=np.where(np.asarray(kinds)==kind)[0]
        if len(ids): chosen.append(ids[np.argmax(area[ids])])
        if len(chosen)==5: break
    uv_cmap=plt.get_cmap("magma").copy(); uv_cmap.set_bad("#d7dde5")
    uv_axes=[fig.add_subplot(2,3,i) for i in range(2,7)]
    uv_values=curv[:,:,:,0][trim&curv_ok]
    uv_lim=safe_limits(uv_values); feature_name=strings(d.get("curvature_feature_names",np.asarray(["curvature_channel_0"])))[0]
    last_im=None
    for a,face_id in zip(uv_axes,chosen):
        data=np.ma.masked_where(~(trim[face_id]&curv_ok[face_id]),curv[face_id,:,:,0]); last_im=a.imshow(data,cmap=uv_cmap,**uv_lim); a.set_title(f"face {face_id}: {kinds[face_id]}\ntrim+curvature-valid UV grid (gray = masked)"); a.set_xticks([]);a.set_yticks([])
    if last_im is not None: fig.colorbar(last_im,ax=uv_axes,shrink=.62,pad=.02,label=f"{feature_name} (1/mm), shared scale")
    fig.suptitle(f"{title}: {len(pos):,} STEP B-Rep faces; {note}",y=.99); save(fig,out)
    return f"{note}; original edge_index is unchanged. UV panels show {feature_name} (1/mm) with a shared numeric scale and actual trim/curvature-valid masks; they are geometric fields, not thermal/stress/defect maps."


def plot_tabular(d: dict[str, np.ndarray], out: Path, title: str) -> str:
    present=d["process_feature_valid"].astype(float); names=strings(d["process_feature_names"]); runs=strings(d["run_id"])
    fig,ax=plt.subplots(figsize=(12,max(7,len(names)*.18))); im=ax.imshow(present.T,cmap="Blues",vmin=0,vmax=1,aspect="auto")
    ax.set_xticks(range(len(runs)),runs); ax.set_yticks(range(len(names)),names,fontsize=6); ax.set_xlabel("run"); ax.set_title("process-channel availability (not magnitude; units differ)")
    fig.colorbar(im,ax=ax,ticks=[0,1],label="missing / observed"); fig.suptitle(f"{title}: {len(names)} raw process channels",y=1.01); save(fig,out)
    return "PNG is an availability grid, intentionally not a raw-value bar chart because process channels have unlike units. The gallery table shows values, units, and missing markers."


def plot_fallback(d: dict[str,np.ndarray], out: Path, title: str) -> str:
    fig,ax=plt.subplots(figsize=(10,4)); ax.axis("off"); lines=[title,"No specialized visualizer applies; exported arrays:"]+[f"{k}: {list(v.shape)} {v.dtype}" for k,v in d.items()]
    ax.text(.02,.96,"\n".join(lines),va="top",family="monospace"); save(fig,out); return "Array inventory preview."


def table_html(d: dict[str,np.ndarray]) -> str:
    if "process_x" not in d: return ""
    names=strings(d["process_feature_names"]); runs=strings(d["run_id"]); x=d["process_x"]; ok=d["process_feature_valid"]
    rows="".join("<tr><td>"+html.escape(name)+"</td>"+"".join(f"<td>{html.escape(format(float(x[i,j]),'.7g')) if ok[i,j] else '<span class=\"muted\">missing</span>'}</td>" for i in range(len(runs)))+"</tr>" for j,name in enumerate(names))
    head="".join(f"<th>{html.escape(run)}</th>" for run in runs)
    geo_names=strings(d["geometry_feature_names"]); geo=d["geometry_x"][0]
    grows="".join(f"<tr><td>{html.escape(n)}</td><td>{float(v):.7g}</td></tr>" for n,v in zip(geo_names,geo))
    return f"<h2>Raw process vectors (channel [unit])</h2><table><thead><tr><th>channel [unit]</th>{head}</tr></thead><tbody>{rows}</tbody></table><h2>Fixed CAD geometry summary</h2><table><thead><tr><th>feature</th><th>value</th></tr></thead><tbody>{grows}</tbody></table><p class=muted>{html.escape(str(scalar(d.get('split_metadata',''))))}</p>"


def render_one(family: str, npz: Path, visual_dir: Path, family_dir: Path) -> tuple[dict[str,Any], str]:
    with np.load(npz,allow_pickle=False) as raw: d={k:raw[k] for k in raw.files}
    preview=visual_dir/f"{npz.stem}.png"; title=f"{family}: {npz.name}"
    if family=="images": description=plot_images(d,preview,title)
    elif family=="points": description=plot_points(d,preview,title)
    elif family=="voxels": description=plot_sparse_voxels(d,preview,title) if "voxel_coords" in d else plot_dense_voxels(d,preview,title)
    elif family=="graphs": description=plot_graph(d,preview,title,graph_requested(family_dir,npz.stem))
    elif family=="brep": description=plot_brep(d,preview,title)
    elif family=="tabular": description=plot_tabular(d,preview,title)
    elif family=="paired" and preview.exists():
        description="STEP/XLS inputs matched to the fused PLY outcome and defect CSV. Includes 2×17 heater maps, provisional CAD overlays, scan residuals, and weak area labels."
    else: description=plot_fallback(d,preview,title)
    schema=npz.with_suffix(".schema.json")
    entry={"source_npz":npz.name,"schema":schema.name if schema.exists() else None,"preview":"visualizations/"+preview.name,"description":description,"arrays":{k:{"shape":list(v.shape),"dtype":str(v.dtype)} for k,v in d.items()}}
    return entry, table_html(d) if family=="tabular" else ""


def family_page(root: Path, family: str, entries: list[dict[str,Any]], tabular_detail: str) -> None:
    directory=root/"data_processed"/family; cards=[]
    for e in entries:
        links=f'<a href="{html.escape(e["source_npz"])}">raw NPZ</a>'
        if e["schema"]: links+=f' · <a href="{html.escape(e["schema"])}">schema</a>'
        cards.append(f'<section class="card"><h2>{html.escape(e["source_npz"])}</h2><p>{html.escape(e["description"])}</p><p>{links}</p><img src="{html.escape(e["preview"])}" alt="preview for {html.escape(e["source_npz"])}"></section>')
    extra = '<p><a href="../../results/graphs/viewer.html">Open existing interactive graph comparison</a></p>' if family=="graphs" else ""
    text=f'<!doctype html><meta charset="utf-8"><title>{family} gallery</title><style>{CSS}</style><main><p><a href="../index.html">← all processed families</a></p><h1>{html.escape(family)} gallery</h1>{extra}{"".join(cards)}{tabular_detail}</main>'
    (directory/"index.html").write_text(text,encoding="utf-8")


def root_page(root: Path, summary: dict[str,int]) -> None:
    preferred={"images":"cad_multiview_256.png","voxels":"cad_surface_sparse_64.png","points":"cad_surface_points.png","graphs":"graph_12000.png","brep":"faces_uv.png","tabular":"process_and_geometry.png","paired":"18h01_cad_ply_targets.png"}
    cards=""
    for family,count in summary.items():
        thumb=preferred.get(family)
        thumb_path=root/"data_processed"/family/"visualizations"/str(thumb)
        image=(f'<a href="{family}/index.html"><img class="thumb" src="{family}/visualizations/{html.escape(str(thumb))}" alt="{html.escape(family)} preview thumbnail"></a>' if thumb_path.exists() else "")
        cards += f'<section class="card"><h2><a href="{family}/index.html">{family}</a></h2><p>{count} NPZ preview(s), each with raw artifact and schema links.</p>{image}</section>'
    text=f'<!doctype html><meta charset="utf-8"><title>Processed data galleries</title><style>{CSS}</style><main><h1>Processed data galleries</h1><p class="muted">Static local scientific artifact reports. Previews are descriptive only: no defect targets, untrained predictions, or fabricated physical fields.</p>{cards}</main>'
    (root/"data_processed"/"index.html").write_text(text,encoding="utf-8")


def generate(root: Path) -> dict[str,int]:
    processed=root/"data_processed"; summary: dict[str,int]={}
    for family in FAMILIES:
        directory=processed/family; directory.mkdir(parents=True,exist_ok=True); visual=directory/"visualizations"; visual.mkdir(exist_ok=True)
        entries=[]; tabular_detail=""
        for npz in sorted(directory.glob("*.npz")):
            entry, detail=render_one(family,npz,visual,directory); entries.append(entry); tabular_detail += detail
        (visual/"manifest.json").write_text(json.dumps({"family":family,"entries":entries},ensure_ascii=False,indent=2),encoding="utf-8")
        family_page(root,family,entries,tabular_detail); summary[family]=len(entries)
    root_page(root,summary); return summary


def main() -> None:
    ap=argparse.ArgumentParser(description="Generate local visual galleries for every data_processed NPZ.")
    ap.add_argument("--root",type=Path,default=root_default(),help="repository root (default: inferred from this script)")
    args=ap.parse_args(); print(json.dumps(generate(args.root),ensure_ascii=False))


if __name__=="__main__": main()