弾丸が壁をすり抜けるとき — 連続衝突検出と包含関数による二分法

一軸を1e-4まで詰めるには半分割が14回要ります。印刷された53箱は14×14の総当たりより少ないです。

60fpsで動くシミュレーションでは、速い弾丸は壁をそのまますり抜けます。あるフレームでは壁の左、次のフレームでは右。どちらの瞬間でも重なっていないので、衝突検出器は何も起きなかったと報告します。本記事はその「あいだ」を検査する連続衝突検出（CCD、二物体の線形軌道全体で最初の接触時刻を求める手法）を扱います。共面条件によって衝突を多項式の求根へと変える過程、そして浮動小数点でも衝突を決して見逃さない包含関数による二分法を実装します。

このテーマは、流体構造連成（FSI）、粒子法、メッシュベースの接触シミュレーションなど、「フレーム間で物体が高速に動く」あらゆるコードに直結します。

二つのフレームの間をすり抜けた弾丸

静的衝突検出は一瞬のスナップショットしか見ません。二つの三角形が今重なっているかを問います。遅い物体には十分です。問題は速度です。

タイムステップ $\Delta t$ のあいだに物体が自身の厚みより遠くへ動くと、スナップショットの間をすり抜けます。この現象をトンネリングと呼びます。さらに厄介なのは復旧できない点です。衝突を一度見逃すと二物体は貫通したまま残り、次のステップではどの面がどちら側かさえ分からなくなります。

point travel を上げると、点が1ステップで横切る距離が伸びます。灰色の棒がその壁です。

point travel: 8.0wall x: 5.0wall drift: 0.0contact radius: 0.18

discrete (endpoints only): MISScontinuous TOI: t = 0.517

point travel を上げると、端点だけを見る discrete 検査はすぐ MISS に落ちます。両端点とも壁から遠いからです。一方 continuous 検査は赤い円で最初の接触時刻 $t^\*$ を正確に指し示します。

端点だけ検査するとなぜ漏れるのか

端点を二つだけ検査するとは、連続軌道を二点で標本化するということです。標本のあいだで起きたことは見えません。

解決の出発点は視点を変えることです。「今重なっているか」ではなく「 $t \in [0,1]$ のいつ重なるか」を問います。各頂点が1ステップのあいだ直線で動くと仮定すれば、位置は時間の1次関数になります。

\mathbf{x}_i(t) = (1-t)\,\mathbf{x}_i^{0} + t\,\mathbf{x}_i^{1}, \qquad t \in [0,1]

ここで $\mathbf{x}_i^{0}$ はステップ開始位置、 $\mathbf{x}_i^{1}$ は終了位置です。これで衝突は時間に関する求根問題になります。

共面条件 — 衝突を多項式の根へ

動く三角形どうしの最初の衝突は二つの場合しかありません。頂点が面に触れるか（vertex–face）、辺が別の辺に触れるか（edge–edge）。どちらの場合も、接触の瞬間に四点が同一平面上に並ぶという幾何的観察が核心です。

頂点 $\mathbf{p}$ と三角形 $(\mathbf{a},\mathbf{b},\mathbf{c})$ の共面条件は、体積（スカラー三重積）が0になることです。

f(t) = \big(\mathbf{p}(t)-\mathbf{a}(t)\big)\cdot \Big[\big(\mathbf{b}(t)-\mathbf{a}(t)\big)\times\big(\mathbf{c}(t)-\mathbf{a}(t)\big)\Big] = 0

各 $\mathbf{x}(t)$ は $t$ の1次式なので、 $f(t)$ は $t$ に関する3次多項式です。これを単変数（univariate）アプローチと呼びます。速いですが落とし穴があります。共面であることは衝突と同じではありません。根 $t^\*$ を求めた後、点が実際に三角形の内部にあるかを別途確認する必要があります。格子へ移すときに内外がひっくり返る件数は、STL をボクセル化する偶奇判定が別に数えています。一マス外すと、連続軌道が壁を抜けたのと同じ症状になります。

より頑健な道は多変数（multivariate）形式です。衝突をギャップベクトルの零点として直接書きます。

\mathbf{F}(t,u,v) = \mathbf{p}(t) - \big[\mathbf{a}(t) + u\,(\mathbf{b}(t)-\mathbf{a}(t)) + v\,(\mathbf{c}(t)-\mathbf{a}(t))\big]

$\mathbf{F}=\mathbf{0}$ かつ $u,v \ge 0$ 、 $u+v \le 1$ となる $(t,u,v)$ が存在すれば衝突です。 $u,v$ は面上の重心座標です。この形式には境界ケースがなく、純粋に代数的です。

数値的求根が見逃す衝突

3次式の根を浮動小数点で解くのが最も速いので、ほとんどのシミュレータがそうしています。しかし2012年以降の大規模ベンチマークが衝撃的な結論を出しました。単変数の数値ソルバの多くが false negative（見逃した衝突）を生むのです。

