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

確率密度関数 (PDF) の変換

このページはレイトレ合宿4!?のアドベントカレンダー2週目の記事として制作しました。
みなさんこんにちは!shockerです。モンテカルロレイトレーシングにおいてPDFという概念は極めて重要です。分散低減のための様々なテクニックにはそれぞれのPDFが存在しますが、異なる単位を持つことが多いので単位を正しく取り扱わないと間違った結果になってしまいます。このページではそれらPDFの変換方法について解説します。

モンテカルロレイトレーシングでは様々な戦略によって確率的に経路をサンプリングしますが、それぞれの戦略で素直に得られる確率密度関数は異なる単位を持っていることも多いので、最終的には単位を揃える必要があります。ここでは確率密度関数の変換方法について解説します。

変数変換とヤコビアン

直交座標系における微小面積の変化
(a) 直交座標系における微小面積の変化
極座標系における微小面積の変化
(b) 極座標系における微小面積の変化
図1. 微小面積の変化
(a)直交座標系においては座標によらず微小面積の大きさは一定だが、(b)極座標系では半径方向で微小面積の大きさが変化する。

PDFの単位の変換には重積分などを勉強する際に登場するヤコビアンが密接に関わってきます。

例として2次元の積分を考えます。図1の左に示すように、直交座標系の場合、微小な変位 $ dx, dy $ によって作られる微小面積は座標中のどこで考えても同じ大きさになります。そのため直交座標では被積分関数だけを積分すれば求めたい積分が求まります。一方で、図1の右に示すような、極座標系の場合、微小な変位 $ dr, d\theta $ によって作られる微小面積は座標によって大きく異なってきます。具体的には2次元直交座標系の場合、微小面積は $ r $ に比例して変化します($ \theta $ には依存しません)。そのため、2次元直交座標系で表された関数を2次元極座標基準で積分する場合には、この微小面積の変化を考慮する必要があり、以下の関係のように表されます。

\begin{equation*} \iint_{A} f(x, y) dx dy = \iint_{A_p} f(r, \theta) r dr d\theta \end{equation*}

極座標系の式では、変数変換された被積分関数に加えて $ r $ が追加でかかっていることがわかります。積分範囲自体は変わりません (座標系が異なるので積分範囲を表す式は変わります)。
こういった変化量を考慮に入れなければならない理由の説明として、例えばある国の人口密度を考えてみます。各県や州ごとの人口密度のデータを既に持っているとして、国の人口密度をそこから求めるには、単なる平均では精度よく求まりません。極めて高い人口密度を持つ場所があったとしても、そこが極めて面積が小さい場所なら、全体の密度に対する影響は大きくなくなるからです。各県や州それぞれの面積に応じたウェイトをかけて平均をとる必要があります。2次元極座標の場合には外側($ r $ が大きい)ほど微小面積が大きくなるので、それに応じたウェイトがかかることなります。(積分の微小面積は無限に小さいですが、無限小でも無限小なりの大きさの違いがあるのです。)

ここでは2次元の直交座標と極座標の積分の関係を表しましたが、この変数変換を一般化すると次のように表されます。

