逆関数法とジオメトリサンプリング実装例

図1. 目標分布 $ f(x) = 1 \,/\, (x + 0.5) $

例題として目標分布を図1に示すような $ [0, 1.5) $ における $ f(x) = 1 \,/\, (x + 0.5) $ としてみます。まずはこれに比例したPDF $ p(x) $ を得るために比例係数(の逆数)を以下のようにして求めます。 \begin{equation*} \int_{0}^{1.5} \frac{1}{x + 0.5} \d x = \Brkb{\ln \Prt{x + 0.5}}_{0}^{1.5} = \ln 4 \end{equation*} したがってサンプルが沿うことになるPDF $ p(x) $ は次のようになります。 \begin{equation*} p(x) = \frac{1}{\ln 4} \cdot \frac{1}{x + 0.5} \end{equation*} 次にこれを積分することでCDF $ P(x) $ を求めます。 \begin{equation*} P(x) = \int p(x) \d x = \frac{1}{\ln 4} \ln (x + 0.5) + A \end{equation*} CDFとしての境界条件 $ P(0) = 0 $ や $ P(1.5) = 1 $ を用いれば積分定数 $ A $ を求めることができて $ P(x) $ は次のようになります。 \begin{equation*} P(x) = \frac{1}{\ln 4} \ln \frac{x + 0.5}{0.5} \end{equation*} これを基に逆関数を計算すると次の式を得ます。 \begin{equation*} P^{-1}(u) = 0.5 \cdot \Prt{4^u - 1} \end{equation*} 図2にこの式を用いて100,000サンプル生成したときのヒストグラムとPDF $ p(x) $ を示します。$ p(x) $ に沿ったサンプリングが実現できていることがわかります。

図2. PDF $ p(x) $ (青線) と実際に達成された分布のヒストグラム

逆関数が解析的に計算できない場合

図3. (a)同時確率密度関数 $ p(x, y) $ と(b)達成された分布のヒストグラム

ここでは例として2次元関数の目標分布 $ f(x, y) = x + 2y $、$ (0 \leq x, y < 1) $ を使います。これを積分してPDF $ p(x, y) $ を求めると次の式を得ます。 \begin{equation*} p(x, y) = \frac{2}{3} (x + 2y) \end{equation*} これを色を使って可視化すると図3(a)のようになります。続いてこれを $ y $ に関して積分することで $ x $ に関する周辺確率密度関数 $ p(x) $ を得ます。 \begin{equation*} p(x) = \frac{2}{3} (x + 1) \end{equation*} このPDFを積分して逆関数を求めると次の式を得ます。 \begin{equation*} P^{-1}(u) = -1 + \sqrt{1 + 3u} \end{equation*} これによって1つ目の一様乱数 $ u_x $ を用いて $ x = P^{-1}(u_x) $ として $ x $ のサンプリングを行うことができます。続いて $ x $ が与えられた場合の $ y $ の条件付き確率密度関数を次のように求めます。 \begin{equation*} p(y|x) = \frac{\displaystyle \frac{2}{3} (x + 2y)}{\displaystyle \frac{2}{3} (x + 1)} = \frac{x + 2y}{x + 1} \end{equation*} この式を積分して $ y $ に関する逆関数を求めると次の式を得ます。 \begin{equation*} P^{-1}(u|x) = \frac{-x + \sqrt{x^2 + 4(x + 1)u}}{2} \end{equation*} 2つ目の一様乱数 $ u_y $ と先にサンプルしておいた $ x $ を用いて$ y = P^{-1}(u_y|x) $ として $ y $ のサンプリングを行うことができます。これら2つの連なった逆関数を用いて1,000,000サンプル生成したときのヒストグラムを可視化したものが図3(b)になります。PDFに沿ったサンプリングができていることがわかります。