根を見逃す道は二つに分かれます。3次式の根は一般に無理数なので、浮動小数点では正確に書けません。それとは別に、二物体が同一平面上を平行に動くと $f(t)$ が恒等的に 0 になります。根が無限にあるこの場合を数値しきい値で落とそうとして、本物の衝突まで捨ててしまいます。

false positive（無い衝突を報告）は少し余分にコストがかかるだけで復旧可能です。しかし false negative はトンネリングへ直結します。衝突応答コードがラインサーチで貫通を防ぐ構造なら、一度の見逃しがシミュレーション全体を崩します。だからCCDに本当に求められる性質は速度ではなく**保守性（conservativeness）**です。見逃すくらいなら幻でも報告せよ、ということです。

包含関数と二分法 — 根を保守的に囲い込む

保守性を保証する最も古く、今なお最も明快な方法は Snyder（1992）の包含関数による二分法です。

核となる道具は包含関数（inclusion function）です。パラメータ箱 $B$ が与えられたとき、 $\mathbf{F}$ が $B$ 上で取りうる全ての値を包む軸平行境界箱 $\Box\mathbf{F}(B)$ を求めます。すると単純な事実が成り立ちます。 $\mathbf{0} \notin \Box\mathbf{F}(B)$ なら $B$ の中に根は決してありません。

幸い $\mathbf{F}$ は各パラメータについて1次（多重線形）です。多重線形関数は直方体上で極値を必ず頂点で取ります。したがって最もタイトな包含箱は、箱の四つ（または八つ）の頂点で $\mathbf{F}$ を評価した値の境界箱です。別途の区間演算は不要です。

アルゴリズムは分割統治です。

候補キューをパラメータ全領域 $[0,1]^d$ で初期化します。
キューから箱 $B$ を取り出し $\Box\mathbf{F}(B)$ を評価します。
$\mathbf{0} \notin \Box\mathbf{F}(B)$ なら捨てます（根なし保証）。
箱の幅が $\delta$ より小さければ根として採用します。
そうでなければ最も広い次元を半分に割り、両方をキューに入れます。
最初の接触時刻が欲しければキューを $t$ 基準で整列します。

パラメータ $(t,s)$ で、包含箱が原点を含まないマスは捨てます。幅が $\delta$ より小さくなったマスだけを根として残します。

Python — 動く点と線分の最初の接触#

2Dで動く点と動く線分の最初の接触時刻を、包含関数による二分法で実装します。ギャップ関数は $\mathbf{F}(t,s) = \mathbf{p}(t) - [\mathbf{a}(t) + s(\mathbf{b}(t)-\mathbf{a}(t))]$ で、根は $(t,s)\in[0,1]^2$ で探します。

import heapq
import numpy as np
 
def gap_F(t, s, query):
    """F(t,s) = P(t) - [A(t) + s (B(t)-A(t))]、直線軌道（t in [0,1]）。"""
    P0, P1, A0, A1, B0, B1 = query
    P = (1 - t) * P0 + t * P1
    A = (1 - t) * A0 + t * A1
    B = (1 - t) * B0 + t * B1
    return P - (A + s * (B - A))
 
def inclusion_aabb(box, query):
    """Fの最もタイトな軸平行境界：四頂点の値のhull（多重線形なので厳密）。"""
    t0, t1, s0, s1 = box
    corners = [gap_F(t0, s0, query), gap_F(t0, s1, query),
               gap_F(t1, s0, query), gap_F(t1, s1, query)]
    lo = np.min(corners, axis=0)
    hi = np.max(corners, axis=0)
    return lo, hi
 
def origin_inside(lo, hi):
    return bool(np.all(lo <= 0) and np.all(hi >= 0))
 
def inclusion_bisection_toi(query, delta=1e-4, max_checks=200000):
    """[0,1]での最初の接触時刻、無ければNone。Snyder/Redonの包含二分法。"""
    # 優先度キューをt0基準で整列 → 最初に収束した箱が最初の接触
    pq = [(0.0, 0.0, 1.0, 0.0, 1.0)]   # (t0, t0, t1, s0, s1)
    checks = 0
    while pq and checks < max_checks:
        _, t0, t1, s0, s1 = heapq.heappop(pq)
        checks += 1
        lo, hi = inclusion_aabb((t0, t1, s0, s1), query)
        if not origin_inside(lo, hi):
            continue                              # この箱に衝突なし
        if max(t1 - t0, s1 - s0) < delta:
            return t0, checks                     # 収束：最初の接触時刻
        if (t1 - t0) >= (s1 - s0):                # 広い次元を半分に
            tm = 0.5 * (t0 + t1)
            heapq.heappush(pq, (t0, t0, tm, s0, s1))
            heapq.heappush(pq, (tm, tm, t1, s0, s1))
        else:
            sm = 0.5 * (s0 + s1)
            heapq.heappush(pq, (t0, t0, t1, s0, sm))
            heapq.heappush(pq, (t0, t0, t1, sm, s1))
    return None, checks
 
