Newer
Older
SmTIAS-LightSim / tests / render / test_physics_validation.py
"""解析解に対する物理量の妥当性検証テスト.

仕様出典:
- docs/06_TEST/TEST_13_物理モデル健全性の一括監査.md

積分球(内面反射率 rho の閉じた球殻)は放射伝達に厳密な閉形式解を持つ.
全放射束 Phi の光源を内部に置くと,壁面の全照度は

    E_total = Phi / (A (1 - rho))          A = 4 pi R^2

となり,直接光が届かない領域の壁面照度は間接成分のみの

    E_indirect = rho * E_total

なので,そこを見たカメラが観測する放射輝度は

    L = rho * E_indirect / pi = rho^2 * Phi / (pi * A * (1 - rho))

これを実測と突き合わせると,**単位系(mm→m)・反射率・発光面の放射束・
センサ校正・多重反射の収束**をまとめて 1 本で検証できる.いずれかが壊れていれば
一致しない.TEST_12 のような「エラーも警告も出ないまま結果だけが歪む」故障を
数値で捕まえるための土台.

本モジュールを `slow` にしない理由: レンダリングを伴うが 32x32 / spp 128 の
極小シーンで合計数秒しかかからない.物理的健全性の検査は通常実行に含めないと
意味がない(実際 `slow` 層はレンダ経路の破損で長期間全滅していたが気づかれなかった).
"""

from __future__ import annotations

import math
from pathlib import Path

import numpy as np
import pytest
import trimesh
import yaml

from sm_tias_lightsim.render.runner import render_scene
from sm_tias_lightsim.scene.schema import load_scene_yaml

# --- 積分球の構成(すべて mm) ---
_SPHERE_RADIUS_MM = 100.0
_REFLECTANCE = 0.8
_EMITTER_RADIANCE = 1000.0  # W/(m^2 sr)
_EMITTER_SIZE_MM = 4.0
# 発光面は -Z を向き z=-50mm に置く.z > -50mm の領域には直接光が一切届かないので,
# +Z 極を見るカメラは純粋に間接成分だけを観測する.
_EMITTER_CENTER_MM = (0.0, 0.0, -50.0)

# 反射率 0.8 では rho^47 ~ 3e-5 なので打ち切り誤差は無視できる
_CONVERGED_MAX_DEPTH = 48
_SPP = 128
_FILM = 32

# 許容誤差 2%: 内訳は icosphere(内接多面体)の面積欠損 ~0.2% と MC ノイズ ~0.3%.
# 実測は解析解比 0.9996〜1.0001 だった(TEST_13).
_TOLERANCE = 0.02

_MM_TO_M = 0.001


def _analytic_shadowed_radiance(radiance: float = _EMITTER_RADIANCE) -> float:
    """直接光が届かない壁面の放射輝度の解析解 [W/(m^2 sr)].

    Args:
        radiance: 発光面の放射輝度 [W/(m^2 sr)].

    Returns:
        壁面の放射輝度 [W/(m^2 sr)].
    """
    area = 4.0 * math.pi * (_SPHERE_RADIUS_MM * _MM_TO_M) ** 2
    emitter_area = (_EMITTER_SIZE_MM * _MM_TO_M) ** 2
    # ランバート発光面の全放射束 Phi = pi * L * A
    flux = math.pi * radiance * emitter_area
    irradiance_total = flux / (area * (1.0 - _REFLECTANCE))
    irradiance_indirect = _REFLECTANCE * irradiance_total
    return _REFLECTANCE * irradiance_indirect / math.pi


def _write_sphere_stl(path: Path, inward: bool) -> None:
    """半径 100mm の球殻 STL を書き出す.

    Args:
        path: 出力先.
        inward: True なら面を反転して法線を内向きにする(空洞の内壁として正しい向き).
    """
    sphere = trimesh.creation.icosphere(subdivisions=4, radius=_SPHERE_RADIUS_MM)
    if inward:
        sphere.faces = sphere.faces[:, ::-1]
    sphere.export(str(path))


def _write_scene_yaml(
    path: Path,
    stl_path: Path,
    max_depth: int = _CONVERGED_MAX_DEPTH,
    radiance: float = _EMITTER_RADIANCE,
) -> Path:
    """積分球シーンの YAML を書き出してパスを返す."""
    cfg = {
        "name": "integrating_sphere",
        "description": "解析解と突き合わせる積分球",
        "geometry": {
            "meshes": [{"id": "shell", "path": str(stl_path), "unit": "mm", "material": "wall"}]
        },
        "materials": {"wall": {"type": "diffuse", "reflectance": [_REFLECTANCE] * 3}},
        "lights": [
            {
                "id": "src",
                "type": "area",
                "shape": "rectangle",
                "center": list(_EMITTER_CENTER_MM),
                "size": [_EMITTER_SIZE_MM, _EMITTER_SIZE_MM],
                "direction": [0.0, 0.0, -1.0],
                "radiance": [radiance] * 3,
            }
        ],
        # +Z 極の壁を見る.発光面は -Z を向いているのでここには直接光が届かない.
        "sensor": {
            "type": "perspective",
            "position": [0.0, 0.0, 50.0],
            "target": [0.0, 0.0, 90.0],
            "up": [0.0, 1.0, 0.0],
            "fov_deg": 30.0,
            "film": {"width": _FILM, "height": _FILM, "rfilter": "box"},
        },
        "render": {
            "variant": "llvm_ad_rgb",
            "integrator": "path",
            "max_depth": max_depth,
            "spp": _SPP,
        },
    }
    path.write_text(yaml.safe_dump(cfg, allow_unicode=True), encoding="utf-8")
    return path


