"""実機 DNG とシミュレーションを，同じノイズ処理で min・放射プロファイル・形について突き合わせる.

目的: 舌撮影デバイスの基準「最大を 1 としたときの最小 ≥ 0.9」に対して，実機を作る前に
シミュレーションの min が参考になるかを確かめる．実機（SmTIAS 定量モード DNG）と
シミュ（校正材質の基準シーン）に**同じ処理**を掛け，

1. min（画像全体・ノイズ補正付き最小/最大比．参考に補正なし・輪帯平均も）
2. 20 分割の放射プロファイル（ROI 平均で正規化した各ビンの値の差）
3. 輝度マップの形（ロバスト最大で正規化したマップの差分・SSIM・NCC・方位フーリエ m=1/2/4）

を数値で比べる．

実機の読み込みは TEST_14 と同じ経路（load_image_dng → lscMap 順適用 → Rec.709 → rot90(-1)
→ 白板 ROI）．実機 ROI（1400×1409 px）はシミュ ROI（464×450 px）へリサイズしてから同じ
処理を掛ける（平滑化 σ を画素で揃えるため．リサイズは約 3×3 画素の平均になりセンサノイズ
も落ちるが，ノイズ補正はリサイズ後の残差から推定するので整合する）．

実行方法:
    uv run python scripts/compare_real_vs_sim_metrics.py
    uv run python scripts/compare_real_vs_sim_metrics.py \
        --sim output/bsdf_grid_fix/r0.85_w0.00_a0.1/luminance.npy --label "ρ=0.85" --tag r085
"""

from __future__ import annotations

import argparse
import csv
import json
import sys
from pathlib import Path

import matplotlib
import numpy as np
from scipy.ndimage import gaussian_filter, median_filter
from skimage.transform import resize

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 MAX_DISTANCE  # noqa: E402

from sm_tias_lightsim.evaluation.dng_loader import load_image_dng  # noqa: E402
from sm_tias_lightsim.evaluation.fullimage import calc_fullimage_min_max_ratio  # noqa: E402
from sm_tias_lightsim.evaluation.radial import calc_radial_profile  # noqa: E402
from sm_tias_lightsim.evaluation.roi_presets import (  # noqa: E402
    REAL_ROI,
    SIM_ROI,
    extract_roi_xywh,
)
from sm_tias_lightsim.evaluation.shape_metrics import calc_ncc, calc_ssim  # noqa: E402
from sm_tias_lightsim.evaluation.zones import compute_normalized_distance_map  # noqa: E402

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

ROOT = Path(__file__).resolve().parents[1]
REAL_DNG = ROOT / "reference_measurement" / "quantitative" / "MiniTIAS_QM_20260601_124500.dng"
REAL_META = REAL_DNG.with_name("MiniTIAS_QM_20260601_124500.meta.json")
DEFAULT_SIM = ROOT / "output" / "bsdf_grid_fix" / "r0.82_w0.00_a0.1" / "luminance.npy"
OUT_DIR = ROOT / "output" / "real_vs_sim"

N_BINS = 20
STRUCTURE_SIGMA = 12.0
DISPLAY_MEDIAN = 5
DISPLAY_SIGMA = 4.0
ROBUST_PERCENTILE = 99.0
ISO_THRESHOLD = 0.9


# ---------------------------------------------------------------------------
# 読み込み
# ---------------------------------------------------------------------------


def _upsample(grid: np.ndarray, shape: tuple[int, int]) -> np.ndarray:
    """(rows, cols) の格子を画像サイズへバイリニア補間する（参照スクリプトと同一）."""
    rows, cols = grid.shape
    h, w = shape
    ys = np.linspace(0, rows - 1, h)
    xs = np.linspace(0, cols - 1, w)
    y0 = np.floor(ys).astype(int)
    y1 = np.minimum(y0 + 1, rows - 1)
    x0 = np.floor(xs).astype(int)
    x1 = np.minimum(x0 + 1, cols - 1)
    wy = (ys - y0)[:, None]
    wx = (xs - x0)[None, :]
    top = grid[y0][:, x0] * (1 - wx) + grid[y0][:, x1] * wx
    bot = grid[y1][:, x0] * (1 - wx) + grid[y1][:, x1] * wx
    return top * (1 - wy) + bot * wy


def lsc_gain_rgb(meta_path: Path, shape: tuple[int, int]) -> np.ndarray:
    """meta.json の lscMap（rows, cols, 4=[R, G_even, G_odd, B]）→ (H, W, 3) の RGB ゲイン."""
    meta = json.loads(meta_path.read_text(encoding="utf-8"))
    rows, cols = meta["lscMapRowCount"], meta["lscMapColumnCount"]
    m = np.array(meta["lscMap"], dtype=np.float64).reshape(rows, cols, 4)
    return np.stack([_upsample(m[:, :, 0], shape),
                     _upsample((m[:, :, 1] + m[:, :, 2]) / 2.0, shape),
                     _upsample(m[:, :, 3], shape)], axis=2)


