この資料では 物理ベースレンダリング(Physically Based Rendering) を実現するために必要な数学・理論・実装について説明しています。

次のような事柄をこの資料では説明しています。

• 確率論
• レイトレでの確率論
• モンテカルロ積分
• レイトレでのモンテカルロ積分
• レンダリング方程式の数値計算

確率論

モンテカルロ積分を行うにあたって、確率論の基礎的な部分をある程度抑えている必要があります。ここではそれらの事柄について説明を行います。

確率空間

試行の結果得られる可能性のある標本全体を 標本空間(Sample Space) といい、$\Omega$で表します。

サイコロを一回投げるという例の場合、それぞれの目に対応する$\Omega = \{\omega_1, \omega_2, \cdots, \omega_6\}$を標本空間として取ります。

標本空間から生成されるσ加法族を$\mathcal{F}$で表し、$\mathcal{F}$の元を 事象(Event) と呼びます。特に$\Omega \in \mathcal{F}$を全事象、$\phi \in \mathcal{F}$を空事象と呼びます。

$\mathcal{F}$の元(つまり事象)に対し、確率測度(Probability Measure) $P:\mathcal{F} \rightarrow [0, 1]$が次のように定義されます。

1. $P$は$\mathcal{F}$上の測度
2. $P(\Omega) = 1$

確率測度が通常の測度と違うところは、$P(\Omega) = 1$が指定されているところです。このため任意の事象$A \in \mathcal{F}$に対し、$0 \le P(A) \le 1$を満たします。

３つ組$(\Omega, \mathcal{F}, P)$を 確率空間(Probability Space) と呼びます。

確率変数

可測関数$X:\Omega \rightarrow \mathbb{R}$を 実数値確率変数(Real Valued Random Variable) といいます。以降確率変数と言う場合は実数値確率変数を指すものとします。

この関数は標本空間の元$\omega \in \Omega$に対応する実数値を定めます。サイコロを1回投げる例では、サイコロの目の数という確率変数$X$を

のように定義できます。

分布測度

確率変数$X$がある値$B$を取る確率$P_X(B)$を計算することを考えます。この確率は確率変数の逆写像を用いて元の確率空間$(\Omega, \mathcal{F}, P)$に引き戻し、そこで$P$を用いることで計算できます。

$P_X$は 分布測度(Distribution Measure) と呼ばれます。

$P_X$は$\mathbb{R}$上の確率測度を定めます。したがって、確率変数$X$により別の確率空間$(\mathbb{R}, \mathcal{B}, P_X)$が誘導されます。($\mathcal{B}$はボレル集合族)

分布関数(累積分布関数)

以下のような関数$F_X(x)$を 分布関数(Distribution Function) または 累積分布関数(Cumulative Distribution Function - cdf) といいます。

分布関数は確率変数から導かれる分布測度を用いて定義されていますが、逆に分布関数を指定することで分布測度が一意に定まります。

分布関数を用いることで様々な確率が計算できます。例えば確率変数$X$が$[a, b]$に値をとる確率$P_X([a, b])$は

と計算できます。

確率密度関数

次のような条件を満たす$f_X(x)$を確率密度関数(Probability Density Function - pdf) といいます。 1. $\int_{-\infty}^{\infty} f(x)dx = 1$ 2. $F_X(x) = \int_{-\infty}^x f(x)dx$

この定義より$\frac{dF_X}{dx} = f(x)$となっていることに注意してください。（$\frac{dF_X}{dx} = f(x)$はRadon-Nikodym微分と呼ばれる。このような$f$が存在することはRadon-Nikodymの定理より言える)

確率密度関数を指定することで、分布関数が定まり、さらに分布測度が定まるので、確率密度関数を指定することで確率変数を構成することがよく行われます。

例: 一様分布

確率変数$X$が一様分布$U(a, b)$に従うとは、$X$の確率密度関数$f_X(x)$が

