Pheemaの学習帳

軽率なアウトプット

レイと三角形の交差判定

この記事は レイトレアドベントカレンダー2018 20日目の記事です。

導入

初めて自作レンダラーを書いたとき、初めて画面に球が表示された感動を今でも忘れることができません。 皆さんも初めてレンダリングしたのは球ではありませんでしたか?[要出典]
しかし、現実は非情です。悲しいことに世の中の多くの3Dモデルは球ではなく三角形ポリゴンの集合として表現されます。したがって3Dモデルをレンダリングするためには、レイと三角形の交差判定は避けては通れない道というわけです。
そこで今回は、レイと三角形の交差判定をおこなう手法のひとつである Möller-Trumbore intersection algorithm [Möller97] について調べてみようと思います。また、各OSSレンダラーやライブラリではどのような交差判定手法を採用しているかについてもオマケ程度にまとめてみることにします。

Möller-Trumbore intersection algorithm

日本語の記事では Tomas Möllerのアルゴリズム と書かれていることが多いようです。しかし、海外の記事等では Möller-Trumbore intersection algorithm と書かれることが多いようなので、本記事ではそちらの表記に合わせることとします。

f:id:Pheema:20181218203846p:plain:w300
図1: 交差するレイと三角形

レイの表現

 \boldsymbol{o} から向き  \boldsymbol{d} に向かって長さ  t  \geq  0 だけ進んだレイ上の点  \boldsymbol{R} (t) は次のように表現することができます。

\begin{align} \boldsymbol{R} (t) = \boldsymbol{o} + t \boldsymbol{d} \label{eq:ray} \end{align}

ただし、向きを示すベクトル  \boldsymbol{d} は単位ベクトルとします。

三角形の内部に存在する点の表現

 \boldsymbol{v}_0, \boldsymbol{v}_1, \boldsymbol{v}_2 を頂点に持つ三角形を考えます。 このとき、三角形の内部に存在する任意の点  \boldsymbol{T} (u, v)

\begin{align} u \geq 0, v \geq 0, u + v \leq 1 \label{eq:uv} \end{align}

を満たす実数  u, v を用いて、次のように表現することができます。このことは「重心座標系」について調べるとより理解が深まるかと思います。

\begin{align} \boldsymbol{T} (u, v) = (1 - u - v) \boldsymbol{v}_0 + u \boldsymbol{v}_1 + v \boldsymbol{v}_2 \label{eq:tri} \end{align}

交差判定の式を導く

レイと三角形の交差判定をするためには、式  (\ref{eq:ray}) と式  (\ref{eq:tri}) から

\begin{align} \boldsymbol{R} (t) = \boldsymbol{T} (u, v) \label{eq:r_and_t} \end{align}

を変数  (t, u, v) について解けばよいです。 解いた結果が  t \geq 0 と式  (\ref{eq:uv}) を満たしているならば、レイと三角形は交差するといえます。

解くべき方程式を得るために、式  (\ref{eq:r_and_t}) を変形していきます。

\begin{align} &\boldsymbol{R} (t) = \boldsymbol{T} (u, v) \notag \\ ~ \notag \\ \Rightarrow& \boldsymbol{o} + t \boldsymbol{d} = (1 - u - v) \boldsymbol{v}_0 + u \boldsymbol{v}_1 + v \boldsymbol{v}_2 \notag \\ ~ \notag \\ \Rightarrow& \begin{bmatrix} -\boldsymbol{d} & \boldsymbol{v}_1 - \boldsymbol{v}_0 & \boldsymbol{v}_2 - \boldsymbol{v}_0 \end{bmatrix} \begin{bmatrix} t \\ u \\ v \end{bmatrix} = \boldsymbol{o} - \boldsymbol{v}_0 \notag \\ \Rightarrow& \begin{bmatrix} -\boldsymbol{d} & \boldsymbol{e}_1 & \boldsymbol{e}_2 \end{bmatrix} \begin{bmatrix} t \\ u \\ v \end{bmatrix} = \boldsymbol{r} \label{eq:intersection} \end{align}

