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

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

• レイトレの数学的モデル
• 光の物理量
• BSDF(BRDF, BTDF)
• レンダリング方程式
• レンダリング方程式の作用素表現
• カメラのモデル化
• 経路積分形式のレンダリング方程式
• レンダリング方程式の計算について

レイトレの数学的モデル

まずレイトレーシングで行っていることを数学的に表現することを考えます。

方向

方向ベクトル$\vec{\omega}$は3次元の単位球面$\mathcal{S}^2 = \{ ||\vec{v}|| = 1 | \vec{v} \in \mathbb{R}^3 \}$上の点として表現します。こうすることで方向ベクトルは大きさが1に正規化されたもののみになります。

物体表面

シーンにある物体表面上の点全体を表す集合を$\mathcal{M}$とします。物体表面上で微積分が行えないと理論が展開できないので、$\mathcal{M}$は 区分的に微分可能な2次元多様体(2-dimensional Piecewise Differencial Manifold) として定義しておきます。

一般の空間上での積分

レイトレでは方向ベクトルによる積分や、物体表面上での積分などが頻繁に出てきます。これらの厳密な定義を知らなくても議論を進めることはできるのですが、ある程度知らないと何故？ってなると思うので、簡単に紹介します。

一般の空間$X$上での積分を定義するためには、測度(Measure)$\mu$ というものを導入し、測度空間(Measure Space)$(X, \mathcal{F}, \mu)$を構成する必要があります。

$\sigma$加法族

$\mathcal{F}$は $\sigma$加法族 といい、$X$の部分集合を持つ集合(集合族)です。$\sigma$加法族は次のような定義になっています。

1. $\phi \in \mathcal{F}$
2. $D \in \mathcal{F} \Longrightarrow D^c \in \mathcal{F}$
3. $D_1, D_2, ... \in \mathcal{F} \Longrightarrow \cup_{i=1}^{\infty} D_i \in \mathcal{F}$

この定義から $\mathcal{F}$に属する集合の和・差・積を高々加算無限回行って得られる集合はまた$\mathcal{F}$に属している ことが分かります。

この$\mathcal{F}$に属している集合を 可測集合(Measurable Set) といいます。後で見るように測度は$\mathcal{F}$に含まれる集合に対して定義されているので、可測集合でない集合では測度が計算できません。

測度

$\mathcal{F}$に含まれる集合$D_1, D_2, ... \in \mathcal{F}$に対して次のような集合関数$\mu$を定義します。

1. $0 \le \mu(D_i) \le {\infty}$ (非負性)
2. $\mu(\phi) = 0$
3. $D_i \cap D_j = \phi (i \neq j)$のとき、$\mu(\cup_{i=1}^{\infty}D_i) = \sum_{i=1}^{\infty}\mu(D_i)$ (可算加法性)

この集合関数を 測度(Measure) といいます。測度は与えられた集合の大きさを測ってくれる関数です。例えば$\mathbb{R}$上なら長さ、$\mathbb{R}^2$なら面積が測度の1つになります。

測度による積分

単関数の積分

$\mathcal{F}$に属し、互いに交わらない集合$D_1, D_2, ..., D_n$が与えられているとします。$D = \cup_{i=1}^{n} D_i$とします。$D_i$上で一定値$a_i \in \mathbb{R}$を取る、$D$上で定義された関数$f$を考えます。この関数は集合の定義関数$\boldsymbol{1}_{D_i}(x)$を用いて

として表現できます。このような関数$f$を 単関数(Simple Function) といいます。

$a_i \ge 0$となる単関数を$f \ge 0$で表すことにします。

単関数の積分を考えましょう。リーマン積分では、一変数関数の積分は$x$軸を一定間隔で分割し、その区間幅と関数の値をかけ合わせたものを足すことで積分を定義していました。

これと同様に、正の単関数$f \ge 0$の測度$\mu$による積分を

として定義します。（これは$D$の表し方によらない)

可測関数

$D$上の一般の関数で上のような積分を定義しようとすると、$\mu(D_i)$が計算できる必要があることが分かります。つまり、$f(x) = a_i$であるような$x$の集合を$f = a_i$と表すとすると、$(f = a_i) \in \mathcal{F}$である必要があります。

$\sigma$加法族の定義から、これは$(f < a_i) \in \mathcal{F}$と同値であることが分かります。これを満たす関数を 可測関数(Measurable Function) といいます。

一般の可測関数の積分

$D$上の任意の可測関数$f$の積分を考えます。$f$を

のように、値が正の部分と負の部分に分割します。このとき$f^+, f^- \ge 0$です。

