Newer
Older
SmTIAS-LightSim / scripts / build_bsdf_grid_panel.py
"""ρ×ω×α 再スイープ(面法線修正後)の結果を可視化する.

`sweep_bsdf_grid_fix.py` が出す output/bsdf_grid_fix/summary.csv を読み,

  1. 応答曲面: 各指標を ρ の関数として (ω, α) 系列ごとに描く.
     実機の値を水平線で引き,どの ρ で交差するかを目で追えるようにする.
  2. 形比較パネル: スコア上位 N 条件と実機を並べる(他の TEST パネルと同じ流儀.
     mean 正規化 + 独立 1-99% 表示 + 同一平滑化).

出力:
    output/_slide_panels/bsdf_grid_response.png   -- 応答曲面(指標 × ρ)
    output/_slide_panels/bsdf_grid_top_shape.png  -- 上位条件 + 実機の形比較

実行方法:
    uv run python scripts/build_bsdf_grid_panel.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

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"

# 実機(全灯)の指標: sweep_bsdf_grid_fix.py と同じ出典
REAL_VALUES = {
    "cov": 0.0695,
    "center_periphery_ratio": 1.0911,
    "gradient_magnitude": 8.35,
    "radial_min_max_ratio": 0.860,
}
METRIC_LABELS = {
    "cov": "CoV",
    "center_periphery_ratio": "中心/周辺比",
    "gradient_magnitude": "勾配 [%]",
    "radial_min_max_ratio": "動径 min/max",
    "ncc_smooth": "NCC 対実機(平滑後)",
    "score": "総合スコア(小さいほど実機に近い)",
}


def load_rows() -> list[dict[str, float]]:
    """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[str, float] = {"job": raw["job"], "stage": raw["stage"]}
            for key, value in raw.items():
                if key in ("job", "stage"):
                    continue
                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_luminance(job: str) -> Path:
    """ジョブ名から luminance.npy のパスを返す(基準点は full_base_fix を指す)."""
    candidate = OUTPUT_BASE / job / "luminance.npy"
    if candidate.exists():
        return candidate
    return REFERENCE_DIR / "luminance.npy"


def plot_response(rows: list[dict[str, float]]) -> Path:
    """指標を ρ の関数として (ω, α) 系列ごとに描き,実機値を水平線で示す."""
    metrics = ["cov", "gradient_magnitude", "center_periphery_ratio",
               "radial_min_max_ratio", "ncc_smooth", "score"]
    series: dict[tuple[float, float], list[dict[str, float]]] = {}
    for row in rows:
        series.setdefault((row["omega"], row["alpha"]), []).append(row)

    fig, axes = plt.subplots(2, 3, figsize=(17, 9), constrained_layout=True)
    for ax, metric in zip(axes.ravel(), metrics, strict=True):
        for (omega, alpha), group in sorted(series.items()):
            group = sorted(group, key=lambda r: r["rho"])
            label = f"ω={omega:g}" if omega == 0 else f"ω={omega:g}, α={alpha:g}"
            ax.plot([r["rho"] for r in group], [r[metric] for r in group],
                    marker="o", markersize=4, linewidth=1.2, label=label)
        if metric in REAL_VALUES:
            ax.axhline(REAL_VALUES[metric], color="k", linestyle="--", linewidth=1.5,
                       label="実機" if metric == "cov" else None)
        ax.set_xlabel("ρ(拡散反射率)")
        ax.set_ylabel(METRIC_LABELS[metric])
        ax.set_title(METRIC_LABELS[metric])
        ax.grid(alpha=0.3)
    axes.ravel()[0].legend(fontsize=7, ncol=2)
    fig.suptitle("面法線修正後の ρ×ω×α 応答曲面 — 破線=実機全灯"
                 "(spp=1024, max_depth=16 固定)", fontsize=13)
    PANEL_DIR.mkdir(parents=True, exist_ok=True)
    path = PANEL_DIR / "bsdf_grid_response.png"
    fig.savefig(path, dpi=150, bbox_inches="tight")
    plt.close(fig)
    return path


def plot_top_shapes(rows: list[dict[str, float]], top: int) -> Path:
    """スコア上位 top 条件と実機の輝度マップを並べる."""
    best = sorted(rows, key=lambda r: r["score"])[:top]
    real = _smooth(np.load(REAL_FULL), False)
    cells = [(f"{r['job']}\nscore={r['score']:.4f}  CoV={r['cov']:.4f}  "
              f"NCC={r['ncc_smooth']:.3f}", _smooth(np.load(_job_luminance(r["job"])), True))
             for r in best]
    real_label = (f"実機 全灯\nCoV={REAL_VALUES['cov']:.4f}  "
                  f"勾配={REAL_VALUES['gradient_magnitude']}%")
    cells.append((real_label, real))

    fig, axes = plt.subplots(1, len(cells), figsize=(4.2 * len(cells), 4.6),
                             constrained_layout=True)
    for ax, (label, img) in zip(np.atleast_1d(axes), cells, strict=True):
        disp = img / img.mean()
        vmin, vmax = np.percentile(disp, [1, 99])
        ax.imshow(disp, cmap="hot", vmin=vmin, vmax=vmax)
        ax.set_title(label, fontsize=10)
        ax.axis("off")
    fig.suptitle("スコア上位条件 vs 実機 — mean 正規化・各セル独立 1–99% 表示(形の比較専用)",
                 fontsize=13)
    PANEL_DIR.mkdir(parents=True, exist_ok=True)
    path = PANEL_DIR / "bsdf_grid_top_shape.png"
    fig.savefig(path, dpi=150, bbox_inches="tight")
    plt.close(fig)
    return path


def main() -> int:
    """応答曲面と上位条件パネルを生成するメインエントリポイント."""
    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"応答曲面: {plot_response(rows)}")
    print(f"形比較:   {plot_top_shapes(rows, args.top)}")

    print("\n=== スコア上位 10 ===")
    header = (f"{'job':22s} {'score':>7s} {'CoV':>7s} {'C/P':>7s} "
              f"{'勾配%':>7s} {'動径':>6s} {'NCC':>6s} {'打切%':>6s}")
    print(header)
    for row in sorted(rows, key=lambda r: r["score"])[:10]:
        print(f"{row['job']:22s} {row['score']:7.4f} {row['cov']:7.4f} "
              f"{row['center_periphery_ratio']:7.4f} {row['gradient_magnitude']:7.2f} "
              f"{row['radial_min_max_ratio']:6.3f} {row['ncc_smooth']:6.3f} "
              f"{row['truncation_loss_pct']:6.1f}")
    print(f"{'実機':22s} {'-':>7s} {REAL_VALUES['cov']:7.4f} "
          f"{REAL_VALUES['center_periphery_ratio']:7.4f} "
          f"{REAL_VALUES['gradient_magnitude']:7.2f} "
          f"{REAL_VALUES['radial_min_max_ratio']:6.3f} {'1.000':>6s} {'-':>6s}")
    return 0


if __name__ == "__main__":
    sys.exit(main())