となっているときをいいます$(a < b)$。これは閉区間$[a, b]$から点を一様に取ることを意味しています。

例: 正規分布

確率変数$X$が正規分布$N(\mu, \sigma^2)$に従うとは、$X$の確率密度関数$f_X(x)$が

になっているときをいいます。確率変数が正規分布に従っていると仮定して、様々な解析を行うことが頻繁にあります。

期待値(平均)

確率変数$X$の 期待値(平均)$E[X]$ を分布測度$P_X$による積分として以下のように定義します。

$D$は$X$の値域です。期待値は$X$から出てくる値が平均的にどのような値を取るかを示しています。

注意としては$X$は常に可積分であるとは限らないので、期待値が計算出来ない場合もあります。

期待値は線形性を持ちます。

分散

確率変数$X$の 分散$V[X]$ を以下で定義します。

分散は$X$から出てくる値が平均周りでどれだけばらついているかを表している量です。分散が大きいことはばらつきが大きいことを示します。

分散も常に計算可能とは限らないことに注意してください。

分散は次のような性質を持ちます。

$\mathrm{Cov}[X, Y]$は共分散といい、

で定義されます。

また、分散はスカラー倍に対しては次のような性質を持ちます。

標準偏差

分散の平方根をとったものを 標準偏差(Standard Deviation) といい、$\sigma$で表します。

分散は元の量を二乗してしまった値のため、元の値の範囲でどれだけばらついているかが分かりにくいという問題点があります。標準偏差は元の値の範囲と一致するように分散を修正したものです。

確率変数変換

確率変数$X$によって誘導された確率空間を$(\mathbb{R}, \mathcal{B}_X, P_X)$、確率変数$Y$によって誘導された確率空間を$(\mathbb{R}, \mathcal{B}_Y, P_Y)$とします。

確率変数$X$と確率変数$Y$の間に、関数$h : \mathbb{R} \rightarrow \mathbb{R}$による関係$Y = h(X)$があるとします。もし$h$が単調増加関数なら

という関係が成り立つため、確率密度関数の間には

という関係が成り立ちます。同様に$h$が単調減少関数の場合もこの式が成り立ちます。つまり、$h$が単調関数ならこのようにして、確率変数変換後の確率密度関数を求めることができます。

多次元確率変数

上では実数値確率変数$X: \Omega \rightarrow \mathbb{R}$を考えましたが、値域を$\mathbb{R^2}$にした$X : \Omega \rightarrow \mathbb{R}^2$を考えます。

このような確率変数に対しても確率密度関数が同様に定義されますが、名前が変わります。

結合確率密度関数

次のような条件を満たす関数$p_{X, Y}(x, y)$を確率変数$X$の 結合確率密度関数(Joint Probability Density Function) といいます。

1. $\int_{\mathbb{R}^2} p_{X, Y}(x, y)dxdy = 1$
2. $P_X(D) = \int_D p_{X, Y}(x, y)dxdy$

周辺確率密度関数

結合確率密度関数$p(x, y)$から次のようにして求まる確率密度関数を 周辺確率密度関数(Marginal Probability Density Function) と呼びます。

これらはそれぞれの次元の実数値確率変数$X$と$Y$の確率密度関数になっています。

条件付き確率密度関数

$Y$の値$y$を観測した上での$X$の確率密度関数$p_{X|Y}(x)$を$X$の 条件付き確率密度関数(Conditional Probability Density Function) といいます。定義は以下の通りです。

確率変数の独立

確率変数$X$と$Y$が 独立(Independent) であるとは

となるときをいいます。

このとき$X$と$Y$の結合確率密度関数は

のように表わされます。

モンテカルロ積分

次のような一変数関数の積分を考えます。

点$x$は確率密度関数$p(x)$が指定され、確率変数$X$になっているとします($[a, b]$に値をとるとする)。$p(x)$により点が$(x_1, x_2, \cdots, x_N)$と$N$個生成されたとします。このとき$I$を以下のような計算で近似することができます。

