ソフトウェアライスタライザ実装 その1 (頂点処理)およびソフトウェアライスタライザ実装 その2 (クリッピング)の続きです。
頂点処理およびクリッピングを終えた後は三角形メッシュが存在するピクセルを塗りつぶす(ラスタライズ)することでシーンの描画が完成します。
今回はラスタライズを実装してソフトウェアラスタライザを完成させました。
主に動画内で省略されていた投影を考慮した線形補間まわりの数学的証明を書いています。
(今回も三葉レイさんの動画を参照・引用しています。)
三角形の内外判定
頂点処理を終えた後、三角形の頂点位置が得られている状況で、各ピクセル(の中心)がその三角形の内部にあるかどうかを判定することが必要となる。
内外判定にはエッジ関数(外積ベクトルの$z$成分)を用いる。外積ベクトルを計算すると、下左図のようにある直線(AB)のどちら側に目的の点(P)が存在するのかどうかを、ベクトルの$z$成分の正負で判定することが出来る。これを下右図のように、三角形の頂点ABCに沿って一周する3つのベクトルについて行うことで、ある点が三角形内部にあるかどうかを判定出来る。全てのエッジ関数が正または全て負の場合のみ内部にあることになる。
内外判定を用いて描画するコードは次のようになる。効率化のためbounding_boxを導入し、色情報(明るさ$I$)は三角メッシュの法線ベクトル$\vec N$と平行光源$\vec d_L$の内積を用いた式
\begin{eqnarray}
I = 0.2 + 0.8*\mathrm{max}(0, \vec N \cdot \vec d_L)
\end{eqnarray}
に基づいて計算している。
// 平行光源
Vector3d light;
light << -1.0, 1.0, 0.0;
light = light/light.norm();
// ビューポート変換 (ラスタライズ)
float pixel_array[512][512][4] = {}; // (nx,ny,RGB色,z値)
for (int i=0;i<512;i++) {
for (int j=0;j<512;j++) {
pixel_array[i][j][3] = 2.0;
}
}
for (int t=0;t<T.size();t++) {
T[t].Bounding_box(); // Bounding_box計算
}
for (int i=0;i<nx;i++) {
for (int j=0;j<ny;j++) {
float x_coord = 1.0/nx + i*(2.0/nx) - 1.0;
float y_coord = 1.0/ny + j*(2.0/ny) - 1.0;
for (int t=0;t<T.size();t++) {
if (T[t].min_x <= x_coord and x_coord <= T[t].max_x) { // BB内ならラスタライズ
if (T[t].min_y <= y_coord and y_coord <= T[t].max_y) {
// Edge関数の計算
float edge_A, edge_B, edge_C;
edge_A = (T[t].C(0)-T[t].B(0))*(y_coord-T[t].B(1)) - (T[t].C(1)-T[t].B(1))*(x_coord-T[t].B(0)); // BC×BP
edge_B = (T[t].A(0)-T[t].C(0))*(y_coord-T[t].C(1)) - (T[t].A(1)-T[t].C(1))*(x_coord-T[t].C(0)); // CA×CP
edge_C = (T[t].B(0)-T[t].A(0))*(y_coord-T[t].A(1)) - (T[t].B(1)-T[t].A(1))*(x_coord-T[t].A(0)); // AB×AP
if (edge_A>=0 and edge_B>=0 and edge_C>=0) { // 三角形の内側(表側)であればラスタライズ
// z-バッファ法 (より手前にあるものを描画)
// ~
// 色情報を計算
// pixel_array[i][j][2] = (色情報);
}
}
}
}
}
}
出力結果は次のようになる。
三角メッシュの形が不自然に浮き彫りになってしまっているが、これは後述する補間の処理を用いることで修正出来る。
重心座標系での線形補間
三角形の内部にあると判定されたピクセルについては、その三角形が持つ色情報やデプス値などに基づいて塗りつぶしを行えばよい。ただし、色情報やデプス値は頂点のみに保存されているものなので、三角形内部に位置するようなピクセルについては3頂点の情報から適切に補間した値を使用する。今回はオーソドックスな線形補間を用いる場合を考える。
近い頂点の影響を良く受けるとして、その影響が線形であるならば、一次元系では上左図のようになる。何らかの量$C_A,C_B$は頂点からの距離の比$s,1-s$を用いて、次のように補間される。
\begin{eqnarray}
C = (1-s)C_A + sC_B
\end{eqnarray}
二次元系では、重心座標系を用いる。対向する三角形の面積の比$\lambda_A:\lambda_B:\lambda_C$を用いて、
\begin{eqnarray}
C = \lambda_A C_A + \lambda_B C_B + \lambda_C C_C
\end{eqnarray}
のように線形補間する。ここで、$\lambda_A,\lambda_B,\lambda_C$等の三角形の面積はエッジ関数(の絶対値)と同じものである。よって内外判定で計算したエッジ関数の値をそのまま再利用することで効率的に計算することが出来る。
if (edge_A>=0 and edge_B>=0 and edge_C>=0) { // 三角形の内側(表側)であればラスタライズ
float lambda_A, lambda_B, lambda_C;
float temp = edge_A + edge_B + edge_C;
lambda_A = edge_A/temp;
lambda_B = edge_B/temp;
lambda_C = edge_C/temp;
float depth = lambda_A*T[t].A(2) + lambda_B*T[t].B(2) + lambda_C*T[t].C(2); // depth
if (depth < pixel_array[i][j][3]) { // depth値が小さい時、
pixel_array[i][j][3] = depth;
for (int temp=0;temp<3;temp++) {
Vector3d N_temp;
N_temp = lambda_A*T[t].N_A + lambda_B*T[t].N_B + lambda_C*T[t].N_C; // 法線ベクトルの補間
N_temp /= N_temp.norm(); // 正規化
pixel_array[i][j][temp] = 0.2 + 0.8*max(0.0, light.dot(N_temp)); // 色の計算
}
}
}
出力画像は次のようになる。補間により三角メッシュの継ぎ目がわからない自然な描画ができている。
投影を考慮した線形補間
実は重心座標系での線形補間は正しくない。前回のクリッピング編でみたように、デバイス座標系は元々3次元空間であるクリッピング座標系を$\omega=1$平面に投影したものなので、奥行き方向($\omega$軸)の情報が失われている(下図)。
簡単のため、まず($x-\omega$)平面で考えると、本来正しい補間の式は
\begin{eqnarray}
C = (1-t)C_A + tC_B
\end{eqnarray}
となる。ここで、投影前の比$t:(1-t)$は次のように$s$を用いて書き直すことが出来る。 $s$の定義から、 ここで、$x$および$\omega$の定義: を代入して、 ----------------------------------------------------------(導出終了)----------------------------------------------------------
(導出過程)
\begin{eqnarray}
\frac{1-s}{s} &=& \frac{\|\frac{x}{\omega}-\frac{x_2}{\omega_2}\|}{\|\frac{x}{\omega}-\frac{x_1}{\omega_1}\|} \\
&=& \frac{\omega \omega_1}{\omega \omega_2} \frac{\|x\omega_2-x_2\omega\|}{\|x\omega_1-x_1\omega\|} \\
&=& \frac{\omega_1}{\omega_2} \frac{\|x\omega_2-x_2\omega\|}{\|x\omega_1-x_1\omega\|}
\end{eqnarray}
\begin{eqnarray}
x &=& tx_2 + (1-t)x_1 \\
\omega &=& t\omega_2 + (1-t)\omega_1 \\
\end{eqnarray}
\begin{eqnarray}
\frac{1-s}{s} &=& \frac{\omega_1}{\omega_2} \frac{\|(tx_2 + (1-t)x_1)\omega_2-x_2 (t\omega_2 + (1-t)\omega_1)\|}{\|(tx_2 + (1-t)x_1)\omega_1-x_1(t\omega_2 + (1-t)\omega_1)\|} \\
&=& \frac{\omega_1}{\omega_2} \frac{\| (1-t) (x_1\omega_2 - x_2\omega_1) \|}{\| t(x_2\omega_1 - x_1\omega_2) \|} \\
&=& \frac{\omega_1}{\omega_2} \frac{1-t}{t}
\end{eqnarray}
これを$t=$の形に書き直して、
\begin{eqnarray}
t = \frac{1}{\frac{\omega_2}{\omega_1}\frac{1-s}{s} + 1}
\end{eqnarray}
補間の式に代入すると、
\begin{eqnarray}
C &=& (1-t)C_A + tC_B \\
&=& C_A + t (C_B - C_A) \\
&=& \frac{C_A \frac{\omega_2}{\omega_1}\frac{1-s}{s}+ C_B}{\frac{\omega_2}{\omega_1}\frac{1-s}{s} + 1} \\
&=& \frac{\frac{C_A}{\omega_1}(1-s)+ \frac{C_B}{\omega_2} s}{\frac{1}{\omega_1}(1-s) + \frac{1}{\omega_2}s}
\end{eqnarray}
ここで、$\frac{1}{\omega} = \frac{1}{\omega_1}(1-s) + \frac{1}{\omega_2}s$の関係式が成立する ----------------------------------------------------------(導出終了)----------------------------------------------------------
(導出過程)
図より、
\begin{eqnarray}
\omega &&= \frac{(1-s)\frac{\omega_2}{\omega_1} \omega_1 + s \omega_2} {(1-s)\frac{\omega_2}{\omega_1} + s} \\
&=& \frac{\omega_2}{(1-s)\frac{\omega_2}{\omega_1} + s} \\
&=& \frac{1}{\frac{1}{\omega_1}(1-s)+\frac{1}{\omega_2}s}
\end{eqnarray}
\begin{eqnarray}
\frac{C}{\omega} = \frac{C_A}{\omega_1}(1-s)+ \frac{C_B}{\omega_2} s
\end{eqnarray}
のようになる。
つまり、$\omega_1,\omega_2$の値を保持しておけば、投影後の比$s:1-s$の値を用いて正しく補間を補正することが可能となる。(ただし今回補正する量は法線ベクトルであり、$\omega$の値は要らない(正規化するので大きさはどうでもよい)ので、投影変換の際に法線ベクトルも同時に$\omega_1,\omega_2$等の値で割る処理をしておけば、少し効率的に出来る。)
これは重心座標系での議論にもそのまま適応出来る。つまり、投影後の重心座標比$\lambda_A,\lambda_B,\lambda_C$を用いるが、その際補正する量を投影前の$\omega$値で割ってやればよい。
// -> デバイス座標系
for (int i=0;i<T.size();i++) {
for (int j=0;j<3;j++) {
T[i].A(j) = T[i].A(j)/T[i].A(3);
T[i].B(j) = T[i].B(j)/T[i].B(3);
T[i].C(j) = T[i].C(j)/T[i].C(3);
// 補間量をomegaで割っておく
T[i].N_A(j) /= T[i].A(3);
T[i].N_B(j) /= T[i].B(3);
T[i].N_C(j) /= T[i].C(3);
}
}
// ビューポート変換 (ラスタライズ)
// ~
出力画像は次のようになった。投影を考慮した補間のするしないでほぼ変わらない画像になった。うーん
Zバッファ法
ラスタライズでは、複数の三角形が重なる時、一番手前側にあるものを描画する。これは、各ピクセル毎にデプス値を保存しておいて、それより近い三角形が現れてた時のみ、ピクセル値を計算・デプス値を更新することで実現出来る。デプス値も三角形頂点の$z$値から線形補間することになるが、実は投影による補正を考える必要はない(投影変換においてそのように$z$値を規格化しているため)。
具体的に式で確かめてみる。クリップ座標系では、先程証明したように、$t$の代わりに$s$で表記した次の補間式
\begin{eqnarray}
\frac{1}{\omega} = \frac{1}{\omega_A}(1-s) + \frac{1}{\omega_B}s
\end{eqnarray}
の関係が成り立つ。カメラ座標系で$z_e$軸($=-\omega$)に戻っても同様の関係は保たれる。つまり、
\begin{eqnarray}
\frac{1}{z_e} = \frac{1}{z_e^A}(1-s) + \frac{1}{z_e^B}s \tag{1}
\end{eqnarray}
である。デバイス座標系で補間を行う場合、補間式がこの$(1)$式と同値なものであれば良い。
デバイス座標系での$z$値を$z_{NDC}$として、変換行列は
\begin{eqnarray}
\begin{pmatrix}
\cdot \\ \cdot \\ z_e \\ 1
\end{pmatrix}
\rightarrow
\begin{pmatrix}
\cdot & 0 & 0 & 0 \\
0 & \cdot & 0 & 0 \\
0 & 0 & -a & -b \\
0 & 0 & -1 & 0
\end{pmatrix}
\begin{pmatrix}
\cdot \\ \cdot \\ z_e \\ 1
\end{pmatrix}
=
\begin{pmatrix}
\cdot \\ \cdot \\ -a z_e - b \\ -z_e
\end{pmatrix}
\rightarrow^{1/z_e}
\begin{pmatrix}
\cdot \\ \cdot \\ a + \frac{b}{z_e} \\ \cdot
\end{pmatrix}
\end{eqnarray}
であったので、$z_{NDC}$は、
\begin{eqnarray}
z_{NDC} = a + \frac{b}{z_e}
\end{eqnarray}
と表せる。ここで、$z_{NDC}$を比$s:1-s$でそのまま線形補間した式は、
\begin{eqnarray}
z_{NDC} = (1-s) z_{NDC}^A + s z_{NDC}^B
\end{eqnarray}
であるが、これを$z_e$で書き直すと、
\begin{eqnarray}
a + \frac{b}{z_e} &=& (1-s) \left(a + \frac{b}{z_e^A}\right) + s \left(a + \frac{b}{z_e^B}\right) \\
\Leftrightarrow \frac{1}{z_e} &=& \frac{1}{z_e^A}(1-s) + \frac{1}{z_e^B}s
\end{eqnarray}
となり、確かにカメラ座標系での補間式$(1)$と一致している。
Z-バッファ法により、複数のオブジェクトが重なるような場合も描画出来るようになった。以下はteapotと球が重なるようなシーンを設定した場合の出力画像である。
GUIアニメーション化
OpenGL+GLFWのフレームワーク上で実装した。1フレームに0.3秒弱の計算時間がかかっているのでカクカク...
ラスタライザ部分の512*512*メッシュ数の3重forループがボトルネックになっているのでもう少し高速化したい。
