diff --git a/CLAUDE.md b/CLAUDE.md index 47f6d2b..a9d43e8 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -23,7 +23,7 @@ ## 開発進捗 -最新: フェーズ 7（実装）— **TEST_18 照明配置依存性の球と立方体の比較**（LED の位置 θ・向き・個数を同内面積ペアで振った 255 条件．配置依存は球が立方体の 1/2〜1/3，対向 2 灯以上で形状差が消える．校正材質 ρ=0.81 でも結論不変（軸上は両形状とも最良になる）．真下向き多灯はバッフル無しで成立するが板向き多灯に劣る）と，**スイープの高速化**（`sm-render`/`sm-evaluate` の本体を `render/pipeline`・`evaluation/pipeline` に関数化してプロセス内実行，`--no-figures`，遮蔽検証を cuda で，`--workers 2`．1 条件 球 40 秒 → 6 秒）．サーバは `sm-remote`（学外からは `SM_REMOTE_HOST=Newton-via-lab`）．**環境**: mitsuba は **3.7.1 に固定**（テストは `uv run python -m pytest`），開発機の **SAC はオフ**（ENV_04 参照）．同一 seed の再レンダは集計指標が一致するので，ノイズフロアは `render.seed` を変えて測る．**spp の指針**（[TEST_19](docs/06_TEST/TEST_19_sppと計算時間ノイズ指標安定性.md)）: GPU 時間は spp に比例（1024 で 5 秒）・固定費 0.8 秒・ノイズは 1/√spp 則．構造 CoV なら 1024，画像全体 min/max 比を暗い条件で比べるなら 2048〜4096．詳細は [docs/PROGRESS.md](docs/PROGRESS.md)． +最新: フェーズ 7（実装）— **TEST_20 実機とシミュの形の突き合わせと配光の検討**（実機 DNG 全灯 + 単灯 3 灯とシミュを同じノイズ処理で比較．同心円方向は ρ≈0.81 で合うが，単灯で明部が反対側へ寄る量と「左上が暗い」傾きはシミュが実機の半分．壁の光沢 ω・α では形全体を崩さずに合わせられず，LED 配光を spot FWHM 45〜60° に狭めると 4 構成の差分 RMS が最良．新しい基準候補は採用材質 ρ=0.85・ω=0.05・α=0.1 + spot 45〜60°．共通の傾きは基板の位置・向きが第一候補で未解決）．前ステップ: TEST_18 照明配置依存性（配置依存は球が立方体の 1/2〜1/3），TEST_19 spp と計算時間（GPU 時間 ∝ spp，構造 CoV は 1024 で十分）．サーバは `sm-remote`（学外からは `SM_REMOTE_HOST=Newton-via-lab`）．**環境**: mitsuba は **3.7.1 に固定**（テストは `uv run python -m pytest`），開発機の **SAC はオフ**（ENV_04 参照）．同一 seed の再レンダは集計指標が一致するので，ノイズフロアは `render.seed` を変えて測る．**spp の指針**（[TEST_19](docs/06_TEST/TEST_19_sppと計算時間ノイズ指標安定性.md)）: GPU 時間は spp に比例（1024 で 5 秒）・固定費 0.8 秒・ノイズは 1/√spp 則．構造 CoV なら 1024，画像全体 min/max 比を暗い条件で比べるなら 2048〜4096．詳細は [docs/PROGRESS.md](docs/PROGRESS.md)． 実装済みモジュール: - `src/sm_tias_lightsim/uniformity.py` — Rec.709 輝度変換 + CoV/std/max-min（SmTIAS-Evaluation コピー） @@ -157,4 +157,5 @@ - 立方体キャビティ導入と球との比較: docs/06_TEST/TEST_16_立方体キャビティ導入と球との比較.md（ロフト・接合・喉の無い立方体を L=120〜360mm でスイープし内面積で球と対応づけ．同面積なら立方体が一貫して均一で，四隅の写り込みは m=4 方位フーリエで検出限界以下と確認） - オフ軸 LED 配置の影響とサイズ依存: docs/06_TEST/TEST_17_オフ軸LED配置の影響とサイズ依存.md（LED を側壁に置いた非対称の定量．サイズとともに減衰するが，頑健性は球 > 立方体で軸上の優劣と逆転．明るくなるのは LED の反対側） - spp と計算時間・ノイズ・指標の安定性: docs/06_TEST/TEST_19_sppと計算時間ノイズ指標安定性.md（spp 256〜8192 × 6 条件．GPU 時間は spp に比例・固定費 0.8 秒，ノイズは 1/√spp 則に一致．構造 CoV は 1024 で十分，画像全体 min/max 比を暗い条件で比べるなら 2048〜4096） +- 実機とシミュの形の突き合わせと配光の検討: docs/06_TEST/TEST_20_実機とシミュの形の突き合わせと配光の検討.md（実機 DNG とシミュを同じノイズ処理で全灯 + 単灯 3 灯の形で比較．お椀は ρ で合うが単灯のピークのずれと傾きはシミュが半分．壁 ω・α のスイープでは形全体を崩さずに合わせられず，LED 配光を spot FWHM 45〜60° に狭めると最良．共通の左上暗の傾きは未解決） - 照明配置依存性の球と立方体の比較: docs/06_TEST/TEST_18_照明配置依存性の球と立方体の比較.md（LED の位置 θ・向き・個数を同内面積ペアで振った 198 条件．配置依存は球が立方体の 1/2〜1/3，対向 2 灯以上で形状差が消える．追加実験: 立方体は面の稜へ寄せると悪化，真下向き多灯はバッフル無しで成立するが板向き多灯には劣る，校正材質 ρ=0.81 でも結論不変だが軸上は両形状とも最良になる） diff --git "a/docs/06_TEST/TEST_02_\343\202\267\343\203\237\343\203\245\343\203\254\343\203\274\343\202\267\343\203\247\343\203\263\345\256\237\350\241\214\350\250\230\351\214\262.md" "b/docs/06_TEST/TEST_02_\343\202\267\343\203\237\343\203\245\343\203\254\343\203\274\343\202\267\343\203\247\343\203\263\345\256\237\350\241\214\350\250\230\351\214\262.md" index 65f1ec3..2b3f1bc 100644 --- "a/docs/06_TEST/TEST_02_\343\202\267\343\203\237\343\203\245\343\203\254\343\203\274\343\202\267\343\203\247\343\203\263\345\256\237\350\241\214\350\250\230\351\214\262.md" +++ "b/docs/06_TEST/TEST_02_\343\202\267\343\203\237\343\203\245\343\203\254\343\203\274\343\202\267\343\203\247\343\203\263\345\256\237\350\241\214\350\250\230\351\214\262.md" @@ -94,6 +94,8 @@ | `output/led_placement_sweep/` (255 ディレクトリ) | TEST_18 | `configs/_tmp_led_placement/scene_*.yaml`（`scripts/sweep_led_placement.py` が生成．同内面積ペア 球 φ ∈ {200,270,340,500}（r ∈ {100,135,170,250}） ↔ 立方体 L ∈ {141,187,234,340} で，LED をキャビティ中心から見た極角 θ=0〜120°（13 段階）に置く位置スイープ・θ=60/90/120° の向き 3 通り（板向き＋バッフル／壁法線向き＋バッフル／真下向き）・θ=90° の輪に N=1〜4 灯（総光束保存）・seed 違いの再レンダ，追加実験の方位角 φ=15〜40°（立方体）と真下向き N=2/4/8 灯．材質・spp=1024・d=16 据え置き．Newton cuda_ad_rgb・mitsuba 3.7.1・`--workers 2` で 1 条件 5〜6 秒） | `configs/eval_baseline_dng_target_shape.yaml` | 配置依存の幅（θ 全域の構造 CoV max−min）は球 0.015〜0.023・立方体 0.025〜0.071．球は θ≈30° が軸上より均一（φ340: 0.0193→0.0057）．対向 2 灯で構造 CoV 0.005〜0.010 と形状差が消える．真下向きはバッフル無しで板向きと同等 — [TEST_18](TEST_18_照明配置依存性の球と立方体の比較.md) 参照 | **配置依存は球が立方体の 1/2〜1/3**（仮説成立）．図なし評価（`--no-figures`）なので各条件に figures/ は無い．`row.json` が条件ごとの集約行（並列・再開用） | | `output/spp_noise_sweep/` (72 ディレクトリ + `summary.csv` / `spp_summary_table.csv` / `spp_noise.png`) | TEST_19 | `configs/_tmp_spp_noise/scene_*.yaml`（`scripts/sweep_spp_noise.py` が TEST_18 の 6 条件の YAML を流用し spp と seed だけ書き換えて生成．球 φ340 軸上/θ=90°・立方体 L=234 軸上/θ=90°・球 φ500 θ=90°・立方体 L=340 θ=90° × spp 256〜8192 × seed 2．Newton GPU 0・cuda・逐次） | `configs/eval_baseline_dng_target_shape.yaml` | GPU 時間は spp に比例（1024 で 5.1 秒），固定費 0.8 秒．ノイズは 1/√spp 則に ±2% で一致．構造 CoV は 1024 で安定，画像全体 min/max は暗い条件で 2048〜4096 が必要 — [TEST_19](TEST_19_sppと計算時間ノイズ指標安定性.md) 参照 | 段階別の壁時計（`mi.render` は `dr.sync_thread` で同期）と seed 差を row.json に記録．`build_spp_noise_comparison.py` で図・表 | | `output/led_placement_sweep_rho081/` (60 ディレクトリ + 解析出力) | TEST_18 | `configs/_tmp_led_placement_rho081/scene_*.yaml`（`sweep_led_placement.py --material calibrated --key-subset` が生成．内壁を TEST_14 の校正材質 純 diffuse ρ=0.81 に置き換えた主要 60 条件＝位置スイープ ペア 1・3，多灯 N=2/4・真下向き ペア 3） | 同上 | 配置依存の幅 球 0.030 vs 立方体 0.077〜0.098．軸上は両形状とも大きく改善（旧材質の光沢が作るリングが消える）— [TEST_18](TEST_18_照明配置依存性の球と立方体の比較.md) 結果 7 | 主結論が材質に依らないことの確認．解析は `build_led_placement_comparison.py --material calibrated` | +| `output/real_vs_sim/` (`sim_rho081_spp4096/`・`sim_single_{bottom,left,right}_rho081_spp4096/` + `real_vs_sim_*.{png,csv}`・`shapes_{legacy,rho081}.{png,csv}`・`bsdf_grid_vs_real.csv`) | TEST_20 | `configs/_tmp_real_vs_sim/scene_*.yaml`（`scene_baseline.yaml` の内壁を TEST_14 の校正材質 純 diffuse ρ=0.81 に置き換え，全灯と単灯 3 灯（`led_bottom/left/right_*` のみ）を spp 4096 でレンダ．Newton GPU 0・cuda） | — （`compare_real_vs_sim_metrics.py` / `compare_real_vs_sim_shapes.py` が luminance.npy と実機 DNG を同じ処理で直接比較） | 全灯の画像全体 min は実機 0.727 vs シミュ 0.787（輪帯平均は 0.873 vs 0.874 で一致）．形は 4 構成とも「左上が暗い」傾きがシミュで半分，単灯のピークのずれも半分以下 — [TEST_20](TEST_20_実機とシミュの形の突き合わせと配光の検討.md) 結果 1・2 | 実機 DNG（`reference_measurement/quantitative/` 全灯・`lightsim/raw/` 単灯）とシミュを同じノイズ処理で突き合わせる基準データ．採用材質の比較は既存 `output/full_base_fix`・`output/singles_fix/` を流用 | +| `output/wall_bsdf_spot/` (228 ディレクトリ + `shape_fit_a.csv` / `shape_fit_b.csv` / `wall_bsdf_heatmap.png` / `spot_fwhm.png`) | TEST_20 | `configs/_tmp_wall_bsdf_spot/scene_*.yaml`（`scripts/sweep_wall_bsdf_spot.py` がサーバ上で生成．フェーズ A: 内壁 blendbsdf の ρ ∈ {0.85, 0.81} × ω ∈ {0〜0.70} × α ∈ {0.02〜0.6} × 4 構成 = 208 条件・spp 2048．フェーズ B: 採用材質のまま LED を全光束一致の spot（FWHM 30〜120°）× 4 構成 = 20 条件・spp 4096．Newton GPU 0・cuda・`--workers 2`） | — （`build_shape_fit_comparison.py` が実機と直接比較） | ω を上げると単灯ピークは実機に届く（ω=0.15）がお椀・傾きが消え ω≥0.3 で形が崩れる．spot FWHM 45〜60° で 4 構成の差分 RMS 平均 0.034（面光源 0.036）と最良 — [TEST_20](TEST_20_実機とシミュの形の突き合わせと配光の検討.md) 結果 3・4 | 単灯ピークのずれは壁の光沢ではなく LED 配光で説明できることの根拠．図なし・`row.json` が条件ごとの記録 | | `output/led_placement_sweep/{summary.csv,analysis.csv,position_summary.csv,position_sensitivity.png,position_maps.png,direction_comparison.png,count_comparison.png,azimuth_comparison.png,ring_comparison.png,analysis_report.txt}`・`_placement_check/cutaway_*.png` | TEST_18 | — | — | 条件別集約 CSV／方位フーリエ m=1〜4・構造 CoV・重心の解析 CSV／位置・向き・灯数・方位角・真下向き多灯の解析図 6 枚／レンダ前 3D 確認の半割り断面 23 枚 | `scripts/build_led_placement_comparison.py`・`show_scene_cutaway.py` の出力．レンダ前 3D 確認（CLAUDE.md 必須ルール）の記録を兼ねる | | `output/cube_size_sweep/board_baffle/{summary.csv,cube_vs_sphere.png,fourfold_component.png}`・`output/cube_size_sweep/_placement_check/cutaway_L*.png` | TEST_16 | — | — | 条件別集約 CSV／球との比較図（輝度マップ 5 枚 + 輪帯プロファイル + 内面積軸の指標カーブ + 直接ペア 3 組）／m=4（cos 4φ）成分の定量図／レンダ前 3D 確認の半割り断面 5 枚 | `scripts/build_cube_sphere_comparison.py`・`analyze_fourfold_component.py`・`show_scene_cutaway.py` の出力．m=4 は全条件 \|a4\| ≤ 0.26%（系統誤差床 ~0.1%）＝四隅の写り込みは実質ゼロ．レンダ前 3D 確認（CLAUDE.md 必須ルール）の記録を兼ねる | | `output/_slide_panels/sphere_loft_disk_before_after.png`・`normal_artifact_explainer.png` | TEST_12 | — | — | 球ロフト円形光源の修正前後 1×2 比較（相対輝度・共通スケール）／法線アーティファクトの機構説明図（模式図＋実物メッシュの補間法線再現＋修正前マップの 4 パネル） | 発表スライド素材（アドホック生成．修正前後比較は CoV 0.1205→0.0628 の視覚版） | @@ -155,6 +157,8 @@ | 球ロフト円形光源版レンダ・評価 (TEST_11 追補) | `$env:DRJIT_LIBLLVM_PATH = "C:\Program Files\LLVM\bin\LLVM-C.dll"; uv run sm-render --scene configs/scene_sphere_r100_loft_disk.yaml --output output/sphere_r100_loft_disk/` → `uv run sm-evaluate --image output/sphere_r100_loft_disk/luminance.npy --config configs/eval_baseline_dng_target_shape.yaml --output output/sphere_r100_loft_disk/` | | 単灯再検証レンダ・集計 (TEST_12 追補) | `$env:DRJIT_LIBLLVM_PATH = "C:\Program Files\LLVM\bin\LLVM-C.dll"; uv run sm-render --scene configs/_tmp_fix_{bottom,left,right}.yaml --output output/singles_fix/{bottom,left,right}/` → `uv run python scripts/build_singles_fix_summary.py` | | spp スイープ (TEST_19) | `uv run python scripts/sweep_spp_noise.py --prepare-only` → `uv run sm-remote push` → `uv run sm-remote run --name spp --gpu 0 --threads 16 -- python scripts/sweep_spp_noise.py` → `uv run sm-remote pull spp_noise_sweep` → `uv run python scripts/build_spp_noise_comparison.py` | +| 実機 vs シミュ 形の突き合わせ (TEST_20) | `uv run python scripts/compare_real_vs_sim_metrics.py --sim output/real_vs_sim/sim_rho081_spp4096/luminance.npy --tag rho081` → `uv run python scripts/compare_real_vs_sim_shapes.py --sim-set {legacy,rho081}` | +| 壁反射特性・spot スイープ (TEST_20) | `uv run sm-remote push --only scripts` → `uv run sm-remote run --name wall --gpu 0 --threads 16 -- python scripts/sweep_wall_bsdf_spot.py --phase all --workers 2` → `uv run sm-remote pull wall_bsdf_spot` → `uv run python scripts/build_shape_fit_comparison.py` | | 照明配置スイープ (TEST_18) | `uv run sm-remote push` → `uv run sm-remote run --name place --gpu 0 --threads 16 -- python scripts/sweep_led_placement.py --workers 2` → `uv run sm-remote pull led_placement_sweep` → `uv run python scripts/build_led_placement_comparison.py`（手元で幾何だけ確認するなら `--prepare-only`） | ## 過去の整理実績 (Cleanup History) diff --git "a/docs/06_TEST/TEST_20_\345\256\237\346\251\237\343\201\250\343\202\267\343\203\237\343\203\245\343\201\256\345\275\242\343\201\256\347\252\201\343\201\215\345\220\210\343\202\217\343\201\233\343\201\250\351\205\215\345\205\211\343\201\256\346\244\234\350\250\216.md" "b/docs/06_TEST/TEST_20_\345\256\237\346\251\237\343\201\250\343\202\267\343\203\237\343\203\245\343\201\256\345\275\242\343\201\256\347\252\201\343\201\215\345\220\210\343\202\217\343\201\233\343\201\250\351\205\215\345\205\211\343\201\256\346\244\234\350\250\216.md" new file mode 100644 index 0000000..0e8a2db --- /dev/null +++ "b/docs/06_TEST/TEST_20_\345\256\237\346\251\237\343\201\250\343\202\267\343\203\237\343\203\245\343\201\256\345\275\242\343\201\256\347\252\201\343\201\215\345\220\210\343\202\217\343\201\233\343\201\250\351\205\215\345\205\211\343\201\256\346\244\234\350\250\216.md" @@ -0,0 +1,212 @@ +# 実機とシミュの形の突き合わせと配光の検討 (Real vs Sim Shape Matching and LED Beam Pattern) + +舌撮影デバイスの基準「最大を 1 としたときの最小 ≥ 0.9」に対し，実機を作る前にシミュレーションが +参考になるかを確かめるため，実機 DNG（全灯 1 枚・単灯 3 枚）とシミュを**同じノイズ処理**で +突き合わせた（2026-09-14〜15）．評価の決め手は集約指標（min 等）ではなく**輝度マップの形** +（全灯 + 単灯 3 灯）とした（ユーザー方針）．結論: **同心円方向（お椀）は材質で合うが，単灯で +明部が点灯基板の反対側へ寄る量がシミュは実機の半分以下**．これは壁の光沢（ω・α）では +形全体を崩さずに合わせられず，**LED の配光を狭める（spot，半値全幅 45〜60°）** ことで説明できる． +全構成に共通する「左上が暗い」傾き成分はどちらでも半分しか出ず，基板の位置・向きか +ISP の LSC 残差が残る課題である． + +関連: + +- 材質パラメータの確定（ρ≈0.81 純 diffuse）と傾き成分の分離: [TEST_14](TEST_14_材質パラメータ再スイープと残差の構造分解.md) +- 単灯の相対結論: [TEST_12 追補](TEST_12_シェーディング法線アーティファクトの発見と修正.md) +- 画像全体 min/max 比の定義: [SPEC_04 評価仕様](../04_SPEC/SPEC_04_評価仕様.md)，[TEST_18](TEST_18_照明配置依存性の球と立方体の比較.md)「均一性の指標」 +- spp の指針: [TEST_19](TEST_19_sppと計算時間ノイズ指標安定性.md) +- 実行記録: [TEST_02_シミュレーション実行記録](TEST_02_シミュレーション実行記録.md) + +## 実験条件 (Setup) + +### 実機データ + +| 構成 | ファイル | 備考 | +| --- | --- | --- | +| 全灯 18 灯 | `reference_measurement/quantitative/MiniTIAS_QM_20260601_124500.dng` | 定量モード・ISO 40・8.3 ms・LSC マップ付き | +| 下 / 左 / 右 単灯 | `lightsim/raw/SmTIAS_QM_20260615_{141420,150300,170948}.dng` | 同一セッション・露出固定 | + +読み込みは TEST_14 と同じ経路（`load_image_dng` → lscMap 順適用 → Rec.709 → rot90(−1) → 白板 ROI +1400×1409 px）．ROI はシミュ ROI（464×450 px）へリサイズしてから同じ処理を掛ける（平滑化 σ を +画素で揃えるため．約 3×3 画素の平均でセンサノイズも落ちるが，ノイズ補正はリサイズ後の残差から +推定するので整合する）． + +### シミュレーション + +基準シーン `scene_baseline.yaml`（実機と同じ筐体・18 灯）から，単灯は `led_bottom_*` / +`led_left_*` / `led_right_*` の 6 灯だけを残して作る．材質は次の 2 つを主に使った． + +- **採用材質**（TEST_15〜18 で使用）: blendbsdf ρ=0.85・ω=0.05・α=0.1（既存レンダ `full_base_fix` / + `singles_fix`，spp 1024） +- **TEST_14 の確定値**: 純 diffuse ρ=0.81（本実験で spp 4096 でレンダ．`output/real_vs_sim/sim_*rho081*`） + +レンダは研究室サーバ Newton GPU 0（cuda_ad_rgb・mitsuba 3.7.1）． + +### 処理と指標（実機・シミュ共通） + +- 形の比較: メディアン 5 px → ガウス σ=4 px → ROI 平均で正規化（=1）．差分（シミュ − 実機）の + RMS・最大，SSIM，NCC． +- 形の記述子: σ=20 平滑化後のピーク位置（ROI 中心からのずれ，半幅比，+x 右・+y 上）， + 上下比・左右比（上 1/3 と下 1/3 等の平均比），1 次平面フィットの傾き（振幅 %・方位）， + 2 次成分（お椀）の振幅，方位フーリエ a1/b1/a2/a4（TEST_16 の機構）． +- min: 画像全体・ノイズ補正付き最小/最大比（σ=12．`evaluation/fullimage.py`）．参考に補正なし・ + 輪帯平均・画素． +- 20 分割放射プロファイル（ROI 平均で正規化）の各ビンの差． + +## 結果 1 — 全灯の min・放射プロファイル・形 (Full-Board: min, Profile, Shape) + +`output/real_vs_sim/real_vs_sim_{legacy,rho081}.png`． + +| 指標 | 実機 | ρ=0.81 純 diffuse（spp 4096） | 採用材質 ρ=0.85 ω=0.05 α=0.1 | +| --- | --- | --- | --- | +| min（画像全体・補正） | **0.727** | 0.787（+0.060） | 0.803（+0.075） | +| min（輪帯平均） | 0.873 | 0.874（+0.001） | 0.901（+0.028） | +| 構造 CoV | 0.067 | 0.046 | 0.039 | +| 20 分割プロファイル RMSE / 最大差 | — | 0.0055 / 0.011 | 0.012 / — | +| 正規化マップ RMS 差 / 最大差 | — | 0.038 / 0.103 | 0.046 / — | +| SSIM / NCC | — | 0.996 / 0.858 | — / 0.858 | +| 方位 m=1 a1 / b1 | +3.8% / +7.1% | −0.3% / +4.0% | −0.2% / +3.6% | + +- **同心円方向はよく合う**: 放射プロファイルの差は最大 0.011，輪帯平均の min は 0.001 差． + ρ=0.81 が TEST_14 どおり妥当． +- **画像全体の min は 0.06 高く出る**．原因は方位方向の非対称で，実機は上（LED と反対側）が暗く + 左右も 6% 非対称なのに，シミュは上下差が半分・左右差ゼロ．TEST_14 で「材質では直せない傾き + 成分」と分離した残差そのもので，輪帯平均では見えず画像全体の min には出る． +- 基準 0.9 に対する判定（不合格）は一致するが，差 0.06 は境目付近では効く． +- シミュのノイズを spp 4096 で半分に下げても min は 0.792 → 0.787 でほぼ不変（補正が機能）． + 実機のノイズ CoV は 0.8% で補正の影響はゼロ． +- TEST_14 のグリッド 59 条件（ρ 0.75〜0.92 × ω × α）を同じ指標で当てても min は 0.774〜0.800 で， + ω・α では差 0.06 は埋まらない（`output/real_vs_sim/bsdf_grid_vs_real.csv`）．ρ を下げると + お椀が深くなりすぎて放射プロファイルが崩れるだけ． + +## 結果 2 — 形の一致: 全灯 + 単灯 3 灯 (Shape Matching, 4 Configurations) + +`output/real_vs_sim/shapes_{legacy,rho081}.png`（上 = 実機・中 = シミュ・下 = 差分）． + +採用材質（ρ=0.85・ω=0.05・α=0.1）: + +| 構成 | 差分 RMS | SSIM | NCC | 傾き 実機 / シミュ | お椀 実機 / シミュ | ピーク位置 実機 / シミュ | +| --- | --- | --- | --- | --- | --- | --- | +| 全灯 | 0.040 | 0.92 | 0.86 | 9.6% / 4.1% | −8.2% / −6.5% | (+0.41, −0.38) / (+0.03, −0.13) | +| 下 単灯 | 0.039 | 0.90 | 0.84 | 8.1% / 3.4% | −10.2% / −6.7% | (+0.14, −0.43) / (−0.07, −0.22) | +| 左 単灯 | 0.045 | 0.90 | 0.92 | 11.9% / 5.3% | −9.4% / −6.2% | (+0.65, −0.52) / (+0.25, −0.31) | +| 右 単灯 | **0.022** | 0.91 | **0.94** | 6.6% / 4.8%（方位一致） | −7.9% / −6.5% | (−1.00, −0.32) / (−0.28, −0.20) | + +- **単灯の明部は点灯基板の反対側へ寄る**（左単灯なら右へ，右単灯なら左へ）．向きはシミュも + 実機も同じだが，**ずれの量はシミュが実機の半分以下**（左 +0.25 vs +0.65，右 −0.28 vs −1.00）． +- 差分は 4 構成すべてで「左上が赤・右下が青」の対角パターン＝実機は左上が暗く右下が明るい傾き + を持ち，シミュはそれを半分しか再現しない（振幅 実機 8〜12% vs シミュ 3〜5%）．右単灯だけは + 方位も振幅もほぼ合う．TEST_14 の「傾きは下・左の基板に固有」と整合するが，右にも弱い同じ向きの + 傾きが残り，共通成分（カメラの LSC 残差）もゼロではない． +- ρ=0.81 純 diffuse にするとお椀の深さは実機と一致し（全灯 −8.3% vs −8.2%），SSIM は 0.92 → 0.97 + に上がるが，**傾きとピークのずれは変わらない**（むしろ拡散が増えてピークのずれが小さくなる）． + +## 結果 3 — 壁の反射特性（ω・α）のスイープ (Wall BSDF Sweep) + +`scripts/sweep_wall_bsdf_spot.py` フェーズ A．ρ ∈ {0.85, 0.81} × ω ∈ {0, 0.05, 0.15, 0.30, 0.50, 0.70} +× α ∈ {0.02, 0.05, 0.1, 0.3, 0.6}（ω=0 は純 diffuse）× 4 構成 = 208 レンダ（spp 2048）． +`output/wall_bsdf_spot/wall_bsdf_heatmap.png`・`shape_fit_a.csv`． + +| ρ=0.85・α=0.1 | ω=0 | 0.05 | 0.15 | 0.30 | 0.50 | +| --- | --- | --- | --- | --- | --- | +| 左単灯 ピーク x（実機 +0.65） | +0.22 | +0.30 | **+0.61** | +0.72 | +0.80 | +| 右単灯 ピーク x（実機 −1.00） | −0.19 | −0.24 | −0.62 | −0.75 | −0.85 | +| 全灯 お椀（実機 −8.2%） | −7.1% | −6.4% | −5.2% | −3.6% | −2.0% | +| 全灯 傾き（実機 9.6%） | 4.2% | 4.0% | 3.4% | 1.9% | 2.3% | +| 4 構成の差分 RMS 平均 | 0.039 | 0.036 | 0.038 | 0.049 | 0.073 | + +- ω を上げると単灯のピークは実機の方向へ大きくずれ，ω=0.15 で実機に届く．**しかし同時にお椀と + 傾きが消えて形全体は崩れる**（差分 RMS は ω=0.30 で 0.049，0.50 で 0.073）． +- α は ω ≥ 0.3 でのみ効く（小さい α ＝鏡面に近いほどピークが端まで飛ぶ）．採用域 ω ≤ 0.15 では + ほぼ無関係． +- 4 構成の差分 RMS 平均の最良は ρ=0.81・ω=0.05・α=0.1（0.0336）．採用材質は 0.0364 で 7 位． + ピークの折り合いも含めた最良は ρ=0.81〜0.85・ω=0.15・α=0.1（ピーク誤差 0.30，RMS 0.036〜0.038）． + 差はいずれも小さく，**ω・α では「ピークのずれ」と「お椀・傾き」を同時には合わせられない**． + +## 結果 4 — LED 配光（spot）のスイープ (LED Beam Pattern Sweep) + +同スクリプト フェーズ B．採用材質のまま LED を Mitsuba の `spot`（面光源と同じ位置・向き・全光束） +にし，半値全幅 FWHM ∈ {30, 45, 60, 90, 120}°（cutoff = 2/3 FWHM，beam = 1/3 FWHM）× 4 構成 += 20 レンダ（spp 4096）．`output/wall_bsdf_spot/spot_fwhm.png`・`shape_fit_b.csv`． + +| FWHM | 30° | 45° | 60° | 90° | 120° | 面光源（ランバート） | +| --- | --- | --- | --- | --- | --- | --- | +| 左単灯 ピーク x（実機 +0.65） | **+0.68** | +0.60 | +0.51 | +0.38 | +0.29 | +0.30 | +| 右単灯 ピーク x（実機 −1.00） | −0.72 | −0.64 | −0.52 | −0.41 | −0.34 | −0.24 | +| 右単灯 傾き（実機 6.6%） | 8.3% | 7.4% | 6.2% | 5.0% | 4.7% | 4.8% | +| 左単灯 傾き（実機 11.9%） | 8.2% | 7.4% | 6.6% | 5.6% | 5.3% | 5.3% | +| 4 構成の差分 RMS 平均 | 0.035 | **0.034** | **0.034** | 0.036 | 0.037 | 0.036 | + +- FWHM を狭めると単灯のピークが実機に近づき，30° で左は実機に一致．傾きも実機側へ寄る． +- ただし FWHM を狭めすぎると全灯・左単灯のピーク y（上下）が実機から離れる． + **FWHM 45〜60° が全体のバランス点**（差分 RMS 平均 0.034 で面光源の 0.036 より良い）． +- 右単灯の差分 RMS は 90° で 0.022 と最良で，面光源と同等． + +## 結論 (Conclusions) + +- 実機との一致は「同心円方向（お椀）」「方位方向（ピークのずれ・傾き）」に分けて見る必要がある． + 材質（ρ）が決めるのは前者だけで，TEST_14 の ρ≈0.81 は妥当． +- 単灯の明部が反対側へ寄る量は，壁の光沢では形全体を崩さずに再現できず，**LED の配光を + ランバートより狭くする（spot FWHM 45〜60°）** ことで説明できる．新しい基準候補は + 「採用材質（ρ=0.85・ω=0.05・α=0.1）+ spot FWHM 45〜60°」． +- 全構成に共通する「左上が暗い」傾き（実機 8〜12%）はどの設定でも半分しか出ない．下・左の + 基板で大きく右で小さいので，基板の位置・向き（または基板ごとの出力差）が第一候補． +- シミュの min（画像全体）は実機より 0.05〜0.07 高く出る．基準 0.9 の合否をシミュ単独で断定する + には，この傾き成分を詰めるか，余裕を持った判定基準（例: シミュで 0.95 以上）にする必要がある． + +## 限界・注意点 (Limitations) + +- 実機は各構成 1 枚（全灯は 2026-06-01，単灯は 2026-06-15 の別セッション）．実機側のノイズ + フロアは撮り直しが無く測れていない（ノイズ CoV は 0.8% で小さい）． +- ISO の min の定義（画像全体か輪帯平均か）は未確認．本実験は画像全体を主，輪帯平均を参考にした． +- spot は点光源で，面光源の 2.8×3.5 mm の広がりは無い．FWHM の cutoff/beam の分け方（2:1）は + 仮置きで，データシートの配光曲線には合わせていない． +- 傾きの共通成分が LSC 残差か幾何かは未切り分け． +- 壁スイープは spp 2048（ノイズ CoV 約 7%）．ピーク位置は σ=20 の平滑化後なので影響は小さいが， + 差分 RMS の 0.002 程度の差は有意ではない． + +## 失敗パターン・方法論の学び (Pitfalls and Methodology) + +- **集約指標（min）だけで材質を選ぶと誤る**．ρ を下げれば min は実機に寄るが，それは + お椀を深くして帳尻を合わせているだけで放射プロファイルは崩れる．形（全灯 + 単灯 3 灯）で + 見ることで「材質で合わせられる成分」と「光源・幾何で合わせる成分」を分離できた． +- **単灯 3 灯は方位方向の診断に必須**．全灯だけでは 3 基板の寄与が混ざって傾きの原因を + 切り分けられない． +- サーバの push は `--only scripts` だけだと新しい src モジュールが無くて落ちる（TEST_19 と同じ）． + `configs/_tmp_*` は push 対象外なので，スイープ用 YAML はサーバ上で生成する設計にした． + +## 次にやること (Next Steps) + +1. 「採用材質 + spot FWHM 45〜60°」で下・左基板の向き（数度）と位置を振り，共通の傾きが + 埋まるかを見る． +2. LED データシートの配光曲線で spot の cutoff/beam を決め直す． +3. 傾きの共通成分の切り分け: 実機で白板を 180° 回転して撮る（傾きが回れば筐体，回らなければ + カメラ側）． +4. 単灯の実機データを全灯と同一セッションで撮り直し，絶対振幅も横並びにする． + +## 再現コマンド (Reproduction) + +```bash +# 全灯: min・放射プロファイル・形（ρ=0.81 は事前にサーバでレンダ） +uv run python scripts/compare_real_vs_sim_metrics.py --sim output/real_vs_sim/sim_rho081_spp4096/luminance.npy \ + --label "ρ=0.81 純 diffuse・spp4096" --tag rho081 +# 形の一致（全灯 + 単灯 3 灯） +uv run python scripts/compare_real_vs_sim_shapes.py --sim-set legacy +uv run python scripts/compare_real_vs_sim_shapes.py --sim-set rho081 +# 壁反射特性・spot スイープ（サーバ．YAML はサーバ上で生成） +uv run sm-remote push --only scripts +uv run sm-remote run --name wall --gpu 0 --threads 16 -- python scripts/sweep_wall_bsdf_spot.py --phase all --workers 2 +uv run sm-remote pull wall_bsdf_spot +uv run python scripts/build_shape_fit_comparison.py +``` + +出力: `output/real_vs_sim/`（`real_vs_sim_*.png/csv`・`shapes_*.png/csv`・`bsdf_grid_vs_real.csv`・ +ρ=0.81 のレンダ 4 本），`output/wall_bsdf_spot/`（228 レンダ + `shape_fit_a.csv` / `shape_fit_b.csv` + +`wall_bsdf_heatmap.png` / `spot_fwhm.png`）． + +## 参照 (References) + +- [TEST_14_材質パラメータ再スイープと残差の構造分解](TEST_14_材質パラメータ再スイープと残差の構造分解.md): ρ の確定と傾き成分 +- [TEST_18_照明配置依存性の球と立方体の比較](TEST_18_照明配置依存性の球と立方体の比較.md): 画像全体 min/max 比の導入 +- [SPEC_04_評価仕様](../04_SPEC/SPEC_04_評価仕様.md): 指標の定義 +- [TEST_02_シミュレーション実行記録](TEST_02_シミュレーション実行記録.md): 出力ディレクトリの索引 diff --git a/docs/PROGRESS.md b/docs/PROGRESS.md index 3acd6e7..aaaf797 100644 --- a/docs/PROGRESS.md +++ b/docs/PROGRESS.md @@ -106,3 +106,4 @@ - **正規化と最小/最大比の定義の見直し — 画像全体・ノイズ補正付きの 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 の文書・図．条件名は半径のまま）． - **TEST_19 spp と計算時間・ノイズ・指標の安定性 (2026-09-14)**: ユーザー発案「GPU で速くなったので spp を上げてノイズと計算時間の関係を測りたい．全体のうち GPU が何秒かも知りたい」．TEST_18 の 6 条件を spp 256〜8192 × seed 2 で回し，段階別の壁時計（`mi.render` は `dr.sync_thread` で同期）を記録した．**結果**: GPU 時間は spp に厳密に比例（1024 で 5.1 秒），CPU 側の固定費は 0.8 秒で一定（PLY 書き出し・評価・保存）＝spp ≥ 1024 で 1 割以下なので CPU 側の最適化はここまで．ノイズは全条件で 1/√spp 則に ±2% で一致（レンダ側に別のノイズ源は無い．ノイズ量は生の輝度で決まり，オフ軸は軸上の 1.5 倍）．構造 CoV は 1024 で安定（8192 との差 ≤0.002）だが，画像全体 min/max 比は暗い大キャビティで 1024 だと 0.01〜0.02 過小評価で 2048〜4096 で収束．補正なしの平滑化 min/max は 8192 でも収束せず，輪帯平均 min/max は spp では改善しない（方位平均の問題）．**推奨**: 構造 CoV の実験は 1024 のまま，min/max 比を暗い条件まで比べるなら 2048〜4096．**失敗パターン**: ① `--only scripts` の push では新しい src モジュール（fullimage.py）が無くて落ちた → 全体を push / ② サンプラー再起動の `pkill -f` が自分自身を殺して nvidia-smi の記録が残らなかった（結論には影響なし）．16384 はユーザー都合で未実施． +- [x] **TEST_20 実機とシミュの形の突き合わせと配光の検討 — 単灯ピークのずれは壁の光沢ではなく LED 配光で説明できる (2026-09-15)**: 実機を作る前にシミュが ISO 基準（min ≥ 0.9）の参考になるかを確かめるため，実機 DNG（全灯 1 枚・単灯 3 枚）とシミュを**同じノイズ処理**（メディアン→ガウス→平均正規化，画像全体 min はノイズ補正付き）で突き合わせた．評価の決め手は min 等の集約指標ではなく**輝度マップの形（全灯 + 単灯 3 灯）**とした（ユーザー方針）．**結果 ①**: 同心円方向（お椀・20 分割プロファイル）は TEST_14 の ρ≈0.81 で最大差 0.011 とよく合い，輪帯平均の min は 0.001 差．しかし画像全体の min はシミュが 0.06 高い＝実機の方位方向の非対称（上が暗い・左右 6%）をシミュが半分しか出さないため．TEST_14 のグリッド 59 条件を当てても ω・α では埋まらない（ρ を下げるとお椀が深くなるだけ）．**結果 ②**: 4 構成すべてで差分は「左上が暗く右下が明るい」対角の傾き（実機 8〜12% vs シミュ 3〜5%），単灯の明部が点灯基板の反対側へ寄る量は向きは合うが実機の半分以下（左 +0.65 vs +0.25，右 −1.00 vs −0.28）．右単灯だけはほぼ一致．**結果 ③（壁 ω・α スイープ 208 条件）**: ω を上げると単灯ピークは ω=0.15 で実機に届くが，同時にお椀と傾きが消えて形全体は崩れる（ω≥0.3 で差分 RMS 1.4〜2 倍）．α は ω≥0.3 でしか効かない．ω・α では「ピークのずれ」と「お椀・傾き」を同時に合わせられない．**結果 ④（spot 配光 20 条件）**: 採用材質のまま LED を全光束一致の spot にして FWHM を狭めると単灯ピークが実機に近づき（30° で左が一致），傾きも実機側へ寄る．**FWHM 45〜60° が全体のバランス点**（差分 RMS 平均 0.034 vs 面光源 0.036）で，新しい基準候補は「採用材質 ρ=0.85・ω=0.05・α=0.1 + spot 45〜60°」．**未解決**: 全構成に共通する左上暗の傾きはどの設定でも半分しか出ない．下・左基板で大きく右で小さいので基板の位置・向き（または出力差）が第一候補．次は spot 45〜60° で基板の向き・位置を振る．**方法論の学び**: min だけで材質を選ぶと ρ を下げて帳尻を合わせる誤りに陥る．形で見ることで「材質で合う成分（お椀）」と「光源・幾何で合わせる成分（方位）」を分離できた．単灯 3 灯は方位方向の診断に必須．**失敗パターン**: ① `--only scripts` の push では新規 src が無い（TEST_19 と同じ）/ ② `configs/_tmp_*` は push 対象外なのでスイープ用 YAML はサーバ上で生成する設計にした / ③ 同名ジョブのログは追記なので `grep exit` の待ち合わせが旧 run に当たる．成果物: `output/real_vs_sim/`・`output/wall_bsdf_spot/`，`scripts/compare_real_vs_sim_{metrics,shapes}.py`・`sweep_wall_bsdf_spot.py`・`build_shape_fit_comparison.py`，記録 [TEST_20](06_TEST/TEST_20_実機とシミュの形の突き合わせと配光の検討.md)． diff --git a/scripts/build_shape_fit_comparison.py b/scripts/build_shape_fit_comparison.py new file mode 100644 index 0000000..96ab4e3 --- /dev/null +++ b/scripts/build_shape_fit_comparison.py @@ -0,0 +1,239 @@ +"""壁の反射特性（ω・α）と spot 配光のスイープを，実機との「形」の一致で評価する（TEST_20）. + +sweep_wall_bsdf_spot.py の出力（output/wall_bsdf_spot/<条件>/luminance.npy）を読み， +compare_real_vs_sim_shapes.py と同じ処理で実機（全灯・下・左・右）と突き合わせる． + +指標（構成ごと）: + rms_diff : ROI 平均正規化マップの差分 RMS（小さいほど形が近い） + peak_dx/dy : σ=20 平滑化後のピーク位置（半幅比．+x 右・+y 上）．実機との差 d_peak + tilt / az : 1 次平面フィットの振幅 % と方位 + bowl : 2 次成分（お椀）の振幅 % + +出力（output/wall_bsdf_spot/）: + shape_fit_a.csv / shape_fit_b.csv : 条件 × 構成の指標 + wall_bsdf_heatmap.png : ρ・構成ごとに ω×α の差分 RMS とピークずれ誤差のヒートマップ + spot_fwhm.png : FWHM に対するピークずれ・傾き・差分 RMS（実機の値を水平線で） + +実行方法: + uv run python scripts/build_shape_fit_comparison.py +""" + +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 + +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 compare_real_vs_sim_shapes import ( # noqa: E402 + CONFIG_LABEL, + CONFIGS, + REAL, + load_real, + load_sim, + normalized_shape, + plane_fit, +) +from sweep_wall_bsdf_spot import ALPHAS, OMEGAS, OUTPUT_BASE, RHOS, SPOT_FWHMS # 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"] +MAX_DISTANCE = 0.95 +PEAK_SIGMA = 20.0 + + +def describe(shape_map: np.ndarray) -> dict[str, float]: + """正規化マップの形の記述子（ピーク位置・傾き・お椀）.""" + h, w = shape_map.shape + heavy = gaussian_filter(shape_map, sigma=PEAK_SIGMA) + iy, ix = np.unravel_index(int(np.argmax(heavy)), heavy.shape) + tilt, az, bowl = plane_fit(shape_map) + return {"peak_dx": (ix - w / 2) / (w / 2), "peak_dy": -(iy - h / 2) / (h / 2), + "tilt": tilt, "tilt_az": az, "bowl": bowl} + + +def compare(sim_roi: np.ndarray, real_map: np.ndarray, real_desc: dict) -> dict[str, float]: + """1 条件の形をシミュ vs 実機で比べる.""" + sim_map = normalized_shape(sim_roi) + inner = compute_normalized_distance_map(*sim_map.shape) <= MAX_DISTANCE + diff = sim_map - real_map + d = describe(sim_map) + d["rms_diff"] = float(np.sqrt(np.mean(diff[inner] ** 2))) + d["d_peak"] = float(np.hypot(d["peak_dx"] - real_desc["peak_dx"], + d["peak_dy"] - real_desc["peak_dy"])) + d["d_tilt"] = float(d["tilt"] - real_desc["tilt"]) + return d + + +def load_all() -> tuple[list[dict], dict[str, dict]]: + """全条件の row.json を読み，実機と比べた指標を付けて返す.""" + real_maps: dict[str, np.ndarray] = {} + real_desc: dict[str, dict] = {} + rows = [] + for path in sorted(OUTPUT_BASE.glob("*/row.json")): + row = json.loads(path.read_text(encoding="utf-8")) + npy = path.parent / "luminance.npy" + if not npy.exists(): + continue + sim_roi = load_sim(npy) + cfg = row["config"] + if cfg not in real_maps: + real_roi = load_real(REAL[cfg], sim_roi.shape) + real_maps[cfg] = normalized_shape(real_roi) + real_desc[cfg] = describe(real_maps[cfg]) + row.update(compare(sim_roi, real_maps[cfg], real_desc[cfg])) + rows.append(row) + return rows, real_desc + + +def write_csv(rows: list[dict], phase: str) -> None: + sel = [r for r in rows if r["phase"] == phase] + if not sel: + return + path = OUTPUT_BASE / f"shape_fit_{phase}.csv" + with path.open("w", newline="", encoding="utf-8") as fh: + writer = csv.DictWriter(fh, fieldnames=list(sel[0].keys())) + writer.writeheader() + writer.writerows(sel) + + +def draw_heatmaps(rows: list[dict], real_desc: dict, output: Path) -> None: + """ρ × 構成 ごとに ω×α の差分 RMS とピークずれ誤差を描く（ω=0 は α 列の先頭に置く）.""" + sel = [r for r in rows if r["phase"] == "a"] + if not sel: + return + omegas = OMEGAS + figsize = (4.4 * len(CONFIGS), 3.6 * 2 * len(RHOS)) + fig, axes = plt.subplots(2 * len(RHOS), len(CONFIGS), figsize=figsize) + for ri, rho in enumerate(RHOS): + for ci, cfg in enumerate(CONFIGS): + for mi, (key, title, cmap) in enumerate((("rms_diff", "差分 RMS", "viridis_r"), + ("d_peak", "ピーク位置の誤差", "magma_r"))): + grid = np.full((len(omegas), len(ALPHAS)), np.nan) + for r in sel: + if r["rho"] != rho or r["config"] != cfg: + continue + oi = omegas.index(r["omega"]) + if r["omega"] == 0.0: + grid[oi, :] = r[key] + else: + grid[oi, ALPHAS.index(r["alpha"])] = r[key] + ax = axes[2 * ri + mi, ci] + im = ax.imshow(grid, cmap=cmap, aspect="auto") + ax.set_xticks(range(len(ALPHAS))) + ax.set_xticklabels([f"{a:g}" for a in ALPHAS]) + ax.set_yticks(range(len(omegas))) + ax.set_yticklabels([f"{o:g}" for o in omegas]) + ax.set_xlabel("α（粗さ）") + ax.set_ylabel("ω（光沢の重み）") + for oi in range(len(omegas)): + for ai in range(len(ALPHAS)): + if not np.isnan(grid[oi, ai]): + ax.text(ai, oi, f"{grid[oi, ai]:.3f}", ha="center", va="center", + fontsize=7, color="white") + ax.set_title(f"ρ={rho:g} {CONFIG_LABEL[cfg]}: {title}", fontsize=9) + fig.colorbar(im, ax=ax, fraction=0.046) + fig.suptitle("フェーズ A: 壁の反射特性（ω×α）と実機との形の一致" + "（小さいほど良い．ω=0 は純 diffuse）", fontsize=12) + fig.tight_layout() + fig.savefig(output, dpi=100) + plt.close(fig) + print(f"出力: {output}") + + +def draw_spot(rows: list[dict], real_desc: dict, output: Path) -> None: + """フェーズ B: FWHM に対するピークずれ・傾き・差分 RMS.""" + sel = [r for r in rows if r["phase"] == "b"] + if not sel: + return + fig, axes = plt.subplots(1, 4, figsize=(20, 5)) + colors = dict(zip(CONFIGS, plt.cm.tab10(np.arange(4)), strict=True)) + # 参照: 採用材質の面光源（フェーズ A の r0.85 w0.05 a0.10） + ref = {r["config"]: r for r in rows + if r["phase"] == "a" and r["rho"] == 0.85 and r["omega"] == 0.05 and r["alpha"] == 0.1} + for ax, key, ylabel in ((axes[0], "peak_dx", "ピーク x（半幅比・右+）"), + (axes[1], "peak_dy", "ピーク y（半幅比・上+）"), + (axes[2], "tilt", "傾き [%]"), (axes[3], "rms_diff", "差分 RMS")): + for cfg in CONFIGS: + pts = sorted([r for r in sel if r["config"] == cfg], key=lambda r: r["fwhm"]) + ax.plot([r["fwhm"] for r in pts], [r[key] for r in pts], marker="o", + color=colors[cfg], label=f"{CONFIG_LABEL[cfg]}（spot）") + if key != "rms_diff": + ax.axhline(real_desc[cfg][key], color=colors[cfg], ls="--", alpha=0.7) + if cfg in ref: + ax.scatter([180.0], [ref[cfg][key]], marker="s", color=colors[cfg], zorder=3) + ax.set_xlabel("spot の半値全幅 FWHM [deg]（180 = ランバート面光源）") + ax.set_ylabel(ylabel) + ax.set_xticks(SPOT_FWHMS + [180.0]) + ax.grid(alpha=0.3) + ax.legend(fontsize=7) + axes[0].set_title("点線 = 実機の値，□ = 採用材質の面光源", fontsize=10) + fig.suptitle("フェーズ B: LED の配光（spot FWHM）と単灯の明部のずれ・傾き", fontsize=12) + fig.tight_layout() + fig.savefig(output, dpi=100) + plt.close(fig) + print(f"出力: {output}") + + +def print_summary(rows: list[dict], real_desc: dict) -> None: + """要点を表で出す.""" + print("\n=== 実機の形 ===") + for cfg in CONFIGS: + d = real_desc.get(cfg) + if d: + print(f" {CONFIG_LABEL[cfg]:>8}: peak=({d['peak_dx']:+.2f},{d['peak_dy']:+.2f}) " + f"傾き {d['tilt']:.1f}%@{d['tilt_az']:.0f}° お椀 {d['bowl']:+.1f}%") + for phase, title in (("a", "フェーズ A（壁）: 4 構成の差分 RMS 平均が小さい順"), + ("b", "フェーズ B（spot）: FWHM ごと")): + sel = [r for r in rows if r["phase"] == phase] + if not sel: + continue + print(f"\n=== {title} ===") + keys = sorted({(r["rho"], r["omega"], r["alpha"], r.get("fwhm")) for r in sel}) + summary = [] + for k in keys: + grp = {r["config"]: r for r in sel + if (r["rho"], r["omega"], r["alpha"], r.get("fwhm")) == k} + if len(grp) < len(CONFIGS): + continue + summary.append((k, float(np.mean([grp[c]["rms_diff"] for c in CONFIGS])), + float(np.mean([grp[c]["d_peak"] for c in CONFIGS])), grp)) + summary.sort(key=lambda s: s[1] if phase == "a" else (s[0][3] or 0)) + print(f"{'ρ':>5} {'ω':>5} {'α':>5} {'FWHM':>5} {'RMS平均':>8} {'ピーク誤差平均':>12} | " + f"{'左 peak x':>9} {'右 peak x':>9} {'全灯 傾き%':>10}") + for k, rms, dpk, grp in summary[: (12 if phase == "a" else 99)]: + fwhm = str(k[3] or "-") + print(f"{k[0]:5.2f} {k[1]:5.2f} {k[2]:5.2f} {fwhm:>5} {rms:8.4f} {dpk:12.3f} | " + f"{grp['left']['peak_dx']:+9.2f} {grp['right']['peak_dx']:+9.2f} " + f"{grp['full']['tilt']:10.1f}") + + +def main() -> int: + parser = argparse.ArgumentParser(description="壁反射特性・spot スイープの形の一致を評価") + parser.parse_args() + rows, real_desc = load_all() + if not rows: + print("エラー: 結果がありません") + return 1 + for phase in ("a", "b"): + write_csv(rows, phase) + draw_heatmaps(rows, real_desc, OUTPUT_BASE / "wall_bsdf_heatmap.png") + draw_spot(rows, real_desc, OUTPUT_BASE / "spot_fwhm.png") + print_summary(rows, real_desc) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/scripts/compare_real_vs_sim_metrics.py b/scripts/compare_real_vs_sim_metrics.py new file mode 100644 index 0000000..d3b151e --- /dev/null +++ b/scripts/compare_real_vs_sim_metrics.py @@ -0,0 +1,330 @@ +"""実機 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()) diff --git a/scripts/compare_real_vs_sim_shapes.py b/scripts/compare_real_vs_sim_shapes.py new file mode 100644 index 0000000..f3ff69d --- /dev/null +++ b/scripts/compare_real_vs_sim_shapes.py @@ -0,0 +1,235 @@ +"""実機とシミュの輝度マップの「形」を，全灯と単灯 3 灯（下・右・左）で突き合わせる. + +評価の決め手は集約指標（min 等）ではなく輝度マップの形そのものである（ユーザー方針）． +本スクリプトは 4 構成それぞれについて，実機（DNG 由来の白板 ROI）とシミュを同じ処理 +（シミュ ROI 解像度へリサイズ → メディアン 5 → ガウス σ=4 → ROI 平均で正規化）で並べ， + +- 実機 / シミュ / 差分（シミュ − 実機）のマップ +- 形の一致度: 差分の RMS・最大，SSIM，NCC，方位フーリエ（a1・b1・a2・a4）， + 上下比・左右比，傾き（1 次平面フィットの振幅と方位），お椀（2 次成分）の大きさ + +を出力する．正規化は ROI 平均（=1）にする．単灯は絶対輝度が実機とシミュで揃っていない +（露出・LED 出力差）ので，形だけを比べるにはこれが自然． + +実行方法: + uv run python scripts/compare_real_vs_sim_shapes.py + uv run python scripts/compare_real_vs_sim_shapes.py --sim-set legacy # ρ=0.85 blend の旧レンダ +""" + +from __future__ import annotations + +import argparse +import csv +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.roi_presets import SIM_ROI, extract_roi_xywh # noqa: E402 +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] +OUT_DIR = ROOT / "output" / "real_vs_sim" +REAL_DIR = ROOT / "lightsim" / "luminance_maps" + +CONFIGS = ("full", "bottom", "left", "right") +CONFIG_LABEL = {"full": "全灯 18 灯", "bottom": "下 単灯", "left": "左 単灯", "right": "右 単灯"} +REAL = { + "full": REAL_DIR / "full_MiniTIAS_QM_20260601_124500.npy", + "bottom": REAL_DIR / "bottom_SmTIAS_QM_20260615_141420.npy", + "left": REAL_DIR / "left_SmTIAS_QM_20260615_150300.npy", + "right": REAL_DIR / "right_SmTIAS_QM_20260615_170948.npy", +} +SIM_SETS = { + # TEST_14 で確定した内壁（純 diffuse ρ=0.81）・spp 4096（本スクリプト用にサーバでレンダ） + "rho081": { + "label": "ρ=0.81 純 diffuse・spp 4096", + "full": OUT_DIR / "sim_rho081_spp4096" / "luminance.npy", + "bottom": OUT_DIR / "sim_single_bottom_rho081_spp4096" / "luminance.npy", + "left": OUT_DIR / "sim_single_left_rho081_spp4096" / "luminance.npy", + "right": OUT_DIR / "sim_single_right_rho081_spp4096" / "luminance.npy", + }, + # TEST_12/14 の旧レンダ（blendbsdf ρ=0.85・ω=0.05・α=0.1・spp 1024） + "legacy": { + "label": "ρ=0.85 光沢 5%（旧材質）・spp 1024", + "full": ROOT / "output" / "full_base_fix" / "luminance.npy", + "bottom": ROOT / "output" / "singles_fix" / "bottom" / "luminance.npy", + "left": ROOT / "output" / "singles_fix" / "left" / "luminance.npy", + "right": ROOT / "output" / "singles_fix" / "right" / "luminance.npy", + }, +} +MEDIAN = 5 +SIGMA = 4.0 +DISPLAY_RANGE = (0.85, 1.15) # ROI 平均 = 1 に対する表示レンジ +DIFF_RANGE = 0.08 + + +def load_sim(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 load_real(path: Path, shape: tuple[int, int]) -> np.ndarray: + """実機 ROI（DNG 由来・LSC 適用済み）をシミュ ROI の解像度へリサイズして返す（生）.""" + arr = np.load(path).astype(np.float64) + return resize(arr, shape, anti_aliasing=True, preserve_range=True) + + +def normalized_shape(roi: np.ndarray) -> np.ndarray: + """平滑化（メディアン 5 → ガウス σ=4）して ROI 平均で正規化した形.""" + smooth = gaussian_filter(median_filter(roi, size=MEDIAN), sigma=SIGMA) + return smooth / float(smooth.mean()) + + +def plane_fit(shape_map: np.ndarray) -> tuple[float, float, float]: + """1 次平面（傾き）を当てはめ，(振幅 %, 方位 deg, 2 次成分の振幅 %) を返す. + + 振幅は ROI の半幅あたりの変化量（TEST_14 の「傾き」と同じ定義）．方位は画像座標で + +X（右）を 0 deg，+Y（画像上）を 90 deg とする．2 次成分は傾きを除いた後の + x²+y² 項（お椀）の振幅． + """ + h, w = shape_map.shape + y, x = np.mgrid[0:h, 0:w] + x = (x - w / 2.0) / (w / 2.0) + y = -(y - h / 2.0) / (h / 2.0) # 画像上を +Y にする + design = np.stack([np.ones(x.size), x.ravel(), y.ravel(), + (x**2 + y**2).ravel(), (x**2 - y**2).ravel(), (x * y).ravel()], axis=1) + coeffs, *_ = np.linalg.lstsq(design, shape_map.ravel(), rcond=None) + tilt = float(np.hypot(coeffs[1], coeffs[2])) * 100.0 + azimuth = float(np.degrees(np.arctan2(coeffs[2], coeffs[1]))) + bowl = float(coeffs[3]) * 100.0 + return tilt, azimuth, bowl + + +def shape_metrics(sim_map: np.ndarray, real_map: np.ndarray, sim_roi: np.ndarray, + real_roi: np.ndarray) -> dict[str, float]: + """形の一致度と，各マップの形の記述子を返す.""" + diff = sim_map - real_map + inner = compute_normalized_distance_map(*sim_map.shape) <= MAX_DISTANCE + out = { + "rms_diff": float(np.sqrt(np.mean(diff[inner] ** 2))), + "max_abs_diff": float(np.abs(diff[inner]).max()), + "ssim": float(calc_ssim(sim_map, real_map, + data_range=float(DISPLAY_RANGE[1] - DISPLAY_RANGE[0]))), + "ncc": float(calc_ncc(sim_map, real_map)), + } + for tag, roi, smap in (("sim", sim_roi, sim_map), ("real", real_roi, real_map)): + rings = ring_harmonics(roi, (1, 2, 4)) + for m in (1, 2, 4): + out[f"{tag}_a{m}"], _ = aggregate(rings, f"a{m}") + out[f"{tag}_b{m}"], _ = aggregate(rings, f"b{m}") + h, w = smap.shape + out[f"{tag}_top_bottom"] = float(smap[: h // 3].mean() / smap[-(h // 3):].mean()) + out[f"{tag}_left_right"] = float(smap[:, : w // 3].mean() / smap[:, -(w // 3):].mean()) + out[f"{tag}_peak_mean"] = float(np.percentile(smap, 99.0)) + out[f"{tag}_tilt"], out[f"{tag}_tilt_az"], out[f"{tag}_bowl"] = plane_fit(smap) + return out + + +def cell_title(who: str, cfg: str, m: dict, tag: str) -> str: + """マップの見出し（上下比・左右比・傾き・お椀）.""" + return (f"{who} {CONFIG_LABEL[cfg]}\n上下比 {m[f'{tag}_top_bottom']:.3f}・" + f"左右比 {m[f'{tag}_left_right']:.3f}\n傾き {m[f'{tag}_tilt']:.1f}%@" + f"{m[f'{tag}_tilt_az']:.0f}°・お椀 {m[f'{tag}_bowl']:+.1f}%") + + +def draw(results: dict[str, dict], maps: dict[str, tuple[np.ndarray, np.ndarray]], + sim_label: str, output: Path) -> None: + """3 行（実機 / シミュ / 差分）× 構成数の列のパネル（results にある構成だけ描く）.""" + configs = [c for c in CONFIGS if c in results] + fig, axes = plt.subplots(3, len(configs), figsize=(4.6 * len(configs), 13.5), squeeze=False) + for c, cfg in enumerate(configs): + real_map, sim_map = maps[cfg] + m = results[cfg] + for r, (data, title) in enumerate(( + (real_map, cell_title("実機", cfg, m, "real")), + (sim_map, cell_title("シミュ", cfg, m, "sim")), + )): + ax = axes[r, c] + im = ax.imshow(data, cmap="hot", vmin=DISPLAY_RANGE[0], vmax=DISPLAY_RANGE[1]) + ax.set_title(title, fontsize=8.5) + ax.axis("off") + if c == len(configs) - 1: + fig.colorbar(im, ax=ax, fraction=0.046, label="輝度 / ROI 平均") + ax = axes[2, c] + im = ax.imshow(sim_map - real_map, cmap="coolwarm", vmin=-DIFF_RANGE, vmax=DIFF_RANGE) + ax.set_title(f"差分（シミュ − 実機）\nRMS {m['rms_diff']:.3f}・最大 " + f"{m['max_abs_diff']:.3f}\nSSIM {m['ssim']:.3f}・NCC {m['ncc']:.3f}", + fontsize=8.5) + ax.axis("off") + if c == len(configs) - 1: + fig.colorbar(im, ax=ax, fraction=0.046) + fig.suptitle(f"輝度マップの形の一致 — 実機 vs シミュ（{sim_label}）．" + "全セル同一処理（シミュ解像度へリサイズ・メディアン 5・ガウス σ=4・ROI 平均 = 1）", + fontsize=12) + 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: + """4 構成の形の比較を図と表にする.""" + parser = argparse.ArgumentParser(description="実機 vs シミュの輝度マップの形の比較") + parser.add_argument("--sim-set", choices=tuple(SIM_SETS), default="rho081") + args = parser.parse_args() + sim_set = SIM_SETS[args.sim_set] + results: dict[str, dict] = {} + maps: dict[str, tuple[np.ndarray, np.ndarray]] = {} + for cfg in CONFIGS: + sim_path, real_path = sim_set[cfg], REAL[cfg] + if not sim_path.exists() or not real_path.exists(): + print(f"スキップ {cfg}: {sim_path if not sim_path.exists() else real_path} が無い") + continue + sim_roi = load_sim(sim_path) + real_roi = load_real(real_path, sim_roi.shape) + sim_map, real_map = normalized_shape(sim_roi), normalized_shape(real_roi) + results[cfg] = shape_metrics(sim_map, real_map, sim_roi, real_roi) + maps[cfg] = (real_map, sim_map) + if not results: + print("エラー: 比較できる構成がありません") + return 1 + draw(results, maps, sim_set["label"], OUT_DIR / f"shapes_{args.sim_set}.png") + with (OUT_DIR / f"shapes_{args.sim_set}.csv").open("w", newline="", encoding="utf-8") as fh: + keys = ["config"] + list(next(iter(results.values())).keys()) + writer = csv.DictWriter(fh, fieldnames=keys) + writer.writeheader() + for cfg, m in results.items(): + writer.writerow({"config": cfg, **m}) + print(f"\n=== 形の一致（{sim_set['label']}）===") + print(f"{'構成':>8} {'RMS差':>6} {'最大差':>6} {'SSIM':>6} {'NCC':>6} | " + f"{'上下比 実/シ':>13} {'左右比 実/シ':>13} {'傾き% 実/シ':>13} {'方位 実/シ':>12} " + f"{'お椀% 実/シ':>12} {'b1% 実/シ':>11}") + for cfg, m in results.items(): + print(f"{CONFIG_LABEL[cfg]:>8} {m['rms_diff']:6.3f} {m['max_abs_diff']:6.3f} " + f"{m['ssim']:6.3f} {m['ncc']:6.3f} | " + f"{m['real_top_bottom']:6.3f}/{m['sim_top_bottom']:6.3f} " + f"{m['real_left_right']:6.3f}/{m['sim_left_right']:6.3f} " + f"{m['real_tilt']:6.1f}/{m['sim_tilt']:6.1f} " + f"{m['real_tilt_az']:5.0f}/{m['sim_tilt_az']:5.0f} " + f"{m['real_bowl']:+5.1f}/{m['sim_bowl']:+5.1f} " + f"{m['real_b1'] * 100:+5.1f}/{m['sim_b1'] * 100:+5.1f}") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/scripts/sweep_wall_bsdf_spot.py b/scripts/sweep_wall_bsdf_spot.py new file mode 100644 index 0000000..e983831 --- /dev/null +++ b/scripts/sweep_wall_bsdf_spot.py @@ -0,0 +1,224 @@ +"""壁の反射特性（ω・α）と LED 配光（spot）のスイープ: 単灯の明部のずれを実機に近づける（TEST_20）. + +実機と採用材質（ρ=0.85・ω=0.05・α=0.1）のシミュを全灯・単灯 3 灯で突き合わせると， +同心円方向（お椀）は合うが，単灯で明部が点灯基板の反対側へ寄る量がシミュは実機の半分以下 +だった（compare_real_vs_sim_shapes.py）．原因の候補は (a) 壁の反射が拡散しすぎ（光沢成分が +一次照射を反対側の壁に集中させる），(b) LED の配光がランバートより狭い，の 2 つ． +本スイープは基準シーン（scene_baseline.yaml）から 4 構成（全灯・下・左・右）を作り， + + フェーズ A: 壁 = blendbsdf(diffuse ρ, roughconductor α) の重み ω と α を振る + フェーズ B: 壁 = 採用材質のまま，LED を spot（半値全幅 FWHM を振る）にする + +をサーバでレンダする（評価は手元で build_shape_fit_comparison.py）． + +実行方法（サーバ）: + uv run sm-remote run --name wall --gpu 0 --threads 16 -- \ + python scripts/sweep_wall_bsdf_spot.py --phase all --workers 2 +""" + +from __future__ import annotations + +import argparse +import copy +import json +import multiprocessing +import sys +import time +from concurrent.futures import ProcessPoolExecutor, as_completed +from dataclasses import dataclass +from pathlib import Path + +import numpy as np +import yaml + +sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src")) + +from sm_tias_lightsim.render.pipeline import render_yaml # noqa: E402 +from sm_tias_lightsim.runtime.overrides import resolve_render_variant # noqa: E402 + +ROOT = Path(__file__).resolve().parents[1] +BASELINE_YAML = ROOT / "configs" / "scene_baseline.yaml" +TMP_DIR = ROOT / "configs" / "_tmp_wall_bsdf_spot" +OUTPUT_BASE = ROOT / "output" / "wall_bsdf_spot" + +CONFIGS = ("full", "bottom", "left", "right") +BOARD_PREFIX = {"bottom": "led_bottom", "left": "led_left", "right": "led_right"} + +# フェーズ A: 壁の反射特性 +RHOS = [0.85, 0.81] +OMEGAS = [0.0, 0.05, 0.15, 0.30, 0.50, 0.70] +ALPHAS = [0.02, 0.05, 0.1, 0.3, 0.6] +SPP_A = 2048 +# フェーズ B: LED 配光（半値全幅 FWHM）．採用材質で +SPOT_FWHMS = [30.0, 45.0, 60.0, 90.0, 120.0] +SPOT_MATERIAL = (0.85, 0.05, 0.1) +SPP_B = 4096 +# 面光源（2.8×3.5 mm の矩形・radiance）と同じ全光束になる spot の軸上強度を作る +AREA_MM2 = 2.8 * 3.5 + + +@dataclass(frozen=True) +class Condition: + """1 レンダ条件.""" + + phase: str + config: str + rho: float + omega: float + alpha: float + fwhm: float | None = None + + @property + def name(self) -> str: + mat = f"r{self.rho:.2f}_w{self.omega:.2f}_a{self.alpha:.2f}" + spot = f"_spot{self.fwhm:g}" if self.fwhm is not None else "" + return f"{self.phase}_{self.config}_{mat}{spot}" + + +def conditions(phase: str) -> list[Condition]: + """条件の一覧（フェーズ A → B の順）.""" + out: list[Condition] = [] + if phase in ("a", "all"): + for rho in RHOS: + for omega in OMEGAS: + for alpha in ALPHAS if omega > 0 else [ALPHAS[2]]: + for cfg in CONFIGS: + out.append(Condition("a", cfg, rho, omega, alpha)) + if phase in ("b", "all"): + rho, omega, alpha = SPOT_MATERIAL + for fwhm in SPOT_FWHMS: + for cfg in CONFIGS: + out.append(Condition("b", cfg, rho, omega, alpha, fwhm)) + return out + + +def material_dict(rho: float, omega: float, alpha: float) -> dict: + """内壁材質の YAML 断片（ω=0 なら純 diffuse）.""" + if omega <= 0.0: + return {"type": "diffuse", "reflectance": [rho, rho, rho]} + return { + "type": "blendbsdf", "weight": omega, + "components": [ + {"type": "diffuse", "reflectance": [rho, rho, rho]}, + {"type": "roughconductor", "alpha": alpha, "distribution": "ggx"}, + ], + } + + +def spot_solid_angle(cutoff_deg: float, beam_deg: float) -> float: + """Mitsuba spot（beam_width まで 1，cutoff まで smoothstep で 0）の実効立体角 [sr].""" + theta = np.linspace(0.0, np.radians(cutoff_deg), 4001) + width = max(np.radians(cutoff_deg - beam_deg), 1e-9) + t = np.clip((np.radians(cutoff_deg) - theta) / width, 0, 1) + falloff = np.where(theta <= np.radians(beam_deg), 1.0, t * t * (3 - 2 * t)) + return float(2 * np.pi * np.trapezoid(falloff * np.sin(theta), theta)) + + +def to_spot(light: dict, fwhm: float) -> dict: + """面光源 1 灯を同じ位置・向き・全光束の spot に置き換える（FWHM = cutoff + beam）.""" + cutoff = fwhm * 2.0 / 3.0 + beam = fwhm / 3.0 + radiance = light["radiance"] + radiance = radiance[0] if isinstance(radiance, list) else radiance + flux = np.pi * radiance * AREA_MM2 * 1e-6 # ランバート面光源の全光束 [W] + intensity = flux / spot_solid_angle(cutoff, beam) + return { + "id": light["id"], "type": "spot", "position": light["center"], + "direction": light["direction"], "cutoff_angle": round(cutoff, 4), + "beam_width": round(beam, 4), "intensity": float(intensity), + } + + +def build_yaml(cond: Condition, base: dict) -> Path: + """条件の YAML を書き出す（メッシュは configs/meshes を相対参照）.""" + cfg = copy.deepcopy(base) + for mesh in cfg["geometry"]["meshes"]: + if not Path(mesh["path"]).is_absolute(): + mesh["path"] = "../" + mesh["path"] + cfg["materials"]["printer_white"] = material_dict(cond.rho, cond.omega, cond.alpha) + lights = cfg["lights"] + if cond.config != "full": + lights = [lt for lt in lights if lt["id"].startswith(BOARD_PREFIX[cond.config])] + if cond.fwhm is not None: + lights = [to_spot(lt, cond.fwhm) for lt in lights] + cfg["lights"] = lights + cfg["render"]["spp"] = SPP_B if cond.phase == "b" else SPP_A + cfg["name"] = cond.name + spot_note = f" / spot FWHM={cond.fwhm:g}°" if cond.fwhm else "" + cfg["description"] = (f"TEST_20 {cond.phase.upper()}: {cond.config} / ρ={cond.rho} " + f"ω={cond.omega} α={cond.alpha}{spot_note}") + TMP_DIR.mkdir(parents=True, exist_ok=True) + dest = TMP_DIR / f"scene_{cond.name}.yaml" + dest.write_text(yaml.safe_dump(cfg, allow_unicode=True, sort_keys=False), encoding="utf-8") + return dest + + +def run_job(cond: Condition, scene_yaml: Path, variant: str) -> tuple[str, float, str]: + """1 条件をレンダする（既存は再利用）．戻り値は (名前, 秒, 状態).""" + out_dir = OUTPUT_BASE / cond.name + if (out_dir / "row.json").exists(): + return cond.name, 0.0, "reuse" + start = time.time() + try: + render_yaml(scene_yaml, out_dir, variant) + except Exception as exc: # noqa: BLE001 - 1 条件の失敗でスイープを止めない + return cond.name, time.time() - start, f"fail: {exc}" + row = {"name": cond.name, "phase": cond.phase, "config": cond.config, "rho": cond.rho, + "omega": cond.omega, "alpha": cond.alpha, "fwhm": cond.fwhm, + "spp": SPP_B if cond.phase == "b" else SPP_A, "variant": variant, + "elapsed_s": round(time.time() - start, 1)} + (out_dir / "row.json").write_text(json.dumps(row, ensure_ascii=False, indent=1), + encoding="utf-8") + return cond.name, time.time() - start, "ok" + + +def main() -> int: + """スイープを実行する.""" + parser = argparse.ArgumentParser(description="壁の反射特性と spot 配光のスイープ") + parser.add_argument("--phase", choices=("a", "b", "all"), default="all") + parser.add_argument("--workers", type=int, default=1) + parser.add_argument("--prepare-only", action="store_true") + parser.add_argument("--list", action="store_true") + args = parser.parse_args() + conds = conditions(args.phase) + if args.list: + for c in conds: + print(c.name) + print(f"{len(conds)} 条件") + return 0 + base = yaml.safe_load(BASELINE_YAML.read_text(encoding="utf-8")) + scenes = {c.name: build_yaml(c, base) for c in conds} + log_path = ROOT / "output" / "_sweep_wall_bsdf_spot.log" + log_path.parent.mkdir(parents=True, exist_ok=True) + log_file = log_path.open("a", encoding="utf-8") + + def log(message: str) -> None: + print(message, flush=True) + log_file.write(message + "\n") + log_file.flush() + + log(f"=== 壁反射特性 / spot スイープ: {len(conds)} 条件（phase={args.phase}）===") + if args.prepare_only: + return 0 + variant = resolve_render_variant(str(base["render"]["variant"])) + done = 0 + if args.workers > 1: + context = multiprocessing.get_context("spawn") + with ProcessPoolExecutor(max_workers=args.workers, mp_context=context) as pool: + futures = {pool.submit(run_job, c, scenes[c.name], variant): c for c in conds} + for fut in as_completed(futures): + name, sec, status = fut.result() + log(f" {name}: {status} ({sec:.0f}s)") + done += status == "ok" + else: + for c in conds: + name, sec, status = run_job(c, scenes[c.name], variant) + log(f" {name}: {status} ({sec:.0f}s)") + done += status == "ok" + log(f"=== 完了: 新規 {done} / {len(conds)} 条件 ===") + log_file.close() + return 0 + + +if __name__ == "__main__": + sys.exit(main())