Newer
Older
SmTIAS-LightSim / scripts / build_shape_fit_comparison.py
"""壁の反射特性(ω・α)と spot 配光のスイープを,実機との「形」の一致で評価する(TEST_20).

sweep_wall_bsdf_spot.py の出力(output/wall_bsdf_spot/<条件>/luminance.npy)を読み,
compare_real_vs_sim_shapes.py と同じ処理で実機(全灯・下・左・右)と突き合わせる.

指標(構成ごと):
    rms_diff    : ROI 平均正規化マップの差分 RMS(小さいほど形が近い)
    peak_dx/dy  : σ=20 平滑化後のピーク位置(半幅比.+x 右・+y 上).実機との差 d_peak
    tilt / az   : 1 次平面フィットの振幅 % と方位
    bowl        : 2 次成分(お椀)の振幅 %

出力(output/wall_bsdf_spot/):
    shape_fit_a.csv / shape_fit_b.csv : 条件 × 構成の指標
    wall_bsdf_heatmap.png  : ρ・構成ごとに ω×α の差分 RMS とピークずれ誤差のヒートマップ
    spot_fwhm.png          : FWHM に対するピークずれ・傾き・差分 RMS(実機の値を水平線で)

実行方法:
    uv run python scripts/build_shape_fit_comparison.py
"""

from __future__ import annotations

import argparse
import csv
import json
import sys
from pathlib import Path

import matplotlib
import numpy as np
from scipy.ndimage import gaussian_filter

matplotlib.use("Agg")
import matplotlib.pyplot as plt  # noqa: E402

sys.path.insert(0, str(Path(__file__).resolve().parent))
sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))

from compare_real_vs_sim_shapes import (  # noqa: E402
    CONFIG_LABEL,
    CONFIGS,
    REAL,
    load_real,
    load_sim,
    normalized_shape,
    plane_fit,
)
from sweep_wall_bsdf_spot import ALPHAS, OMEGAS, OUTPUT_BASE, RHOS, SPOT_FWHMS  # noqa: E402

from sm_tias_lightsim.evaluation.zones import compute_normalized_distance_map  # noqa: E402

plt.rcParams["font.family"] = ["Yu Gothic", "MS Gothic", "DejaVu Sans"]
MAX_DISTANCE = 0.95
PEAK_SIGMA = 20.0


def describe(shape_map: np.ndarray) -> dict[str, float]:
    """正規化マップの形の記述子(ピーク位置・傾き・お椀)."""
    h, w = shape_map.shape
    heavy = gaussian_filter(shape_map, sigma=PEAK_SIGMA)
    iy, ix = np.unravel_index(int(np.argmax(heavy)), heavy.shape)
    tilt, az, bowl = plane_fit(shape_map)
    return {"peak_dx": (ix - w / 2) / (w / 2), "peak_dy": -(iy - h / 2) / (h / 2),
            "tilt": tilt, "tilt_az": az, "bowl": bowl}


def compare(sim_roi: np.ndarray, real_map: np.ndarray, real_desc: dict) -> dict[str, float]:
    """1 条件の形をシミュ vs 実機で比べる."""
    sim_map = normalized_shape(sim_roi)
    inner = compute_normalized_distance_map(*sim_map.shape) <= MAX_DISTANCE
    diff = sim_map - real_map
    d = describe(sim_map)
    d["rms_diff"] = float(np.sqrt(np.mean(diff[inner] ** 2)))
    d["d_peak"] = float(np.hypot(d["peak_dx"] - real_desc["peak_dx"],
                                 d["peak_dy"] - real_desc["peak_dy"]))
    d["d_tilt"] = float(d["tilt"] - real_desc["tilt"])
    return d


def load_all() -> tuple[list[dict], dict[str, dict]]:
    """全条件の row.json を読み,実機と比べた指標を付けて返す."""
    real_maps: dict[str, np.ndarray] = {}
    real_desc: dict[str, dict] = {}
    rows = []
    for path in sorted(OUTPUT_BASE.glob("*/row.json")):
        row = json.loads(path.read_text(encoding="utf-8"))
        npy = path.parent / "luminance.npy"
        if not npy.exists():
            continue
        sim_roi = load_sim(npy)
        cfg = row["config"]
        if cfg not in real_maps:
            real_roi = load_real(REAL[cfg], sim_roi.shape)
            real_maps[cfg] = normalized_shape(real_roi)
            real_desc[cfg] = describe(real_maps[cfg])
        row.update(compare(sim_roi, real_maps[cfg], real_desc[cfg]))
        rows.append(row)
    return rows, real_desc