def load_real_roi() -> np.ndarray:
    """実機 DNG → linear raw → LSC 順適用 → Rec.709 × 255 → portrait → 白板 ROI."""
    rgb = load_image_dng(str(REAL_DNG), output_format="rgb_float32").astype(np.float64)
    rgb = rgb * lsc_gain_rgb(REAL_META, rgb.shape[:2])
    luma = 0.2126 * rgb[..., 0] + 0.7152 * rgb[..., 1] + 0.0722 * rgb[..., 2]
    luma = np.rot90(luma * 255.0, -1)
    return extract_roi_xywh(luma, REAL_ROI)


def load_sim_roi(path: Path) -> np.ndarray:
    """シミュ luminance.npy → ROI."""
    arr = np.load(path).astype(np.float64)
    if arr.ndim == 3:
        arr = 0.2126 * arr[..., 0] + 0.7152 * arr[..., 1] + 0.0722 * arr[..., 2]
    return extract_roi_xywh(arr, SIM_ROI)


# ---------------------------------------------------------------------------
# 指標
# ---------------------------------------------------------------------------


def structure_and_noise_cov(roi: np.ndarray, seed: int = 0) -> tuple[float, float, float]:
    """(生 CoV, 構造 CoV, ノイズ CoV)．build_radius_sweep_visual.decompose_cov と同じ手順."""
    mean = roi.mean()
    residual = roi - gaussian_filter(median_filter(roi, size=DISPLAY_MEDIAN), sigma=DISPLAY_SIGMA)
    smoothed = gaussian_filter(median_filter(roi, size=DISPLAY_MEDIAN), sigma=STRUCTURE_SIGMA)
    flat = residual.ravel().copy()
    np.random.default_rng(seed).shuffle(flat)
    leak = gaussian_filter(median_filter(flat.reshape(roi.shape), size=DISPLAY_MEDIAN),
                           sigma=STRUCTURE_SIGMA).std() / mean
    structure = float(np.sqrt(max((smoothed.std() / smoothed.mean()) ** 2 - leak**2, 0.0)))
    return float(roi.std() / mean), structure, float(residual.std() / mean)


