diff --git a/CLAUDE.md b/CLAUDE.md index 8f96628..0fb74ec 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -23,7 +23,7 @@ ## 開発進捗 -最新: フェーズ 7（実装）— Step 8-14 の物理モデル精緻化と並行し，**[TEST_13](docs/06_TEST/TEST_13_物理モデル健全性の一括監査.md) 物理モデル健全性の一括監査**を実施．TEST_12 と同型の「無警告で結果だけが歪む」故障を横断検査し，**絶対放射量が積分球の解析解と 0.2% 一致**（単位系・反射率・光源パワー・センサ校正は健全）・シェーディング法線と幾何法線の一致（逸脱 0.044 度）・メッシュ法線の向き（裏面ヒット 0）を確認．**修正 3 件**: 片面 BSDF の脆弱性（法線裏返りで無警告の黒化 → `twosided` ラップ．実シーンには無影響）／**LLVM 初期化の順序バグでレンダテスト 8 件が長期間全滅していた**問題／ TEST_08 軸取り違えの回帰テスト欠落．**要判断**: `max_depth=16` は収束値の 90.3% しかエネルギーを拾えていない（CoV への影響は +0.35% で形状指標の過去結論は無影響．設定変更は LED radiance 再校正を伴うため保留，警告ログのみ実装）．**可視化の復旧**: `sm-visualize` の起動不能は VTK の破損ではなく Windows Smart App Control による DLL ブロックと判明し，`vtk==9.4.2` への固定で解消（可視化テスト 54 件が全て通過）．再破損を即検知する import 監視テストを通常実行に常設．詳細な経緯・設計判断・失敗パターンは [docs/PROGRESS.md](docs/PROGRESS.md)． +最新: フェーズ 7（実装）— 面法線修正（TEST_12）後の条件で内壁材質 3 軸を全数スイープした**[TEST_14](docs/06_TEST/TEST_14_材質パラメータ再スイープと残差の構造分解.md) 材質パラメータ再スイープと残差の構造分解**を実施．**内壁材質は ρ ≈ 0.81 の純 diffuse に確定**（振幅系 4 指標の実機交差点 0.80〜0.82 と，残差のお椀成分のゼロ交差 0.81 が独立に一致）．**ω・α はほぼ無効で blendbsdf は不要**＝調整パラメータが 1 つに減った．**最重要**: 全 60 条件で NCC が 0.81〜0.86 に張り付き，**形は材質パラメータでは原理的に動かない**ことを網羅探索で実証．残差を `傾き + お椀` に分解すると ρ が動かすのはお椀だけで，**傾き（約 6%・方位 −140 deg）は不変**，両者を除いた残りは約 1%＝ほぼ MC ノイズ．既存単灯に同じ分解を当てると傾きは bottom 6.00%・left 7.17% に対し **right だけ 1.84%** ＝カメラ由来ではなく **bottom/left のモデル化誤差**と確定し，次の主題に設定．**未実施の前提**: ノイズフロア（seed 違い 2 本）未測定のため上位 4 条件は横並びで主張は「ρ ≈ 0.80〜0.82」まで．詳細な経緯・設計判断・失敗パターンは [docs/PROGRESS.md](docs/PROGRESS.md)． 実装済みモジュール: - `src/sm_tias_lightsim/uniformity.py` — Rec.709 輝度変換 + CoV/std/max-min（SmTIAS-Evaluation コピー） @@ -146,3 +146,4 @@ - シミュレーション実行記録: docs/06_TEST/TEST_02_シミュレーション実行記録.md（output/ 各サブディレクトリの索引．Step との対応・シーン/評価 config・主要指標・再現コマンド・削除推奨） - 輝度マップ形状追跡: docs/06_TEST/TEST_03_輝度マップ形状追跡.md（数値とは独立に「形」を実機に近づける観点．比較パネル画像と LSC モデル化の構造的意義） - 物理モデル健全性の一括監査: docs/06_TEST/TEST_13_物理モデル健全性の一括監査.md（「無警告で結果だけが歪む」故障モードの横断監査．解析解による絶対放射量検証・法線監査・打ち切りバイアス・常設した再発防止テスト） +- 材質パラメータ再スイープと残差の構造分解: docs/06_TEST/TEST_14_材質パラメータ再スイープと残差の構造分解.md（面法線修正後の ρ×ω×α 全 60 条件スイープ．ρ ≈ 0.81 純 diffuse の確定と，材質では直せない傾き成分の分離） 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 e5f049a..96df4da 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" @@ -75,6 +75,8 @@ | `output/sphere_r100_loft_disk/` | TEST_11 追補 | `configs/scene_sphere_r100_loft_disk.yaml`（ロフト一体キャビティ＋同面積円板単灯 r=1.7662mm＝完全軸対称，spp=1024・d=16） | `configs/eval_baseline_dng_target_shape.yaml` | CoV 0.1205 ほか全指標が矩形版と MC ノイズ内で一致 — [TEST_11 §追補](TEST_11_球体筐体実験.md) 参照 | TEST_11 追補（円形光源）．発光面形状の不感を確認．**以後の球実験の正準構成** | | `output/sphere_r100{,_fit,_loft}/geometry/` | TEST_11 | — | — | レンダ前 3D 配置確認のスクリーンショット（`--no-bake` レイアウト・LED 4 アングル・path ベイク） | 幾何変更時のレンダ前 3D 確認ワークフロー（CLAUDE.md 必須ルール）の記録．fit 版はレンダ未実施（ロフト版に置換） | | `output/full_base_fix/` | TEST_12 | `configs/_tmp_fn_full_base.yaml`（scene_baseline_blend の w=0.05/α=0.1/d=16 版＝TEST_10 full_base の再構成．修正前コードで保存済み結果と一致検証済み） | `configs/eval_baseline_dng_target_shape.yaml` | **CoV 0.0610 / 勾配 6.83% / C/P 1.073**（面法線修正後）— [TEST_12](TEST_12_シェーディング法線アーティファクトの発見と修正.md) 参照 | **面法線修正後の現行筐体全灯の新基準**．実機（CoV 0.0695）と数値レベルで整合 | +| `output/bsdf_grid_fix/` (60 ディレクトリ) | TEST_14 | `configs/_tmp_fn_full_base.yaml` の `printer_white` のみを `scripts/sweep_bsdf_grid_fix.py` で ρ ∈ {0.75, 0.78, 0.82, 0.85, 0.88, 0.90} × ω ∈ {0, 0.05, 0.15, 0.30} × α ∈ {0.1, 0.3, 0.6} に書き換えた 60 種（d=16, spp=1024 据え置き．ω=0 は純 diffuse） | `configs/eval_baseline_dng_target.yaml` | 最良 ρ0.82 ω0: **CoV 0.0669 / 勾配 8.38% / C/P 1.0915 / 動径 0.862**（実機 0.0695 / 8.35% / 1.0911 / 0.860）— [TEST_14](TEST_14_材質パラメータ再スイープと残差の構造分解.md) 参照 | **面法線修正後の材質パラメータ確定**．ρ ≈ 0.81 純 diffuse，ω・α は無効．全 60 条件で NCC が 0.81〜0.86 に張り付き形は不動 | +| `output/bsdf_grid_fix/summary.csv` | TEST_14 | — | — | 60 条件集約 CSV（ρ/ω/α + 指標 + 誤差成分 + score + 打ち切り損失） | `scripts/build_bsdf_grid_panel.py`・`build_bsdf_grid_visual.py` の入力．重み付けを変えた再ランキングにも使える | | `output/sphere_r100_loft_disk_fix/` | TEST_12 | `configs/scene_sphere_r100_loft_disk.yaml`（面法線修正後の再レンダ） | 同上 | CoV 0.0628 / angular_ratio 0.005（完全同心円） | 面法線修正後の球体筐体の正準結果．現行筐体とほぼ同等＝球形化の優位性消滅 | | `output/_slide_panels/sphere_fix_{shape_roi,relative}.png` | TEST_12 | — | — | 現行筐体 fix／球ロフト円形 fix／実機の 1×3 パネル（修正後のみで比較） | TEST_12 視覚比較・発表スライド素材（`scripts/build_sphere_panel.py` 現行版が生成） | | `output/singles_fix/{bottom,left,right}/` | TEST_12 追補 | `configs/_tmp_fix_{bottom,left,right}.yaml`（検証済み `_tmp_fn_full_base.yaml` から灯グループ抽出．blend w0.05_a0.1・d=16・spp=1024） | 形指標のみ（`scripts/build_singles_fix_summary.py`） | 灯間 NCC（平滑）が修正前 0.99 均等 → **修正後 0.85/0.93/0.73 の非対称** — [TEST_12 §追補](TEST_12_シェーディング法線アーティファクトの発見と修正.md) 参照 | TEST_09 中核結論「3 灯同形」の撤回根拠．単灯 spp=1024 はノイズ 5.3% > 構造 4% のため raw 指標はノイズ支配（以後の単灯は spp≥4096 推奨） | diff --git "a/docs/06_TEST/TEST_14_\346\235\220\350\263\252\343\203\221\343\203\251\343\203\241\343\203\274\343\202\277\345\206\215\343\202\271\343\202\244\343\203\274\343\203\227\343\201\250\346\256\213\345\267\256\343\201\256\346\247\213\351\200\240\345\210\206\350\247\243.md" "b/docs/06_TEST/TEST_14_\346\235\220\350\263\252\343\203\221\343\203\251\343\203\241\343\203\274\343\202\277\345\206\215\343\202\271\343\202\244\343\203\274\343\203\227\343\201\250\346\256\213\345\267\256\343\201\256\346\247\213\351\200\240\345\210\206\350\247\243.md" new file mode 100644 index 0000000..7493a51 --- /dev/null +++ "b/docs/06_TEST/TEST_14_\346\235\220\350\263\252\343\203\221\343\203\251\343\203\241\343\203\274\343\202\277\345\206\215\343\202\271\343\202\244\343\203\274\343\203\227\343\201\250\346\256\213\345\267\256\343\201\256\346\247\213\351\200\240\345\210\206\350\247\243.md" @@ -0,0 +1,227 @@ +# 材質パラメータ再スイープと残差の構造分解 (Material Parameter Re-sweep and Residual Decomposition) + +面法線アーティファクト（[TEST_12](TEST_12_シェーディング法線アーティファクトの発見と修正.md)）の修正後， +内壁の材質パラメータ 3 軸 (ρ / ω / α) を全数スイープして実機に最も近い条件を確定し， +さらに残差を低次成分に分解して「材質で直せる誤差」と「直せない誤差」を分離した（2026-08-30〜31）． + +関連: + +- 前提となる修正: [TEST_12_シェーディング法線アーティファクトの発見と修正](TEST_12_シェーディング法線アーティファクトの発見と修正.md) +- 形の追跡方針: [TEST_03_輝度マップ形状追跡](TEST_03_輝度マップ形状追跡.md) +- 実行記録: [TEST_02_シミュレーション実行記録](TEST_02_シミュレーション実行記録.md) + +## 背景・動機 (Motivation) + +反射率スイープ（`output/sweep_reflectance/`）と blend BSDF スイープ（`output/blend_sweep/`）は +いずれも面法線修正**前**に実施しており，当時はシミュの不均一性の 6 割超がシェーディング法線 +アーティファクトだった．したがって「どの材質パラメータが実機に近いか」という結論は +**アーティファクトへのフィッティングになっていた可能性が高く，引き直しが必要**だった． + +加えて TEST_12 で実機とシミュが数値レベルで整合するようになったため，**指標を使った比較が +初めて意味を持つ**状態になった．ここで材質パラメータを確定させれば，以後の設計検討 +（未製作筐体の予測）の土台になる． + +## 実験条件 (Setup) + +固定条件は既存レンダと完全に同一に保ち，`printer_white` の材質定義だけを書き換えた +（テンプレート: `configs/_tmp_fn_full_base.yaml`）． + +| 項目 | 値 | +| --- | --- | +| レンダ | spp=1024 / max_depth=16 / 512² / llvm_ad_rgb / path / orthographic | +| 幾何・光源・センサ | 変更なし（現行筐体・18 灯・`ortho_scale` 50×52） | +| 評価 | `configs/eval_baseline_dng_target.yaml`（linear + mean アンカー 207.82） | +| 比較対象 | 実機全灯 `lightsim/luminance_maps/full_MiniTIAS_QM_20260601_124500.npy` | + +スイープ軸（ω=0 は α に依存しないため α は 1 値のみ．計 60 条件）: + +| 軸 | 意味 | 値 | +| --- | --- | --- | +| ρ | diffuse 成分の反射率 | 0.75 / 0.78 / 0.82 / 0.85 / 0.88 / 0.90 | +| ω | blendbsdf の roughconductor 混合比（0 なら純 diffuse） | 0 / 0.05 / 0.15 / 0.30 | +| α | roughconductor の GGX 粗さ | 0.1 / 0.3 / 0.6 | + +実行は無人 6 時間枠で行い，58 ジョブを 3.98 時間で完走した（1 ジョブ 179〜332 秒． +基準点 ρ=0.85/ω=0.05/α=0.1 は `output/full_base_fix/` を再利用し再レンダしていない）． + +### 指標の扱い + +形の一致度は **平滑化（median 5px + gaussian σ=4）後の NCC 対実機**を用いる． +TEST_12 追補で確認したとおり生マップの NCC は MC ノイズに支配される（単灯でノイズ 5.3% +> 構造振幅 3.8〜4.4%）ためである．振幅系は CoV / 中心周辺比 / 勾配 / 動径 min-max の +実機との相対乖離を用い，総合スコアは形と振幅を半々で重み付けした（`sweep_bsdf_grid_fix.py`）． + +## 結果 1 — ρ が支配的，ω・α はほぼ無効 (Results: ρ Dominates) + +スコア上位 5 条件: + +| 条件 | score | CoV | C/P | 勾配 [%] | 動径 min/max | NCC | 打切損失 [%] | +| --- | --- | --- | --- | --- | --- | --- | --- | +| **ρ0.82 ω0（純 diffuse）** | **0.0789** | 0.0669 | 1.0915 | 8.38 | 0.862 | 0.857 | 5.1 | +| ρ0.78 ω0.05 α0.1 | 0.0801 | 0.0701 | 1.0925 | 8.47 | 0.858 | 0.853 | 3.0 | +| ρ0.82 ω0.05 α0.3 | 0.0854 | 0.0665 | 1.0899 | 8.25 | 0.866 | 0.856 | 6.0 | +| ρ0.82 ω0.05 α0.6 | 0.0873 | 0.0681 | 1.0934 | 8.54 | 0.866 | 0.854 | 6.0 | +| ρ0.85 ω0.15 α0.6 | 0.0900 | 0.0672 | 1.0902 | 8.27 | 0.867 | 0.845 | 12.9 | +| **実機 全灯** | — | **0.0695** | **1.0911** | **8.35** | **0.860** | 1.000 | — | + +純 diffuse（ω=0）系列の ρ 応答: + +| ρ | CoV | 勾配 [%] | C/P | 動径 min/max | NCC | score | +| --- | --- | --- | --- | --- | --- | --- | +| 0.75 | 0.0761 | 10.33 | 1.1152 | 0.831 | 0.852 | 0.1749 | +| 0.78 | 0.0723 | 9.52 | 1.1053 | 0.844 | 0.852 | 0.1301 | +| **0.82** | **0.0669** | **8.38** | **1.0915** | **0.862** | 0.857 | **0.0789** | +| 0.85 | 0.0626 | 7.54 | 1.0816 | 0.874 | 0.860 | 0.1204 | +| 0.88 | 0.0584 | 6.79 | 1.0728 | 0.888 | 0.856 | 0.1655 | +| 0.90 | 0.0552 | 6.19 | 1.0660 | 0.898 | 0.858 | 0.1977 | +| 実機 | 0.0695 | 8.35 | 1.0911 | 0.860 | 1.000 | — | + +読み取れること: + +1. **振幅系 4 指標が単調かつ独立に同じ ρ を指す**．実機値との交差点は CoV → 0.80， + 勾配 → 0.815，C/P → 0.82，動径 min-max → 0.815．**偶然では揃わない一致**であり， + ρ が真に効いているパラメータであることの裏付けになる． +2. **ω・α はほぼ効かない**．ρ を固定して ω・α を全振りしても CoV の変化幅は 0.008〜0.011 + にとどまるのに対し，ρ 単独では 0.0542〜0.0797 と 3 倍動く． +3. **例外は ω=0.30・α=0.1（強い準鏡面）のみ**で，これだけ全指標で外れ値になる + （勾配 3.5〜4.8%，NCC 0.61〜0.70）．壁が鏡に近づく別レジームであり，実機の白色樹脂とは + かけ離れている． +4. **結論として blendbsdf は不要**．純 diffuse でよく，調整対象パラメータが ρ 1 つに減る． + +図: `output/_slide_panels/bsdf_grid_response.png`（指標 × ρ の応答曲面）， +`output/_slide_panels/bsdf_grid_contact_sheet.png`（全 60 条件の輝度マップ一覧）． + +## 結果 2 — 形は材質パラメータでは動かない (Shape Is Insensitive) + +**60 条件すべてで NCC は 0.6118〜0.8599**，外れ値の ω=0.30/α=0.1 系列を除けば +**0.8135〜0.8599** に収まり，基準点（0.8544）からほとんど改善しない． +応答曲面の NCC パネルはほぼ水平線になる． + +コンタクトシートで目視しても **60 セルが実質すべて同じ絵**（左右対称の中央ブロブ＋上端の +暗い帯）であり，実機の「右下が明るく左上隅が暗い斜め構造」はどの条件でも現れない． + +振幅系が ρ で 3 倍動くのに形の相関が動かない ── **残る形の差は材質パラメータでは +原理的に埋まらない**ことが，網羅探索によって実証された． + +これは [TEST_03](TEST_03_輝度マップ形状追跡.md) の主題「数値の縮小 ≠ 形の改善」が， +面法線修正後の新しい条件でも成立することを意味する． + +## 結果 3 — 残差の構造分解 (Residual Decomposition) + +指標だけでは最終判断にならないため，実機との残差マップ +（mean 正規化後の `(sim - real) / mean(real)` [%]）を直接観察し， +さらに `定数 + 傾き(X, Y) + お椀(R²)` で最小二乗分解した． + +純 diffuse 系列: + +| ρ | 残差 RMS [%] | 傾き [%] | 方位 [deg] | お椀 [%] | 分解後の残り [%] | +| --- | --- | --- | --- | --- | --- | +| 0.75 | 3.58 | 5.79 | −146.1 | **−1.77** | 1.01 | +| 0.78 | 3.63 | 6.00 | −143.7 | **−0.90** | 0.99 | +| 0.82 | 3.72 | 6.19 | −140.7 | **+0.27** | 0.99 | +| 0.85 | 3.83 | 6.36 | −138.2 | +1.12 | 0.97 | +| 0.88 | 4.00 | 6.56 | −136.3 | +1.91 | 0.98 | +| 0.90 | 4.12 | 6.69 | −134.9 | +2.46 | 0.94 | + +全 60 条件での広がり: 傾き 5.79〜8.74%，方位 −146.1〜−127.1 deg，お椀 −2.12〜+4.92%． + +読み取れること: + +1. **ρ が制御しているのはお椀（動径）成分だけ**であり，**ρ ≈ 0.81 でゼロ交差**する． + 振幅系 4 指標から読んだ 0.80〜0.82 と独立に一致する． +2. **傾き成分（約 6%・方位 −140 deg）はどの ρ でもほぼ一定**．光沢を混ぜると + 最大 8.74% まで**悪化**する．材質パラメータでは手が届かない． +3. **傾きとお椀を除くと残りは約 1%** ＝ ほぼ MC ノイズ．モデル誤差はこの 2 成分で + ほぼ説明し尽くされている． +4. したがって **残差 RMS で ρ を選んではならない**．RMS は支配項である傾き（直せない誤差）を + 測ってしまうため ρ=0.75 が最良に見えるが，これは誤り．**お椀成分がゼロになる ρ** が + 正しい選択基準である． + +図: `output/_slide_panels/bsdf_grid_rho_series.png`（上段 輝度マップ / 下段 残差マップ）， +`output/_slide_panels/bsdf_grid_top_residual.png`（上位条件の残差マップ）． + +### 単灯での確認 — 傾きは基板ごとに違う (Per-Board Check) + +既存の単灯レンダ（`output/singles_fix/`，ρ=0.85・ω=0.05・α=0.1．全灯と同一条件）に +同じ分解を適用した（追加レンダなし）: + +| 灯 | 残差 RMS [%] | 傾き [%] | 方位 [deg] | お椀 [%] | 分解後の残り [%] | +| --- | --- | --- | --- | --- | --- | +| 全灯 | 3.99 | 6.55 | −137.4 | 1.77 | 0.99 | +| bottom | 3.98 | 6.00 | −143.0 | 3.56 | 1.23 | +| left | 4.55 | 7.17 | −151.7 | 3.16 | 1.30 | +| **right** | **2.30** | **1.84** | **−62.0** | 1.43 | 1.94 | + +**傾きは全灯共通の性質ではなく，基板ごとに異なる**．bottom と left は 6〜7% で方位も近い +（−143 / −152 deg）のに対し，**right だけ 1.84% とほぼ傾きがない**（傾きが小さいため方位 +−62 deg 自体はほとんど意味を持たない）． + +これは原因の切り分けに直結する．**仮に傾きがカメラ由来（LSC 補正残差）であれば，同じカメラ・ +同じ ROI を使う以上 4 条件すべてに同一の傾きが乗るはずである**．そうなっていない以上， +**傾きの主要部分はカメラのアーティファクトではなく，bottom と left の基板のモデル化誤差**と +判断できる．right は既にほぼ整合している． + +なお，実機単灯の ROI mean は bottom 57.75 / left 61.46 / right 63.55（露出・ISO 固定の +同一セッション）で **右 > 左 > 下 の 10% の広がり**があるのに対し，シミュ（`singles_fix`）は +6.539 / 6.559 / 6.575 と **0.5% 以内でほぼフラット**であり，基板間の出力差が未モデル化である． +ただし一律のゲイン差だけでは「right だけ傾きが小さい」ことを説明できないため， +位置・向き・配光のいずれかが bottom / left で外れている可能性が高い． + +## 結論 (Conclusions) + +1. **内壁材質は ρ ≈ 0.81 の純 diffuse（ω=0）を採用する**．振幅系 4 指標と残差のお椀成分が + 独立に同じ値を指しており，材質パラメータの探索はここで決着とみなしてよい． + blendbsdf は不要． +2. **残る誤差は約 6% の傾き成分**（方位 −140 deg ＝ 左上が明るすぎ・右下が暗すぎ）であり， + 材質では直せない．これが次の主題になる． +3. **傾きは bottom / left の基板に固有**であり，right はほぼ整合している．カメラ由来ではなく + モデル化誤差である． + +## 限界・注意点 (Limitations) + +- **ノイズフロアを測っていない**．上位 4 条件（score 0.0789〜0.0873）は実質横並びであり， + 「ρ=0.82 が ρ=0.78 より良い」とは言えない．正しい主張は「ρ ≈ 0.80〜0.82」まで． + 同一条件を異なる seed で 2 回レンダし，シミュ同士の指標差を測れば各指標の到達可能上限が + 出る（未実施）． +- **`max_depth=16` 据え置きによる打ち切り損失**が ρ とともに増える（ρ=0.82・ω=0 で 5.1%， + ρ=0.90 で 20.6%）．採用条件は許容範囲だが，高 ρ 側の条件は打ち切り誤差と交絡している． + なお blendbsdf の roughconductor 成分は `material="none"`（反射率 1.0 の完全反射体）のため， + ω が大きいほど実効反射率が上がり打ち切り損失も増える（ρ0.88・ω0.30 で 26.8%）． +- **総合スコアの重み付け（形と振幅で半々）は恣意的**である．全成分を + `output/bsdf_grid_fix/summary.csv` に残してあるため，別の重みで並べ替え直せる． +- **実機データのセッションが揃っていない**．全灯は 2026-06-01（MiniTIAS），単灯 3 灯は + 2026-06-15（SmTIAS）で撮影セッションが別．単灯どうしの比較は露出・ISO 固定の同一 + セッションなので有効だが，全灯と単灯の絶対振幅を横並びにするには露出正規化が要る． +- **比較対象が実機 1 台・全灯 1 枚**である．ρ の妥当性を単灯 3 条件で独立検証していない． + +## 次にやること (Next Steps) + +1. **単灯 3 本を ρ=0.81・純 diffuse で引き直す**（3 レンダ約 10 分）．ρ の独立検証と， + bottom / left の傾きが ρ 非依存で残ることの確認が同時に取れる． +2. **傾き成分の原因切り分け**: LED 基板ごとの出力・位置・向き・配光のどれが効くかを + 単灯で個別に検証する． +3. **ノイズフロアの測定**（seed 違い 2 本のレンダ）．以後すべての一致度指標の解釈の前提になる． +4. `configs/scene_baseline.yaml` の既定材質を ρ=0.81 純 diffuse に更新するかの判断． + +## 再現コマンド (Reproduction) + +```powershell +# スイープ（無人実行．時間予算で自動打ち切り，1 ジョブ完了ごとに summary.csv を更新） +uv run python scripts/sweep_bsdf_grid_fix.py --budget-hours 5.5 + +# 指標の可視化（応答曲面 + 上位条件の形比較） +uv run python scripts/build_bsdf_grid_panel.py --top 3 + +# 輝度マップそのものの可視化（一覧 / ρ 系列 + 残差 / 上位条件の残差） +uv run python scripts/build_bsdf_grid_visual.py --top 4 +``` + +出力: `output/bsdf_grid_fix/summary.csv`（60 条件）・各条件の `luminance.npy` と `metrics.csv`・ +ログ `output/_sweep_bsdf_grid_fix.log`・図 `output/_slide_panels/bsdf_grid_*.png`． + +## 参照 (References) + +- [TEST_12_シェーディング法線アーティファクトの発見と修正](TEST_12_シェーディング法線アーティファクトの発見と修正.md): 本再スイープの前提となる修正 +- [TEST_13_物理モデル健全性の一括監査](TEST_13_物理モデル健全性の一括監査.md): 打ち切りバイアスの見積もり式と絶対放射量の検証 +- [TEST_03_輝度マップ形状追跡](TEST_03_輝度マップ形状追跡.md): 「数値と形は別軸」という評価方針 +- [TEST_02_シミュレーション実行記録](TEST_02_シミュレーション実行記録.md): 出力ディレクトリの索引 diff --git a/docs/PROGRESS.md b/docs/PROGRESS.md index a9eb8e6..433338f 100644 --- a/docs/PROGRESS.md +++ b/docs/PROGRESS.md @@ -80,3 +80,4 @@ - **TEST_12 追補 2 — ロフト w1.4 の修正後再検証で TEST_10 の判定を確定 (2026-08-18)**: ロフト w1.4（全灯 1＋単灯 3）を修正版で再レンダ．**④「全灯に実質不感」は維持**（CoV 0.0610→0.0605，全指標差 ≈1%＝実機 w1.4 測定への予測「従来筐体とほぼ同じ」も不変）．**③「対称拡大は 3 灯を同形化」も維持・明確化**（灯間 NCC 平滑 0.85/0.93/0.73→0.91/0.96/0.84 で実機の非対称から遠ざかる）．一方 **①「bottom 単灯の平ら化＝幾何レバー初実証」は撤回相当に弱まる**（平滑 p/m 変化 -0.9% のみ．TEST_10 の raw -10% はアーティファクト期の見かけが大半．修正後 raw p/m はヘビーテイルノイズの最大値ゆらぎで不採用）．②「left/right 不動」維持．TEST_10 冒頭に BEFORE FIX 注記（結論の生死を明記）．成果物: `output/full_w1p4_fix/`・`output/singles_fix/*_w1p4/`，summary/pairwise CSV と 3×4 パネル更新，記録 [TEST_12 §ロフト再検証](06_TEST/TEST_12_シェーディング法線アーティファクトの発見と修正.md)． - **TEST_13 物理モデル健全性の一括監査 — 潜在バグ 3 件を修正し再発防止テストを常設 (2026-08-30)**: TEST_12 と同じ「エラーも警告も出ないまま結果だけが歪む」故障モードが他に残っていないかを横断監査（動機: 同型のバグが TEST_08 でも起きており，そのとき TEST_07 が誤出力を「物理的に正しい」と一度正当化している＝3 度目を防ぐ）．**健全と確認**: シェーディング法線 == 幾何法線（実測逸脱 0.044 度）・メッシュ法線の向き（裏面ヒット 0・全 15 STL 巻き順一貫）・**絶対放射量（積分球の解析解と 0.2% 一致＝単位系/反射率/光源パワー/センサ校正が正しい）**・面光源の軸マッピング．**修正した問題 3 件**: ① 片面 BSDF の脆弱性（法線が裏返ると無警告で真っ黒．`twosided` ラップで恒久修正．実シーンでは mean +0.015%＝無影響を確認）／ ② **LLVM 初期化の順序バグで `slow` レンダテスト 8 件が長期間全滅**（`scene/lights.py` がモジュール import 時に mitsuba を読むため `runner.py` の自動探索が手遅れ．開発機は環境変数が永続設定済みで症状が隠れていた）／ ③ TEST_08 の軸取り違えに実効的な回帰テストが無かった．**要判断の発見**: `max_depth=16` は収束値の **90.3%** しかエネルギーを拾えていない（間接照明のみの構成のため．ただし CoV への影響は +0.35% で形状指標の過去結論は無影響．設定変更は LED radiance 再校正を伴うため未実施，警告ログのみ実装）．**未解決**: VTK 破損で `sm-visualize` と可視化テスト 39 件が実行不能（依存関係の修復が必要）．設計判断: 対称性テストは負のコントロールが発火せず（立方体では平滑法線でも象限比 1.002 vs 1.003）検出力を実証できなかったため**常設テストに含めない**（検出できない検査は「検査済み」という誤った安心を与える）．積分球テストは 1.5 秒と軽量なので `slow` を付けず通常実行に含めた（`slow` 隔離が ② の全滅を見逃した原因のため）．記録 [TEST_13](06_TEST/TEST_13_物理モデル健全性の一括監査.md)． - **sm-visualize の起動不能を解消 — 原因は VTK 破損ではなく Smart App Control (2026-08-30)**: TEST_13 で唯一未解決だった 3D 可視化の全滅（可視化テスト 39 件・`sm-visualize` 実行不能）を解決．**当初「VTK のインストール破損」と誤診していたが，真因は Windows Smart App Control（強制モード）が vtk 9.6.1 の DLL 209 個のうち 14 個をレピュテーション判定でブロックしていたこと**．署名の有無は判定基準ではなく（通った DLL も弾かれた DLL も未署名），`.pyd` は全て正常に存在し `uv.lock` も 5/17 から不変だったため，`uv sync` や同一バージョン再インストールでは直らない種類だった．誤診の原因は最初に観測したメッセージが `No module named vtkmodules.vtkCommonCore` という二次症状だったこと＝**ネイティブ拡張の ImportError は表示メッセージを額面どおり受け取らない**という教訓．対処は `vtk==9.4.2` への固定（より長く配布されレピュテーションが確立したバージョン．pyvista 0.48.0 の制約 `>=9.2.2,<9.7.0,!=9.4.1` 内）で，全 147 モジュールがブロック 0 件で読み込め可視化テスト 54 件が全通過．設計判断: Smart App Control の無効化は**一度オフにすると Windows 再インストールなしに戻せない**ため不採用（可視化ツール 1 本のために恒久的にマルウェア保護を手放すのは釣り合わない）．WSL2 への切り出しは環境二重管理の負債を負うため次善策として保留．**再発防止**: 破損が 8/17〜8/30 の約 2 週間気づかれなかった原因は可視化テストが全て `slow` 隔離されていたことで，TEST_13 検査 7（LLVM 初期化バグでレンダテストが長期間全滅）とまったく同じ構図．import 可否だけを見る軽量テスト（1 秒未満）を `slow` 無しで常設し，`pyproject.toml` の固定にも理由をコメントで残した（Dependabot / `/deps-update` による静かな再破損を防ぐため）．判断は `/consult`（advisor: Fable）に諮り，A（ダウングレード）→ B（WSL2）の順・C（SAC 無効化）不採用・import 監視テスト常設という助言に沿って実施．記録 [TEST_13 検査 8](06_TEST/TEST_13_物理モデル健全性の一括監査.md)． +- **TEST_14 材質パラメータ再スイープと残差の構造分解 — ρ を確定し「材質では直せない誤差」を分離 (2026-08-31)**: 面法線修正（TEST_12）前に実施した反射率・blend BSDF スイープはアーティファクトへのフィッティングだった可能性が高いため，ρ×ω×α の全 60 条件を修正後条件（spp=1024・d=16 据え置き）で引き直した（無人 6 時間枠・58 ジョブ 3.98h）．**結論は ρ ≈ 0.81 の純 diffuse**で，振幅系 4 指標（CoV/勾配/C-P/動径）の実機交差点が 0.80〜0.82 に揃い，残差のお椀成分のゼロ交差 0.81 とも独立に一致した．**ω・α はほぼ無効**（ρ 固定で CoV 変化幅 0.008〜0.011 に対し ρ 単独では 3 倍動く）で **blendbsdf は不要**＝調整パラメータが 1 つに減る．**最重要の発見**: 60 条件すべてで NCC が 0.81〜0.86 に張り付き（外れ値の準鏡面系列を除く），**形は材質パラメータでは原理的に動かない**ことを網羅探索で実証した．残差を `傾き + お椀` に分解すると **ρ が動かすのはお椀だけ・傾き（約 6%・方位 −140 deg）は不変**で，両者を除いた残りは約 1%＝ほぼ MC ノイズ．**方法論上の学び**: 残差 RMS は支配項である傾き（直せない誤差）を測ってしまうため ρ の選択基準に使ってはならない（RMS だけ見ると ρ=0.75 が最良に見える）．**原因の切り分け**: 既存単灯に同じ分解を当てると傾きは bottom 6.00%・left 7.17% に対し **right だけ 1.84%** ＝基板ごとに違う．カメラ由来（LSC 残差）なら 4 条件に同一の傾きが乗るはずなので，**傾きは bottom/left のモデル化誤差**と確定．実機単灯の振幅も右>左>下で 10% ばらつくのにシミュは 0.5% 以内でフラット（未モデル化）だが，一律ゲイン差では right だけ整合する事実を説明できないため位置・向き・配光の疑いが濃い．**未実施の前提**: ノイズフロア（seed 違い 2 本）を測っていないため上位 4 条件（score 0.079〜0.087）は横並びで「ρ=0.82 > 0.78」とは言えず，主張は「ρ ≈ 0.80〜0.82」まで．また d=16 据え置きにより打ち切り損失が ρ とともに増える（採用条件で 5.1%，ρ=0.90 で 20.6%）点は高 ρ 側との交絡として残る．成果物: `output/bsdf_grid_fix/`（60 条件 + summary.csv），図 `output/_slide_panels/bsdf_grid_{response,contact_sheet,rho_series,top_residual,top_shape}.png`，`scripts/sweep_bsdf_grid_fix.py`・`build_bsdf_grid_panel.py`・`build_bsdf_grid_visual.py`，記録 [TEST_14](06_TEST/TEST_14_材質パラメータ再スイープと残差の構造分解.md)． diff --git a/scripts/build_bsdf_grid_panel.py b/scripts/build_bsdf_grid_panel.py new file mode 100644 index 0000000..5f00638 --- /dev/null +++ b/scripts/build_bsdf_grid_panel.py @@ -0,0 +1,182 @@ +"""ρ×ω×α 再スイープ（面法線修正後）の結果を可視化する. + +`sweep_bsdf_grid_fix.py` が出す output/bsdf_grid_fix/summary.csv を読み， + + 1. 応答曲面: 各指標を ρ の関数として (ω, α) 系列ごとに描く． + 実機の値を水平線で引き，どの ρ で交差するかを目で追えるようにする． + 2. 形比較パネル: スコア上位 N 条件と実機を並べる（他の TEST パネルと同じ流儀． + mean 正規化 + 独立 1-99% 表示 + 同一平滑化）． + +出力: + output/_slide_panels/bsdf_grid_response.png -- 応答曲面（指標 × ρ） + output/_slide_panels/bsdf_grid_top_shape.png -- 上位条件 + 実機の形比較 + +実行方法: + uv run python scripts/build_bsdf_grid_panel.py [--top 4] +""" + +from __future__ import annotations + +import argparse +import csv +import sys +from pathlib import Path + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +from scipy.ndimage import gaussian_filter, median_filter + +sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src")) + +from sm_tias_lightsim.evaluation.roi_presets import SIM_ROI, extract_roi_xywh # noqa: E402 + +matplotlib.rcParams["font.family"] = ["Meiryo", "Yu Gothic", "DejaVu Sans"] + +ROOT = Path(__file__).resolve().parents[1] +SUMMARY_CSV = ROOT / "output" / "bsdf_grid_fix" / "summary.csv" +OUTPUT_BASE = ROOT / "output" / "bsdf_grid_fix" +REFERENCE_DIR = ROOT / "output" / "full_base_fix" +REAL_FULL = ROOT / "lightsim" / "luminance_maps" / "full_MiniTIAS_QM_20260601_124500.npy" +PANEL_DIR = ROOT / "output" / "_slide_panels" + +# 実機（全灯）の指標: sweep_bsdf_grid_fix.py と同じ出典 +REAL_VALUES = { + "cov": 0.0695, + "center_periphery_ratio": 1.0911, + "gradient_magnitude": 8.35, + "radial_min_max_ratio": 0.860, +} +METRIC_LABELS = { + "cov": "CoV", + "center_periphery_ratio": "中心/周辺比", + "gradient_magnitude": "勾配 [%]", + "radial_min_max_ratio": "動径 min/max", + "ncc_smooth": "NCC 対実機（平滑後）", + "score": "総合スコア（小さいほど実機に近い）", +} + + +def load_rows() -> list[dict[str, float]]: + """summary.csv を数値化して読み込む.""" + if not SUMMARY_CSV.exists(): + raise SystemExit(f"エラー: {SUMMARY_CSV} がありません（スイープ未実行）") + rows = [] + with SUMMARY_CSV.open(encoding="utf-8") as fh: + for raw in csv.DictReader(fh): + row: dict[str, float] = {"job": raw["job"], "stage": raw["stage"]} + for key, value in raw.items(): + if key in ("job", "stage"): + continue + row[key] = float(value) if value not in ("", None) else float("nan") + rows.append(row) + return rows + + +def _smooth(arr: np.ndarray, is_sim: bool) -> np.ndarray: + """他パネルと同一の平滑化（median5 → gauss σ4）と ROI 抽出.""" + if arr.ndim == 3: + arr = 0.2126 * arr[..., 0] + 0.7152 * arr[..., 1] + 0.0722 * arr[..., 2] + arr = arr.astype(np.float64) + if is_sim: + arr = extract_roi_xywh(arr, SIM_ROI) + return gaussian_filter(median_filter(arr, size=5), sigma=4) + + +def _job_luminance(job: str) -> Path: + """ジョブ名から luminance.npy のパスを返す（基準点は full_base_fix を指す）.""" + candidate = OUTPUT_BASE / job / "luminance.npy" + if candidate.exists(): + return candidate + return REFERENCE_DIR / "luminance.npy" + + +def plot_response(rows: list[dict[str, float]]) -> Path: + """指標を ρ の関数として (ω, α) 系列ごとに描き，実機値を水平線で示す.""" + metrics = ["cov", "gradient_magnitude", "center_periphery_ratio", + "radial_min_max_ratio", "ncc_smooth", "score"] + series: dict[tuple[float, float], list[dict[str, float]]] = {} + for row in rows: + series.setdefault((row["omega"], row["alpha"]), []).append(row) + + fig, axes = plt.subplots(2, 3, figsize=(17, 9), constrained_layout=True) + for ax, metric in zip(axes.ravel(), metrics, strict=True): + for (omega, alpha), group in sorted(series.items()): + group = sorted(group, key=lambda r: r["rho"]) + label = f"ω={omega:g}" if omega == 0 else f"ω={omega:g}, α={alpha:g}" + ax.plot([r["rho"] for r in group], [r[metric] for r in group], + marker="o", markersize=4, linewidth=1.2, label=label) + if metric in REAL_VALUES: + ax.axhline(REAL_VALUES[metric], color="k", linestyle="--", linewidth=1.5, + label="実機" if metric == "cov" else None) + ax.set_xlabel("ρ（拡散反射率）") + ax.set_ylabel(METRIC_LABELS[metric]) + ax.set_title(METRIC_LABELS[metric]) + ax.grid(alpha=0.3) + axes.ravel()[0].legend(fontsize=7, ncol=2) + fig.suptitle("面法線修正後の ρ×ω×α 応答曲面 — 破線＝実機全灯" + "（spp=1024, max_depth=16 固定）", fontsize=13) + PANEL_DIR.mkdir(parents=True, exist_ok=True) + path = PANEL_DIR / "bsdf_grid_response.png" + fig.savefig(path, dpi=150, bbox_inches="tight") + plt.close(fig) + return path + + +def plot_top_shapes(rows: list[dict[str, float]], top: int) -> Path: + """スコア上位 top 条件と実機の輝度マップを並べる.""" + best = sorted(rows, key=lambda r: r["score"])[:top] + real = _smooth(np.load(REAL_FULL), False) + cells = [(f"{r['job']}\nscore={r['score']:.4f} CoV={r['cov']:.4f} " + f"NCC={r['ncc_smooth']:.3f}", _smooth(np.load(_job_luminance(r["job"])), True)) + for r in best] + real_label = (f"実機 全灯\nCoV={REAL_VALUES['cov']:.4f} " + f"勾配={REAL_VALUES['gradient_magnitude']}%") + cells.append((real_label, real)) + + fig, axes = plt.subplots(1, len(cells), figsize=(4.2 * len(cells), 4.6), + constrained_layout=True) + for ax, (label, img) in zip(np.atleast_1d(axes), cells, strict=True): + disp = img / img.mean() + vmin, vmax = np.percentile(disp, [1, 99]) + ax.imshow(disp, cmap="hot", vmin=vmin, vmax=vmax) + ax.set_title(label, fontsize=10) + ax.axis("off") + fig.suptitle("スコア上位条件 vs 実機 — mean 正規化・各セル独立 1–99% 表示（形の比較専用）", + fontsize=13) + PANEL_DIR.mkdir(parents=True, exist_ok=True) + path = PANEL_DIR / "bsdf_grid_top_shape.png" + fig.savefig(path, dpi=150, bbox_inches="tight") + plt.close(fig) + return path + + +def main() -> int: + """応答曲面と上位条件パネルを生成するメインエントリポイント.""" + parser = argparse.ArgumentParser(description="ρ×ω×α スイープ結果の可視化") + parser.add_argument("--top", type=int, default=4, help="形比較に並べる上位条件数") + args = parser.parse_args() + + rows = load_rows() + print(f"読み込み: {len(rows)} 条件") + print(f"応答曲面: {plot_response(rows)}") + print(f"形比較: {plot_top_shapes(rows, args.top)}") + + print("\n=== スコア上位 10 ===") + header = (f"{'job':22s} {'score':>7s} {'CoV':>7s} {'C/P':>7s} " + f"{'勾配%':>7s} {'動径':>6s} {'NCC':>6s} {'打切%':>6s}") + print(header) + for row in sorted(rows, key=lambda r: r["score"])[:10]: + print(f"{row['job']:22s} {row['score']:7.4f} {row['cov']:7.4f} " + f"{row['center_periphery_ratio']:7.4f} {row['gradient_magnitude']:7.2f} " + f"{row['radial_min_max_ratio']:6.3f} {row['ncc_smooth']:6.3f} " + f"{row['truncation_loss_pct']:6.1f}") + print(f"{'実機':22s} {'-':>7s} {REAL_VALUES['cov']:7.4f} " + f"{REAL_VALUES['center_periphery_ratio']:7.4f} " + f"{REAL_VALUES['gradient_magnitude']:7.2f} " + f"{REAL_VALUES['radial_min_max_ratio']:6.3f} {'1.000':>6s} {'-':>6s}") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/scripts/build_bsdf_grid_visual.py b/scripts/build_bsdf_grid_visual.py new file mode 100644 index 0000000..45a2252 --- /dev/null +++ b/scripts/build_bsdf_grid_visual.py @@ -0,0 +1,225 @@ +"""ρ×ω×α スイープの輝度マップを「目で判断する」ための可視化. + +指標（CoV / NCC 等）は分布を集約した量であり，形そのものの判断材料にはならない +（TEST_03 の主題）．本スクリプトは summary.csv の順位ではなく **輝度マップそのもの** +を並べ，どのパラメータで形がどう変わるか／変わらないかを目視できるようにする． + +生成物: + 1. bsdf_grid_contact_sheet.png -- 全 60 条件 + 実機のコンタクトシート + （行 = ρ，列 = (ω, α)．全セル共通カラースケールなので差がそのまま見える） + 2. bsdf_grid_rho_series.png -- ω=0 純 diffuse の ρ 系列 + 実機． + 上段 = 輝度マップ，下段 = 実機との相対残差マップ [%] + 3. bsdf_grid_top_residual.png -- スコア上位条件の残差マップ（発散カラーマップ） + +残差の定義: 平均で正規化したうえで (sim - real) / mean(real) [%]． +正（赤）= シミュが明るすぎる領域，負（青）= 暗すぎる領域． + +実行方法: + uv run python scripts/build_bsdf_grid_visual.py [--top 4] +""" + +from __future__ import annotations + +import argparse +import csv +import sys +from pathlib import Path + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +from scipy.ndimage import gaussian_filter, median_filter +from skimage.transform import resize + +sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src")) + +from sm_tias_lightsim.evaluation.roi_presets import SIM_ROI, extract_roi_xywh # noqa: E402 + +matplotlib.rcParams["font.family"] = ["Meiryo", "Yu Gothic", "DejaVu Sans"] + +ROOT = Path(__file__).resolve().parents[1] +SUMMARY_CSV = ROOT / "output" / "bsdf_grid_fix" / "summary.csv" +OUTPUT_BASE = ROOT / "output" / "bsdf_grid_fix" +REFERENCE_DIR = ROOT / "output" / "full_base_fix" +REAL_FULL = ROOT / "lightsim" / "luminance_maps" / "full_MiniTIAS_QM_20260601_124500.npy" +PANEL_DIR = ROOT / "output" / "_slide_panels" + +RESIDUAL_LIMIT = 12.0 # 残差カラーバーの上下限 [%] + + +def load_rows() -> list[dict]: + """summary.csv を数値化して読み込む.""" + if not SUMMARY_CSV.exists(): + raise SystemExit(f"エラー: {SUMMARY_CSV} がありません（スイープ未実行）") + rows = [] + with SUMMARY_CSV.open(encoding="utf-8") as fh: + for raw in csv.DictReader(fh): + row: dict = {"job": raw["job"], "stage": raw["stage"]} + for key, value in raw.items(): + if key not in ("job", "stage"): + row[key] = float(value) if value not in ("", None) else float("nan") + rows.append(row) + return rows + + +def _smooth(arr: np.ndarray, is_sim: bool) -> np.ndarray: + """他パネルと同一の平滑化（median5 → gauss σ4）と ROI 抽出.""" + if arr.ndim == 3: + arr = 0.2126 * arr[..., 0] + 0.7152 * arr[..., 1] + 0.0722 * arr[..., 2] + arr = arr.astype(np.float64) + if is_sim: + arr = extract_roi_xywh(arr, SIM_ROI) + return gaussian_filter(median_filter(arr, size=5), sigma=4) + + +def _job_map(job: str) -> np.ndarray: + """ジョブ名から平滑化・mean 正規化済みの輝度マップを返す.""" + path = OUTPUT_BASE / job / "luminance.npy" + if not path.exists(): + path = REFERENCE_DIR / "luminance.npy" + img = _smooth(np.load(path), True) + return img / img.mean() + + +def _real_map(shape: tuple[int, int]) -> np.ndarray: + """実機マップを指定 shape に合わせ mean 正規化して返す.""" + real = _smooth(np.load(REAL_FULL), False) + real = resize(real, shape, anti_aliasing=True) + return real / real.mean() + + +def _residual(sim_norm: np.ndarray, real_norm: np.ndarray) -> np.ndarray: + """mean 正規化済みどうしの相対残差 [%]（正＝シミュが明るすぎ）.""" + return 100.0 * (sim_norm - real_norm) + + +def contact_sheet(rows: list[dict]) -> Path: + """全条件を 行=ρ × 列=(ω, α) で並べた一覧（共通カラースケール）.""" + rho_values = sorted({r["rho"] for r in rows}) + combos = sorted({(r["omega"], r["alpha"]) for r in rows}) + lookup = {(r["rho"], r["omega"], r["alpha"]): r for r in rows} + + ncols = len(combos) + 1 # 右端に実機列 + fig, axes = plt.subplots(len(rho_values), ncols, + figsize=(1.75 * ncols, 1.95 * len(rho_values)), + constrained_layout=True) + sample = _job_map(rows[0]["job"]) + real = _real_map(sample.shape) + vmin, vmax = np.percentile(real, [1, 99]) + + for row_idx, rho in enumerate(rho_values): + for col_idx, (omega, alpha) in enumerate(combos): + ax = axes[row_idx][col_idx] + entry = lookup.get((rho, omega, alpha)) + if entry is None: + ax.axis("off") + continue + ax.imshow(_job_map(entry["job"]), cmap="hot", vmin=vmin, vmax=vmax) + ax.set_title(f"NCC={entry['ncc_smooth']:.3f}", fontsize=6.5, pad=2) + ax.set_xticks([]) + ax.set_yticks([]) + if col_idx == 0: + ax.set_ylabel(f"ρ={rho:g}", fontsize=9) + if row_idx == 0: + label = "ω=0" if omega == 0 else f"ω={omega:g}\nα={alpha:g}" + ax.text(0.5, 1.35, label, transform=ax.transAxes, + ha="center", va="bottom", fontsize=8) + ax_real = axes[row_idx][-1] + im = ax_real.imshow(real, cmap="hot", vmin=vmin, vmax=vmax) + ax_real.set_xticks([]) + ax_real.set_yticks([]) + ax_real.set_title("実機", fontsize=7, pad=2) + if row_idx == 0: + ax_real.text(0.5, 1.35, "参照", transform=ax_real.transAxes, + ha="center", va="bottom", fontsize=8) + + fig.colorbar(im, ax=axes, shrink=0.5, label="mean 正規化輝度（全セル共通スケール）") + fig.suptitle("全 60 条件の輝度マップ一覧 — 行=ρ，列=(ω, α)，右端=実機" + "（全セル共通カラースケール・mean 正規化）", fontsize=13) + PANEL_DIR.mkdir(parents=True, exist_ok=True) + path = PANEL_DIR / "bsdf_grid_contact_sheet.png" + fig.savefig(path, dpi=140, bbox_inches="tight") + plt.close(fig) + return path + + +def rho_series(rows: list[dict]) -> Path: + """ω=0 の ρ 系列を輝度マップ（上段）と残差マップ（下段）で並べる.""" + series = sorted([r for r in rows if r["omega"] == 0], key=lambda r: r["rho"]) + sample = _job_map(series[0]["job"]) + real = _real_map(sample.shape) + vmin, vmax = np.percentile(real, [1, 99]) + + ncols = len(series) + 1 + fig, axes = plt.subplots(2, ncols, figsize=(2.55 * ncols, 5.9), + constrained_layout=True) + for col, entry in enumerate(series): + sim = _job_map(entry["job"]) + axes[0][col].imshow(sim, cmap="hot", vmin=vmin, vmax=vmax) + axes[0][col].set_title(f"ρ={entry['rho']:g}\nCoV={entry['cov']:.4f} " + f"NCC={entry['ncc_smooth']:.3f}", fontsize=9) + res = axes[1][col].imshow(_residual(sim, real), cmap="RdBu_r", + vmin=-RESIDUAL_LIMIT, vmax=RESIDUAL_LIMIT) + rms = float(np.sqrt(np.mean(_residual(sim, real) ** 2))) + axes[1][col].set_title(f"残差 RMS={rms:.2f}%", fontsize=9) + top = axes[0][-1].imshow(real, cmap="hot", vmin=vmin, vmax=vmax) + axes[0][-1].set_title("実機 全灯\nCoV=0.0695", fontsize=9) + axes[1][-1].axis("off") + for ax in axes.ravel(): + if ax.axison: + ax.set_xticks([]) + ax.set_yticks([]) + + fig.colorbar(top, ax=axes[0], shrink=0.8, label="mean 正規化輝度") + fig.colorbar(res, ax=axes[1], shrink=0.8, label="相対残差 [%]（赤=シミュが明るい）") + fig.suptitle("純 diffuse（ω=0）の ρ 系列 — 上段: 輝度マップ / 下段: 実機との残差", + fontsize=13) + PANEL_DIR.mkdir(parents=True, exist_ok=True) + path = PANEL_DIR / "bsdf_grid_rho_series.png" + fig.savefig(path, dpi=150, bbox_inches="tight") + plt.close(fig) + return path + + +def top_residual(rows: list[dict], top: int) -> Path: + """スコア上位条件の残差マップを並べる（形の判断はここが主戦場）.""" + best = sorted(rows, key=lambda r: r["score"])[:top] + sample = _job_map(best[0]["job"]) + real = _real_map(sample.shape) + + fig, axes = plt.subplots(1, top, figsize=(3.7 * top, 4.5), constrained_layout=True) + for ax, entry in zip(np.atleast_1d(axes), best, strict=True): + residual = _residual(_job_map(entry["job"]), real) + im = ax.imshow(residual, cmap="RdBu_r", vmin=-RESIDUAL_LIMIT, vmax=RESIDUAL_LIMIT) + rms = float(np.sqrt(np.mean(residual ** 2))) + p95 = float(np.percentile(np.abs(residual), 95)) + ax.set_title(f"{entry['job']}\nRMS={rms:.2f}% P95={p95:.2f}% " + f"NCC={entry['ncc_smooth']:.3f}", fontsize=10) + ax.set_xticks([]) + ax.set_yticks([]) + fig.colorbar(im, ax=axes, shrink=0.8, label="相対残差 [%]（赤=シミュが明るい / 青=暗い）") + fig.suptitle("スコア上位条件の残差マップ — 赤青の構造が「材質では埋まらない形の差」", + fontsize=13) + PANEL_DIR.mkdir(parents=True, exist_ok=True) + path = PANEL_DIR / "bsdf_grid_top_residual.png" + fig.savefig(path, dpi=150, bbox_inches="tight") + plt.close(fig) + return path + + +def main() -> int: + """3 種の可視化を生成するメインエントリポイント.""" + parser = argparse.ArgumentParser(description="ρ×ω×α スイープの輝度マップ可視化") + parser.add_argument("--top", type=int, default=4, help="残差マップに並べる上位条件数") + args = parser.parse_args() + + rows = load_rows() + print(f"読み込み: {len(rows)} 条件") + print(f"一覧: {contact_sheet(rows)}") + print(f"ρ 系列: {rho_series(rows)}") + print(f"残差: {top_residual(rows, args.top)}") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/scripts/sweep_bsdf_grid_fix.py b/scripts/sweep_bsdf_grid_fix.py new file mode 100644 index 0000000..10d2049 --- /dev/null +++ b/scripts/sweep_bsdf_grid_fix.py @@ -0,0 +1,344 @@ +"""面法線修正後の ρ × ω × α 再スイープ（TEST_12 後の条件引き直し）. + +面法線バグ（TEST_12）修正前に実施した `sweep_reflectance.py` / `sweep_blend_bsdf.py` +の結論は，評価板のシェーディング法線アーティファクトが不均一性の 6 割を占めていた +状態で得たものであり無効．本スクリプトは同じ 3 軸を **修正後の条件で引き直し**， +実機（全灯）に最も近い組合せを探す． + +スイープ軸: + ρ : diffuse reflectance（内壁 printer_white の拡散反射率） + ω : blendbsdf weight（roughconductor 成分の混合比．0 なら純 diffuse） + α : roughconductor の GGX 粗さ + +固定条件（既存レンダと完全に同一に保つ）: + spp=1024 / max_depth=16 / 512x512 / llvm_ad_rgb / path / orthographic + → configs/_tmp_fn_full_base.yaml をテンプレートとし materials のみ書き換える + +探索戦略（無人実行を想定）: + ρ×ω×α の全数グリッドを，粗い ρ（COARSE_RHO）→細かい ρ（FINE_RHO）の 2 パスで流す． + 時間予算で打ち切られても ρ 軸が粗いなりに全域カバーされている状態を保つための順序． + 1 ジョブ完了ごとに summary.csv を書き直し，既にレンダ済みのジョブは再利用する + （基準点 ρ=0.85/ω=0.05/α=0.1 は output/full_base_fix/ を再利用し再レンダしない）. + +出力（1 ジョブ完了ごとに逐次追記するので途中終了しても結果は残る）: + output/bsdf_grid_fix//luminance.npy, metrics.csv + output/bsdf_grid_fix/summary.csv + configs/_tmp_bsdf_grid/scene_.yaml + +実行方法: + uv run python scripts/sweep_bsdf_grid_fix.py [--budget-hours 5.5] +""" + +from __future__ import annotations + +import argparse +import copy +import csv +import subprocess +import sys +import time +from pathlib import Path + +import numpy as np +import yaml +from scipy.ndimage import gaussian_filter, median_filter +from skimage.transform import resize + +sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src")) + +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 # noqa: E402 + +ROOT = Path(__file__).resolve().parents[1] + +TEMPLATE_YAML = ROOT / "configs" / "_tmp_fn_full_base.yaml" +EVAL_CONFIG_YAML = ROOT / "configs" / "eval_baseline_dng_target.yaml" +TMP_DIR = ROOT / "configs" / "_tmp_bsdf_grid" +OUTPUT_BASE = ROOT / "output" / "bsdf_grid_fix" +SUMMARY_CSV = OUTPUT_BASE / "summary.csv" + +REAL_FULL = ROOT / "lightsim" / "luminance_maps" / "full_MiniTIAS_QM_20260601_124500.npy" +REFERENCE_JOB_DIR = ROOT / "output" / "full_base_fix" + +# 基準点（= output/full_base_fix/ を生成した条件） +REF_RHO, REF_OMEGA, REF_ALPHA = 0.85, 0.05, 0.1 + +# 全数グリッド．1 ジョブ約 200 秒（実測）なので 5.5h でおよそ 90 ジョブ入る． +# 粗パス（COARSE_RHO）を先に流し切ってから細パス（FINE_RHO）で ρ を埋める順序にし， +# 時間予算で打ち切られても ρ 軸が粗いなりに全域カバーされている状態を保つ． +COARSE_RHO = [0.75, 0.85, 0.90] +FINE_RHO = [0.78, 0.82, 0.88] +OMEGA_VALUES = [0.0, 0.05, 0.15, 0.30] +ALPHA_VALUES = [0.1, 0.3, 0.6] + +# 実機（全灯）の指標: configs/eval_baseline_dng_target.yaml ヘッダの記録値 +REAL_COV = 0.0695 +REAL_CP = 1.0911 +REAL_GRADIENT = 8.35 +REAL_RADIAL = 0.860 + +SUMMARY_COLUMNS = [ + "job", "rho", "omega", "alpha", "stage", "elapsed_s", + "cov", "center_periphery_ratio", "gradient_magnitude", "radial_min_max_ratio", + "ncc_smooth", "peak_mean_smooth", "top_bot_smooth", + "e_cov", "e_cp", "e_grad", "e_radial", "e_ncc", "score", + "truncation_loss_pct", +] + + +def job_name(rho: float, omega: float, alpha: float) -> str: + """ρ/ω/α からジョブ名（ディレクトリ名）を作る.""" + return f"r{rho:.2f}_w{omega:.2f}_a{alpha:.1f}" + + +def truncation_loss_pct(rho: float, omega: float) -> float: + """max_depth=16 で取りこぼす相対エネルギー [%] を実効反射率から見積もる. + + blendbsdf の roughconductor 成分は material="none"（反射率 1.0 の完全反射体）なので + 実効的な壁の反射率は rho*(1-omega) + omega となり，rho 単体より高くなる． + TEST_13 の check_truncation_bias と同じ rho_eff^(max_depth-1) を用いる． + """ + rho_eff = rho * (1.0 - omega) + omega + return round(100.0 * rho_eff ** 15, 2) + + +def build_scene_yaml(rho: float, omega: float, alpha: float, dest: Path) -> None: + """テンプレートの printer_white だけを差し替えたシーン YAML を書き出す. + + ω=0 のときは blendbsdf を使わず純 diffuse にする（weight=0 の縮退を避ける）. + 他のセクション（geometry / lights / sensor / render）は一切変更しない. + """ + config = yaml.safe_load(TEMPLATE_YAML.read_text(encoding="utf-8")) + config = copy.deepcopy(config) + if omega <= 0.0: + material = {"type": "diffuse", "reflectance": [rho, rho, rho]} + else: + material = { + "type": "blendbsdf", + "weight": omega, + "components": [ + {"type": "diffuse", "reflectance": [rho, rho, rho]}, + {"type": "roughconductor", "alpha": alpha, "distribution": "ggx"}, + ], + } + config["materials"]["printer_white"] = material + config["name"] = f"bsdf_grid_{job_name(rho, omega, alpha)}" + # メッシュの相対パスは YAML 自身のディレクトリ基準で解決されるため， + # テンプレート（configs/ 直下）から別ディレクトリへ移す時点で絶対パス化する． + for mesh in config["geometry"]["meshes"]: + mesh["path"] = str((TEMPLATE_YAML.parent / mesh["path"]).resolve()) + dest.parent.mkdir(parents=True, exist_ok=True) + dest.write_text(yaml.safe_dump(config, allow_unicode=True, sort_keys=False), encoding="utf-8") + + +def _smooth(arr: np.ndarray, is_sim: bool) -> np.ndarray: + """build_full_vs_singles_panel.py と同一の平滑化（median5 → gauss σ4）.""" + if arr.ndim == 3: + arr = 0.2126 * arr[..., 0] + 0.7152 * arr[..., 1] + 0.0722 * arr[..., 2] + arr = arr.astype(np.float64) + if is_sim: + arr = extract_roi_xywh(arr, SIM_ROI) + return gaussian_filter(median_filter(arr, size=5), sigma=4) + + +def _shape_stats(sim_path: Path, real_smooth: np.ndarray) -> dict[str, float]: + """平滑化後の NCC 対実機・peak/mean・top/bot を返す（MC ノイズ支配を避けるため平滑後）.""" + sim = _smooth(np.load(sim_path), True) + resized = resize(sim, real_smooth.shape, anti_aliasing=True) + half = sim.shape[0] // 2 + return { + "ncc_smooth": float(calc_ncc(resized, real_smooth)), + "peak_mean_smooth": float(sim.max() / sim.mean()), + "top_bot_smooth": float(sim[:half].mean() / sim[half:].mean()), + } + + +def _score(metrics: dict[str, float]) -> dict[str, float]: + """実機との相対乖離から総合スコアを組む（小さいほど実機に近い）. + + 形（NCC）と振幅系 4 指標（CoV / C-P / 勾配 / 動径 min-max）を半々で重み付けする． + C/P と動径 min/max は「1 からの超過分」で相対化する（1 近傍で比が鈍るため）. + """ + e_cov = abs(metrics["cov"] - REAL_COV) / REAL_COV + e_cp = abs((metrics["center_periphery_ratio"] - 1.0) - (REAL_CP - 1.0)) / (REAL_CP - 1.0) + e_grad = abs(metrics["gradient_magnitude"] - REAL_GRADIENT) / REAL_GRADIENT + e_radial = (abs((1.0 - metrics["radial_min_max_ratio"]) - (1.0 - REAL_RADIAL)) + / (1.0 - REAL_RADIAL)) + e_ncc = 1.0 - metrics["ncc_smooth"] + amplitude = (e_cov + e_cp + e_grad + e_radial) / 4.0 + return { + "e_cov": e_cov, "e_cp": e_cp, "e_grad": e_grad, "e_radial": e_radial, + "e_ncc": e_ncc, "score": 0.5 * amplitude + 0.5 * e_ncc, + } + + +def read_metrics_csv(path: Path) -> dict[str, float]: + """sm-evaluate が出す metrics.csv の 1 行目を辞書で返す.""" + with path.open(encoding="utf-8") as fh: + row = next(csv.DictReader(fh)) + out = {} + for key in ("cov", "center_periphery_ratio", "gradient_magnitude", "radial_min_max_ratio"): + out[key] = float(row[key]) + return out + + +def run_job(rho: float, omega: float, alpha: float, stage: str, + real_smooth: np.ndarray, log) -> dict[str, float] | None: + """1 条件をレンダ→評価→指標算出まで実行して集約行を返す（失敗時 None）. + + 既に metrics.csv と luminance.npy が揃っているジョブは再レンダせず再利用する + （中断・再開や，グリッドを広げて流し直すときに無駄な計算をしないため）. + """ + name = job_name(rho, omega, alpha) + out_dir = OUTPUT_BASE / name + scene_yaml = TMP_DIR / f"scene_{name}.yaml" + if (out_dir / "metrics.csv").exists() and (out_dir / "luminance.npy").exists(): + row = {"job": name, "rho": rho, "omega": omega, "alpha": alpha, + "stage": stage + "(再利用)", "elapsed_s": 0.0} + row.update(read_metrics_csv(out_dir / "metrics.csv")) + row.update(_shape_stats(out_dir / "luminance.npy", real_smooth)) + row.update(_score(row)) + row["truncation_loss_pct"] = truncation_loss_pct(rho, omega) + log(f"--- {stage}: {name} → 既存結果を再利用 (score={row['score']:.4f})") + return row + log(f"--- {stage}: {name} (rho={rho}, omega={omega}, alpha={alpha}) ---") + build_scene_yaml(rho, omega, alpha, scene_yaml) + out_dir.mkdir(parents=True, exist_ok=True) + start = time.time() + for cmd in ( + ["uv", "run", "sm-render", "--scene", str(scene_yaml), "--output", str(out_dir)], + ["uv", "run", "sm-evaluate", "--image", str(out_dir / "luminance.npy"), + "--config", str(EVAL_CONFIG_YAML), "--output", str(out_dir)], + ): + result = subprocess.run(cmd, cwd=ROOT, capture_output=True, text=True) + if result.returncode != 0: + log(f" !! 失敗 ({cmd[2]}): {result.stderr[-600:]}") + return None + elapsed = time.time() - start + row = {"job": name, "rho": rho, "omega": omega, "alpha": alpha, + "stage": stage, "elapsed_s": round(elapsed, 1)} + row.update(read_metrics_csv(out_dir / "metrics.csv")) + row.update(_shape_stats(out_dir / "luminance.npy", real_smooth)) + row.update(_score(row)) + row["truncation_loss_pct"] = truncation_loss_pct(rho, omega) + log(f" -> {elapsed:.0f}s CoV={row['cov']:.4f} (実機 {REAL_COV}) " + f"NCC={row['ncc_smooth']:.4f} 勾配={row['gradient_magnitude']:.2f}% " + f"score={row['score']:.4f}") + return row + + +def append_summary(rows: list[dict[str, float]]) -> None: + """集約 CSV を毎回まるごと書き直す（途中終了でも最新状態が残る）.""" + OUTPUT_BASE.mkdir(parents=True, exist_ok=True) + ordered = sorted(rows, key=lambda r: r["score"]) + with SUMMARY_CSV.open("w", newline="", encoding="utf-8") as fh: + writer = csv.DictWriter(fh, fieldnames=SUMMARY_COLUMNS) + writer.writeheader() + for row in ordered: + writer.writerow({k: row.get(k, "") for k in SUMMARY_COLUMNS}) + + +def reference_row(real_smooth: np.ndarray, log) -> dict[str, float] | None: + """基準点は既存の output/full_base_fix/ を再利用する（再レンダしない）.""" + metrics_csv = REFERENCE_JOB_DIR / "metrics.csv" + luminance = REFERENCE_JOB_DIR / "luminance.npy" + if not (metrics_csv.exists() and luminance.exists()): + log(" 基準点 output/full_base_fix/ が無いため再レンダします") + return run_job(REF_RHO, REF_OMEGA, REF_ALPHA, "stage1", real_smooth, log) + row = {"job": job_name(REF_RHO, REF_OMEGA, REF_ALPHA), "rho": REF_RHO, + "omega": REF_OMEGA, "alpha": REF_ALPHA, "stage": "reference", "elapsed_s": 0.0} + row.update(read_metrics_csv(metrics_csv)) + row.update(_shape_stats(luminance, real_smooth)) + row.update(_score(row)) + row["truncation_loss_pct"] = truncation_loss_pct(REF_RHO, REF_OMEGA) + log(f" 基準点（既存 full_base_fix 再利用）: CoV={row['cov']:.4f} " + f"NCC={row['ncc_smooth']:.4f} score={row['score']:.4f}") + return row + + +def build_grid(rho_values: list[float]) -> list[tuple[float, float, float]]: + """指定した ρ 値集合について ω×α の全組合せを返す. + + ω=0（純 diffuse）は α に依存しないので α は 1 値だけ回す. + """ + combos = [] + for rho in rho_values: + for omega in OMEGA_VALUES: + for alpha in (ALPHA_VALUES[:1] if omega <= 0.0 else ALPHA_VALUES): + combos.append((rho, omega, alpha)) + return combos + + +def main() -> int: + """段階スイープを時間予算内で実行するメインエントリポイント.""" + parser = argparse.ArgumentParser(description="面法線修正後の ρ×ω×α 再スイープ") + parser.add_argument("--budget-hours", type=float, default=5.5, + help="全体の時間予算（時間）．超過が見込まれる時点で打ち切る") + args = parser.parse_args() + + deadline = time.time() + args.budget_hours * 3600.0 + log_path = ROOT / "output" / "_sweep_bsdf_grid_fix.log" + log_path.parent.mkdir(parents=True, exist_ok=True) + log_file = log_path.open("w", encoding="utf-8") + + def log(message: str) -> None: + print(message, flush=True) + log_file.write(message + "\n") + log_file.flush() + + if not REAL_FULL.exists(): + log(f"エラー: 実機マップ {REAL_FULL} が見つかりません") + return 1 + real_smooth = _smooth(np.load(REAL_FULL), False) + + log(f"=== ρ×ω×α 再スイープ（面法線修正後）予算 {args.budget_hours}h ===") + log(f" 固定条件: spp=1024 / max_depth=16 / 512^2（{TEMPLATE_YAML.name} 準拠）") + log(f" 基準点: rho={REF_RHO}, omega={REF_OMEGA}, alpha={REF_ALPHA}") + + rows: list[dict[str, float]] = [] + ref = reference_row(real_smooth, log) + if ref is not None: + rows.append(ref) + append_summary(rows) + + durations: list[float] = [] + + def budget_ok() -> bool: + """残り時間が想定 1 ジョブ分を確保できるか（実測の最大所要で判定）.""" + expected = max(durations) if durations else 900.0 + return time.time() + expected < deadline + + done = {(r["rho"], r["omega"], r["alpha"]) for r in rows} + passes = (("粗パス", build_grid(COARSE_RHO)), ("細パス", build_grid(FINE_RHO))) + for pass_name, combos in passes: + pending = [c for c in combos if c not in done] + log(f"\n=== {pass_name}: {len(pending)} ジョブ（ρ×ω×α 全組合せ） ===") + for rho, omega, alpha in pending: + if not budget_ok(): + log(f" 時間予算のため {pass_name} を打ち切ります") + break + row = run_job(rho, omega, alpha, pass_name, real_smooth, log) + if row is not None: + rows.append(row) + done.add((rho, omega, alpha)) + if row["elapsed_s"] > 0: + durations.append(row["elapsed_s"]) + append_summary(rows) + else: + continue + break + + log(f"\n=== 完了: {len(rows)} 条件 ===") + log(f" 実機: CoV={REAL_COV:.4f}, C/P={REAL_CP:.4f}, " + f"勾配={REAL_GRADIENT:.2f}%, 動径 min/max={REAL_RADIAL:.3f}") + log(" スコア上位 5:") + for row in sorted(rows, key=lambda r: r["score"])[:5]: + log(f" {row['job']:22s} score={row['score']:.4f} CoV={row['cov']:.4f} " + f"NCC={row['ncc_smooth']:.4f} 勾配={row['gradient_magnitude']:.2f}%") + log(f" 集約 CSV: {SUMMARY_CSV}") + log_file.close() + return 0 + + +if __name__ == "__main__": + sys.exit(main())