Newer
Older
SmTIAS-LightSim / scripts / generate_cube_housing.py
# TEST_16(立方体キャビティ実験)用: 評価板の矩形輪郭に整合する立方体筐体メッシュ生成.
#
# 球(generate_sphere_loft_housing.py)では「矩形の評価板 → 円形の球断面」を繋ぐ
# ロフトが必要で,その接合の幾何(板直下の喉・曲率不連続)が TEST_15 の主要な
# 交絡要因だった.立方体は評価板もキャビティ断面も矩形なのでロフトが構造的に不要になり,
# 上面はそのまま「板の縁から外へ広がる水平な天井」=TEST_15 で最良だった radial ロフトが
# 目指した形になる.
"""立方体キャビティ筐体(評価板輪郭に整合する閉ソリッド)STL を生成する.

構成(上から下へ,単位 mm):
    1. リム矩形: 評価板輪郭から rim_inset だけ内側の矩形.z=rim_z(板の実体に
       食い込む高さ)に置き,板とのオーバーラップで隙間なくシールする
    2. 天井: z=rim_z の水平面.リム矩形から内壁正方形(L×L)まで広がる
    3. 側壁 4 枚: z=rim_z から z=rim_z-L まで
    4. 底面: z=rim_z-L の水平面

内面と外面(thickness オフセット.外面の天板は rim_z+thickness に持ち上げて
内面天井との共面を避ける)を生成し,リム帯で縫合した閉ソリッドを出力する.

球版と違い全ての面が平面なので,三角形分割は物理(法線・交差)に影響しない.
分割解像度はレイの数値的なすり抜け回避と可視化の見やすさのためだけにある.
周方向サンプルは半ステップずらし,メッシュエッジが XZ / YZ 対称面に乗らないようにする
(球スイープの遮蔽検証で実測したすり抜けの再発防止).

使用例:
    uv run python scripts/generate_cube_housing.py --size 200 \
        --output configs/meshes/立方体_L200.stl
"""

from __future__ import annotations

import argparse
from pathlib import Path

import numpy as np
import trimesh
from generate_sphere_loft_housing import _rect_radius, _stitch_rings


def _corner_angles(half_x: float, half_y: float) -> np.ndarray:
    """矩形(半幅 half_x, half_y)の 4 隅の方位角を返す."""
    t = float(np.arctan2(half_y, half_x))
    return np.array([t, np.pi - t, np.pi + t, 2.0 * np.pi - t])


def _build_theta(n_theta: int, rect_halves: list[tuple[float, float]]) -> np.ndarray:
    """周方向サンプル角を返す(均等分割 + 全矩形の隅の厳密角).

    均等分割は半ステップずらして 0 / 90 / 180 / 270 度を避ける(対称面上の
    メッシュエッジをレイが数値的にすり抜ける現象の回避).隅の角度を明示的に
    含めることで,矩形の角が theta 分割で面取りされず厳密に表現される.
    """
    base = (np.arange(n_theta) + 0.5) * (2.0 * np.pi / n_theta)
    corners = np.concatenate([_corner_angles(hx, hy) for hx, hy in rect_halves])
    return np.unique(np.concatenate([base, corners]))


