diff --git a/scripts/analyze_fourfold_component.py b/scripts/analyze_fourfold_component.py index 59c72c1..a5ac3ea 100644 --- a/scripts/analyze_fourfold_component.py +++ b/scripts/analyze_fourfold_component.py @@ -69,9 +69,15 @@ return extract_roi_xywh(arr, SIM_ROI) -def ring_harmonics(roi: np.ndarray) -> list[dict[str, float]]: +def ring_harmonics( + roi: np.ndarray, harmonics: tuple[int, ...] = HARMONICS +) -> list[dict[str, float]]: """全方位カバーの各輪帯について方位フーリエ係数（ring mean 比）を返す. + Args: + roi: ROI の輝度 2 次元配列． + harmonics: 計算する次数 m の組（既定は HARMONICS．多灯配置の TEST_18 では m=3 も使う）． + Returns: 輪帯ごとの dict のリスト．キーは d_center / n_pixels と，各 m について am / bm / se_m @@ -112,7 +118,7 @@ "d_center": float((ring_edges[ring - 1] + ring_edges[ring]) / 2.0), "n_pixels": float(mask.sum()), } - for m in HARMONICS: + for m in harmonics: cos_m, sin_m = np.cos(m * sector_phi), np.sin(m * sector_phi) row[f"a{m}"] = float(2.0 / N_SECTORS * np.sum(means * cos_m) / ring_mean) row[f"b{m}"] = float(2.0 / N_SECTORS * np.sum(means * sin_m) / ring_mean) diff --git a/scripts/build_led_placement_comparison.py b/scripts/build_led_placement_comparison.py new file mode 100644 index 0000000..ad0c252 --- /dev/null +++ b/scripts/build_led_placement_comparison.py @@ -0,0 +1,478 @@ +"""照明配置スイープの解析: 位置・向き・個数への依存を球と立方体で比較する（TEST_18）. + +sweep_led_placement.py の出力（output/led_placement_sweep/<条件>/）を読み，条件ごとに + +1. **方位フーリエ成分**（輪帯ごとの逆分散重み集約）: 1 灯目が +X 方位にあるので + m=1 の符号付き係数 a1 が「LED 側がどれだけ明るい/暗いか」，N 灯の等方位配置では + m=N の係数 aN が「灯数の痕跡」を直接与える．直交成分 b_m はノイズ対照 +2. **構造 CoV**（MC ノイズを差し引いた CoV）と輪帯 min/max +3. **輝度重心**（metrics.csv） + +を求め，次の 3 つの図と集約表を出力する． + +- position_sensitivity.png : θ に対する感度曲線（形状 × サイズ）と，配置依存の + 大きさをまとめた指標（θ 全域での構造 CoV の幅・|a1| の最大）の内面積依存 +- direction_comparison.png : 向き（板向き / 壁法線向き / 真下向き）の比較 +- count_comparison.png : 灯数の比較（構造 CoV・|m1|・|mN|） +- position_maps.png : 代表ペアの輝度マップ（θ を横に並べる） + +実行方法: + uv run python scripts/build_led_placement_comparison.py +""" + +from __future__ import annotations + +import argparse +import csv +import json +import sys +from pathlib import Path + +import matplotlib +import numpy as np + +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 analyze_fourfold_component import aggregate, load_roi, ring_harmonics # noqa: E402 +from build_radial_comparison import radial_stats, ring_min_max_ratio # noqa: E402 +from build_radius_sweep_visual import decompose_cov, normalized_map # noqa: E402 +from sweep_led_placement import ( # noqa: E402 + AIMS, + COUNT_N, + COUNT_THETA, + DIRECTION_THETAS, + OUTPUT_BASE, + SIZE_PAIRS, + THETA_GRID, + collect_rows, + write_summary, +) + +plt.rcParams["font.family"] = ["Yu Gothic", "MS Gothic", "DejaVu Sans"] + +HARMONICS = (1, 2, 3, 4) +SHAPE_COLOR = {"sphere": "tab:blue", "cube": "tab:red"} +SHAPE_LABEL = {"sphere": "球", "cube": "立方体"} +SHAPE_STYLE = {"sphere": dict(marker="o", ls="-"), "cube": dict(marker="s", ls="--")} +AIM_LABEL = {"board": "板中心向き", "normal": "壁法線向き", "down": "真下向き"} +MAP_THETAS = [0.0, 30.0, 60.0, 90.0, 120.0] +MAP_PAIR_INDEX = 2 # 輝度マップに使う代表ペア（r=170 / L=234） + +ANALYSIS_COLUMNS = [ + "name", "group", "shape", "size_mm", "pair", "theta_deg", "aim", "n_led", "rep", + "cavity_area_mm2", "baffle_radius_mm", "baffle_ratio", "occlusion_blocked_frac", + "raw_mean_roi", "cov", "structure_cov", "noise_cov", "ring_min_max", + "centroid_x", "centroid_y", "centroid_offset", + "a1", "b1", "se1", "m1_abs", "a2", "b2", "se2", "a3", "b3", "se3", "a4", "b4", "se4", + "aN", "mN_abs", +] + + +def size_label(shape: str, size: float) -> str: + """条件ラベル用のサイズ表記.""" + return f"{'r' if shape == 'sphere' else 'L'}={size:g}" + + +def pair_index(shape: str, size: float) -> int: + """サイズペアの番号（0=最小）を返す（該当なしは -1）.""" + idx = 0 if shape == "sphere" else 1 + for i, pair in enumerate(SIZE_PAIRS): + if pair[idx] == size: + return i + return -1 + + +def analyze_condition(row: dict) -> dict | None: + """1 条件の輝度マップから非対称・均一性の指標を集める.""" + directory = OUTPUT_BASE / row["name"] + lum = directory / "luminance.npy" + if not lum.exists(): + return None + result = {k: row.get(k) for k in ("name", "group", "shape", "size_mm", "theta_deg", "aim", + "n_led", "rep", "cavity_area_mm2", "baffle_radius_mm", + "occlusion_blocked_frac", "raw_mean_roi", "cov")} + result["pair"] = pair_index(row["shape"], row["size_mm"]) + result["baffle_ratio"] = row["baffle_radius_mm"] / row["size_mm"] + rings = ring_harmonics(load_roi(lum), HARMONICS) + for m in HARMONICS: + result[f"a{m}"], result[f"se{m}"] = aggregate(rings, f"a{m}") + result[f"b{m}"], _ = aggregate(rings, f"b{m}") + result["m1_abs"] = float(np.hypot(result["a1"], result["b1"])) + n = int(row["n_led"]) + result["aN"] = result[f"a{n}"] + result["mN_abs"] = float(np.hypot(result[f"a{n}"], result[f"b{n}"])) + _, structure, noise = decompose_cov(lum) + centers, means, _ = radial_stats(lum, 20) + result["structure_cov"] = structure + result["noise_cov"] = noise + result["ring_min_max"] = ring_min_max_ratio(centers, means) + result["centroid_x"] = row.get("intensity_centroid_x_norm", float("nan")) + result["centroid_y"] = row.get("intensity_centroid_y_norm", float("nan")) + result["centroid_offset"] = row.get("intensity_centroid_offset_norm", float("nan")) + return result + + +def select(results: list[dict], **conds) -> list[dict]: + """条件で絞り込む（値が None のキーは無視）.""" + out = [] + for r in results: + if all(v is None or r.get(k) == v for k, v in conds.items()): + out.append(r) + return out + + +def size_colors(n: int): + """サイズペアごとの色.""" + return plt.cm.viridis(np.linspace(0.0, 0.85, n)) + + +# --------------------------------------------------------------------------- +# 図 1: 位置の感度曲線 +# --------------------------------------------------------------------------- + + +def robustness_summary(results: list[dict]) -> list[dict]: + """位置スイープについて，形状 × ペアごとに配置依存の大きさを集約する.""" + rows = [] + for shape in ("sphere", "cube"): + for i, pair in enumerate(SIZE_PAIRS): + size = pair[0 if shape == "sphere" else 1] + pos = sorted(select(results, group="position", shape=shape, size_mm=size), + key=lambda r: r["theta_deg"]) + if len(pos) < 3: + continue + struct = np.array([r["structure_cov"] for r in pos]) + a1 = np.array([r["a1"] for r in pos]) + on_axis = next((r for r in pos if r["theta_deg"] == 0.0), None) + rows.append({ + "shape": shape, "pair": i, "size_mm": size, + "area": pos[0]["cavity_area_mm2"], "n_theta": len(pos), + "struct_onaxis": on_axis["structure_cov"] if on_axis else float("nan"), + "struct_mean": float(struct.mean()), + "struct_max": float(struct.max()), + "struct_range": float(struct.max() - struct.min()), + "struct_std": float(struct.std(ddof=1)), + "a1_absmax": float(np.abs(a1).max()), + "a1_absmean": float(np.abs(a1).mean()), + "theta_worst": float(pos[int(np.argmax(struct))]["theta_deg"]), + }) + return rows + + +def draw_position(results: list[dict], summary: list[dict], output: Path) -> None: + """θ に対する感度曲線と，配置依存の集約指標の内面積依存を描く.""" + fig, axes = plt.subplots(2, 3, figsize=(20, 11)) + colors = size_colors(len(SIZE_PAIRS)) + for shape in ("sphere", "cube"): + for i, pair in enumerate(SIZE_PAIRS): + size = pair[0 if shape == "sphere" else 1] + pos = sorted(select(results, group="position", shape=shape, size_mm=size), + key=lambda r: r["theta_deg"]) + if not pos: + continue + theta = [r["theta_deg"] for r in pos] + label = f"{SHAPE_LABEL[shape]} {size_label(shape, size)}" + style = dict(SHAPE_STYLE[shape], color=colors[i], ms=4, lw=1.4, label=label) + axes[0, 0].errorbar(theta, [r["a1"] * 100 for r in pos], + yerr=[r["se1"] * 100 for r in pos], capsize=2, **style) + axes[0, 1].plot(theta, [r["structure_cov"] for r in pos], **style) + axes[0, 2].plot(theta, [r["centroid_x"] for r in pos], **style) + axes[0, 0].set_ylabel("a1（輪帯平均比 %・正=LED 側が明るい）") + axes[0, 0].set_title("m=1 非対称 vs LED 極角 θ", fontsize=11) + axes[0, 1].set_ylabel("構造 CoV") + axes[0, 1].set_title("構造 CoV vs θ（低いほど均一）", fontsize=11) + axes[0, 2].set_ylabel("輝度重心 X（正規化）") + axes[0, 2].set_title("輝度重心の X オフセット vs θ", fontsize=11) + for ax in axes[0]: + ax.axhline(0.0, color="gray", lw=0.8, ls="--") + ax.set_xlabel("LED 極角 θ [deg]（0=底の中心・90=側壁の中段・120=板の直下）") + ax.set_xticks(THETA_GRID[::2]) + ax.grid(alpha=0.3) + ax.legend(fontsize=7, ncol=2) + + for key, ax, title, ylabel in ( + ("struct_range", axes[1, 0], "配置依存の幅: 構造 CoV の max−min（θ 全域）", + "構造 CoV の幅"), + ("a1_absmax", axes[1, 1], "非対称の最大: max|a1|（θ 全域）", "max|a1| [%]"), + ("struct_max", axes[1, 2], "最悪配置の構造 CoV（点線は軸上 θ=0）", "構造 CoV"), + ): + for shape in ("sphere", "cube"): + rows = [r for r in summary if r["shape"] == shape] + scale = 100.0 if key == "a1_absmax" else 1.0 + ax.plot([r["area"] for r in rows], [r[key] * scale for r in rows], + color=SHAPE_COLOR[shape], label=f"{SHAPE_LABEL[shape]}", + **{k: v for k, v in SHAPE_STYLE[shape].items()}) + if key == "struct_max": + ax.plot([r["area"] for r in rows], [r["struct_onaxis"] for r in rows], + color=SHAPE_COLOR[shape], marker=SHAPE_STYLE[shape]["marker"], + ls=":", alpha=0.6, label=f"{SHAPE_LABEL[shape]}（軸上）") + ax.set_xscale("log") + ax.set_xlabel("キャビティ内面積 (mm²・対数)") + ax.set_ylabel(ylabel) + ax.set_title(title, fontsize=11) + ax.grid(alpha=0.3, which="both") + ax.legend(fontsize=8) + fig.suptitle( + "照明位置への依存（1 灯・板中心向き・バッフル付き）: 球 vs 同内面積の立方体", + fontsize=15, + ) + fig.tight_layout() + fig.savefig(output, dpi=100) + plt.close(fig) + print(f"出力: {output}") + + +def draw_position_maps(results: list[dict], output: Path) -> None: + """代表ペアの輝度マップを θ ごとに並べる（輪帯ピーク=1.0・共通スケール）.""" + fig, axes = plt.subplots(2, len(MAP_THETAS), figsize=(3.6 * len(MAP_THETAS), 7.4)) + maps: dict[tuple[str, float], np.ndarray] = {} + for shape in ("sphere", "cube"): + size = SIZE_PAIRS[MAP_PAIR_INDEX][0 if shape == "sphere" else 1] + for theta in MAP_THETAS: + rows = select(results, group="position", shape=shape, size_mm=size, theta_deg=theta) + if rows: + maps[(shape, theta)] = normalized_map(OUTPUT_BASE / rows[0]["name"] + / "luminance.npy") + if not maps: + plt.close(fig) + return + vmin = float(np.percentile(np.concatenate([m.ravel() for m in maps.values()]), 0.5)) + image = None + for row_idx, shape in enumerate(("sphere", "cube")): + size = SIZE_PAIRS[MAP_PAIR_INDEX][0 if shape == "sphere" else 1] + for col, theta in enumerate(MAP_THETAS): + ax = axes[row_idx, col] + ax.axis("off") + data = maps.get((shape, theta)) + if data is None: + continue + image = ax.imshow(data, cmap="hot", vmin=vmin, vmax=1.0) + rows = select(results, group="position", shape=shape, size_mm=size, theta_deg=theta) + ax.set_title(f"{SHAPE_LABEL[shape]} {size_label(shape, size)} θ={theta:g}°\n" + f"a1={rows[0]['a1'] * 100:+.2f}% 構造CoV={rows[0]['structure_cov']:.4f}", + fontsize=9) + if image is not None: + fig.colorbar(image, ax=axes.ravel().tolist(), fraction=0.02, label="輪帯ピーク=1.0") + fig.suptitle("LED 極角 θ による輝度マップの変化（LED は +X＝画像右側．輪帯ピーク正規化）", + fontsize=13) + fig.savefig(output, dpi=100) + plt.close(fig) + print(f"出力: {output}") + + +# --------------------------------------------------------------------------- +# 図 2: 向き +# --------------------------------------------------------------------------- + + +def draw_direction(results: list[dict], output: Path) -> None: + """向き（板中心 / 壁法線 / 真下）の比較を θ × サイズで描く.""" + fig, axes = plt.subplots(2, len(DIRECTION_THETAS), figsize=(6.2 * len(DIRECTION_THETAS), 10)) + colors = size_colors(len(SIZE_PAIRS)) + width = 0.09 + for col, theta in enumerate(DIRECTION_THETAS): + for metric_row, (key, ylabel, scale) in enumerate(( + ("structure_cov", "構造 CoV", 1.0), ("m1_abs", "|m1|（輪帯平均比 %）", 100.0), + )): + ax = axes[metric_row, col] + for a_idx, aim in enumerate(AIMS): + for i, pair in enumerate(SIZE_PAIRS): + for s_idx, shape in enumerate(("sphere", "cube")): + size = pair[s_idx] + rows = select(results, shape=shape, size_mm=size, theta_deg=theta, + aim=aim, n_led=1, rep=1) + if not rows: + continue + x = a_idx + (i * 2 + s_idx - 3.5) * width + hatch = "" if shape == "sphere" else "//" + ax.bar(x, rows[0][key] * scale, width=width, color=colors[i], + edgecolor=SHAPE_COLOR[shape], hatch=hatch, lw=1.0) + if rows[0]["baffle_ratio"] > 0.3: + ax.text(x, rows[0][key] * scale, "†", ha="center", va="bottom", + fontsize=8) + ax.set_xticks(range(len(AIMS))) + ax.set_xticklabels([AIM_LABEL[a] for a in AIMS]) + ax.set_ylabel(ylabel) + ax.set_title(f"θ={theta:g}°: {ylabel}", fontsize=11) + ax.grid(alpha=0.3, axis="y") + handles = [plt.Rectangle((0, 0), 1, 1, color=colors[i], + label=f"ペア{i + 1}: r={p[0]:g} / L={p[1]:g}") + for i, p in enumerate(SIZE_PAIRS)] + handles += [ + plt.Rectangle((0, 0), 1, 1, fc="white", ec=SHAPE_COLOR["sphere"], label="球（枠青）"), + plt.Rectangle((0, 0), 1, 1, fc="white", ec=SHAPE_COLOR["cube"], hatch="//", + label="立方体（枠赤・斜線）"), + ] + axes[0, 0].legend(handles=handles, fontsize=8, loc="upper left") + fig.suptitle("LED の向きの影響（1 灯）．† = バッフル半径がキャビティ寸法の 30% 超" + "（板向き・法線向きはバッフル付き，真下向きは板が背後でバッフル無し）", + fontsize=13) + fig.tight_layout() + fig.savefig(output, dpi=100) + plt.close(fig) + print(f"出力: {output}") + + +# --------------------------------------------------------------------------- +# 図 3: 灯数 +# --------------------------------------------------------------------------- + + +def draw_count(results: list[dict], output: Path) -> None: + """灯数 N の効果（θ=90°）と N=4 の θ 依存を描く.""" + fig, axes = plt.subplots(1, 3, figsize=(20, 6)) + colors = size_colors(len(SIZE_PAIRS)) + ns = [1] + COUNT_N + for shape in ("sphere", "cube"): + for i, pair in enumerate(SIZE_PAIRS): + size = pair[0 if shape == "sphere" else 1] + rows = [select(results, shape=shape, size_mm=size, theta_deg=COUNT_THETA, + aim="board", n_led=n, rep=1) for n in ns] + xs = [n for n, r in zip(ns, rows, strict=True) if r] + picked = [r[0] for r in rows if r] + if not picked: + continue + label = f"{SHAPE_LABEL[shape]} {size_label(shape, size)}" + style = dict(SHAPE_STYLE[shape], color=colors[i], ms=5, lw=1.4, label=label) + axes[0].plot(xs, [r["structure_cov"] for r in picked], **style) + axes[1].plot(xs, [r["m1_abs"] * 100 for r in picked], **style) + axes[1].plot(xs, [r["mN_abs"] * 100 for r in picked], marker="x", ls=":", + color=colors[i], lw=1.0, alpha=0.8) + n4 = sorted(select(results, shape=shape, size_mm=size, aim="board", n_led=4, rep=1), + key=lambda r: r["theta_deg"]) + if n4: + axes[2].plot([r["theta_deg"] for r in n4], [r["structure_cov"] for r in n4], + **style) + axes[0].set_xlabel("灯数 N（θ=90° の輪に等方位・総光束一定）") + axes[0].set_ylabel("構造 CoV") + axes[0].set_title("灯数 vs 構造 CoV", fontsize=11) + axes[1].set_xlabel("灯数 N") + axes[1].set_ylabel("方位成分（輪帯平均比 %）") + axes[1].set_title("|m1|（実線）と |m=N|（点線 ×）vs 灯数", fontsize=11) + axes[2].set_xlabel("LED 極角 θ [deg]") + axes[2].set_ylabel("構造 CoV") + axes[2].set_title("N=4 の θ 依存（θ=60/90/120）", fontsize=11) + for ax in axes: + ax.grid(alpha=0.3) + ax.legend(fontsize=7, ncol=2) + axes[0].set_xticks(ns) + axes[1].set_xticks(ns) + fig.suptitle("灯数の効果: 多灯化で配置の痕跡はどれだけ消えるか（板中心向き・各灯にバッフル）", + fontsize=14) + fig.tight_layout() + fig.savefig(output, dpi=100) + plt.close(fig) + print(f"出力: {output}") + + +# --------------------------------------------------------------------------- +# 表 +# --------------------------------------------------------------------------- + + +def print_tables(results: list[dict], summary: list[dict]) -> None: + """集約表を標準出力に書く.""" + print("\n=== 位置スイープの集約（θ 全域）===") + print(f"{'形状':>6} {'サイズ':>7} {'内面積':>9} {'θ数':>3} {'構造CoV軸上':>11} {'平均':>8} " + f"{'最大':>8} {'幅':>8} {'max|a1|%':>9} {'最悪θ':>6}") + for r in summary: + print(f"{SHAPE_LABEL[r['shape']]:>6} {r['size_mm']:7g} {r['area']:9.0f} {r['n_theta']:3d} " + f"{r['struct_onaxis']:11.4f} {r['struct_mean']:8.4f} {r['struct_max']:8.4f} " + f"{r['struct_range']:8.4f} {r['a1_absmax'] * 100:9.2f} {r['theta_worst']:6.0f}") + + print("\n=== 位置スイープの全条件 ===") + print(f"{'条件':>30} {'a1 [%]':>14} {'b1 [%]':>8} {'構造CoV':>8} {'ノイズ':>7} " + f"{'輪帯min/max':>11} {'重心X':>8} {'生mean':>7} {'遮蔽':>6}") + for r in sorted(select(results, group="position"), + key=lambda r: (r["shape"], r["size_mm"], r["theta_deg"])): + occl = r["occlusion_blocked_frac"] + occl = float(occl) if occl not in ("", None) else float("nan") + print(f"{r['name']:>30} {r['a1'] * 100:7.3f}±{r['se1'] * 100:5.3f} {r['b1'] * 100:8.3f} " + f"{r['structure_cov']:8.4f} {r['noise_cov']:7.4f} {r['ring_min_max']:11.4f} " + f"{r['centroid_x']:8.4f} {r['raw_mean_roi']:7.3f} {occl:6.3f}") + + print("\n=== 向き（1 灯）===") + print(f"{'条件':>30} {'バッフルR/寸法':>13} {'a1 [%]':>8} {'|m1| [%]':>9} {'構造CoV':>8} " + f"{'生mean':>7}") + for r in sorted(select(results, group="direction") + select(results, group="position", + n_led=1), + key=lambda r: (r["shape"], r["size_mm"], r["theta_deg"], r["aim"])): + if r["theta_deg"] not in DIRECTION_THETAS: + continue + print(f"{r['name']:>30} {r['baffle_ratio']:13.3f} {r['a1'] * 100:8.3f} " + f"{r['m1_abs'] * 100:9.3f} {r['structure_cov']:8.4f} {r['raw_mean_roi']:7.3f}") + + print("\n=== 灯数 ===") + print(f"{'条件':>30} {'|m1| [%]':>9} {'aN [%]':>8} {'|mN| [%]':>9} {'構造CoV':>8} " + f"{'生mean':>7}") + for r in sorted(select(results, aim="board", rep=1), + key=lambda r: (r["shape"], r["size_mm"], r["theta_deg"], r["n_led"])): + if r["n_led"] == 1 and r["theta_deg"] not in (COUNT_THETA, 60.0, 120.0): + continue + print(f"{r['name']:>30} {r['m1_abs'] * 100:9.3f} {r['aN'] * 100:8.3f} " + f"{r['mN_abs'] * 100:9.3f} {r['structure_cov']:8.4f} {r['raw_mean_roi']:7.3f}") + + print("\n=== 再現性（同一条件の再レンダ）===") + for shape in ("sphere", "cube"): + reps = [r for r in results if r["shape"] == shape and r["theta_deg"] == COUNT_THETA + and r["aim"] == "board" and r["n_led"] == 1 + and pair_index(shape, r["size_mm"]) == 2] + if len(reps) < 2: + continue + struct = np.array([r["structure_cov"] for r in reps]) + a1 = np.array([r["a1"] for r in reps]) + print(f" {SHAPE_LABEL[shape]} {size_label(shape, reps[0]['size_mm'])} n={len(reps)}: " + f"構造CoV {struct.mean():.4f} ± {struct.std(ddof=1):.4f}, " + f"a1 {a1.mean() * 100:.3f} ± {a1.std(ddof=1) * 100:.3f} %") + + +def write_analysis_csv(results: list[dict], summary: list[dict]) -> None: + """条件別の解析値と位置スイープの集約を CSV に書く.""" + path = OUTPUT_BASE / "analysis.csv" + with path.open("w", newline="", encoding="utf-8") as fh: + writer = csv.DictWriter(fh, fieldnames=ANALYSIS_COLUMNS) + writer.writeheader() + for r in results: + writer.writerow({k: r.get(k, "") for k in ANALYSIS_COLUMNS}) + path = OUTPUT_BASE / "position_summary.csv" + if summary: + with path.open("w", newline="", encoding="utf-8") as fh: + writer = csv.DictWriter(fh, fieldnames=list(summary[0].keys())) + writer.writeheader() + writer.writerows(summary) + (OUTPUT_BASE / "analysis.json").write_text( + json.dumps({"results": results, "position_summary": summary}, ensure_ascii=False, + indent=1), + encoding="utf-8", + ) + + +def main() -> int: + """照明配置スイープの解析図と集約表を生成する.""" + parser = argparse.ArgumentParser(description="照明配置スイープの解析図を生成") + parser.add_argument("--no-figures", action="store_true", help="表と CSV だけ出す") + args = parser.parse_args() + + write_summary() + rows = collect_rows() + results = [res for res in (analyze_condition(r) for r in rows) if res is not None] + if not results: + print("エラー: スイープの結果がありません（先にレンダを実行してください）") + return 1 + summary = robustness_summary(results) + write_analysis_csv(results, summary) + if not args.no_figures: + draw_position(results, summary, OUTPUT_BASE / "position_sensitivity.png") + draw_position_maps(results, OUTPUT_BASE / "position_maps.png") + draw_direction(results, OUTPUT_BASE / "direction_comparison.png") + draw_count(results, OUTPUT_BASE / "count_comparison.png") + print_tables(results, summary) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/scripts/profile_sweep_step.py b/scripts/profile_sweep_step.py new file mode 100644 index 0000000..18a69a7 --- /dev/null +++ b/scripts/profile_sweep_step.py @@ -0,0 +1,142 @@ +"""スイープ 1 条件の所要時間の内訳を測る（高速化の当たりを付けるための計測）. + +sm-render / sm-evaluate をサブプロセスで呼ぶ現行方式と，同じ処理をプロセス内で +直接呼ぶ方式の各段階（起動・import・PLY 書き出し・シーン読込・レンダ・評価・図・遮蔽検証） +を秒で出す． + +実行方法: + uv run python scripts/profile_sweep_step.py \ + configs/_tmp_led_placement/scene_sphere_r170_th090_board_n1.yaml +""" + +from __future__ import annotations + +import argparse +import os +import shutil +import subprocess +import sys +import tempfile +import time +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parent)) +sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src")) + +ROOT = Path(__file__).resolve().parents[1] +EVAL_YAML = ROOT / "configs" / "eval_baseline_dng_target_shape.yaml" + + +class Timer: + """段階ごとの経過秒を記録する.""" + + def __init__(self) -> None: + self.rows: list[tuple[str, float]] = [] + self._t = time.perf_counter() + + def lap(self, label: str) -> None: + now = time.perf_counter() + self.rows.append((label, now - self._t)) + self._t = now + + def report(self, title: str) -> None: + print(f"\n=== {title} ===") + for label, sec in self.rows: + print(f" {label:<40} {sec:8.2f} s") + print(f" {'合計':<40} {sum(s for _, s in self.rows):8.2f} s") + + +def profile_subprocess(scene: Path, out_dir: Path) -> None: + """現行方式: uv run sm-render → uv run sm-evaluate をサブプロセスで.""" + t = Timer() + subprocess.run(["uv", "run", "python", "-c", "pass"], cwd=ROOT, check=True, + capture_output=True) + t.lap("uv run python -c pass（起動のみ）") + subprocess.run(["uv", "run", "python", "-c", "import mitsuba"], cwd=ROOT, check=True, + capture_output=True) + t.lap("uv run python -c 'import mitsuba'") + subprocess.run(["uv", "run", "sm-render", "--scene", str(scene), "--output", str(out_dir)], + cwd=ROOT, check=True, capture_output=True) + t.lap("uv run sm-render（全体）") + subprocess.run(["uv", "run", "sm-evaluate", "--image", str(out_dir / "luminance.npy"), + "--config", str(EVAL_YAML), "--output", str(out_dir)], + cwd=ROOT, check=True, capture_output=True) + t.lap("uv run sm-evaluate（全体）") + t.report("現行方式（サブプロセス）") + + +def profile_inprocess(scene: Path, out_dir: Path, variant: str) -> None: + """プロセス内で各段階を直接呼ぶ.""" + t = Timer() + import numpy as np + + from sm_tias_lightsim.runtime.llvm_init import ensure_drjit_libllvm_path + + if variant.startswith("llvm_"): + ensure_drjit_libllvm_path() + import mitsuba as mi + + t.lap("import mitsuba") + mi.set_variant(variant) + t.lap(f"set_variant({variant})") + + from sm_tias_lightsim.scene.builder import build_scene + from sm_tias_lightsim.scene.schema import load_scene_yaml + + config = load_scene_yaml(scene) + config = config.model_copy( + update={"render": config.render.model_copy(update={"variant": variant})} + ) + t.lap("load_scene_yaml") + scene_dict, tempdir = build_scene(config, project_root=scene.parent) + t.lap("build_scene（STL 読込 + PLY 書き出し）") + try: + mscene = mi.load_dict(scene_dict) + t.lap("mi.load_dict（PLY 読込 + 加速構造）") + for i in range(2): + img = mi.render(mscene, spp=config.render.spp) + arr = np.array(img) + t.lap(f"mi.render spp={config.render.spp}（{i + 1} 回目）") + finally: + shutil.rmtree(tempdir, ignore_errors=True) + from sm_tias_lightsim.render.output import save_luminance, save_png + + save_png(arr, out_dir / "render.png") + save_luminance(arr, out_dir / "luminance.npy") + t.lap("save_png + save_luminance") + + from sm_tias_lightsim.evaluation.pipeline import evaluate_image + + image = out_dir / "luminance.npy" + evaluate_image(image, EVAL_YAML, out_dir, save_figures=True) + t.lap("evaluate_image（図あり）") + evaluate_image(image, EVAL_YAML, out_dir, save_figures=False) + t.lap("evaluate_image（図なし）") + + from verify_baffle_occlusion import check_occlusion + + result = check_occlusion(scene, board_grid=41, led_count=32) + t.lap(f"check_occlusion（variant={mi.variant()}・{result['blocked']:.0f}/{result['total']:.0f}）") + result = check_occlusion(scene, board_grid=41, led_count=32) + t.lap("check_occlusion（2 回目）") + t.report(f"プロセス内方式（variant={variant}）") + + +def main() -> int: + parser = argparse.ArgumentParser(description="スイープ 1 条件の所要時間の内訳") + parser.add_argument("scene", type=Path) + parser.add_argument("--variant", default=os.environ.get("SM_RENDER_VARIANT", "llvm_ad_rgb")) + parser.add_argument("--skip-subprocess", action="store_true") + args = parser.parse_args() + out_dir = Path(tempfile.mkdtemp(prefix="sm_profile_")) + try: + if not args.skip_subprocess: + profile_subprocess(args.scene.resolve(), out_dir) + profile_inprocess(args.scene.resolve(), out_dir, args.variant) + finally: + shutil.rmtree(out_dir, ignore_errors=True) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/scripts/show_scene_cutaway.py b/scripts/show_scene_cutaway.py index b16d6b8..41f6d4e 100644 --- a/scripts/show_scene_cutaway.py +++ b/scripts/show_scene_cutaway.py @@ -79,7 +79,7 @@ shown = mesh.clip(normal=(0, 1, 0), origin=(0, 0, 0)) if clipped else mesh if clipped: color, opacity = HOUSING_COLOR, 1.0 - elif entry["id"] == "baffle": + elif entry["id"].startswith("baffle"): color, opacity = BAFFLE_COLOR, 1.0 else: color, opacity = BOARD_COLOR, 0.55 diff --git a/scripts/sweep_led_placement.py b/scripts/sweep_led_placement.py new file mode 100644 index 0000000..e0016b5 --- /dev/null +++ b/scripts/sweep_led_placement.py @@ -0,0 +1,776 @@ +"""照明配置スイープ: LED の位置・向き・個数への依存を球と立方体で比較する（TEST_18）. + +TEST_16 で「軸上配置なら同内面積の立方体が球より均一」，TEST_17 で「側壁 1 点への +オフ軸配置では球のほうが頑健」と分かった．本スイープは「積分球が使われる理由は +照明配置への依存が小さいからではないか」という仮説を，配置を系統的に振って検証する． + +条件の軸（同内面積の球・立方体ペアごと）: + A. 位置 : LED をキャビティ中心から見た極角 θ（0°=底の中心＝軸上，90°=中心と + 同じ高さの側壁，120°=板のすぐ下）で壁面上を動かす．板中心向き＋バッフル + B. 向き : θ=60/90/120° で「板中心向き」「壁法線向き（キャビティ中心向き）」 + 「真下向き（板が発光面の背後＝バッフル不要）」を比べる + C. 個数 : θ の輪に N 灯を等方位に置き総光束を保存する（radiance / N） + D. 再現性 : 同一条件の再レンダでノイズフロアを測る + +θ の定義（両形状共通）: + キャビティ中心 c（球: 球中心，立方体: 立方体中心）から方向 + u = (sinθ cosφ, sinθ sinφ, -cosθ) に伸ばした半直線が内壁と交わる点に LED を置く． + 球ではこれが「LED から見た板の立体角」を一意に決める自然なパラメータで， + 立方体には同じ角度をそのまま適用する（θ<45° は底面，45°<θ<135° は側壁）． + 立方体の稜（θ=45°）は LED を置けないので θ グリッドから外す． + +バッフルは TEST_17 と同じ「LED から BAFFLE_GAP 先に置く傾いた円板」を，発光面の前方に +ある板の全点を隠す最小半径 × 安全率で数値的に決める．TEST_17 では外殻の貫通を +エラーにしたが，板の近く（θ が大きい）では円板が壁を大きく貫かざるを得ないので +本スイープでは許容し，食い込み量を記録するだけにする（外側に出た部分はキャビティ +内から見えないため光輸送には影響しない）． + +固定条件は TEST_15/16/17 と同一（材質・レンダ・評価・評価板・LED 面積・総光束）． + +実行方法: + uv run python scripts/sweep_led_placement.py --prepare-only # 幾何 + YAML のみ + uv run python scripts/sweep_led_placement.py --shapes sphere # 球だけレンダ + uv run python scripts/sweep_led_placement.py --groups position count # グループ指定 + uv run python scripts/sweep_led_placement.py --workers 3 # 3 条件同時に実行 +""" + +from __future__ import annotations + +import argparse +import copy +import csv +import io +import json +import multiprocessing +import sys +import time +from concurrent.futures import ProcessPoolExecutor, as_completed +from contextlib import redirect_stdout +from dataclasses import dataclass +from pathlib import Path + +import numpy as np +import trimesh +import yaml + +sys.path.insert(0, str(Path(__file__).resolve().parent)) +sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src")) + +from generate_cube_housing import build_cube_housing # noqa: E402 +from sweep_cube_size import cavity_geometry as cube_cavity_geometry # noqa: E402 +from sweep_offaxis_led import STANDOFF_MARGIN # noqa: E402 +from sweep_sphere_radius import ( # noqa: E402 + BAFFLE_GAP, + BAFFLE_MARGIN, + BAFFLE_THICKNESS, + BOARD_HALF, + BOARD_SURFACE_Z, + EVAL_CONFIG_YAML, + LED_RADIUS, + RIM_INSET, + RIM_Z, + TEMPLATE_YAML, + THICKNESS, + raw_stats, + read_metrics_csv, +) +from sweep_sphere_radius import ( + build_mesh as build_sphere_mesh, +) +from sweep_sphere_radius import ( + cavity_geometry as sphere_cavity_geometry, +) +from sweep_sphere_radius import ( + center_z as sphere_center_z, +) +from sweep_sphere_radius import ( + cut_z as sphere_cut_z, +) + +from sm_tias_lightsim.evaluation.pipeline import evaluate_image # noqa: E402 +from sm_tias_lightsim.render.pipeline import render_yaml # noqa: E402 +from sm_tias_lightsim.runtime.overrides import resolve_render_variant # noqa: E402 +from sm_tias_lightsim.scene.paths import relative_mesh_path # noqa: E402 + +ROOT = Path(__file__).resolve().parents[1] + +TMP_DIR = ROOT / "configs" / "_tmp_led_placement" +MESH_DIR = TMP_DIR / "meshes" +OUTPUT_BASE = ROOT / "output" / "led_placement_sweep" + +SPHERE_CAVITY = "similar_radial" # TEST_15 最良・TEST_17 と同じ + +# 同内面積ペア（球の内半径, 立方体の内寸）．立方体は 6L^2+594 = 球の内面積 から丸めた +SIZE_PAIRS: list[tuple[float, float]] = [ + (100.0, 141.0), + (135.0, 187.0), + (170.0, 234.0), + (250.0, 340.0), +] +REPEAT_PAIR_INDEX = 2 # 再現性チェックに使うペア（r=170 / L=234） + +# A. 位置: 極角グリッド（立方体の稜 45° を避ける．120° は球の loft 接続（約 121°）直下） +THETA_GRID = [0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 110.0, 120.0] +# B. 向き +DIRECTION_THETAS = [60.0, 90.0, 120.0] +AIMS = ("board", "normal", "down") +# C. 個数 +COUNT_THETA = 90.0 +COUNT_N = [2, 3, 4] +COUNT_N4_THETAS = [60.0, 120.0] +# D. 再現性（同一条件の追加レンダ回数） +REPEATS = 3 + +GROUPS = ("position", "direction", "count", "repeat") + +# 評価板の中心（キャビティ側の面）＝ 板向き LED の照準点 +BOARD_CENTER = np.array([0.0, 0.0, BOARD_SURFACE_Z]) +# LED 円板を立方体の稜に近づけすぎない余裕（mm） +EDGE_MARGIN = 3.0 +# バッフル半径の上限（キャビティ寸法比）．超える条件は物理的に無理があるので実行しない +BAFFLE_RADIUS_CAP_RATIO = 0.6 + +SUMMARY_COLUMNS = [ + "name", "group", "shape", "size_mm", "theta_deg", "aim", "n_led", "rep", + "cavity_area_mm2", "led_x_mm", "led_y_mm", "led_z_mm", "led_standoff_mm", "tilt_deg", + "baffled", "baffle_radius_mm", "baffle_embed_mm", "board_front_frac", + "occlusion_blocked_frac", "occlusion_total", + "elapsed_s", "render_variant", "raw_mean_roi", "raw_mean_full", "raw_rel_std_roi", + "cov", "center_periphery_ratio", "gradient_magnitude", "radial_min_max_ratio", + "angular_ratio", "pca_eccentricity", "intensity_centroid_offset_norm", + "intensity_centroid_x_norm", "intensity_centroid_y_norm", "max_min_ratio", +] + + +@dataclass(frozen=True) +class Condition: + """1 レンダ条件（配置の指定）.""" + + group: str + shape: str + size: float + theta_deg: float + aim: str + n_led: int + rep: int = 1 + + @property + def base_name(self) -> str: + """幾何（メッシュ・バッフル）を共有する名前（再レンダ条件は seed だけが違う）.""" + prefix = "r" if self.shape == "sphere" else "L" + return (f"{self.shape}_{prefix}{self.size:g}_th{self.theta_deg:03.0f}" + f"_{self.aim}_n{self.n_led}") + + @property + def name(self) -> str: + """出力ディレクトリ名.""" + return self.base_name + (f"_rep{self.rep}" if self.rep > 1 else "") + + @property + def azimuths_deg(self) -> list[float]: + """各 LED の方位角（1 灯目は +X 側，以降等間隔）.""" + return [360.0 * k / self.n_led for k in range(self.n_led)] + + +def conditions() -> list[Condition]: + """全条件を「重要な順」（位置 → 個数 → 向き → 再現性）に並べて返す. + + 途中で止まっても主実験（位置）が先に揃うようにする．同じ条件が複数グループから + 生まれる場合（θ=90 板向き 1 灯は位置・向き・個数の共通アンカー）は最初の 1 つだけ残す． + """ + result: list[Condition] = [] + seen: set[str] = set() + + def add(cond: Condition) -> None: + if cond.name not in seen: + seen.add(cond.name) + result.append(cond) + + for shape_index, shape in enumerate(("sphere", "cube")): + for pair in SIZE_PAIRS: + size = pair[shape_index] + for theta in THETA_GRID: + add(Condition("position", shape, size, theta, "board", 1)) + for shape_index, shape in enumerate(("sphere", "cube")): + for pair in SIZE_PAIRS: + size = pair[shape_index] + for n in COUNT_N: + add(Condition("count", shape, size, COUNT_THETA, "board", n)) + for theta in COUNT_N4_THETAS: + add(Condition("count", shape, size, theta, "board", 4)) + for shape_index, shape in enumerate(("sphere", "cube")): + for pair in SIZE_PAIRS: + size = pair[shape_index] + for theta in DIRECTION_THETAS: + for aim in AIMS: + add(Condition("direction", shape, size, theta, aim, 1)) + for shape_index, shape in enumerate(("sphere", "cube")): + size = SIZE_PAIRS[REPEAT_PAIR_INDEX][shape_index] + for rep in range(2, REPEATS + 2): + add(Condition("repeat", shape, size, COUNT_THETA, "board", 1, rep)) + return result + + +# --------------------------------------------------------------------------- +# 幾何 +# --------------------------------------------------------------------------- + + +def cavity_center(shape: str, size: float) -> np.ndarray: + """θ の基準になるキャビティ中心（mm）を返す.""" + if shape == "cube": + return np.array([0.0, 0.0, RIM_Z - size / 2.0]) + return np.array([0.0, 0.0, sphere_center_z(size, SPHERE_CAVITY)]) + + +def cavity_area(shape: str, size: float) -> float: + """比較の横軸に使う内面積（mm²）を返す（各スイープと同じ定義）.""" + if shape == "cube": + return cube_cavity_geometry(size)["cavity_area_mm2"] + return sphere_cavity_geometry(size, SPHERE_CAVITY)["cavity_area_mm2"] + + +def polar_direction(theta_deg: float, azimuth_deg: float) -> np.ndarray: + """キャビティ中心から見た単位方向 u（θ=0 が -Z＝底）を返す.""" + theta, phi = np.radians(theta_deg), np.radians(azimuth_deg) + return np.array([np.sin(theta) * np.cos(phi), np.sin(theta) * np.sin(phi), -np.cos(theta)]) + + +def wall_anchor( + shape: str, size: float, theta_deg: float, azimuth_deg: float +) -> tuple[np.ndarray, np.ndarray]: + """LED を貼り付ける内壁上の点と，そこでの内向き法線を返す（mm）. + + Raises: + ValueError: 立方体の稜・球の loft 接続に近すぎて LED 円板を置けない場合． + """ + center = cavity_center(shape, size) + u = polar_direction(theta_deg, azimuth_deg) + if shape == "sphere": + point = center + size * u + # 球の壁は loft 接続平面（cut_z）より下だけ．接続の上では法線の定義が変わる + if point[2] > sphere_cut_z(size, SPHERE_CAVITY): + raise ValueError(f"θ={theta_deg:g}: 球の loft 接続より上に出ます") + return point, -u + half = size / 2.0 + with np.errstate(divide="ignore"): + t_axes = np.where(np.abs(u) > 1e-12, half / np.abs(u), np.inf) + axis = int(np.argmin(t_axes)) + point = center + t_axes[axis] * u + normal = np.zeros(3) + normal[axis] = -np.sign(u[axis]) + others = [i for i in range(3) if i != axis] + margin = float(min(half - abs(point[i] - center[i]) for i in others)) + if margin < LED_RADIUS + EDGE_MARGIN: + raise ValueError(f"θ={theta_deg:g}: 立方体の稜に近すぎます（余裕 {margin:.1f}mm）") + return point, normal + + +def _board_points(grid: int = 41) -> np.ndarray: + """評価板のキャビティ側の面の格子点（四隅を含む）を返す.""" + gx = np.linspace(-BOARD_HALF[0], BOARD_HALF[0], grid) + gy = np.linspace(-BOARD_HALF[1], BOARD_HALF[1], grid) + mesh_x, mesh_y = np.meshgrid(gx, gy) + z = np.full(mesh_x.size, BOARD_SURFACE_Z) + return np.stack([mesh_x.ravel(), mesh_y.ravel(), z], axis=1) + + +def _disk_points(center: np.ndarray, direction: np.ndarray, radius: float, n: int) -> np.ndarray: + """direction を法線とする円板の外周点 n 個 + 中心を返す.""" + axis = np.array([0.0, 1.0, 0.0]) if abs(direction[1]) < 0.9 else np.array([0.0, 0.0, 1.0]) + u = np.cross(direction, axis) + u /= np.linalg.norm(u) + v = np.cross(direction, u) + angles = 0.1234 + np.linspace(0.0, 2.0 * np.pi, n, endpoint=False) + ring = center + radius * (np.outer(np.cos(angles), u) + np.outer(np.sin(angles), v)) + return np.vstack([center[None, :], ring]) + + +def led_geometry( + shape: str, size: float, theta_deg: float, azimuth_deg: float, aim: str +) -> dict: + """LED 円板の中心・発光方向・浮かせ量・傾き角と，板への直接光の有無を返す. + + 浮かせ量は「傾いた円板の縁が壁面へ食い込む量（LED_RADIUS×sinτ）+ マージン」 + （τ は発光方向と壁法線のなす角）．板向きは浮かせた後の中心から向け直す． + """ + wall, normal = wall_anchor(shape, size, theta_deg, azimuth_deg) + center = wall.copy() + for _ in range(2): + if aim == "board": + direction = BOARD_CENTER - center + elif aim == "normal": + direction = normal.copy() + elif aim == "down": + direction = np.array([0.0, 0.0, -1.0]) + else: + raise ValueError(f"未知の向き指定: {aim}") + direction = direction / np.linalg.norm(direction) + sin_tilt = float(np.linalg.norm(np.cross(direction, normal))) + standoff = LED_RADIUS * sin_tilt + STANDOFF_MARGIN + center = wall + standoff * normal + tilt_deg = float(np.degrees(np.arctan2(sin_tilt, float(direction @ normal)))) + front = ((_board_points() - center) @ direction) > 0.0 + return { + "center": center, "direction": direction, "normal": normal, + "standoff": standoff, "tilt_deg": tilt_deg, + "board_front_frac": float(front.mean()), + } + + +def baffle_radius(led: dict) -> float: + """発光面の前方にある板の全点から LED 全点を隠すバッフル半径（mm）を返す. + + LED 円板上の点 Q × 板上の点 P の直線がバッフル平面（LED から BAFFLE_GAP 先・ + 法線 = 発光方向）を横切る位置の，バッフル中心からの面内距離の最大値 × 安全率． + 発光面の背後にある P（直接光が届かない）は除外する． + """ + center, direction = led["center"], led["direction"] + b_center = center + BAFFLE_GAP * direction + q = _disk_points(center, direction, LED_RADIUS, 64) + p = _board_points() + delta = p[None, :, :] - q[:, None, :] # (M, N, 3) + denom = delta @ direction + front = denom > 1e-9 + if not front.any(): + return 0.0 + t = np.where(front, ((b_center - q) @ direction)[:, None] / np.where(front, denom, 1.0), 0.0) + x = q[:, None, :] + t[..., None] * delta + offset = x - b_center + radial = offset - (offset @ direction)[..., None] * direction[None, None, :] + return float(BAFFLE_MARGIN * np.linalg.norm(radial, axis=-1)[front].max()) + + +def _wall_distance(shape: str, size: float, points: np.ndarray) -> float: + """点群から解析的なキャビティ内壁までの最小距離（mm）を返す（負なら壁の外）.""" + if shape == "cube": + margins = np.concatenate([ + (size / 2.0 - np.abs(points[:, 0:2])).ravel(), + (points[:, 2] - (RIM_Z - size)), + (RIM_Z - points[:, 2]), + ]) + return float(margins.min()) + return float((size - np.linalg.norm(points - cavity_center(shape, size), axis=1)).min()) + + +def build_baffle_mesh( + shape: str, size: float, led: dict, dest: Path +) -> tuple[trimesh.Trimesh, float, float]: + """LED の視線上に置く傾いたバッフル円板を生成して書き出す. + + Returns: + (メッシュ, 半径 mm, 内壁への食い込み深さ mm)．食い込みは外殻の貫通も含めて許容する + （壁の外側に出た部分はキャビティ内から見えず，光輸送に影響しない）． + + Raises: + ValueError: 半径がキャビティ寸法に対して大きすぎる（物理的に成立しない）場合． + """ + sections = 128 + radius = baffle_radius(led) + if radius > BAFFLE_RADIUS_CAP_RATIO * size: + raise ValueError( + f"バッフル半径 {radius:.1f}mm がキャビティ寸法 {size:g}mm の " + f"{BAFFLE_RADIUS_CAP_RATIO:.0%} を超えます（板が発光面の前方に掠めて入る配置）" + ) + mesh = trimesh.creation.cylinder(radius=radius, height=BAFFLE_THICKNESS, sections=sections) + mesh.apply_transform(trimesh.transformations.rotation_matrix(np.pi / sections, [0, 0, 1])) + direction = led["direction"] + mesh.apply_transform(trimesh.geometry.align_vectors([0.0, 0.0, 1.0], direction)) + mesh.apply_translation(led["center"] + BAFFLE_GAP * direction) + embed = max(0.0, -_wall_distance(shape, size, mesh.vertices)) + dest.parent.mkdir(parents=True, exist_ok=True) + mesh.export(dest) + return mesh, radius, embed + + +def build_housing_mesh(shape: str, size: float, dest: Path) -> trimesh.Trimesh: + """筐体メッシュを生成して書き出す（軸上スイープと同一のビルダ）.""" + if dest.exists(): + return trimesh.load(dest, force="mesh") + if shape == "sphere": + return build_sphere_mesh(size, dest, SPHERE_CAVITY) + mesh = build_cube_housing( + size=size, thickness=THICKNESS, + board_half_x=BOARD_HALF[0], board_half_y=BOARD_HALF[1], + rim_inset=RIM_INSET, rim_z=RIM_Z, + ) + dest.parent.mkdir(parents=True, exist_ok=True) + mesh.export(dest) + return mesh + + +def housing_mesh_path(shape: str, size: float) -> Path: + """筐体 STL の置き場所（同じ形状・サイズの条件で共有）.""" + prefix = "r" if shape == "sphere" else "L" + return MESH_DIR / f"筐体_{shape}_{prefix}{size:g}.stl" + + +def build_scene_yaml( + cond: Condition, mesh_path: Path, leds: list[dict], baffles: list[Path | None], dest: Path +) -> None: + """テンプレートから当該条件のシーン YAML を書き出す. + + 材質・センサ・レンダ設定・LED の半径は変更しない．LED の位置・向き・個数と + バッフルだけを変え，多灯では radiance を 1/N にして総光束を保存する． + """ + config = copy.deepcopy(yaml.safe_load(TEMPLATE_YAML.read_text(encoding="utf-8"))) + config["name"] = f"led_placement_{cond.name}" + config["render"]["seed"] = cond.rep - 1 + config["description"] = ( + f"照明配置スイープ [{cond.shape} / {cond.size:g}mm]: θ={cond.theta_deg:g}° の壁面に " + f"{cond.n_led} 灯（向き={cond.aim}）．" + + ("バッフルで直接光を遮断．" if any(b is not None for b in baffles) + else "板は発光面の背後（直接光なし・バッフル無し）．") + ) + yaml_dir = dest.parent + for mesh in config["geometry"]["meshes"]: + if mesh["id"] == "sphere_housing": + if cond.shape == "cube": + mesh["id"] = "cube_housing" + mesh["path"] = relative_mesh_path(mesh_path, yaml_dir) + else: + mesh["path"] = relative_mesh_path(TEMPLATE_YAML.parent / mesh["path"], yaml_dir) + for k, baffle_path in enumerate(baffles): + if baffle_path is None: + continue + config["geometry"]["meshes"].append({ + "id": f"baffle_{k}", + "path": relative_mesh_path(baffle_path, yaml_dir), + "unit": "mm", + "material": "printer_white", + }) + template_light = config["lights"][0] + config["lights"] = [] + for k, led in enumerate(leds): + light = copy.deepcopy(template_light) + light["id"] = f"led_{k}" + light["center"] = [round(float(c), 4) for c in led["center"]] + light["direction"] = [round(float(d), 6) for d in led["direction"]] + light["radiance"] = [float(r) / cond.n_led for r in template_light["radiance"]] + config["lights"].append(light) + dest.parent.mkdir(parents=True, exist_ok=True) + dest.write_text(yaml.safe_dump(config, allow_unicode=True, sort_keys=False), encoding="utf-8") + + +def prepare_condition(cond: Condition) -> tuple[Path, dict]: + """1 条件の筐体・LED・バッフルを生成し，シーン YAML のパスと幾何の記録を返す. + + Raises: + ValueError: LED を置けない（稜・接続に近い），円板が壁を貫く，バッフルが大きすぎる． + """ + mesh_path = housing_mesh_path(cond.shape, cond.size) + housing = build_housing_mesh(cond.shape, cond.size, mesh_path) + if not housing.is_watertight: + raise ValueError(f"{cond.name}: watertight でないメッシュが生成されました") + leds: list[dict] = [] + baffles: list[Path | None] = [] + record: dict = { + "cavity_area_mm2": cavity_area(cond.shape, cond.size), + "baffled": 0, "baffle_radius_mm": 0.0, "baffle_embed_mm": 0.0, + "board_front_frac": 0.0, + } + for k, azimuth in enumerate(cond.azimuths_deg): + led = led_geometry(cond.shape, cond.size, cond.theta_deg, azimuth, cond.aim) + rim = _disk_points(led["center"], led["direction"], LED_RADIUS, 64) + # 球は解析的な球面に対して測る．θ=120° でも円板は loft 接続平面より下に収まり， + # 接続直上 12mm では loft が球の延長より外側（r=100）か 1.7mm 以内の内側なので + # 解析球面での判定が実質そのまま成り立つ（生成メッシュで確認済み） + clearance = _wall_distance(cond.shape, cond.size, rim) + if clearance < 0.0: + raise ValueError(f"{cond.name}: LED 円板が壁を貫いています（余裕 {clearance:.3f}mm）") + leds.append(led) + if led["board_front_frac"] > 0.0: + baffle_path = MESH_DIR / f"バッフル_{cond.base_name}_{k}.stl" + _, radius, embed = build_baffle_mesh(cond.shape, cond.size, led, baffle_path) + baffles.append(baffle_path) + record["baffled"] = 1 + record["baffle_radius_mm"] = round(radius, 3) + record["baffle_embed_mm"] = round(max(record["baffle_embed_mm"], embed), 3) + else: + baffles.append(None) + if k == 0: + record.update({ + "led_x_mm": round(float(led["center"][0]), 3), + "led_y_mm": round(float(led["center"][1]), 3), + "led_z_mm": round(float(led["center"][2]), 3), + "led_standoff_mm": round(float(led["standoff"]), 3), + "tilt_deg": round(led["tilt_deg"], 2), + "board_front_frac": round(led["board_front_frac"], 4), + "led_clearance_mm": round(clearance, 3), + }) + # 再レンダ条件（rep>1）は幾何が同じでも seed を変えた別 YAML にする．seed が同じだと + # シーンを読み直しても指標が 4 桁まで一致し，ノイズフロアにならない（サーバで実測） + scene_yaml = TMP_DIR / f"scene_{cond.name}.yaml" + build_scene_yaml(cond, mesh_path, leds, baffles, scene_yaml) + return scene_yaml, record + + +def prepare( + conds: list[Condition], log, report_suffix: str = "" +) -> dict[str, tuple[Path, dict]]: + """条件群の幾何を生成し，置けない条件はスキップして記録する. + + report_suffix は形状別に並列で走らせるときにレポートを分けるための識別子． + """ + prepared: dict[str, tuple[Path, dict]] = {} + skipped: list[tuple[str, str]] = [] + log(f"=== 幾何生成: {len(conds)} 条件 ===") + for cond in conds: + try: + scene_yaml, record = prepare_condition(cond) + except ValueError as exc: + skipped.append((cond.name, str(exc))) + log(f" SKIP {cond.name}: {exc}") + continue + prepared[cond.name] = (scene_yaml, record) + baffle_note = (f"バッフル R={record['baffle_radius_mm']:6.2f} " + f"食い込み={record['baffle_embed_mm']:5.2f}" if record["baffled"] + else "バッフル無し（板は背後）") + log(f" {cond.name:>34} LED=({record['led_x_mm']:7.2f},{record['led_y_mm']:7.2f}," + f"{record['led_z_mm']:8.2f}) 傾き={record['tilt_deg']:5.1f}° " + f"余裕={record['led_clearance_mm']:4.2f} 板前方={record['board_front_frac']:.2f} " + f"{baffle_note}") + report = TMP_DIR / f"prepare_report{report_suffix}.csv" + keys = ["name", "status", "note", "cavity_area_mm2", "led_x_mm", "led_y_mm", "led_z_mm", + "led_standoff_mm", "tilt_deg", "led_clearance_mm", "board_front_frac", + "baffled", "baffle_radius_mm", "baffle_embed_mm"] + with report.open("w", newline="", encoding="utf-8") as fh: + writer = csv.DictWriter(fh, fieldnames=keys) + writer.writeheader() + for cond in conds: + if cond.name in prepared: + row = {"name": cond.name, "status": "ok", "note": ""} + row.update(prepared[cond.name][1]) + else: + note = next(msg for name, msg in skipped if name == cond.name) + row = {"name": cond.name, "status": "skip", "note": note} + writer.writerow({k: row.get(k, "") for k in keys}) + log(f" 生成 {len(prepared)} / スキップ {len(skipped)}（{report}）") + return prepared + + +# --------------------------------------------------------------------------- +# 実行 +# --------------------------------------------------------------------------- + + +# 解決済みの Mitsuba variant（プロセス内で 1 回だけ決める．ワーカープロセスでも各 1 回） +_render_variant: str | None = None + + +def render_variant() -> str: + """このスイープが使う variant を返す（環境変数 SM_RENDER_VARIANT > テンプレート YAML）. + + 条件ごとにシーン YAML を読み直すと 1 条件あたり 1 秒近く掛かるので，全条件で共通の + テンプレート（TEMPLATE_YAML）の render.variant を 1 回だけ解決して使い回す． + """ + global _render_variant # noqa: PLW0603 + if _render_variant is None: + template = yaml.safe_load(TEMPLATE_YAML.read_text(encoding="utf-8")) + _render_variant = resolve_render_variant(str(template["render"]["variant"])) + return _render_variant + + +def ensure_render_variant() -> str: + """レンダ前にプロセスの variant を確定させる（遮蔽検証も同じ variant で走らせるため）. + + check_occlusion は mi.variant() が設定済みならそれを使うので，先に呼んでおくと + YAML の llvm_ad_rgb ではなくサーバの cuda_ad_rgb でレイキャストが走る + （球 1 条件で 25 秒 → 1 秒未満）．import mitsuba を伴うので， + 既存結果を再利用するだけの条件では呼ばない． + """ + # mitsuba を読み込む import なので，モジュールレベルではなくここで行う + from sm_tias_lightsim.render.runner import ensure_variant + + variant = render_variant() + ensure_variant(variant) + return variant + + +def verify_occlusion(scene_yaml: Path, messages: list[str]) -> dict[str, float]: + """バッフル付き条件について板→LED の遮蔽率をレイキャストで測る（失敗しても続行）. + + Dr.Jit のスレッド数は ensure_render_variant 内の apply_drjit_thread_count が + SM_DRJIT_THREADS を見て制限するので，ここでは触らない． + """ + try: + from verify_baffle_occlusion import check_occlusion + + result = check_occlusion(scene_yaml, board_grid=41, led_count=32) + return { + "occlusion_blocked_frac": round(result["blocked"] / result["total"], 6), + "occlusion_total": int(result["total"]), + } + except Exception as exc: # noqa: BLE001 - 検証の失敗でスイープを止めない + messages.append(f" !! 遮蔽検証に失敗（続行）: {exc}") + return {"occlusion_blocked_frac": float("nan"), "occlusion_total": 0} + + +def run_job(cond: Condition, scene_yaml: Path, record: dict) -> tuple[dict | None, list[str]]: + """1 条件をレンダ→評価まで実行し row.json を書いて集約行とログ行を返す. + + ワーカープロセスからも呼べるよう，引数・戻り値は pickle 可能なものだけにする + （主プロセスの log クロージャは渡せないので，ログはメッセージとして返して + 主プロセスが書き出す）． + + Returns: + (集約行, ログ行のリスト)．失敗した条件の集約行は None． + """ + messages: list[str] = [] + out_dir = OUTPUT_BASE / cond.name + row_json = out_dir / "row.json" + if row_json.exists(): + messages.append(f"--- {cond.name}: 既存結果を再利用") + return json.loads(row_json.read_text(encoding="utf-8")), messages + variant = ensure_render_variant() + row: dict = { + "name": cond.name, "group": cond.group, "shape": cond.shape, "size_mm": cond.size, + "theta_deg": cond.theta_deg, "aim": cond.aim, "n_led": cond.n_led, "rep": cond.rep, + } + row.update(record) + row["render_variant"] = variant + if record["baffled"]: + row.update(verify_occlusion(scene_yaml, messages)) + messages.append(f"--- {cond.name}: レンダ開始") + out_dir.mkdir(parents=True, exist_ok=True) + start = time.time() + # sm-render / sm-evaluate をサブプロセスで呼ぶ代わりにプロセス内で実行する + # （起動 + import mitsuba + 図の描画で 1 条件あたり 4 秒前後を使っていたため）． + # 標準出力は条件ごとの評価サマリーで埋まるので取り込んでおき，失敗時だけ末尾を残す． + captured = io.StringIO() + try: + with redirect_stdout(captured): + render_yaml(scene_yaml, out_dir, variant) + evaluate_image( + out_dir / "luminance.npy", EVAL_CONFIG_YAML, out_dir, save_figures=False, + ) + except Exception as exc: # noqa: BLE001 - 1 条件の失敗でスイープ全体を止めない + messages.append(f" !! 失敗: {exc}") + tail = captured.getvalue().strip()[-800:] + if tail: + messages.append(f" 出力: {tail}") + return None, messages + row["elapsed_s"] = round(time.time() - start, 1) + row.update(read_metrics_csv(out_dir / "metrics.csv")) + with (out_dir / "metrics.csv").open(encoding="utf-8") as fh: + metrics = next(csv.DictReader(fh)) + for key in ("intensity_centroid_x_norm", "intensity_centroid_y_norm"): + row[key] = float(metrics[key]) if metrics.get(key) else float("nan") + row.update(raw_stats(out_dir / "luminance.npy")) + row_json.write_text(json.dumps(row, ensure_ascii=False, indent=1), encoding="utf-8") + occl = row.get("occlusion_blocked_frac", float("nan")) + cx, cy = row["intensity_centroid_x_norm"], row["intensity_centroid_y_norm"] + messages.append( + f" -> {row['elapsed_s']:.0f}s CoV={row['cov']:.4f} centroid=({cx:+.4f},{cy:+.4f}) " + f"angular_ratio={row['angular_ratio']:.4f} 生mean={row['raw_mean_roi']:.4f} " + f"遮蔽={occl:.4f}") + return row, messages + + +def collect_rows() -> list[dict]: + """出力ディレクトリ配下の row.json をすべて読む（並列ジョブの結果を合流させる）.""" + rows = [] + for path in sorted(OUTPUT_BASE.glob("*/row.json")): + rows.append(json.loads(path.read_text(encoding="utf-8"))) + return rows + + +def write_summary() -> Path: + """全 row.json から集約 CSV を書き直す（並列に走る形状別ジョブでも壊れない）.""" + order = {name: i for i, name in enumerate(c.name for c in conditions())} + rows = sorted(collect_rows(), key=lambda r: order.get(r["name"], len(order))) + summary_csv = OUTPUT_BASE / "summary.csv" + summary_csv.parent.mkdir(parents=True, exist_ok=True) + with summary_csv.open("w", newline="", encoding="utf-8") as fh: + writer = csv.DictWriter(fh, fieldnames=SUMMARY_COLUMNS) + writer.writeheader() + for row in rows: + writer.writerow({k: row.get(k, "") for k in SUMMARY_COLUMNS}) + return summary_csv + + +def main() -> int: + """照明配置スイープを実行するメインエントリポイント.""" + parser = argparse.ArgumentParser(description="照明配置スイープ（球 + 立方体）") + parser.add_argument("--prepare-only", action="store_true", + help="メッシュとシーン YAML の生成だけ行う（3D 配置確認用）") + parser.add_argument("--shapes", choices=("sphere", "cube"), nargs="*", default=None, + help="実行する形状（省略時は両方）") + parser.add_argument("--groups", choices=GROUPS, nargs="*", default=None, + help="実行する条件グループ（省略時は全部）") + parser.add_argument("--limit", type=int, default=None, + help="先頭 N 条件だけ実行する（動作確認用）") + parser.add_argument("--workers", type=int, default=1, + help="同時に実行する条件数（既定 1＝逐次．GPU 1 枚なら 2〜3 が目安）") + parser.add_argument("--list", action="store_true", help="条件一覧を表示して終了") + args = parser.parse_args() + + conds = [c for c in conditions() + if (not args.shapes or c.shape in args.shapes) + and (not args.groups or c.group in args.groups)] + if args.limit: + conds = conds[: args.limit] + if args.list: + for c in conds: + print(c.name) + print(f"{len(conds)} 条件") + return 0 + + log_path = ROOT / "output" / "_sweep_led_placement.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"=== 照明配置スイープ（{len(conds)} 条件）===") + suffix = "_" + "_".join(args.shapes) if args.shapes else "" + prepared = prepare(conds, log, suffix) + if args.prepare_only: + log(f"--prepare-only のため終了．シーン YAML: {TMP_DIR}") + log_file.close() + return 0 + + jobs = [(cond, *prepared[cond.name]) for cond in conds if cond.name in prepared] + done = 0 + if args.workers > 1: + # ワーカーは spawn（Windows 互換）で立ち上げ，各プロセスが自前で variant を + # 1 回設定してから条件を処理する．集約 CSV の書き出しは主プロセスだけが行う． + log(f"=== {args.workers} 並列で実行（1 ワーカー 1 プロセス）===") + context = multiprocessing.get_context("spawn") + with ProcessPoolExecutor(max_workers=args.workers, mp_context=context) as pool: + futures = { + pool.submit(run_job, cond, scene_yaml, record): cond + for cond, scene_yaml, record in jobs + } + for future in as_completed(futures): + cond = futures[future] + try: + row, messages = future.result() + except Exception as exc: # noqa: BLE001 - 1 条件の異常で全体を止めない + log(f" !! {cond.name}: ワーカーが異常終了（続行）: {exc}") + continue + for message in messages: + log(message) + if row is not None: + done += 1 + write_summary() + else: + for cond, scene_yaml, record in jobs: + row, messages = run_job(cond, scene_yaml, record) + for message in messages: + log(message) + if row is not None: + done += 1 + write_summary() + + log(f"\n=== 完了: {done} / {len(prepared)} 条件 ===") + log(f" 集約 CSV: {write_summary()}") + log_file.close() + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/sm_tias_lightsim/cli/evaluate.py b/src/sm_tias_lightsim/cli/evaluate.py index 5a5d4ba..c642f35 100644 --- a/src/sm_tias_lightsim/cli/evaluate.py +++ b/src/sm_tias_lightsim/cli/evaluate.py @@ -3,6 +3,10 @@ PNG または luminance.npy を入力として均一性評価を実行し， metrics.csv と可視化 PNG を出力ディレクトリに保存する． +処理本体は `evaluation/pipeline.py` の `evaluate_image()` にある（スイープ等から +プロセス内で直接呼べるようにするため）．本モジュールは引数のパースと， +例外 → 「エラー: …」表示 + 終了コード 1 への変換だけを担う． + matplotlib のバックエンドは Agg を強制する（headless 環境対応）． evaluate は Mitsuba を使わないため，render.py と異なり遅延 import は不要． ただし Agg バックエンドの設定順序を保つため，matplotlib.use() は pyplot import より前に置く． @@ -17,54 +21,37 @@ import matplotlib # headless 環境（サーバ・CI 等）でも動くよう Agg バックエンドを強制する． -# visualize.py より先に設定することで，pyplot import 時に Agg が有効になる． +# pipeline（→ visualize.py）より先に設定することで，pyplot import 時に Agg が有効になる． matplotlib.use("Agg") -import numpy # noqa: E402 -import pandas # noqa: E402 -from PIL import Image # noqa: E402 - -from sm_tias_lightsim.evaluation.normalize import apply_mean_anchor # noqa: E402 -from sm_tias_lightsim.evaluation.oetf import apply_srgb_oetf, solve_exposure # noqa: E402 -from sm_tias_lightsim.evaluation.radial import ( # noqa: E402 +# evaluate_image / EvaluationStepError 以外は，補助関数が本モジュールにあった頃の +# 参照互換のための再エクスポート（テスト・スクリプトが cli.evaluate.<名前> を参照している）． +from sm_tias_lightsim.evaluation.pipeline import ( # noqa: E402, F401 + _FIGURES_DIRNAME, + _HISTOGRAM_FILENAME, + _LUMINANCE_MAP_FILENAME, + _METRICS_FILENAME, + _RADIAL_PROFILE_FILENAME, + _ZONE_MAP_FILENAME, + EvaluationStepError, + _apply_normalization_to_luminance, + _apply_oetf_to_luminance, + _load_eval_config, + _load_luminance, + _load_reference_luminance, + _print_summary, + _run_radial_min_max, + _run_radial_profile, + _run_shape_metrics, + _run_zone_analysis, + _save_figures, + _save_metrics_csv, + evaluate_image, +) +from sm_tias_lightsim.evaluation.radial import ( # noqa: E402, F401 calc_radial_min_max_ratio, calc_radial_profile, ) -from sm_tias_lightsim.evaluation.roi import extract_roi # noqa: E402 -from sm_tias_lightsim.evaluation.shape_metrics import ( # noqa: E402 - calc_angular_radial_decomposition, - calc_intensity_centroid, - calc_intensity_pca, - calc_ncc, - calc_ssim, -) -from sm_tias_lightsim.evaluation.visualize import ( # noqa: E402 - save_histogram, - save_luminance_map, - save_radial_profile, - save_zone_map, -) -from sm_tias_lightsim.evaluation.zones import ( # noqa: E402 - calc_zone_stats, - classify_zones, - compute_normalized_distance_map, -) -from sm_tias_lightsim.scene.schema import ( # noqa: E402 - EvaluationConfig, - MetricsConfig, - NormalizationConfig, - OetfConfig, - ZonesConfig, -) -from sm_tias_lightsim.uniformity import calc_uniformity, to_grayscale # noqa: E402 - -# 出力ファイル名を一元管理する（変更時の影響範囲を最小化するため） -_METRICS_FILENAME = "metrics.csv" -_FIGURES_DIRNAME = "figures" -_LUMINANCE_MAP_FILENAME = "luminance_map.png" -_HISTOGRAM_FILENAME = "histogram.png" -_RADIAL_PROFILE_FILENAME = "radial_profile.png" -_ZONE_MAP_FILENAME = "zone_map.png" def _parse_args(argv: list[str] | None = None) -> argparse.Namespace: @@ -106,431 +93,19 @@ metavar="OUTPUT_DIR", help="出力先ディレクトリ（存在しない場合は親ごと作成する）．", ) + parser.add_argument( + "--no-figures", + action="store_true", + help="figures/ の可視化 PNG を作らない（指標だけ欲しいスイープ用．描画分だけ速い）．", + ) return parser.parse_args(argv) -def _apply_oetf_to_luminance( - roi_luminance: numpy.ndarray, - oetf_cfg: OetfConfig, -) -> numpy.ndarray: - """OETF 順方向適用：ROI linear luminance を sRGB 符号化空間へ変換する． - - 無彩色シーン（シミュ LED・材料が RGB 等価）では luminance への直接 OETF 適用は - 「チャネル別 sRGB OETF → Rec.709 luma 加重和」と厳密一致する（SPEC_03 参照）． - - 露出ソルバは ROI 上で解く（target_mean は実測 ROI mean のため ROI 上で収束させる）． - - Args: - roi_luminance: ROI 切り出し済みの linear 輝度配列（H' x W', float64）． - oetf_cfg: OetfConfig（enabled=True が前提）． - - Returns: - OETF 適用後の ROI 輝度配列（H' x W', float64, quantize_uint8=True なら [0,255]）． - - Raises: - ValueError: oetf_cfg.exposure が None の場合（anchor_mean モードで必須）． - """ - if oetf_cfg.exposure is None: - msg = "oetf.enabled=True かつ anchor_mean モードでは exposure 設定が必須です．" - raise ValueError(msg) - - roi_lin = roi_luminance.astype(numpy.float64) - target_mean = oetf_cfg.exposure.target_mean - # ROI linear 値で露出スケールを決定（実測 mean が ROI mean のため ROI 上で解く） - scale, clipped = solve_exposure(roi_lin, target_mean) - clip_pct = clipped * 100.0 - print(f" OETF exposure : scale={scale:.6g} clip={clip_pct:.1f}%") - - # ROI に対してのみ OETF を適用する（指標算出は ROI 上で行うため） - return apply_srgb_oetf(roi_lin, scale, quantize_uint8=oetf_cfg.quantize_uint8) - - -def _apply_normalization_to_luminance( - roi_luminance: numpy.ndarray, - norm_cfg: NormalizationConfig, -) -> numpy.ndarray: - """linear mean anchor 正規化を ROI 輝度配列に適用する． - - ROI mean を norm_cfg.target_mean に揃える線形ゲインを適用する． - 形状指標（CoV / 勾配 / C/P 比 / max-min 比）は線形ゲインに不変なので， - 絶対値スケールだけが変わる． - - Args: - roi_luminance: ROI 切り出し済みの linear 輝度配列（H' x W', float64）． - norm_cfg: NormalizationConfig（enabled=True が前提）． - - Returns: - 正規化後の ROI 輝度配列（H' x W', float64）． - - Raises: - ValueError: ROI mean が 0 の場合（apply_mean_anchor から伝播）． - """ - roi_lin = roi_luminance.astype(numpy.float64) - target_mean = norm_cfg.target_mean - normalized, alpha = apply_mean_anchor(roi_lin, target_mean) - print(f" normalization : alpha={alpha:.6g} target_mean={target_mean}") - return normalized - - -def _load_luminance(image_path: Path) -> numpy.ndarray: - """画像ファイルを読み込み，輝度配列（2D float64）を返す． - - 拡張子が .npy の場合は numpy.load で読み込む． - それ以外（PNG 等）は PIL で RGB 画像として読み込み，to_grayscale で変換する． - .npy が 3D（H, W, 3）だった場合も to_grayscale を通す． - - Args: - image_path: 読み込む画像ファイルのパス． - - Returns: - 輝度配列（H x W, float64）． - - Raises: - FileNotFoundError: ファイルが存在しない場合． - ValueError: ファイルを輝度配列に変換できない場合． - """ - if not image_path.exists(): - msg = f"画像ファイルが見つかりません: {image_path}" - raise FileNotFoundError(msg) - - if image_path.suffix.lower() == ".npy": - data = numpy.load(image_path) - if data.ndim == 3: - return to_grayscale(data.astype(numpy.float64)) - if data.ndim == 2: - return data.astype(numpy.float64) - msg = f".npy の次元数が不正です（ndim={data.ndim}）．2D または 3D を期待します．" - raise ValueError(msg) - - pil_image = Image.open(image_path).convert("RGB") - return to_grayscale(numpy.array(pil_image, dtype=numpy.float64)) - - -def _load_eval_config(config_path: Path) -> EvaluationConfig: - """YAML を読み込んで EvaluationConfig を返す． - - SceneConfig の evaluation サブツリー，または EvaluationConfig 互換の - トップレベル YAML の両方をサポートする． - - Args: - config_path: YAML ファイルのパス． - - Returns: - EvaluationConfig インスタンス． - - Raises: - FileNotFoundError: ファイルが存在しない場合． - pydantic.ValidationError: スキーマ違反． - yaml.YAMLError: YAML 構文エラー． - """ - import yaml - - if not config_path.exists(): - msg = f"設定 YAML ファイルが見つかりません: {config_path}" - raise FileNotFoundError(msg) - - with config_path.open(encoding="utf-8") as f: - data = yaml.safe_load(f) - - # "evaluation" キーがあれば SceneConfig 埋め込みとみなしてサブツリーを取得する - if isinstance(data, dict) and "evaluation" in data: - eval_data = data["evaluation"] - else: - eval_data = data - - return EvaluationConfig.model_validate(eval_data) - - -def _run_zone_analysis( - roi_luminance: numpy.ndarray, - zones_config: ZonesConfig | None, -) -> tuple[numpy.ndarray, dict] | tuple[None, dict]: - """ゾーン分析を実行し，zone_map とゾーン統計を返す． - - Args: - roi_luminance: ROI 輝度配列（H x W, float64）． - zones_config: ZonesConfig インスタンス（enabled / boundaries を持つ）． - None の場合は無効扱い． - - Returns: - (zone_map, zone_stats) のタプル． - zones が無効なら (None, {}) を返す． - - Raises: - ValueError: ゾーン分析の計算に失敗した場合． - """ - if zones_config is None or not zones_config.enabled: - return None, {} - - h, w = roi_luminance.shape - distance_map = compute_normalized_distance_map(h, w) - zone_map = classify_zones(distance_map, zones_config.boundaries) - zone_stats = calc_zone_stats(roi_luminance, zone_map) - return zone_map, zone_stats - - -def _run_radial_profile( - roi_luminance: numpy.ndarray, - metrics_config: MetricsConfig | None, -) -> numpy.ndarray | None: - """放射状プロファイルを算出して返す． - - Args: - roi_luminance: ROI 輝度配列（H x W, float64）． - metrics_config: MetricsConfig インスタンス（radial_profile フラグを持つ）． - None の場合は無効扱い． - - Returns: - 放射状プロファイル配列（(n_bins, 2)），無効時は None． - """ - if metrics_config is None or not metrics_config.radial_profile: - return None - return calc_radial_profile(roi_luminance) - - -def _run_radial_min_max( - roi_luminance: numpy.ndarray, - metrics_config: MetricsConfig | None, - radial_profile: numpy.ndarray | None, -) -> dict: - """動径 min/max 比を算出して返す． - - radial_min_max が True のとき，radial_profile が既に算出されていればそれを再利用し， - 未算出（None）の場合はその場で calc_radial_profile を呼び出して算出する． - 重複算出を避けるため，radial_profile の結果を優先的に使用する． - - Args: - roi_luminance: ROI 輝度配列（H x W, float64）． - metrics_config: MetricsConfig インスタンス（radial_min_max フラグを持つ）． - None の場合は無効扱い． - radial_profile: _run_radial_profile で算出済みのプロファイル配列． - None の場合はここで算出する． - - Returns: - 動径 min/max 比指標の辞書（radial_min_max_ratio / radial_min_distance / - radial_min / radial_max の 4 キー）．無効時は空辞書． - """ - if metrics_config is None or not metrics_config.radial_min_max: - return {} - - # radial_profile が未算出の場合はここで算出する（重複算出を避けるため既存結果を優先） - profile = radial_profile if radial_profile is not None else calc_radial_profile(roi_luminance) - return calc_radial_min_max_ratio(profile) - - -def _load_reference_luminance(reference_path: Path) -> numpy.ndarray: - """形状指標用の参照画像を読み込み，輝度配列（2D float64）を返す． - - SSIM / NCC の参照画像として使用する．_load_luminance と同一の変換ロジックを用い， - .npy 形式（2D/3D）と PNG 形式の両方に対応する． - 呼び出し元（_run_shape_metrics）で shape を確認してから calc_ssim / calc_ncc に渡すこと． - - Args: - reference_path: 参照画像ファイルのパス（.npy または .png）． - - Returns: - 輝度配列（H x W, float64）． - - Raises: - FileNotFoundError: ファイルが存在しない場合． - ValueError: ファイルを輝度配列に変換できない場合． - """ - return _load_luminance(reference_path) - - -def _run_shape_metrics( - roi_luminance: numpy.ndarray, - metrics_config: MetricsConfig | None, - config_dir: Path, -) -> dict: - """形状指標（輝度重心 / 角度-放射分解 / PCA 主軸 / SSIM / NCC）を算出する． - - metrics.shape_metrics=True のとき，以下の指標を算出して辞書として返す: - - 輝度重心 5 値 (calc_intensity_centroid) - - 角度-放射分解 3 値 (calc_angular_radial_decomposition) - - PCA 主軸 4 値 (calc_intensity_pca) - - SSIM / NCC（shape_reference が指定されたときのみ，それ以外は float("nan")） - - shape_reference パスは config YAML のディレクトリを基準に解決する． - 絶対パスの場合はそのまま使用する． - - Args: - roi_luminance: ROI 輝度配列（H x W, float64）． - metrics_config: MetricsConfig インスタンス（shape_metrics フラグを持つ）． - None の場合は無効扱い． - config_dir: 評価 YAML ファイルのディレクトリ（shape_reference の相対パス解決用）． - - Returns: - 形状指標の辞書（12 キー）．無効時は空辞書． - SSIM / NCC が計算不能な場合は float("nan") を格納する． - """ - if metrics_config is None or not metrics_config.shape_metrics: - return {} - - result: dict = {} - - # 輝度重心 - centroid = calc_intensity_centroid(roi_luminance) - result["intensity_centroid_x_px"] = centroid["centroid_x_px"] - result["intensity_centroid_y_px"] = centroid["centroid_y_px"] - result["intensity_centroid_x_norm"] = centroid["centroid_x_norm"] - result["intensity_centroid_y_norm"] = centroid["centroid_y_norm"] - result["intensity_centroid_offset_norm"] = centroid["centroid_offset_norm"] - - # 角度-放射分解 - ard = calc_angular_radial_decomposition(roi_luminance) - result["angular_variance"] = ard["angular_variance"] - result["radial_variance"] = ard["radial_variance"] - result["angular_ratio"] = ard["angular_ratio"] - - # PCA 主軸 - pca = calc_intensity_pca(roi_luminance) - result["pca_primary_axis_angle_deg"] = pca["primary_axis_angle_deg"] - result["pca_eccentricity"] = pca["eccentricity"] - - # SSIM / NCC（shape_reference が指定されている場合のみ） - reference_path_str = metrics_config.shape_reference - if reference_path_str is not None: - ref_path = Path(reference_path_str) - if not ref_path.is_absolute(): - ref_path = config_dir / ref_path - - try: - ref_luminance = _load_reference_luminance(ref_path) - # shape が一致するか確認（一致しない場合は警告してスキップ） - if ref_luminance.shape != roi_luminance.shape: - print( - f" 警告: shape_reference の shape {ref_luminance.shape} が " - f"ROI の shape {roi_luminance.shape} と異なるため SSIM/NCC をスキップします．" - ) - result["ssim"] = float("nan") - result["ncc"] = float("nan") - else: - result["ssim"] = calc_ssim(roi_luminance, ref_luminance) - result["ncc"] = calc_ncc(roi_luminance, ref_luminance) - except FileNotFoundError: - print( - f" 警告: shape_reference ファイルが見つかりません: {ref_path}．" - "SSIM/NCC をスキップします．" - ) - result["ssim"] = float("nan") - result["ncc"] = float("nan") - except Exception as exc: # noqa: BLE001 - print( - f" 警告: shape_reference の読み込みに失敗しました — {exc}．" - "SSIM/NCC をスキップします．" - ) - result["ssim"] = float("nan") - result["ncc"] = float("nan") - else: - result["ssim"] = float("nan") - result["ncc"] = float("nan") - - return result - - -def _save_metrics_csv(metrics: dict, csv_path: Path) -> None: - """均一性指標を CSV ファイルとして保存する． - - Args: - metrics: 指標名 → 値の辞書． - csv_path: 保存先 CSV パス． - """ - df = pandas.DataFrame([metrics]) - df.to_csv(csv_path, index=False) - - -def _save_figures( - roi_luminance: numpy.ndarray, - figures_dir: Path, - radial_profile: numpy.ndarray | None, - zone_map: numpy.ndarray | None, -) -> None: - """可視化 PNG を figures/ ディレクトリに保存する． - - Args: - roi_luminance: ROI 輝度配列（H x W, float64）． - figures_dir: 保存先ディレクトリ（事前に mkdir 済みであること）． - radial_profile: 放射状プロファイル配列，None の場合は保存しない． - zone_map: ゾーンマップ配列，None の場合は保存しない． - """ - save_luminance_map(roi_luminance, figures_dir / _LUMINANCE_MAP_FILENAME) - save_histogram(roi_luminance, figures_dir / _HISTOGRAM_FILENAME) - - if radial_profile is not None: - save_radial_profile(radial_profile, figures_dir / _RADIAL_PROFILE_FILENAME) - - if zone_map is not None: - save_zone_map(zone_map, figures_dir / _ZONE_MAP_FILENAME, luminance=roi_luminance) - - -def _print_summary(metrics: dict, csv_path: Path, figures_dir: Path) -> None: - """評価サマリーを標準出力に表示する． - - Args: - metrics: 指標名 → 値の辞書． - csv_path: metrics.csv の保存先パス（表示用）． - figures_dir: figures/ ディレクトリのパス（表示用）． - """ - print("=== 均一性評価サマリー ===") - print(f" CoV : {metrics['cov']:.6f}") - print(f" std : {metrics['std']:.4f}") - print(f" max/min ratio : {metrics['max_min_ratio']:.4f}") - print(f" mean : {metrics['mean']:.4f}") - print(f" max : {metrics['max']:.4f}") - print(f" min : {metrics['min']:.4f}") - - if "center_periphery_ratio" in metrics: - print(f" C/P ratio : {metrics['center_periphery_ratio']:.4f}") - print(f" gradient [%] : {metrics['gradient_magnitude']:.2f}") - - if "radial_min_max_ratio" in metrics: - print(f" radial min/max : {metrics['radial_min_max_ratio']:.4f}") - print(f" radial min dist: {metrics['radial_min_distance']:.4f}") - - if "intensity_centroid_offset_norm" in metrics: - print(f" centroid offset: {metrics['intensity_centroid_offset_norm']:.4f} (norm)") - print( - f" centroid (norm): x={metrics['intensity_centroid_x_norm']:.4f}" - f" y={metrics['intensity_centroid_y_norm']:.4f}" - ) - ang_ratio = metrics.get("angular_ratio", float("nan")) - if ang_ratio != ang_ratio: # NaN check - print(" angular_ratio : NaN (uniform image)") - else: - print(f" angular_ratio : {ang_ratio:.4f} (0=concentric, 1=directional)") - print(f" pca axis [deg] : {metrics['pca_primary_axis_angle_deg']:.1f}") - print(f" pca eccentricit: {metrics['pca_eccentricity']:.4f}") - ssim_val = metrics.get("ssim", float("nan")) - ncc_val = metrics.get("ncc", float("nan")) - if ssim_val == ssim_val: # not NaN - print(f" SSIM : {ssim_val:.4f}") - print(f" NCC : {ncc_val:.4f}") - - print(f"\nmetrics.csv : {csv_path}") - print(f"figures/ : {figures_dir}") - - def main(argv: list[str] | None = None) -> int: """sm-evaluate コマンドのメインエントリポイント． - 処理フロー: - 1. コマンドライン引数をパースする． - 2. 画像ファイルを読み込み，輝度配列を得る． - 3. YAML 設定ファイルを読み込み，EvaluationConfig を得る． - 4. 出力ディレクトリと figures/ サブディレクトリを作成する． - 5. extract_roi で ROI を切り出す． - 6. normalization.enabled のとき linear mean anchor 正規化を適用する． - 7. oetf.enabled のとき sRGB OETF を適用する． - 8. calc_uniformity で均一性指標を算出する． - 9. zones.enabled のときゾーン分析を実行する． - 10. metrics.radial_profile のとき放射状プロファイルを算出する． - 11. metrics.radial_min_max のとき動径 min/max 比を算出する（radial_profile 再利用）． - 12. metrics.shape_metrics のとき形状指標を算出する - （重心 / 角度-放射分解 / PCA / SSIM / NCC，metrics.csv に 12 列末尾追加）． - 13. metrics.csv を保存する． - 14. 可視化 PNG を figures/ に保存する． - 15. 算出指標のサマリーを標準出力に表示する． + 引数をパースして `evaluate_image()` を呼び，段階別のエラーを終了コードに変換する． + 評価そのものの処理フローは `evaluation/pipeline.py` の `evaluate_image()` を参照． Args: argv: 引数リスト．None のとき sys.argv[1:] を使用． @@ -542,112 +117,15 @@ # ここでは ValueError / FileNotFoundError 等の実行時エラーのみ捕捉する args = _parse_args(argv) - image_path: Path = args.image.resolve() - config_path: Path = args.config.resolve() - output_dir: Path = args.output.resolve() - try: - luminance = _load_luminance(image_path) - except FileNotFoundError as exc: - print(f"エラー: {exc}", file=sys.stderr) - return 1 - except Exception as exc: # noqa: BLE001 - print(f"エラー: 画像の読み込みに失敗しました — {exc}", file=sys.stderr) - return 1 - - try: - eval_config = _load_eval_config(config_path) - except FileNotFoundError as exc: - print(f"エラー: {exc}", file=sys.stderr) - return 1 - except Exception as exc: # noqa: BLE001 - print(f"エラー: 設定 YAML の読み込みに失敗しました — {exc}", file=sys.stderr) - return 1 - - figures_dir = output_dir / _FIGURES_DIRNAME - try: - output_dir.mkdir(parents=True, exist_ok=True) - figures_dir.mkdir(parents=True, exist_ok=True) - except OSError as exc: - print(f"エラー: 出力ディレクトリを作成できません — {exc}", file=sys.stderr) - return 1 - - try: - roi_luminance = extract_roi(luminance, eval_config.roi) - except (ValueError, NotImplementedError) as exc: - print(f"エラー: ROI 抽出に失敗しました — {exc}", file=sys.stderr) - return 1 - - # linear mean anchor 正規化（enabled のとき ROI 抽出後・OETF 前に挟む）． - # OETF との併用禁止は EvaluationConfig の model_validator で保証されている． - norm_cfg = eval_config.normalization - if norm_cfg is not None and norm_cfg.enabled: - try: - roi_luminance = _apply_normalization_to_luminance(roi_luminance, norm_cfg) - except ValueError as exc: - print(f"エラー: linear 正規化に失敗しました — {exc}", file=sys.stderr) - return 1 - - # OETF 順方向適用（enabled のとき ROI 抽出後・指標算出前に挟む）． - # eval_real.yaml（実測 PNG）は oetf を設定しないこと： - # 実測 PNG は既に sRGB 符号化済みのため，ここで OETF を掛けると二重適用になる． - oetf_cfg = eval_config.oetf - if oetf_cfg is not None and oetf_cfg.enabled: - try: - roi_luminance = _apply_oetf_to_luminance(roi_luminance, oetf_cfg) - except ValueError as exc: - print(f"エラー: OETF 適用に失敗しました — {exc}", file=sys.stderr) - return 1 - - try: - uniformity_stats = calc_uniformity(roi_luminance) - except ValueError as exc: - print(f"エラー: 均一性指標の算出に失敗しました — {exc}", file=sys.stderr) - return 1 - - metrics: dict = dict(uniformity_stats) - - try: - zone_map, zone_stats = _run_zone_analysis(roi_luminance, eval_config.zones) - metrics.update(zone_stats) - except ValueError as exc: - print(f"エラー: ゾーン分析に失敗しました — {exc}", file=sys.stderr) - return 1 - - try: - radial_profile = _run_radial_profile(roi_luminance, eval_config.metrics) - except Exception as exc: # noqa: BLE001 - print(f"エラー: 放射状プロファイルの算出に失敗しました — {exc}", file=sys.stderr) - return 1 - - try: - radial_min_max_stats = _run_radial_min_max( - roi_luminance, eval_config.metrics, radial_profile + evaluate_image( + args.image.resolve(), + args.config.resolve(), + args.output.resolve(), + save_figures=not args.no_figures, ) - metrics.update(radial_min_max_stats) - except Exception as exc: # noqa: BLE001 - print(f"エラー: 動径 min/max 比の算出に失敗しました — {exc}", file=sys.stderr) + except EvaluationStepError as exc: + print(f"エラー: {exc}", file=sys.stderr) return 1 - try: - shape_stats = _run_shape_metrics(roi_luminance, eval_config.metrics, config_path.parent) - metrics.update(shape_stats) - except Exception as exc: # noqa: BLE001 - print(f"エラー: 形状指標の算出に失敗しました — {exc}", file=sys.stderr) - return 1 - - csv_path = output_dir / _METRICS_FILENAME - try: - _save_metrics_csv(metrics, csv_path) - except Exception as exc: # noqa: BLE001 - print(f"エラー: metrics.csv の保存に失敗しました — {exc}", file=sys.stderr) - return 1 - - try: - _save_figures(roi_luminance, figures_dir, radial_profile, zone_map) - except Exception as exc: # noqa: BLE001 - print(f"エラー: 可視化 PNG の保存に失敗しました — {exc}", file=sys.stderr) - return 1 - - _print_summary(metrics, csv_path, figures_dir) return 0 diff --git a/src/sm_tias_lightsim/cli/render.py b/src/sm_tias_lightsim/cli/render.py index 5512dde..8f4e8e5 100644 --- a/src/sm_tias_lightsim/cli/render.py +++ b/src/sm_tias_lightsim/cli/render.py @@ -1,13 +1,18 @@ # pyright: reportCallIssue=none # Mitsuba 3 は variant ベースの動的属性解決を行うため Pylance/Pyright が # mi.* を object 型と推論し reportCallIssue を誤検知する（coding-style.md の方針）． -# 下記の Mitsuba 依存 import（load_scene_yaml / render_scene / save_png / save_luminance）を -# main() 内で遅延 import しているのは，--help だけ使う場面（Mitsuba 不在の軽量環境）でも -# 起動できるようにするためである． +# 処理本体（render_yaml）を main() 内で遅延 import しているのは，--help だけ使う場面 +# （Mitsuba 不在の軽量環境）でも起動できるようにするためである． +# render/pipeline.py 自体も mitsuba 依存 import を関数内に閉じ込めており， +# llvm_* variant の DRJIT_LIBLLVM_PATH 設定が mitsuba import より前に走ることを保証する． """sm-render コマンドラインエントリポイント． YAML シーン定義を読み込み，Mitsuba 3 でレンダリングを実行して PNG と luminance.npy を出力ディレクトリに保存する． + +処理本体は `render/pipeline.py` の `render_yaml()` にある（スイープ等から +プロセス内で直接呼べるようにするため）．本モジュールは引数のパースと， +例外 → 「エラー: …」表示 + 終了コード 1 への変換だけを担う． """ from __future__ import annotations @@ -16,10 +21,6 @@ import sys from pathlib import Path -# Step 6 の evaluate CLI も同じファイルを読むため，名前を定数として一元管理する -_PNG_FILENAME = "render.png" -_NPY_FILENAME = "luminance.npy" - def _parse_args(argv: list[str] | None = None) -> argparse.Namespace: """コマンドライン引数をパースして返す． @@ -67,13 +68,8 @@ def main(argv: list[str] | None = None) -> int: """sm-render コマンドのメインエントリポイント． - 処理フロー: - 1. コマンドライン引数をパースする． - 2. YAML シーン定義を load_scene_yaml で読み込み，SceneConfig を得る． - 3. 出力ディレクトリを作成する（親ごと，既存なら何もしない）． - 4. render_scene で Mitsuba レンダリングを実行し，HDR 配列を得る． - 5. save_png で render.png を保存する． - 6. save_luminance で luminance.npy を保存する． + 引数をパースして `render_yaml()` を呼び，段階別のエラーを終了コードに変換する． + レンダそのものの処理フローは `render/pipeline.py` の `render_yaml()` を参照． Args: argv: 引数リスト．None のとき sys.argv[1:] を使用． @@ -85,64 +81,17 @@ # ここでは ValueError / FileNotFoundError 等の実行時エラーのみ捕捉する args = _parse_args(argv) - scene_path: Path = args.scene.resolve() - output_dir: Path = args.output.resolve() + # mitsuba 依存のモジュールを読み込む前に引数をパースする（--help を軽量に保つ） + from sm_tias_lightsim.render.pipeline import RenderStepError, render_yaml try: - from sm_tias_lightsim.scene.schema import load_scene_yaml - - config = load_scene_yaml(scene_path) - except FileNotFoundError as exc: - print(f"エラー: シーン YAML ファイルが見つかりません — {exc}", file=sys.stderr) - return 1 - except Exception as exc: # noqa: BLE001 - print(f"エラー: シーン YAML の読み込みに失敗しました — {exc}", file=sys.stderr) - return 1 - - # 同じ YAML を手元（llvm）とサーバ（cuda）で使い回すため，variant だけ実行側から差し替える - # 優先順位: --variant > 環境変数 SM_RENDER_VARIANT > YAML の render.variant - from sm_tias_lightsim.runtime.overrides import resolve_render_variant - - variant = resolve_render_variant(config.render.variant, args.variant) - if variant != config.render.variant: - print(f"variant を上書き: {config.render.variant} → {variant}") - config = config.model_copy( - update={"render": config.render.model_copy(update={"variant": variant})} + png_path, npy_path = render_yaml( + args.scene.resolve(), + args.output.resolve(), + args.variant, ) - - try: - output_dir.mkdir(parents=True, exist_ok=True) - except OSError as exc: - print(f"エラー: 出力ディレクトリを作成できません — {exc}", file=sys.stderr) - return 1 - - # project_root を scene YAML の親ディレクトリとすることで， - # YAML 内の STL 相対パスが YAML ファイルを基点に解決される - project_root = scene_path.parent - try: - # llvm_* variant の場合，mitsuba の import チェーン（sensors.py / lights.py）が - # 走る前に DRJIT_LIBLLVM_PATH を設定しておく必要がある． - # runner.py の import（build_scene → sensors.py / lights.py → mitsuba）より先に呼ぶ． - if config.render.variant.startswith("llvm_"): - from sm_tias_lightsim.runtime.llvm_init import ensure_drjit_libllvm_path - - ensure_drjit_libllvm_path() - from sm_tias_lightsim.render.runner import render_scene - - hdr_image = render_scene(config, project_root=project_root) - except Exception as exc: # noqa: BLE001 - print(f"エラー: レンダリングに失敗しました — {exc}", file=sys.stderr) - return 1 - - png_path = output_dir / _PNG_FILENAME - npy_path = output_dir / _NPY_FILENAME - try: - from sm_tias_lightsim.render.output import save_luminance, save_png - - save_png(hdr_image, png_path) - save_luminance(hdr_image, npy_path) - except Exception as exc: # noqa: BLE001 - print(f"エラー: 出力ファイルの保存に失敗しました — {exc}", file=sys.stderr) + except RenderStepError as exc: + print(f"エラー: {exc}", file=sys.stderr) return 1 print(f"レンダリング完了: {png_path}") diff --git a/src/sm_tias_lightsim/evaluation/pipeline.py b/src/sm_tias_lightsim/evaluation/pipeline.py new file mode 100644 index 0000000..9ee3f52 --- /dev/null +++ b/src/sm_tias_lightsim/evaluation/pipeline.py @@ -0,0 +1,631 @@ +"""評価パイプラインの本体（sm-evaluate の処理を関数として提供する）． + +`evaluate_image()` が「画像読込 → 評価設定読込 → ROI 抽出 → 正規化 / OETF → +各種指標 → metrics.csv 保存 → 図保存」の一連の処理を行う．CLI（`cli/evaluate.py`）は +引数のパースと例外→終了コードの変換だけを担い，本モジュールを呼ぶ． + +スイープスクリプトのようにプロセス内から繰り返し評価する用途では +`save_figures=False` を渡すと figures/ の描画（1 条件あたり約 1.8 秒）を丸ごと省ける． + +matplotlib のバックエンドは Agg を強制する（headless 環境対応）． +evaluate は Mitsuba を使わないため，render 側と異なり遅延 import は不要． +ただし Agg バックエンドの設定順序を保つため，matplotlib.use() は pyplot import より前に置く． +""" + +from __future__ import annotations + +from pathlib import Path + +import matplotlib + +# headless 環境（サーバ・CI 等）でも動くよう Agg バックエンドを強制する． +# visualize.py より先に設定することで，pyplot import 時に Agg が有効になる． +matplotlib.use("Agg") + +import numpy # noqa: E402 +import pandas # noqa: E402 +from PIL import Image # noqa: E402 + +from sm_tias_lightsim.evaluation.normalize import apply_mean_anchor # noqa: E402 +from sm_tias_lightsim.evaluation.oetf import apply_srgb_oetf, solve_exposure # noqa: E402 +from sm_tias_lightsim.evaluation.radial import ( # noqa: E402 + calc_radial_min_max_ratio, + calc_radial_profile, +) +from sm_tias_lightsim.evaluation.roi import extract_roi # noqa: E402 +from sm_tias_lightsim.evaluation.shape_metrics import ( # noqa: E402 + calc_angular_radial_decomposition, + calc_intensity_centroid, + calc_intensity_pca, + calc_ncc, + calc_ssim, +) +from sm_tias_lightsim.evaluation.visualize import ( # noqa: E402 + save_histogram, + save_luminance_map, + save_radial_profile, + save_zone_map, +) +from sm_tias_lightsim.evaluation.zones import ( # noqa: E402 + calc_zone_stats, + classify_zones, + compute_normalized_distance_map, +) +from sm_tias_lightsim.scene.schema import ( # noqa: E402 + EvaluationConfig, + MetricsConfig, + NormalizationConfig, + OetfConfig, + ZonesConfig, +) +from sm_tias_lightsim.uniformity import calc_uniformity, to_grayscale # noqa: E402 + +# 出力ファイル名を一元管理する（変更時の影響範囲を最小化するため） +_METRICS_FILENAME = "metrics.csv" +_FIGURES_DIRNAME = "figures" +_LUMINANCE_MAP_FILENAME = "luminance_map.png" +_HISTOGRAM_FILENAME = "histogram.png" +_RADIAL_PROFILE_FILENAME = "radial_profile.png" +_ZONE_MAP_FILENAME = "zone_map.png" + + +class EvaluationStepError(RuntimeError): + """評価パイプラインのどの段階で失敗したかを文言として保持する例外． + + CLI は `print(f"エラー: {exc}")` するだけで従来どおりの段階別メッセージを出せる． + error-handling.md の「独自例外は最小限」に沿って階層は作らず，この 1 クラスのみとする + （段階名を CLI 側で分岐せずにメッセージを保つための最小限の器）． + """ + + +def _apply_oetf_to_luminance( + roi_luminance: numpy.ndarray, + oetf_cfg: OetfConfig, +) -> numpy.ndarray: + """OETF 順方向適用：ROI linear luminance を sRGB 符号化空間へ変換する． + + 無彩色シーン（シミュ LED・材料が RGB 等価）では luminance への直接 OETF 適用は + 「チャネル別 sRGB OETF → Rec.709 luma 加重和」と厳密一致する（SPEC_03 参照）． + + 露出ソルバは ROI 上で解く（target_mean は実測 ROI mean のため ROI 上で収束させる）． + + Args: + roi_luminance: ROI 切り出し済みの linear 輝度配列（H' x W', float64）． + oetf_cfg: OetfConfig（enabled=True が前提）． + + Returns: + OETF 適用後の ROI 輝度配列（H' x W', float64, quantize_uint8=True なら [0,255]）． + + Raises: + ValueError: oetf_cfg.exposure が None の場合（anchor_mean モードで必須）． + """ + if oetf_cfg.exposure is None: + msg = "oetf.enabled=True かつ anchor_mean モードでは exposure 設定が必須です．" + raise ValueError(msg) + + roi_lin = roi_luminance.astype(numpy.float64) + target_mean = oetf_cfg.exposure.target_mean + # ROI linear 値で露出スケールを決定（実測 mean が ROI mean のため ROI 上で解く） + scale, clipped = solve_exposure(roi_lin, target_mean) + clip_pct = clipped * 100.0 + print(f" OETF exposure : scale={scale:.6g} clip={clip_pct:.1f}%") + + # ROI に対してのみ OETF を適用する（指標算出は ROI 上で行うため） + return apply_srgb_oetf(roi_lin, scale, quantize_uint8=oetf_cfg.quantize_uint8) + + +def _apply_normalization_to_luminance( + roi_luminance: numpy.ndarray, + norm_cfg: NormalizationConfig, +) -> numpy.ndarray: + """linear mean anchor 正規化を ROI 輝度配列に適用する． + + ROI mean を norm_cfg.target_mean に揃える線形ゲインを適用する． + 形状指標（CoV / 勾配 / C/P 比 / max-min 比）は線形ゲインに不変なので， + 絶対値スケールだけが変わる． + + Args: + roi_luminance: ROI 切り出し済みの linear 輝度配列（H' x W', float64）． + norm_cfg: NormalizationConfig（enabled=True が前提）． + + Returns: + 正規化後の ROI 輝度配列（H' x W', float64）． + + Raises: + ValueError: ROI mean が 0 の場合（apply_mean_anchor から伝播）． + """ + roi_lin = roi_luminance.astype(numpy.float64) + target_mean = norm_cfg.target_mean + normalized, alpha = apply_mean_anchor(roi_lin, target_mean) + print(f" normalization : alpha={alpha:.6g} target_mean={target_mean}") + return normalized + + +def _load_luminance(image_path: Path) -> numpy.ndarray: + """画像ファイルを読み込み，輝度配列（2D float64）を返す． + + 拡張子が .npy の場合は numpy.load で読み込む． + それ以外（PNG 等）は PIL で RGB 画像として読み込み，to_grayscale で変換する． + .npy が 3D（H, W, 3）だった場合も to_grayscale を通す． + + Args: + image_path: 読み込む画像ファイルのパス． + + Returns: + 輝度配列（H x W, float64）． + + Raises: + FileNotFoundError: ファイルが存在しない場合． + ValueError: ファイルを輝度配列に変換できない場合． + """ + if not image_path.exists(): + msg = f"画像ファイルが見つかりません: {image_path}" + raise FileNotFoundError(msg) + + if image_path.suffix.lower() == ".npy": + data = numpy.load(image_path) + if data.ndim == 3: + return to_grayscale(data.astype(numpy.float64)) + if data.ndim == 2: + return data.astype(numpy.float64) + msg = f".npy の次元数が不正です（ndim={data.ndim}）．2D または 3D を期待します．" + raise ValueError(msg) + + pil_image = Image.open(image_path).convert("RGB") + return to_grayscale(numpy.array(pil_image, dtype=numpy.float64)) + + +def _load_eval_config(config_path: Path) -> EvaluationConfig: + """YAML を読み込んで EvaluationConfig を返す． + + SceneConfig の evaluation サブツリー，または EvaluationConfig 互換の + トップレベル YAML の両方をサポートする． + + Args: + config_path: YAML ファイルのパス． + + Returns: + EvaluationConfig インスタンス． + + Raises: + FileNotFoundError: ファイルが存在しない場合． + pydantic.ValidationError: スキーマ違反． + yaml.YAMLError: YAML 構文エラー． + """ + import yaml + + if not config_path.exists(): + msg = f"設定 YAML ファイルが見つかりません: {config_path}" + raise FileNotFoundError(msg) + + with config_path.open(encoding="utf-8") as f: + data = yaml.safe_load(f) + + # "evaluation" キーがあれば SceneConfig 埋め込みとみなしてサブツリーを取得する + if isinstance(data, dict) and "evaluation" in data: + eval_data = data["evaluation"] + else: + eval_data = data + + return EvaluationConfig.model_validate(eval_data) + + +def _run_zone_analysis( + roi_luminance: numpy.ndarray, + zones_config: ZonesConfig | None, +) -> tuple[numpy.ndarray, dict] | tuple[None, dict]: + """ゾーン分析を実行し，zone_map とゾーン統計を返す． + + Args: + roi_luminance: ROI 輝度配列（H x W, float64）． + zones_config: ZonesConfig インスタンス（enabled / boundaries を持つ）． + None の場合は無効扱い． + + Returns: + (zone_map, zone_stats) のタプル． + zones が無効なら (None, {}) を返す． + + Raises: + ValueError: ゾーン分析の計算に失敗した場合． + """ + if zones_config is None or not zones_config.enabled: + return None, {} + + h, w = roi_luminance.shape + distance_map = compute_normalized_distance_map(h, w) + zone_map = classify_zones(distance_map, zones_config.boundaries) + zone_stats = calc_zone_stats(roi_luminance, zone_map) + return zone_map, zone_stats + + +def _run_radial_profile( + roi_luminance: numpy.ndarray, + metrics_config: MetricsConfig | None, +) -> numpy.ndarray | None: + """放射状プロファイルを算出して返す． + + Args: + roi_luminance: ROI 輝度配列（H x W, float64）． + metrics_config: MetricsConfig インスタンス（radial_profile フラグを持つ）． + None の場合は無効扱い． + + Returns: + 放射状プロファイル配列（(n_bins, 2)），無効時は None． + """ + if metrics_config is None or not metrics_config.radial_profile: + return None + return calc_radial_profile(roi_luminance) + + +def _run_radial_min_max( + roi_luminance: numpy.ndarray, + metrics_config: MetricsConfig | None, + radial_profile: numpy.ndarray | None, +) -> dict: + """動径 min/max 比を算出して返す． + + radial_min_max が True のとき，radial_profile が既に算出されていればそれを再利用し， + 未算出（None）の場合はその場で calc_radial_profile を呼び出して算出する． + 重複算出を避けるため，radial_profile の結果を優先的に使用する． + + Args: + roi_luminance: ROI 輝度配列（H x W, float64）． + metrics_config: MetricsConfig インスタンス（radial_min_max フラグを持つ）． + None の場合は無効扱い． + radial_profile: _run_radial_profile で算出済みのプロファイル配列． + None の場合はここで算出する． + + Returns: + 動径 min/max 比指標の辞書（radial_min_max_ratio / radial_min_distance / + radial_min / radial_max の 4 キー）．無効時は空辞書． + """ + if metrics_config is None or not metrics_config.radial_min_max: + return {} + + # radial_profile が未算出の場合はここで算出する（重複算出を避けるため既存結果を優先） + profile = radial_profile if radial_profile is not None else calc_radial_profile(roi_luminance) + return calc_radial_min_max_ratio(profile) + + +def _load_reference_luminance(reference_path: Path) -> numpy.ndarray: + """形状指標用の参照画像を読み込み，輝度配列（2D float64）を返す． + + SSIM / NCC の参照画像として使用する．_load_luminance と同一の変換ロジックを用い， + .npy 形式（2D/3D）と PNG 形式の両方に対応する． + 呼び出し元（_run_shape_metrics）で shape を確認してから calc_ssim / calc_ncc に渡すこと． + + Args: + reference_path: 参照画像ファイルのパス（.npy または .png）． + + Returns: + 輝度配列（H x W, float64）． + + Raises: + FileNotFoundError: ファイルが存在しない場合． + ValueError: ファイルを輝度配列に変換できない場合． + """ + return _load_luminance(reference_path) + + +def _run_shape_metrics( + roi_luminance: numpy.ndarray, + metrics_config: MetricsConfig | None, + config_dir: Path, +) -> dict: + """形状指標（輝度重心 / 角度-放射分解 / PCA 主軸 / SSIM / NCC）を算出する． + + metrics.shape_metrics=True のとき，以下の指標を算出して辞書として返す: + - 輝度重心 5 値 (calc_intensity_centroid) + - 角度-放射分解 3 値 (calc_angular_radial_decomposition) + - PCA 主軸 4 値 (calc_intensity_pca) + - SSIM / NCC（shape_reference が指定されたときのみ，それ以外は float("nan")） + + shape_reference パスは config YAML のディレクトリを基準に解決する． + 絶対パスの場合はそのまま使用する． + + Args: + roi_luminance: ROI 輝度配列（H x W, float64）． + metrics_config: MetricsConfig インスタンス（shape_metrics フラグを持つ）． + None の場合は無効扱い． + config_dir: 評価 YAML ファイルのディレクトリ（shape_reference の相対パス解決用）． + + Returns: + 形状指標の辞書（12 キー）．無効時は空辞書． + SSIM / NCC が計算不能な場合は float("nan") を格納する． + """ + if metrics_config is None or not metrics_config.shape_metrics: + return {} + + result: dict = {} + + # 輝度重心 + centroid = calc_intensity_centroid(roi_luminance) + result["intensity_centroid_x_px"] = centroid["centroid_x_px"] + result["intensity_centroid_y_px"] = centroid["centroid_y_px"] + result["intensity_centroid_x_norm"] = centroid["centroid_x_norm"] + result["intensity_centroid_y_norm"] = centroid["centroid_y_norm"] + result["intensity_centroid_offset_norm"] = centroid["centroid_offset_norm"] + + # 角度-放射分解 + ard = calc_angular_radial_decomposition(roi_luminance) + result["angular_variance"] = ard["angular_variance"] + result["radial_variance"] = ard["radial_variance"] + result["angular_ratio"] = ard["angular_ratio"] + + # PCA 主軸 + pca = calc_intensity_pca(roi_luminance) + result["pca_primary_axis_angle_deg"] = pca["primary_axis_angle_deg"] + result["pca_eccentricity"] = pca["eccentricity"] + + # SSIM / NCC（shape_reference が指定されている場合のみ） + reference_path_str = metrics_config.shape_reference + if reference_path_str is not None: + ref_path = Path(reference_path_str) + if not ref_path.is_absolute(): + ref_path = config_dir / ref_path + + try: + ref_luminance = _load_reference_luminance(ref_path) + # shape が一致するか確認（一致しない場合は警告してスキップ） + if ref_luminance.shape != roi_luminance.shape: + print( + f" 警告: shape_reference の shape {ref_luminance.shape} が " + f"ROI の shape {roi_luminance.shape} と異なるため SSIM/NCC をスキップします．" + ) + result["ssim"] = float("nan") + result["ncc"] = float("nan") + else: + result["ssim"] = calc_ssim(roi_luminance, ref_luminance) + result["ncc"] = calc_ncc(roi_luminance, ref_luminance) + except FileNotFoundError: + print( + f" 警告: shape_reference ファイルが見つかりません: {ref_path}．" + "SSIM/NCC をスキップします．" + ) + result["ssim"] = float("nan") + result["ncc"] = float("nan") + except Exception as exc: # noqa: BLE001 + print( + f" 警告: shape_reference の読み込みに失敗しました — {exc}．" + "SSIM/NCC をスキップします．" + ) + result["ssim"] = float("nan") + result["ncc"] = float("nan") + else: + result["ssim"] = float("nan") + result["ncc"] = float("nan") + + return result + + +def _save_metrics_csv(metrics: dict, csv_path: Path) -> None: + """均一性指標を CSV ファイルとして保存する． + + Args: + metrics: 指標名 → 値の辞書． + csv_path: 保存先 CSV パス． + """ + df = pandas.DataFrame([metrics]) + df.to_csv(csv_path, index=False) + + +def _save_figures( + roi_luminance: numpy.ndarray, + figures_dir: Path, + radial_profile: numpy.ndarray | None, + zone_map: numpy.ndarray | None, +) -> None: + """可視化 PNG を figures/ ディレクトリに保存する． + + Args: + roi_luminance: ROI 輝度配列（H x W, float64）． + figures_dir: 保存先ディレクトリ（事前に mkdir 済みであること）． + radial_profile: 放射状プロファイル配列，None の場合は保存しない． + zone_map: ゾーンマップ配列，None の場合は保存しない． + """ + save_luminance_map(roi_luminance, figures_dir / _LUMINANCE_MAP_FILENAME) + save_histogram(roi_luminance, figures_dir / _HISTOGRAM_FILENAME) + + if radial_profile is not None: + save_radial_profile(radial_profile, figures_dir / _RADIAL_PROFILE_FILENAME) + + if zone_map is not None: + save_zone_map(zone_map, figures_dir / _ZONE_MAP_FILENAME, luminance=roi_luminance) + + +def _print_summary(metrics: dict, csv_path: Path, figures_dir: Path | None) -> None: + """評価サマリーを標準出力に表示する． + + Args: + metrics: 指標名 → 値の辞書． + csv_path: metrics.csv の保存先パス（表示用）． + figures_dir: figures/ ディレクトリのパス（表示用）． + 図を保存しない場合は None（figures/ の行を出力しない）． + """ + print("=== 均一性評価サマリー ===") + print(f" CoV : {metrics['cov']:.6f}") + print(f" std : {metrics['std']:.4f}") + print(f" max/min ratio : {metrics['max_min_ratio']:.4f}") + print(f" mean : {metrics['mean']:.4f}") + print(f" max : {metrics['max']:.4f}") + print(f" min : {metrics['min']:.4f}") + + if "center_periphery_ratio" in metrics: + print(f" C/P ratio : {metrics['center_periphery_ratio']:.4f}") + print(f" gradient [%] : {metrics['gradient_magnitude']:.2f}") + + if "radial_min_max_ratio" in metrics: + print(f" radial min/max : {metrics['radial_min_max_ratio']:.4f}") + print(f" radial min dist: {metrics['radial_min_distance']:.4f}") + + if "intensity_centroid_offset_norm" in metrics: + print(f" centroid offset: {metrics['intensity_centroid_offset_norm']:.4f} (norm)") + print( + f" centroid (norm): x={metrics['intensity_centroid_x_norm']:.4f}" + f" y={metrics['intensity_centroid_y_norm']:.4f}" + ) + ang_ratio = metrics.get("angular_ratio", float("nan")) + if ang_ratio != ang_ratio: # NaN check + print(" angular_ratio : NaN (uniform image)") + else: + print(f" angular_ratio : {ang_ratio:.4f} (0=concentric, 1=directional)") + print(f" pca axis [deg] : {metrics['pca_primary_axis_angle_deg']:.1f}") + print(f" pca eccentricit: {metrics['pca_eccentricity']:.4f}") + ssim_val = metrics.get("ssim", float("nan")) + ncc_val = metrics.get("ncc", float("nan")) + if ssim_val == ssim_val: # not NaN + print(f" SSIM : {ssim_val:.4f}") + print(f" NCC : {ncc_val:.4f}") + + print(f"\nmetrics.csv : {csv_path}") + if figures_dir is not None: + print(f"figures/ : {figures_dir}") + + + +def evaluate_image( + image_path: Path, + config_path: Path, + output_dir: Path, + save_figures: bool = True, +) -> dict: + """画像 1 枚を評価し，metrics.csv（と図）を出力して指標の辞書を返す． + + 処理フロー: + 1. 画像ファイルを読み込み，輝度配列を得る． + 2. YAML 設定ファイルを読み込み，EvaluationConfig を得る． + 3. 出力ディレクトリ（save_figures のとき figures/ も）を作成する． + 4. extract_roi で ROI を切り出す． + 5. normalization.enabled のとき linear mean anchor 正規化を適用する． + 6. oetf.enabled のとき sRGB OETF を適用する． + 7. calc_uniformity で均一性指標を算出する． + 8. zones.enabled のときゾーン分析を実行する． + 9. metrics.radial_profile のとき放射状プロファイルを算出する． + 10. metrics.radial_min_max のとき動径 min/max 比を算出する（radial_profile 再利用）． + 11. metrics.shape_metrics のとき形状指標を算出する． + 12. metrics.csv を保存する． + 13. save_figures のとき可視化 PNG を figures/ に保存する． + 14. 算出指標のサマリーを標準出力に表示する． + + Args: + image_path: 評価対象の画像ファイル（PNG または .npy）． + config_path: 評価設定 YAML（SceneConfig または EvaluationConfig 互換）． + output_dir: 出力先ディレクトリ（存在しなければ親ごと作成する）． + save_figures: False のとき figures/ を作らず図も保存しない + （スイープ等で図が不要なときに描画時間を省くため）． + + Returns: + 指標名 → 値の辞書（metrics.csv に書いた内容と同じ）． + + Raises: + EvaluationStepError: いずれかの段階で失敗した場合（失敗した段階を示す文言を持つ）． + """ + image_path = Path(image_path).resolve() + config_path = Path(config_path).resolve() + output_dir = Path(output_dir).resolve() + + try: + luminance = _load_luminance(image_path) + except FileNotFoundError as exc: + raise EvaluationStepError(str(exc)) from exc + except Exception as exc: + msg = f"画像の読み込みに失敗しました — {exc}" + raise EvaluationStepError(msg) from exc + + try: + eval_config = _load_eval_config(config_path) + except FileNotFoundError as exc: + raise EvaluationStepError(str(exc)) from exc + except Exception as exc: + msg = f"設定 YAML の読み込みに失敗しました — {exc}" + raise EvaluationStepError(msg) from exc + + figures_dir = output_dir / _FIGURES_DIRNAME if save_figures else None + try: + output_dir.mkdir(parents=True, exist_ok=True) + if figures_dir is not None: + figures_dir.mkdir(parents=True, exist_ok=True) + except OSError as exc: + msg = f"出力ディレクトリを作成できません — {exc}" + raise EvaluationStepError(msg) from exc + + try: + roi_luminance = extract_roi(luminance, eval_config.roi) + except (ValueError, NotImplementedError) as exc: + msg = f"ROI 抽出に失敗しました — {exc}" + raise EvaluationStepError(msg) from exc + + # linear mean anchor 正規化（enabled のとき ROI 抽出後・OETF 前に挟む）． + # OETF との併用禁止は EvaluationConfig の model_validator で保証されている． + norm_cfg = eval_config.normalization + if norm_cfg is not None and norm_cfg.enabled: + try: + roi_luminance = _apply_normalization_to_luminance(roi_luminance, norm_cfg) + except ValueError as exc: + msg = f"linear 正規化に失敗しました — {exc}" + raise EvaluationStepError(msg) from exc + + # OETF 順方向適用（enabled のとき ROI 抽出後・指標算出前に挟む）． + # eval_real.yaml（実測 PNG）は oetf を設定しないこと： + # 実測 PNG は既に sRGB 符号化済みのため，ここで OETF を掛けると二重適用になる． + oetf_cfg = eval_config.oetf + if oetf_cfg is not None and oetf_cfg.enabled: + try: + roi_luminance = _apply_oetf_to_luminance(roi_luminance, oetf_cfg) + except ValueError as exc: + msg = f"OETF 適用に失敗しました — {exc}" + raise EvaluationStepError(msg) from exc + + try: + uniformity_stats = calc_uniformity(roi_luminance) + except ValueError as exc: + msg = f"均一性指標の算出に失敗しました — {exc}" + raise EvaluationStepError(msg) from exc + + metrics: dict = dict(uniformity_stats) + + try: + zone_map, zone_stats = _run_zone_analysis(roi_luminance, eval_config.zones) + metrics.update(zone_stats) + except ValueError as exc: + msg = f"ゾーン分析に失敗しました — {exc}" + raise EvaluationStepError(msg) from exc + + try: + radial_profile = _run_radial_profile(roi_luminance, eval_config.metrics) + except Exception as exc: + msg = f"放射状プロファイルの算出に失敗しました — {exc}" + raise EvaluationStepError(msg) from exc + + try: + radial_min_max_stats = _run_radial_min_max( + roi_luminance, eval_config.metrics, radial_profile + ) + metrics.update(radial_min_max_stats) + except Exception as exc: + msg = f"動径 min/max 比の算出に失敗しました — {exc}" + raise EvaluationStepError(msg) from exc + + try: + shape_stats = _run_shape_metrics(roi_luminance, eval_config.metrics, config_path.parent) + metrics.update(shape_stats) + except Exception as exc: + msg = f"形状指標の算出に失敗しました — {exc}" + raise EvaluationStepError(msg) from exc + + csv_path = output_dir / _METRICS_FILENAME + try: + _save_metrics_csv(metrics, csv_path) + except Exception as exc: + msg = f"metrics.csv の保存に失敗しました — {exc}" + raise EvaluationStepError(msg) from exc + + if figures_dir is not None: + try: + _save_figures(roi_luminance, figures_dir, radial_profile, zone_map) + except Exception as exc: + msg = f"可視化 PNG の保存に失敗しました — {exc}" + raise EvaluationStepError(msg) from exc + + _print_summary(metrics, csv_path, figures_dir) + return metrics diff --git a/src/sm_tias_lightsim/render/pipeline.py b/src/sm_tias_lightsim/render/pipeline.py new file mode 100644 index 0000000..2c34542 --- /dev/null +++ b/src/sm_tias_lightsim/render/pipeline.py @@ -0,0 +1,118 @@ +"""レンダパイプラインの本体（sm-render の処理を関数として提供する）． + +`render_yaml()` が「シーン YAML 読込 → variant 解決 → レンダ → render.png / +luminance.npy 保存」を行う．CLI（`cli/render.py`）は引数のパースと， +例外 → 「エラー: …」表示 + 終了コード 1 への変換だけを担う． + +スイープスクリプトのように何十条件も回す用途では，サブプロセス（`uv run sm-render`）の +起動と `import mitsuba` が 1 条件あたり 2 秒前後かかるため，本関数をプロセス内で +繰り返し呼ぶことで丸ごと省ける． + +mitsuba に触れる import（runner / builder / schema）はすべて関数内で行う． +モジュールレベルで import すると，llvm_* variant 用の DRJIT_LIBLLVM_PATH 設定より前に +LLVM 初期化が走ってしまうため（cli/render.py と同じ遅延 import 方針）． +""" + +from __future__ import annotations + +from pathlib import Path + +# Step 6 の evaluate CLI も同じファイルを読むため，名前を定数として一元管理する +_PNG_FILENAME = "render.png" +_NPY_FILENAME = "luminance.npy" + + +class RenderStepError(RuntimeError): + """レンダパイプラインのどの段階で失敗したかを文言として保持する例外． + + CLI は `print(f"エラー: {exc}")` するだけで従来どおりの段階別メッセージを出せる． + error-handling.md の「独自例外は最小限」に沿って階層は作らず，この 1 クラスのみとする． + """ + + +def render_yaml( + scene_yaml: Path, + output_dir: Path, + variant: str | None = None, +) -> tuple[Path, Path]: + """シーン YAML をレンダリングし，render.png と luminance.npy を出力する． + + 処理フロー: + 1. YAML シーン定義を load_scene_yaml で読み込み，SceneConfig を得る． + 2. variant を解決する（引数 > 環境変数 SM_RENDER_VARIANT > YAML の render.variant）． + 3. 出力ディレクトリを作成する（親ごと，既存なら何もしない）． + 4. llvm_* のとき DRJIT_LIBLLVM_PATH を設定してから render_scene を実行する． + 5. render.png と luminance.npy を保存する． + + Args: + scene_yaml: シーン定義 YAML ファイルのパス． + output_dir: 出力先ディレクトリ（存在しない場合は親ごと作成する）． + variant: Mitsuba variant の明示指定．None なら環境変数 / YAML から解決する． + + Returns: + (render.png のパス, luminance.npy のパス)． + + Raises: + RenderStepError: いずれかの段階で失敗した場合（失敗した段階を示す文言を持つ）． + """ + scene_path = Path(scene_yaml).resolve() + out_dir = Path(output_dir).resolve() + + try: + from sm_tias_lightsim.scene.schema import load_scene_yaml + + config = load_scene_yaml(scene_path) + except FileNotFoundError as exc: + msg = f"シーン YAML ファイルが見つかりません — {exc}" + raise RenderStepError(msg) from exc + except Exception as exc: + msg = f"シーン YAML の読み込みに失敗しました — {exc}" + raise RenderStepError(msg) from exc + + # 同じ YAML を手元（llvm）とサーバ（cuda）で使い回すため，variant だけ実行側から差し替える + # 優先順位: 引数 variant > 環境変数 SM_RENDER_VARIANT > YAML の render.variant + from sm_tias_lightsim.runtime.overrides import resolve_render_variant + + resolved = resolve_render_variant(config.render.variant, variant) + if resolved != config.render.variant: + print(f"variant を上書き: {config.render.variant} → {resolved}") + config = config.model_copy( + update={"render": config.render.model_copy(update={"variant": resolved})} + ) + + try: + out_dir.mkdir(parents=True, exist_ok=True) + except OSError as exc: + msg = f"出力ディレクトリを作成できません — {exc}" + raise RenderStepError(msg) from exc + + # project_root を scene YAML の親ディレクトリとすることで， + # YAML 内の STL 相対パスが YAML ファイルを基点に解決される + project_root = scene_path.parent + try: + # llvm_* variant の場合，mitsuba の import チェーン（sensors.py / lights.py）が + # 走る前に DRJIT_LIBLLVM_PATH を設定しておく必要がある． + # runner.py の import（build_scene → sensors.py / lights.py → mitsuba）より先に呼ぶ． + if config.render.variant.startswith("llvm_"): + from sm_tias_lightsim.runtime.llvm_init import ensure_drjit_libllvm_path + + ensure_drjit_libllvm_path() + from sm_tias_lightsim.render.runner import render_scene + + hdr_image = render_scene(config, project_root=project_root) + except Exception as exc: + msg = f"レンダリングに失敗しました — {exc}" + raise RenderStepError(msg) from exc + + png_path = out_dir / _PNG_FILENAME + npy_path = out_dir / _NPY_FILENAME + try: + from sm_tias_lightsim.render.output import save_luminance, save_png + + save_png(hdr_image, png_path) + save_luminance(hdr_image, npy_path) + except Exception as exc: + msg = f"出力ファイルの保存に失敗しました — {exc}" + raise RenderStepError(msg) from exc + + return png_path, npy_path diff --git a/src/sm_tias_lightsim/render/runner.py b/src/sm_tias_lightsim/render/runner.py index 49914fc..8eb362d 100644 --- a/src/sm_tias_lightsim/render/runner.py +++ b/src/sm_tias_lightsim/render/runner.py @@ -26,6 +26,70 @@ _current_variant: str | None = None +def ensure_variant(variant: str) -> None: + """Mitsuba の variant をプロセス内で 1 度だけ設定する（2 回目以降は何もしない）． + + llvm_* variant 要求時は mi.set_variant 前に DRJIT_LIBLLVM_PATH を自動探索して設定する． + mitsuba は関数内でインポートするため，モジュール import 時点では LLVM 初期化が走らない． + set_variant 成功後に variant 状態を INFO ログで出力する． + 既に別の variant が設定されている場合は警告ログを出して続行する + （Mitsuba は同一プロセスで variant を切り替えられないため）． + + レンダ前だけでなく，レイキャストによる検証（遮蔽チェック等）を目的の variant で + 走らせたいときにも先に呼べるよう，render_scene から切り出してある． + + Args: + variant: 設定したい Mitsuba variant（例: "llvm_ad_rgb" / "cuda_ad_rgb"）． + + Raises: + ImportError: JIT の初期化に失敗した場合（LLVM-C.dll 未検出等）． + """ + global _current_variant # noqa: PLW0603 + + if _current_variant is None: + # llvm_* variant の場合のみ DRJIT_LIBLLVM_PATH を自動探索する + # scalar_* / cuda_* は LLVM-C.dll 不要なので探索しない + # mitsuba の import 前に os.environ['DRJIT_LIBLLVM_PATH'] を設定することが必須． + # mitsuba を先に import すると LLVM 初期化が走り，LLVM-C.dll 未設定で ImportError になる． + if variant.startswith("llvm_"): + ensure_drjit_libllvm_path() + + # mitsuba は DRJIT_LIBLLVM_PATH 設定後にインポートする． + # モジュールレベルで import すると，ensure_drjit_libllvm_path() の前に LLVM 初期化が走り + # ImportError になるため，この位置に遅延 import する． + import mitsuba as mi + + # 共用サーバで LLVM バックエンドが全コアを占有しないよう，SM_DRJIT_THREADS があれば制限する + apply_drjit_thread_count() + + if _current_variant is None: + try: + mi.set_variant(variant) + except ImportError as e: + if variant.startswith("llvm_") and "LLVM" in str(e): + _logger.error( + "Mitsuba variant=%s の JIT 初期化に失敗しました．" + "LLVM-C.dll が見つかりません．" + "Windows: C:\\Program Files\\LLVM\\bin\\LLVM-C.dll を配置するか，" + "DRJIT_LIBLLVM_PATH 環境変数を設定してください．" + "詳細は docs/02_ENV/ENV_02_環境構築手順.md を参照．", + variant, + ) + raise + + _current_variant = variant + log_variant_status(_logger, variant, mi.variant()) + elif _current_variant != variant: + _logger.warning( + "Mitsuba variant is already set to %r; " + "requested variant %r will be ignored. " + "Restart the process to change the variant.", + _current_variant, + variant, + ) + # 同じ variant が既に設定されている場合は何もしない + + def render_scene(config: SceneConfig, project_root: Path | None = None) -> numpy.ndarray: """SceneConfig からシーンを構築してレンダリングし，HDR 画像（H, W, 3）を返す． @@ -43,57 +107,15 @@ Returns: HDR 画像 (H, W, 3) の float32 配列．値域は物理量（0〜∞ も許容）． """ - global _current_variant # noqa: PLW0603 + ensure_variant(config.render.variant) - desired_variant = config.render.variant - - if _current_variant is None: - # llvm_* variant の場合のみ DRJIT_LIBLLVM_PATH を自動探索する - # scalar_* / cuda_* は LLVM-C.dll 不要なので探索しない - # mitsuba の import 前に os.environ['DRJIT_LIBLLVM_PATH'] を設定することが必須． - # mitsuba を先に import すると LLVM 初期化が走り，LLVM-C.dll 未設定で ImportError になる． - if desired_variant.startswith("llvm_"): - ensure_drjit_libllvm_path() - - # mitsuba は DRJIT_LIBLLVM_PATH 設定後にインポートする． - # モジュールレベルで import すると，ensure_drjit_libllvm_path() の前に LLVM 初期化が走り - # ImportError になるため，この位置に遅延 import する． + # ensure_variant で DRJIT_LIBLLVM_PATH 設定済み・variant 設定済みのため，ここで import する import mitsuba as mi - # 共用サーバで LLVM バックエンドが全コアを占有しないよう，SM_DRJIT_THREADS があれば制限する - apply_drjit_thread_count() - - if _current_variant is None: - try: - mi.set_variant(desired_variant) - except ImportError as e: - if desired_variant.startswith("llvm_") and "LLVM" in str(e): - _logger.error( - "Mitsuba variant=%s の JIT 初期化に失敗しました．" - "LLVM-C.dll が見つかりません．" - "Windows: C:\\Program Files\\LLVM\\bin\\LLVM-C.dll を配置するか，" - "DRJIT_LIBLLVM_PATH 環境変数を設定してください．" - "詳細は docs/02_ENV/ENV_02_環境構築手順.md を参照．", - desired_variant, - ) - raise - - _current_variant = desired_variant - log_variant_status(_logger, desired_variant, mi.variant()) - elif _current_variant != desired_variant: - _logger.warning( - "Mitsuba variant is already set to %r; " - "requested variant %r will be ignored. " - "Restart the process to change the variant.", - _current_variant, - desired_variant, - ) - # 同じ variant が既に設定されている場合は何もしない - scene_dict, tempdir = build_scene(config, project_root=project_root) try: scene = mi.load_dict(scene_dict) - image_tensor = mi.render(scene, spp=config.render.spp) + image_tensor = mi.render(scene, spp=config.render.spp, seed=config.render.seed) arr = numpy.array(image_tensor) finally: shutil.rmtree(tempdir, ignore_errors=True) diff --git a/src/sm_tias_lightsim/scene/schema.py b/src/sm_tias_lightsim/scene/schema.py index 6eee4fa..ec9996f 100644 --- a/src/sm_tias_lightsim/scene/schema.py +++ b/src/sm_tias_lightsim/scene/schema.py @@ -98,6 +98,10 @@ ) max_depth: int = Field(default=8, ge=1, description="最大光線追跡深度") spp: int = Field(default=256, ge=1, description="ピクセルあたりのサンプル数") + seed: int = Field( + default=0, ge=0, + description="乱数シード（mi.render の seed）．再レンダでノイズの実現を変えるときに使う", + ) # ---------------------------------------------------------------------------