diff --git a/scripts/build_blend_bottom_panel.py b/scripts/build_blend_bottom_panel.py new file mode 100644 index 0000000..e55e21b --- /dev/null +++ b/scripts/build_blend_bottom_panel.py @@ -0,0 +1,295 @@ +"""TEST_09 実験5: bottom 単灯 blendbsdf スイープ（バグ修正後再検証）の解析とパネル生成. + +sweep レンダ（output/blend_bottom/w{W}_a{A}/luminance.npy）の完了後に実行すること. +旧 blend 実験（scripts/sweep_blend_bsdf.py，3 灯フルシーン・絶対指標）は面光源スケール軸 +バグ込みの結果のため，本スクリプトは TEST_09 の「形」評価に統一して再検証する: + +- peak/mean(ROI) = max / mean（実機 bottom ≈ 1.218） +- top/bot(ROI) = 上半分 mean / 下半分 mean（実機 bottom ≈ 0.931） +- NCC(ROI, 対実機) = シミュ ROI を実機解像度へ resize して正規化相関（area diffuse 基準 0.720） + +生成するファイル: + output/blend_bottom/summary.csv -- 全ジョブの形指標 + output/_slide_panels/blend_bottom_shape_roi.png -- 輝度マップ + 指標パネル（ROI 揃え） + +使い方: + uv run python scripts/build_blend_bottom_panel.py + +前提: + - output/blend_bottom/w{W}_a{A}/luminance.npy が存在すること + - output/depth_bottom/d16/luminance.npy（area diffuse 基準）が存在すること + - lightsim/luminance_maps/bottom_*.npy（実機単灯実測，ROI 適用済み）が存在すること +""" + +from __future__ import annotations + +import csv +import sys +from pathlib import Path + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +from skimage.transform import resize + +from sm_tias_lightsim.evaluation.roi_presets import SIM_ROI, extract_roi_xywh +from sm_tias_lightsim.evaluation.shape_metrics import calc_ncc + +# 日本語ラベル用フォント（Windows 標準．無ければ DejaVu にフォールバック） +matplotlib.rcParams["font.family"] = ["Meiryo", "Yu Gothic", "DejaVu Sans"] + +# --------------------------------------------------------------------------- +# 定数 +# --------------------------------------------------------------------------- + +ROOT = Path(__file__).resolve().parents[1] + +SWEEP_BASE = ROOT / "output" / "blend_bottom" +BASELINE_NPY = ROOT / "output" / "depth_bottom" / "d16" / "luminance.npy" +REAL_NPY = ROOT / "lightsim" / "luminance_maps" / "bottom_SmTIAS_QM_20260615_141420.npy" +PANEL_PATH = ROOT / "output" / "_slide_panels" / "blend_bottom_shape_roi.png" +SUMMARY_CSV = SWEEP_BASE / "summary.csv" + +# スイープグリッド（configs/_tmp_blend_bottom/ と整合） +WEIGHTS = [0.05, 0.15, 0.3, 0.6] +ALPHAS = [0.1, 0.3, 0.6] + +# 共有 colorbar のパーセンタイル（mean=1 正規化後） +VMIN_PERCENTILE = 1.0 +VMAX_PERCENTILE = 99.0 + + +# --------------------------------------------------------------------------- +# 指標計算（TEST_09 実験1〜3 と同一定義） +# --------------------------------------------------------------------------- + + +def load_sim_roi(npy_path: Path) -> np.ndarray: + """シミュ luminance.npy を読み込み SIM_ROI で切り出す. + + Args: + npy_path: luminance.npy のパス（2D 輝度 or 3D RGB）. + + Returns: + ROI 切り出し済み輝度配列（H × W, float64）. + + Raises: + FileNotFoundError: ファイルが存在しない場合. + """ + if not npy_path.exists(): + msg = f"luminance.npy が見つかりません: {npy_path}" + raise FileNotFoundError(msg) + arr = np.load(npy_path) + if arr.ndim == 3: + arr = 0.2126 * arr[..., 0] + 0.7152 * arr[..., 1] + 0.0722 * arr[..., 2] + return extract_roi_xywh(arr, SIM_ROI) + + +def calc_shape_metrics(roi: np.ndarray, real: np.ndarray) -> dict[str, float]: + """TEST_09 の 3 形指標（peak/mean・top/bot・NCC 対実機）を算出する. + + NCC はシミュ ROI を実機解像度へ resize してから計算する + （TEST_09 掲載値と同一の手順: d16 で 0.7207 を再現確認済み）. + + Args: + roi: シミュ ROI 輝度配列（H × W）. + real: 実機実測輝度配列（ROI 適用済み，H' × W'）. + + Returns: + {"peak_mean", "top_bot", "ncc"} の辞書. + """ + h = roi.shape[0] + sim_resized = resize(roi, real.shape, anti_aliasing=True) + return { + "peak_mean": float(roi.max() / roi.mean()), + "top_bot": float(roi[: h // 2].mean() / roi[h // 2 :].mean()), + "ncc": float(calc_ncc(sim_resized, real)), + } + + +# --------------------------------------------------------------------------- +# パネル描画 +# --------------------------------------------------------------------------- + + +def _build_panel( + job_rois: dict[str, np.ndarray], + job_metrics: dict[str, dict[str, float]], + baseline_roi: np.ndarray, + baseline_metrics: dict[str, float], + real: np.ndarray, + real_metrics: dict[str, float], + output_path: Path, +) -> None: + """(1 + len(WEIGHTS)) × 3 グリッドの形比較パネルを生成・保存する. + + 先頭行に基準（area diffuse d16）と実機を置き，以下 weight × alpha グリッド. + 形の比較のため全画像を mean=1 に正規化し共有スケールで描画する. + + Args: + job_rois: ジョブ名 → ROI 輝度配列. + job_metrics: ジョブ名 → 形指標辞書. + baseline_roi: area diffuse 基準の ROI 輝度配列. + baseline_metrics: 基準の形指標辞書. + real: 実機実測輝度配列. + real_metrics: 実機の形指標辞書（NCC は 1.0）. + output_path: 保存先 PNG パス. + """ + normed = {k: v / v.mean() for k, v in job_rois.items()} + baseline_n = baseline_roi / baseline_roi.mean() + real_n = real / real.mean() + + all_vals = np.concatenate( + [baseline_n.ravel(), real_n.ravel()] + [v.ravel() for v in normed.values()] + ) + vmin = float(np.percentile(all_vals, VMIN_PERCENTILE)) + vmax = float(np.percentile(all_vals, VMAX_PERCENTILE)) + + n_rows = 1 + len(WEIGHTS) + n_cols = len(ALPHAS) + fig, axes = plt.subplots( + n_rows, n_cols, figsize=(n_cols * 4.2, n_rows * 4.0), constrained_layout=True + ) + + def _show(ax: plt.Axes, img: np.ndarray | None, title: str) -> None: + if img is not None: + ax.imshow(img, cmap="hot", vmin=vmin, vmax=vmax, aspect="auto") + else: + ax.set_facecolor("0.15") + ax.set_title(title, fontsize=9, pad=4) + ax.set_axis_off() + + def _fmt(m: dict[str, float]) -> str: + return f"p/m={m['peak_mean']:.3f} t/b={m['top_bot']:.3f} NCC={m['ncc']:.3f}" + + # 先頭行: 基準 area diffuse / 実機 / 空 + # 実機はシミュよりレンジが狭く共有スケールでは形が潰れるため， + # 実機セルのみ自身の 1–99 パーセンタイルで正規化して形を可視化する + # （厳密 min–max では裾 2% の画素にカラーレンジの約 4 割を食われる． + # 共有スケール側と同じパーセンタイル流儀に揃える） + _show(axes[0][0], baseline_n, f"area diffuse (d16 基準)\n{_fmt(baseline_metrics)}") + real_vmin = float(np.percentile(real_n, VMIN_PERCENTILE)) + real_vmax = float(np.percentile(real_n, VMAX_PERCENTILE)) + axes[0][1].imshow(real_n, cmap="hot", vmin=real_vmin, vmax=real_vmax, aspect="auto") + axes[0][1].set_title( + f"実機 bottom（1–99% 正規化・独立スケール）\n{_fmt(real_metrics)}", fontsize=9, pad=4 + ) + axes[0][1].set_axis_off() + _show(axes[0][2], None, "") + + for row_idx, weight in enumerate(WEIGHTS, start=1): + for col_idx, alpha in enumerate(ALPHAS): + job = f"w{weight}_a{alpha}" + if job in normed: + _show( + axes[row_idx][col_idx], + normed[job], + f"w={weight}, α={alpha}\n{_fmt(job_metrics[job])}", + ) + else: + _show(axes[row_idx][col_idx], None, f"w={weight}, α={alpha}\n[MISSING]") + + fig.colorbar( + plt.cm.ScalarMappable(norm=plt.Normalize(vmin=vmin, vmax=vmax), cmap="hot"), + ax=axes, + fraction=0.015, + pad=0.02, + label="Luminance (mean=1 normalized)", + ) + fig.suptitle( + "TEST_09 実験5 — bottom 単灯 blend BSDF 再検証（バグ修正後, spp=1024, d=16, ROI 揃え）\n" + f"実機目標: peak/mean={real_metrics['peak_mean']:.3f}, " + f"top/bot={real_metrics['top_bot']:.3f} / " + f"area diffuse 基準 NCC={baseline_metrics['ncc']:.3f}", + fontsize=11, + y=1.03, + ) + + output_path.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(output_path, dpi=150, bbox_inches="tight") + plt.close(fig) + print(f"パネル保存: {output_path}") + + +# --------------------------------------------------------------------------- +# メイン処理 +# --------------------------------------------------------------------------- + + +def main() -> int: + """解析とパネル生成のメインエントリポイント. + + Returns: + 終了コード（成功なら 0，入力不足なら 1）. + """ + if not REAL_NPY.exists(): + print(f"エラー: 実機 npy が見つかりません: {REAL_NPY}", file=sys.stderr) + return 1 + real = np.load(REAL_NPY).astype(np.float64) + real_metrics = { + "peak_mean": float(real.max() / real.mean()), + "top_bot": float(real[: real.shape[0] // 2].mean() / real[real.shape[0] // 2 :].mean()), + "ncc": 1.0, + } + + baseline_roi = load_sim_roi(BASELINE_NPY) + baseline_metrics = calc_shape_metrics(baseline_roi, real) + + job_rois: dict[str, np.ndarray] = {} + job_metrics: dict[str, dict[str, float]] = {} + for weight in WEIGHTS: + for alpha in ALPHAS: + job = f"w{weight}_a{alpha}" + try: + roi = load_sim_roi(SWEEP_BASE / job / "luminance.npy") + except FileNotFoundError: + print(f" 警告: {job} 未完了（スキップ）") + continue + job_rois[job] = roi + job_metrics[job] = calc_shape_metrics(roi, real) + + if not job_rois: + print("エラー: 有効なジョブが 1 件もありません.", file=sys.stderr) + return 1 + + # summary.csv + SUMMARY_CSV.parent.mkdir(parents=True, exist_ok=True) + with SUMMARY_CSV.open("w", encoding="utf-8", newline="") as f: + writer = csv.writer(f) + writer.writerow(["job", "weight", "alpha", "peak_mean", "top_bot", "ncc"]) + writer.writerow( + ["area_diffuse_d16", "", "", *(f"{baseline_metrics[k]:.4f}" for k in + ("peak_mean", "top_bot", "ncc"))] + ) + for weight in WEIGHTS: + for alpha in ALPHAS: + job = f"w{weight}_a{alpha}" + if job in job_metrics: + m = job_metrics[job] + writer.writerow( + [job, weight, alpha, f"{m['peak_mean']:.4f}", + f"{m['top_bot']:.4f}", f"{m['ncc']:.4f}"] + ) + print(f"summary 保存: {SUMMARY_CSV}") + + # コンソール表 + print(f"\n{'job':<16}{'peak/mean':>10}{'top/bot':>10}{'NCC':>8}") + print(f"{'area diffuse':<16}{baseline_metrics['peak_mean']:>10.3f}" + f"{baseline_metrics['top_bot']:>10.3f}{baseline_metrics['ncc']:>8.3f}") + for weight in WEIGHTS: + for alpha in ALPHAS: + job = f"w{weight}_a{alpha}" + if job in job_metrics: + m = job_metrics[job] + print(f"{job:<16}{m['peak_mean']:>10.3f}{m['top_bot']:>10.3f}{m['ncc']:>8.3f}") + print(f"{'実機 bottom':<16}{real_metrics['peak_mean']:>10.3f}" + f"{real_metrics['top_bot']:>10.3f}{'-':>8}") + + _build_panel( + job_rois, job_metrics, baseline_roi, baseline_metrics, real, real_metrics, PANEL_PATH + ) + return 0 + + +if __name__ == "__main__": + sys.exit(main())