"""球体筐体の半径スイープ（キャビティサイズが照明均一性に与える効果の定量化）.

TEST_11 / TEST_12 の球ロフト筐体について，**評価板のサイズとレンダ条件を据え置いたまま
球の内半径だけを変える**．半径が大きいほど球面積に対する評価板の占有率が下がるため，
理想積分球（開口比ゼロ）に漸近するはずで，その効きを CoV / angular_ratio 等で測る．

幾何のパラメータ化（重要）:
    ロフト接続円（z = CUT_Z の断面円 CUT_R）を**固定**し，球中心 z を動かして球を深くする．

        center_z = CUT_Z - sqrt(radius^2 - CUT_R^2)
        pole_z   = center_z - radius
        led_z    = pole_z + LED_POLE_OFFSET

    球中心を固定して半径を変える方式は採らない．断面円まで一緒に大きくなり，評価板まわりの
    スカート（ロフト）形状が半径ごとに変わって「球の大きさ」の効果と交絡するため．
    radius = 100 のとき center_z = -92 となり，既存 configs/scene_sphere_r100_loft_disk.yaml
    と厳密に一致する（パイプラインの検証点）．

固定条件（既存レンダと完全に同一に保つ）:
    材質  : 既存の blendbsdf rho=0.85 / omega=0.05 / alpha=0.1
            （TEST_14 の校正値 rho=0.81 純 diffuse では**なく**，TEST_11/12 と直接比較するため）
    光源  : 軸上の円板 radius=1.7662mm（面積 9.8mm^2 保存）．評価板と反対側（-Z 極側）
    レンダ: spp=1024 / max_depth=16 / 512^2 / llvm_ad_rgb / path / orthographic (50x52)
    評価  : configs/eval_baseline_dng_target_shape.yaml（shape_metrics 込み）

LED 構成（--variant）:
    board_baffle（既定）: LED を **-Z 極の壁面に密着**させ **評価板向き (0,0,+1)** にし，
        直前にバッフル円板を置いて
        評価板への直接光を幾何的にゼロにする．バッフルは板の全点から LED を隠す最小半径
        （baffle_radius）に安全率を掛けたもので，半径が大きいほど板の見込み角が小さくなり
        バッフルも小さくて済む．
    axial_wall: LED を壁向き (0,0,-1) にした従来構成（TEST_11 / TEST_12 と地続き）．
        バッフル無しでも板は発光面の裏側にあるため直接光はゼロ．

絶対輝度について:
    評価設定は mean アンカー正規化（target_mean=207.82）をかけるので metrics.csv の mean は
    常に 207.82 になり絶対的な明るさの情報が消える．球が大きいほど暗くなる（おおむね 1/r^2）
    のは実機設計上の重要なトレードオフなので，luminance.npy の**生の mean** を別途記録する．

実行方法:
    uv run python scripts/sweep_sphere_radius.py --prepare-only   # メッシュ + シーン YAML のみ
    uv run python scripts/sweep_sphere_radius.py                  # 上記 + レンダ + 評価
"""

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 trimesh
import yaml

sys.path.insert(0, str(Path(__file__).resolve().parent))
sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))

from generate_sphere_loft_housing import build_loft_housing  # noqa: E402

from sm_tias_lightsim.evaluation.roi_presets import SIM_ROI, extract_roi_xywh  # noqa: E402

ROOT = Path(__file__).resolve().parents[1]

TEMPLATE_YAML = ROOT / "configs" / "scene_sphere_r100_loft_disk.yaml"
EVAL_CONFIG_YAML = ROOT / "configs" / "eval_baseline_dng_target_shape.yaml"
TMP_DIR = ROOT / "configs" / "_tmp_sphere_radius"
MESH_DIR = TMP_DIR / "meshes"
OUTPUT_BASE = ROOT / "output" / "sphere_radius_sweep"

# 幾何定数（既存 r=100 / center_z=-92 構成から導出．変更しないこと）
CUT_Z = -40.0
CUT_R = float(np.sqrt(100.0**2 - 52.0**2))  # = 85.4166 ロフト接続円の半径