$f(x)$はこの2つの関数を用いて$f(x) = f^+(x) - f^-(x)$と表せます。

$f$が可測であることは、$f^+$と$f^-$が両方とも可測であることと同値です。

まず$f^+$の積分を考えます。$f^+$に対しては、$D$上の単関数$f_i(x) \ge 0$の単調増加列で、$f$に各点収束するものが存在します。この単調増加列を$\{f_i\}$としたとき、$f^+$の測度$\mu$による積分を

として定義します(これは$\{f_i\}$の表し方によらない)。$f^-$の積分も同様に定義されます。

次に$f$の積分を考えます。$f = f^+ - f^-$と表されているので、$\int_E f^+d\mu$と$\int_E f^-d\mu$のどちらか一方が有限の値を取るとき、$f$の測度$\mu$による積分を

と定義します。この値が有限のとき、$f$は$D$上可積分であるといいます。

$\mathcal{S}^2$上での積分

$\mathcal{S}^2$上の表面積を表す測度を 立体角測度(Solid Angle Measure) といい、$\sigma$で表します。

これは平面の角度の概念を3次元空間に一般化したようなものです。$\Omega \in \mathcal{S}^2$上の可測関数$f$の立体角測度による積分を

と表します。(多様体論の微分形式からこのような測度が導かれるようです)

具体的な計算方法

この積分は半径1の球面座標系$(\theta, \phi)$を用いて表すことができます。立体角と球面座標系には

という関係があります。($(x, y, z) \rightarrow (r, \theta, \phi)$のヤコビアンを求めて$r=1$と置くと$\sin{\theta}$が得られる)

これを用いて立体角による積分は$\Omega$を区間$[\theta_1, \theta_2] \times [\phi_1, \phi_2]$として表現することで

と計算することができます。これは後で頻繁に使います。

$\mathcal{M}$上での積分

$\mathcal{M}$上の表面積を表す測度を 面積測度(Area Measure) といい、$A$で表します。

$D \in \mathcal{M}$上の可測関数$f$の面積測度による積分を

と表します。(これも立体角測度と同様に、多様体論からこのような測度が存在することが言えます)

レイ

レイは始点$\vec{o} \in \mathbb{R}^3$と、方向ベクトル$\vec{d} \in \mathcal{S}^2$を定義することで表現できます。レイ上にある全ての点集合$R$は

と表現できます。ここで$\mathbb{R}^+$は正の実数全体を表す集合です。

レイキャスト関数

シーン$\mathcal{M}$が与えられたとき、レイを飛ばして$\mathcal{M}$との衝突点$x'$を求める関数を レイキャスト関数(Raycasting Function) $r(x, \vec{d})$として定義します。

$t_{close}$は$\vec{o}$から$\vec{d}$にレイを飛ばしたときの衝突距離を表します。

光の物理量

物理ベースレンダリングを行うためには、まず光のどういった物理量を計算すれば良いのかを知らなくてはなりません。物理ベースレンダリングでは 放射測定学(Radiometry) と呼ばれる物理の分野で定義される物理量が計算の基礎になります。

光のエネルギー

波長$\lambda$の光子(Photon) 1つのエネルギー$e_{\lambda}$は

として定義されます。ここで$h$はプランク定数、$c$は光速です。

放射エネルギー

波長$\lambda$の光子が$n$個あるとき、スペクトル放射エネルギー$Q_{\lambda}$は

として定義されます。これを全波長に渡って積分することで、光の放射エネルギー$Q$が

として定義されます。これはそれぞれの波長で存在する光子の数を合計したものと見ることができます。

放射束

放射束(Radiant Flux) $\Phi$は単位時間あたりの放射エネルギーで、以下のように定義されます。

これは単位時間あたりの光子の数を表していると見ることができます。

放射照度

放射照度(Irradiance)$E$ は$\mathcal{M}$上の単位面積あたりの放射束を表し、

と定義されます。これは$\mathcal{M}$上のある点$x$を単位時間あたりに通過する光子の数と見ることができます。

放射輝度

放射輝度(Radiance)$L$は$\mathcal{M}$上の単位面積あたり、単位立体角あたりの放射束として定義され

と表されます。$\theta$は点$x$における法線$\vec{N(x)}$と方向ベクトル$\vec{\omega}$のなす角度です。$|\cos{\theta}|$は

として、2つのベクトルの内積をとることで計算できます。$|\cos{\theta}|dA$は物体表面上の単位面積$dA$を方向ベクトル$\vec{\omega}$の方向に投影した単位面積を表しています。