これを$I$の モンテカルロ積分(Montecarlo Integration) といいます。

各項$\frac{f(X)}{p(X)}$の期待値を計算してみると

となり、各項の期待値は$I$に一致していることが分かります。したがって$E[\bar{I}]$も$I$に一致しているので、$\bar{I}$は$I$の 不偏推定量(Unbiased Estimator) になっています。

ここで 大数の法則(Law of large numbers) と呼ばれる次の定理を紹介します。

大数の法則

$X_1, X_2, \cdots, X_N$を平均$\mu$を持つ$n$個の確率変数とする。$\bar{X_n} = \frac{1}{n}\sum_{i=1}^N X_i$とすると、$n \to \infty$としたときに$\bar{X}$は確率1で$\mu$に収束する。

この定理より、モンテカルロ積分$\bar{I}$はNを無限大にすると確率1で求めたい積分$\bar{I}$に一致することが分かります。これがモンテカルロ積分で$I$が計算できる理由です。

モンテカルロ積分の誤差

一方で$\bar{I}$の分散を計算してみると

したがって分散は$O(N)$のオーダーで減少することが分かります。

標準偏差は推定誤差を表していますが、これは$O(\sqrt{N})$のオーダーで減少します。これはアルゴリズム論的にはあまり良くない性質です。例えば誤差を$\frac{1}{100}$にするためには、サンプル数を$10000$倍にしないといけないことを意味しています。

モンテカルロ積分の収束の遅さを改善するための手法として、重点的サンプリング(Importance Sampling) と呼ばれるものがあります。

高次元積分のモンテカルロ積分

モンテカルロ積分の利点は高次元積分を少ないサンプル数で計算が可能なことです。

以下のような高次元積分を考えます。

それぞれの次元の点を確率変数$X_1, \cdots, X_k$で表します。それぞれの次元の確率変数を互いに独立になるように設定すれば、$p_{X_1, \cdots, X_k} = p_{X_1} \cdots p_{X_k}$となるため、以下のようにモンテカルロ積分で$I$を求めることができます。

ここで$x_j^i$は次元$j$の$i$番目のサンプルを表します。

確率密度関数から点をサンプリングする方法

モンテカルロ積分を行うためには指定した確率密度関数からサンプリング点を生成する必要があります。次のような方法が有名です。

• 逆関数法(Inverse Transform Sampling)
• マルコフ連鎖モンテカルロ法(Markov Chain Monte Carlo - MCMC)

逆関数法

サンプリング点を生成したい確率変数を$X$(ここでは実数値確率変数とします)とし、その確率密度関数を$p_X(x)$とします。まず$X$の分布関数$F_X(x)$を求めます。

$F_X(x)$が狭義単調増加関数と仮定します。すると逆関数$F_X^{-1}(y)$が存在します。

ここで$U$を一様分布$U(0, 1)$に従う確率変数とすると、$F_X^{-1}(U)$の分布関数は

となり、$X$の分布関数$F_X(x)$に一致することが分かります。

したがって$[0, 1]$の一様乱数$u$を生成し、$x = F_X^{-1}(u)$とすることで確率変数$X$のサンプル$x$が得られます。

例1

以下のような積分をモンテカルロ積分で計算してみましょう。

$[0, \pi]$の一様分布を使ってサンプリング点$X$を生成することにします。$X$の確率密度関数は$p_X(x) = \frac{1}{\pi}$です。

逆関数法を使うと、$[0, 1]$の一様乱数$u$に対して$x = \pi u$とすれば$X$のサンプルが生成できます。

以下はC++でこれを実装したコードです。

#include <iostream>
#include <cmath>
#include <random>


std::random_device rnd_dev;
std::mt19937 mt(rnd_dev());
std::uniform_real_distribution<double> dist(0, 1);
inline double rnd() {
  return dist(mt);
}


