物理ベースレンダリング2 #CG

qiitaでは数式をひとつの記事にたくさん入れると重くてムリだったので、２つに分けました。

前半はこちら

物理ベースレンダリング1

こちらを参照して自分の理解のためにまとめました。基本的に理論的な部分をまとめています。

物理ベースレンダラedupt解説 - http://www.slideshare.net/h013/edupt-kaisetsu-22852235
双方向パストレーシングレンダラedubpt解説 - http://www.slideshare.net/h013/edubpt-v100
Monte Carlo Ray Tracing - http://www.cs.rutgers.edu/~decarlo/readings/mcrt-sg03c.pdf
Bi-directional Reflectance Distribution Function (BRDF) - Ruigang Yang - http://www.vis.uky.edu/~ryang/Teaching/CS684-fall05/lectures/lec10-BRDF.pdf

実装はこちら。rustやってみたかったのでなんとなくrustで書きました。まだおそい。

パストレーシング

レンダリング方程式をモンテカルロ積分することが目的。モンテカルロ積分で再帰するたびに計算量が指数関数的に増加してしまう。それを防ぐためにパストレーシングでは反射の時に1方向のみサンプリングすることで解決する。

必要なものは

$PDF$ - モンテカルロ積分するときの確率分布関数
$\boldsymbol{d}$ - 単位半球面上のサンプル(上のPDFに準じたもの)
$p_{RR}$ - 再帰打ち切り用のロシアンルーレットの乱数と確率
$f_r$ - BRDF

これらを計算し、レンダリング方程式に代入してモンテカルロ法を適用する。

微小立体角に関して関数をモンテカルロ法を用いて積分する。
モンテカルロ積分を行うときには、何らかの確率密度関数$P_\sigma$に基づいたサンプルが必要
この時、$P_\sigma(\Omega) = 1$を満たせば良い。

$\Omega$ - 単位半球面
$f(\boldsymbol{\omega})$ - 単位半球面に定義された任意の関数

一様分布によるサンプリング

サンプルが一様分布の場合、確率密度関数は定数になる

P_\sigma(\boldsymbol{\omega}) = C

\int_\Omega P_\sigma(\boldsymbol{\omega}) {\rm d} \sigma(\boldsymbol{\omega}) = 1

これを満たせば良い。単位半球面の立体角は$2 \pi$

\int_{\Omega} P_\sigma(\boldsymbol{\omega}) {\rm d} \sigma(\boldsymbol{\omega}) = \int_{\Omega} C {\rm d} \sigma(\boldsymbol{\omega}) = 2 \pi C = 1

C = \frac{1}{2 \pi}

P_\sigma(\boldsymbol{\omega}) = \frac{1}{2 \pi}

確率密度関数

P_\sigma(\boldsymbol{\omega}) = PDF(\theta, \phi) = \frac{1}{2 \pi}

累積分布関数

\begin{array}{l}
F(\theta, \phi) &=& \int P_\sigma(\boldsymbol{\omega}) {\rm d} \boldsymbol{\omega} \\
&=& \frac{1}{2 \pi} \int^\phi_0 \int^\theta_0 {\rm sin} \theta {\rm d} \theta {\rm d} \phi \\
&=& \frac{\phi}{2 \pi} \int^\theta_0 {\rm sin} \theta {\rm d} \theta \\
&=& \frac{\phi}{2 \pi}(1 - {\rm cos} \theta)
\end{array}

パラメータを分離すると

F(\phi) = \frac{\phi}{2 \pi} \\
F(\theta) = 1 - {\rm cos} \theta

確率密度関数に基づいたサンプルを生成するため、累積分布関数の逆関数を考える

$u_1$, $u_2$ - $(0, 1]$の範囲の乱数

\phi_i = 2 \pi u_1 \\
\theta_i = {\rm arccos} (u_2)

これを用いて単位半球面上でのサンプリングを考え、サンプル点のベクトルを求める。

$\boldsymbol{d}$ - サンプル点
$\boldsymbol{u}, \boldsymbol{v}, \boldsymbol{w}$ - 正規直交基底

\boldsymbol{d} = \boldsymbol{d_u} + \boldsymbol{d_v} + \boldsymbol{d_w}

\begin{array}{l}
 \boldsymbol{d_u} &=& {\rm sin} \theta_i {\rm cos} \phi_i \boldsymbol{u} \\
&=& {\rm sin}({\rm arccos} (u_2)) {\rm cos}(2 \pi u_1) \boldsymbol{u} \\
&=& \sqrt{1- u_2^2} {\rm cos}(2 \pi u_1) \boldsymbol{u}
 \end{array}

\begin{array}{l}
 \boldsymbol{d_v} &=& {\rm sin} \theta_i {\rm sin} \phi_i \boldsymbol{v} \\
&=& {\rm sin}({\rm arccos} (u_2)) {\rm sin}(2 \pi u_1) \boldsymbol{v} \\
&=& \sqrt{1 - u_2^2} {\rm sin}(2 \pi u_1) \boldsymbol{v}
 \end{array}

\begin{array}{l}
 \boldsymbol{d_w} &=& {\rm cos} \theta_i \boldsymbol{w} \\