放射輝度は$\mathcal{M}$上のある点$x$にある方向$\vec{\omega}$からやってくる、単位時間あたりに通過する光子の数を表していると考えることができます。

物理ベースレンダリングでは放射輝度が最も重要な物理量です。レイトレではカメラセンサー上の点に、ある方向からやってくる光の強さを計算することになりますが、この光の強さというのがまさに放射輝度になっています。

放射輝度の重要な性質として、真空中では光が進行しても放射輝度の値が変化しない ことが挙げられます。真空のシーンでは物体表面上でしか放射輝度の変化が発生しないので、レイが空気中を進行している最中の光の作用を全て無視することができます。

放射輝度から放射照度を求める

放射輝度$L(x, \vec{\omega})$に$|\cos{\theta}|$をかけたものを立体角測度で積分することで放射照度が求まります。

放射輝度から放射束を求める

上で求めた$E$をさらに面積測度で積分することで放射束が求まります。

BRDF

光が物体に当たると反射・屈折が発生します。

点$x$から方向$\vec{\omega}$に出射していく放射輝度を$L_o(x, \vec{\omega})$、点$x$に方向$\vec{\omega}$から入射してくる放射輝度を$L_i(x, \vec{\omega})$と表すことにします。

BRDF(Bidirectional Reflectance Distribution Function - 双方向反射率分布関数) は物体表面上のある点$x$における、光の入射方向$\vec{\omega_i}$と反射方向$\vec{\omega_o}$が与えられたときの反射率を表す関数$f$で、次のように定義されます。

BRDFは入射してくる単位放射照度あたりの、反射される放射輝度の値を表しています。

このような定義になっているのは、実験結果から$E_i(x) \propto L_o(x, \vec{\omega_o})$という関係が知られているためです。

光の入射と反射方向は通常物体表面より上側に限られるため、$\vec{\omega_i}, \vec{\omega_o}$は半球面$\Omega = \{ \vec{\omega} | \vec{N(x)}\cdot \vec{\omega} \ge 0, \vec{\omega} \in \mathcal{S}^2\}$に属しています。

BRDFがエネルギー保存の法則を満たすために

となります。また、入射と反射の方向を入れ替えても同じなので

という関係があります。

$f$を指定することで、様々な反射特性を表現することができます。

拡散反射のBRDF

拡散反射(Diffuse Reflection) では光がすべての方向に一様に反射します。

これはBRDFの値が一定であることを意味しており、拡散反射のBRDFは一定の反射率を$\rho$として

と表されます。($\pi$が出てくるのはエネルギー保存則を満たすため)

完全鏡面のBRDF

完全鏡面の場合、入射した光は必ず一意の反射方向に反射するので、デルタ関数$\delta(\vec{\omega_i}, \vec{\omega_o})$を用いて

と表されます。

反射方程式

入射してくる放射輝度$L_i(x, \vec{\omega_i})$とBRDFが与えられている時、反射する放射輝度$L_o(x, \vec{\omega_o})$の値を計算してみましょう。

BRDFの定義から

両辺を立体角測度に関して全領域(半球面$\Omega$)で積分することで

この方程式を 反射方程式(Reflection Equation) といいます。

レンダリング方程式

物体が点$x$で発光し、方向$\vec{\omega_o}$に放射している放射輝度の値を$L_e(x, \vec{\omega_o})$で表すことにします。すると反射方程式と組み合わせることで、$x$から$\vec{\omega_o}$に反射する放射輝度の値は

として表されます。この方程式を レンダリング方程式(The Rendering Equation) といいます。(他にも光輸送方程式(LTE - Light Transport Equation)などの名前があります)

右辺第一項は発光を表しているので発光項、右辺第二項は反射を表しているので反射項と呼ばれます。

レンダリング方程式の重要性

レンダリング方程式は物理ベースレンダリングで最も重要な基礎方程式になっています。画像を生成するためにはカメラセンサー上の点$x$に、ある方向$\vec{\omega_i}$からやってくる放射輝度の値$L_i(x, \vec{\omega_i})$を知る必要があるのですが、この値をレンダリング方程式を用いて計算することができます。

つまり、物理ベースレンダリングとはレンダリング方程式を計算することに帰着されます。世の中には物理ベースレンダリングを達成するための様々なアルゴリズムが存在しますが、それらは全てこのレンダリング方程式をいかにして計算するかということを工夫したものです。

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

この方程式を計算するためには$L_e$の値、BRDFの値、そして$L_i$を知っている必要があります。$L_e$とBRDFはシーンを作った時点で既知ですが、$L_i$は未知なので計算する必要があります。