球面や半球面上のサンプリングは方向のサンプリングなどに用いられます。PDFを基に逆関数法を使ってサンプリング方法を導出することができますが、注意点として球面座標における積分では「確率密度関数の変換」で示したようにヤコビアンを考える必要があります。立体角に関するPDF $ p_\sigma(\vomega) $ は球面座標上におけるPDF $ p(\theta, \phi) $ と次の関係を持ちます。 \begin{eqnarray*} p_\sigma(\vomega) \d \vomega &=& p(\theta, \phi) \d\theta \d\phi \\ p(\theta, \phi) &=& \sin\theta \cdot p_\sigma(\vomega) \end{eqnarray*}

図4. 球面・半球面上のサンプリング

(a)球面上一様分布、(b)半球面上一様分布、(c)半球面上cos分布

球面内一様分布

半球面内一様分布

半球面内cos分布

平面内のサンプリングはレンダリングにおいては、光源上の点のサンプリングなどに用いられます。円内のサンプリングは光源に加えて簡単なレンズモデルにおけるレンズ上の点のサンプリング、そして一様分布ではない可能性もありますが、ある面上の点の周辺領域のサンプリングにも応用されます。

長方形(平行四辺形)内一様分布

図5. 平行四辺形内一様分布

平行四辺形内一様分布のPDFは次のような形で表されます。 \begin{equation*} p_A(\vx) = \frac{1}{\absb{\vec{a} \times \vec{b}}} \hspace{3mm} \Brk{\mathrm{m}^{-2}} \end{equation*} ここで $ \vec{a} $ と $ \vec{b} $ は平行四辺形の形状を決める2つのベクトルです(2次元でも3次元でも構いません)。$ \absb{\vec{a} \times \vec{b}} $ は平行四辺形の面積に等しくなるので、PDFは面積の逆数に等しくなります。平行四辺形中の位置 $ \vec{x} $ のサンプリングは次のようになります。 \begin{equation*} \vec{x} = \vec{o} + u_1 \vec{a} + u_2 \vec{b} \end{equation*} ここで $ \vec{o} $ は平行四辺形の基準点を表しています。ここでは逆関数法の計算手順の詳細は示しませんでしたが、$ \vec{a} $ と $ \vec{b} $ を基底ベクトルとする座標系における座標 $ a, b $, $ (0 \leq a, b < 1) $ を考えることで求めることができます。図5に $ \vec{o} = (2, 1) $, $ \vec{a} = (4, 1) $, $ \vec{b} = (1, 4) $ とした場合のサンプリング結果を示しています。

円内一様分布

図6. 円内一様分布 ($ R = 2 $)

円内一様分布のPDFは次のように円の面積の逆数として表されます。 \begin{equation*} p_A(\vx) = \frac{1}{\pi R^2} \hspace{3mm} \Brk{\mathrm{m}^{-2}} \end{equation*} ここで $ R $ は円の半径です。

単純な実装

単純に円座標(極座標)内でヤコビアンに気をつけつつサンプリングを導出すると次の形を得ます。 \begin{eqnarray*} r &=& R \sqrt{u_1} \\ \theta &=& 2 \pi u_2 \end{eqnarray*} 図6に $ R = 2 $ とした場合のサンプリング結果を示しています。

歪みの少ない実装

図7. 円へのマッピング手法比較

(a)単純な手法、(b)concentric mapping

上の単純な実装でも円内一様分布としては間違いないのですが、入力となる2つの一様乱数が入る $ [0, 1)^2 $ の領域を多数の小領域に区切って考えたとき、各小領域の円内へのマップのされ方は場所によって大きく異なり、非常に歪みの大きいところが存在します。層化サンプリングや準モンテカルロ法、ブルーノイズサンプリングなどの適用を考える場合には、元のサンプル列が持つ優れた特性を破壊してしまうため好ましくありません。[Shirley1997]では歪みが比較的均等なマッピングであるconcentric mappingが提案されています。理論の紹介は省きますが、単位円のサンプリングに関して次に示すような手順を実行します。

