Newer
Older
SmTIAS-LightSim / scripts / build_minmax_asymmetry_comparison.py
"""輪帯平均と画像全体の最小/最大比を,非対称な照明配置(TEST_18)で比較する.

輪帯平均の min/max(`ring_min_max`)は方位方向に平均するので,LED をオフ軸に置いて
生じる「片側だけ明るい」非対称を打ち消してしまう.画像全体のノイズ補正付き min/max
(`fullimage_min_max`,`evaluation/fullimage.py`)は方位方向に平均しないので,
その非対称が値に乗るはずである.本スクリプトはその差が実際に非対称の大きさ |a1|
(m=1 方位フーリエ成分)とともに開くことを,位置スイープ(4 ペア × 13 θ)で確かめる.

入力は `build_led_placement_comparison.py` が書いた analysis.csv(両定義の列を含む).
luminance.npy の読み直しはしない.

出力(output/led_placement_sweep/):
    minmax_asymmetry.png / minmax_asymmetry.csv

実行方法:
    uv run python scripts/build_minmax_asymmetry_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"))

import build_led_placement_comparison as placement  # noqa: E402
import sweep_led_placement as sweep  # noqa: E402
from sweep_led_placement import MATERIALS, SHAPES, SIZE_PAIRS, set_material  # noqa: E402

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

CSV_COLUMNS = [
    "name", "shape", "size_mm", "theta_deg", "m1_abs",
    "ring_min_max", "fullimage_min_max", "diff", "structure_cov",
]

# θ=0 は軸上(対称)・θ=120 は板の直下(最も非対称).標準出力で対比する
THETA_SYMMETRIC = 0.0
THETA_OFFAXIS = 120.0


def load_rows(analysis_csv: Path) -> list[dict]:
    """analysis.csv から位置スイープの行を読み,両定義と差を持つ dict のリストを返す."""
    rows: list[dict] = []
    with analysis_csv.open(encoding="utf-8", newline="") as fh:
        for raw in csv.DictReader(fh):
            if raw.get("group") != "position":
                continue
            if not raw.get("fullimage_min_max"):
                continue
            ring = float(raw["ring_min_max"])
            full = float(raw["fullimage_min_max"])
            rows.append({
                "name": raw["name"],
                "shape": raw["shape"],
                "size_mm": float(raw["size_mm"]),
                "theta_deg": float(raw["theta_deg"]),
                "m1_abs": float(raw["m1_abs"]),
                "ring_min_max": ring,
                "fullimage_min_max": full,
                "diff": ring - full,
                "structure_cov": float(raw["structure_cov"]),
            })
    return sorted(rows, key=lambda r: (r["shape"], r["size_mm"], r["theta_deg"]))


def series(rows: list[dict], shape: str, size: float) -> list[dict]:
    """同じ形状・サイズの系列を θ 昇順で取り出す."""
    return [r for r in rows if r["shape"] == shape and r["size_mm"] == size]


def draw(rows: list[dict], output: Path) -> None:
    """θ に対する両定義の比較(形状ごと)と,差 vs 非対称の大きさの散布図を描く."""
    fig, axes = plt.subplots(1, 3, figsize=(19, 5.6))
    colors = placement.size_colors(len(SIZE_PAIRS))

    for ax, shape in zip(axes[:2], SHAPES, strict=False):
        for i, pair in enumerate(SIZE_PAIRS):
            size = placement.pair_size(pair, shape)
            sel = series(rows, shape, size)
            if not sel:
                continue
            theta = [r["theta_deg"] for r in sel]
            label = placement.size_label(shape, size)
            ax.plot(theta, [r["fullimage_min_max"] for r in sel], color=colors[i],
                    marker="o", ms=4, lw=1.6, ls="-", label=f"{label} 画像全体")
            ax.plot(theta, [r["ring_min_max"] for r in sel], color=colors[i],
                    marker="x", ms=4, lw=1.2, ls=":", alpha=0.75, label=f"{label} 輪帯平均")
        ax.set_xlabel("LED 極角 θ [deg](0=底の中心・90=側壁の中段・120=板の直下)")
        ax.set_ylabel("最小/最大比(1 が理想)")
        ax.set_title(f"{placement.SHAPE_LABEL[shape]}:実線=画像全体(ノイズ補正)"
                     "/点線=輪帯平均", fontsize=11)
        ax.grid(alpha=0.3)
        ax.legend(fontsize=7, ncol=2)

    ax = axes[2]
    for shape in SHAPES:
        for i, pair in enumerate(SIZE_PAIRS):
            size = placement.pair_size(pair, shape)
            sel = series(rows, shape, size)
            if not sel:
                continue
            ax.scatter([r["m1_abs"] * 100 for r in sel], [r["diff"] for r in sel],
                       color=colors[i], marker=placement.SHAPE_STYLE[shape]["marker"],
                       s=26, alpha=0.8,
                       label=f"{placement.SHAPE_LABEL[shape]} {placement.size_label(shape, size)}")
    ax.axhline(0.0, color="gray", lw=0.8, ls="--")
    ax.set_xlabel("非対称の大きさ |a1|(輪帯平均比 %)")
    ax.set_ylabel("輪帯平均 − 画像全体(差)")
    ax.set_title("方位非対称が大きいほど輪帯平均は均一に見える(差が開く)", fontsize=11)
    ax.grid(alpha=0.3)
    ax.legend(fontsize=7, ncol=2)

    fig.suptitle("最小/最大比の定義比較 — 輪帯平均は方位方向に平均するため"
                 "オフ軸 LED の非対称を過小評価する(TEST_18 位置スイープ)", fontsize=14)
    fig.tight_layout()
    output.parent.mkdir(parents=True, exist_ok=True)
    fig.savefig(output, dpi=110)
    plt.close(fig)
    print(f"出力: {output}")


def report(rows: list[dict]) -> None:
    """形状ごとに軸上/オフ軸での差と,差と |a1| の相関係数を表示する."""
    print("\n=== 輪帯平均 − 画像全体(差)===")
    print(f"{'形状':>6} {'サイズ':>7} {'θ=0 の差':>10} {'θ=120 の差':>11} "
          f"{'差の最大':>9} {'相関 r(差, |a1|)':>16}")
    for shape in SHAPES:
        for pair in SIZE_PAIRS:
            size = placement.pair_size(pair, shape)
            sel = series(rows, shape, size)
            if len(sel) < 3:
                continue
            diffs = np.array([r["diff"] for r in sel])
            m1 = np.array([r["m1_abs"] for r in sel])
            at_sym = next((r["diff"] for r in sel if r["theta_deg"] == THETA_SYMMETRIC), np.nan)
            at_off = next((r["diff"] for r in sel if r["theta_deg"] == THETA_OFFAXIS), np.nan)
            corr = float(np.corrcoef(diffs, m1)[0, 1]) if diffs.std() > 0 else float("nan")
            print(f"{placement.SHAPE_LABEL[shape]:>6} {size:7g} {at_sym:10.4f} {at_off:11.4f} "
                  f"{diffs.max():9.4f} {corr:16.3f}")

    for shape in SHAPES:
        sel = [r for r in rows if r["shape"] == shape]
        if len(sel) < 3:
            continue
        diffs = np.array([r["diff"] for r in sel])
        m1 = np.array([r["m1_abs"] for r in sel])
        corr = float(np.corrcoef(diffs, m1)[0, 1])
        print(f"{placement.SHAPE_LABEL[shape]}(全サイズ {len(sel)} 条件): "
              f"差と |a1| の相関係数 r={corr:.3f}")


def write_csv(rows: list[dict], output: Path) -> None:
    """比較表を CSV に書き出す."""
    output.parent.mkdir(parents=True, exist_ok=True)
    with output.open("w", newline="", encoding="utf-8") as fh:
        writer = csv.DictWriter(fh, fieldnames=CSV_COLUMNS)
        writer.writeheader()
        for r in rows:
            writer.writerow({k: r[k] for k in CSV_COLUMNS})
    print(f"出力: {output}")


def main() -> int:
    """位置スイープについて両定義の min/max 比を比較する図と表を出す."""
    parser = argparse.ArgumentParser(description="輪帯平均と画像全体の min/max 比を比較")
    parser.add_argument("--material", choices=tuple(MATERIALS), default="legacy",
                        help="解析する材質系列(スイープの --material と同じ)")
    args = parser.parse_args()
    set_material(args.material)
    output_base = sweep.OUTPUT_BASE

    analysis_csv = output_base / "analysis.csv"
    if not analysis_csv.exists():
        print(f"エラー: {analysis_csv} がありません."
              "先に build_led_placement_comparison.py を実行してください")
        return 1

    rows = load_rows(analysis_csv)
    if not rows:
        print(f"エラー: {analysis_csv} に fullimage_min_max 列の入った位置スイープの行が"
              "ありません.build_led_placement_comparison.py を再実行してください")
        return 1

    write_csv(rows, output_base / "minmax_asymmetry.csv")
    draw(rows, output_base / "minmax_asymmetry.png")
    report(rows)
    return 0


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