inline double f(double x) {
  double v = std::sin(x);
  return v*v;
}


int main() {
  const int N = 10000;

  double sum = 0;
  for(int i = 0; i < N; i++) {
    double x = M_PI * rnd();
    sum += M_PI*f(x);
  }
  sum /= N;

  std::cout << sum << std::endl;

  return 0;
}

収束する様子は以下のようになります。

重点的サンプリング

$p(x)$を被積分関数$f(x)$に比例させるように取る方法です。

もし完全に比例させるように取ることができて、$p(x) = kf(x)$とできたとします。 このとき$\int_a^b p(x)dx = 1$という条件があるので$k = \frac{1}{I}$です。

この状況でモンテカルロ積分の分散を求めると

となり、分散が0なので確率1で$\bar{I}$は$I$に一致することが分かります。

実際には$p(x)$を$f(x)$に完全に比例するように取ることはできないですが、ある程度似るように取ることができれば収束が早くなることが期待できます。

例2

例1のモンテカルロ積分に重点的サンプリングを適用してみましょう。 被積分関数は$\sin^2(x)$なので、比例定数$k$を用いて$p_X(x) = k\sin x$と置きます。

確率密度関数は$\int_0^{\pi} p_X(x)dx = 1$を満たさないといけないので、これより$k = \frac{1}{2}$と求まります。したがって$p_X(x) = \frac{1}{2}\sin x$です。

$X$のサンプルを生成するために逆関数法を使用します。$X$の分布関数$F_X(x)$を求めると

となります。したがって$F_X^{-1}(y)$は

と求まります。$[0, 1]$の一様乱数$u$を用いて$x = F_X^{-1}(u)$として$X$のサンプルが生成できます。

以下はこれを実装したC++のコードです。

#include <iostream>
#include <cmath>
#include <random>


std::random_device rnd_dev;
std::mt19937 mt(rnd_dev());
std::uniform_real_distribution<double> dist(0, 1);
inline double rnd() {
  return dist(mt);
}


inline double f(double x) {
  double v = std::sin(x);
  return v*v;
}
inline double p(double x) {
  return 0.5 * std::sin(x);
}


int main() {
  const int N = 10000;

  double sum = 0;
  for(int i = 0; i < N; i++) {
    double u = rnd();
    double x = std::acos(1 - 2*u);
    sum += f(x)/p(x);
  }
  sum /= N;

  std::cout << sum << std::endl;

  return 0;
}

以下は収束の様子を重点的サンプリングを行わない場合を比較したものです。 重点的サンプリングを行ったほうが早く理論値の周辺に収束していることが分かります。

レイトレでのモンテカルロ積分

モンテカルロ積分をレンダリング方程式に適用しようとする場合、確率変数になるのは点の位置$x \in \mathcal{M}$や、方向ベクトル$\vec{\omega} \in \mathcal{S^2}$になります。つまり確率変数は実数値確率変数ではなく、$\mathcal{M}$-確率変数や$\mathcal{S^2}$-確率変数となっています。

確率空間の設定

確率空間は明示的には与えず、抽象的に$(\Omega, \mathcal{F}, P)$と置くことにします。

$\mathcal{S^2}$-確率変数

$\mathcal{S}^2$-確率変数$X_{\mathcal{S^2}} : \Omega \rightarrow \mathcal{S^2}$は方向ベクトルを生成する確率変数です。この確率変数から新たに確率空間$(\mathcal{S^2}, \mathcal{F}_{\mathcal{S^2}}, P_{\mathcal{S^2}})$が誘導されます。

$X_{\mathcal{S^2}}$の確率密度関数を$p_{\mathcal{S^2}}(\vec{\omega}) = \frac{dP_{\mathcal{S^2}}}{d\sigma}$と表すことにします。方向ベクトルの集合$D \in \mathcal{S^2}$が与えられたとき、この集合(事象)に対する確率$P_{\mathcal{S^2}}$は