def _build_surface(
    theta: np.ndarray,
    rim_r: np.ndarray,
    wall_r: np.ndarray,
    rim_z: float,
    top_z: float,
    bottom_z: float,
    n_ceiling: int,
    n_wall: int,
    n_bottom: int,
) -> tuple[np.ndarray, np.ndarray]:
    """リム矩形 →(襟)→ 天井 → 側壁 → 底面の 1 枚曲面を生成する.

    top_z > rim_z のとき,リムから top_z まで垂直に立ち上がる襟を先頭に挿入する
    (外面用.天板を持ち上げて内面天井 z=rim_z との共面=ゼロ厚を避ける).

    Returns:
        (vertices, rings): vertices は (N, 3),rings は (n_ring, n_theta) の頂点 index.
        リングは上(リム)から下(底面中心直前)の順.末尾に底面中心の極頂点を 1 つ追加.
    """
    rings_rho: list[np.ndarray] = []
    rings_z: list[float] = []
    if top_z > rim_z:
        for z in np.linspace(rim_z, top_z, 2, endpoint=False):
            rings_rho.append(rim_r)
            rings_z.append(float(z))
    for t in np.linspace(0.0, 1.0, n_ceiling, endpoint=False):
        rings_rho.append((1.0 - t) * rim_r + t * wall_r)
        rings_z.append(top_z)
    for z in np.linspace(top_z, bottom_z, n_wall, endpoint=False):
        rings_rho.append(wall_r)
        rings_z.append(float(z))
    for t in np.linspace(0.0, 1.0, n_bottom, endpoint=False):
        rings_rho.append((1.0 - t) * wall_r)
        rings_z.append(bottom_z)

    n_ring, n_theta = len(rings_rho), len(theta)
    vertices = np.empty((n_ring * n_theta + 1, 3), dtype=np.float64)
    rings = np.arange(n_ring * n_theta).reshape(n_ring, n_theta)
    for i, (rho, z) in enumerate(zip(rings_rho, rings_z, strict=True)):
        vertices[rings[i], 0] = rho * np.cos(theta)
        vertices[rings[i], 1] = rho * np.sin(theta)
        vertices[rings[i], 2] = z
    vertices[-1] = (0.0, 0.0, bottom_z)
    return vertices, rings


def build_cube_housing(
    size: float,
    thickness: float,
    board_half_x: float,
    board_half_y: float,
    rim_inset: float,
    rim_z: float,
    n_theta: int = 256,
    n_ceiling: int = 24,
    n_wall: int = 64,
    n_bottom: int = 32,
) -> trimesh.Trimesh:
    """立方体キャビティ筐体の閉ソリッド(内面+外面+リム帯)を生成する.

    Args:
        size: キャビティ内寸の一辺 L(mm).内壁は x, y = ±L/2,深さは rim_z から L.
        thickness: 殻厚(mm).外面は内面をおおむね thickness オフセットした形.
        board_half_x: 評価板輪郭の X 半幅(mm).
        board_half_y: 評価板輪郭の Y 半幅(mm).
        rim_inset: リム矩形を板輪郭から内側へ寄せる量(mm).板とのオーバーラップ代.
        rim_z: リム矩形と天井の Z 座標(mm).板の実体に食い込む高さに置いてシールする.
        n_theta: 周方向の均等分割数(隅の厳密角が追加されるため実分割数は少し増える).
        n_ceiling: 天井(リム矩形 → 内壁正方形)の分割リング数.
        n_wall: 側壁の Z 方向リング数.
        n_bottom: 底面(内壁正方形 → 中心)の分割リング数.

    Returns:
        watertight な筐体ソリッドメッシュ(mm 単位).

    Raises:
        ValueError: 評価板が上面に収まらない size が指定された場合,または
            生成メッシュが watertight にならなかった場合.
    """
    if size / 2.0 <= max(board_half_x, board_half_y):
        msg = (
            f"size={size}mm では評価板({2 * board_half_x:g}×{2 * board_half_y:g}mm)が"
            f"上面に収まりません(size > {2 * max(board_half_x, board_half_y):g} が必要)"
        )
        raise ValueError(msg)

    hx_in = board_half_x - rim_inset
    hy_in = board_half_y - rim_inset
    half_in = size / 2.0
    half_out = half_in + thickness
    theta = _build_theta(
        n_theta,
        [
            (hx_in, hy_in),
            (half_in, half_in),
            (hx_in + thickness, hy_in + thickness),
            (half_out, half_out),
        ],
    )

    v_in, rings_in = _build_surface(
        theta,
        _rect_radius(theta, hx_in, hy_in),
        _rect_radius(theta, half_in, half_in),
        rim_z,
        rim_z,
        rim_z - size,
        n_ceiling,
        n_wall,
        n_bottom,
    )
    v_out, rings_out = _build_surface(
        theta,
        _rect_radius(theta, hx_in + thickness, hy_in + thickness),
        _rect_radius(theta, half_out, half_out),
        rim_z,
        rim_z + thickness,
        rim_z - size - thickness,
        n_ceiling,
        n_wall,
        n_bottom,
    )

    f_in = _stitch_rings(rings_in, len(v_in) - 1)
    f_out = _stitch_rings(rings_out, len(v_out) - 1)

    offset = len(v_in)
    vertices = np.vstack([v_in, v_out])
    # リム帯: 内面上端リングと外面上端リングを z=rim_z の平面内で接続して閉じる
    top_in, top_out = rings_in[0], rings_out[0] + offset
    n_ring_theta = rings_in.shape[1]
    band: list[tuple[int, int, int]] = []
    for j in range(n_ring_theta):
        k = (j + 1) % n_ring_theta
        band.append((top_in[j], top_out[j], top_out[k]))
        band.append((top_in[j], top_out[k], top_in[k]))
    faces = np.vstack([f_in, f_out + offset, np.asarray(band, dtype=np.int64)])

    mesh = trimesh.Trimesh(vertices=vertices, faces=faces, process=False)
    mesh.fix_normals()  # 閉ソリッドとして法線を外向きに統一(内面は自動的にキャビティ向き)
    if not mesh.is_watertight:
        raise ValueError("生成メッシュが watertight になりませんでした(縫合パラメータを確認)")
    return mesh


