diff --git a/scripts/build_minmax_definition_comparison.py b/scripts/build_minmax_definition_comparison.py new file mode 100644 index 0000000..2ac2a3d --- /dev/null +++ b/scripts/build_minmax_definition_comparison.py @@ -0,0 +1,166 @@ +"""最小/最大比の「最大・最小の定義」による違いを軸上スイープ（TEST_15/16）で比較する. + +均一性の指標として使っている最小/最大比は，最大・最小をどう取るかで値が変わる． +本スクリプトは同じ luminance.npy に対して次の定義を計算し，サイズ依存の傾向と定義間の差を +表と図にする（TEST_18 の非対称な分布に進む前に，対称な分布で各定義の基準を押さえる）． + +定義: + pixel : 画素の min / max（評価パイプラインの max_min_ratio と同じ．MC ノイズに支配される） + ring : 20 輪帯平均の min / max（評価パイプラインの radial_min_max_ratio と同じ） + ring_inner : 輪帯平均のうち四隅だけの最外輪帯（d > 0.95）を除いた min / max + （scripts/build_radial_comparison.py の ring_min_max_ratio と同じ） + smooth : 平滑化（メディアン 5 → ガウス σ=4）後の画素の min / max + robust : 平滑化後の 1 / 99 パーセンタイル（方位に依らないロバスト定義．TEST_18 の + マップ表示と同じ「最大」） + +出力（output/_normalization_comparison/）: + minmax_definitions.csv : 条件 × 定義の値 + minmax_definitions.png : サイズ依存を定義ごとに重ね描き（球・立方体） + +実行方法: + uv run python scripts/build_minmax_definition_comparison.py +""" + +from __future__ import annotations + +import argparse +import csv +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 build_led_placement_maps import ROBUST_PERCENTILE, smoothed_roi # noqa: E402 +from build_normalization_comparison import OUT_DIR, conditions # noqa: E402 +from build_radial_comparison import MAX_DISTANCE # noqa: E402 + +from sm_tias_lightsim.evaluation.radial import calc_radial_profile # noqa: E402 +from sm_tias_lightsim.evaluation.roi_presets import SIM_ROI, extract_roi_xywh # noqa: E402 + +plt.rcParams["font.family"] = ["Yu Gothic", "MS Gothic", "DejaVu Sans"] + +DEFINITIONS = ("pixel", "ring", "ring_inner", "smooth", "robust") +DEF_LABEL = { + "pixel": "画素 min/max", + "ring": "輪帯平均 min/max（20 ビン）", + "ring_inner": "輪帯平均 min/max（最外ビン除外）", + "smooth": "平滑化後 画素 min/max", + "robust": f"平滑化後 {100 - ROBUST_PERCENTILE:g}/{ROBUST_PERCENTILE:g} %ile", +} +DEF_STYLE = { + "pixel": dict(marker="x", ls=":", color="gray"), + "ring": dict(marker="o", ls="-", color="tab:blue"), + "ring_inner": dict(marker="o", ls="--", color="tab:cyan"), + "smooth": dict(marker="s", ls="-.", color="tab:orange"), + "robust": dict(marker="D", ls="-", color="tab:red"), +} +LOW_PERCENTILE = 100.0 - ROBUST_PERCENTILE + + +def min_max_ratios(path: Path) -> dict[str, float]: + """1 条件の輝度マップについて全定義の最小/最大比を返す.""" + arr = np.load(path).astype(np.float64) + if arr.ndim == 3: + arr = 0.2126 * arr[..., 0] + 0.7152 * arr[..., 1] + 0.0722 * arr[..., 2] + roi = extract_roi_xywh(arr, SIM_ROI) + profile = calc_radial_profile(roi, n_bins=20) + centers, means = profile[:, 0], profile[:, 1] + valid = ~np.isnan(means) + inner = valid & (centers <= MAX_DISTANCE) + smooth = smoothed_roi(path) + return { + "pixel": float(roi.min() / roi.max()), + "ring": float(means[valid].min() / means[valid].max()), + "ring_inner": float(means[inner].min() / means[inner].max()), + "smooth": float(smooth.min() / smooth.max()), + "robust": float(np.percentile(smooth, LOW_PERCENTILE) + / np.percentile(smooth, ROBUST_PERCENTILE)), + } + + +def draw(rows: list[dict], output: Path) -> None: + """定義ごとの最小/最大比をサイズに対して描く（球・立方体で 1 枚ずつ + 定義間の差）.""" + fig, axes = plt.subplots(1, 3, figsize=(19, 5.5)) + for ax, shape, xlabel in ((axes[0], "sphere", "球 内半径 r [mm]"), + (axes[1], "cube", "立方体 内寸 L [mm]")): + sel = [r for r in rows if r["shape"] == shape] + sizes = [r["size_mm"] for r in sel] + for d in DEFINITIONS: + ax.plot(sizes, [r[d] for r in sel], label=DEF_LABEL[d], ms=5, lw=1.4, **DEF_STYLE[d]) + ax.set_xlabel(xlabel) + ax.set_ylabel("最小/最大比（1 が理想）") + ax.set_title("球 similar_radial（TEST_15）" if shape == "sphere" else "立方体（TEST_16）", + fontsize=11) + ax.grid(alpha=0.3) + ax.legend(fontsize=8) + ax = axes[2] + for shape, marker in (("sphere", "o"), ("cube", "s")): + sel = [r for r in rows if r["shape"] == shape] + areas = [r["area"] for r in sel] + label = "球" if shape == "sphere" else "立方体" + ax.plot(areas, [r["ring_inner"] - r["robust"] for r in sel], marker=marker, ls="-", + color="tab:red", label=f"{label}: 輪帯（最外除外）− ロバスト") + ax.plot(areas, [r["ring"] - r["ring_inner"] for r in sel], marker=marker, ls="--", + color="tab:blue", label=f"{label}: 輪帯（全） − 輪帯（最外除外）") + ax.plot(areas, [r["smooth"] - r["robust"] for r in sel], marker=marker, ls=":", + color="tab:orange", label=f"{label}: 平滑化 min/max − ロバスト") + ax.axhline(0.0, color="gray", lw=0.8) + ax.set_xscale("log") + ax.set_xlabel("キャビティ内面積 (mm²・対数)") + ax.set_ylabel("定義間の差") + ax.set_title("定義間の差（画素 min/max は桁が違うので除外）", fontsize=11) + ax.grid(alpha=0.3, which="both") + ax.legend(fontsize=7) + fig.suptitle("最小/最大比の定義による違い — 軸上 1 灯・板向き・バッフル付き（対称な分布）", + fontsize=13) + fig.tight_layout() + output.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(output, dpi=110) + plt.close(fig) + print(f"出力: {output}") + + +def main() -> int: + """全定義の最小/最大比を計算して表と図を出す.""" + parser = argparse.ArgumentParser(description="最小/最大比の定義による違いを比較") + parser.parse_args() + from sweep_cube_size import cavity_geometry as cube_geo + from sweep_sphere_radius import cavity_geometry as sphere_geo + + rows = [] + for shape, items in conditions().items(): + for label, path in items: + size = float(label.split("=")[1]) + area = (sphere_geo(size, "similar_radial") if shape == "sphere" + else cube_geo(size))["cavity_area_mm2"] + row = {"shape": shape, "size_mm": size, "area": area} + row.update(min_max_ratios(path)) + rows.append(row) + if not rows: + print("エラー: TEST_15/16 の結果が見つかりません") + return 1 + OUT_DIR.mkdir(parents=True, exist_ok=True) + with (OUT_DIR / "minmax_definitions.csv").open("w", newline="", encoding="utf-8") as fh: + writer = csv.DictWriter(fh, fieldnames=list(rows[0].keys())) + writer.writeheader() + writer.writerows(rows) + draw(rows, OUT_DIR / "minmax_definitions.png") + + print(f"\n{'形状':>7} {'サイズ':>6} " + " ".join(f"{d:>11}" for d in DEFINITIONS) + + f" {'輪帯内−ロバスト':>14}") + for r in rows: + print(f"{r['shape']:>7} {r['size_mm']:6g} " + + " ".join(f"{r[d]:11.4f}" for d in DEFINITIONS) + + f" {r['ring_inner'] - r['robust']:+14.4f}") + return 0 + + +if __name__ == "__main__": + sys.exit(main())