Newer
Older
SmTIAS-Evaluation / scripts / compare_composite_vs_sum.py
"""全灯点灯の実測と単灯 3 枚の加算合成を比較するスクリプト.

同一セッションで撮影した「全灯点灯」1 枚と「下・左・右の単灯」3 枚(いずれも白板 DNG)
について,単灯輝度マップの画素ごとの加算(下+左+右)が全灯実測を再現するかを検証する.
照明が線形加算で合成される(灯間の相互作用がない)なら 実測 ≈ 加算 となるはずで,
「3 灯の相補で合成均一性が成立する」仮説(TEST_04)の定量検証にあたる.

構成:
- 上段: 全灯実測 / 単灯加算 の輝度マップ(共通絶対カラースケール)+ 比マップ(実測/加算)
- 下段左: 動径プロファイル重ね描き(実測・加算の絶対値)
- 下段右: 数値表(平均・CoV・max/min 比・中心/周辺比・勾配・動径 min/max 比 + 比の統計)

すべて linear + LSC 補正・ROI=smtias.whiteboard で揃えて算出する.

使い方:
    python scripts/compare_composite_vs_sum.py --all <dng> --bottom <dng> --left <dng> --right <dng>
"""

import argparse
import sys
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np

PROJECT_ROOT = Path(__file__).resolve().parent.parent
sys.path.insert(0, str(PROJECT_ROOT))

from src.analysis.spatial import (  # noqa: E402
    calc_distance_map,
    calc_radial_min_max_ratio,
    calc_radial_profile,
    calc_spatial_uniformity,
)
from src.analysis.uniformity import calc_uniformity  # noqa: E402
from src.config import load_roi_config  # noqa: E402
from src.io.dng_pipeline import dng_luminance  # noqa: E402
from src.io.loader import extract_roi  # noqa: E402

ROI_CONFIG = PROJECT_ROOT / "config" / "roi_config.json"
ROI_KEY = "smtias.whiteboard"
DEFAULT_OUTPUT = PROJECT_ROOT / "output" / "smtias" / "figures" / "compare_composite_vs_sum.png"

# 共通絶対カラースケールの設定(外れ値に強いパーセンタイルで vmin/vmax を決める)
SCALE_PERCENTILES = (1.0, 99.0)


def _luminance(dng_path: Path, roi: dict) -> np.ndarray:
    """DNG を linear + LSC 補正で輝度化し ROI を切り出す(portrait)."""
    luma_full = dng_luminance(str(dng_path), apply_lsc=True)
    return extract_roi(luma_full, roi["x"], roi["y"], roi["width"], roi["height"])


def _metrics(luminance: np.ndarray) -> dict:
    """均一性・空間・動径の指標をまとめて算出する."""
    distance_map = calc_distance_map(*luminance.shape)
    profile = calc_radial_profile(luminance, distance_map)
    radial = calc_radial_min_max_ratio(profile)
    return {
        **calc_uniformity(luminance),
        **calc_spatial_uniformity(luminance)["zone_stats"],
        "radial_min_max_ratio": radial["radial_min_max_ratio"],
        "profile": profile,
    }


def _parse_args() -> argparse.Namespace:
    """CLI 引数を解釈する."""
    parser = argparse.ArgumentParser(
        description="全灯実測と単灯 3 枚の加算合成(下+左+右)を比較する."
    )
    parser.add_argument("--all", required=True, dest="all_on", help="全灯点灯の DNG パス")
    parser.add_argument("--bottom", required=True, help="下のみ点灯の DNG パス")
    parser.add_argument("--left", required=True, help="左のみ点灯の DNG パス")
    parser.add_argument("--right", required=True, help="右のみ点灯の DNG パス")
    parser.add_argument("--output", default=str(DEFAULT_OUTPUT), help="出力する比較図のパス")
    parser.add_argument(
        "--title-suffix", default="", help="図タイトルの末尾に付ける文字列(例: '(w1p4)')"
    )
    return parser.parse_args()