def _parse_args(argv: list[str] | None = None) -> argparse.Namespace:
    """コマンドライン引数をパースして返す."""
    parser = argparse.ArgumentParser(
        prog="generate_cube_housing",
        description="評価板輪郭に整合する立方体キャビティ筐体 STL(mm 単位)を生成する.",
    )
    parser.add_argument(
        "--size", type=float, default=200.0, help="キャビティ内寸 L(mm).既定 200"
    )
    parser.add_argument("--thickness", type=float, default=3.0, help="殻厚(mm).既定 3")
    parser.add_argument(
        "--board-half",
        type=float,
        nargs=2,
        default=(49.5, 51.0),
        metavar=("HX", "HY"),
        help="評価板輪郭の半幅 X Y(mm).既定 49.5 51(99×102mm 板)",
    )
    parser.add_argument(
        "--rim-inset",
        type=float,
        default=1.5,
        help="リムを板輪郭から内側へ寄せる量(mm).既定 1.5",
    )
    parser.add_argument(
        "--rim-z",
        type=float,
        default=0.0,
        help="リムと天井の Z 座標(mm).既定 0(板裏面 -4.5 より 4.5mm 食い込みシール)",
    )
    parser.add_argument(
        "--output",
        type=Path,
        default=Path("configs/meshes/立方体_L200.stl"),
        help="出力 STL パス.既定 configs/meshes/立方体_L200.stl",
    )
    return parser.parse_args(argv)


def main(argv: list[str] | None = None) -> int:
    """筐体を生成して STL 出力し,検証情報を表示する."""
    args = _parse_args(argv)
    mesh = build_cube_housing(
        size=args.size,
        thickness=args.thickness,
        board_half_x=args.board_half[0],
        board_half_y=args.board_half[1],
        rim_inset=args.rim_inset,
        rim_z=args.rim_z,
    )
    args.output.parent.mkdir(parents=True, exist_ok=True)
    mesh.export(args.output)
    print(f"出力: {args.output}")
    print(f"faces: {len(mesh.faces)}, watertight: {mesh.is_watertight}")
    print(f"bounds (mm): {mesh.bounds.round(2).tolist()}")
    print(f"体積 (mm^3): {mesh.volume:.0f}(正値なら法線外向きで正常)")
    return 0


if __name__ == "__main__":
    raise SystemExit(main())