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


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

前半はこちら

[物理ベースレンダリング1](http://qiita.com/pnlybubbles/items/416a3273cd245d1a1c04)

---

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

* 物理ベースレンダラ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で書きました。まだおそい。

https://github.com/pnlybubbles/path-trace

## パストレーシング

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

必要なものは

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

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

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

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

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

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

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

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

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

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

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

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

確率密度関数

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

累積分布関数

```math
\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}
```

パラメータを分離すると

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

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

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

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

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

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

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

```math
\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}
```

```math
\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}
```

```math
\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}
```

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

```math
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$に近づく。

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

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

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

これらを代入すると

```math
\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}
```

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

```math
{\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$にしたがって重点的にサンプリングしたい場合、確率密度関数を次のようにして求める。

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

```math
\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}
```

```math
C = \frac{1}{\pi}
```

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

確率密度関数

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

累積分布関数

```math
\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}
```

パラメータを分離すると

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

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

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

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

```math
\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}
```

```math
\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}
```

```math
\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}
```

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

```math
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})
```

```math
\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}
```

```math
{\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
![image_16.png](https://qiita-image-store.s3.amazonaws.com/0/14939/e3a0e152-ccf6-91e3-f3f8-d571b9fe3e66.png)

64 sample
![image_64.png](https://qiita-image-store.s3.amazonaws.com/0/14939/e9126f09-5165-15ab-28ce-2cbe5c39e904.png)

256 sample
![image_256.png](https://qiita-image-store.s3.amazonaws.com/0/14939/9022504e-a961-ffbe-b669-6545636a9301.png)

1024 sample
![image_1024.png](https://qiita-image-store.s3.amazonaws.com/0/14939/1641991f-5c9a-c122-e4b6-4667175f36e4.png)

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