と計算できます。

$\mathcal{M}$-確率変数

$\mathcal{M}$-確率変数$X_{\mathcal{M}} : \Omega \rightarrow \mathcal{M}$は物体表面上の点を生成する確率変数です。この確率変数から新たに確率空間$(\mathcal{M}, \mathcal{F}_{\mathcal{M}}, P_{\mathcal{M}})$が誘導されます。

$X_{\mathcal{M}}$の確率密度関数を$p_{\mathcal{M}}(x) = \frac{dP_{\mathcal{M}}}{dA}$と表すことにします。物体表面上の点集合$D \in \mathcal{M}$が与えられたとき、この集合(事象)に対する確率$P_{\mathcal{M}}$は

と計算できます。

$\mathcal{S^2}$-確率変数と$\mathcal{M}$-確率変数の変換

立体角測度と面積測度の間には次のような関係がありました。

確率変数変換の式から、$p_\mathcal{S^2}$と$p_\mathcal{M}$の間には次のような関係が成り立ちます。

この式を使うと、経路積分形式で表されたレンダリングの方程式において、方向のサンプリングと物体表面上の点のサンプリングを組み合わせて使うことができるようになります。

レンダリング方程式のモンテカルロ積分

経路積分形式のレンダリング方程式は以下のように表されていました。

これにモンテカルロ積分を適用することを考えます。経路点$x_0, \cdots, x_k$に対応する確率変数を$X_0, \cdots, X_k$とし、それぞれの確率変数は互いに独立であると仮定します。$i$番目のサンプルを$(x_0^i, \cdots, x_k^i)$と表し、それぞれの確率変数の確率密度関数を$p_{X_0}(x_0), \cdots, p_{X_k}(x_k)$とすると

のようにしてモンテカルロ積分することができます。無限級数の各項が積分になっているので、項ごとにモンテカルロ積分を適用することでこの形を得ます。

この形では各項の評価に毎回違うサンプルを使っていますが、これは計算コストが高くなるためにあまり良い方法ではありません。後で見るようにPath Tracingでは$x_0, x_1, \cdots$を逐次的に生成し、各項でサンプルの使い回しを行うことで計算コストを抑えています。

無限級数の評価

レンダリング方程式を数値計算するためには、上記の式のように無限級数を評価する必要があります。

次のような無限級数を考えます。

コンピューターで無限個の項を足し合わせることはできないので、有限個の項の総和で無限級数を近似する必要があります。

単純に足し合わせる最大の項数を設定して近似を行っても良いのですが、この場合それより後ろにある項が評価されないために、わずかですがバイアスが残ってしまいます。また、どこまで足し合わせれば良いかは必ずしも明確ではないです。

そこで使用されるのが ロシアンルーレット(Russian Roulette) と呼ばれる確率論的な評価方法です。

ロシアンルーレットによる無限級数の評価

ロシアンルーレットでは次のような手順で無限級数の評価を行います。

1. 次の項を足す確率$p \in [0, 1]$を設定します。
2. $i$回目の項を足すかどうか決めているとします。$[0, 1]$の一様分布に従う乱数$u$を生成します。
3. $u < p$なら$\frac{S_i}{p^i}$を足し合わせます。そうでないなら級数の評価を打ち切ります。

$S_i$ではなくて$\frac{S_i}{p^i}$を足していることがポイントです。$i$番目の項が足される確率は$p^i$なので、$\frac{S_i}{p^i}$で足せば期待値が元の無限級数と一致することになります。

次の項を足す確率$p$の設定

$p$を小さな値に設定してしまうと分散が大きくなってしまうので、最初は確率を1にし、徐々に確率を下げていく方法がよく取られます。

レンダリング方程式の数値計算

ここではレンダリング方程式のモンテカルロ積分を利用して、レンダリング方程式の数値解を計算するためのアルゴリズムについて考えます。

