last modified: 2021-12-09 03:21:16
このページではMathJaxによる数式表示を行っています。対応しているブラウザーであれば少し待てば数式が表示されます。
(画面左下にMathJaxのタイプセットステータスが表示されます。)
環境によっては表示が崩れているかもしれません。その際はブラウザーを変更してみてください。

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

このページはレイトレ Advent Calendar 2021の記事として制作しました。皆さんこんにちは、@Shocker_0x15です。このウェブサイトでまとめてる情報もパストレーシングといった基礎に至るまではだいぶ埋まったよなぁ、何書こうかなぁと思ってたらそういえばサンプリング周り書いてなかったなと気が付きました。

このページでは逆関数法の理論とレンダラーの実装にあたってよく使用するジオメトリのサンプリング実装例について紹介します。

逆関数法

ある関数 $ f(x) $ に比例した確率密度分布でサンプリングを行いたい場合、逆関数法(Inversion Method, Inverse Transform Method)が使える場合には解析的なサンプリングを行うことができます。逆関数法では関数 $ f(x) $ に比例した確率密度関数 (PDF) $ p(x) $ の累積分布関数 (CDF) $ P(X \leq x) $ を求め、$ P(x) $ の逆関数 $ P^{-1}(u) $ に $ [0, 1) $ の一様乱数を入力することで $ p(x) $ に沿った、つまり $ f(x) $ に比例した分布でサンプルを得ることができます。

一次元関数の例

目標分布(1D)
図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) $ に沿ったサンプリングが実現できていることがわかります。

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

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

CDFによっては逆関数が解析的に計算できない場合があります。その場合でもニュートン法などの球根アルゴリズムによって実用上問題ないレベルでのサンプリングを行うことができますが、どうしても処理としては重くなりがちです。逆関数法とは異なるアプローチで解析的なサンプリングを行える手法Triangle-cut Parameterization[Heitz2020]も提案されています。

逆関数法は2次元以上の関数にも適用することができます。サンプル分布を比例させたい2次元関数の目標分布 $ f(x, y) $ を考えます。まずは1次元の場合と同様に積分を経て比例係数を求めることによって(同時)確率密度関数 $ p(x, y) $ を求めます。このとき $ p(x, y) $ が各次元に綺麗に分離できる、すなわち $ p(x, y) = p_x(x) p_y(y) $ として表せるかつ、それぞれの次元の値の範囲が独立している場合は、上に示した1次元の手順を用いて各次元のサンプルを「個別に」得るだけです。一般的な分離できない場合は少しだけ導出が面倒になります。その場合の手順を次に示します。

まずは1つ目の次元 $ x $ に関する周辺確率密度関数 $ p(x) $ を求めます。これは $ x $ 以外の次元全てに関して積分を行うことで求めることができます。 \begin{equation} p(x) = \int p(x, y) \d y \end{equation} $ x $ に関する周辺確率密度関数 $ p(x) $ が求まれば、ある $ x $ が与えられた場合の $ y $ の条件付き確率密度関数 $ p(y | x) $ を求めることができます。 \begin{equation} p(y | x) = \frac{p(x, y)}{p(x)} \end{equation} まずは $ p(x) $ に関して通常の1次元の逆関数サンプリングを行います。そして得られた $ x $ を用いれば $ p(y|x) $ もまた $ y $ に関する1次元の関数になるので通常の逆関数サンプリングを行います。3次元以上の場合も同じ考え方を再帰的に適用することで逆関数法を使用することができます(解析的に逆関数が計算できるかは別として)。

二次元関数の例

2Dの同時確率密度関数と達成された分布
図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分布

球面内一様分布

球面上の一様分布のPDFは次のように全球の立体角 $ 4\pi $ の逆数として定数で表されます。 \begin{equation*} p_\sigma(\vomega) = \frac{1}{4\pi} \hspace{3mm} \Brk{\mathrm{sr}^{-1}} \end{equation*} これを上記で示したように逆関数法でサンプリング方法を求めると $ \theta $, $ \phi $ はそれぞれ次のように求められます。 \begin{eqnarray*} \theta &=& \cos^{-1} \Prt{1 - 2 u_1} \\ \phi &=& 2 \pi u_2 \end{eqnarray*} $ u_1 $, $ u_2 $ はいずれも $ [0, 1) $ の一様乱数です。これによって達成された分布を図4(a)に示しています。3次元の図なのでいまいち把握しづらいかもしれませんが、球面に一様に分布していそうなことはわかるでしょうか。

