"""壁の反射特性(ω・α)と 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())