Newer
Older
SmTIAS-LightSim / scripts / build_blend_bottom_panel.py
"""TEST_09 実験5: 単灯 blendbsdf スイープ(バグ修正後再検証)の解析とパネル生成.

sweep レンダ(output/blend_bottom/{job}/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)

対象灯は --side で選ぶ(既定 bottom=フルグリッド 4×3.left / right は時間短縮のため
縮小グリッド w{0.05, 0.3}×α{0.1, 0.6} + Phase 2 の w0.6_a0.1).

生成するファイル:
    output/blend_bottom/summary{_side}.csv                  -- 全ジョブの形指標
    output/_slide_panels/blend_{side}_shape_roi.png         -- 輝度マップ + 指標パネル(ROI 揃え)

使い方:
    uv run python scripts/build_blend_bottom_panel.py [--side {bottom,left,right}]

前提:
    - output/blend_bottom/ に対象ジョブの luminance.npy が存在すること
    - area diffuse 基準(bottom: output/depth_bottom/d16,left/right: output/sanity_fix_1024/)
    - lightsim/luminance_maps/{side}_*.npy(実機単灯実測,ROI 適用済み)が存在すること
"""

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 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"
PANEL_DIR = ROOT / "output" / "_slide_panels"
REAL_DIR = ROOT / "lightsim" / "luminance_maps"

# 灯ごとの設定(configs/_tmp_blend_bottom/ のジョブ命名・グリッドと整合)
# baseline: area diffuse(純 Lambertian)の同条件出力.
# grid: (weights, alphas).left/right は縮小グリッド+Phase 2 の w0.6_a0.1 のみ存在し,
#       欠けセルはパネル上 [MISSING] 表示になる.
SIDES: dict[str, dict] = {
    "bottom": {
        "baseline": ROOT / "output" / "depth_bottom" / "d16" / "luminance.npy",
        "real": REAL_DIR / "bottom_SmTIAS_QM_20260615_141420.npy",
        "job_prefix": "",
        "weights": [0.05, 0.15, 0.3, 0.6],
        "alphas": [0.1, 0.3, 0.6],
    },
    "left": {
        "baseline": ROOT / "output" / "sanity_fix_1024" / "left" / "luminance.npy",
        "real": REAL_DIR / "left_SmTIAS_QM_20260615_150300.npy",
        "job_prefix": "left_",
        "weights": [0.05, 0.3, 0.6],
        "alphas": [0.1, 0.6],
    },
    "right": {
        "baseline": ROOT / "output" / "sanity_fix_1024" / "right" / "luminance.npy",
        "real": REAL_DIR / "right_SmTIAS_QM_20260615_170948.npy",
        "job_prefix": "right_",
        "weights": [0.05, 0.3, 0.6],
        "alphas": [0.1, 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(
    side: str,
    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)) × len(alphas) グリッドの形比較パネルを生成・保存する.

    先頭行に基準(area diffuse)と実機を置き,以下 weight × alpha グリッド.
    形の比較のため全画像を mean=1 に正規化し共有スケールで描画する.

    Args:
        side: 対象灯("bottom" / "left" / "right").グリッド定義は SIDES から引く.
        job_rois: ジョブ名 → ROI 輝度配列.
        job_metrics: ジョブ名 → 形指標辞書.
        baseline_roi: area diffuse 基準の ROI 輝度配列.
        baseline_metrics: 基準の形指標辞書.
        real: 実機実測輝度配列.
        real_metrics: 実機の形指標辞書(NCC は 1.0).
        output_path: 保存先 PNG パス.
    """
    cfg = SIDES[side]
    weights, alphas, prefix = cfg["weights"], cfg["alphas"], cfg["job_prefix"]

    normed = {k: v / v.mean() for k, v in job_rois.items()}
    baseline_n = baseline_roi / baseline_roi.mean()
    real_n = real / real.mean()

    n_rows = 1 + len(weights)
    n_cols = max(len(alphas), 2)
    fig, axes = plt.subplots(
        n_rows, n_cols, figsize=(n_cols * 4.2, n_rows * 4.0), constrained_layout=True
    )

    # 各セルを自身の 1–99 パーセンタイルで独立ストレッチする(形の比較専用パネル).
    # 共有スケールだとレンジの広いセルに引きずられ個々の形のコントラストが薄まる.
    # 厳密 min–max でなくパーセンタイルなのは,裾の外れ画素(実機は裾 2% が
    # カラーレンジの約 4 割を占有,シミュは MC ノイズ)にレンジを食われないため.
    def _show(ax: plt.Axes, img: np.ndarray | None, title: str) -> None:
        if img is not None:
            cell_vmin = float(np.percentile(img, VMIN_PERCENTILE))
            cell_vmax = float(np.percentile(img, VMAX_PERCENTILE))
            ax.imshow(img, cmap="hot", vmin=cell_vmin, vmax=cell_vmax, aspect="equal")
        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 / 実機 / 空
    _show(axes[0][0], baseline_n, f"area diffuse(基準)\n{_fmt(baseline_metrics)}")
    _show(axes[0][1], real_n, f"実機 {side}\n{_fmt(real_metrics)}")
    for col_idx in range(2, n_cols):
        _show(axes[0][col_idx], None, "")

    for row_idx, weight in enumerate(weights, start=1):
        for col_idx, alpha in enumerate(alphas):
            job = f"{prefix}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]")
        for col_idx in range(len(alphas), n_cols):
            _show(axes[row_idx][col_idx], None, "")

    fig.suptitle(
        f"TEST_09 実験5 — {side} 単灯 blend BSDF 再検証(バグ修正後, spp=1024, d=16, ROI 揃え)\n"
        f"各セル独立 1–99% 正規化(形の比較専用・明るさの絶対比較は不可)/ "
        f"実機目標: peak/mean={real_metrics['peak_mean']:.3f}, "
        f"top/bot={real_metrics['top_bot']:.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).
    """
    parser = argparse.ArgumentParser(description="単灯 blend BSDF スイープの形指標解析とパネル生成")
    parser.add_argument(
        "--side",
        choices=sorted(SIDES),
        default="bottom",
        help="対象灯(既定: bottom)",
    )
    args = parser.parse_args()
    side = args.side
    cfg = SIDES[side]
    weights, alphas, prefix = cfg["weights"], cfg["alphas"], cfg["job_prefix"]

    real_npy: Path = cfg["real"]
    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(cfg["baseline"])
    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"{prefix}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(bottom は既存名 summary.csv を維持)
    summary_csv = SWEEP_BASE / ("summary.csv" if side == "bottom" else f"summary_{side}.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", "", "", *(f"{baseline_metrics[k]:.4f}" for k in
                                       ("peak_mean", "top_bot", "ncc"))]
        )
        for weight in weights:
            for alpha in alphas:
                job = f"{prefix}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':<22}{'peak/mean':>10}{'top/bot':>10}{'NCC':>8}")
    print(f"{'area diffuse':<22}{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"{prefix}w{weight}_a{alpha}"
            if job in job_metrics:
                m = job_metrics[job]
                print(f"{job:<22}{m['peak_mean']:>10.3f}{m['top_bot']:>10.3f}{m['ncc']:>8.3f}")
    print(f"{'実機 ' + side:<22}{real_metrics['peak_mean']:>10.3f}"
          f"{real_metrics['top_bot']:>10.3f}{'-':>8}")

    panel_path = PANEL_DIR / f"blend_{side}_shape_roi.png"
    _build_panel(
        side, job_rois, job_metrics, baseline_roi, baseline_metrics, real, real_metrics,
        panel_path,
    )
    return 0


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