半球面内一様分布

半球面上の一様分布のPDFは全球と同様に半球の立体角 $ 2\pi $ の逆数として定数で表されます。 \begin{equation*} p_\sigma(\vomega) = \frac{1}{2\pi} \hspace{3mm} \Brk{\mathrm{sr}^{-1}} \end{equation*} サンプリングは次のように表されます。 \begin{eqnarray*} \theta &=& \cos^{-1} \Prt{1 - u_1} \\ \phi &=& 2 \pi u_2 \end{eqnarray*} $ 1 - u_1 $ は実質 $ u_1 $ として表しても問題ありません。これによって達成された分布を図4(b)に示しています。

半球面内cos分布

半球内で天頂角 $ \theta $ が小さいほど密になるような分布、具体的には $ \cos\theta $ に比例する場合の分布をcos分布と呼んだりします。cos分布のPDFは次のように単純にcosに比例するはずです。 \begin{equation*} p_\sigma(\vomega) \propto \cos\theta \end{equation*} せっかくなので天頂角周りへの密具合を指数 $ n $ で制御できるような一般的なcos分布を考えてみます。 \begin{equation*} p_\sigma(\vomega) \propto \cos^{n} \theta \end{equation*} PDFを求めると次のようになります。 \begin{equation*} p_\sigma(\vomega) = \frac{n + 1}{2\pi} \cos^{n} \theta \hspace{3mm} \Brk{\mathrm{sr}^{-1}} \end{equation*} サンプリングを導出すると次の形を得ます。 \begin{eqnarray*} \theta &=& \cos^{-1} \Prt{1 - u_1}^{\frac{1}{n + 1}} \\ \phi &=& 2 \pi u_2 \end{eqnarray*} $ n = 0 $ の場合には半球面上の一様分布と等しくなります。通常のcos分布($ n = 1 $)によって達成された分布を図4(c)に示しています。図4(b)の半球一様分布と比較して天頂角周りにサンプルが集中していることがわかります。cos分布は完全拡散面(ランバート面)におけるBRDFとcos項の積に対する重点的サンプリングなどに用いられます。また指数 $ n $ 付きの一般的な分布に関しても適当な方向に傾けることで簡単・軽負荷な光沢面のサンプリング実装などに使えるかもしれません。

平面内のサンプリング

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

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

平行四辺形内一様分布
図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にそれぞれのサンプリング方法を用いてレンダリングした例を示します。このシーンでは床に対して垂直な光源が奥に見えますが、平面中の一様サンプリングの場合は光源が接地している箇所の周辺の分散が大きくなっているのに対し、球面サンプリングの場合は分散が画面全体で比較的一様になっていることがわかります。
光源の放射発散度(や配光特性)が一様なケースでは球面三角形によるサンプリングのほうが、サンプリングコストは高くとも生成されるサンプルの質が高いため総合的には高い効率となるでしょう。光源がテクスチャーによって制御される放射発散度の分布を持っている場合、単純な平面中のサンプリングだとテクスチャーの値に応じた重点的サンプリングを簡単に行えますが、球面三角形だと難しいかもしれません。そのような場合には平面中のサンプリングのほうが優れる可能性があります。

参考文献

  • [Arvo1995] James Arvo - "Stratified Sampling of Spherical Triangles", 1995
  • [Heitz2019] Eric Heitz - "A Low-Distortion Map Between Triangle and Square", 2019
  • [Heitz2020] Eric Heitz - "Can't Invert the CDF? The Triangle-Cut Parameterization of the Region under the Curve", 2020
  • [Pharr2016] Matt Pharr, Wenzel Jakob, Greg Humphreys - "PHYSICALLY BASED RENDERING: From Theory to Implementation Third Edition", 2016, Morgan Kaufmann
  • [Shirley1997] Peter Shirley, Kenneth Chiu - "A Low Distortion Map Between Disk and Square", 1997
  • [Ureña2013] Carlos Ureña, Marcos Fajardo, Alan King - "An Area-Preserving Parameterization for Spherical Rectangles", 2013

Today : 00000025 Total : 00209109
Copyright © 2022 shocker.