&=& {\rm cos}({\rm arccos} (u_2)) \boldsymbol{w} \\
&=& u_2 \boldsymbol{w}
 \end{array}

これらをレンタリング方程式に当てはめる。

L_o(\boldsymbol{x}, \boldsymbol{\omega}) = L_e(\boldsymbol{x}, \boldsymbol{\omega}) + \int_{\Omega_x} f_r(\boldsymbol{x}, \boldsymbol{\omega}, \boldsymbol{\omega^\prime}) L_i(\boldsymbol{x}, \boldsymbol{\omega^\prime}) {\rm cos} \theta {\rm d} \sigma(\boldsymbol{\omega^\prime})

レンダリング方程式にモンテカルロ積分を適用して、推定値を求める。このとき$\hat L_o$はサンプルする方向を1方向のみとすることで、再帰が深くなっても計算量は指数関数的に増加しない。この$\hat L_o$を何度も評価することで真値$L_o$に近づく。

L_o = \frac{1}{N} \sum_{i = 1}^N {\hat L_o}

今回、BRDF項($f_r$)は完全拡散面を考えるため

f_r = \frac{\rho}{\pi}

これらを代入すると

\begin{array}{l}
{\hat L_o(\boldsymbol{x}}, \boldsymbol{\omega}) &=& L_e(\boldsymbol{x}, \boldsymbol{\omega}) + \frac{f_r(\boldsymbol{x}, \boldsymbol{\omega}, \boldsymbol{\omega^\prime_i}) L_i(\boldsymbol{x}, \boldsymbol{\omega^\prime_i}) {\rm cos} \theta}{PDF(\boldsymbol{\omega^\prime_i})} \\
&=& L_e(\boldsymbol{x}, \boldsymbol{\omega}) + \frac{f_r(\boldsymbol{x}, \boldsymbol{\omega}, \boldsymbol{\omega^\prime_i}) L_o(\boldsymbol{r}(\boldsymbol{x}, \boldsymbol{\omega^\prime}), - \boldsymbol{\omega^\prime}) {\rm cos} \theta}{PDF(\boldsymbol{\omega^\prime_i})} \\
&=& L_e(\boldsymbol{x}, \boldsymbol{\omega}) + \frac{(\rho / \pi) L_o(\boldsymbol{r}(\boldsymbol{x}, \boldsymbol{\omega^\prime}), - \boldsymbol{\omega^\prime}) {\rm cos} \theta}{1/ 2 \pi} \\
&=& L_e(\boldsymbol{x}, \boldsymbol{\omega}) + 2 \rho L_o(\boldsymbol{r}(\boldsymbol{x}, \boldsymbol{\omega^\prime}), - \boldsymbol{\omega^\prime}) {\rm cos} \theta
\end{array}

再帰的な関数になっているため、実装ではロシアンルーレットで再帰の中断を行う。

{\hat L_o(\boldsymbol{x}}, \boldsymbol{\omega}) =
\begin{cases}
L_e(\boldsymbol{x}, \boldsymbol{\omega}) + 2 \rho L_o(\boldsymbol{r}(\boldsymbol{x}, \boldsymbol{\omega^\prime}), - \boldsymbol{\omega^\prime}) {\rm cos} \theta \frac{1}{p_{RR}}  & (r_{RR} \lt p_{RR}) \\
L_e(\boldsymbol{x}, \boldsymbol{\omega}) & (r_{RR} \geq p_{RR})
\end{cases}

余弦によるインポータンスサンプリング

レンダリング方程式のcos項に注目すると、真上からの寄与を大きく余弦にしたがって取るようにすると分散を抑えられる(ランバートの余弦則)。

${\rm cos} \theta$にしたがって重点的にサンプリングしたい場合、確率密度関数を次のようにして求める。

P_\sigma(\boldsymbol{\omega}) = C {\rm cos} \theta

\begin{array}{l}
\int_{\Omega} P_\sigma(\boldsymbol{\omega}) {\rm d} \sigma(\boldsymbol{\omega}) &=& \int_{\Omega} C {\rm cos} \theta {\rm d} \sigma(\boldsymbol{\omega}) \\
&=& C \int^{2 \pi}_0 \int^{\pi/2}_0 {\rm cos} \theta {\rm sin} \theta {\rm d} \theta {\rm d} \phi \\
&=& \pi C \int^{\pi / 2}_0 {\rm sin} 2 \theta {\rm d} \theta \\
&=& \pi C
\end{array}

C = \frac{1}{\pi}

P_\sigma(\boldsymbol{\omega}) = \frac{{\rm cos} \theta}{\pi}

確率密度関数

P_\sigma(\boldsymbol{\omega}) = PDF(\theta, \phi) = \frac{{\rm cos} \theta}{\pi}

累積分布関数