# キャビティ形状の指定（--cavity）
#   fixed_cut : 接続円（CUT_R, CUT_Z）を固定して球中心だけ動かす．3 次ロフト．
#       半径を上げると接続面が球の頂点に近づき，r>=120 では**球の頂点が評価板より下**に
#       もぐる（＝ほぼ閉じた玉に細い首が生えた形）．接合の傾きは 0.25→2.75 と 11 倍に
#       急峻化し，接合部の曲率比も 9.3→0.11→1.34 と暴れる．
#   similar_g2: キャビティ全体を半径に相似させ（center_z=-0.92r, cut_z=-0.40r），
#       ロフトを 5 次 Hermite にして接合部で**傾きと曲率の両方**を球に一致させる．
#       接合の傾きは全半径で 0.609 一定，曲率比は厳密に 1.00．球の頂点は常に板より上
#       （z=+7〜+20mm）に留まり，どの半径でも「お椀にフタ」の自然な形になる．
#       比は r=100 の既存構成（center_z=-92, cut_z=-40）をそのまま無次元化した値なので，
#       r=100 では接続円・球中心が既存と一致する（ロフトの曲線だけ 5 次になる）．
CENTER_Z_RATIO = -0.92
CUT_Z_RATIO = -0.40
CAVITIES = ("fixed_cut", "similar_g2")

# LED を -Z 極から何 mm 手前に置くか（構成ごとに異なる）
# axial_wall: TEST_11 の分散低減実測（NEE 可視領域が球面積の 25%）に基づく 50mm
# board_baffle: 壁面に密着させる．厳密に 0 にすると平らな発光円板が湾曲した壁面を貫くため
#   （円板半径 1.7662mm でのサジッタは r=88 でも 0.018mm）最小限だけ浮かせる．
#   極の壁面法線は +Z＝評価板向きなので，密着と「評価板を向く」が同時に成り立つ
LED_WALL_STANDOFF = 0.5
LED_POLE_OFFSET = {"axial_wall": 50.0, "board_baffle": LED_WALL_STANDOFF}

# 筐体生成の既定値（generate_sphere_loft_housing の CLI 既定と同一）
THICKNESS = 3.0
BOARD_HALF = (49.5, 51.0)
RIM_INSET = 1.5
RIM_Z = 0.0

# LED 円板の半径（テンプレート YAML と一致させること）
LED_RADIUS = 1.7662
# 評価板の最遠点（四隅）までの軸からの距離．バッフル寸法の律速になる
BOARD_CORNER = float(np.hypot(*BOARD_HALF))
# 評価板のキャビティ側の面の Z（照明ムラ補正ボード.stl の最下面）．
# センサは z=-50 から +Z を向いてこの面を撮るので，遮蔽すべき対象はこの面である
BOARD_SURFACE_Z = -4.5

# バッフル（board_baffle 構成）: LED の直前に置く円板で評価板への直接光を完全に塞ぐ
BAFFLE_GAP = 10.0  # LED 発光面からバッフル中心面までの距離（mm）
BAFFLE_THICKNESS = 1.0  # 円板の厚み（mm）
BAFFLE_MARGIN = 1.05  # 幾何的な最小半径に対する安全率

# LED の向き．axial_wall: 壁（-Z 極）向き＝従来 / board_baffle: 評価板（+Z）向き
VARIANTS = ("axial_wall", "board_baffle")

# r=170 は元の TIAS のサイズ，r=100 は既存結果との突合用アンカー．
RADII = [88.0, 95.0, 100.0, 110.0, 120.0, 135.0, 150.0, 170.0, 200.0, 250.0]

SUMMARY_COLUMNS = [
    "radius_mm", "center_z_mm", "pole_z_mm", "led_z_mm",
    "cavity_area_mm2", "board_area_ratio", "elapsed_s",
    "raw_mean_roi", "raw_mean_full", "raw_rel_std_roi",
    "cov", "center_periphery_ratio", "gradient_magnitude", "radial_min_max_ratio",
    "angular_ratio", "pca_eccentricity", "intensity_centroid_offset_norm",
    "max_min_ratio",
]


