Newer
Older
SmTIAS-LightSim / scripts / build_offaxis_comparison.py
"""オフ軸 LED スイープの解析: 非対称の量とサイズによる減衰を球・立方体で比較する(TEST_17).

LED を +X 側の側壁に置いたことで生じる非対称を,次の 3 層で測る.

1. **m=1 符号付き係数 a1**(輪帯ごとの方位フーリエ・逆分散重み集約): ズレ方向が +X と
   分かっているので,cos φ 成分の符号付き振幅が「LED 側がどれだけ明るい/暗いか」を直接
   与える.構成は y=0 の鏡映対称なので直交成分 b1 は対称性からゼロのはず=ノイズ対照
   (m=4 解析と同じ設計).軸上(on-axis)の同条件の |m1| が方法の背景値を与える
2. **輝度重心の X オフセット**(metrics.csv の intensity_centroid_x_norm)
3. **構造 CoV / 輪帯 min/max**: 非対称による総合的な均一性の悪化量(軸上との差)

仮説「大きいほど光源位置の影響が薄まる」が正しければ,a1・重心オフセット・
Δ構造CoV がキャビティ内面積とともに減衰するはずである.

実行方法:
    uv run python scripts/build_offaxis_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 analyze_fourfold_component import aggregate, 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_cube_size import OUTPUT_BASE as CUBE_OUTPUT_BASE  # noqa: E402
from sweep_cube_size import RUN_DIR as CUBE_RUN_DIR  # noqa: E402
from sweep_cube_size import job_name as cube_job_name  # noqa: E402
from sweep_offaxis_led import (  # noqa: E402
    CUBE_SIZES,
    OUTPUT_BASE,
    SPHERE_RADII,
    job_name,
)
from sweep_sphere_radius import OUTPUT_BASE as SPHERE_OUTPUT_BASE  # noqa: E402
from sweep_sphere_radius import job_name as sphere_job_name  # noqa: E402

plt.rcParams["font.family"] = ["Yu Gothic", "MS Gothic", "DejaVu Sans"]

SPHERE_ONAXIS_DIR = SPHERE_OUTPUT_BASE / "board_baffle_similar_radial"

SHAPE_COLOR = {"sphere": "tab:blue", "cube": "tab:red"}
SHAPE_LABEL = {"sphere": "球 similar_radial", "cube": "立方体"}


def onaxis_dir(shape: str, size: float) -> Path:
    """軸上(on-axis)スイープの同条件ディレクトリを返す."""
    if shape == "sphere":
        return SPHERE_ONAXIS_DIR / sphere_job_name(size)
    return CUBE_OUTPUT_BASE / CUBE_RUN_DIR / cube_job_name(size)


def read_summary_rows() -> list[dict]:
    """オフ軸スイープの summary.csv を読み,欠損条件を除いて返す."""
    path = OUTPUT_BASE / "summary.csv"
    with path.open(encoding="utf-8") as fh:
        rows = []
        for row in csv.DictReader(fh):
            parsed: dict = {"shape": row["shape"]}
            parsed.update(
                {k: (float(v) if v else float("nan")) for k, v in row.items() if k != "shape"}
            )
            rows.append(parsed)
    return sorted(rows, key=lambda r: (r["shape"], r["size_mm"]))


def load_roi_for(path: Path) -> np.ndarray:
    """luminance.npy から ROI を読み出す(analyze_fourfold_component と同じ経路)."""
    from analyze_fourfold_component import load_roi

    return load_roi(path)


def analyze_condition(shape: str, size: float) -> dict | None:
    """1 条件のオフ軸・軸上の非対称指標と均一性指標を集める."""
    off_dir = OUTPUT_BASE / job_name(shape, size)
    on_dir = onaxis_dir(shape, size)
    if not (off_dir / "luminance.npy").exists() or not (on_dir / "luminance.npy").exists():
        return None

    result: dict = {"shape": shape, "size_mm": size}
    for tag, directory in (("off", off_dir), ("on", on_dir)):
        lum = directory / "luminance.npy"
        rings = ring_harmonics(load_roi_for(lum))
        for key in ("a1", "b1"):
            result[f"{tag}_{key}"], result[f"{tag}_{key}_se"] = aggregate(rings, key)
        _, structure, _ = decompose_cov(lum)
        centers, means, _ = radial_stats(lum, 20)
        result[f"{tag}_structure_cov"] = structure
        result[f"{tag}_ring_min_max"] = ring_min_max_ratio(centers, means)
        with (directory / "metrics.csv").open(encoding="utf-8") as fh:
            row = next(csv.DictReader(fh))
        result[f"{tag}_centroid_x"] = float(row["intensity_centroid_x_norm"])
        result[f"{tag}_gradient"] = float(row["gradient_magnitude"])
    return result


def draw_maps(fig, gs, rows: list[dict]) -> None:
    """オフ軸の輝度マップを形状別に 1 行ずつ並べる(輪帯ピーク=1.0・共通スケール)."""
    for row_idx, shape in enumerate(("sphere", "cube")):
        selected = [r for r in rows if r["shape"] == shape]
        maps = [
            (r["size_mm"], normalized_map(OUTPUT_BASE / job_name(shape, r["size_mm"])
                                          / "luminance.npy"))
            for r in selected
        ]
        if not maps:
            continue
        stacked = np.concatenate([m.ravel() for _, m in maps])
        vmin = float(np.percentile(stacked, 0.5))
        for i, (size, data) in enumerate(maps[:5]):
            ax = fig.add_subplot(gs[row_idx, i])
            image = ax.imshow(data, cmap="hot", vmin=vmin, vmax=1.0)
            prefix = "r" if shape == "sphere" else "L"
            ax.set_title(f"{SHAPE_LABEL[shape]} {prefix}={size:g}", fontsize=9)
            ax.axis("off")
            if i == len(maps[:5]) - 1:
                fig.colorbar(image, ax=ax, fraction=0.046, label="最大輪帯=1.0")


def draw_metric(ax, results: list[dict], key_off: str, key_on: str | None,
                ylabel: str, title: str, se_key: str | None = None,
                scale: float = 1.0) -> None:
    """指標を内面積(対数軸)に対して球・立方体×オフ軸/軸上で描く."""
    for shape in ("sphere", "cube"):
        rows = [r for r in results if r["shape"] == shape]
        areas = [r["area"] for r in rows]
        values = [r[key_off] * scale for r in rows]
        if se_key:
            errors = [r[se_key] * scale for r in rows]
            ax.errorbar(areas, values, yerr=errors, marker="o", capsize=3,
                        color=SHAPE_COLOR[shape], label=f"{SHAPE_LABEL[shape]}(オフ軸)")
        else:
            ax.plot(areas, values, marker="o", color=SHAPE_COLOR[shape],
                    label=f"{SHAPE_LABEL[shape]}(オフ軸)")
        if key_on:
            ax.plot(areas, [r[key_on] * scale for r in rows], marker="s", ls="--", ms=4,
                    color=SHAPE_COLOR[shape], alpha=0.55,
                    label=f"{SHAPE_LABEL[shape]}(軸上)")
    ax.axhline(0.0, color="gray", lw=0.8, ls="--")
    ax.set_xscale("log")
    ax.set_xlabel("キャビティ内面積 (mm²・対数)")
    ax.set_ylabel(ylabel)
    ax.set_title(title, fontsize=11)
    ax.legend(fontsize=7)
    ax.grid(alpha=0.3, which="both")


def print_table(results: list[dict]) -> None:
    """非対称指標の集約表を出力する."""
    print("\n=== 非対称の指標(a1 > 0 = LED 側(+X)が明るい・輪帯平均比 %)===")
    print(f"{'条件':>14} {'a1 オフ軸 [%]':>15} {'b1(対照)[%]':>14} {'a1 軸上 [%]':>12} "
          f"{'重心X(off)':>10} {'構造CoV off':>11} {'on':>8} {'Δ':>8}")
    for r in results:
        name = job_name(r["shape"], r["size_mm"])
        print(f"{name:>14} {r['off_a1'] * 100:8.3f}±{r['off_a1_se'] * 100:5.3f} "
              f"{r['off_b1'] * 100:8.3f}±{r['off_b1_se'] * 100:4.3f} "
              f"{r['on_a1'] * 100:12.3f} {r['off_centroid_x']:10.5f} "
              f"{r['off_structure_cov']:11.4f} {r['on_structure_cov']:8.4f} "
              f"{r['off_structure_cov'] - r['on_structure_cov']:+8.4f}")


def main() -> int:
    """オフ軸 LED スイープの比較図を生成する."""
    parser = argparse.ArgumentParser(description="オフ軸 LED スイープの解析図を生成")
    parser.add_argument("--output", type=Path, default=None, help="出力 PNG パス")
    args = parser.parse_args()

    summary = read_summary_rows()
    results = []
    for shape, sizes in (("sphere", SPHERE_RADII), ("cube", CUBE_SIZES)):
        for size in sizes:
            res = analyze_condition(shape, size)
            if res is None:
                continue
            match = [r for r in summary if r["shape"] == shape and r["size_mm"] == size]
            res["area"] = match[0]["cavity_area_mm2"] if match else float("nan")
            results.append(res)
    if not results:
        print("エラー: オフ軸スイープの結果がありません(先にレンダを実行してください)")
        return 1

    output = args.output or OUTPUT_BASE / "offaxis_comparison.png"
    fig = plt.figure(figsize=(20, 14))
    gs = fig.add_gridspec(3, 6, height_ratios=[1.1, 1.1, 1])
    draw_maps(fig, gs, results)
    draw_metric(
        fig.add_subplot(gs[2, 0:2]), results, "off_a1", "on_a1",
        "a1(輪帯平均比 %・正=LED 側が明るい)",
        "m=1 非対称成分 vs 内面積", se_key="off_a1_se", scale=100.0,
    )
    draw_metric(
        fig.add_subplot(gs[2, 2:4]), results, "off_centroid_x", "on_centroid_x",
        "輝度重心 X オフセット(正規化)", "重心の LED 側への寄り vs 内面積",
    )
    draw_metric(
        fig.add_subplot(gs[2, 4:6]), results, "off_structure_cov", "on_structure_cov",
        "構造 CoV", "総合均一性: オフ軸 vs 軸上",
    )
    fig.suptitle(
        "オフ軸 LED(側壁・板中心向き・バッフル付き)の影響とサイズによる減衰"
        " 輝度マップは輪帯ピーク正規化・LED は +X 側", fontsize=15,
    )
    fig.tight_layout()
    output.parent.mkdir(parents=True, exist_ok=True)
    fig.savefig(output, dpi=100)
    print(f"出力: {output}")
    print_table(results)
    return 0


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