ここ数年，レイトレ合宿で最底辺のあたりの順位をうろついているので，流石にちょっとヤバいと思い始めてきたので，本格的に勉強をしようかなと思いました。最初にレイトレで参考した本がPeter Shirleyの”Realstic Ray Tracing”だったのですが，ペーパーバッグ版を買ったのですが，印字などが恐ろしく汚く画像も黒塗り状態で全く分からない状態で，知りたいところがあまり良くわからないような状態でした。今回，ハードカバー版を入手したので，長年よくわからなかったところが，こういう画像だったのか！と知れるようになりました。読めるようになり俄然やる気が湧いてきたので，学習した内容をここにまとめておこうと思います。1章はベクトルクラスやRGBクラス，乱数など基本内容なので飛ばして，2章から始めることにします。
2章はレイとオブジェクトの交差判定について記述されています。
直線，線分，レイを総称して「レイ」と記述することにします。レイはパラメトリックな直線として記述することができ，以下の式で表します。
\begin{eqnarray} {\mathbf p}(t) = {\mathbf o} + t {\mathbf d} \tag{1} \end{eqnarray}
ここで，\({\mathbf o}\)は起点で，\({\mathbf d}\)は方向ベクトル，\(t\)は区間に対応するものになります。
サーフェイスは陰方程式(implicit equations)またはパラメトリック方程式(parametric equations)の2つのどちらかでよく記述されます。
陰方程式はサーフェイス上の点群を暗黙的に次のように定義します。
\begin{eqnarray} f(x, y, z) = 0 \tag{2} \end{eqnarray}
\((x, y, z)\)は任意の点を表す。これを\({\mathbf p}\)とすれば，
\begin{eqnarray} f({\mathbf p}) = 0 \tag{3} \end{eqnarray}
となる。ここで，式(1)を代入すれば
\begin{eqnarray} f({\mathbf p}(t)) &=& 0 \\ f({\mathbf o} + t {\mathbf d}) &=& 0 \tag{4} \end{eqnarray}
式(4)を満たすような例として，法線ベクトル\({\mathbf n}\)を持ち点\({\mathbf a}\)を通る無限平面について考えると，この平面は陰方程式として次のように記述できます。
\begin{eqnarray} ({\mathbf p} - {\mathbf a}) \cdot {\mathbf n} = 0 \tag{5} \end{eqnarray}
\({\mathbf a}\)と\({\mathbf n}\)は既知の値とすることに注意。この方程式は”平面法線に対して\({\mathbf a}\)から\({\mathbf p}\)へのベクトルが垂直である”ということを意味しています。\({\mathbf p}\)が平面でないときは，\({\mathbf p} - {\mathbf a}\)は\({\mathbf n}\)と直角になりません。 式(1)を式(5)に代入すれば，
\begin{eqnarray} ({\mathbf o} + t {\mathbf d} - {\mathbf a}) \cdot {\mathbf n} = 0 \tag{6} \end{eqnarray}
上式において，未知のものは\(t\)のみなので，\(t\)について解き
\begin{eqnarray} & & ({\mathbf o} + t {\mathbf d} - {\mathbf a}) \cdot {\mathbf n} = 0 \\ & \iff & {\mathbf o} \cdot {\mathbf n} + t {\mathbf d} \cdot {\mathbf n} - {\mathbf a} \cdot {\mathbf n} = 0 \\ & \iff& t {\mathbf d} \cdot {\mathbf n} = {\mathbf a} \cdot {\mathbf n} - {\mathbf o} \cdot {\mathbf n} \\ t &=& \frac{({\mathbf a} - {\mathbf o}) \cdot {\mathbf n} }{ {\mathbf d} \cdot {\mathbf n} } \tag{7} \end{eqnarray}
今，興味があるのは直線上のある区間で交差があるかということなので，\(t\)がこの範囲内であるかを判定すればよいことになります。分母がゼロとなるケースは平面に対して平行であることを表し，交差は定義されない。分子と分母がともにゼロになる場合は，レイは平面上にあること示します。
ライティングの計算のために，サーフェイス法線が必要となります。交差点\({\mathbf p}\)におけるサーフェイス法線は，陰関数の勾配で与えられます。つまり，
\begin{eqnarray} {\mathbf n} = {\nabla} f({\mathbf p}) = \left( \frac{\partial f({\mathbf p})}{\partial x}, \frac{\partial f({\mathbf p})}{\partial y}, \frac{\partial f({\mathbf p})}{\partial z} \right) \tag{8} \end{eqnarray}
となります。勾配ベクトルはサーフェイス内側あるいは外側になるかもしれないので注意してください。
３次元のサーフェイスを指定する別の方法は，2次元パラメータで表現する方法です。このサーフェイスは次の形式を持ちます。
\begin{eqnarray} x &=& f(u, v) \\ y &=& g(u, v) \tag{9} \\ z &=& h(u, v) \\ \end{eqnarray}
例えば，z-upとする原点を中心とする単位球を極座標系で表すと，次のパラメトリック方程式が得られます。
\begin{eqnarray} x &=& \cos \phi \sin \theta \\ y &=& \sin \phi \sin \theta \tag{10} \\ z &=& \cos \theta \end{eqnarray}
パラメトリックサーフェイスとレイが交差するためには，デカルト座標系ですべて一致するように以下のように設定します。
\begin{eqnarray} o_x + t d_x = f(u, v) \\ o_y + t d_y = g(u, v) \tag{11} \\ o_z + t d_z = h(u, v) \end{eqnarray}
上記を３つの未知数(\(t, u, v\))について数値的に解けばよいです。
パラメトリックサーフェイスについての法線ベクトルは次で与えられます。
\begin{eqnarray} {\mathbf n}(u, v) = \left( \frac{ \partial f }{ \partial u }, \frac{ \partial g }{ \partial u }, \frac{ \partial h }{ \partial u } \right) \times \left( \frac{ \partial f }{ \partial v}, \frac{ \partial g }{ \partial v }, \frac{ \partial h }{ \partial v } \right) \tag{12} \end{eqnarray}
半径\(R\),中心\( {\mathbf c} = (c_x, c_y, c_z) \)を持つ球を考えます。陰方程式で交差を表せば
\begin{eqnarray} (x - c_x)^2 + (y - c_y)^2 + (z - c_z)^2 - R^2 = 0 \tag{13} \end{eqnarray}
となり，これをベクトルで表現すれば
\begin{eqnarray} ({\mathbf p} - {\mathbf c}) \cdot ({\mathbf p} - {\mathbf c}) - R^2 = 0 \tag{14} \end{eqnarray}
となります。ここで，式(1)を代入して，２次方程式の解の公式を適用して\(t\)について解けば
\begin{eqnarray} t = \frac{ -{\mathbf d} \cdot ({\mathbf o} - {\mathbf c}) \pm \sqrt{ ({\mathbf d} \cdot ({\mathbf o} - {\mathbf c})^2 ) - ({\mathbf d} \cdot {\mathbf d}) ( ({\mathbf o} - {\mathbf c}) \cdot ({\mathbf o} - {\mathbf c}) - R^2 ) } }{ ({\mathbf d} \cdot {\mathbf d}) } \tag{15} \end{eqnarray}
が得られます。単位法線は\( ({\mathbf p} - {\mathbf c}) / R\)で求められます。
三角形は3点\( {\mathbf a} \)，\( {\mathbf b} \)，\( {\mathbf c} \)によって定義され，同一直線上にない場合に，これらは平面を定義します。この平面を記述する一般的な方法として重心座標系(barycentric coordinates)があります:
\begin{eqnarray} {\mathbf p}(\alpha, \beta, \gamma) = \alpha {\mathbf a} + \beta {\mathbf b} + \gamma {\mathbf c} \tag{16} \end{eqnarray}
\(\alpha\)，\(\beta\)，\(\gamma\)は次の制限があります。
\begin{eqnarray} \alpha + \beta + \gamma = 1 \tag{17} \end{eqnarray}
点\({\mathbf p}\)が三角形内部にあるためには，
\begin{eqnarray} 0 \lt \alpha \lt 1 \\ 0 \lt \beta \lt 1 \tag{18} \\ 0 \lt \gamma \lt 1 \end{eqnarray}
を満たす必要があります。変数を１つ消すために式(17)を\(\alpha = 1 - \beta - \gamma\)と変形し，式(16)に代入すると
\begin{eqnarray} {\mathbf p}(\beta, \gamma) = {\mathbf a} + \beta ({\mathbf b} - {\mathbf a}) + \gamma ( {\mathbf c} - {\mathbf a}) \tag{19} \end{eqnarray}
が得られます。点が三角形内にある条件は
\begin{eqnarray} \beta + \gamma & \lt & 1 \\ 0 & \lt & \beta \tag{20} \\ 0 & \lt & \gamma \end{eqnarray}
となります。式(19)に式(1)を代入し，成分表示すると
\begin{eqnarray} o_x + t d_x &=& a_x + \beta (b_x - a_x) + \gamma( c_x - a_x ) \\ o_y + t d_y &=& a_y + \beta (b_y - a_y) + \gamma( c_y - a_y ) \tag{21} \\ o_z + t d_z &=& a_z + \beta (b_z - a_z) + \gamma( c_z - a_z) \end{eqnarray}
が満たすべき式で，これらの式を用いて\(t\)，\(\beta\)，\(\gamma\)を求めればよいことになります。
ここで，式(21)を行列で書き直し，クラメルの公式を用いることで解が求まります。
行列\({\mathbf A}\)を
\begin{eqnarray} {\mathbf A} = \begin{bmatrix} a_x - b_x & a_x - c_x & d_x \\ a_y - b_y & a_y - c_y & d_y \\ a_z - b_z & a_z - c_z & d_z \end{bmatrix} \tag{22} \end{eqnarray}
とすると，\(\beta\)，\(\gamma\)，\(t\)は行列式\(| {\mathbf A} | \)を用いて
\begin{eqnarray} \beta &=& \frac{ \begin{vmatrix} a_x - o_x & a_x - c_x & d_x \\ a_y - o_y & a_y - c_y & d_y \\ a_z - o_z & a_z - c_z & d_z \end{vmatrix} }{ | {\mathbf A} |}, \tag{23} \\ \gamma &=& \frac{ \begin{vmatrix} a_x - b_x & a_x - o_x & d_x \\ a_y - b_y & a_y - o_y & d_y \\ a_z - b_z & a_z - o_z & d_z \end{vmatrix} }{ | {\mathbf A} | }, \tag{24} \\ t &=& \frac{ \begin{vmatrix} a_x - b_x & a_x - c_x & a_x - o_x \\ a_y - b_y & a_y - c_y & a_y - o_y \\ a_z - b_z & a_z - c_z & a_z - o_z \end{vmatrix}}{ | {\mathbf A} |} \tag{25} \end{eqnarray}
と求めることができます。
三角形の点が反時計回りに格納されている場合，正規化されていない法線ベクトルは
\begin{eqnarray} {\mathbf n} = ({\mathbf b} - {\mathbf a}) \times ({\mathbf c} - {\mathbf c}) \tag{26} \end{eqnarray}
で計算できます。
オブジェクトが１個以上ある場合は，レイに最初のヒットするオブジェクトを見つければよいです。オブジェクト指向型の言語の場合は抽象クラスやインタフェース，スーパークラスでヒットするかどうかを判定するhitルーチンを作っておき，これを継承して呼び出し，交差だったらtrueを返すようにし，\(t\)のようなあるデータを埋めるように実装します。疑似コードは次のようになります。
function surfacelist.hit(ray o + t * d, t0, t1, prim)
// prim はヒットオブジェクトへのポインタ
// t と prim は交差があった場合に格納されます
bool hitone ← false
t ← t1
for all surf in list do
if surf → hit(ray o + t * d, t0, t1, prim) then
hitone ← true
return hitone 実装コードをGithubにまとめておきました。
今回は2章を読みました。既知の内容が多かったので，特に躓くこともなく読めました。