式の見通しを良くするために三角形の辺を示すベクトル  \boldsymbol{e}_1 \equiv \boldsymbol{v}_1 - \boldsymbol{v}_0, \boldsymbol{e}_2 \equiv \boldsymbol{v}_2 - \boldsymbol{v}_0 と、右辺をまとめたベクトル  \boldsymbol{r} \equiv \boldsymbol{o} - \boldsymbol{v}_0 を定義したことに注意してください。

なお、列ベクトル  (t, u, v) の左にかかっている行列は列ベクトルを3つ並べてつくった  3 \times 3 正方行列です。 行列の各要素は次のように表されます。

\begin{align*} \begin{bmatrix} -\boldsymbol{d} & \boldsymbol{e}_1 & \boldsymbol{e}_0 \end{bmatrix} = \begin{bmatrix} -d_x & e_{1, x} & e_{0, x} \\ -d_y & e_{1, y} & e_{0, y} \\ -d_z & e_{1, z} & e_{0, z} \end{bmatrix} \end{align*}

以上の結果から、レイと三角形の交差判定は式  (\ref{eq:intersection}) (t, u, v) について解く問題へ帰着します。(3元1次の連立方程式ですね!)

交差判定の式を解く

 (\ref{eq:intersection})の解は、クラメルの公式を用いることで次のように書き下すことができます。

\begin{align} \begin{bmatrix} t \\ u \\ v \end{bmatrix} &= \frac{1}{ \begin{vmatrix} -\boldsymbol{d} & \boldsymbol{e}_1 & \boldsymbol{e}_2 \end{vmatrix} } \begin{bmatrix} \begin{vmatrix} \boldsymbol{r} & \boldsymbol{e}_1 & \boldsymbol{e}_2 \end{vmatrix} \\ \begin{vmatrix} -\boldsymbol{d} & \boldsymbol{r} & \boldsymbol{e}_2 \end{vmatrix} \\ \begin{vmatrix} -\boldsymbol{d} & \boldsymbol{e}_1 & \boldsymbol{r} \end{vmatrix} \end{bmatrix} \notag \\ ~ \notag \\ &= \frac{1}{ \left( \boldsymbol{d} \times \boldsymbol{e}_2 \right) \cdot \boldsymbol{e}_1} \begin{bmatrix} \left( \boldsymbol{r} \times \boldsymbol{e}_1 \right) \cdot \boldsymbol{e}_2 \\ \left( \boldsymbol{d} \times \boldsymbol{e}_2 \right) \cdot \boldsymbol{r} \\ \left( \boldsymbol{r} \times \boldsymbol{e}_1 \right) \cdot \boldsymbol{d} \label{eq:solved} \end{bmatrix} \end{align}

上記の式変形については、次の行列式の性質を利用しました。


  1. 行列式 i 列と  j 列を入れ替えると  -1 が掛かる
     \boldsymbol{b}_k (1 \leq k \leq n)行列式 k 列目の列ベクトルであるとします。このとき元の行列式 i 列目と  j を入れ替えた行列式の結果は、元の行列式の結果に  -1 を掛けたものとなります。すなわち、次の式が成り立ちます。 \begin{align*} \begin{vmatrix} \boldsymbol{b}_1 & \cdots & \boldsymbol{b}_i & \cdots & \boldsymbol{b}_j & \cdots & \boldsymbol{b}_n \end{vmatrix} \\ = -\begin{vmatrix} \boldsymbol{b}_1 & \cdots & \boldsymbol{b}_j & \cdots & \boldsymbol{b}_i & \cdots & \boldsymbol{b}_n \end{vmatrix} \end{align*}
  2.  3 \times 3 行列の行列式スカラー3重積で表すことができる
     3 \times 1 ベクトル  \boldsymbol{c}_0, \boldsymbol{c}_1, \boldsymbol{c}_2 について次の式が成り立ちます。 \begin{align*} \begin{vmatrix} \boldsymbol{c}_0 & \boldsymbol{c}_1 & \boldsymbol{c}_2 \end{vmatrix} = \boldsymbol{c}_0 \cdot \left( \boldsymbol{c}_1 \times \boldsymbol{c}_2 \right) \end{align*}

 (\ref{eq:solved}) において、 \boldsymbol{\alpha} \equiv \boldsymbol{d} \times \boldsymbol{e}_2, \boldsymbol{\beta} \equiv \boldsymbol{r} \times \boldsymbol{e}_1 とすれば、最終的に計算するべき式は次のように表すことができます。

