ホーム＜ゲームつくろー！＜衝突判定編＜ 移動する２つの球の衝突場所と時刻を得る

3D衝突編
その９　移動する２つの球の衝突場所と時刻を得る

　ゲームなどのシミュレーションにおける球体（パーティクル）の衝突は、その移動が離散的であるため、非常に高速で動いている物同士の「すり抜け」問題が起こる事があります。これを防ぐには移動前と移動後の位置からなる「球体スウィープ」を用いて衝突を判定する必要があります。

　この章では、2つの移動している球の衝突位置とその時刻を得る方法を模索してみます。これがあれば、どれだけ高速で動いている球同士でも衝突を検知できます！

�@ 球体スウィープボリューム（SSV）

　球を移動させることで出来る薬のカプセルのような境界ボリュームの事を「球体スウィープボリューム（sphere-swept volume:SSV）」と言います。今回考えるパーティクルのような球がある点から別の点まで直線的に移動するとSSVができあがります。SSVの交差はそれほど難しくありませんが、その応用範囲はかなりに広いため非常に有用です。

�A パーティクルの位置を時間で表現する

　2つのパーティクルをA、Bとしましょう。その半径はrA、rBとします。また移動前のパーティクル位置をA(t=0)、移動後をA(t=1)と時間tを用いて表すとします：

　上図でパーティクルAはt=0からt=1まで等速直線運動で移動するとします。これはパーティクルBも同じです。そうして出来るSSV同士が接しているか否かは、ずばり上の緑色のベクトルの長さが2つのパーティクルの半径の合計よりも長いが短いかで決まります。緑色のベクトルは同時刻におけるパーティクルAとBの中心点を線分で結んだものです。このベクトルは時間tと共に変化しますので、時間の関数でその長さを表現できれば衝突時刻を算出できます。この章の中心はその関数の導出です。ちょっと面倒な点もありますが詳しくじっくり見ていきましょう。

　まず、ある時刻tにおける緑色のベクトルを、

と表わすことにします。これはAからBへ伸びるベクトルであるため、

となります。ここまでは何も難しいことは無いですね。A(t)及びB(t)は等速で直線的に移動すると考えていますので、次のように表わすことができます：

t=0で最初の位置、tが増えるほど前の位置から次の位置へ移動して、t=1で最終位置に到達します。これを上のC(t)に戻すと：

ここから緑色のベクトルC(t)はその初期位置C(0)とC(1)と時間tで表せることが分かりました。

　ベクトルC(t)の長さは、各成分を2乗して足しそのルートを取れば算出できます。ただルートは計算上扱いがやや面倒なので、長さの2乗で進めていきましょう：

　さて、この式を丁寧に展開していきます：

　ベクトルC(t)の長さのべき乗は最終的に時間tの2次関数になることがわかりました。P、Q、Rはそれぞれ、

です。C(0)もC(1)もパーティクルA、Bの初期位置及び到達位置から算出出来ますので、上のP,Q,Rはすべて定数になります。この式をベースとして衝突時刻を算出してみましょう。

�B パーティクルの衝突時刻・位置の算出

　パーティクルAとBが衝突する時、両者の距離はその半径の合計と一緒になります。�Aで時刻tから両者の距離（のべき乗）を算出する2次式を求めました。よって�Aの式の左辺がパーティクルAとBの半径の合計のべき乗になるtがあれば、それが衝突時刻になります。すなわち：

このような2次方程式になります。これをtについて解く（解の公式）と、

となります。

　解の公式はルートの中がマイナスになると実数解が存在しないことになります。これは「衝突しない」事を表わします。ルートの中がちょうど0になる場合、これはt=-Q/Pという時刻に両者が接することを表わします。そしてルートの中がプラスになる時、tは2つの解を持ちます。これは両者が接して、めり込んで、抜けるという状態がある事を表わします。抜ける時に接する瞬間がまたあるため答えが2つになるのです。P、Q、Rの値を求めておけば、上の式からたちどころに衝突時間tを算出でき、解の存在から衝突も判定出来ます。ただし一つ注意することがあります。解の公式で得られるtは0未満もしくは1以上になる事があります。それは過去に戻れば衝突する、もしくは未来のどこかで衝突する事を表しますが、今回の移動範囲内では衝突しない事も述べています。

　衝突位置（衝突接点）は衝突した時刻tでの緑色のベクトル上にあってパーティクルAから半径rAだけ進んだ位置になります。ただ衝突時のベクトルCの長さは双方のパーティクルの半径の和になりますから、その内分比を用いてもOKです：

�C P、Q、Rの意味合い

　�Aの式でまとめた定数P、Q、R。これらを表す式は意味合いがあります：

ベクトルDをC(1)-C(0)と置くと、PはベクトルDの長さのべき乗、QはC(0)とベクトルDの内積、そしてRはベクトルC(0)の長さのべき乗である事がわかります。ですからDを最初に求めておくと計算の見通しが良くなりますね。

�D 例外処理