def center_z(radius: float, cavity: str = "fixed_cut") -> float:
    """球中心の Z 座標を返す（球は -Z 側に伸びる）."""
    if cavity == "similar_g2":
        return CENTER_Z_RATIO * radius
    if radius <= CUT_R:
        raise ValueError(f"radius は CUT_R={CUT_R:.3f}mm より大きい必要があります（指定 {radius}）")
    return CUT_Z - float(np.sqrt(radius**2 - CUT_R**2))


def cut_z(radius: float, cavity: str = "fixed_cut") -> float:
    """ロフトと球の接続平面の Z 座標を返す."""
    return CUT_Z_RATIO * radius if cavity == "similar_g2" else CUT_Z


def cut_radius(radius: float, cavity: str = "fixed_cut") -> float:
    """接続円の半径を返す（球と接続平面の交線）."""
    return float(np.sqrt(radius**2 - (cut_z(radius, cavity) - center_z(radius, cavity)) ** 2))


def loft_blend(cavity: str) -> str:
    """キャビティ形状に対応するロフトのブレンド次数を返す."""
    return "quintic" if cavity == "similar_g2" else "cubic"


def led_z(radius: float, variant: str, cavity: str = "fixed_cut") -> float:
    """LED 円板の Z 座標（-Z 極から構成ごとのオフセットだけ手前）を返す."""
    return center_z(radius, cavity) - radius + LED_POLE_OFFSET[variant]


def job_name(radius: float) -> str:
    """半径からジョブ名（ディレクトリ名）を作る."""
    return f"r{radius:g}"


def baffle_radius(
    radius: float, variant: str = "board_baffle", cavity: str = "fixed_cut"
) -> float:
    """評価板の全点から LED を隠すのに必要なバッフル円板の半径（mm）を返す.

    LED（半径 LED_RADIUS の円板，z = led_z）から評価板のキャビティ側の面
    （z = BOARD_SURFACE_Z）へ向かう直線が z = led_z + BAFFLE_GAP の平面を横切る位置は，
    LED 上の点 Q と板上の点 P の内分点なので，
    f = BAFFLE_GAP / (BOARD_SURFACE_Z - led_z) として

        |C| <= (1 - f) * LED_RADIUS + f * BOARD_CORNER

    で厳密に上から押さえられる（等号は Q と P が同じ方位にあるとき）．
    これに安全率 BAFFLE_MARGIN を掛けたものを採用する．
    """
    f = BAFFLE_GAP / (BOARD_SURFACE_Z - led_z(radius, variant, cavity))
    return BAFFLE_MARGIN * ((1.0 - f) * LED_RADIUS + f * BOARD_CORNER)


def build_baffle_mesh(
    radius: float, dest: Path, variant: str = "board_baffle", cavity: str = "fixed_cut"
) -> trimesh.Trimesh:
    """LED 直前のバッフル円板（薄い円柱ソリッド）を生成して書き出す.

    円板は半セクタだけ回転させる．無回転だと蓋の三角形ファンのエッジが XZ / YZ 平面に
    完全一致し，その平面内を通るレイがエッジ上の交差になって数値的にすり抜ける
    （verify_baffle_occlusion.py で実測した現象）．物理的な影響は無いが，
    遮蔽の検証を明確にするため対称面上にエッジを置かない．
    """
    sections = 128
    mesh = trimesh.creation.cylinder(
        radius=baffle_radius(radius, variant, cavity), height=BAFFLE_THICKNESS, sections=sections
    )
    mesh.apply_transform(trimesh.transformations.rotation_matrix(np.pi / sections, [0, 0, 1]))
    mesh.apply_translation([0.0, 0.0, led_z(radius, variant, cavity) + BAFFLE_GAP])
    dest.parent.mkdir(parents=True, exist_ok=True)
    mesh.export(dest)
    return mesh