$L_i(x, \vec{\omega_i})$を計算するには再びレンダリング方程式を適用することができます。$x' = r(x, \vec{\omega_i})$とすると、$L_i(x, \vec{\omega_i}) = L_o(x', -\vec{\omega_i})$なので

として、$L_i(x, \vec{\omega_i})$が計算できます。

ここで$L_i(x', \vec{\omega_i'})$の計算に再びレンダリング方程式を適用します。このように再帰的にレンダリング方程式を使うことで、$L_o(x, \vec{\omega_o})$を計算することができます。これを理論的にきれいな形としてまとめたのが次の内容です。

レンダリング方程式の作用素表現

レンダリング方程式は再帰的な方程式になっていますが、これを関数解析で出てくる線形作用素を使用することで、簡単な代数方程式として表現することができます。

ここでは放射輝度の関数空間$\mathcal{L}$を考えます。この空間上で定義される線形写像$F:\mathcal{L} \rightarrow \mathcal{L}$を線形作用素(Linear Operator) といいます。線形作用素は関数を別の関数に移す線形写像です。身近な例としては微分・積分が線形作用素の1つになっています。

局所散乱作用素

レンダリング方程式の再帰をきれいな形として書くために、次のような局所散乱作用素(Local Scattering Operator)$S$を定義します。

局所散乱作用素$S$を使用することで、ある点の入射放射輝度$L_i$が与えられているときに、その点で反射される放射輝度$L_o$を

として表現できます。

進行作用素

さらに放射輝度の真空中の伝播を表すために、進行作用素$G$を定義します。

進行作用素$G$を利用することで、ある点から出射した放射輝度$L_o$が別の点に入射するとき、その点における入射放射輝度$L_i$を

として表現できます。

レンダリング方程式の作用素表現

局所散乱作用素$S$と進行作用素$G$を組み合わせることで、レンダリング方程式を

として表現することができます。ここでさらに$T = SG$を光輸送作用素として定義すれば

となります。これを レンダリング方程式の作用素表現 といいます。これは放射輝度の関数空間上で成り立つ代数的関係を表しています。

レンダリング方程式のノイマン級数展開

作用素表現のレンダリング方程式

は単純な代数方程式になっています。この方程式は関数方程式ですが、関数のあらゆる引数に対してこの関係が成り立っているので、シーン中の放射輝度の均衡関係を表していると見ることができます。(光子が光源から発射されて反射を繰り返していくと、シーン中のあらゆる地点で光子の単位時間あたりの通過量が一定になる)

これを変形することでレンダリング方程式の解を明示的に表現することが可能です。

右辺の$TL_o$を左辺に移行することで

ここで$I$は恒等作用素であり、$Ih = h$です。 $(I - T)$の逆作用素$(I - T)^{-1}$を用いることで

と書くことができます(逆作用素はBRDFがエネルギー保存の法則を満たすなら存在することが示されます)。右辺には$L_e$しかないので、これがレンダリング方程式の明示的な解になっています。

$(I - T)^{-1}$はノイマン級数展開することが可能で

と表されます。したがってレンダリング方程式は

と表現されます。これは レンダリング方程式のノイマン級数展開 と呼ばれます。この形では、出射放射輝度を求めるためには、シーン中の発光放射輝度を0回、1回、2回、$\cdots$, $\infty$回反射させたものを足し合わせれば良いということを言っています。

カメラのモデル化

光の伝達はレンダリング方程式としてモデル化することができました。次にカメラをモデル化することを考えてみます。

カメラはセンサー上に大量の画素を持ち、画素1つ1つが入射してくる光子の数を記録することによって光の強さを記録しています。これはつまり、非常に短い時間では放射束を記録しているということです。

カメラは入射してくる光をセンサー上に結像させるためにレンズを前面に持っています。レンズは一般的に複数枚の凹凸レンズを組み合わせて構成されます。

カメラからレイを出射することを考えます。物理的にはセンサーからレイが出発して、多くのレンズで屈折を繰り返してからシーンに出ていくようにモデル化した方が良いですが、ここでは話をシンプルにするために、途中の屈折過程はブラックボックス化し、レンズ最前面からいきなりレイが出射すると考えます。(画素の範囲外に出るレイを除けば画素上から出るレイと1対1対応する)

