このページではレンダリング方程式を数値計算する方法の一つであるモンテカルロ積分について説明する。
レンダリング方程式を再掲する。
L o ( x , ω o ⃗ ) = L e ( x , ω o ⃗ ) + ∫ Ω f ( x , ω i ⃗ , ω o ⃗ ) L i ( x , ω i ⃗ ) ( ω i ⃗ ⋅ n ⃗ ) d ω i ⃗ L_o(x, \vec{\omega_o}) = L_e(x, \vec{\omega_o}) + \int_{\Omega} f(x, \vec{\omega_i}, \vec{\omega_o})L_i(x, \vec{\omega_i})(\vec{\omega_i}\cdot\vec{n})d\vec{\omega_i} L o ( x , ω o ) = L e ( x , ω o ) + ∫ Ω f ( x , ω i , ω o ) L i ( x , ω i ) ( ω i ⋅ n ) d ω i 積分の項は数値積分を用いて計算することができる。L i ( x , ω i ⃗ ) L_i(x, \vec{\omega_i}) L i ( x , ω i ) 次元の呪い )
次元の呪い そこで評価点を乱数を用いてランダムに発生させ、その点で値を評価した結果を足し合わせていく方法がある。これが モンテカルロ積分(Monte Carlo Integration) である。
以下のような積分を評価することを考える。
I = ∫ a b f ( x ) d x I = \int_a^b f(x)dx I = ∫ a b f ( x ) d x サンプリング点{ x 1 , x 2 , … , x N } \{ x_1, x_2, \dots, x_N \} { x 1 , x 2 , … , x N } [ a , b ] [a, b] [ a , b ] p ( x i ) = 1 b − a p(x_i) = \frac{1}{b - a} p ( x i ) = b − a 1
I ≈ 1 N ∑ i = 1 N f ( x i ) p ( x i ) = b − a N ∑ i = 1 N f ( x i ) I \approx \frac{1}{N}\sum_{i = 1}^N \frac{f(x_i)}{p(x_i)} = \frac{b - a}{N}\sum_{i = 1}^N f(x_i) I ≈ N 1 i = 1 ∑ N p ( x i ) f ( x i ) = N b − a i = 1 ∑ N f ( x i ) で計算することができる。
E [ 1 N ∑ i = 1 N f ( x i ) p ( x i ) ] = 1 N ∑ i = 1 N E [ f ( x ) p ( x ) ] = 1 N ∑ i = 1 N ∫ f ( x ) p ( x ) p ( x ) d x = ∫ f ( x ) d x \begin{aligned}
E[\frac{1}{N}\sum_{i = 1}^N\frac{f(x_i)}{p(x_i)}] &= \frac{1}{N}\sum_{i = 1}^N E[\frac{f(x)}{p(x)}] \\
&= \frac{1}{N}\sum_{i=1}^N\int \frac{f(x)}{p(x)}p(x)dx \\
&= \int f(x)dx
\end{aligned} E [ N 1 i = 1 ∑ N p ( x i ) f ( x i ) ] = N 1 i = 1 ∑ N E [ p ( x ) f ( x ) ] = N 1 i = 1 ∑ N ∫ p ( x ) f ( x ) p ( x ) d x = ∫ f ( x ) d x 期待値が積分と一致することが分かる。
V [ 1 N ∑ i = 1 N f ( x i ) p ( x i ) ] = 1 N 2 ∑ i = 1 N V [ f ( x ) p ( x ) ] = 1 N 2 ∑ i = 1 N ∫ ( f ( x ) p ( x ) − I ) 2 p ( x ) d x = 1 N ∫ ( f ( x ) p ( x ) − I ) 2 p ( x ) d x \begin{aligned}
V[\frac{1}{N}\sum_{i=1}^N\frac{f(x_i)}{p(x_i)}] &= \frac{1}{N^2}\sum_{i=1}^N V[\frac{f(x)}{p(x)}] \\
&= \frac{1}{N^2}\sum_{i=1}^N\int (\frac{f(x)}{p(x)} - I)^2p(x)dx \\
&= \frac{1}{N}\int (\frac{f(x)}{p(x)} - I)^2p(x)dx
\end{aligned} V [ N 1 i = 1 ∑ N p ( x i ) f ( x i ) ] = N 2 1 i = 1 ∑ N V [ p ( x ) f ( x ) ] = N 2 1 i = 1 ∑ N ∫ ( p ( x ) f ( x ) − I ) 2 p ( x ) d x = N 1 ∫ ( p ( x ) f ( x ) − I ) 2 p ( x ) d x 分散はO ( N ) O(N) O ( N ) O ( N ) O(\sqrt{N}) O ( N ) 1 2 \frac{1}{2} 2 1
モンテカルロ積分の収束の様子 π \pi π
レンダリング方程式は無限次元積分になるが、モンテカルロ積分を使うことで数値計算することができる。
L o ( x , ω o ⃗ ) = L e ( x , ω o ⃗ ) + ∫ Ω f ( x , ω i ⃗ , ω o ⃗ ) L i ( x , ω i ⃗ ) ( ω i ⃗ ⋅ n ⃗ ) d ω i ⃗ L_o(x, \vec{\omega_o}) = L_e(x, \vec{\omega_o}) + \int_{\Omega} f(x, \vec{\omega_i}, \vec{\omega_o})L_i(x, \vec{\omega_i})(\vec{\omega_i}\cdot\vec{n})d\vec{\omega_i} L o ( x , ω o ) = L e ( x , ω o ) + ∫ Ω f ( x , ω i , ω o ) L i ( x , ω i ) ( ω i ⋅ n ) d ω i これをモンテカルロ積分の形式に直す。{ ω 1 ⃗ , ω 2 ⃗ , … , ω N ⃗ } \{\vec{\omega_1}, \vec{\omega_2}, \dots, \vec{\omega_N} \} { ω 1 , ω 2 , … , ω N } p ( ω i ⃗ ) p(\vec{\omega_i}) p ( ω i )
L o ( x , ω o ⃗ ) ≈ L e ( x , ω o ⃗ ) + 1 N ∑ i = 1 N f ( x , ω i ⃗ , ω o ⃗ ) L i ( x , ω i ⃗ ) ( ω i ⃗ ⋅ n ⃗ ) p ( ω i ⃗ ) L_o(x, \vec{\omega_o}) \approx L_e(x, \vec{\omega_o}) + \frac{1}{N}\sum_{i=1}^N\frac{f(x, \vec{\omega_i}, \vec{\omega_o})L_i(x, \vec{\omega_i})(\vec{\omega_i}\cdot\vec{n})}{p(\vec{\omega_i})} L o ( x , ω o ) ≈ L e ( x , ω o ) + N 1 i = 1 ∑ N p ( ω i ) f ( x , ω i , ω o ) L i ( x , ω i ) ( ω i ⋅ n ) で計算することができる。L i ( x , ω i ⃗ ) L_i(x, \vec{\omega_i}) L i ( x , ω i )
コンピューターで無限回の反射を計算することはできないので、一定の回数の反射を計算したところで計算を終了させる必要がある。計算を終える反射回数の決定方法として ロシアンルーレット(Russian Roulette) という方法がよく使われる。
次のような無限級数を評価することを考える。
S = L 1 + L 2 + … S = L_1 + L_2 + \dots S = L 1 + L 2 + … これを有限の項で打ち切って
S ≈ L 1 + L 2 + ⋯ + L N S \approx L_1 + L_2 + \dots + L_N S ≈ L 1 + L 2 + ⋯ + L N とすると誤差が残ってしまう。そこで次の項を足すかどうかを確率的に選択することで 統計的に偏りのない 計算をすることができる。
u u u [ 0 , 1 ] [0, 1] [ 0 , 1 ] p p p
u < p u < p u < p u ≥ p u \ge p u ≥ p
というルールで計算するのが ロシアンルーレット(Russian Roulette) である。
モンテカルロ積分とロシアンルーレットを組み合わせればレンダリング方程式を統計的に偏りなく計算することができる。
Path TracingではL i ( x , ω i ⃗ ) L_i(x, \vec{\omega_i}) L i ( x , ω i ) 2 10 2^{10} 2 1 0
Path Tracingのアルゴリズムは以下のようになる。
RGB Li ( Ray ray) {
if ( rand ( ) >= p) return RGB ( 0 )
レイを飛ばして物体と衝突するか調べる
if ( 衝突しなかった) return 空の色
if ( 物体が光っている) return Le
nextRay = 反射させたレイを生成
pdf = 生成したレイの確率密度
f = BRDF項
cos = コサイン項
return 1 / ( p* pdf) * f * cos * Li ( nextRay)
} この関数をN N N
通常のモンテカルロ積分は収束が遅いという問題を抱えている。収束を早めるための方法として以下のようなものがある。
重点的サンプリング(Importance Sampling)
層化サンプリング(Stratified Sampling)
準モンテカルロ積分(Quasi Monte Carlo Integration)
以下の積分を評価することを考える。
I = ∫ f ( x ) d x I = \int f(x)dx I = ∫ f ( x ) d x モンテカルロ積分の分散Vは
V = 1 N ∫ ( f ( x ) p ( x ) − I ) 2 p ( x ) d x V = \frac{1}{N}\int (\frac{f(x)}{p(x)} - I)^2p(x)dx V = N 1 ∫ ( p ( x ) f ( x ) − I ) 2 p ( x ) d x で与えられる。仮に確率密度関数をf ( x ) f(x) f ( x ) p ( x ) = k f ( x ) p(x) = kf(x) p ( x ) = k f ( x ) ∫ p ( x ) d x = 1 \int p(x)dx = 1 ∫ p ( x ) d x = 1
1 = k ∫ f ( x ) d x k = 1 I 1 = k\int f(x)dx \\
k = \frac{1}{I} 1 = k ∫ f ( x ) d x k = I 1 よって分散V V V
V = 1 N ∫ ( I − I ) 2 1 I f ( x ) d x = 0 V = \frac{1}{N}\int (I - I)^2 \frac{1}{I}f(x)dx = 0 V = N 1 ∫ ( I − I ) 2 I 1 f ( x ) d x = 0 したがって、確率密度関数を被積分関数f ( x ) f(x) f ( x )
重点的サンプリング このように確率密度関数p ( x ) p(x) p ( x ) f ( x ) f(x) f ( x ) 重点的サンプリング(Importance Sampling) という。
重点的サンプリングによる収束の様子
モンテカルロ積分ではサンプリング点をランダムに取るが、もしサンプリングされた点に偏りがあると分散が増大する原因になる。そこでサンプリングする領域を幾つかの小領域に分けて、必ず各小領域からサンプリング点を一つずつ取ることで、全体的に偏りなくサンプリング点を得ることができる。
このようにして分散を減少させる方法を 層化サンプリング(Stratified Sampling) という。
層化サンプリング法では小領域から乱数を生成することで分散の低減を図ったが、次元数が高くなると次元の呪いで小領域に分割するのが困難になる。
そこで、そもそも使用する乱数がある程度の規則性を持って分布させることを考える。このような規則性を持った乱数のような数のことを 超一様分布列(Low Discrepancy Sequence) または 準乱数 といい、超一様分布列を用いてモンテカルロ積分を行う方法を 準モンテカルロ積分(Quasi Monte Carlo Integration) という。
準モンテカルロ積分の誤差の標準偏差の収束オーダーはO ( N ) O(N) O ( N ) O ( N ) O(\sqrt{N}) O ( N )
超一様分布列として
ハルトン列(Halton Sequence)
ソボル列(Sobol Sequence)
といったものがある。
Halton列 [ 0 , 1 ] × [ 0 , 1 ] [0, 1]\times[0,1] [ 0 , 1 ] × [ 0 , 1 ]
Halton列によるモンテカルロ積分
拡散反射面ではf ( x , ω i ⃗ , ω o ⃗ ) = ρ π f(x, \vec{\omega_i}, \vec{\omega_o}) = \frac{\rho}{\pi} f ( x , ω i , ω o ) = π ρ
L o ( x , ω o ⃗ ) = L e ( x , ω o ⃗ ) + ρ π ∫ Ω L i ( x , ω i ⃗ ) ( ω i ⃗ ⋅ n ⃗ ) d ω i ⃗ L_o(x, \vec{\omega_o}) = L_e(x, \vec{\omega_o}) + \frac{\rho}{\pi}\int_{\Omega}L_i(x, \vec{\omega_i})(\vec{\omega_i}\cdot\vec{n})d\vec{\omega_i} L o ( x , ω o ) = L e ( x , ω o ) + π ρ ∫ Ω L i ( x , ω i ) ( ω i ⋅ n ) d ω i となる。積分の中身に注目するとコサイン項( ω i ⃗ ⋅ n ⃗ ) = cos θ (\vec{\omega_i}\cdot\vec{n}) = \cos{\theta} ( ω i ⋅ n ) = cos θ ω i ⃗ \vec{\omega_i} ω i p ( ω i ⃗ ) p(\vec{\omega_i}) p ( ω i ) cos θ \cos{\theta} cos θ
一般的に拡散反射面以外でも、BRDFf ( x , ω i ⃗ , ω o ⃗ ) f(x, \vec{\omega_i}, \vec{\omega_o}) f ( x , ω i , ω o ) ( ω i ⃗ ⋅ ω o ⃗ ) (\vec{\omega_i}\cdot\vec{\omega_o}) ( ω i ⋅ ω o ) p ( ω i ⃗ ) p(\vec{\omega_i}) p ( ω i )
cos θ \cos{\theta} cos θ p ( ω i ⃗ ) p(\vec{\omega_i}) p ( ω i ) 比例定数k k k
p ( ω i ⃗ ) = k cos θ p(\vec{\omega_i}) = k\cos{\theta} p ( ω i ) = k cos θ とおく。∫ Ω p ( ω i ⃗ ) d ω i ⃗ = 1 \int_{\Omega} p(\vec{\omega_i})d\vec{\omega_i} = 1 ∫ Ω p ( ω i ) d ω i = 1
k ∫ Ω cos θ d ω i ⃗ = 1 k\int_{\Omega} \cos{\theta}d\vec{\omega_i} = 1 k ∫ Ω cos θ d ω i = 1 立体角を球面座標系に変換するとd ω i ⃗ = sin θ d θ d ϕ d\vec{\omega_i} = \sin{\theta}d\theta d\phi d ω i = sin θ d θ d ϕ
k ∫ 0 2 π ∫ 0 π 2 cos θ sin θ d θ d ϕ = 1 k\int_0^{2\pi}\int_0^{\frac{\pi}{2}}\cos{\theta}\sin{\theta}d\theta d\phi = 1 k ∫ 0 2 π ∫ 0 2 π cos θ sin θ d θ d ϕ = 1 計算していくと
k = 1 π k = \frac{1}{\pi} k = π 1 よって
p ( ω i ⃗ ) = cos θ π p(\vec{\omega_i}) = \frac{\cos{\theta}}{\pi} p ( ω i ) = π cos θ と求まる。
p ( ω i ⃗ ) p(\vec{\omega_i}) p ( ω i ) 確率密度関数p ( ω i ⃗ ) p(\vec{\omega_i}) p ( ω i )
逆関数法(Inverse Transform Method)
棄却サンプリング法(Rejection Sampling)
マルコフ連鎖モンテカルロ法(Markov Chain Monte Carlo Method)
逆関数法は以下のような手順で確率密度関数p ( x ) p(x) p ( x )
[ 0 , 1 ] [0, 1] [ 0 , 1 ] u u u 累積分布関数F ( x ) F(x) F ( x ) F − 1 ( x ) F^{-1}(x) F − 1 ( x ) u u u
方向ω i ⃗ \vec{\omega_i} ω i
p i = p ( ω i ⃗ ) p_i = p(\vec{\omega_i}) p i = p ( ω i ) [ 0 , 1 ] [0, 1] [ 0 , 1 ] u u u u < p u < p u < p ω i ⃗ \vec{\omega_i} ω i u ≥ p u \ge p u ≥ p
省略
累積分布関数F ( x ) F(x) F ( x )
p ( ω i ⃗ ) = cos θ π p(\vec{\omega_i}) = \frac{\cos{\theta}}{\pi} p ( ω i ) = π c o s θ 逆関数法を用いてp ( ω i ⃗ ) = cos θ π p(\vec{\omega_i}) = \frac{\cos{\theta}}{\pi} p ( ω i ) = π c o s θ ω i ⃗ \vec{\omega_i} ω i
立体角のままでは扱いずらいので、球面座標系に変換することを考える。立体角による確率密度関数p ( ω ⃗ ) p(\vec{\omega}) p ( ω ) p ( θ , ϕ ) p(\theta, \phi) p ( θ , ϕ )
p ( ω ⃗ ) sin θ = p ( θ , ϕ ) p(\vec{\omega})\sin{\theta} = p(\theta, \phi) p ( ω ) sin θ = p ( θ , ϕ ) という関係がある。
これを利用しp ( ω i ⃗ ) = cos θ π p(\vec{\omega_i}) = \frac{\cos{\theta}}{\pi} p ( ω i ) = π c o s θ
p ( θ , ϕ ) = cos θ sin θ π p(\theta, \phi) = \frac{\cos{\theta}\sin{\theta}}{\pi} p ( θ , ϕ ) = π cos θ sin θ となる
これはθ \theta θ ϕ \phi ϕ θ \theta θ p ( θ ) p(\theta) p ( θ )
θ \theta θ p ( θ ) p(\theta) p ( θ )
p ( θ ) = ∫ 0 2 π cos θ sin θ π d ϕ = 2 cos θ sin θ = sin 2 θ p(\theta) = \int_0^{2\pi} \frac{\cos{\theta}\sin{\theta}}{\pi}d\phi = 2\cos{\theta}\sin{\theta} = \sin{2\theta} p ( θ ) = ∫ 0 2 π π cos θ sin θ d ϕ = 2 cos θ sin θ = sin 2 θ よってϕ \phi ϕ p ( ϕ ∣ θ ) p(\phi | \theta) p ( ϕ ∣ θ )
p ( ϕ ∣ θ ) = p ( θ , ϕ ) p ( θ ) = 1 2 π p(\phi | \theta) = \frac{p(\theta, \phi)}{p(\theta)} = \frac{1}{2\pi} p ( ϕ ∣ θ ) = p ( θ ) p ( θ , ϕ ) = 2 π 1 と求まる。
θ \theta θ P ( θ ) P(\theta) P ( θ )
P ( θ ) = ∫ 0 θ sin 2 t d t = 1 − cos 2 θ 2 P(\theta) = \int_0^{\theta}\sin{2t}dt = \frac{1 - \cos{2\theta}}{2} P ( θ ) = ∫ 0 θ sin 2 t d t = 2 1 − cos 2 θ 同様に
P ( ϕ ∣ θ ) = ∫ 0 ϕ 1 2 π d t = 1 2 π ϕ P(\phi | \theta) = \int_0^{\phi}\frac{1}{2\pi}dt = \frac{1}{2\pi}\phi P ( ϕ ∣ θ ) = ∫ 0 ϕ 2 π 1 d t = 2 π 1 ϕ P ( θ ) P(\theta) P ( θ )
P θ − 1 ( u ) = 1 2 cos − 1 ( 1 − 2 u ) P_{\theta}^{-1}(u) = \frac{1}{2}\cos^{-1}(1 - 2u) \\ P θ − 1 ( u ) = 2 1 cos − 1 ( 1 − 2 u ) P ( ϕ ∣ θ ) P(\phi | \theta) P ( ϕ ∣ θ )
P ϕ ∣ θ − 1 ( v ) = 2 π v P_{\phi | \theta}^{-1}(v) = 2\pi v P ϕ ∣ θ − 1 ( v ) = 2 π v と求まる。
あとはu , v u, v u , v [ 0 , 1 ] [0, 1] [ 0 , 1 ]
θ = 1 2 cos − 1 ( 1 − 2 u ) ϕ = 2 π v \theta = \frac{1}{2}\cos^{-1}(1 - 2u) \\
\phi = 2\pi v θ = 2 1 cos − 1 ( 1 − 2 u ) ϕ = 2 π v とすれば良い。
あとは球面座標系から直交座標系( x , y , z ) (x, y, z) ( x , y , z )
x = cos ϕ sin θ y = cos θ z = sin ϕ sin θ \begin{aligned}
x &= \cos{\phi}\sin{\theta} \\
y &= \cos{\theta} \\
z &= \sin{\phi}\sin{\theta}
\end{aligned} x y z = cos ϕ sin θ = cos θ = sin ϕ sin θ プログラムにすると次のようになる。
Vec3 randomCosineHemisphere ( double & pdf, const Vec3& n) {
double u = rnd ( ) ;
double v = rnd ( ) ;
double theta = 0.5 * std:: acos ( 1 - 2 * u) ;
double phi = 2 * M_PI* v;
pdf = 1 / M_PI * std:: cos ( theta) ;
double x = std:: cos ( phi) * std:: sin ( theta) ;
double y = std:: cos ( theta) ;
double z = std:: sin ( phi) * std:: sin ( theta) ;
Vec3 xv, zv;
orthonormalBasis ( n, xv, zv) ;
return x* xv + y* n + z* zv;
} pdfにはサンプリングされた方向ω i ⃗ \vec{\omega_i} ω i p ( ω i ⃗ ) p(\vec{\omega_i}) p ( ω i )
これを使うと拡散反射面の重点的サンプリングは次のように書ける。
Vec3 getRadiance ( const Ray& ray, int depth = 0 ) {
. . .
if ( res. hitSphere-> type == "diffuse" ) {
double pdf;
Vec3 nextDir = randomCosineHemisphere ( pdf, res. hitNormal) ;
Ray nextRay ( res. hitPos + 0.001 * res. hitNormal, nextDir) ;
double cos_term = std:: max ( dot ( nextDir, res. hitNormal) , 0.0f ) ;
return 1 / roulette * 1 / pdf * res. hitSphere-> color/ M_PI * cos_term * getRadiance ( nextRay, depth + 1 , roulette) ;
}
. . .
} 半球一様サンプリング 重点的サンプリング 上の画像は同じ100サンプル数で一様サンプリングと重点的サンプリングの結果を比較したものである。重点的サンプリングを使ったほうがノイズが減少しているのが分かる。