def write_csv(rows: list[dict], phase: str) -> None:
    sel = [r for r in rows if r["phase"] == phase]
    if not sel:
        return
    path = OUTPUT_BASE / f"shape_fit_{phase}.csv"
    with path.open("w", newline="", encoding="utf-8") as fh:
        writer = csv.DictWriter(fh, fieldnames=list(sel[0].keys()))
        writer.writeheader()
        writer.writerows(sel)


def draw_heatmaps(rows: list[dict], real_desc: dict, output: Path) -> None:
    """ρ × 構成 ごとに ω×α の差分 RMS とピークずれ誤差を描く(ω=0 は α 列の先頭に置く)."""
    sel = [r for r in rows if r["phase"] == "a"]
    if not sel:
        return
    omegas = OMEGAS
    figsize = (4.4 * len(CONFIGS), 3.6 * 2 * len(RHOS))
    fig, axes = plt.subplots(2 * len(RHOS), len(CONFIGS), figsize=figsize)
    for ri, rho in enumerate(RHOS):
        for ci, cfg in enumerate(CONFIGS):
            for mi, (key, title, cmap) in enumerate((("rms_diff", "差分 RMS", "viridis_r"),
                                                     ("d_peak", "ピーク位置の誤差", "magma_r"))):
                grid = np.full((len(omegas), len(ALPHAS)), np.nan)
                for r in sel:
                    if r["rho"] != rho or r["config"] != cfg:
                        continue
                    oi = omegas.index(r["omega"])
                    if r["omega"] == 0.0:
                        grid[oi, :] = r[key]
                    else:
                        grid[oi, ALPHAS.index(r["alpha"])] = r[key]
                ax = axes[2 * ri + mi, ci]
                im = ax.imshow(grid, cmap=cmap, aspect="auto")
                ax.set_xticks(range(len(ALPHAS)))
                ax.set_xticklabels([f"{a:g}" for a in ALPHAS])
                ax.set_yticks(range(len(omegas)))
                ax.set_yticklabels([f"{o:g}" for o in omegas])
                ax.set_xlabel("α(粗さ)")
                ax.set_ylabel("ω(光沢の重み)")
                for oi in range(len(omegas)):
                    for ai in range(len(ALPHAS)):
                        if not np.isnan(grid[oi, ai]):
                            ax.text(ai, oi, f"{grid[oi, ai]:.3f}", ha="center", va="center",
                                    fontsize=7, color="white")
                ax.set_title(f"ρ={rho:g} {CONFIG_LABEL[cfg]}: {title}", fontsize=9)
                fig.colorbar(im, ax=ax, fraction=0.046)
    fig.suptitle("フェーズ A: 壁の反射特性(ω×α)と実機との形の一致"
                 "(小さいほど良い.ω=0 は純 diffuse)", fontsize=12)
    fig.tight_layout()
    fig.savefig(output, dpi=100)
    plt.close(fig)
    print(f"出力: {output}")


