"""ρ×ω×α スイープの輝度マップを「目で判断する」ための可視化. 指標（CoV / NCC 等）は分布を集約した量であり，形そのものの判断材料にはならない （TEST_03 の主題）．本スクリプトは summary.csv の順位ではなく **輝度マップそのもの** を並べ，どのパラメータで形がどう変わるか／変わらないかを目視できるようにする． 生成物: 1. bsdf_grid_contact_sheet.png -- 全 60 条件 + 実機のコンタクトシート （行 = ρ，列 = (ω, α)．全セル共通カラースケールなので差がそのまま見える） 2. bsdf_grid_rho_series.png -- ω=0 純 diffuse の ρ 系列 + 実機． 上段 = 輝度マップ，下段 = 実機との相対残差マップ [%] 3. bsdf_grid_top_residual.png -- スコア上位条件の残差マップ（発散カラーマップ） 残差の定義: 平均で正規化したうえで (sim - real) / mean(real) [%]． 正（赤）= シミュが明るすぎる領域，負（青）= 暗すぎる領域． 実行方法: uv run python scripts/build_bsdf_grid_visual.py [--top 4] """ from __future__ import annotations import argparse import csv import sys from pathlib import Path import matplotlib import matplotlib.pyplot as plt import numpy as np from scipy.ndimage import gaussian_filter, median_filter from skimage.transform import resize sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src")) from sm_tias_lightsim.evaluation.roi_presets import SIM_ROI, extract_roi_xywh # noqa: E402 matplotlib.rcParams["font.family"] = ["Meiryo", "Yu Gothic", "DejaVu Sans"] ROOT = Path(__file__).resolve().parents[1] SUMMARY_CSV = ROOT / "output" / "bsdf_grid_fix" / "summary.csv" OUTPUT_BASE = ROOT / "output" / "bsdf_grid_fix" REFERENCE_DIR = ROOT / "output" / "full_base_fix" REAL_FULL = ROOT / "lightsim" / "luminance_maps" / "full_MiniTIAS_QM_20260601_124500.npy" PANEL_DIR = ROOT / "output" / "_slide_panels" RESIDUAL_LIMIT = 12.0 # 残差カラーバーの上下限 [%] def load_rows() -> list[dict]: """summary.csv を数値化して読み込む.""" if not SUMMARY_CSV.exists(): raise SystemExit(f"エラー: {SUMMARY_CSV} がありません（スイープ未実行）") rows = [] with SUMMARY_CSV.open(encoding="utf-8") as fh: for raw in csv.DictReader(fh): row: dict = {"job": raw["job"], "stage": raw["stage"]} for key, value in raw.items(): if key not in ("job", "stage"): row[key] = float(value) if value not in ("", None) else float("nan") rows.append(row) return rows def _smooth(arr: np.ndarray, is_sim: bool) -> np.ndarray: """他パネルと同一の平滑化（median5 → gauss σ4）と ROI 抽出.""" if arr.ndim == 3: arr = 0.2126 * arr[..., 0] + 0.7152 * arr[..., 1] + 0.0722 * arr[..., 2] arr = arr.astype(np.float64) if is_sim: arr = extract_roi_xywh(arr, SIM_ROI) return gaussian_filter(median_filter(arr, size=5), sigma=4) def _job_map(job: str) -> np.ndarray: """ジョブ名から平滑化・mean 正規化済みの輝度マップを返す.""" path = OUTPUT_BASE / job / "luminance.npy" if not path.exists(): path = REFERENCE_DIR / "luminance.npy" img = _smooth(np.load(path), True) return img / img.mean() def _real_map(shape: tuple[int, int]) -> np.ndarray: """実機マップを指定 shape に合わせ mean 正規化して返す.""" real = _smooth(np.load(REAL_FULL), False) real = resize(real, shape, anti_aliasing=True) return real / real.mean() def _residual(sim_norm: np.ndarray, real_norm: np.ndarray) -> np.ndarray: """mean 正規化済みどうしの相対残差 [%]（正＝シミュが明るすぎ）.""" return 100.0 * (sim_norm - real_norm) def contact_sheet(rows: list[dict]) -> Path: """全条件を 行=ρ × 列=(ω, α) で並べた一覧（共通カラースケール）.""" rho_values = sorted({r["rho"] for r in rows}) combos = sorted({(r["omega"], r["alpha"]) for r in rows}) lookup = {(r["rho"], r["omega"], r["alpha"]): r for r in rows} ncols = len(combos) + 1 # 右端に実機列 fig, axes = plt.subplots(len(rho_values), ncols, figsize=(1.75 * ncols, 1.95 * len(rho_values)), constrained_layout=True) sample = _job_map(rows[0]["job"]) real = _real_map(sample.shape) vmin, vmax = np.percentile(real, [1, 99]) for row_idx, rho in enumerate(rho_values): for col_idx, (omega, alpha) in enumerate(combos): ax = axes[row_idx][col_idx] entry = lookup.get((rho, omega, alpha)) if entry is None: ax.axis("off") continue ax.imshow(_job_map(entry["job"]), cmap="hot", vmin=vmin, vmax=vmax) ax.set_title(f"NCC={entry['ncc_smooth']:.3f}", fontsize=6.5, pad=2) ax.set_xticks([]) ax.set_yticks([]) if col_idx == 0: ax.set_ylabel(f"ρ={rho:g}", fontsize=9) if row_idx == 0: label = "ω=0" if omega == 0 else f"ω={omega:g}\nα={alpha:g}" ax.text(0.5, 1.35, label, transform=ax.transAxes, ha="center", va="bottom", fontsize=8) ax_real = axes[row_idx][-1] im = ax_real.imshow(real, cmap="hot", vmin=vmin, vmax=vmax) ax_real.set_xticks([]) ax_real.set_yticks([]) ax_real.set_title("実機", fontsize=7, pad=2) if row_idx == 0: ax_real.text(0.5, 1.35, "参照", transform=ax_real.transAxes, ha="center", va="bottom", fontsize=8) fig.colorbar(im, ax=axes, shrink=0.5, label="mean 正規化輝度（全セル共通スケール）") fig.suptitle("全 60 条件の輝度マップ一覧 — 行=ρ，列=(ω, α)，右端=実機" "（全セル共通カラースケール・mean 正規化）", fontsize=13) PANEL_DIR.mkdir(parents=True, exist_ok=True) path = PANEL_DIR / "bsdf_grid_contact_sheet.png" fig.savefig(path, dpi=140, bbox_inches="tight") plt.close(fig) return path def rho_series(rows: list[dict]) -> Path: """ω=0 の ρ 系列を輝度マップ（上段）と残差マップ（下段）で並べる.""" series = sorted([r for r in rows if r["omega"] == 0], key=lambda r: r["rho"]) sample = _job_map(series[0]["job"]) real = _real_map(sample.shape) vmin, vmax = np.percentile(real, [1, 99]) ncols = len(series) + 1 fig, axes = plt.subplots(2, ncols, figsize=(2.55 * ncols, 5.9), constrained_layout=True) for col, entry in enumerate(series): sim = _job_map(entry["job"]) axes[0][col].imshow(sim, cmap="hot", vmin=vmin, vmax=vmax) axes[0][col].set_title(f"ρ={entry['rho']:g}\nCoV={entry['cov']:.4f} " f"NCC={entry['ncc_smooth']:.3f}", fontsize=9) res = axes[1][col].imshow(_residual(sim, real), cmap="RdBu_r", vmin=-RESIDUAL_LIMIT, vmax=RESIDUAL_LIMIT) rms = float(np.sqrt(np.mean(_residual(sim, real) ** 2))) axes[1][col].set_title(f"残差 RMS={rms:.2f}%", fontsize=9) top = axes[0][-1].imshow(real, cmap="hot", vmin=vmin, vmax=vmax) axes[0][-1].set_title("実機 全灯\nCoV=0.0695", fontsize=9) axes[1][-1].axis("off") for ax in axes.ravel(): if ax.axison: ax.set_xticks([]) ax.set_yticks([]) fig.colorbar(top, ax=axes[0], shrink=0.8, label="mean 正規化輝度") fig.colorbar(res, ax=axes[1], shrink=0.8, label="相対残差 [%]（赤=シミュが明るい）") fig.suptitle("純 diffuse（ω=0）の ρ 系列 — 上段: 輝度マップ / 下段: 実機との残差", fontsize=13) PANEL_DIR.mkdir(parents=True, exist_ok=True) path = PANEL_DIR / "bsdf_grid_rho_series.png" fig.savefig(path, dpi=150, bbox_inches="tight") plt.close(fig) return path def top_residual(rows: list[dict], top: int) -> Path: """スコア上位条件の残差マップを並べる（形の判断はここが主戦場）.""" best = sorted(rows, key=lambda r: r["score"])[:top] sample = _job_map(best[0]["job"]) real = _real_map(sample.shape) fig, axes = plt.subplots(1, top, figsize=(3.7 * top, 4.5), constrained_layout=True) for ax, entry in zip(np.atleast_1d(axes), best, strict=True): residual = _residual(_job_map(entry["job"]), real) im = ax.imshow(residual, cmap="RdBu_r", vmin=-RESIDUAL_LIMIT, vmax=RESIDUAL_LIMIT) rms = float(np.sqrt(np.mean(residual ** 2))) p95 = float(np.percentile(np.abs(residual), 95)) ax.set_title(f"{entry['job']}\nRMS={rms:.2f}% P95={p95:.2f}% " f"NCC={entry['ncc_smooth']:.3f}", fontsize=10) ax.set_xticks([]) ax.set_yticks([]) fig.colorbar(im, ax=axes, shrink=0.8, label="相対残差 [%]（赤=シミュが明るい / 青=暗い）") fig.suptitle("スコア上位条件の残差マップ — 赤青の構造が「材質では埋まらない形の差」", fontsize=13) PANEL_DIR.mkdir(parents=True, exist_ok=True) path = PANEL_DIR / "bsdf_grid_top_residual.png" fig.savefig(path, dpi=150, bbox_inches="tight") plt.close(fig) return path def main() -> int: """3 種の可視化を生成するメインエントリポイント.""" parser = argparse.ArgumentParser(description="ρ×ω×α スイープの輝度マップ可視化") parser.add_argument("--top", type=int, default=4, help="残差マップに並べる上位条件数") args = parser.parse_args() rows = load_rows() print(f"読み込み: {len(rows)} 条件") print(f"一覧: {contact_sheet(rows)}") print(f"ρ 系列: {rho_series(rows)}") print(f"残差: {top_residual(rows, args.top)}") return 0 if __name__ == "__main__": sys.exit(main())