diff --git a/scripts/build_shape_fit_comparison.py b/scripts/build_shape_fit_comparison.py new file mode 100644 index 0000000..96ab4e3 --- /dev/null +++ b/scripts/build_shape_fit_comparison.py @@ -0,0 +1,239 @@ +"""壁の反射特性（ω・α）と 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()) diff --git a/scripts/sweep_wall_bsdf_spot.py b/scripts/sweep_wall_bsdf_spot.py new file mode 100644 index 0000000..e983831 --- /dev/null +++ b/scripts/sweep_wall_bsdf_spot.py @@ -0,0 +1,224 @@ +"""壁の反射特性（ω・α）と LED 配光（spot）のスイープ: 単灯の明部のずれを実機に近づける（TEST_20）. + +実機と採用材質（ρ=0.85・ω=0.05・α=0.1）のシミュを全灯・単灯 3 灯で突き合わせると， +同心円方向（お椀）は合うが，単灯で明部が点灯基板の反対側へ寄る量がシミュは実機の半分以下 +だった（compare_real_vs_sim_shapes.py）．原因の候補は (a) 壁の反射が拡散しすぎ（光沢成分が +一次照射を反対側の壁に集中させる），(b) LED の配光がランバートより狭い，の 2 つ． +本スイープは基準シーン（scene_baseline.yaml）から 4 構成（全灯・下・左・右）を作り， + + フェーズ A: 壁 = blendbsdf(diffuse ρ, roughconductor α) の重み ω と α を振る + フェーズ B: 壁 = 採用材質のまま，LED を spot（半値全幅 FWHM を振る）にする + +をサーバでレンダする（評価は手元で build_shape_fit_comparison.py）． + +実行方法（サーバ）: + uv run sm-remote run --name wall --gpu 0 --threads 16 -- \ + python scripts/sweep_wall_bsdf_spot.py --phase all --workers 2 +""" + +from __future__ import annotations + +import argparse +import copy +import json +import multiprocessing +import sys +import time +from concurrent.futures import ProcessPoolExecutor, as_completed +from dataclasses import dataclass +from pathlib import Path + +import numpy as np +import yaml + +sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src")) + +from sm_tias_lightsim.render.pipeline import render_yaml # noqa: E402 +from sm_tias_lightsim.runtime.overrides import resolve_render_variant # noqa: E402 + +ROOT = Path(__file__).resolve().parents[1] +BASELINE_YAML = ROOT / "configs" / "scene_baseline.yaml" +TMP_DIR = ROOT / "configs" / "_tmp_wall_bsdf_spot" +OUTPUT_BASE = ROOT / "output" / "wall_bsdf_spot" + +CONFIGS = ("full", "bottom", "left", "right") +BOARD_PREFIX = {"bottom": "led_bottom", "left": "led_left", "right": "led_right"} + +# フェーズ A: 壁の反射特性 +RHOS = [0.85, 0.81] +OMEGAS = [0.0, 0.05, 0.15, 0.30, 0.50, 0.70] +ALPHAS = [0.02, 0.05, 0.1, 0.3, 0.6] +SPP_A = 2048 +# フェーズ B: LED 配光（半値全幅 FWHM）．採用材質で +SPOT_FWHMS = [30.0, 45.0, 60.0, 90.0, 120.0] +SPOT_MATERIAL = (0.85, 0.05, 0.1) +SPP_B = 4096 +# 面光源（2.8×3.5 mm の矩形・radiance）と同じ全光束になる spot の軸上強度を作る +AREA_MM2 = 2.8 * 3.5 + + +@dataclass(frozen=True) +class Condition: + """1 レンダ条件.""" + + phase: str + config: str + rho: float + omega: float + alpha: float + fwhm: float | None = None + + @property + def name(self) -> str: + mat = f"r{self.rho:.2f}_w{self.omega:.2f}_a{self.alpha:.2f}" + spot = f"_spot{self.fwhm:g}" if self.fwhm is not None else "" + return f"{self.phase}_{self.config}_{mat}{spot}" + + +def conditions(phase: str) -> list[Condition]: + """条件の一覧（フェーズ A → B の順）.""" + out: list[Condition] = [] + if phase in ("a", "all"): + for rho in RHOS: + for omega in OMEGAS: + for alpha in ALPHAS if omega > 0 else [ALPHAS[2]]: + for cfg in CONFIGS: + out.append(Condition("a", cfg, rho, omega, alpha)) + if phase in ("b", "all"): + rho, omega, alpha = SPOT_MATERIAL + for fwhm in SPOT_FWHMS: + for cfg in CONFIGS: + out.append(Condition("b", cfg, rho, omega, alpha, fwhm)) + return out + + +def material_dict(rho: float, omega: float, alpha: float) -> dict: + """内壁材質の YAML 断片（ω=0 なら純 diffuse）.""" + if omega <= 0.0: + return {"type": "diffuse", "reflectance": [rho, rho, rho]} + return { + "type": "blendbsdf", "weight": omega, + "components": [ + {"type": "diffuse", "reflectance": [rho, rho, rho]}, + {"type": "roughconductor", "alpha": alpha, "distribution": "ggx"}, + ], + } + + +def spot_solid_angle(cutoff_deg: float, beam_deg: float) -> float: + """Mitsuba spot（beam_width まで 1，cutoff まで smoothstep で 0）の実効立体角 [sr].""" + theta = np.linspace(0.0, np.radians(cutoff_deg), 4001) + width = max(np.radians(cutoff_deg - beam_deg), 1e-9) + t = np.clip((np.radians(cutoff_deg) - theta) / width, 0, 1) + falloff = np.where(theta <= np.radians(beam_deg), 1.0, t * t * (3 - 2 * t)) + return float(2 * np.pi * np.trapezoid(falloff * np.sin(theta), theta)) + + +def to_spot(light: dict, fwhm: float) -> dict: + """面光源 1 灯を同じ位置・向き・全光束の spot に置き換える（FWHM = cutoff + beam）.""" + cutoff = fwhm * 2.0 / 3.0 + beam = fwhm / 3.0 + radiance = light["radiance"] + radiance = radiance[0] if isinstance(radiance, list) else radiance + flux = np.pi * radiance * AREA_MM2 * 1e-6 # ランバート面光源の全光束 [W] + intensity = flux / spot_solid_angle(cutoff, beam) + return { + "id": light["id"], "type": "spot", "position": light["center"], + "direction": light["direction"], "cutoff_angle": round(cutoff, 4), + "beam_width": round(beam, 4), "intensity": float(intensity), + } + + +def build_yaml(cond: Condition, base: dict) -> Path: + """条件の YAML を書き出す（メッシュは configs/meshes を相対参照）.""" + cfg = copy.deepcopy(base) + for mesh in cfg["geometry"]["meshes"]: + if not Path(mesh["path"]).is_absolute(): + mesh["path"] = "../" + mesh["path"] + cfg["materials"]["printer_white"] = material_dict(cond.rho, cond.omega, cond.alpha) + lights = cfg["lights"] + if cond.config != "full": + lights = [lt for lt in lights if lt["id"].startswith(BOARD_PREFIX[cond.config])] + if cond.fwhm is not None: + lights = [to_spot(lt, cond.fwhm) for lt in lights] + cfg["lights"] = lights + cfg["render"]["spp"] = SPP_B if cond.phase == "b" else SPP_A + cfg["name"] = cond.name + spot_note = f" / spot FWHM={cond.fwhm:g}°" if cond.fwhm else "" + cfg["description"] = (f"TEST_20 {cond.phase.upper()}: {cond.config} / ρ={cond.rho} " + f"ω={cond.omega} α={cond.alpha}{spot_note}") + TMP_DIR.mkdir(parents=True, exist_ok=True) + dest = TMP_DIR / f"scene_{cond.name}.yaml" + dest.write_text(yaml.safe_dump(cfg, allow_unicode=True, sort_keys=False), encoding="utf-8") + return dest + + +def run_job(cond: Condition, scene_yaml: Path, variant: str) -> tuple[str, float, str]: + """1 条件をレンダする（既存は再利用）．戻り値は (名前, 秒, 状態).""" + out_dir = OUTPUT_BASE / cond.name + if (out_dir / "row.json").exists(): + return cond.name, 0.0, "reuse" + start = time.time() + try: + render_yaml(scene_yaml, out_dir, variant) + except Exception as exc: # noqa: BLE001 - 1 条件の失敗でスイープを止めない + return cond.name, time.time() - start, f"fail: {exc}" + row = {"name": cond.name, "phase": cond.phase, "config": cond.config, "rho": cond.rho, + "omega": cond.omega, "alpha": cond.alpha, "fwhm": cond.fwhm, + "spp": SPP_B if cond.phase == "b" else SPP_A, "variant": variant, + "elapsed_s": round(time.time() - start, 1)} + (out_dir / "row.json").write_text(json.dumps(row, ensure_ascii=False, indent=1), + encoding="utf-8") + return cond.name, time.time() - start, "ok" + + +def main() -> int: + """スイープを実行する.""" + parser = argparse.ArgumentParser(description="壁の反射特性と spot 配光のスイープ") + parser.add_argument("--phase", choices=("a", "b", "all"), default="all") + parser.add_argument("--workers", type=int, default=1) + parser.add_argument("--prepare-only", action="store_true") + parser.add_argument("--list", action="store_true") + args = parser.parse_args() + conds = conditions(args.phase) + if args.list: + for c in conds: + print(c.name) + print(f"{len(conds)} 条件") + return 0 + base = yaml.safe_load(BASELINE_YAML.read_text(encoding="utf-8")) + scenes = {c.name: build_yaml(c, base) for c in conds} + log_path = ROOT / "output" / "_sweep_wall_bsdf_spot.log" + log_path.parent.mkdir(parents=True, exist_ok=True) + log_file = log_path.open("a", encoding="utf-8") + + def log(message: str) -> None: + print(message, flush=True) + log_file.write(message + "\n") + log_file.flush() + + log(f"=== 壁反射特性 / spot スイープ: {len(conds)} 条件（phase={args.phase}）===") + if args.prepare_only: + return 0 + variant = resolve_render_variant(str(base["render"]["variant"])) + done = 0 + if args.workers > 1: + context = multiprocessing.get_context("spawn") + with ProcessPoolExecutor(max_workers=args.workers, mp_context=context) as pool: + futures = {pool.submit(run_job, c, scenes[c.name], variant): c for c in conds} + for fut in as_completed(futures): + name, sec, status = fut.result() + log(f" {name}: {status} ({sec:.0f}s)") + done += status == "ok" + else: + for c in conds: + name, sec, status = run_job(c, scenes[c.name], variant) + log(f" {name}: {status} ({sec:.0f}s)") + done += status == "ok" + log(f"=== 完了: 新規 {done} / {len(conds)} 条件 ===") + log_file.close() + return 0 + + +if __name__ == "__main__": + sys.exit(main())