\begin{align} \begin{bmatrix} t \\ u \\ v \end{bmatrix} = \frac{1}{\boldsymbol{\alpha} \cdot \boldsymbol{e}_1} \begin{bmatrix} \boldsymbol{\beta} \cdot \boldsymbol{e}_2 \\ \boldsymbol{\alpha} \cdot \boldsymbol{r} \\ \boldsymbol{\beta} \cdot \boldsymbol{d} \end{bmatrix} \end{align}

導出途中では行列や行列式をゴニョゴニョしていましたが、最終的にはベクトルの内積外積のみのシンプルな式になりました。 (t, u, v) さえ求めてしまえば、式  (\ref{eq:ray}) や 式  (\ref{eq:tri}) から衝突した点の座標も求めることができます。

実装について

コードが長くなっても見づらいと思うので、今回は レイが三角形の表側・裏側のどちらからでも交差する 場合のみを扱います。元論文 [Möller97] では、三角形の裏面からの交差を無視する場合の実装についても触れられています。もし必要であればそちらを参照してください。

必要のない計算を避けるために

  • 交差しないことが分かったらその先の交差判定処理を行わない(早期リターン)
  • すぐに必要とならない計算はできるだけ関数の後ろに持っていく

といった工夫が高速化につながるかと思います。

以下の実装例は [Möller97] の論文に基づいたものです。記事上ではconst参照等は省いていたり、なんとなくstd::optionalを使ってみたりしていますが、実際に動かす際は各自のコードに合う形に直してください。

実装例を見る

// Dot(Vector3, Vector3): Vector3同士の内積
// Cross(Vector3, Vector3): Vector3同士の外積

// レイと三角形の交差判定
//   交差する場合:   Vector3(t, u, v) 
//   交差しない場合: std::nullopt 
// が戻り値となる。
std::optional<Vector3> Intersect(
    Vector3 o,  // レイの原点
    Vector3 d,  // レイの向き
    Vector3 v0, // 三角形の頂点0
    Vector3 v1, // 三角形の頂点1
    Vector3 v2, // 三角形の頂点2
) {
    // 微小な定数([Möller97] での値)
    constexpr float kEpsilon = 1e-6f;

    Vector3 e1 = v1 - v0;
    Vector3 e2 = v2 - v0;

    Vector3 alpha = Cross(d, e2);
    float det = Dot(e1, alpha);

    // 三角形に対して、レイが平行に入射するような場合 det = 0 となる。
    // det が小さすぎると 1/det が大きくなりすぎて数値的に不安定になるので
    // det ≈ 0 の場合は交差しないこととする。
    if (-kEpsilon < det && det < kEpsilon) {
        return std::nullopt;
    }

    float invDet = 1.0f / det;
    Vector3 r = o - v0;
    
    // u が 0 <= u <= 1 を満たしているかを調べる。
    float u = Dot(alpha, r) * invDet;
    if (u < 0.0f || u > 1.0f) {
        return std::nullopt;
    }

    Vector3 beta = Cross(r, e1);

    // v が 0 <= v <= 1 かつ u + v <= 1 を満たすことを調べる。
    // すなわち、v が 0 <= v <= 1 - u をみたしているかを調べればOK。
    float v = Dot(d, beta) * invDet;
    if (v < 0.0f || u + v > 1.0f) {
        return std::nullopt;
    }

    // t が 0 <= t を満たすことを調べる。
    float t = Dot(e2, beta) * invDet;
    if (t < 0.0f) {
        return std::nullopt;
    }

    // 交差した!!!!
    return Vector3(t, u, v);
}

各レンダラー、ライブラリにおける交差判定手法

ここまでは Möller-Trumbore intersection algorithm について扱ってきました。では、世の中に公開されているレンダラーやライブラリではどのような手法が採用されているのでしょうか。個人的になんとなく気になったリポジトリ達を探ってみました。

※下記内容に関してはあくまでも 2018/12/20時点 での話です。

Cycles