# 1フレームで壁を横切る弾丸。端点だけ見ると両方とも壁から遠い。
query = (np.array([0.5, 3.5]), np.array([9.5, 3.5]),   # 点 P0 -> P1
         np.array([5.0, 1.0]), np.array([5.0, 1.0]),   # 壁端 A（静止）
         np.array([5.0, 6.0]), np.array([5.0, 6.0]))   # 壁端 B（静止）
 
toi, checks = inclusion_bisection_toi(query, delta=1e-4)
print(f"time-of-impact t* = {toi:.5f}   (検査した箱: {checks})")

出力：

time-of-impact t* = 0.50000   (検査した箱: 53)

step は分割を1レベル進め、 $\delta$ は赤いマスが止まってよい幅です。

point travel: 8.0wall drift: 1.5tolerance δ: 0.060

level: 0
active boxes: 1
earliest TOI ≈ 0.000

赤いマスが $\delta$ より狭くなった根領域です。 $\delta$ を上げると 53 回より早く止まり、マスは大きくなります。出力の $t^{*}=0.50000$ は、その幅で左端 $t_0$ が五桁に乗った値です。

解析解と正確に一致します。点は $x=5$ に $t=0.5$ で到達し、そのとき $y=3.5$ は線分の範囲 $[1,6]$ の中にあります。端点だけ見る検査が MISS と答えたまさにその衝突を、53回の箱評価で囲い込みました。保守的でありながら安価です。

衝突位置は割り算一行です。 $x$ は 0.5 から 9.5 へ動くので距離は 9、壁 $x=5$ までは 4.5 です。 $t=(5-0.5)/9=0.5$ が出力の 0.50000 と一致します。 $y=3.5$ は固定で、線分 $[1,6]$ の両端から 2.5 ずつ離れています。 $\delta=10^{-4}$ まで一軸だけ狭めると $2^{13}=8192$ 、 $2^{14}=16384$ なので 14 回の半分割です。検査数は 53 です。 $14\times 14=196$ の総当たりより少ないです。頂点での包含判定が空の箱を捨てるからです。返す値は箱の左端 $t_0$ です。印刷が 0.50000 なのは、その $t_0$ が五桁の中で 0.5 と一致するということです。

δ ひとつで精度と速度を交換する

包含二分法の最も実用的な美点は、ユーザのつまみがただ一つ、 $\delta$ である点です。

$\delta$ を上げると箱の分割が減り早く止まります。代わりに根領域が大きくなり false positive が増えます。下げると精密になりますが検査回数が増えます。リアルタイムシミュレーションでは検査回数の上限 $m_I$ を併せて置き、最悪実行時間を固定します。上限に達したら最も最近記録した衝突区間を保守的に返します — それでも見逃しません。

ここに最小分離（minimum separation）を加えると、製造公差や丸めによる貫通まで防げます。距離をユークリッドの代わりにチェビシェフ（L∞）に変える小さなトリックが問題を単純化します。 $\mathbf{F}=\mathbf{0}$ の代わりに $\lVert\mathbf{F}\rVert_\infty \le d$ を解けばよく、これは原点の代わりに辺の長さ $2d$ の立方体が包含箱と重なるかを見ることと同じです。コードでは箱を $d$ だけ大きくする一行で済みます。

コードに移す前に押さえること

常に保守的な向きに丸める。 頂点での $\mathbf{F}$ 評価に前進誤差限界を足して包含箱を少し大きくします。この一行が false negative を防ぎます。
単変数3次式の誘惑を警戒する。 速いですが共面平行運動で無限の根に引っかかります。保守性が必要なら多変数＋二分法を使います。
広域フェーズ（broad phase）を先に敷く。 空間ハッシュやAABB木で候補ペアをふるってからCCDを呼びます。CCDは狭域フェーズの高価な最後の関門です。
$\delta$ と $m_I$ をシーンスケールで正規化する。 絶対長ではなくドメインの大きさに合わせれば、速度やスケールが変わっても動きます。

δ を上げるとマスは膨らみ、見逃しは残らない

$\delta=10^{-4}$ は、一軸を 14 回割れという要求です。その幅で出力は $t^{*}=0.50000$ 、検査 53 回です。 $\delta$ を上げると赤いマスは大きくなり false positive は増えますが、原点を包含箱の外へ押し出しはしません。見逃すのは幅ではなく、共面の 3 次式を数値しきい値で落とすときです。

端点二つだけ見ると、この弾丸は MISS です。 $t=0.5$ で $x=5$ 、 $y=3.5$ が線分 $[1,6]$ の中にあるのに、両スナップショットは壁から遠いです。流体が格子の上にあり固体がセルのあいだを通るときは、Peskin の埋め込み境界が接触を格子点に打ち込みません。フレームのあいだに境界がセルを横切ると、端点検査はまた MISS を返します。