def build_mesh(radius: float, dest: Path, cavity: str = "fixed_cut") -> trimesh.Trimesh:
    """指定半径・指定キャビティ形状の筐体 STL を生成して書き出す."""
    mesh = build_loft_housing(
        radius=radius,
        thickness=THICKNESS,
        center_z=center_z(radius, cavity),
        board_half_x=BOARD_HALF[0],
        board_half_y=BOARD_HALF[1],
        rim_inset=RIM_INSET,
        rim_z=RIM_Z,
        cut_z=cut_z(radius, cavity),
        blend=loft_blend(cavity),
    )
    dest.parent.mkdir(parents=True, exist_ok=True)
    mesh.export(dest)
    return mesh


def build_scene_yaml(
    radius: float, mesh_path: Path, dest: Path, variant: str, cavity: str,
    baffle_path: Path | None = None,
) -> None:
    """テンプレートから当該半径・当該構成のシーン YAML を書き出す.

    材質・センサ・レンダ設定・LED の半径/radiance は一切変更しない．変えるのは
    筐体メッシュ・LED の Z・LED の向き（構成による）・バッフルの有無だけ．
    メッシュの相対パスは YAML 自身のディレクトリ基準で解決されるため絶対パス化する．
    """
    config = copy.deepcopy(yaml.safe_load(TEMPLATE_YAML.read_text(encoding="utf-8")))
    config["name"] = f"sphere_radius_{variant}_{cavity}_{job_name(radius)}"
    if variant == "board_baffle":
        aim = "壁面に密着し評価板向き（+Z）．バッフルで直接光を遮断"
    else:
        aim = "壁向き（-Z）で評価板への直接光が幾何的にゼロ"
    config["description"] = (
        f"球体筐体の半径スイープ [{variant} / {cavity}]: 内半径 {radius:g}mm"
        f"（球中心 z={center_z(radius, cavity):.3f}・接続円 {cut_radius(radius, cavity):.3f}mm"
        f" @ z={cut_z(radius, cavity):.3f}・ロフト {loft_blend(cavity)}）．"
        f"LED は極から {LED_POLE_OFFSET[variant]:g}mm 手前，{aim}．"
    )
    for mesh in config["geometry"]["meshes"]:
        if mesh["id"] == "sphere_housing":
            mesh["path"] = str(mesh_path.resolve())
        else:
            mesh["path"] = str((TEMPLATE_YAML.parent / mesh["path"]).resolve())
    if baffle_path is not None:
        config["geometry"]["meshes"].append({
            "id": "baffle",
            "path": str(baffle_path.resolve()),
            "unit": "mm",
            "material": "printer_white",
        })
    for light in config["lights"]:
        light["center"] = [0.0, 0.0, led_z(radius, variant, cavity)]
        light["direction"] = [0.0, 0.0, 1.0 if variant == "board_baffle" else -1.0]
    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 cavity_geometry(radius: float, cavity: str = "fixed_cut") -> dict[str, float]:
    """キャビティ内面積と，評価板がそれに占める割合（開口比）を返す.

    内面積 = 球のお椀（接続面より下の球冠）+ ロフト部側面（錐台近似）+ 評価板．
    開口比は積分球としての理想度の目安（小さいほど理想に近い）．
    """
    z_cut = cut_z(radius, cavity)
    r_cut = cut_radius(radius, cavity)
    cap_height = z_cut - (center_z(radius, cavity) - radius)  # 極から接続面までの高さ
    sphere_area = 2.0 * np.pi * radius * cap_height
    rect_perimeter = 2.0 * (2 * BOARD_HALF[0] + 2 * BOARD_HALF[1])
    slant = float(np.hypot(abs(z_cut - RIM_Z), r_cut - float(np.hypot(*BOARD_HALF)) * 0.5))
    loft_area = 0.5 * (2 * np.pi * r_cut + rect_perimeter) / 2.0 * slant
    board_area = (2 * BOARD_HALF[0]) * (2 * BOARD_HALF[1])
    total = sphere_area + loft_area + board_area
    return {
        "cavity_area_mm2": round(total, 1),
        "board_area_ratio": round(board_area / total, 5),
    }


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))
    keys = ("cov", "center_periphery_ratio", "gradient_magnitude", "radial_min_max_ratio",
            "angular_ratio", "pca_eccentricity", "intensity_centroid_offset_norm",
            "max_min_ratio")
    return {k: (float(row[k]) if row.get(k) else float("nan")) for k in keys}