def radial_profile_norm(roi: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """20 ビン放射プロファイル（ROI 平均で正規化）の (ビン中心距離, 値) を返す."""
    profile = calc_radial_profile(roi, n_bins=N_BINS)
    return profile[:, 0], profile[:, 1] / roi.mean()


def ring_min_max(roi: np.ndarray) -> float:
    """輪帯平均 min/max（最外ビン除外．TEST_18 と同じ）."""
    profile = calc_radial_profile(roi, n_bins=N_BINS)
    valid = (~np.isnan(profile[:, 1])) & (profile[:, 0] <= MAX_DISTANCE)
    return float(profile[valid, 1].min() / profile[valid, 1].max())


def display_map(roi: np.ndarray) -> np.ndarray:
    """表示用: メディアン 5 → ガウス σ=4 で平滑化し，99 パーセンタイルを 1 に正規化."""
    smooth = gaussian_filter(median_filter(roi, size=DISPLAY_MEDIAN), sigma=DISPLAY_SIGMA)
    return smooth / float(np.percentile(smooth, ROBUST_PERCENTILE))


def metrics_for(roi: np.ndarray) -> dict[str, float]:
    """1 枚の ROI について全指標を計算する."""
    pixel_cov, structure, noise = structure_and_noise_cov(roi)
    full = calc_fullimage_min_max_ratio(roi, sigma=STRUCTURE_SIGMA)
    smooth = gaussian_filter(roi, sigma=STRUCTURE_SIGMA)
    rings = ring_harmonics(roi, (1, 2, 4))
    out = {
        "mean": float(roi.mean()),
        "pixel_cov": pixel_cov,
        "structure_cov": structure,
        "noise_cov": noise,
        "min_fullimage": full["fullimage_min_max_ratio"],
        "min_fullimage_raw": float(smooth.min() / smooth.max()),
        "min_ring": ring_min_max(roi),
        "min_pixel": float(roi.min() / roi.max()),
        "noise_span_over_mean": float((full["fullimage_noise_max"] - full["fullimage_noise_min"])
                                      / roi.mean()),
    }
    for m in (1, 2, 4):
        out[f"a{m}"], out[f"se{m}"] = aggregate(rings, f"a{m}")
        out[f"b{m}"], _ = aggregate(rings, f"b{m}")
    return out


def shape_agreement(sim_roi: np.ndarray, real_roi: np.ndarray) -> dict[str, float]:
    """形の一致度: 正規化マップの差分・SSIM・NCC・放射プロファイルの差."""
    sim_map, real_map = display_map(sim_roi), display_map(real_roi)
    diff = sim_map - real_map
    dist = compute_normalized_distance_map(*sim_roi.shape)
    inner = dist <= MAX_DISTANCE
    _, sim_prof = radial_profile_norm(sim_roi)
    _, real_prof = radial_profile_norm(real_roi)
    valid = ~(np.isnan(sim_prof) | np.isnan(real_prof))
    prof_diff = (sim_prof - real_prof)[valid]
    return {
        "map_rms_diff": float(np.sqrt(np.mean(diff[inner] ** 2))),
        "map_max_abs_diff": float(np.abs(diff[inner]).max()),
        "map_mean_diff": float(diff[inner].mean()),
        "ssim": float(calc_ssim(sim_map, real_map, data_range=1.0)),
        "ncc": float(calc_ncc(sim_map, real_map)),
        "profile_rmse": float(np.sqrt(np.mean(prof_diff**2))),
        "profile_max_abs": float(np.abs(prof_diff).max()),
        "profile_bins": int(valid.sum()),
    }


# ---------------------------------------------------------------------------
# 図
# ---------------------------------------------------------------------------


def draw(sim_roi: np.ndarray, real_roi: np.ndarray, sim_m: dict, real_m: dict, agree: dict,
         label: str, output: Path) -> None:
    """実機・シミュ・差分のマップ，放射プロファイル，min の比較を 1 枚に描く."""
    fig, axes = plt.subplots(2, 3, figsize=(18, 11))
    sim_map, real_map = display_map(sim_roi), display_map(real_roi)
    vmin = 0.80
    for ax, data, title, m in ((axes[0, 0], real_map, "実機（DNG・LSC 適用・リサイズ）", real_m),
                               (axes[0, 1], sim_map, f"シミュ（{label}）", sim_m)):
        im = ax.imshow(data, cmap="hot", vmin=vmin, vmax=1.0)
        ax.set_title(f"{title}\nmin={m['min_fullimage']:.3f}"
                     f"（補正なし {m['min_fullimage_raw']:.3f}・"
                     f"輪帯 {m['min_ring']:.3f}）\n構造 CoV {m['structure_cov']:.4f}・"
                     f"ノイズ CoV {m['noise_cov']:.3f}", fontsize=9)
        ax.axis("off")
        fig.colorbar(im, ax=ax, fraction=0.046, label="輝度（99 パーセンタイル = 1）")
    ax = axes[0, 2]
    diff = sim_map - real_map
    im = ax.imshow(diff, cmap="coolwarm", vmin=-0.1, vmax=0.1)
    ax.set_title(f"差分（シミュ − 実機）\nRMS {agree['map_rms_diff']:.3f}・最大 "
                 f"{agree['map_max_abs_diff']:.3f}・SSIM {agree['ssim']:.3f}・"
                 f"NCC {agree['ncc']:.3f}", fontsize=9)
    ax.axis("off")
    fig.colorbar(im, ax=ax, fraction=0.046)

    ax = axes[1, 0]
    centers, sim_prof = radial_profile_norm(sim_roi)
    _, real_prof = radial_profile_norm(real_roi)
    ax.plot(centers, real_prof, marker="o", color="black", label="実機")
    ax.plot(centers, sim_prof, marker="s", color="tab:red", label=f"シミュ（{label}）")
    ax.axvline(MAX_DISTANCE, color="gray", lw=0.8, ls=":")
    ax.set_xlabel("正規化楕円距離（ビン中心）")
    ax.set_ylabel("輪帯平均 / ROI 平均")
    ax.set_title(f"20 分割 放射プロファイル（RMSE {agree['profile_rmse']:.4f}・最大差 "
                 f"{agree['profile_max_abs']:.4f}）", fontsize=10)
    ax.grid(alpha=0.3)
    ax.legend()
    ax = axes[1, 1]
    ax.bar(centers, sim_prof - real_prof, width=0.04, color="tab:red")
    ax.axhline(0, color="gray", lw=0.8)
    ax.set_xlabel("正規化楕円距離（ビン中心）")
    ax.set_ylabel("シミュ − 実機")
    ax.set_title("放射プロファイルの差（ビンごと）", fontsize=10)
    ax.grid(alpha=0.3)

    ax = axes[1, 2]
    keys = [("min_fullimage", "画像全体\n（補正）"),
            ("min_fullimage_raw", "画像全体\n（補正なし）"),
            ("min_ring", "輪帯平均"), ("min_pixel", "画素")]
    x = np.arange(len(keys))
    ax.bar(x - 0.18, [real_m[k] for k, _ in keys], width=0.36, color="black", label="実機")
    ax.bar(x + 0.18, [sim_m[k] for k, _ in keys], width=0.36, color="tab:red", label="シミュ")
    ax.axhline(ISO_THRESHOLD, color="tab:green", ls="--", label=f"基準 {ISO_THRESHOLD}")
    ax.set_xticks(x)
    ax.set_xticklabels([lbl for _, lbl in keys], fontsize=9)
    ax.set_ylim(0, 1.05)
    ax.set_ylabel("最小/最大比")
    ax.set_title("min の定義ごとの実機 vs シミュ", fontsize=10)
    for i, (k, _) in enumerate(keys):
        ax.text(i - 0.18, real_m[k] + 0.01, f"{real_m[k]:.3f}", ha="center", fontsize=8)
        ax.text(i + 0.18, sim_m[k] + 0.01, f"{sim_m[k]:.3f}", ha="center", fontsize=8)
    ax.legend(fontsize=8, loc="lower left")
    ax.grid(alpha=0.3, axis="y")
    fig.suptitle("実機 DNG vs シミュレーション — 同じノイズ処理での "
                 "min・放射プロファイル・形の一致", 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="実機 DNG とシミュの min・プロファイル・形の比較")
    parser.add_argument("--sim", type=Path, default=DEFAULT_SIM, help="シミュの luminance.npy")
    parser.add_argument("--label", default="ρ=0.82 純 diffuse", help="図のラベル")
    parser.add_argument("--tag", default="r082", help="出力ファイル名の識別子")
    args = parser.parse_args()

    real_full = load_real_roi()
    sim_roi = load_sim_roi(args.sim)
    real_roi = resize(real_full, sim_roi.shape, anti_aliasing=True, preserve_range=True)
    real_m, sim_m = metrics_for(real_roi), metrics_for(sim_roi)
    agree = shape_agreement(sim_roi, real_roi)
    # 実機は元解像度でも計算し，リサイズの影響を確認する（σ は画素比で拡大）
    scale = real_full.shape[0] / sim_roi.shape[0]
    native = calc_fullimage_min_max_ratio(real_full, sigma=STRUCTURE_SIGMA * scale)
    real_m["min_fullimage_native"] = native["fullimage_min_max_ratio"]

    OUT_DIR.mkdir(parents=True, exist_ok=True)
    draw(sim_roi, real_roi, sim_m, real_m, agree, args.label,
         OUT_DIR / f"real_vs_sim_{args.tag}.png")
    with (OUT_DIR / f"real_vs_sim_{args.tag}.csv").open("w", newline="", encoding="utf-8") as fh:
        writer = csv.writer(fh)
        writer.writerow(["metric", "real", "sim", "sim_minus_real"])
        for k in sim_m:
            writer.writerow([k, real_m.get(k, ""), sim_m[k],
                             sim_m[k] - real_m[k] if k in real_m else ""])
        writer.writerow(["min_fullimage_native", real_m["min_fullimage_native"], "", ""])
        for k, v in agree.items():
            writer.writerow([k, "", v, ""])

    print(f"\n=== 実機 vs シミュ（{args.label}）===")
    print(f"{'指標':>24} {'実機':>9} {'シミュ':>9} {'差':>9}")
    for k in ("mean", "pixel_cov", "noise_cov", "structure_cov", "min_fullimage",
              "min_fullimage_raw", "min_ring", "min_pixel", "a1", "a2", "a4", "b1"):
        print(f"{k:>24} {real_m[k]:9.4f} {sim_m[k]:9.4f} {sim_m[k] - real_m[k]:+9.4f}")
    print(f"{'min_fullimage_native':>24} {real_m['min_fullimage_native']:9.4f}   "
          f"（実機を元解像度・σ×{scale:.2f} で）")
    print("\n=== 形の一致 ===")
    for k, v in agree.items():
        print(f"{k:>24} {v:9.4f}" if isinstance(v, float) else f"{k:>24} {v:>9}")
    print(f"\n基準 {ISO_THRESHOLD}: 実機 min={real_m['min_fullimage']:.3f} → "
          f"{'合格' if real_m['min_fullimage'] >= ISO_THRESHOLD else '不合格'}，"
          f"シミュ min={sim_m['min_fullimage']:.3f} → "
          f"{'合格' if sim_m['min_fullimage'] >= ISO_THRESHOLD else '不合格'}")
    return 0


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