こちらを参照して自分の理解のためにまとめました。基本的に理論的な部分をまとめています。
物理ベースレンダラ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で書きました。まだおそい。
放射束 (Flux)
定義: 単位時間あたり、ある領域を通過するフォトン数 (エネルギー)
\Phi \left[ \frac{\rm J}{\rm S} \right] = \Phi [{\rm W}]
放射照度 (Irradiance)
定義: 放射束の面積密度 (単位面積あたりの放射束)
E(x) \left[ \frac{\rm W}{\rm m^2} \right] = \frac{{\rm d} \Phi}{{\rm d} A}
放射輝度 (Radiance)
定義: 放射照度の立体角密度 (単位立体角あたりの照射照度)
L(x, \boldsymbol{\omega}) = \frac{{\rm d} E_\omega (x)}{{\rm d} \omega} \left[ \frac{\rm W}{\rm m^2 \cdot sr} \right]
ランバート反射 (Lambertian)
領域$ A $を$ \Phi [{\rm W}] $の光が通過
放射照度$ E $
E = \frac{\Phi}{A}
領域$ A_\omega$を$ \Phi [{\rm W}] $の光が通過
放射照度$ E_\omega $
A_\omega = A {\rm cos} \theta
E_\omega = \frac{\Phi}{A_\omega} = \frac{\Phi}{A {\rm cos} \theta}
ゆえに
E = E_\omega {\rm cos} \theta
BRDF (Bi-directional Reflectance Distribution Function)
双方反射分布関数。
入射する放射照度(Irradiance)に対して出射する放射輝度(Radiance)の比として表される。放射輝度の式を使って次のように変形できる。
f_r(\boldsymbol{x}, \boldsymbol{\omega}, \boldsymbol{\omega^\prime}) = \frac{{\rm d} L(\boldsymbol{x}, \boldsymbol{\omega})}{{\rm d} E_\omega(\boldsymbol{x}, \boldsymbol{\omega^\prime})} = \frac{{\rm d} L(\boldsymbol{x}, \boldsymbol{\omega})}{L_i(\boldsymbol{x}, \boldsymbol{\omega^\prime}) {\rm cos} \theta {\rm d} \boldsymbol{\omega^\prime}}
- BRDFは[0, 1]の範囲に縛られない。
- BRDFはunit-lessではない。($sr^{-1}$)
例
半球面で完全一様散乱(ランバート反射)の場合でのBRDFの正規化を行う
完全一様散乱の時、$f_r$は定数であるので次のようにおく
f_r = k_d
- $\Omega$ - 単位半球面
- $\rho$ - 反射率(albedo) (実装ではRGB)
\begin{array}{r}
\int_\Omega f_r {\rm cos} \theta {\rm d}\omega &=& \rho \\
\int_0^{2 \pi} \int_0^{\pi / 2} k_d {\rm cos} \theta {\rm sin} \theta {\rm d} \theta {\rm d} \phi &=& \rho \\
\pi k_d \int_0^{\pi / 2} {\rm sin} (2 \theta) {\rm d} \theta &=& \rho \\
\pi k_d &=& \rho
\end{array}
これより
f_r = \frac{\rho}{\pi}
レンダリング方程式
L_o(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})
- $ \boldsymbol{x} $ - 画像の各ピクセルについてカメラからレイを飛ばしたときのシーンとの交差点
- $ L_o $ - $ x $から$ \boldsymbol{\omega} $への放射輝度
- $ L_e $ - $ x $から$ \boldsymbol{\omega} $への自己発光による放射輝度
- $ f_r $ - $ x $におけるBRDF
- $ L_i $ - $ \boldsymbol{\omega^\prime} $方向から$ x $への放射輝度
- $ {\rm cos} \theta $ - 入射角度による放射輝度値の影響
この時$\boldsymbol{r}(\boldsymbol{x}, \boldsymbol{\omega})$を$x$から$\boldsymbol{\omega}$方向へのレイとシーンの交差点とすると
L_i(\boldsymbol{x}, \boldsymbol{\omega^\prime}) = L_o(\boldsymbol{r}(\boldsymbol{x}, \boldsymbol{\omega^\prime}), - \boldsymbol{\omega^\prime})
よって次式の様に再帰的に表せる
L_o(x, \boldsymbol{\omega}) = L_e(\boldsymbol{x}, \boldsymbol{\omega}) + \int_{\Omega_x} f_r(\boldsymbol{x}, \boldsymbol{\omega}, \boldsymbol{\omega^\prime}) L_o(\boldsymbol{r}(\boldsymbol{x}, \boldsymbol{\omega^\prime}), - \boldsymbol{\omega^\prime}) {\rm cos} \theta {\rm d} \sigma(\boldsymbol{\omega^\prime})
ロシアンルーレット
計算機は無限級数を利用することができないため、再帰の終了条件を作る必要性がある。その際に、統計的にUnbiasedなアルゴリズムとしてロシアンルーレットが知られている。
統計的にUnbiased = 期待値が真値と等価
C = C_0 + C_1 + C_2 + \cdots
*$\hat C$- 推定値
E[{\hat C}] = C
次のように$D_i$を定義すると
D_i = C_i + C_{i + 1} + C_{i+2} + \cdots
\begin{array}{l}
C &=& D_0 \\
&=& C_0 + D_1 \\
&=& C_0 + C_1 + D_2 \\
&=& \cdots
\end{array}
再帰的な表記にすると
C = D_0 \\
D_i = C_i + D_{i + 1}
処理を乱数によって分岐させる
*$r_0$は$(0, 1]$の乱数
*$p_0$は$(0, 1)$の定数
{\hat C} =
\begin{cases}
\frac{D_0}{p_0} & (r_0 \lt p_0) \\
0 & (r_0 \geq p_0)
\end{cases}
この時、期待値は
E[{\hat C}] = D_0 = C
- 実際の実装では$r_0$を生成し、$p_0$未満なら$D_0$の評価を行い、$p_0$以上なら処理を打ち切り評価値を$0$とする。 *$p_0$は基本的に任意の値で良いが、$D_0$の値の大小に比例させたほうが分散$V[{\hat C}]$は小さい。
次にロシアンルーレットによって$D_0$の評価を行うことになったとすると
{\hat D_0} = C_0 +
\begin{cases}
\frac{D_1}{p_1} & (r_1 \lt p_1) \\
0 & (r_1 \geq p_1)
\end{cases}
これを再帰的な表記にすると
{\hat C} =
\begin{cases}
\frac{D_0}{p_0} & (r_0 \lt p_0) \\
0 & (r_0 \geq p_0)
\end{cases}
{\hat D_i} = {C_i} +
\begin{cases}
\frac{D_{i + 1}}{p_{i + 1}} & (r_{i + 1} \lt p_{i + 1}) \\
0 & (r_{i + 1} \geq p_{i+ 1})
\end{cases}
このように確率的に処理を打ち切ることを選択することで、トータルの計算量を現実的な量に抑えつつ、期待値$E[{\hat C}] = C$と真値と等価(Unbiased)にすることができる。
* 最初の数項で$p_0 = 1$を与え、ロシアンルーレットを行わないことで、その数項でのロシアンルーレットのよる分散を抑えることができる。(最初の数項の寄与が大きいう特徴を知っているため)
* 実装上の注意点として、ある一定の$p_0$を利用するとたまたま再帰が深くなった時にスタックオーバーフローを起こしてしまうため、ある一定の深さ以上においては$p_0$を極端に小さくすることで深くなりすぎることを避けるようにする。
確率について
確率分布(確率測度)$P(D)$は$\Omega$上に定義される測度。確率変数$X$が集合$D$に含まれる確率とする。
対応する確率密度関数$p$は$\Omega$上の測度$\mu$を使ってRadon-Nikodym微分を行うと$p(x) = \frac{{\rm d} P}{{\rm d} \mu}$となりこれは$P = \int_D p(x){\rm d} \mu (x)$を満たすものである。
一般に、測度$\mu$に関する確率密度関数$p(x)$を$P_\mu (x)$と表記する。
例
- $S^2$ - 単位球面全体
- $\sigma(D)$ - 立体角測度。単位球面上のある領域を与えると、その領域の立体角を返す関数。(単位: ステラジアン)
- $P_\sigma (\boldsymbol{\omega})$ - $\sigma$に関する確率密度関数。$S^2$上のある一方向についての確率密度関数。
全球から一方向をサンプリングするとする。
一様にサンプリングするということは確率密度が全方向で定数なので
P_\sigma(\boldsymbol{\omega}) = C
確率分布の定義より$P(S^2) = 1$となるため
\int_{S^2} P_\sigma(\boldsymbol{\omega}) {\rm d} \sigma(\boldsymbol{\omega}) = \int_{S^2} C {\rm d} \sigma(\boldsymbol{\omega}) = 4 \pi C = 1
C = \frac{1}{4 \pi}
P_\sigma(\boldsymbol{\omega}) = \frac{1}{4 \pi}
モンテカルロ積分
$f(x)$が高次元、再帰的、非連続の場合、解析的に解を求めることは困難。
I = \int_D f(x) d \mu (x)
$f(x)$- 任意の関数
モンテカルロ推定器
{\hat I} = \frac{1}{N} \sum_{i=1}^N \frac{f(X_i)}{PDF(X_i)}
- $PDF(x)$- 確率密度関数
- $N$- サンプル数
- $X_i$- 確率変数(PDFによってサンプリングされる$X_i \in D$)
期待値$E[{\hat I}] = I$となる。Nを増加させていくとモンテカルロ推定値$\hat I$は真値$I$に対する誤差が$- O(\frac{1}{\sqrt{N}})$の速度で減少していく。$N$を∞に飛ばすと$I$に収束する。
- Consistentな推定器
- Nを∞に飛ばした時に$\hat I$が真値$I$に収束する。
- Unbiasedな推定器
- $\hat I$の期待値が真値$I$と等しい.
qiitaでは数式をひとつの記事にたくさん入れると重くてムリだったので、2つに分けました。
後半はこちら