Newer
Older
SmTIAS-LightSim / scripts / build_loft_w1p4_panel.py
"""TEST_10 後側ロフト変形 w=1.4 実験の形比較パネル生成.

output/loft_w1p4/ のレンダ完了後に実行すること.
行 = bottom / left / right / full,列 = 元メッシュ基準 / 後側_w1p4 / 実機 の 4×3 グリッド.
シミュ ROI と実機マップは解像度が異なる(450×464 vs 約 1400×1409)ため,
全セルを共通解像度へ resize してから描画し,行・列の枠を完全に揃える.
各セルは独立 1–99 パーセンタイル正規化(形の比較専用・明るさの絶対比較は不可).

生成するファイル:
    output/loft_w1p4/summary.csv                    -- 全ジョブの形指標
    output/_slide_panels/loft_w1p4_shape_roi.png    -- 4×3 形比較パネル(ROI 揃え)

使い方:
    uv run python scripts/build_loft_w1p4_panel.py

前提:
    - output/loft_w1p4/{bottom,left,right,full,full_base}/luminance.npy が存在すること
    - 元メッシュ単灯基準 output/blend_bottom/{,left_,right_}w0.05_a0.1/luminance.npy
    - 実機マップ lightsim/luminance_maps/{side}_*.npy(単灯 3 種 + full.full は
      reference_measurement/quantitative/MiniTIAS_QM_20260601_124500.dng から
      単灯と同一手順(linear raw + LSC 順適用 + Rec.709 + REAL_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

matplotlib.rcParams["font.family"] = ["Meiryo", "Yu Gothic", "DejaVu Sans"]

ROOT = Path(__file__).resolve().parents[1]
PANEL_PATH = ROOT / "output" / "_slide_panels" / "loft_w1p4_shape_roi.png"
SUMMARY_CSV = ROOT / "output" / "loft_w1p4" / "summary.csv"

REAL = {
    "bottom": ROOT / "lightsim" / "luminance_maps" / "bottom_SmTIAS_QM_20260615_141420.npy",
    "left": ROOT / "lightsim" / "luminance_maps" / "left_SmTIAS_QM_20260615_150300.npy",
    "right": ROOT / "lightsim" / "luminance_maps" / "right_SmTIAS_QM_20260615_170948.npy",
    "full": ROOT / "lightsim" / "luminance_maps" / "full_MiniTIAS_QM_20260601_124500.npy",
}
BASE = {
    "bottom": ROOT / "output" / "blend_bottom" / "w0.05_a0.1" / "luminance.npy",
    "left": ROOT / "output" / "blend_bottom" / "left_w0.05_a0.1" / "luminance.npy",
    "right": ROOT / "output" / "blend_bottom" / "right_w0.05_a0.1" / "luminance.npy",
    "full": ROOT / "output" / "loft_w1p4" / "full_base" / "luminance.npy",
}

SIDES = ["bottom", "left", "right", "full"]
VERSIONS = ["base", "w1p4", "real"]
VERSION_LABELS = {"base": "元メッシュ基準", "w1p4": "後側_w1p4", "real": "実機"}

# 全灯ロフト重みスイープ(w2p0/w4p0 拡張分.ラベル → luminance.npy のディレクトリ)
FULL_SWEEP = [
    ("w0.5(元)", ROOT / "output" / "loft_w1p4" / "full_base"),
    ("w1.4", ROOT / "output" / "loft_w1p4" / "full"),
    ("w2.0", ROOT / "output" / "loft_w1p4" / "full_w2p0"),
    ("w4.0", ROOT / "output" / "loft_w1p4" / "full_w4p0"),
]
SWEEP_PANEL_PATH = ROOT / "output" / "_slide_panels" / "loft_full_sweep_shape_roi.png"
SWEEP_REL_PANEL_PATH = ROOT / "output" / "_slide_panels" / "loft_full_sweep_relative.png"

# 全セルの共通描画解像度(シミュ ROI と同程度.実機はダウンサンプル)
CELL_SHAPE = (450, 464)
VMIN_PERCENTILE = 1.0
VMAX_PERCENTILE = 99.0


def _load_luminance(path: Path) -> np.ndarray:
    """luminance.npy を読み込み 2D 輝度にする(3D RGB は Rec.709 変換)."""
    arr = np.load(path)
    if arr.ndim == 3:
        arr = 0.2126 * arr[..., 0] + 0.7152 * arr[..., 1] + 0.0722 * arr[..., 2]
    return arr.astype(np.float64)


def _metrics(img: np.ndarray, real: np.ndarray | None) -> dict[str, float]:
    """TEST_09/10 の形指標(peak/mean・top/bot・NCC 対実機)を算出する."""
    h = img.shape[0]
    m = {
        "peak_mean": float(img.max() / img.mean()),
        "top_bot": float(img[: h // 2].mean() / img[h // 2 :].mean()),
    }
    if real is not None:
        m["ncc"] = float(calc_ncc(resize(img, real.shape, anti_aliasing=True), real))
    else:
        m["ncc"] = float("nan")
    return m


def main() -> int:
    """パネルと summary CSV を生成するメインエントリポイント.

    Returns:
        終了コード(成功なら 0,入力不足なら 1).
    """
    cells: dict[tuple[str, str], np.ndarray] = {}
    rows: list[list[str]] = []
    for side in SIDES:
        real = _load_luminance(REAL[side]) if side in REAL else None
        for ver in ["base", "w1p4"]:
            if ver == "base":
                path = BASE[side]
            else:
                path = ROOT / "output" / "loft_w1p4" / side / "luminance.npy"
            if not path.exists():
                print(f"エラー: {path} が見つかりません", file=sys.stderr)
                return 1
            roi = extract_roi_xywh(_load_luminance(path), SIM_ROI)
            cells[(side, ver)] = roi
            m = _metrics(roi, real)
            rows.append(
                [side, ver, f"{m['peak_mean']:.4f}", f"{m['top_bot']:.4f}", f"{m['ncc']:.4f}"]
            )
        if real is not None:
            cells[(side, "real")] = real
            rm = _metrics(real, real)
            rows.append([side, "real", f"{rm['peak_mean']:.4f}", f"{rm['top_bot']:.4f}", "1.0000"])

    # 全灯ロフト重みスイープ(w2p0/w4p0 は 4×3 パネル外のため summary にのみ追記)
    real_full = _load_luminance(REAL["full"])
    sweep_cells: list[tuple[str, np.ndarray, dict[str, float]]] = []
    for label, job_dir in FULL_SWEEP:
        npy = job_dir / "luminance.npy"
        if not npy.exists():
            print(f"  警告: {npy} 未完了(sweep パネルからスキップ)")
            continue
        roi = extract_roi_xywh(_load_luminance(npy), SIM_ROI)
        m = _metrics(roi, real_full)
        sweep_cells.append((label, roi, m))
        if label.startswith(("w2.0", "w4.0")):
            rows.append(
                ["full", label, f"{m['peak_mean']:.4f}", f"{m['top_bot']:.4f}", f"{m['ncc']:.4f}"]
            )
    sweep_cells.append(("実機(w0.5)", real_full, _metrics(real_full, real_full)))

    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(["side", "version", "peak_mean", "top_bot", "ncc_vs_real"])
        writer.writerows(rows)
    print(f"summary 保存: {SUMMARY_CSV}")

    # 全セルを共通解像度に resize して行・列の枠を完全に揃える
    fig, axes = plt.subplots(
        len(SIDES),
        len(VERSIONS),
        figsize=(len(VERSIONS) * 4.2, len(SIDES) * 4.2),
        constrained_layout=True,
    )
    for r, side in enumerate(SIDES):
        for c, ver in enumerate(VERSIONS):
            ax = axes[r][c]
            ax.set_axis_off()
            if (side, ver) not in cells:
                ax.set_facecolor("0.15")
                if r == 0:
                    ax.set_title(VERSION_LABELS[ver], fontsize=11, pad=6)
                ax.text(
                    0.5,
                    0.5,
                    "実機全灯 npy 無し",
                    ha="center",
                    va="center",
                    transform=ax.transAxes,
                    fontsize=10,
                    color="0.4",
                )
                continue
            img = cells[(side, ver)]
            img = resize(img / img.mean(), CELL_SHAPE, anti_aliasing=True)
            ax.imshow(
                img,
                cmap="hot",
                vmin=float(np.percentile(img, VMIN_PERCENTILE)),
                vmax=float(np.percentile(img, VMAX_PERCENTILE)),
                aspect="equal",
            )
            m = _metrics(cells[(side, ver)], None)
            title = f"p/m={m['peak_mean']:.3f}  t/b={m['top_bot']:.3f}"
            if r == 0:
                title = f"{VERSION_LABELS[ver]}\n{title}"
            ax.set_title(title, fontsize=10, pad=6)
        # 行ラベル(左端セルの外側)
        axes[r][0].text(
            -0.06,
            0.5,
            side,
            ha="right",
            va="center",
            rotation=90,
            transform=axes[r][0].transAxes,
            fontsize=13,
            fontweight="bold",
        )

    fig.suptitle(
        "TEST_10 後側ロフト変形 w=1.4 — blend w0.05_a0.1, spp=1024, d=16, ROI 揃え\n"
        "各セル独立 1–99% 正規化(形の比較専用)・全セル共通解像度",
        fontsize=12,
    )
    PANEL_PATH.parent.mkdir(parents=True, exist_ok=True)
    fig.savefig(PANEL_PATH, dpi=150, bbox_inches="tight")
    plt.close(fig)
    print(f"パネル保存: {PANEL_PATH}")

    # 全灯ロフト重みスイープの 1×N パネル
    n = len(sweep_cells)
    fig2, axes2 = plt.subplots(1, n, figsize=(n * 4.2, 4.6), constrained_layout=True)
    for ax, (label, img, m) in zip(axes2, sweep_cells, strict=True):
        disp = resize(img / img.mean(), CELL_SHAPE, anti_aliasing=True)
        ax.imshow(
            disp,
            cmap="hot",
            vmin=float(np.percentile(disp, VMIN_PERCENTILE)),
            vmax=float(np.percentile(disp, VMAX_PERCENTILE)),
            aspect="equal",
        )
        ax.set_title(
            f"{label}\np/m={m['peak_mean']:.3f}  t/b={m['top_bot']:.3f}  NCC={m['ncc']:.3f}",
            fontsize=10,
            pad=6,
        )
        ax.set_axis_off()
    fig2.suptitle(
        "TEST_10 拡張 — 全灯ロフト重みスイープ(blend w0.05_a0.1, spp=1024, d=16, ROI 揃え)\n"
        "各セル独立 1–99% 正規化(形の比較専用)・全セル共通解像度",
        fontsize=12,
    )
    fig2.savefig(SWEEP_PANEL_PATH, dpi=150, bbox_inches="tight")
    plt.close(fig2)
    print(f"パネル保存: {SWEEP_PANEL_PATH}")

    # 相対輝度(各マップの max=1.0)版スイープパネル.全セル共通スケールで
    # 「全体がどれだけ 1 に近いか=均一性の絶対比較」を可能にする.
    # max はシミュ MC ノイズの外れ画素に引きずられないよう resize 後(共通解像度・
    # 実機相当に平滑化された状態)の最大値を用いる.
    rel_cells = []
    for label, img, _ in sweep_cells:
        if label.startswith("実機"):
            continue
        disp = resize(img.astype(np.float64), CELL_SHAPE, anti_aliasing=True)
        rel = disp / disp.max()
        rel_cells.append((label, rel))
    vmin = min(float(r.min()) for _, r in rel_cells)
    n_rel = len(rel_cells)
    fig3, axes3 = plt.subplots(1, n_rel, figsize=(n_rel * 4.4, 5.0), constrained_layout=True)
    for ax, (label, rel) in zip(axes3, rel_cells, strict=True):
        im = ax.imshow(rel, cmap="hot", vmin=vmin, vmax=1.0, aspect="equal")
        ax.set_title(
            f"{label}\nmean/max={rel.mean():.3f}  min/max={rel.min():.3f}",
            fontsize=10,
            pad=6,
        )
        ax.set_axis_off()
    fig3.colorbar(im, ax=axes3, shrink=0.85, label="相対輝度(max=1.0)")
    fig3.suptitle(
        "TEST_10 拡張 — 全灯ロフト重みスイープ・相対輝度(各セル max=1.0,全セル共通スケール)\n"
        "完全均一なら全面 1.0.mean/max・min/max が 1 に近いほど均一",
        fontsize=12,
    )
    fig3.savefig(SWEEP_REL_PANEL_PATH, dpi=150, bbox_inches="tight")
    plt.close(fig3)
    print(f"パネル保存: {SWEEP_REL_PANEL_PATH}")
    return 0


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