def draw_spot(rows: list[dict], real_desc: dict, output: Path) -> None:
    """フェーズ B: FWHM に対するピークずれ・傾き・差分 RMS."""
    sel = [r for r in rows if r["phase"] == "b"]
    if not sel:
        return
    fig, axes = plt.subplots(1, 4, figsize=(20, 5))
    colors = dict(zip(CONFIGS, plt.cm.tab10(np.arange(4)), strict=True))
    # 参照: 採用材質の面光源(フェーズ A の r0.85 w0.05 a0.10)
    ref = {r["config"]: r for r in rows
           if r["phase"] == "a" and r["rho"] == 0.85 and r["omega"] == 0.05 and r["alpha"] == 0.1}
    for ax, key, ylabel in ((axes[0], "peak_dx", "ピーク x(半幅比・右+)"),
                            (axes[1], "peak_dy", "ピーク y(半幅比・上+)"),
                            (axes[2], "tilt", "傾き [%]"), (axes[3], "rms_diff", "差分 RMS")):
        for cfg in CONFIGS:
            pts = sorted([r for r in sel if r["config"] == cfg], key=lambda r: r["fwhm"])
            ax.plot([r["fwhm"] for r in pts], [r[key] for r in pts], marker="o",
                    color=colors[cfg], label=f"{CONFIG_LABEL[cfg]}(spot)")
            if key != "rms_diff":
                ax.axhline(real_desc[cfg][key], color=colors[cfg], ls="--", alpha=0.7)
            if cfg in ref:
                ax.scatter([180.0], [ref[cfg][key]], marker="s", color=colors[cfg], zorder=3)
        ax.set_xlabel("spot の半値全幅 FWHM [deg](180 = ランバート面光源)")
        ax.set_ylabel(ylabel)
        ax.set_xticks(SPOT_FWHMS + [180.0])
        ax.grid(alpha=0.3)
        ax.legend(fontsize=7)
    axes[0].set_title("点線 = 実機の値,□ = 採用材質の面光源", fontsize=10)
    fig.suptitle("フェーズ B: LED の配光(spot FWHM)と単灯の明部のずれ・傾き", fontsize=12)
    fig.tight_layout()
    fig.savefig(output, dpi=100)
    plt.close(fig)
    print(f"出力: {output}")


def print_summary(rows: list[dict], real_desc: dict) -> None:
    """要点を表で出す."""
    print("\n=== 実機の形 ===")
    for cfg in CONFIGS:
        d = real_desc.get(cfg)
        if d:
            print(f"  {CONFIG_LABEL[cfg]:>8}: peak=({d['peak_dx']:+.2f},{d['peak_dy']:+.2f}) "
                  f"傾き {d['tilt']:.1f}%@{d['tilt_az']:.0f}° お椀 {d['bowl']:+.1f}%")
    for phase, title in (("a", "フェーズ A(壁): 4 構成の差分 RMS 平均が小さい順"),
                         ("b", "フェーズ B(spot): FWHM ごと")):
        sel = [r for r in rows if r["phase"] == phase]
        if not sel:
            continue
        print(f"\n=== {title} ===")
        keys = sorted({(r["rho"], r["omega"], r["alpha"], r.get("fwhm")) for r in sel})
        summary = []
        for k in keys:
            grp = {r["config"]: r for r in sel
                   if (r["rho"], r["omega"], r["alpha"], r.get("fwhm")) == k}
            if len(grp) < len(CONFIGS):
                continue
            summary.append((k, float(np.mean([grp[c]["rms_diff"] for c in CONFIGS])),
                            float(np.mean([grp[c]["d_peak"] for c in CONFIGS])), grp))
        summary.sort(key=lambda s: s[1] if phase == "a" else (s[0][3] or 0))
        print(f"{'ρ':>5} {'ω':>5} {'α':>5} {'FWHM':>5} {'RMS平均':>8} {'ピーク誤差平均':>12} | "
              f"{'左 peak x':>9} {'右 peak x':>9} {'全灯 傾き%':>10}")
        for k, rms, dpk, grp in summary[: (12 if phase == "a" else 99)]:
            fwhm = str(k[3] or "-")
            print(f"{k[0]:5.2f} {k[1]:5.2f} {k[2]:5.2f} {fwhm:>5} {rms:8.4f} {dpk:12.3f} | "
                  f"{grp['left']['peak_dx']:+9.2f} {grp['right']['peak_dx']:+9.2f} "
                  f"{grp['full']['tilt']:10.1f}")


def main() -> int:
    parser = argparse.ArgumentParser(description="壁反射特性・spot スイープの形の一致を評価")
    parser.parse_args()
    rows, real_desc = load_all()
    if not rows:
        print("エラー: 結果がありません")
        return 1
    for phase in ("a", "b"):
        write_csv(rows, phase)
    draw_heatmaps(rows, real_desc, OUTPUT_BASE / "wall_bsdf_heatmap.png")
    draw_spot(rows, real_desc, OUTPUT_BASE / "spot_fwhm.png")
    print_summary(rows, real_desc)
    return 0


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