Quad
四角形との交差を考えていく.
平面はあくまで無限に伸びていく感じだったけど,今回のは大きさがちゃんとある.
以下のような感じで有限.
起点となるQ Q Q があり、v ⃗ \vec{v} v とu ⃗ \vec{u} u のベクトルで平面を張る感じ.
このベクトルが無限遠だったらPlaneと同じなんだけど、有限なのでそこがちょっと面倒ポイント.
さて、そしたらまずは平面の方程式を考えてみる.
こいつは以下のような形だった.
A x + B y + C z = D
\begin{equation}
\begin{split}
Ax+By+Cz=D
\end{split}
\end{equation}
A x + B y + C z = D
これをn ⃗ = ( A , B , C ) , w ⃗ = ( x , y , z ) \vec{n}=(A,B,C),\vec{w}=(x,y,z) n = ( A , B , C ) , w = ( x , y , z ) とすると、以下のように書き換えられる.
n ⃗ ⋅ w ⃗ = D
\begin{equation}
\begin{split}
\vec{n} \cdot \vec{w}=D
\end{split}
\end{equation}
n ⋅ w = D
このn ⃗ \vec{n} n は法線ベクトルとなる.
平面に対して垂直となるのが大事.
さて、そしたらいつも通り光線r ⃗ = o ⃗ + t d ⃗ \vec{r} = \vec{o} +t \vec{d} r = o + t d を代入する.
これはr ⃗ = w ⃗ \vec{r} = \vec{w} r = w となる位置に代入すればよいので、以下のようになる.
n ⃗ ⋅ ( o ⃗ + t d ⃗ ) = D
\begin{equation}
\begin{split}
\vec{n} \cdot (\vec{o}+t\vec{d})=D
\end{split}
\end{equation}
n ⋅ ( o + t d ) = D
さて、変換をしていこう.
n ⃗ ⋅ ( o ⃗ + t d ⃗ ) = D n ⃗ ⋅ o ⃗ + t n ⃗ ⋅ d ⃗ = D t n ⃗ ⋅ d ⃗ = D − n ⃗ ⋅ o ⃗ t = D − n ⃗ ⋅ o ⃗ n ⃗ ⋅ d ⃗
\begin{equation}
\begin{split}
& \vec{n} \cdot (\vec{o}+t\vec{d})=D \\
& \vec{n} \cdot \vec{o} + t \vec{n} \cdot \vec{d} = D \\
& t \vec{n} \cdot \vec{d} = D - \vec{n} \cdot \vec{o} \\
& t = \frac{D - \vec{n} \cdot \vec{o}}{\vec{n} \cdot \vec{d}}
\end{split}
\end{equation}
n ⋅ ( o + t d ) = D n ⋅ o + t n ⋅ d = D t n ⋅ d = D − n ⋅ o t = n ⋅ d D − n ⋅ o
次はこのt t t を計算するために、わかってないパラメータn ⃗ , D \vec{n},D n , D を考える.
まずはn ⃗ \vec{n} n を考え直す.
これは法線ということはu ⃗ , v ⃗ \vec{u},\vec{v} u , v で張る平面に垂直ということになるため、外積で表せる.
n ⃗ = u ⃗ × v ⃗
\begin{equation}
\begin{split}
\vec{n} = \vec{u} \times \vec{v}
\end{split}
\end{equation}
n = u × v
これで法線はOK.
次にD D D 、これは(2)を解けばよい.
(2)でわからないのはw ⃗ \vec{w} w .
これはどこかの点であればいいんだけど、図の中でわかってるのは起点Q Q Q である.
つまりw ⃗ = Q \vec{w}=Q w = Q としてしまえば以下のように言える.
n ⃗ ⋅ Q = D
\begin{equation}
\begin{split}
\vec{n} \cdot Q = D
\end{split}
\end{equation}
n ⋅ Q = D
あとはこれをコードに落とせばOK.
まず変数はm_q,m_u,m_vの3つを用意.
法線を計算した後に正規化しておく.
あとはQ Q Q とn ⃗ \vec{n} n の内積でD D D を計算するだけ.
これで役者はそろった.
Quad ( const Vec3 & q , const Vec3 & u , const Vec3 & v )
: m_q ( q )
, m_u ( u )
, m_v ( v )
{
Vec3 n = Cross ( m_u , m_v );
m_normal = n . normalized ();
m_d = Dot ( m_normal , m_q );
}
あとはt t t を計算してあげればOK.
double denom = Dot ( m_normal , ray . direction . xyz ());
if ( Abs ( denom ) < EPSILON_VALUE ) { return false ; }
double t = ( m_d - Dot ( m_normal , ray . origin . xyz ())) / denom ;
if ( ! interval . Contains ( t )) { return false ; }
これで結果を見てみる.
見た目が平面と変わらない...
範囲指定をしていないので当たり前ではある.
平面上の矩形外はtをはじく必要があるが、現段階では省いてないのが問題.
次はこれをやっていく.
今回の四角形はQ Q Q ,Q + u ⃗ Q+\vec{u} Q + u ,Q + v ⃗ Q+\vec{v} Q + v ,Q + u ⃗ + v ⃗ Q+\vec{u}+\vec{v} Q + u + v の4点が張る面である.
ここで平面の点をP P P とすると,
P = Q + α u ⃗ + β v ⃗ 0 < = α < = 1 , 0 < = β < = 1
\begin{equation}
\begin{split}
& P = Q + \alpha \vec{u} + \beta \vec{v} \\
& 0 <= \alpha <=1, 0 <= \beta <=1
\end{split}
\end{equation}
P = Q + α u + β v 0 <= α <= 1 , 0 <= β <= 1
と表せる.
ここで新しくp ⃗ \vec{p} p を以下のように表してみる.
p ⃗ = P − Q = α u ⃗ + β v ⃗
\begin{equation}
\begin{split}
\vec{p} = P-Q = \alpha \vec{u} + \beta \vec{v}
\end{split}
\end{equation}
p = P − Q = α u + β v
そしてこのp ⃗ \vec{p} p とu ⃗ \vec{u} u の外積を取ると,
u ⃗ × p ⃗ = u ⃗ × ( α u ⃗ + β v ⃗ ) = α u ⃗ × u ⃗ + β u ⃗ × v ⃗ = β u ⃗ × v ⃗
\begin{equation}
\begin{split}
& \vec{u} \times \vec{p} \\
& = \vec{u} \times (\alpha \vec{u} + \beta \vec{v}) \\
& = \alpha \vec{u} \times \vec{u} + \beta \vec{u} \times \vec{v} \\
& = \beta \vec{u} \times \vec{v}
\end{split}
\end{equation}
u × p = u × ( α u + β v ) = α u × u + β u × v = β u × v
となる.
同様にv ⃗ \vec{v} v に関しても計算すると,
v ⃗ × p ⃗ = v ⃗ × ( α u ⃗ + β v ⃗ ) = α v ⃗ × u ⃗ + β v ⃗ × v ⃗ = α v ⃗ × u ⃗
\begin{equation}
\begin{split}
& \vec{v} \times \vec{p} \\
& = \vec{v} \times (\alpha \vec{u} + \beta \vec{v}) \\
& = \alpha \vec{v} \times \vec{u} + \beta \vec{v} \times \vec{v} \\
& = \alpha \vec{v} \times \vec{u}
\end{split}
\end{equation}
v × p = v × ( α u + β v ) = α v × u + β v × v = α v × u
さて、この後が妙な方法なんだけど、外積だと移項ができないので、内積にしてしまおうという考えでn ⃗ \vec{n} n を掛ける.
n ⃗ ⋅ ( v ⃗ × p ⃗ ) = α n ⃗ ⋅ ( v ⃗ × u ⃗ ) n ⃗ ⋅ ( u ⃗ × p ⃗ ) = β n ⃗ ⋅ ( u ⃗ × v ⃗ )
\begin{equation}
\begin{split}
& \vec{n} \cdot (\vec{v} \times \vec{p}) = \alpha \vec{n} \cdot (\vec{v} \times \vec{u}) \\
& \vec{n} \cdot (\vec{u} \times \vec{p}) = \beta \vec{n} \cdot (\vec{u} \times \vec{v})
\end{split}
\end{equation}
n ⋅ ( v × p ) = α n ⋅ ( v × u ) n ⋅ ( u × p ) = β n ⋅ ( u × v )
内積にしたので塊で移行ができる.
α , β \alpha,\beta α , β はスカラーなのでこれもOKなはず.
α = n ⃗ ⋅ ( p ⃗ × v ⃗ ) n ⃗ ⋅ ( u ⃗ × v ⃗ ) β = n ⃗ ⋅ ( u ⃗ × p ⃗ ) n ⃗ ⋅ ( u ⃗ × v ⃗ )
\begin{equation}
\begin{split}
& \alpha = \frac{\vec{n} \cdot (\vec{p} \times \vec{v})}{\vec{n} \cdot (\vec{u} \times \vec{v})} \\
& \beta = \frac{\vec{n} \cdot (\vec{u} \times \vec{p})}{\vec{n} \cdot (\vec{u} \times \vec{v})}
\end{split}
\end{equation}
α = n ⋅ ( u × v ) n ⋅ ( p × v ) β = n ⋅ ( u × v ) n ⋅ ( u × p )
さて、ここで
w ⃗ = n ⃗ n ⃗ ⋅ ( u ⃗ × v ⃗ ) = n ⃗ n ⃗ ⋅ n ⃗
\begin{equation}
\begin{split}
\vec{w} = \frac{\vec{n}}{\vec{n} \cdot (\vec{u} \times \vec{v})} = \frac{\vec{n}}{\vec{n} \cdot \vec{n}}
\end{split}
\end{equation}
w = n ⋅ ( u × v ) n = n ⋅ n n
という定義をしておくと、次のように綺麗にまとまる.
α = w ⃗ ⋅ ( p ⃗ × v ⃗ ) β = w ⃗ ⋅ ( u ⃗ × p ⃗ )
\begin{equation}
\begin{split}
& \alpha = \vec{w} \cdot (\vec{p} \times \vec{v}) \\
& \beta = \vec{w} \cdot (\vec{u} \times \vec{p})
\end{split}
\end{equation}
α = w ⋅ ( p × v ) β = w ⋅ ( u × p )
これであとはα , β \alpha,\beta α , β が求まったので、これが範囲に収まってるかを調べるだけ!!
これをコードに落とし込もう!
まずはw ⃗ \vec{w} w を計算しておく.
これはコンストラクタで行っておく.
Vec3 n = Cross ( m_u , m_v );
m_normal = n . normalized ();
m_d = Dot ( m_normal , m_q );
m_w = n / Dot ( n , n );
次にα , β \alpha,\beta α , β を計算しておく.
Vec3 p = ray . point_at ( t );
Vec3 planarHitPointVector = p - m_q ;
double alpha = Dot ( m_w , Cross ( planarHitPointVector , m_v ));
double beta = Dot ( m_w , Cross ( m_u , planarHitPointVector ));
最後にこれが範囲内かを判定すればOK!!
// 範囲内?
Interval unit { 0.0 , 1.0 };
if ( ! unit . Contains ( alpha ) || ! unit . Contains ( beta )) { return false ; }
ここまで出来たら結果を見てみよう.
描画できてるね.
この平面が6個あればBoxを作ることができる.
RayTracing in One Weekendシリーズなんかはこれを6個並べてボックスを作って、並べることでコーネルボックスを構成してるね.
これだけで出来ちゃうのが面白いところ、といったところで今回はここまで.