レンダリングにおけるimportance samplingの基礎

Importance Samplingは、モンテカルロ積分の期待値に近い値を、より少ないサンプリング数で得ようとする手法です。ここでは特に、レンダリングに関連した取り扱いについて重点的に記していきたいと思います。
2016-09-19 InlineのLaTexを入れ替え。ついでに編集。

モンテカルロ積分: Monte Carlo Integraion

背後にある基本的なアイデア

モンテカルロ積分の背後にある基本的なアイデアは、xの関数f(x)の領域Ωでの積分を、ランダムなサンプリングXi=(x1,x2,…xN)を用いて近似することです。
\begin{aligned} I &=\int_\Omega f(x) \textrm{d}x \\ \hat{I} &= \frac{\Omega}{N}\sum_{i=1}^{N}{f(Xi)} \end{aligned}

上記の\hat{I}のことをMonte Carlo Estimatorと呼ぶそうです。なぜこの式で近似できるかは直感的に分かる気がしますが、まずは、これの期待値を考えます。
\begin{aligned} E[\hat{I}] &= E[\frac{\Omega}{N}\sum_{i=1}^{N}{f(Xi)}] \\ &= \frac{\Omega}{N} \sum_{i=1}^{N}{E[f(Xi)]} \\ &= \Omega E[f(Xi)] \end{aligned}

ここでE[f(Xi)]は領域Ω内での離散的なサンプリングXiの関数f(Xi)の期待値です。これを連続な関数f(x)を用いて書き直すならば、f(x)に、xの確率密度関数:probability density function:PDF p(x)を乗算して積分することになります。
\begin{aligned} \Omega E[f(x)] = \Omega \int_{\Omega}f(x)p(x) \textrm{d}x \end{aligned}

例えばPDFが一定、つまりXiの分布が一様であると仮定すると、以下の様になります。
\begin{aligned} p(x) &= \frac{1}{\Omega} \\ E[\hat{I}] &= \int_{\Omega}f(x)\textrm{d}x = I \end{aligned}
つまりMonte Carlo Estimatorの期待値を求めることで、元来求めたかった積分の解である、Iを求めることができるわけです。

積分したい関数f(x)について

f(x)について何も前提条件がなければ、上記のように一様なサンプリングを行うしかありません。ただ、f(Xi)の分散が大きい場合、得られる解の分散も大きくなります。つまり近似した積分の値の誤差が大きくなる傾向にあります。
もし、関数f(x)に関して何らかの前提条件、たとえば、領域Ω内のある領域では値が大きくなり、またある領域では値が0に近くなるなどが分かれば、値が大きくなるところを重点的にサンプリングすることで、f(Xi) の分散が小さくなるようにサンプリング出来るでしょう。ただし、サンプリングの分布は領域Ωにおいて一様ではなくなり、サンプリングのPDF p(x)は定数ではなくなるため、上記の式は成り立ちません。

ここで、f(x)を以下のように分解した式のg(x)の積分について考えます。
\begin{aligned} f(x) = g(x)p(x) \end{aligned}

\begin{aligned} I_g &=\int_\Omega g(x) \textrm{d}x \\ \hat{I_g} &= \frac{\Omega}{N}\sum_{i=1}^{N}{g(Xi)} \end{aligned}

\hat{I_g}の期待値は、以下のようになります。
\begin{aligned} E[\hat{I_g}] &= \Omega E[g(x)] &= \Omega \int_{\Omega} f(x) \textrm{d}x \end{aligned}

これを、元のIの式に代入すると、以下のようになります。
\begin{aligned} I&=\int_{\Omega} f(x) \textrm{d}x = \frac{1}{\Omega}E[\hat{I_g}] &= E[\frac{1}{N} \sum_{i=1}^{N} \frac{f(Xi)}{p(Xi)}] \end{aligned}

つまり、サンプリングの確率密度を考慮したMonte Carlo Estimatorを定義すると、以下のようになります。
\begin{aligned} \hat{I} = \frac{1}{N} \sum_{i=1}^{N} \frac{f(Xi)}{p(Xi)} \end{aligned}

Lambert Diffuseを、一様ランダムベクトルで積分してみる

まずはじめに、Lambert Diffuseモデルの反射を考えてみたいと思います。
\begin{aligned} L_o = \int_{\Omega} \frac{kd}{\pi} L_i \cos\theta \textrm{d}\omega \end{aligned}
θは微小面の法線と入射光の成す角です。 Lambert DiffuseのBRDFはkd/πで一定です。Liが入射光のradianceで、これを見かけ上の微小面積に投影してBRDFを乗算すればLo : radianceが求まると思います。
さて、先のモンテカルロ積分の式に上記積分を当てはめていきます。まずLiのサンプリングを決めます。サンプリングは半球の全領域で一様とします。したがって、PDFは以下のように一定になります。
\begin{aligned} p(\omega) = \frac{1}{2\pi} \end{aligned}