def main() -> None:
    """全灯実測 vs 単灯加算の比較図を生成して保存する."""
    args = _parse_args()
    roi = load_roi_config(str(ROI_CONFIG), ROI_KEY)
    if roi is None:
        raise SystemExit(f"ROI '{ROI_KEY}' が config に見つかりません")

    measured = _luminance(Path(args.all_on), roi)
    singles = [_luminance(Path(p), roi) for p in (args.bottom, args.left, args.right)]
    summed = np.sum(singles, axis=0)

    labels = ["Measured (all LEDs)", "Sum (bottom+left+right)"]
    lumas = [measured, summed]
    metrics = [_metrics(luma) for luma in lumas]

    # 比マップ(実測 / 加算)とその統計
    ratio = measured / np.maximum(summed, 1e-6)
    ratio_stats = {
        "mean": float(np.mean(ratio)),
        "median": float(np.median(ratio)),
        "p5": float(np.percentile(ratio, 5)),
        "p95": float(np.percentile(ratio, 95)),
    }

    flat_all = np.concatenate([luma.ravel() for luma in lumas])
    vmin, vmax = np.percentile(flat_all, SCALE_PERCENTILES)

    fig = plt.figure(figsize=(15, 9))
    gs = fig.add_gridspec(2, 3, height_ratios=[1.0, 0.85], hspace=0.28, wspace=0.22)

    # 上段: 実測・加算(共通絶対スケール)
    map_axes = []
    im = None
    for col, (label, luma) in enumerate(zip(labels, lumas)):
        ax = fig.add_subplot(gs[0, col])
        im = ax.imshow(luma, cmap="hot", interpolation="nearest", vmin=vmin, vmax=vmax)
        ax.set_title(f"{label}  (mean={metrics[col]['mean']:.1f})", fontsize=11)
        ax.set_xlabel("X (px)")
        if col == 0:
            ax.set_ylabel("Y (px)")
        map_axes.append(ax)
    cbar = fig.colorbar(im, ax=map_axes, fraction=0.020, pad=0.02)
    cbar.set_label(f"Luminance (common scale, p{SCALE_PERCENTILES[0]:g}-p{SCALE_PERCENTILES[1]:g})")

    # 上段右: 比マップ(実測/加算.1.0 中心の発散スケール)
    ax_ratio = fig.add_subplot(gs[0, 2])
    half_range = max(abs(ratio_stats["p5"] - 1.0), abs(ratio_stats["p95"] - 1.0))
    im_r = ax_ratio.imshow(
        ratio, cmap="RdBu_r", interpolation="nearest",
        vmin=1.0 - half_range, vmax=1.0 + half_range,
    )
    ax_ratio.set_title(
        f"Ratio measured/sum  (median={ratio_stats['median']:.3f})", fontsize=11
    )
    ax_ratio.set_xlabel("X (px)")
    fig.colorbar(im_r, ax=ax_ratio, fraction=0.040, pad=0.02).set_label("Ratio")

    # 下段左: 動径プロファイル重ね描き(絶対値)
    ax_prof = fig.add_subplot(gs[1, 0:2])
    for label, color, marker, m in zip(labels, ["darkorange", "steelblue"], ["o", "s"], metrics):
        profile = m["profile"]
        ax_prof.plot(
            profile[:, 0], profile[:, 1], marker=marker, ms=4, color=color,
            label=f"{label}  CoV={m['cov']:.3f}",
        )
    ax_prof.set_title("Radial profile (absolute luminance)", fontsize=11)
    ax_prof.set_xlabel("Normalized distance from center")
    ax_prof.set_ylabel("Mean luminance")
    ax_prof.grid(True, linestyle="--", alpha=0.5)
    ax_prof.legend(loc="lower left", fontsize=9)

    # 下段右: 数値表
    ax_tbl = fig.add_subplot(gs[1, 2])
    ax_tbl.axis("off")
    metric_rows = [
        ("Mean luminance", "mean", "{:.2f}"),
        ("CoV", "cov", "{:.3f}"),
        ("Max/min ratio", "max_min_ratio", "{:.3f}"),
        ("Center/periphery", "center_periphery_ratio", "{:.3f}"),
        ("Gradient (%)", "gradient_magnitude", "{:.2f}"),
        ("Radial min/max", "radial_min_max_ratio", "{:.3f}"),
    ]
    rows = [
        [name] + [fmt.format(m[key]) for m in metrics]
        for name, key, fmt in metric_rows
    ]
    rows.append(["Ratio mean [p5, p95]",
                 f"{ratio_stats['mean']:.3f}",
                 f"[{ratio_stats['p5']:.3f}, {ratio_stats['p95']:.3f}]"])
    table = ax_tbl.table(
        cellText=rows,
        colLabels=["Metric", "Measured", "Sum"],
        cellLoc="center",
        loc="center",
    )
    table.auto_set_font_size(False)
    table.set_fontsize(9)
    table.scale(1.0, 1.6)
    for j in range(3):
        table[0, j].set_facecolor("#dfe6ee")
        table[0, j].set_text_props(weight="bold")

    title = "All-LEDs measured vs single-LED sum (bottom+left+right)"
    if args.title_suffix:
        title += f" {args.title_suffix}"
    fig.suptitle(title, fontsize=14, y=0.98)

    output_path = Path(args.output)
    output_path.parent.mkdir(parents=True, exist_ok=True)
    fig.savefig(str(output_path), dpi=150, bbox_inches="tight")
    plt.close(fig)

    print(f"全灯実測 vs 単灯加算の比較図を保存しました: {output_path}")
    for label, m in zip(labels, metrics):
        print(f"  {label}: mean={m['mean']:.2f}, CoV={m['cov']:.4f}, "
              f"radial min/max={m['radial_min_max_ratio']:.4f}")
    print(f"  Ratio measured/sum: mean={ratio_stats['mean']:.4f}, "
          f"median={ratio_stats['median']:.4f}, "
          f"p5={ratio_stats['p5']:.4f}, p95={ratio_stats['p95']:.4f}")


if __name__ == "__main__":
    main()