　ここまでのお話にはいくつか例外が生じます。例えばP=0になる事があります。その時解の公式の分母が0になってしまうためtを求められません。P=0になるのはt=0からt=1まで双方のパーティクルが同じ速度で平行に移動した時です。速度が同じなのでいつまでも衝突する事がありません。もしくはt=0の時点で衝突していた場合はどのtでも衝突し続けます。よって、P=0の場合はt=0での双方のパーティクルの衝突を判断することになります。衝突している場合の衝突時刻はt=0です。問題は衝突位置で、めり込んでいる状態の場合はいつまでもめり込んだままで離れてくれません。この場合は仕方が無いのでパーティクルA(0)からB(0)へベクトルを伸ばし、半径rA分の長さの位置を衝突点としましょう。

　さらに問題なのはP=0（同速度の平行移動）でパーティクルの中心点が一緒になっている場合です。この場合あらゆる所が衝突位置になってしまいます。こうなるともうその位置の数学的な定義はできませんので、中心点を衝突点とするしかないかもしれません。

　例外その2はt=0の時にすでに衝突している状態です。この状態、案外良くあります（2つの球が接している状態でゲームがスタートするとか）。この場合衝突しているのは明らかですが、めり込んでいる場合の衝突位置がやはり問題になります。接している時の事を考慮するなら、パーティクルA(0)からB(0)へベクトルを伸ばし、半径rA分の長さの位置を衝突点とした方が無難かもしれません。

�D 関数公開

　では今回の理屈を踏まえた2つの移動するパーティクル（球）の衝突時刻と衝突位置を算出する関数を公開します。引数等で使用されているSphere等は3D衝突編その０「3D基本プリミティブの型定義」で定義されています：

2つの球の衝突までの時間と位置を取得

// 2つの球の衝突までの時間と位置を取得 // 双方の球は指定のベクトル方向に等速で進むと仮定 // stepSec : 双方のパーティクルが進む時間 // s1 : パーティクルA // v1 : パーティクルAの速度ベクトル // s2 : パーティクルB // v2 : パーティクルBの速度ベクトル // outSec : 衝突時刻 // outColPos : 衝突位置（接点） // outColPosS1(option) : 衝突時のパーティクルAの中心座標 // outColPosS2(option) : 衝突時のパーティクルBの中心座標 bool calcIntervalSphereSphere(     float stepSec,     const Sphere &s1,     const Vec3 &v1,     const Sphere &s2,     const Vec3 &v2,     float &outSec,     Point &outColPos,     Point *outColPosS1 = 0,     Point *outColPosS2 = 0 ) {     // 前位置及び到達位置におけるパーティクル間のベクトルを算出     Vec3 C0 = s2.p - s1.p;     Vec3 A1 = s1.p + v1 * stepSec;     Vec3 B1 = s2.p + v2 * stepSec;     Vec3 C1 = B1 - A1;     Vec3 D = C1 - C0;     float rAB = s1.r + s2.r;     float rABSq = rAB * rAB;     float P = D.lengthSq();     // 衝突判定に解の公式を使えるか？     if ( P == 0 ) {         // 平行移動 //         // t = 0で衝突しているか？         if ( ( s2.p - s1.p ).lengthSq() > rABSq ) {             return false;         }         outSec = 0.0f;         if ( outColPosS1 != 0 )             *outColPosS1 = s1.p;         if ( outColPosS2 != 0 )             *outColPosS2 = s2.p;         if ( s2.p == s1.p ) {             // 中心点も一緒なので中心点を衝突点として返す             outColPos = s1.p;             return true;         }         // 球A->Bのベクトル方向に長さrAの所を衝突点とする         outColPos = s1.p + ( s1.r / rAB ) * C0;         return true; // 同じ方向に移動     }     // 衝突検知可能 //     // 最初から衝突している？     if ( ( s2.p - s1.p ).lengthSq() <= rABSq ) {         outSec = 0.0f;         outColPos = s1.p + s1.r / rAB * C0;         if ( outColPosS1 != 0 )             *outColPosS1 = s1.p;         if ( outColPosS2 != 0 )             *outColPosS2 = s2.p;         return true;     }     float Q = C0.dot( D );     float R = C0.lengthSq();     // 衝突判定式     float judge = Q * Q - P * ( R - rAB * rAB );     if ( judge < 0 ) {         // 衝突していない         return false;     }     // 衝突時間の算出     float judge_rt = sqrtf( judge );     float t_plus = ( -Q + judge_rt ) / P;     float t_minus = ( -Q - judge_rt ) / P;     if ( t_minus > t_plus ) {         // t_minusを小さい方に         float tmp = t_minus;         t_minus = t_plus;         t_plus = tmp;     }     // 時間外衝突か？     if ( t_minus < 0.0f || t_minus > 1.0f ) {         return false;     }     // 衝突位置の決定     outSec = t_minus * stepSec;     Vec3 Atc = s1.p + v1 * stepSec * t_minus;     Vec3 Btc = s2.p + v2 * stepSec * t_minus;     outColPos = Atc + s1.r / rAB * ( Btc - Atc );     if ( outColPosS1 != 0 )     *outColPosS1 = Atc;     if ( outColPosS2 != 0 )     *outColPosS2 = Btc;     return true; // 衝突 }