一方f(x)は、以下のようになります。
\begin{aligned} f(x) = \frac{kd}{\pi} Li \cos\theta \end{aligned}

Loの近似値は以下のように表すことができます。
\begin{aligned} \hat{L_o} = \frac{2kd}{N} \sum_{i=1}^{N} L_{ii} \cos\theta_i \end{aligned}

半球上に分布する一様ランダムベクトルの生成

モンテカルロ積分は、式内にサンプリングのPDFがあるので、上記のモンテカルロ積分の式は、下記のようなサンプリングの関数と対になって使用しなければなりません。今回のサンプリングは一様ランダムベクトルなのでベクトルの生成方法はいろいろ考えられますが、一例として下記のように行うこととします。使用するのは2つの0~1の一様乱数, u1, u2です。これを用いて半球状に一様なベクトルを生成します。z軸を天頂方向（つまり0~1）としています。

vec3 SampleUniform(float u1, float u2)
{
  float r = sqrt(1.0f - u1 * u1);
  float phi = 2 * PI * u2;
  return vec3(r * cos(phi) , r * sin(phi), u1);
}

半球上に分布する一様ランダムベクトルの生成(2)

上記のサンプリング関数をPDFから求めることも出来ます。逆変換法:inversion methodと呼ばれる方法です。
まず、PDFは、p(ω) = 1/2πで一定ですが、具体的なベクトルのPDFをを求めるため、極座標系の関数p(θ,φ)が必要になります。単位球の微小面積は以下のように求まるので、これを利用し、p(θ,φ)を求めます。
\begin{aligned} \textrm{d}\omega = \sin\theta \textrm{d}\theta \textrm{d}\phi \end{aligned}
\begin{aligned} p(\theta, \phi) = \frac{\sin\theta}{2\pi} \end{aligned}
次に、θの周辺密度関数:marginal density function:MDFを求めます。
\begin{aligned} p(\theta) = \int_{0}^{2\pi}p(\theta, \phi)\textrm{d}\phi = \sin\theta \end{aligned}
次に、φの条件付密度関数:conditional density function p(φ|θ)を求めます。
\begin{aligned} p(\phi | \theta) = \frac{p(\theta, \phi)}{p(\theta)} = \frac{1}{2\pi} \end{aligned}
ここで、p(θ)とp(φ|θ)の累積分布関数:cumulative distribution function:CDF P(θ)とP(φ|θ)を求めます。
\begin{aligned} P(\theta) &= \int_{0}^{\theta} \sin\theta' \textrm{d}\theta' = -\cos\theta + 1\\ P(\phi | \theta) &= \int_{0}^{\phi}\frac{1}{2\pi}\textrm{d}\phi' = \frac{\phi}{2\pi} \end{aligned}

さて。ここで上記のCDFについて考えてみます。CDFはサンプリングの領域において0~1へと増加していきます。たとえば、P(θ)は、θが0~π/2の値をとるのにしたがって、0~1へと増加していきます。この増加の傾きは、θのPDFによって与えられます。
ここで、CDFの逆関数を考えてみます。P(θ)をξ1とした場合、P(θ)の逆関数は以下のようになります。
\begin{aligned} \theta = \arccos(1-\xi_1) \end{aligned}
この関数は、ξ1が0~1の値をとるのにしたがって増加してゆき、θは0~π/2の値をとります。この関数の傾きはθのPDFの逆数になるはずです。ここで、ξ1に0~1の一様乱数を与えたときに得られるθの分布は、θのPDFに従うことになります。
同様にP(φ|θ)の逆関数は以下のようになり、こちらもφのMDFの分布に従います。
\begin{aligned} \phi = 2\pi \xi_2 \end{aligned}

このことを利用して、2つの0~1の一様乱数より、サンプリングベクトルの生成を行うことが出来ます。