$j$番目の画素が記録する光の強さ$I_j$を考えます。レンズ最前面上の点を$x \in D$、ここから出るレイの方向を$\vec{\omega}$とします。このレイに対応する画素のセンサー感度を$W_e^{(j)}(x, \vec{\omega})$(これはMeasure Contribution Functionと呼ばれる)、レンズ最前面にレイの方向から入射してくる放射輝度を$L_i(x, \vec{\omega})$とします。すると$\theta$を$x$における法線ベクトルと$\vec{\omega}$のなす角として

となります。

$W_e$の関数形を変えることで、様々なカメラを表現することが可能です。

ピンホールカメラ

ピンホールカメラの場合、ピンホールの位置を$x_p$として、$D = \{x_p\}$です。また画素$j$に対応してレイの出射方向が一意に定まるため、$W_e$はデルタ関数で表現されることになります。したがって$I_j$を計算する式は積分をはずせて

のように簡略化されます。

経路積分形式のレンダリング方程式

レンダリング方程式は方向に関する積分方程式として表されていましたが、これを光の経路$(x_0, x_1, \cdots, x_k)$に関する積分として表現することも可能です。ここではそのような表現について見ていきます。

ある点$x$から別の点$x'$に向かう方向ベクトルを$x \rightarrow x'$として表すことにします。

立体角測度と面積測度の関係

物体表面上の点$x'$における微小単位面積$dA$を点$x$に置かれた単位球面上に投影することを考えます。すると

という関係があることが分かります。ここで$\cos{\theta'} = \vec{N(x')}\cdot (x' \rightarrow x)$です。$V(x, x')$は 可視関数(Visible Function) と呼ばれ、$x$と$x'$がお互いに見えるときに1、それ以外のときに0をとる関数です。

この関係式を用いて立体角測度と面積測度に関する積分を互いに積分変数変換することが可能になります。

3点形式のレンダリング方程式

ある点$x_1$からレンズ上の点$x_0$に出射する放射輝度$L_o(x_1 \rightarrow x_0)$を考えます。レンダリング方程式を立体角測度に関する積分から面積測度に関する積分に積分変数変換することによって、

と書くことができます。

$G(x_1, x_2)$は 幾何項(Geometry Term) と呼ばれ、積分変換変換による項と$|\cos{\theta}|$をかけたもので

と表されます。 このように視点$x_0$, 注目点$x_1$, 反射方向の点$x_2$を用いた形のレンダリング方程式を 3点形式のレンダリング方程式 といいます。

経路積分形式のレンダリング方程式

3点形式のレンダリング方程式をノイマン級数展開することで、次のような経路積分形式のレンダリング方程式が得られます。

カメラモデルの導入

カメラモデルの式

を経路積分形式のレンダリング方程式に組み込むことを考えます。立体角測度による積分を面積測度による積分に変換することで

となります。($D$を$\mathcal{M}$に変えていることに注意。このようにしても$W_e$を$D$以外の点で0にすることで値は変化しない)

$L_i(x_0 \rightarrow x_1) = L_o(x_1 \rightarrow x_0)$なので、経路積分形式のレンダリング方程式を代入することで

のように書けます。ここで$f(x_{-1} \rightarrow x_0 \rightarrow x_1) = 1$としています。

$\bar{x_k} = (x_0, x_1, \cdots, x_k)$とおきます。これを長さ$k+1$の光の 経路(Path) と呼びます。

$T_k(\bar{x_k}) = \prod_{i=1}^{k}\{ f(x_{i-2} \rightarrow x_{i-1} \rightarrow x_i)V(x_{i-1}, x_i)G(x_{i-1}, x_i) \}$を スループット(Throughput) として定義すると

のようになります。これが 経路積分形式(Path Integral Form) と呼ばれるレンダリング方程式の形です。

レンダリング方程式の計算について

再帰的な形のレンダリング方程式を経路積分形式で表現することで、レンダリング方程式は$n$次積分の無限級数になっていることが分かりました。数値積分を利用すればこれは数値計算が可能な形になっています。

しかし、10回の反射を経た光は10次元の積分、100回の反射を経た光は100次元の積分、$\cdots$というように、高次元の積分を数値的に計算する必要があります。

高次元の積分は通常の数値積分の手法だと計算量が爆発的に多くなります。例えばそれぞれの軸を10個のサンプル点に分割する場合、100次元積分を計算しようとすると$10^{100}$のサンプル点を生成することになります。これは現代のコンピューターを持ってしても計算が不可能な値です。

これを解決するための手法として、確率論的に積分の値を求める モンテカルロ積分(Montecarlo Integration) というものがあります。レイトレではモンテカルロ積分を用いてレンダリング方程式を解くことが一般的です。