diff --git a/CLAUDE.md b/CLAUDE.md index d93c0a6..ef5f3b8 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -58,6 +58,7 @@ - `src/sm_tias_lightsim/remote/{commands,files,transfer}.py` + `src/sm_tias_lightsim/cli/remote.py` — `sm-remote` CLI（研究室サーバ Newton での実行）．push（tar over ssh．`.venv`/`output`/`.git`/`configs/_tmp_*` 除外）/ run（tmux でデタッチ起動．既定 `cuda_ad_rgb`・`--gpu`・`--threads`）/ status / logs / attach / kill / pull（`output/` 配下のみ．要求パス外・不正な名前・リンクを拒否）．サーバ上で実行すると止まって案内．使い方は ENV_04「研究室サーバでの実行」 - `src/sm_tias_lightsim/evaluation/dng_loader.py` — DNG 読み込みモジュール．SmTIAS-Evaluation `src/io/dng_loader.py` を SmTIAS-LightSim にコピーしたもの (オリジナルは SmTIAS-Evaluation 側で管理．バグ修正・機能追加はそちらで反映)．主要 API: `load_image_dng(path, demosaic, black_level, output_format)`．`output_format` は `rgb_float32` / `luma_float32` / `rgb_uint8`．linear 系出力 (`rgb_float32` / `luma_float32`) は `output_color=raw` でカラーマトリクス非適用・per-channel 比例性保持．`gamma=(1,1) + no_auto_bright=True + use_camera_wb=False + user_wb=[1,1,1,1]` で linear raw を保つ．依存: `rawpy>=0.21.0` - `src/sm_tias_lightsim/evaluation/radial.py` — 20 ビン放射状輝度プロファイル．既存 `calc_radial_profile` に加えて SmTIAS-Evaluation `src/analysis/spatial.py:calc_radial_min_max_ratio` をマージ．動径 min/max 比 (20 点輪帯平均) は画素単位 max/min よりノイズに対しロバスト．nan ビン除外を移植時に補強 (既存 `calc_radial_profile` は空ビンを nan で返す) +- `src/sm_tias_lightsim/evaluation/fullimage.py` — 画像全体で評価するノイズ補正付き最小/最大比（`calc_fullimage_min_max_ratio`）．σ=12 平滑化の min/max からノイズ場（残差の空間シャッフルを同じ平滑化）の極値を差し引く．動径 min/max 比が方位平均で非対称を見逃す問題（TEST_18）と，画素・パーセンタイルがノイズの裾を拾う問題の両方を避ける．`MetricsConfig.fullimage_min_max` / `fullimage_sigma` で有効化（既定 False）．偏りとノイズ床の注意は SPEC_04 - `src/sm_tias_lightsim/evaluation/shape_metrics.py` — 輝度分布「形」評価指標モジュール（新規）．現状 CoV/勾配/C/P/動径 min/max は不均一性の絶対量しか測れず，Y 軸単調勾配と中央ピーク楕円を構造的に区別できない問題を解決する 5 つの純関数: `calc_intensity_centroid`（輝度重心，5 値）/ `calc_angular_radial_decomposition`（極座標分解，angular_ratio が 0=同心円〜1=方向依存，3 値）/ `calc_ssim`（構造的類似度，scikit-image 依存，shape_reference 指定時のみ）/ `calc_ncc`（正規化相関係数，絶対輝度不変）/ `calc_intensity_pca`（weighted PCA で主軸角度と eccentricity，4 値）．mask 省略可能，完全均一画像で angular_ratio/NCC は NaN を返す．`MetricsConfig.shape_metrics=true` / `shape_reference` フィールドで有効化，後方互換保持．依存追加: `scikit-image>=0.21.0`．詳細: [TECH_11 形状指標の見直し](docs/05_TECH/TECH_11_形状指標の見直しと新指標導入.md) ※ 本欄は**最新ステップ 1 行のみ上書き更新**．詳細な進捗履歴（動機・設計判断・失敗パターン）は [docs/PROGRESS.md](docs/PROGRESS.md) に追記する．運用ルールは `.claude/rules/progress-log.md` を参照（該当ファイル編集時に自動ロードされる）． diff --git a/configs/eval_baseline_dng_target_shape.yaml b/configs/eval_baseline_dng_target_shape.yaml index 5cea129..e35cbdc 100644 --- a/configs/eval_baseline_dng_target_shape.yaml +++ b/configs/eval_baseline_dng_target_shape.yaml @@ -37,6 +37,7 @@ max_min: true radial_profile: true radial_min_max: true + fullimage_min_max: true # 画像全体のノイズ補正付き min/max（方位依存の不均一も測る） shape_metrics: true shape_reference: null # 実機 DNG 由来の参照輝度マップを指定すると SSIM/NCC が有効化される normalization: diff --git "a/docs/04_SPEC/SPEC_01_\343\202\242\343\203\274\343\202\255\343\203\206\343\202\257\343\203\201\343\203\243\350\250\255\350\250\210.md" "b/docs/04_SPEC/SPEC_01_\343\202\242\343\203\274\343\202\255\343\203\206\343\202\257\343\203\201\343\203\243\350\250\255\350\250\210.md" index ee3c6cc..b54e6e0 100644 --- "a/docs/04_SPEC/SPEC_01_\343\202\242\343\203\274\343\202\255\343\203\206\343\202\257\343\203\201\343\203\243\350\250\255\350\250\210.md" +++ "b/docs/04_SPEC/SPEC_01_\343\202\242\343\203\274\343\202\255\343\203\206\343\202\257\343\203\201\343\203\243\350\250\255\350\250\210.md" @@ -124,6 +124,7 @@ | `evaluation/normalize` | linear 輝度 → mean を target_mean に揃えた輝度 | 線形ゲイン．形状指標は不変．OETF と排他 | | `evaluation/dng_loader` | DNG（RAW）→ linear raw 配列 | 実測リファレンスの linear 取得 | | `evaluation/shape_metrics` | 輝度配列（+ 任意の参照画像）→ 形状指標 | 輝度重心 / 角度-放射分解 / PCA / SSIM / NCC | +| `evaluation/fullimage` | 輝度配列 → 画像全体のノイズ補正付き min/max 比 | 方位に依らない．σ 平滑化 + ノイズ場の極値で補正．純関数 | | `evaluation/visualize` | 各種統計 → matplotlib 図 PNG | 可視化のみ．データ加工はしない | | `evaluation/pipeline` | 画像 + 評価 YAML + 出力先 → 指標 dict（metrics.csv を書く） | `save_figures=False` で図を省略．段階別の `EvaluationStepError` を投げる | | `visualize/baker` | シーン + 光源 → 各 mesh face 中心の照度（W/m²） | direct（直接光）/ path（多重反射込み）モード | diff --git "a/docs/04_SPEC/SPEC_02_\343\202\267\343\203\274\343\203\263\345\256\232\347\276\251\344\273\225\346\247\230.md" "b/docs/04_SPEC/SPEC_02_\343\202\267\343\203\274\343\203\263\345\256\232\347\276\251\344\273\225\346\247\230.md" index b28e727..12ebb4d 100644 --- "a/docs/04_SPEC/SPEC_02_\343\202\267\343\203\274\343\203\263\345\256\232\347\276\251\344\273\225\346\247\230.md" +++ "b/docs/04_SPEC/SPEC_02_\343\202\267\343\203\274\343\203\263\345\256\232\347\276\251\344\273\225\346\247\230.md" @@ -212,6 +212,8 @@ max_min: true radial_profile: true radial_min_max: false # 動径 min/max 比の 4 列を metrics.csv に追加（既定 false） + fullimage_min_max: false # 画像全体のノイズ補正付き min/max 比の 5 列を追加（既定 false） + fullimage_sigma: 12.0 # 上記の平滑化 σ [px]（> 0．ROI の短辺が 4σ 以上必要） shape_metrics: false # 形状指標（輝度重心 / 角度-放射分解 / PCA 主軸）を追加（既定 false） shape_reference: null # SSIM / NCC 算出用の参照画像パス（.npy / .png）．shape_metrics=true 時のみ有効 oetf: # 任意．省略時 OETF 無効（linear のまま評価＝後方互換） diff --git "a/docs/04_SPEC/SPEC_04_\350\251\225\344\276\241\344\273\225\346\247\230.md" "b/docs/04_SPEC/SPEC_04_\350\251\225\344\276\241\344\273\225\346\247\230.md" index 5518beb..9a7155f 100644 --- "a/docs/04_SPEC/SPEC_04_\350\251\225\344\276\241\344\273\225\346\247\230.md" +++ "b/docs/04_SPEC/SPEC_04_\350\251\225\344\276\241\344\273\225\346\247\230.md" @@ -246,6 +246,47 @@ `metrics.csv` に `radial_min_max_ratio` / `radial_min_distance` / `radial_min` / `radial_max` の 4 列を末尾に追加．既定 False（既存挙動互換）．`metrics.radial_profile` の有無に関わらず内部で動径プロファイルを算出するため独立して有効化可能．実装: `evaluation/radial.py:calc_radial_min_max_ratio`（SmTIAS-Evaluation の `src/analysis/spatial.py` からマージ）． +### 画像全体 min/max 比（オプション，2026-09-14 導入） + +`metrics.fullimage_min_max: true` で有効化．**画像全体で評価する，ノイズ補正付きの最小/最大比**． +動径 min/max 比は方位方向に平均するため，同心円状でない輝度分布（オフ軸 LED 配置など）では +明るい側と暗い側が同じ輪帯で混ざり非対称を過小評価する（TEST_18 では非対称が強い条件で +0.97〜0.99 に張り付き，差は非対称の大きさ \|a1\| と相関 0.93〜0.95）．一方，画素単位や平滑化後の +パーセンタイルは MC ノイズの裾を拾い，生の輝度が低い（相対ノイズが大きい）条件で頭打ちになる． +本指標はその両方を避ける． + +定義（`evaluation/fullimage.py:calc_fullimage_min_max_ratio`）: + +1. ROI 輝度を σ = `fullimage_sigma`（既定 12 px）のガウス平滑化して `smooth` を得る． +2. ノイズ場を推定する: 残差（ROI − ガウス(σ=4)(メディアン 5 px)(ROI)）を空間シャッフルして同じ + σ で平滑化し，その最小値 `noise_min`・最大値 `noise_max` を 4 回の平均で求める（構造を壊して + ノイズの統計だけ残す．構造 CoV と同じ発想）． +3. `fullimage_min = min(smooth) − noise_min`，`fullimage_max = max(smooth) − noise_max`， + `fullimage_min_max_ratio = fullimage_min / fullimage_max`． + +| 指標 | 定義 | 値域 | +| --- | --- | --- | +| **画像全体 min/max 比** | ノイズ補正後の `min / max` | 0〜1（1.0 で均一．小さいほどムラ大） | + +`metrics.csv` に `fullimage_min_max_ratio` / `fullimage_min` / `fullimage_max` / `fullimage_noise_min` / +`fullimage_noise_max` の 5 列を追加（動径 min/max の後・形状指標の前）．既定 False（既存挙動互換）． +線形ゲインに不変（mean anchor 正規化の前後で同じ値）．同じ画像・同じ設定なら常に同じ値（乱数は +関数内で seed 固定）．入力が 2D でない・σ ≤ 0・`min(H, W) < 4σ`・補正後の max ≤ 0 は ValueError． + +検証（`scripts/build_fullimage_minmax.py`，軸上スイープ 20 条件）: 対称な分布では動径 min/max 比と +±0.01 以内で一致し，一様像＋同量のノイズで測ったノイズ床は ±0.004．ブロック平均・多項式フィット +・補正なし平滑化の各候補より偏りとノイズ床の両方で優れた． + +**使用上の注意**: + +- **平滑化による偏り**: σ の平滑化が構造をなますため，比は真値より少し大きく（均一寄りに）出る + （256² の線形勾配で真値 0.800 → 0.813，ガウス山で 0.857 → 0.865）．偏りは σ と ROI サイズの比で + 決まるので，**σ を固定して条件間を比較する**用途に使い，絶対値を他の定義と直接比べない． +- **ノイズ床はノイズ量に比例**: 一様像で相対ノイズ 5% → 0.003，10% → 0.006，20% → 0.012． + 構造のある像ではノイズが増えると比が均一寄りに上振れするので，**同じ spp（同じノイズ水準） + の条件同士で比較する**． +- ROI が小さいときは σ を下げる（`fullimage_sigma`）．4σ 未満の辺は受け付けない． + ### 形状指標 (Shape Metrics — オプション、2026-06-03 導入) **動機**: 既存の 4 指標 (CoV/勾配/C/P/動径 min/max) は **不均一性の絶対量しか測れず、「形」を構造的に区別できない**．シミュ「Y 軸単調勾配」と実機「中央ピーク楕円」のように形が違う 2 つの画像でも、数値だけは似通うことがある．詳細は [TECH_11 形状指標の見直しと新指標導入](../05_TECH/TECH_11_形状指標の見直しと新指標導入.md)． diff --git a/docs/PROGRESS.md b/docs/PROGRESS.md index 848bef0..39c2d5c 100644 --- a/docs/PROGRESS.md +++ b/docs/PROGRESS.md @@ -102,3 +102,5 @@ - **スイープの高速化 — レンダ/評価の関数化・プロセス内実行・cuda 遮蔽検証・2 並列 (2026-09-14)**: ユーザー指示「効率化を先に」を受け，`scripts/profile_sweep_step.py` で 1 条件の内訳を実測してから手を入れた．**計測結果**: GPU のレンダ本体（spp=1024・512²）が 4.4〜5.2 秒で床．その上に `uv run sm-render`/`sm-evaluate` の起動・import・YAML・PLY 書き出しで約 2.5 秒，図の描画で約 2 秒，さらに遮蔽検証を YAML の llvm で回していたため球で約 25 秒（頂点 21 万個のメッシュを CPU で読み直す）が乗り，球 1 条件 40 秒・立方体 13 秒だった．**対策**: ① `sm-evaluate`/`sm-render` の本体を `evaluation/pipeline.py:evaluate_image`・`render/pipeline.py:render_yaml` に関数化し CLI は引数パースと例外変換だけに（CLI 挙動は不変）/ ② `--no-figures`（`save_figures=False`）/ ③ `runner.ensure_variant` を切り出し，遮蔽検証の前に cuda を確定（25 秒 → 0.5 秒）/ ④ スイープはサブプロセスをやめプロセス内で呼ぶ / ⑤ `--workers N`（spawn のプロセスプール）．**実測**: 立方体 13 → 5.5 秒，球 40 → 約 6 秒．`--workers` は GPU 律速のため 1→5.5・2→5.0・3→4.5〜5.8 秒でほぼ変わらず，**2 を推奨**（CPU 側処理を隠す分だけ）．残り 75 条件が 8 分で終わった．**発見**: 同一 seed の再レンダは集計指標が 4 桁まで一致し（画素はノイズ程度に変わるが集計で消える）ノイズフロアにならない → `RenderConfig.seed` を追加し `mi.render` に渡す（再現性条件は rep ごとに seed を変えた YAML）．`/implement` の 4 フェーズ（coder → tester 101 件追加 → refactorer → docs）を，ユーザー就寝中の指示に基づき Phase 1 の動作確認をサーバ実走で代行して完走．GPU は 0 番のみ使用（ユーザー指示）． - **TEST_18 照明配置依存性の球と立方体の比較 — 配置依存は球が立方体の 1/2〜1/3 (2026-09-14)**: ユーザーの仮説「積分球が使われる理由は照明配置への依存が小さいから」を，同内面積 4 ペア（球 φ200/270/340/500（r=100/135/170/250）↔ 立方体 L=141/187/234/340）で **LED の位置（キャビティ中心から見た極角 θ=0〜120° を 13 段階）・向き（板向き＋バッフル／壁法線向き＋バッフル／真下向き）・個数（θ=90° の輪に N=1〜4）・再現性（seed 違い）** の 198 条件で検証した（Newton cuda・1 条件 5〜6 秒）．**結果 ①**: θ 全域での構造 CoV の幅は球 0.015〜0.023・立方体 0.025〜0.071（1.2〜3.1 倍），max\|a1\| は 3〜4% vs 4.5〜11%．立方体は θ とともに単調に悪化し板の直下で軸上の 3.4〜5 倍．**球は θ≈30° が軸上より均一**（φ340: 0.0193 → 0.0057．軸上のリング構造が崩れる）で「どこに置いてもそこそこ」．立方体の優位は軸上の 1 点だけ．**結果 ②**: 真下向き（板が背後＝バッフル無し）は板向き＋バッフルと同等で 5〜30% 明るい．壁法線向き（キャビティ中心向き）は最悪（2〜5 倍）．**結果 ③**: 対向 2 灯で m=1 は ≤0.3% に消え構造 CoV は 0.005〜0.010＝単灯軸上より均一で形状差が消える．立方体は N=3 で崩れる（3 回対称が 4 回対称の壁と噛み合わない）．**ノイズフロア**: 構造 CoV ±0.0007・a1 ±0.08%．設計への含意は「単灯なら球，多灯にできるなら側壁の対向 2〜4 灯・真下向き・バッフル無しで立方体でも球と同等」．**方法論**: 極角 θ を両形状共通のパラメータにし，立方体の稜（45°）は外す．板の近くのバッフルは外殻の貫通を許容し食い込み量を記録（壁法線向きは半径がキャビティ寸法の 0.6 を超えたら無効）．解析は TEST_16/17 の方位フーリエ（m=1〜4）＋構造 CoV＋重心．**追加実験（同日夜）**: ① 方位角: 立方体は面の中央が最良で稜へ寄せるほど構造 CoV が 1.3〜1.75 倍に悪化（球は方位に依らない）/ ② 真下向き多灯（バッフル無し）: 2 灯以上で非対称は消え単灯軸上より均一だが，板向き＋バッフルの多灯には 1.5〜2.5 倍劣り N を増やしても頭打ち（一次照射が壁の下半分に集中し動径構造が残る）．③ 校正材質（純 diffuse ρ=0.81）で主要 60 条件を引き直し: 配置依存の幅は球 0.030 vs 立方体 0.077〜0.098 と結論は不変（むしろ差が広がる）．ただし旧材質の光沢 5% が作っていた軸上のリング構造が消え，**校正材質では両形状とも軸上が最良**（「球は θ≈30° が最良」は旧材質の産物）．全 255 + 60 条件（TEST_18 結果 5〜7）．**失敗パターン**: spawn のワーカーには主プロセスで切り替えた出力ディレクトリ（材質）が伝わらず，旧材質の row.json を再利用して 60 条件が空振りした → 材質を run_job の引数で渡す． + +- **正規化と最小/最大比の定義の見直し — 画像全体・ノイズ補正付きの min/max 比を指標に追加 (2026-09-14)**: ユーザー指摘「輪帯ピーク正規化は同心円状でない分布に向かない」から，輝度マップの正規化と min/max 比の定義を軸上スイープ（TEST_15 球 10・TEST_16 立方体 10）で系統的に比べた（`scripts/build_normalization_comparison.py`・`build_minmax_definition_comparison.py`・`build_fullimage_minmax.py`，出力 `output/_normalization_comparison/`）．**発見**: ① 表示の正規化は対称な分布なら輪帯ピークとロバスト最大（平滑化後 99 パーセンタイル）が 1% 以内で同じ．非対称では輪帯ピークが破綻するので表示は後者に統一 / ② min/max 比は「最小」側がノイズの裾に敏感で，平滑化後のパーセンタイルは大きいキャビティ（生輝度が低くノイズ 26%）で 0.91〜0.92 に頭打ち．輪帯平均は数千画素の平均でノイズに勝つが方位を平均する / ③ 画像全体で使える定義 7 候補のうち「σ=12 平滑化の min/max をノイズ場（残差の空間シャッフル）の極値で補正」が対称分布で輪帯平均と ±0.01 で一致しノイズ床 ±0.004 と最良．ブロック平均は粗さと床のトレードオフ，多項式フィットは形を追えない．**実装**: `/implement` で `evaluation/fullimage.py` と `MetricsConfig.fullimage_min_max` を追加（既定 False．TEST_18 の評価 YAML で有効化）．TEST_18 に適用すると輪帯平均は非対称条件で 0.97〜0.99 に張り付き，新指標は球 0.87〜0.91・立方体 0.65 まで下がる（差と \|a1\| の相関 0.93〜0.95）．TEST_18 の全結果表と輝度マップパネルに min（最大=1 のときの最小値）を併記．**注意（テスターの指摘）**: 平滑化で真値より均一寄りに出る（線形勾配で +0.013）・ノイズ床はノイズ量に比例（相対 20% で 0.012）ため，σ と spp を揃えた条件間比較に使う．球のサイズ表記は内径 φ に統一（TEST_18 の文書・図．条件名は半径のまま）． diff --git a/scripts/build_led_placement_comparison.py b/scripts/build_led_placement_comparison.py index ec57ea1..0afc739 100644 --- a/scripts/build_led_placement_comparison.py +++ b/scripts/build_led_placement_comparison.py @@ -123,7 +123,8 @@ "occlusion_blocked_frac", "raw_mean_roi", "cov")} result["pair"] = pair_index(row["shape"], row["size_mm"]) result["baffle_ratio"] = row["baffle_radius_mm"] / row["size_mm"] - rings = ring_harmonics(load_roi(lum), HARMONICS) + roi = load_roi(lum) + rings = ring_harmonics(roi, HARMONICS) for m in HARMONICS: result[f"a{m}"], result[f"se{m}"] = aggregate(rings, f"a{m}") result[f"b{m}"], _ = aggregate(rings, f"b{m}") @@ -141,8 +142,7 @@ result["noise_cov"] = noise result["ring_min_max"] = ring_min_max_ratio(centers, means) # 輪帯平均は方位方向に平均するのでオフ軸の非対称を打ち消す．画像全体の定義を併記する - result["fullimage_min_max"] = calc_fullimage_min_max_ratio( - load_roi(lum))["fullimage_min_max_ratio"] + result["fullimage_min_max"] = calc_fullimage_min_max_ratio(roi)["fullimage_min_max_ratio"] result["centroid_x"] = row.get("intensity_centroid_x_norm", float("nan")) result["centroid_y"] = row.get("intensity_centroid_y_norm", float("nan")) result["centroid_offset"] = row.get("intensity_centroid_offset_norm", float("nan")) diff --git a/src/sm_tias_lightsim/cli/evaluate.py b/src/sm_tias_lightsim/cli/evaluate.py index c642f35..0ac07eb 100644 --- a/src/sm_tias_lightsim/cli/evaluate.py +++ b/src/sm_tias_lightsim/cli/evaluate.py @@ -40,6 +40,7 @@ _load_luminance, _load_reference_luminance, _print_summary, + _run_fullimage_min_max, _run_radial_min_max, _run_radial_profile, _run_shape_metrics, diff --git a/src/sm_tias_lightsim/evaluation/fullimage.py b/src/sm_tias_lightsim/evaluation/fullimage.py new file mode 100644 index 0000000..fae4894 --- /dev/null +++ b/src/sm_tias_lightsim/evaluation/fullimage.py @@ -0,0 +1,174 @@ +"""画像全体で評価する，ノイズ補正付き最小/最大比のモジュール． + +輪帯平均の min/max（``radial.calc_radial_min_max_ratio``）は方位方向に平均するため， +同心円状でない輝度分布（TEST_18 のオフ軸 LED 配置）では不均一を平均で打ち消してしまう． +一方，画素単位の min/max や平滑化後のパーセンタイルはモンテカルロ（MC）ノイズの裾を拾い， +生の輝度が低い大きなキャビティで値が頭打ちになる． + +本モジュールはその中間として「ガウス平滑化した像の min/max を，同じ平滑化を掛けた +ノイズ場の極値で補正する」定義を提供する．平滑化で方位方向に平均せずノイズだけを落とし， +残ったノイズの寄与（平滑化後でも残る極値のバイアス）を統計的に差し引く． + +補正に使うノイズ場は，残差（生 − メディアン/ガウスで均した像）を空間シャッフルして +同じ σ で平滑化したもの．シャッフルで構造を壊し，ノイズの振幅分布だけを残す +（``evaluation`` 外の ``scripts/build_radius_sweep_visual.py:decompose_cov`` と同じ発想）． + +仮定: + ノイズが空間的に独立（画素間で無相関）であること．シャッフルは画素の並べ替えなので， + もともと空間相関を持つノイズ（デノイザ適用後の像等）では補正量を過小評価する． + MC レンダの生出力（デノイズ無し）を前提とする． + +検証: + 軸上スイープ（TEST_15 球 / TEST_16 立方体の計 20 条件，分布が同心円対称なので + 輪帯平均が基準として使える）で，本定義は輪帯平均と ±0.01 以内で一致し， + 一様像にノイズだけを載せたときの見かけの不均一（ノイズ床）は ±0.004 だった． + 検証スクリプト: ``scripts/build_fullimage_minmax.py``（定義名 ``smooth_corr``）． +""" + +from __future__ import annotations + +import numpy +from scipy.ndimage import gaussian_filter, median_filter + + +def estimate_noise_extremes( + image: numpy.ndarray, + sigma: float, + residual_sigma: float = 4.0, + median_size: int = 5, + n_shuffle: int = 4, + seed: int = 0, +) -> tuple[float, float]: + """ノイズだけの場に σ=``sigma`` の平滑化を掛けたときの (min, max) の期待値を返す． + + 手順: + 1. 残差 $r = I - G_{\\sigma_r}(M_k(I))$ を取る（$M_k$ = サイズ $k$ のメディアン + フィルタ，$G_{\\sigma_r}$ = σ=``residual_sigma`` のガウス平滑化）． + メディアンで firefly を落としてからガウスで均すことで，残差に構造を残さない． + 2. 残差を空間シャッフルして構造を壊し，σ=``sigma`` のガウス平滑化を掛ける． + 3. その場の min / max を ``n_shuffle`` 回平均する． + + 乱数は ``numpy.random.default_rng(seed)`` をこの関数内で生成するため， + 同じ入力・同じ ``seed`` なら常に同じ値を返す（呼び出し順に依存しない）． + + Args: + image: 2D 輝度配列（H x W）． + sigma: 本体の平滑化と同じガウス σ（画素単位）． + residual_sigma: 残差を取るときのガウス σ（画素単位）． + median_size: 残差を取るときのメディアンフィルタのサイズ（画素）． + n_shuffle: シャッフルの反復回数（多いほど期待値の推定が安定する）． + seed: 乱数シード． + + Returns: + (ノイズ場の最小値の期待値, ノイズ場の最大値の期待値)． + それぞれ平滑化後の像の min / max から差し引くための量（min 側は負，max 側は正）． + """ + img = numpy.asarray(image, dtype=numpy.float64) + rng = numpy.random.default_rng(seed) + + residual = img - gaussian_filter(median_filter(img, size=median_size), sigma=residual_sigma) + + mins: list[float] = [] + maxs: list[float] = [] + for _ in range(n_shuffle): + flat = residual.ravel().copy() + rng.shuffle(flat) + field = gaussian_filter(flat.reshape(img.shape), sigma=sigma) + mins.append(float(field.min())) + maxs.append(float(field.max())) + + return float(numpy.mean(mins)), float(numpy.mean(maxs)) + + +def calc_fullimage_min_max_ratio( + image: numpy.ndarray, + sigma: float = 12.0, + residual_sigma: float = 4.0, + median_size: int = 5, + n_shuffle: int = 4, + seed: int = 0, +) -> dict[str, float]: + """画像全体で評価する，ノイズ補正付き最小/最大比を算出する． + + $\\mathrm{ratio} = \\dfrac{\\min G_\\sigma(I) - n_{\\min}}{\\max G_\\sigma(I) - n_{\\max}}$ + + ここで $G_\\sigma$ は σ=``sigma`` のガウス平滑化，$(n_{\\min}, n_{\\max})$ は + ``estimate_noise_extremes`` が返すノイズ場の極値の期待値． + 平滑化で MC ノイズを落としつつ方位方向には平均しないので， + オフ軸 LED のような方位依存の不均一も min/max に反映される． + + 検証（``scripts/build_fullimage_minmax.py``）: 同心円対称な軸上スイープ 20 条件で + 輪帯平均の min/max と ±0.01 以内で一致し，ノイズ床（一様像 + ノイズでの 1 からのずれ）は + ±0.004．ノイズが空間的に独立であることを仮定する（モジュール docstring 参照）． + + Args: + image: 2D 輝度配列（H x W）．ROI 抽出後の輝度を渡すこと． + sigma: ガウス平滑化の σ（画素単位）．既定 12.0 は 464x450 の ROI で検証した値． + residual_sigma: ノイズ場推定の残差を取るガウス σ（画素単位）． + median_size: ノイズ場推定の残差を取るメディアンフィルタのサイズ（画素）． + n_shuffle: ノイズ場の極値を平均する回数． + seed: ノイズ場推定の乱数シード（同じ入力・同じ seed なら決定論的）． + + Returns: + 以下のキーを持つ辞書: + - ``"fullimage_min_max_ratio"`` (float): 補正後の最小/最大比（1 が理想）． + - ``"fullimage_min"`` (float): 補正後の最小値（平滑化後 min − ノイズ場 min）． + - ``"fullimage_max"`` (float): 補正後の最大値（平滑化後 max − ノイズ場 max）． + - ``"fullimage_noise_min"`` (float): ノイズ場の最小値の期待値． + - ``"fullimage_noise_max"`` (float): ノイズ場の最大値の期待値． + + Raises: + ValueError: image が 2D の数値配列でない場合． + ValueError: sigma が 0 以下の場合． + ValueError: 画像が平滑化に対して小さすぎる場合（min(H, W) < 4 * sigma）． + ValueError: 補正後の最大値が 0 以下の場合（画像がノイズだけとみなせる）． + """ + img = numpy.asarray(image) + if img.ndim != 2: + msg = f"image は 2D 配列を期待します（ndim={img.ndim}）" + raise ValueError(msg) + if not numpy.issubdtype(img.dtype, numpy.number): + msg = f"image は数値配列を期待します（dtype={img.dtype}）" + raise ValueError(msg) + img = img.astype(numpy.float64) + + if sigma <= 0.0: + msg = f"sigma は正の値を期待します（sigma={sigma}）" + raise ValueError(msg) + + # 平滑化カーネルが画像からはみ出すと端の値が境界条件に支配されるため下限を設ける + min_side = min(img.shape) + if min_side < 4.0 * sigma: + msg = ( + f"画像が平滑化 σ={sigma:g} に対して小さすぎます" + f"（min(H, W)={min_side} < 4σ={4.0 * sigma:g}）．sigma を小さくしてください．" + ) + raise ValueError(msg) + + smooth = gaussian_filter(img, sigma=sigma) + noise_min, noise_max = estimate_noise_extremes( + img, + sigma=sigma, + residual_sigma=residual_sigma, + median_size=median_size, + n_shuffle=n_shuffle, + seed=seed, + ) + + corrected_min = float(smooth.min()) - noise_min + corrected_max = float(smooth.max()) - noise_max + if corrected_max <= 0.0: + msg = ( + "ノイズ補正後の最大値が 0 以下です" + f"（平滑化後 max={smooth.max():.6g}, ノイズ場 max={noise_max:.6g}）．" + "画像が実質ノイズのみの可能性があります．" + ) + raise ValueError(msg) + + return { + "fullimage_min_max_ratio": corrected_min / corrected_max, + "fullimage_min": corrected_min, + "fullimage_max": corrected_max, + "fullimage_noise_min": noise_min, + "fullimage_noise_max": noise_max, + } diff --git a/src/sm_tias_lightsim/evaluation/pipeline.py b/src/sm_tias_lightsim/evaluation/pipeline.py index f8e89fa..21fa584 100644 --- a/src/sm_tias_lightsim/evaluation/pipeline.py +++ b/src/sm_tias_lightsim/evaluation/pipeline.py @@ -28,6 +28,7 @@ import pandas # noqa: E402 from PIL import Image # noqa: E402 +from sm_tias_lightsim.evaluation.fullimage import calc_fullimage_min_max_ratio # noqa: E402 from sm_tias_lightsim.evaluation.normalize import apply_mean_anchor # noqa: E402 from sm_tias_lightsim.evaluation.oetf import apply_srgb_oetf, solve_exposure # noqa: E402 from sm_tias_lightsim.evaluation.radial import ( # noqa: E402 @@ -321,6 +322,31 @@ return calc_radial_min_max_ratio(profile) +def _run_fullimage_min_max( + roi_luminance: numpy.ndarray, + metrics_config: MetricsConfig | None, +) -> dict: + """画像全体のノイズ補正付き min/max 比を算出して返す． + + 輪帯平均の min/max（`_run_radial_min_max`）は方位方向に平均するため， + オフ軸 LED のような方位依存の不均一を打ち消してしまう．本指標は平滑化で + MC ノイズだけを落とし，残るノイズの寄与を差し引いた min/max を返す． + + Args: + roi_luminance: ROI 輝度配列（H x W, float64）． + metrics_config: MetricsConfig インスタンス（fullimage_min_max フラグと + fullimage_sigma を持つ）．None の場合は無効扱い． + + Returns: + 画像全体 min/max 比指標の辞書（fullimage_min_max_ratio / fullimage_min / + fullimage_max / fullimage_noise_min / fullimage_noise_max の 5 キー）．無効時は空辞書． + """ + if metrics_config is None or not metrics_config.fullimage_min_max: + return {} + + return calc_fullimage_min_max_ratio(roi_luminance, sigma=metrics_config.fullimage_sigma) + + def _load_reference_luminance(reference_path: Path) -> numpy.ndarray: """形状指標用の参照画像を読み込み，輝度配列（2D float64）を返す． @@ -512,6 +538,9 @@ print(f" radial min/max : {metrics['radial_min_max_ratio']:.4f}") print(f" radial min dist: {metrics['radial_min_distance']:.4f}") + if "fullimage_min_max_ratio" in metrics: + print(f" fullimage min/max: {metrics['fullimage_min_max_ratio']:.4f}") + if "intensity_centroid_offset_norm" in metrics: print(f" centroid offset: {metrics['intensity_centroid_offset_norm']:.4f} (norm)") print( @@ -555,10 +584,11 @@ 8. zones.enabled のときゾーン分析を実行する． 9. metrics.radial_profile のとき放射状プロファイルを算出する． 10. metrics.radial_min_max のとき動径 min/max 比を算出する（radial_profile 再利用）． - 11. metrics.shape_metrics のとき形状指標を算出する． - 12. metrics.csv を保存する． - 13. save_figures のとき可視化 PNG を figures/ に保存する． - 14. 算出指標のサマリーを標準出力に表示する． + 11. metrics.fullimage_min_max のとき画像全体のノイズ補正付き min/max 比を算出する． + 12. metrics.shape_metrics のとき形状指標を算出する． + 13. metrics.csv を保存する． + 14. save_figures のとき可視化 PNG を figures/ に保存する． + 15. 算出指標のサマリーを標準出力に表示する． Args: image_path: 評価対象の画像ファイル（PNG または .npy）． @@ -622,6 +652,9 @@ _run_radial_min_max(roi_luminance, eval_config.metrics, radial_profile) ) + with _step("画像全体 min/max 比の算出に失敗しました"): + metrics.update(_run_fullimage_min_max(roi_luminance, eval_config.metrics)) + with _step("形状指標の算出に失敗しました"): metrics.update( _run_shape_metrics(roi_luminance, eval_config.metrics, config_path.parent) diff --git a/src/sm_tias_lightsim/scene/schema.py b/src/sm_tias_lightsim/scene/schema.py index ec9996f..1ff8078 100644 --- a/src/sm_tias_lightsim/scene/schema.py +++ b/src/sm_tias_lightsim/scene/schema.py @@ -452,6 +452,23 @@ "True のとき radial_profile の有無に関わらず内部で radial_profile を算出する．" ), ) + fullimage_min_max: bool = Field( + default=False, + description=( + "画像全体のノイズ補正付き min/max 比（fullimage_min_max_ratio / fullimage_min / " + "fullimage_max / fullimage_noise_min / fullimage_noise_max の 5 列を " + "metrics.csv に追加）．輪帯平均と違い方位方向に平均しないため，" + "オフ軸 LED のような方位依存の不均一も評価できる．" + ), + ) + fullimage_sigma: float = Field( + default=12.0, + gt=0, + description=( + "fullimage_min_max のガウス平滑化 σ（画素単位）．" + "既定 12.0 は 464x450 の ROI で検証した値．" + ), + ) shape_metrics: bool = Field( default=False, description=(