みんな大好きBlenderに付属している物理ベースレンダラーです。Cyclesのソースコードはgitで管理されているので ここ を参考にCloneしてくるなり、Web から見てもよいかと思います。 現在、交差判定が書かれているファイルはBlenderリポジトリ内の /intern/cycles/util/util_math_intersect.h のようです。

この実装についてのコミットメッセージはこちらから見ることができます。このコミットメッセージによると元々 [Woop13] の手法を使っていたようですが、条件によっては計算精度の問題でアーティファクトが出てしまうことがあり、Embreeの実装をベースとしたPlücker 座標を用いた手法に変更したようです。

Embree

Intelが開発しているレイトレ用のOSSライブラリです。レイトレーシングに使用するような空間構造や交差判定などの処理が一通り揃っています。SIMDを使ってごりごりとチューニングしている印象です。

交差判定の処理が書かれているファイルは /embree/tree/master/kernels/geometry 下にまとまっています。

ファイル 手法
triangle_intersector_moeller.h [Möller97] の手法
triangle_intersector_pluecker.h Plücker 座標を用いた手法(元論文あるのか不明)
triangle_intersector_woop.h [Woop13] の手法

LuxCoreRender

自分が昔Blenderで使ってみたときは"LuxRender"という名前だった気がするのですが、いつの間にか名前が変わってますね。 ソースコードGithub から見ることができます。交差判定の処理が書かれている箇所は /include/luxrays/core/geometry/triangle.h で、手法としては [Möller97] のようですね。

LuxCore側のBVHを使用する場合はこの Intersect 関数が使用されるようですが、Embree側のBVHを利用する場合はそちらの実装が使われるようです。

Mitsuba

研究向けのレンダラーです。論文等で「新しいレンダリングの手法を提案したのでMitsubaに組み込んでみたよ!」みたいな使われ方をしている印象があります。交差判定の処理が書かれているファイルは /include/mitsuba/core/triangle.h で、手法としては [Möller97] のようです。

pbrt-v3

みんな大好き "Physically Based Rendering: From Theory To Implementation 3rd Edition" のレンダラーです。書籍の方については10月に オンライン版 が出たことで話題になりましたね。

こちらもソースコードGithub から見ることができます。交差判定の書かれているファイルは /src/shapes/triangle.cpp で、手法としては [Woop13] のようです。

PhysX

NVIDIA社が開発している物理エンジンで、最近(2018/12/03)オープンソース化するとの発表があり、いつの間にかGithub上にリポジトリが作られていました。レイトレからは話が逸れますが、物理エンジンでも衝突判定や可視判定で交差判定の処理をするはずなのでこちらも調べてみました。

交差判定の書かれているファイルは PhysX_3.4/Source/GeomUtils/src/intersection/GuIntersectionRayTriangle.h のようです。

実装はほぼ [Möller97] です。変更箇所としては、除算 1/det の計算が関数のほぼ最後に来ていますね。実装例では途中でこの除算を計算していますが、交差しない場合はこの除算の計算が無駄になるのでこの方が元の論文のコードよりも速そうな気もします。

まとめ

今回は Möller-Trumbore intersection algorithm [Möller97] の導出と実装例について調べてみました。実装も簡単で比較的高速な手法なようなので、初めてレンダラーを作るときなどには十分な手法かと思います。また、あまり精度を気にしないようなゲーム上でのレイキャスト処理などにも使用できそうです。

たまに趣味で取り組んでいるレンダラーはよくよく見てみると [Möller97] っぽい交差判定になっていたのですが、文献を調査してみるといかに雰囲気で計算してかを実感しました…。今回の調査を踏まえて実装しなおして、精度に問題があるようであれば別の手法も考慮に入れていきたいと思います。

なお、今回の記事では他手法との比較については 自分の気力の 紙面の都合上議論をしませんでした。計算誤差を考慮した手法などもあるようなので、また機会があれば調べてみようと思います!

参考文献

  • [Möller97] Tomas Möller, Ben Trumbore. 1997. Fast, Minimum Storage Ray/Triangle Intersection. Journal of Graphics Tools (JGT) 2, 1, 21-28
  • [Woop13] Sven Woop, Carsten Benthin, Ingo Wald. 2013. Watertight Ray/Triangle Intersection. Journal of Computer Graphics Techniques (JCGT) 2, 1, 65-82.