Newer
Older
SmTIAS-LightSim / scripts / build_blend_bottom_panel.py
"""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())