まず2つの一様乱数を $ [-1, 1) $ にマップします。 \begin{eqnarray*} o_x &=& 2 * u_1 - 1 \\ o_y &=& 2 * u_2 - 1 \end{eqnarray*} 続いて $ o_x $, $ o_y $ 間の関係に応じて場合分け処理を行い $ r $ と $ \theta $ を得ます。ここで得る $ r $ と $ \theta $ は負の値にもなり得ることに注意が必要です。

$ o_x = 0 $, $ o_y = 0 $ の場合	\begin{eqnarray*} r = 0 \\ \theta = 0 \end{eqnarray*}
$ \abs{o_x} > \abs{o_y} $ の場合	\begin{eqnarray*} r &=& o_x \\ \theta &=& \frac{\pi}{4} \frac{o_y}{o_x} \end{eqnarray*}
それ以外($ \abs{o_x} \leq \abs{o_y} $)の場合	\begin{eqnarray*} r &=& o_y \\ \theta &=& \frac{\pi}{2} - \frac{\pi}{4} \frac{o_x}{o_y} \end{eqnarray*}

最後に直交座標への変換を行います。 \begin{eqnarray*} x &=& r \cos\theta \\ y &=& r \sin\theta \end{eqnarray*} 図7に単純なマッピングとconcentric mappingの分布の比較を示しています。入力として一様乱数の代わりに $ [0, 1)^2 $ の空間内で比較的均等に分布したサンプルを用いました。どちらも全体としては一様分布ですが、(b) concentric mappingは(a)単純なマッピングに比べて少しサンプルの粗密のゆらぎが抑えられていることがわかるでしょうか。

三角形内の一様分布

図8. 三角形内一様分布

三角形内一様分布のPDFは次のような形で表されます。 \begin{equation*} p_A(\vx) = \frac{1}{0.5 \cdot \absb{\vec{a} \times \vec{b}}} \hspace{3mm} \Brk{\mathrm{m}^{-2}} \end{equation*} ここで $ \vec{a} $ と $ \vec{b} $ は平行四辺形の場合と同様に、三角形の形状を決める2つのベクトルです(2次元でも3次元でも構いません)。$ 0.5 \cdot \absb{\vec{a} \times \vec{b}} $ は三角形の面積に等しくなるので、PDFは面積の逆数に等しくなります。

単純な実装

三角形中の位置 $ \vec{x} $ のサンプリングは単純な実装の場合次のようになります。 \begin{equation*} t_a = 1 - \sqrt{u_1} \\ t_b = (1 - t_a) u_2 \\ \vec{x} = \vec{o} + t_a \vec{a} + t_b \vec{b} \end{equation*} ここで $ \vec{o} $ は三角形の基準点を表しています。一旦入力サンプルを $ \vec{a} $, $ \vec{b} $ を基底とするローカル空間中のパラメター $ t_a $, $ t_b $ に変換し、それらを用いて $ \vec{x} $ を表します。図8に $ \vec{o} = (1, 1) $, $ \vec{a} = (4, 1) $, $ \vec{b} = (1, 4) $ とした場合のサンプリング結果を示しています。

歪みの少ない実装

前述の円内一様分布の場合と同様に、三角形に関してもマッピングの歪みが少ない実装が[Heitz2019]で提案されています。サンプリング手順は次のようになります。

function $ \mathrm{squareToTriangle} $ (単位矩形 $ [0, 1)^2 $ 中のサンプル $ u_1 $, $ u_2 $) {
	$ t_1 \LAR 0.5 u_1 $
	$ t_2 \LAR 0.5 u_2 $
	$ t_\mathrm{offset} = t_2 - t_1 $
	if $ t_\mathrm{offset} > 0 $
		$ t_2 \LAR t_2 + t_\mathrm{offset} $
	else
		$ t_1 \LAR t_1 - t_\mathrm{offset} $
	return $ t_1 $, $ t_2 $
}