def raw_stats(luminance_npy: Path) -> dict[str, float]:
    """正規化前の生の輝度統計（絶対的な明るさ）を返す.

    metrics.csv は mean アンカー正規化後なので mean が常に 207.82 になり，
    「球が大きいほど暗くなる」という設計上重要な情報が失われる．ここで拾っておく．
    """
    arr = np.load(luminance_npy).astype(np.float64)
    if arr.ndim == 3:
        arr = 0.2126 * arr[..., 0] + 0.7152 * arr[..., 1] + 0.0722 * arr[..., 2]
    roi = extract_roi_xywh(arr, SIM_ROI)
    return {
        "raw_mean_roi": float(roi.mean()),
        "raw_mean_full": float(arr.mean()),
        "raw_rel_std_roi": float(roi.std() / roi.mean()),
    }


def prepare(radii: list[float], variant: str, cavity: str, log) -> dict[float, Path]:
    """全半径について筐体・バッフルのメッシュとシーン YAML を生成し，watertight を検証する."""
    scenes: dict[float, Path] = {}
    log(f"=== 幾何生成（構成 {variant} / キャビティ {cavity}）===")
    for radius in radii:
        mesh_path = MESH_DIR / f"球ロフト_{cavity}_{job_name(radius)}.stl"
        mesh = build_mesh(radius, mesh_path, cavity)
        if not mesh.is_watertight:
            raise ValueError(f"r={radius}: watertight でないメッシュが生成されました")
        baffle_path = None
        baffle_note = ""
        if variant == "board_baffle":
            baffle_path = MESH_DIR / f"バッフル_{cavity}_{job_name(radius)}.stl"
            baffle = build_baffle_mesh(radius, baffle_path, variant, cavity)
            if not baffle.is_watertight:
                raise ValueError(f"r={radius}: バッフルが watertight になりませんでした")
            baffle_note = (f" バッフル R={baffle_radius(radius, variant, cavity):5.2f}mm "
                           f"z={led_z(radius, variant, cavity) + BAFFLE_GAP:8.2f}")
        scene_yaml = TMP_DIR / f"scene_{variant}_{cavity}_{job_name(radius)}.yaml"
        build_scene_yaml(radius, mesh_path, scene_yaml, variant, cavity, baffle_path)
        scenes[radius] = scene_yaml
        geo = cavity_geometry(radius, cavity)
        log(f"  r={radius:6.1f} center_z={center_z(radius, cavity):9.3f} "
            f"cut_z={cut_z(radius, cavity):8.2f} 接続円={cut_radius(radius, cavity):7.2f} "
            f"LED_z={led_z(radius, variant, cavity):9.3f} "
            f"開口比={geo['board_area_ratio'] * 100:5.2f}%{baffle_note}")
    return scenes


def run_dir_name(variant: str, cavity: str) -> str:
    """出力ディレクトリ名を返す.

    既定キャビティ（fixed_cut）は変更前に実施した board_baffle 10 条件の結果が
    output/sphere_radius_sweep/board_baffle/ に既にあるため，そのまま再利用できるよう
    構成名だけを使う．新しいキャビティはサフィックスを付けて分ける．
    """
    return variant if cavity == "fixed_cut" else f"{variant}_{cavity}"