\begin{equation*} \int \cdots \int_D f(x_1, \ldots, x_n) dx_1, \ldots, dx_n = \int \cdots \int_{D'} f(\theta_1, \ldots, \theta_n) \abs{\mathrm{det}(J_{\Vec{x} \, \RAR \, \Vec{\theta}})} d\theta_1, \ldots, d\theta_n \end{equation*} \begin{equation*} J_{\Vec{x} \, \RAR \, \Vec{\theta}} = \begin{pmatrix} \displaystyle\frac{dx_1}{d\theta_1} & \ldots & \displaystyle\frac{dx_1}{d\theta_n} \\ \vdots & \ddots & \vdots \\ \displaystyle\frac{dx_n}{d\theta_1} & \ldots & \displaystyle\frac{dx_n}{d\theta_n} \\ \end{pmatrix} \end{equation*}

ここで、$ J_{\Vec{x} \, \RAR \, \Vec{\theta}} $ が所謂ヤコビアン(Jacobian, ヤコビ行列)と呼ばれるものです。積分の変数変換ではヤコビアンの行列式の絶対値を使います。ちなみにヤコビアンの行列式を単にヤコビアンと呼んだりもするので、このへんは文脈次第でしょう。

ヤコビアンは変数変換の際に必要となることが多いので、変数変換という操作と結びつけがちですが、最初から極座標で表された関数をそのまま極座標で積分したい、という際にも必要なので気をつけましょう。

PDFの単位の変換

PDFの単位の変換は結局変数変換と同じ考え方になります。ただし、PDFの単位は積分変数と逆数の関係にあるので、ヤコビアンの向きは逆向きに考えます。

\begin{equation*} p_{\Vec{\theta}}(\theta_1, \ldots, \theta_n) = \frac{1}{\abs{\mathrm{det}(J_{\Vec{\theta} \, \RAR \, \Vec{x}})}} p_{\Vec{x}}(x_1, ... x_n) \end{equation*}

つまり変数 $ \Vec{x} $ の関数であるPDF $ p_{\Vec{x}} $ を変数 $ \Vec{\theta} $ で表されるPDFに変換したい場合は $ J_{\Vec{x} \, \RAR \, \Vec{\theta}} $ ではなく、$ J_{\Vec{\theta} \, \RAR \, \Vec{x}} $ を考えます。

変換の例

IBLのサンプリングにおけるuv座標と立体角

IBLによる無限遠光源のサンプリング時には、多くの場合まずは2次元の画像上でサンプリングを行いますが、レンダリング方程式が立体角(や物体表面の面積)に関する積分になっているので、画像上の確率密度から立体角に関する確率密度に変換する必要があります。サンプリングを順序にそって説明すると、まず2次元の$ [0, 1)^2 $のテクスチャー座標をPDF $ p_{tc}(u, v) $ に沿ってサンプリング、テクスチャー座標を球面座標の$ (\theta, \phi) $に変換、球面座標を方向ベクトル$ (x, y, z) $に変換します。一連の流れは以下のようになります。

\begin{equation*} (u, v) \xrightarrow[\displaystyle J_{\Vec{\theta} \, \RAR \, tc}]{} (\phi, \theta) \xrightarrow[\displaystyle J_{\Vec{x} \, \RAR \, \Vec{\theta}}]{} (x, y, z) \end{equation*}

前述のPDFの変換式を当てはめると以下の関係を得ます。

\begin{eqnarray*} p_{\Vec{\theta}}(\theta, \phi) &=& \frac{1}{\abs{\mathrm{det}(J_{\Vec{\theta} \, \RAR \, tc})}} p_{tc}(u, v) \\ p_\sigma(\vomega) = p_{\Vec{x}}(x, y, z) &=& \frac{1}{\abs{\mathrm{det}(J_{\Vec{x} \, \RAR \, \Vec{\theta}})}} p_{\Vec{\theta}}(\theta, \phi) \\ &=& \frac{1}{\abs{\mathrm{det}(J_{\Vec{x} \, \RAR \, \Vec{\theta}})}} \frac{1}{\abs{\mathrm{det}(J_{\Vec{\theta} \, \RAR \, tc})}} p_{tc}(u, v) \end{eqnarray*}

$ (\phi, \theta) \rightarrow (x, y, z) $ で次元が変わってしまっているように見えますが、実際には単位球 $ x^2 + y^2 + z^2 = 1 $ もしくは $ r = 1 $ という制限があるので2次元のままです。ただし、ヤコビアンに関しては $ x = r \cos\phi \sin\theta,\, y = r \sin\phi \sin\theta,\, z = r \cos\theta $ という3次元の関係を考える必要があります。

それぞれのヤコビアンの行列式の絶対値は次のように計算できます。

\begin{equation*} \abs{\mathrm{det}(J_{\Vec{\theta} \, \RAR \, tc})} = \abs{\begin{vmatrix} 2\pi & 0 \\ 0 & \pi \\ \end{vmatrix}} = 2\pi^2 \end{equation*} \begin{equation*} \abs{\mathrm{det}(J_{\Vec{x} \, \RAR \, \Vec{\theta}})} = \abs{\begin{vmatrix} \cos\phi \sin\theta & r \cos\phi \cos\theta & -r \sin\phi \sin\theta \\ \sin\phi \sin\theta & r \sin\phi \cos\theta & r \cos\phi \sin\theta \\ \cos\theta & -r \sin\theta & 0 \\ \end{vmatrix}} = r^2 \sin\theta \end{equation*}

したがってIBLのある方向からサンプリングするPDF $ p_\sigma(\vomega) $ は次のように表すことができます。

\begin{equation*} p_\sigma(\vomega) = \frac{1}{2\pi^2 \sin\theta} p_{tc}(u, v) \end{equation*}

参考文献

  • [Pharr2010] Matt Pharr, Greg Humphreys - "PHYSICALLY BASED RENDERING: From Theory to Implementation Second Edition", 2010, Morgan Kaufmann

Today : 00000071 Total : 00061010
Copyright © 2017 shocker.