この処理で得た $ t_1 $, $ t_2 $ は三角形の1つ目と2つ目の頂点に対応する重心座標(barycentric coordinates)になっています。これより任意の三角形の頂点座標 $ \vec{v}_1 $, $ \vec{v}_2 $, $ \vec{v}_3 $ を用いてサンプル点 $ x $ を次のように構成できます。 \begin{equation*} \vec{x} = t_1 \vec{v}_1 + t_2 \vec{v}_2 + (1 - t_1 - t_2) \vec{v}_3 \end{equation*}

任意の多角形内の一様分布

任意の多角形は三角形の集合に分割することができます。多角形を構成する三角形の中から一つを面積に比例した確率で選択し、選択された三角形中で一様分布からのサンプリングを行うことで任意の多角形内の一様分布サンプリングを行うことができます。

図9. (a)平面上の三角形と球面三角形の対応、(b)球面三角形内の一様分布

少し難しいですが、球面に投影した三角形中で一様サンプリングする手法が[Arvo1995]によって提案されています。球面上の三角形のことを球面三角形(Spherical Triangle)と呼びます。サンプリング自体の紹介に入る前に3つの頂点 $ \vec{v}_a $, $ \vec{v}_b $, $ \vec{v}_c $ を持つ通常の三角形がどのようにして球面上の三角形に対応付けられるかについて説明します(図9(a))。球面三角形を構成する各頂点 $ \vec{A} $, $ \vec{B} $, $ \vec{C} $ は基準点 $ \vec{o} $ に対する相対的な位置ベクトルを正規化することで得られます。すなわち次の関係になります。 \begin{equation*} \vec{A} = \frac{\vec{v}_A - \vec{o}}{\norm{\vec{v}_A - \vec{o}}} \hspace{5mm} \vec{B} = \frac{\vec{v}_B - \vec{o}}{\norm{\vec{v}_B - \vec{o}}} \hspace{5mm} \vec{C} = \frac{\vec{v}_C - \vec{o}}{\norm{\vec{v}_C - \vec{o}}} \end{equation*} 球面上の頂点 $ \vec{A} $, $ \vec{B} $, $ \vec{C} $ それぞれに対応する"球面上"の内角 $ \alpha $, $ \beta $, $ \gamma $ に関するcosは次のように表されます。 \begin{equation*} \cos\alpha = -\vec{n}_{AB} \cdot \vec{n}_{CA} \hspace{5mm} \cos\beta = -\vec{n}_{BC} \cdot \vec{n}_{AB} \hspace{5mm} \cos\gamma = -\vec{n}_{CA} \cdot \vec{n}_{BC} \end{equation*} ここで $ \vec{n}_{AB} $, $ \vec{n}_{BC} $, $ \vec{n}_{CA} $ は頂点間の弧が含まれる面の(正規化された)法線ベクトルで次のように表されます。 \begin{equation*} \vec{n}_{AB} = \frac{\vec{A} \times \vec{B}}{\norm{\vec{A} \times \vec{B}}} \hspace{5mm} \vec{n}_{BC} = \frac{\vec{B} \times \vec{C}}{\norm{\vec{B} \times \vec{C}}} \hspace{5mm} \vec{n}_{CA} = \frac{\vec{C} \times \vec{A}}{\norm{\vec{C} \times \vec{A}}} \end{equation*} この球面三角形の面積、単位球に投影しているので同時に立体角 $ \Omega $ はジラール(Girard)の式より次のようになることが知られています。 \begin{equation*} \Omega = \alpha + \beta + \gamma - \pi \end{equation*} それでは必要な球面三角形の情報が揃ったところでサンプリング手順の紹介です。理論について興味がある方は論文を参照してください。アルゴリズムを簡単に示すために2つのベクトル間の演算 $ [\vec{x} | \vec{y}] $ を次のように定義しておきます。 \begin{equation*} [\vec{x} | \vec{y}] \equiv \frac{\vec{x} - (\vec{x} \cdot \vec{y}) \vec{y}}{\norm{\vec{x} - (\vec{x} \cdot \vec{y}) \vec{y}}} \end{equation*} 球面三角形中の点 $ \vec{P} $ をサンプリングするアルゴリズムを以下に示します。