\begin{array}{l}
F(\theta, \phi) &=& \int P_\sigma(\boldsymbol{\omega}) {\rm d} \boldsymbol{\omega} \\
&=& \frac{1}{\pi} \int {\rm cos} \theta {\rm d} \boldsymbol{\omega} \\
&=& \frac{1}{\pi} \int^\phi_0 \int^\theta_0 {\rm cos} \theta {\rm sin} \theta {\rm d} \theta {\rm d} \phi \\
&=& \frac{\phi}{2 \pi} \int^\theta_0 {\rm sin} (2 \theta) {\rm d} \theta \\
&=& \frac{\phi}{4 \pi}(1 - {\rm cos} (2 \theta)) \\
&=& \frac{\phi}{2 \pi}(1 - {\rm cos^2} \theta)
\end{array}

パラメータを分離すると

F(\phi) = \frac{\phi}{2 \pi} \\
F(\theta) = 1 - {\rm cos^2} \theta

確率密度関数に基づいたサンプルを生成するため、累積分布関数の逆関数を考える

\phi_i = 2 \pi u_1 \\
\theta_i = {\rm arccos} (\sqrt{1 - u_2})

これを用いて単位半球面上でのサンプリングを考えるとベクトルの各成分は

\begin{array}{l}
 \boldsymbol{d_u} &=& {\rm sin} \theta_i {\rm cos} \phi_i \boldsymbol{u} \\
&=& {\rm sin}({\rm arccos} (\sqrt{1 - u_2})) {\rm cos}(2 \pi u_1) \boldsymbol{u} \\
&=& \sqrt{u_2} {\rm cos}(2 \pi u_1) \boldsymbol{u}
 \end{array}

\begin{array}{l}
 \boldsymbol{d_v} &=& {\rm sin} \theta_i {\rm sin} \phi_i \boldsymbol{v} \\
&=& {\rm sin}({\rm arccos} (\sqrt{1 - u_2})) {\rm sin}(2 \pi u_1) \boldsymbol{v} \\
&=& \sqrt{u_2} {\rm sin}(2 \pi u_1) \boldsymbol{v}
 \end{array}

\begin{array}{l}
 \boldsymbol{d_w} &=& {\rm cos} \theta_i \boldsymbol{w} \\
&=& {\rm cos}({\rm arccos} (\sqrt{1 - u_2})) \boldsymbol{w} \\
&=& \sqrt{1 - u_2} \boldsymbol{w}
 \end{array}

レンダリング方程式にモンテカルロ積分を適用して、推定値を求める。

L_o(\boldsymbol{x}, \boldsymbol{\omega}) = L_e(\boldsymbol{x}, \boldsymbol{\omega}) + \int_{\Omega_x} f_r(\boldsymbol{x}, \boldsymbol{\omega}, \boldsymbol{\omega^\prime}) L_i(\boldsymbol{x}, \boldsymbol{\omega^\prime}) {\rm cos} \theta {\rm d} \sigma(\boldsymbol{\omega^\prime})

\begin{array}{l}
{\hat L_o(\boldsymbol{x}}, \boldsymbol{\omega}) &=& L_e(\boldsymbol{x}, \boldsymbol{\omega}) + \frac{f_r(\boldsymbol{x}, \boldsymbol{\omega}, \boldsymbol{\omega^\prime_i}) L_i(\boldsymbol{x}, \boldsymbol{\omega^\prime_i}) {\rm cos} \theta}{PDF(\boldsymbol{\omega^\prime_i})} \\
&=& L_e(\boldsymbol{x}, \boldsymbol{\omega}) + \frac{f_r(\boldsymbol{x}, \boldsymbol{\omega}, \boldsymbol{\omega^\prime_i}) L_o(\boldsymbol{r}(\boldsymbol{x}, \boldsymbol{\omega^\prime}), - \boldsymbol{\omega^\prime}) {\rm cos} \theta}{PDF(\boldsymbol{\omega^\prime_i})} \\
&=& L_e(\boldsymbol{x}, \boldsymbol{\omega}) + \frac{(\rho / \pi) L_o(\boldsymbol{r}(\boldsymbol{x}, \boldsymbol{\omega^\prime}), - \boldsymbol{\omega^\prime}) {\rm cos} \theta}{{\rm cos} \theta / \pi} \\
&=& L_e(\boldsymbol{x}, \boldsymbol{\omega}) + \rho L_o(\boldsymbol{r}(\boldsymbol{x}, \boldsymbol{\omega^\prime}), - \boldsymbol{\omega^\prime})
\end{array}

{\hat L_o(\boldsymbol{x}}, \boldsymbol{\omega}) =
\begin{cases}
L_e(\boldsymbol{x}, \boldsymbol{\omega}) + \rho L_o(\boldsymbol{r}(\boldsymbol{x}, \boldsymbol{\omega^\prime}), - \boldsymbol{\omega^\prime}) \frac{1}{p_{RR}}  & (r_{RR} \lt p_{RR}) \\
L_e(\boldsymbol{x}, \boldsymbol{\omega}) & (r_{RR} \geq p_{RR})
\end{cases}

レンダリング結果

512x512で一様サンプルとcos重点的サンプルでレンダリングしてみました。左が一様サンプル、右がcos重点的サンプルです。

16 sample

64 sample

256 sample

1024 sample

サンプル数を多くするとだんだんと収束し、cos重点的サンプルの方が分散が小さいのがわかります。