def _render_mean_radiance(
    tmp_path: Path,
    inward: bool,
    max_depth: int,
    radiance: float = _EMITTER_RADIANCE,
) -> float:
    """積分球をレンダリングし,観測された平均放射輝度を返す."""
    tmp_path.mkdir(parents=True, exist_ok=True)
    stl_path = tmp_path / ("sphere_inward.stl" if inward else "sphere_outward.stl")
    _write_sphere_stl(stl_path, inward=inward)
    scene_yaml = _write_scene_yaml(
        tmp_path / "scene.yaml", stl_path, max_depth=max_depth, radiance=radiance
    )
    image = render_scene(load_scene_yaml(scene_yaml), project_root=tmp_path)
    return float(np.mean(image))


# ---------------------------------------------------------------------------
# 絶対放射量の検証
# ---------------------------------------------------------------------------


def test_integrating_sphere_radiance_matches_analytic(tmp_path):
    """積分球の壁面放射輝度が解析解と一致すること.

    単位系・反射率・発光面の放射束・センサ校正・多重反射の収束を一括で検証する.
    """
    measured = _render_mean_radiance(tmp_path, inward=True, max_depth=_CONVERGED_MAX_DEPTH)
    assert measured == pytest.approx(_analytic_shadowed_radiance(), rel=_TOLERANCE)


def test_rendered_radiance_scales_linearly_with_emitter_radiance(tmp_path):
    """発光面の radiance を 2 倍にすると観測値も 2 倍になること.

    輸送が線形であること(露出・トーンマップ等の非線形処理が紛れ込んでいないこと)を
    実レンダで確認する.シミュ値を実測にアンカーする校正はこの線形性に依存している.
    """
    base = _render_mean_radiance(tmp_path / "base", inward=True, max_depth=_CONVERGED_MAX_DEPTH)
    doubled = _render_mean_radiance(
        tmp_path / "doubled",
        inward=True,
        max_depth=_CONVERGED_MAX_DEPTH,
        radiance=_EMITTER_RADIANCE * 2.0,
    )
    assert doubled == pytest.approx(2.0 * base, rel=_TOLERANCE)


# ---------------------------------------------------------------------------
# 片面 BSDF による無警告の黒化に対する耐性
# ---------------------------------------------------------------------------


def test_flipped_normals_do_not_black_out_the_cavity(tmp_path):
    """法線が裏返った空洞でも黒化せず解析解に一致すること.

    Mitsuba の diffuse は片面 BSDF なので,twosided ラップが無いと
    ここは 0.0 になる(エラーも警告も出ない).
    materials.build_scene_bsdf_dict のラップが効いていることの実証.
    """
    measured = _render_mean_radiance(tmp_path, inward=False, max_depth=_CONVERGED_MAX_DEPTH)
    assert measured > 0.0, "空洞が真っ黒になっている(片面 BSDF に逆戻りしている)"
    assert measured == pytest.approx(_analytic_shadowed_radiance(), rel=_TOLERANCE)


def test_flipped_and_correct_normals_agree(tmp_path):
    """法線の向きによらず同じ結果になること.

    不透明な壁の反射は面の向き付けという描画上の規約に依存してはならない.
    """
    inward = _render_mean_radiance(tmp_path / "a", inward=True, max_depth=_CONVERGED_MAX_DEPTH)
    outward = _render_mean_radiance(tmp_path / "b", inward=False, max_depth=_CONVERGED_MAX_DEPTH)
    assert outward == pytest.approx(inward, rel=_TOLERANCE)


# ---------------------------------------------------------------------------
# 多重反射の打ち切りバイアス
# ---------------------------------------------------------------------------


def test_shallow_max_depth_underestimates_radiance(tmp_path):
    """max_depth が浅いと解析解を有意に下回ること.

    scene/audit.py の打ち切り警告が実在の問題に対する警告であることの裏付け.
    rho=0.8・max_depth=8 では収束値の約 74% しか拾えない(TEST_13 実測).
    """
    shallow = _render_mean_radiance(tmp_path, inward=True, max_depth=8)
    analytic = _analytic_shadowed_radiance()
    assert shallow < analytic * 0.9


def test_deeper_max_depth_converges_towards_analytic(tmp_path):
    """max_depth を深くすると解析解に近づくこと(単調な収束)."""
    analytic = _analytic_shadowed_radiance()
    shallow = _render_mean_radiance(tmp_path / "a", inward=True, max_depth=8)
    deep = _render_mean_radiance(tmp_path / "b", inward=True, max_depth=_CONVERGED_MAX_DEPTH)
    assert abs(deep - analytic) < abs(shallow - analytic)