function $ \mathrm{sampleSphericalTriangle} $ (単位矩形 $ [0, 1)^2 $ 中のサンプル $ u_1 $, $ u_2 $) {
	$ \hat{\Omega} = \Omega u_1 $
	$ s = \sin\Prt{\hat{\Omega} - \alpha} $
	$ t = \cos\Prt{\hat{\Omega} - \alpha} $
	$ u = t - \cos\alpha $
	$ v = s + \Prt{\vec{A} \cdot \vec{B}} \sin\alpha $
	$ \displaystyle q = \frac{(vt - us) \cos\alpha - v}{(vs + ut) \sin\alpha} $
	$ \vec{C}' = q \vec{A} + \sqrt{1 - q^2} [\vec{C} | \vec{A}] $
	$ z = 1 - u_2 (1 - \vec{C}' \cdot \vec{B}) $
	$ \vec{P} = z \vec{B} + \sqrt{1 - z^2} [\vec{C}' | \vec{B}] $
	return $ \vec{P} $
}

このアルゴリズムによってサンプリングした結果を図9(b)に示しています。このアルゴリズムは球面三角形の立体角中で一様サンプリングを実現するので、そのPDFは立体角の逆数、すなわち \begin{equation*} p_\sigma(\vomega) = \frac{1}{\Omega} \hspace{3mm} \Brk{\mathrm{sr}^{-1}} \end{equation*} となります。このPDFを元の三角形中の点 $ \vec{x} $ を得る、面積に関するPDFに変換すると次のようになります。 \begin{equation*} p_A(\vec{x}) = \frac{1}{\Omega} \frac{\absb{\vec{P} \cdot \vec{n}}}{\norm{\vec{x} - \vec{o}}^2} \hspace{3mm} \Brk{\mathrm{m}^{-2}} \end{equation*} ここで $ \vec{n} $ は元の平面三角形の(正規化)法線ベクトルを表しています。なお、$ \vec{x} $ は基準点 $ \vec{o} $ と方向 $ \vec{P} $ を用いてレイと三角形の交差判定から求めることができます。(交差は確実に起こることがわかっているのでいくつかの処理は省略できます。平面三角形と球面三角形の関係からもっと直接的に求める方法がありそうな気もする...)
ちなみにここでは球面三角形のサンプリングに関して紹介しましたが、球面長方形に関する同様のサンプリングも提案されています[Ureña2013]。2つの球面三角形を用いて球面長方形を表現するよりも微妙に特性が優れているようです。

レンダリングにおける利点

光源平面中の一様サンプリング 1spp	光源の各三角形を球面投影した立体角中で 一様サンプリング 1spp	リファレンス
図10. 光輸送問題への球面三角形サンプリング適用例

球面上のサンプリングも平面中のサンプリングと同じく、結局は光源上の点のサンプリングのために用いたりするのですが、モンテカルロ積分的に球面上のサンプリングは利点を持っています。点 $ \vx $ において $ \vz $ 方向に反射(透過)する放射輝度について考えてみます。ここでは簡単のために多重反射などは考えず、ある光源からの直接照明だけを考えます。解きたい積分は次で表されます(「レンダリング方程式2」参照)。 \begin{equation*} L_r(\vx \RAR \vz) = \int_{\mathcal{M}_L} f_s(\vy \RAR \vx \RAR \vz) L_e(\vy \RAR \vx) G(\vy \LRAR \vx) V(\vy \LRAR \vx) \d A(\vy) \end{equation*} ここで $ \vy $ は光源上の点、$ \mathcal{M}_L $ は光源の面を表しています。これをそのままモンテカルロ積分に変形すると推定量は次の形になります。(簡単のためサンプル数は1とします。) \begin{equation*} \Abk{L_r(\vx \RAR \vz)} = \frac{f_s(\vy \RAR \vx \RAR \vz) L_e(\vy \RAR \vx) G(\vy \LRAR \vx) V(\vy \LRAR \vx)}{p_A(\vy)} \end{equation*} 簡単のため光源の形状が単一の三角形だとすると平面中の一様サンプリングの場合の推定量は次の形になります。$ A_{\mathrm{tri}} $ は三角形の面積を表しています。 \begin{equation*} \Abk{L_r(\vx \RAR \vz)} = \frac{f_s(\vy \RAR \vx \RAR \vz) L_e(\vy \RAR \vx) G(\vy \LRAR \vx) V(\vy \LRAR \vx)}{\frac{1}{A_{\mathrm{tri}}}} = f_s(\vy \RAR \vx \RAR \vz) L_e(\vy \RAR \vx) G(\vy \LRAR \vx) V(\vy \LRAR \vx) A_{\mathrm{tri}} \end{equation*} 一方で、球面上の三角形中の一様サンプリングを行った場合の推定量は次の形になります。 \begin{equation*} \Abk{L_r(\vx \RAR \vz)} = \frac{f_s(\vy \RAR \vx \RAR \vz) L_e(\vy \RAR \vx) G(\vy \LRAR \vx) V(\vy \LRAR \vx)}{\frac{1}{\Omega_{\mathrm{tri}}(\vx)} \frac{\absb{\vomega_{\vx \vy} \cdot \vec{n}_{\vy}}}{\norm{\vy - \vx}^2}} = f_s(\vy \RAR \vx \RAR \vz) L_e(\vy \RAR \vx) \absb{\vomega_{\vx \vy} \cdot \vec{n}_{\vx}} V(\vy \LRAR \vx) \Omega_{\mathrm{tri}}(\vx) \end{equation*} $ \Omega_{\mathrm{tri}}(\vx) $ は球面上に投影した三角形の立体角を表しています。球面三角形の立体角は基準となる点 $ \vx $ によって変わるので $ \vx $ の関数となっています。2つのサンプリングによるそれぞれの推定量を見比べてみると後者では $ G $ 項が部分的にキャンセルされていることがわかります。ここで重要なのが $ G $ 項の分母に含まれる距離の2乗の項がキャンセルされていることです。シェーディング点 $ \vx $ が光源の面に非常に近い場合、光源面上のどこをサンプリングするかによって $ G $ 項の値は非常にばらつきます。すなわち前者のモンテカルロ推定量の分散が大きくなる可能性があります。一方で後者では距離の2乗項がキャンセルされているので相対的に分散は小さくなります。図10にそれぞれのサンプリング方法を用いてレンダリングした例を示します。このシーンでは床に対して垂直な光源が奥に見えますが、平面中の一様サンプリングの場合は光源が接地している箇所の周辺の分散が大きくなっているのに対し、球面サンプリングの場合は分散が画面全体で比較的一様になっていることがわかります。ただし球面サンプリングの場合でも光源を構成する2つの三角形の選択は一様に行っているため、接地箇所の分散は依然わずかに大きくなっています。

光源の放射発散度(や配光特性)が一様なケースでは球面三角形によるサンプリングのほうが、サンプリングコストは高くとも生成されるサンプルの質が高いため総合的には高い効率となるでしょう。光源がテクスチャーによって制御される放射発散度の分布を持っている場合、単純な平面中のサンプリングだとテクスチャーの値に応じた重点的サンプリングを簡単に行えますが、球面三角形だと難しいかもしれません。そのような場合には平面中のサンプリングのほうが優れる可能性があります。