vec3 SampleUniform(float u1, float u2)
{
  float theta = arccos(1-u1);
  float phi = 2 * pi * u2;

  return vec3(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
}

少し計算を整理すれば、以下のようになります。章のはじめに書いた関数と、u1の向きが逆になっているだけで、同様の関数が得られました。

vec3 SampleUniform(float u1, float u2)
{
  z = 1 - u1
  r = sqrt(1- z*z);
  float phi = 2 * pi * u2;

  return vec3(r * cos(phi), r * sin(phi),  z);
}

Lambert Diffuseをcosine weighted samplingで積分してみる

Lambert Diffuseモデルにおけるサンプリングベクトルの分布について

ここで、Lambert Diffueseをモンテカルロ積分をする上で、もっと効率の良いサンプリングベクトルの分布について考えてみたいと思います。簡単なことですが、cos(θ)が0に近い領域ではLiの値にかかわらずLoの値が小さくなります。ということは、その領域でLiの不正確なサンプリングをしても、全体の誤差に与える影響が小さいということになります。言い換えれば、サンプルの数を減らしても全体の誤差に与える影響が小さいということです。一方でcos(θ)が大きな値になる領域では、全体の誤差に与える影響が大きいので、なるべく正確なサンプリングをしたほうが良好な結果が得られるはずです。つまり、積分したい関数の値とサンプリングのPDFは比例関係であるのが望ましいということです。
\begin{aligned} f(x) \propto p(x) \end{aligned}

ただし、この条件は一般的に非常に困難です。なぜなら、PDFは領域全体で積分した時に1になる必要があり、上式の比例定数を求めるためにはf(x)の積分を求める必要があります。
\begin{aligned} p(x) = \frac{1}{\int_{\Omega}f(x)\textrm{d}x}f(x) \end{aligned}
この式内の積分を求めることは、モンテカルロ積分で求めたいものを求めることになり本末転倒になるというわけです。したがって、通常p(x)はf(x)の係数の一部であったり、f(x)をテイラー展開をして低次の項を用いたりするようです。

cosine weighted sampling

上記をふまえると、LambertDiffueseをモンテカルロ積分するならば、PDF:p(ω)はcos(θ)に比例するのが好ましいといえます。PDFは領域全体で1にならなくてはならないことを用いて比例定数を求めます。
\begin{aligned} 1 &= c \int_{\Omega} p(\omega) \textrm{d}\omega \\ &=c \int_{\phi=0}^{2\pi} \int_{\theta=0}^{\frac{1}{2\pi}} \cos(\theta)\sin(\theta) \textrm{d}\theta \textrm{d}\phi \\ &= c 2 \pi \int_{0}^{\frac{1}{2\pi}} \cos(\theta)\sin(\theta) \textrm{d}\theta \\ &= c \pi \int_{0}^{\frac{1}{2\pi}} \sin(2\theta) \textrm{d}\theta = \frac{c}{\pi} \\ c &= \frac{1}{\pi} \end{aligned}

したがって、PDFは以下のようになります。
\begin{aligned} p(\theta, \phi) = \frac{1}{\pi}\sin(\theta)\cos(\theta) \end{aligned}

上記のPDFにinversion methodを適用して、サンプリング用の関数を作ります。
θのMDFは、以下のようになります。
\begin{aligned} p(\theta) = \int_{0}^{2\pi}p(\theta, \phi) \textrm{d}\phi = 2 \sin\theta\cos\theta \end{aligned}

一方、φのconditional density functionは以下のようになります。
\begin{aligned} p(\phi|\theta) = \frac{p(\theta, \phi)}{p(\theta)} = \frac{1}{2\pi} \end{aligned}

CDFはそれぞれ以下のようになります。
\begin{aligned} P(\theta) &= 2 \int_{0}^{\theta} \sin\theta'\cos\theta'\textrm{d}\theta' \\ &= \int_{0}^{\theta} \sin(2\theta')\textrm{d}\theta' \\ &= \frac{1}{2}(\cos(2\theta)-1) \\ &= \sin2\theta \\ P(\phi|\theta) &= \frac{\phi}{2\pi} \end{aligned}

CDFの逆関数はそれぞれ以下のようになります。
\begin{aligned} \theta &= \arcsin\sqrt{\xi_1}\\ \phi &= \frac{2}{\pi} \xi_2 \end{aligned}

上記をふまえ、サンプリング関数は以下のようになります。

vec3 SampleCosineWeighted(float u1, float u2)
{
  float theta = asin(sqrt(u1));
  float phi = 2 * pi * u2;

  return vec3(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
}

上記を少し整理すると、以下のようになります。

vec3 SampleCosineWeighted(float u1, float u2)
{
  float r = sqrt(u1);
  float phi = 2 * pi * u2;

  return vec3(r * cos(phi), r * sin(phi), sqrt(1-u1));
}

一方で、モンテカルロ積分の式は、下記の通りとなります。
\begin{aligned} p(\omega) = \frac{cos(\theta)}{\pi} \end{aligned}
\begin{aligned} \hat{L_o} &= \frac{kd}{\pi}\frac{1}{N} \sum_{i=1}^{N} \frac{L_{ii}\cos\theta_i}{\frac{cos(\theta_i)}{\pi}} \\ &= \frac{kd}{N} \sum_{i=1}^{N} L_{ii} \end{aligned}