レンダリング方程式を数値計算するためのアルゴリズムとしては以下のようなものがあります。

• Path Tracing(PT)
• Bidirectional Path Tracing(BDPT)
• Metropolis Light Transport(MLT)
• Stochastic Progressive Photon Mapping(SPPM)

これらのアルゴリズムのどこに相違点があるかと言うと、大まかに言えばサンプルの生成方法に違いがあると言えます。例えば次に見るようにPath Tracingでは視点から逐次的に経路点のサンプルを生成していきますが、BDPTでは視点と光源の両方向から経路点のサンプルを生成していきます。

Path Tracing

パストレーシング(Path Tracing - 経路追跡法) は視点から経路点の生成を行っていく方法です。

ここでレンダリング方程式のモンテカルロ積分の式を再掲します。

この形では$(x_0, x_1, \cdots, x_k)$の各経路点を自由に生成することができます。しかし、もし本当に経路点を自由に生成したとすると、サンプリングされる経路の大半は可視関数$V(x_0, x_1)\cdots$のせいで0に評価されてしまうでしょう。これは明らかに無駄な計算を行っています。

Path Tracingでは経路点を直接サンプリングするのではなく、方向ベクトルをサンプリングすることで間接的に経路点をサンプリングします。

経路の始点$x_0$をサンプリングし、そこから先は方向ベクトル$\vec{\omega_0}, \cdots$をサンプリングすることで経路を以下のように生成することができます。

上のモンテカルロ積分を確率密度関数の変換公式を用いると

となり、$V$と$G$が消えてBRDFとコサイン項のみが残ることが分かります。($p_{X_0}(x_0)$が消えていることに注意。一般的なカメラの場合、$x_0$を定めると$\vec{\omega_0}$が定まるため、逆に$\vec{\omega_0}$を定めると$x_0$が定まってしまう。このとき$X_0$の確率密度関数にデルタ関数が含まれるので、1とみなしてよい)

ここでさらに、方向ベクトルの確率密度関数をBRDFとコサイン項に比例するように取って重点的サンプリングを行えば、分散を大幅に減少させることができます。このようなサンプリング方法をBRDF Samplingと呼びます。

さらにPath Tracingでは経路の生成を逐次的に行うことで、計算コストを抑えて効率よくモンテカルロ積分を行います。どういうことかというと、$N_1 = N_2 = \cdots = N$として、総和の順番を入れ替えることで(厳密には無限級数が絶対収束することを示さないといけない)、$\bar{I_j}$の計算が

のように置き換えられます。この形では長さ$k = 1$で$\vec{\omega_0}^i$をサンプリングすると、それが長さ$k \ge 2$の経路の評価においても使われるようになっています。つまり、方向ベクトルのサンプリングを繰り返すことで、経路を1つずつ長くしていくということです。

具体的には$i$番目のサンプル$(\vec{\omega_0}^i, \vec{\omega_1}^i, \cdots)$の評価において $k = 1$では

を評価して加算し、$k = 2$では

を評価して加算していきます。$k = 3$以降も同様です。

どこで経路を打ち切るかが問題になりますが、ここでロシアンルーレットを使用します。ロシアンルーレットの確率$p$を設定し、各$k$で項を評価する前に確率$1 - p$で項を評価を取りやめます。もし項を評価する場合は評価した値を$p^k$で割る必要があることに注意が必要です。

Path Tracingのアルゴリズム

まとめると、Path Tracingは以下のようなアルゴリズムです。

1. レンズ上の点をサンプリングし、始点$x_0$を決める。このとき$\vec{\omega_0}$も定まる。
2. ロシアンルーレットを行う。打ち切りなら終了。
3. $x_0$から$\vec{\omega_0}$にレイを飛ばし、衝突点や法線を求める。
4. 長さ$k$の項の評価を行う
5. BRDFを用いて$\vec{\omega_1}$の重点的サンプリングを行う。
6. 2に戻り、同様の手順を繰り返す。