diff --git a/scripts/compare_real_vs_sim_metrics.py b/scripts/compare_real_vs_sim_metrics.py new file mode 100644 index 0000000..595459e --- /dev/null +++ b/scripts/compare_real_vs_sim_metrics.py @@ -0,0 +1,327 @@ +"""実機 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}（補正なし {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}・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} （実機を元解像度・σ×{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())