def run_job(
    radius: float, scene_yaml: Path, variant: str, cavity: str, log
) -> dict[str, float] | None:
    """1 半径をレンダ→評価まで実行して集約行を返す（失敗時 None）.

    既に metrics.csv と luminance.npy が揃っているジョブは再レンダせず再利用する．
    """
    name = job_name(radius)
    out_dir = OUTPUT_BASE / run_dir_name(variant, cavity) / name
    row: dict[str, float] = {
        "radius_mm": radius,
        "center_z_mm": round(center_z(radius, cavity), 3),
        "pole_z_mm": round(center_z(radius, cavity) - radius, 3),
        "led_z_mm": round(led_z(radius, variant, cavity), 3),
        "elapsed_s": 0.0,
    }
    row.update(cavity_geometry(radius, cavity))
    if (out_dir / "metrics.csv").exists() and (out_dir / "luminance.npy").exists():
        log(f"--- {name}: 既存結果を再利用")
    else:
        log(f"--- {name}: レンダ開始（center_z={center_z(radius, cavity):.2f} "
            f"LED_z={led_z(radius, variant, cavity):.2f}）")
        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)],
        ):
            # encoding を明示する．既定では Windows の cp932 で復号され，
            # sm-render が出す日本語ログで UnicodeDecodeError になって
            # 失敗時のエラー本文が読めなくなる（実測）．
            result = subprocess.run(
                cmd, cwd=ROOT, capture_output=True, text=True,
                encoding="utf-8", errors="replace",
            )
            if result.returncode != 0:
                log(f"  !! 失敗 ({cmd[2]}): {result.stderr[-800:]}")
                return None
        row["elapsed_s"] = round(time.time() - start, 1)
    row.update(read_metrics_csv(out_dir / "metrics.csv"))
    row.update(raw_stats(out_dir / "luminance.npy"))
    log(f"  -> {row['elapsed_s']:.0f}s CoV={row['cov']:.4f} "
        f"angular_ratio={row['angular_ratio']:.4f} C/P={row['center_periphery_ratio']:.4f} "
        f"勾配={row['gradient_magnitude']:.2f}% 生mean={row['raw_mean_roi']:.4f}")
    return row


def write_summary(rows: list[dict[str, float]], variant: str, cavity: str) -> None:
    """集約 CSV を半径昇順で書き直す（途中終了でも最新状態が残る）."""
    summary_csv = OUTPUT_BASE / run_dir_name(variant, cavity) / "summary.csv"
    summary_csv.parent.mkdir(parents=True, exist_ok=True)
    with summary_csv.open("w", newline="", encoding="utf-8") as fh:
        writer = csv.DictWriter(fh, fieldnames=SUMMARY_COLUMNS)
        writer.writeheader()
        for row in sorted(rows, key=lambda r: r["radius_mm"]):
            writer.writerow({k: row.get(k, "") for k in SUMMARY_COLUMNS})


def main() -> int:
    """半径スイープを実行するメインエントリポイント."""
    parser = argparse.ArgumentParser(description="球体筐体の半径スイープ")
    parser.add_argument("--prepare-only", action="store_true",
                        help="メッシュとシーン YAML の生成だけ行う（3D 配置確認用）")
    parser.add_argument("--radii", type=float, nargs="*", default=None,
                        help="半径リスト（mm）．省略時は既定グリッド")
    parser.add_argument("--variant", choices=VARIANTS, default="board_baffle",
                        help="LED 構成．board_baffle=評価板向き+バッフル / axial_wall=壁向き")
    parser.add_argument("--cavity", choices=CAVITIES, default="similar_g2",
                        help="キャビティ形状．similar_g2=相似+5次G2接合 / fixed_cut=接続円固定")
    args = parser.parse_args()

    radii = args.radii if args.radii else RADII
    variant, cavity = args.variant, args.cavity
    log_path = ROOT / "output" / "_sweep_sphere_radius.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"=== 球体筐体 半径スイープ（構成 {variant} / キャビティ {cavity}）===")
    scenes = prepare(radii, variant, cavity, log)
    if args.prepare_only:
        log(f"--prepare-only のため終了．シーン YAML: {TMP_DIR}")
        log_file.close()
        return 0

    rows: list[dict[str, float]] = []
    for radius in radii:
        row = run_job(radius, scenes[radius], variant, cavity, log)
        if row is not None:
            rows.append(row)
            write_summary(rows, variant, cavity)

    log(f"\n=== 完了: {len(rows)} 条件 ===")
    log(f"  集約 CSV: {OUTPUT_BASE / run_dir_name(variant, cavity) / 'summary.csv'}")
    log_file.close()
    return 0


if __name__ == "__main__":
    sys.exit(main())