レンダリングにおけるImportance Samplingの基礎(2)

今回も前回に引き続きImportanceSamplingについてです。

PhongのBRDF

PhongのBRDFは、視線の反射ベクトルRとL_iの関数になっています。残りは正規化項(n+1)/(2π）やSpecular係数のk_sといった要素になっています。

L_o = \int_{\Omega} \frac{(n+1)k_sL_i}{2\pi} (\vec{R} \cdot \vec{L_i})^{n} \textrm{d}\omega

今回も極座標系を使いますが、反射ベクトルRを天頂とした極座標を用いて、以下のように式を書きます。

L_o = \frac{(n+1)k_s}{2\pi} \int_{\phi=0}^{2\pi}\int_{\theta=0}^{\frac{\pi}{2}} L_i \cos^{n}\theta \sin\theta \textrm{d}\theta \textrm{d}\phi

上記の式をモンテカルロ積分するのですが、PDFには正規化Phongの部分をそのまま使えば半球で積分すれば1になるので都合が良いです。なぜ下記の式を積分すると1になるかは、以前に書いた記事「Phongの正規化」を見ていただければ分かると思います。

p(\theta,\phi)=\frac{n+1}{2\pi}\cos^{n}\theta\sin\theta

あとは、前回の記事でも使用した、inversion methodを使って式を導出していきます。MDFとconditional density functionはそれぞれ以下のようになります。

p(\theta) = (n+1)\cos^{n}\theta \sin\theta
p(\phi|\theta) = \frac{1}{2\pi}

さらに、それぞれのCDFは以下のようになります。

P(\theta) = (n+1)\int_{0}^{\theta}\cos^{n}\theta'\sin\theta' \textrm{d}\theta

昔の記事でPhongの正規化を行った際と同様の漸化式と部分積分を使って積分します。

\\ \begin{aligned} I_n &= \int_{0}^{\theta}\cos^{n}\theta'\sin\theta' \textrm{d}\theta \\ f &= \cos^{n}\theta \\ g' &= \sin\theta \\ I_n &= (1 - \cos^{n+1}\theta) - nI_n \\ I_n &= \frac{1 - \cos^{n+1}\theta}{n+1} \\ \end{aligned}
P(\theta) = 1-\cos^{n+1}\theta

φの方は自明なので、導出は省略です。

P(\phi | \theta) = \frac{\phi}{2\pi}

P(θ)をξ1とし、P(φ|θ)をξ2として、逆関数を求めます。
\begin{aligned} \theta &= \arccos((1-\xi_1)^{\frac{1}{n+1}}) \\ \phi &= 2\pi\xi_2 \end{aligned}

これを用いたサンプリング用の関数は以下のようになります。
2016-08-21修正。
上記関数で求まるθ,φは、反射ベクトルの摂動です。Samplingの関数は接空間を基準にしています。反射ベクトルを求めるのは簡単ですが、反射ベクトルは単なるベクトルで、空間の概念が無いので、そのままでは摂動させることができません。サンプリングのベクトルを求めるには、反射ベクトル空間が必要になります。下記の例では、入射ベクトルと反射ベクトルの外積のベクトルを補助的に使って、任意軸の回転を用いて、反射ベクトルを摂動させます。特異点がある上に、計算コストが高いのであまり良い方法とは言えないと思います。

vec3 SamplePhong(float u1, float u2, vec3 V)
{
  float theta = acos(pow(1-u1, 1/(n+1)));
  float phi = 2 * pi * u2;

  vec3 R = vec3(-V.x, -V.y, V.z);
  vec3 Ro = cross(V, R);
　if (length(Ro) < epsilon) {
    // in the case when V and R are aligned with Z axis.
    return vec3(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta))
  }
  Ro = normalize(Ro);
  vec3 Rphi = AxisRotation(R, Ro, phi);
  vec3 Rtheta = AxisRotation(Rphi, R, theta);
  return Rtheta;
}

一方でモンテカルロ積分の式は、以下のようになります。

\hat{L_o}=\frac{k_s}{N}\sum_{1}^{N} L_{ii}

ここで一点気をつけなくてはならない点があります。サンプリング用のベクトルを返す関数は、視線ベクトルを反射させたものを摂動させているので、法線との角度が90度以上のベクトルを返すことがあります。このベクトルの場合はPhongの値は明らかにゼロですが、モンテカルロ積分の式に視線ベクトルや法線の関数が残っていないので、Liを加算するときには、Liのベクトルを別途チェックする必要があります。

話は逸れて、正規分布とBox Muller transformについて

（レンダリングと一旦離れますが最終的な目標はexp(x)が含まれるBRDFをImportance Samplingすることにあります。ただし、今回はそこまで到達しません。）
ここで、Importance Samplingの確率密度関数:PDFに正規分布が使えないかを考えてみます。 標準正規分布の式は以下のとおりで、全区間の積分は１であることが知られています。したがって、この関数はPDFとしてそのまま使用可能です。

f(x) = \frac{1}{\sqrt{2\pi}} e^{-\frac{x^2}{2}}

さらに2変数の正規分布について考えます。まずは、これの積分を考えます。

f(x,y) = \frac{1}{2\pi} e^{-\frac{x^2}{2}} e^{-\frac{y^2}{2}}
\int_{-\infty}^{\infty} f(x,y) = \frac{1}{2\pi} \int_{x=-\infty}^{\infty} \int_{y=-\infty}^{\infty} e^{-\frac{x^2+y^2}{2}} dxdy

この式を、極座標系に書き換えて積分します。

\begin{aligned} x &= r \cos{\theta}\\ y &= r sin{\theta}\\ dxdy &= rdrd{\theta}\\ \end{aligned}

上記のように置換すると、以下のようになります。

\begin{aligned} \int_{x=-\infty}^{\infty} \int_{y=-\infty}^{\infty} f(x,y) &=\frac{1}{2\pi} \int_{r=0}^{\infty} \int_{\theta=0}^{2\pi} r e^{-\frac{r^2}{2}} drd\theta \\ &= \int_{r=0}^{\infty} r e^{-\frac{r^2}{2}} dr \end{aligned}

さらに以下のように変数変換を行います。

\begin{aligned} s&=\frac{r^2}{2} \\ ds&= rdr \end{aligned}

\begin{aligned} \int_{r=0}^{\infty} r e^{-\frac{r^2}{2}} dr = \int_{s=0}^{\infty} e^{-s} ds = 1 \end{aligned}

さて、積分が1だということが分かったので、先ほどの式をそのままPDFとして使用します。x,yが標準正規分布に従うならば、その確率密度は以下のようになるはずです。
p(x,y) = \frac{1}{2\pi} e^{-\frac{x^2+y^2}{2}}

MDFとconditional density functionを求めますが、その前に極座標系に変数変換します。

\begin{aligned} x &= \sqrt{R}\cos\theta \\ y &= \sqrt{R}\sin\theta \end{aligned}

この変数変換のヤコビアンは以下のようになります。確率密度関数p(x,y)は、x,yがある微小領域dx,dyをとる場合の確率ですので、変数変換する際は、ヤコビアンを掛ける必要があります。

\begin{aligned} det(J) &= \begin{array}{|cc|} \frac{\sigma x}{\sigma R} & \frac{\sigma x}{\sigma \theta} \\ \frac{\sigma y}{\sigma R} & \frac{\sigma y}{\sigma \theta} \end{array} = \begin{array}{|cc|} \frac{1}{2\sqrt{R}}\cos\theta & -\sqrt{R}\sin\theta \\ \frac{1}{2\sqrt{R}}\sin\theta & \sqrt{R}\cos\theta \end{array} = \frac{1}{2} \end{aligned}

したがって、先の確率密度関数は以下のようになります。
p(R, \theta) = \frac{1}{4\pi} e^{-\frac{R}{2}}

inversion methodを使うための積分は簡単なので割愛します。各CDFは以下のようになります。
\begin{aligned} P(\theta) &= \frac{1}{2\pi}\theta \\ P(R|\theta) &= 1 - e^{-\frac{R}{2}} \end{aligned}

P(θ)をξ1とし、P(R|θ)をξ2として、逆関数を求めます。
\begin{aligned} \theta &= 2\pi\xi_1 \\ R &= -2\ln(1-\xi_2) \end{aligned}

Rは極座標系における半径の平方根だったので、実際の半径:rは、Rの平方根となります。
\begin{aligned} \theta &= 2\pi\xi_1 \\ r &= \sqrt{-2\ln(1-\xi_2)} \end{aligned}

ξ2が0~1の一様乱数であるとするならば、1-ξ2はξ2と同義であると考えられるので、以下のように書き換えられます。
\begin{aligned} \theta &= 2\pi\xi_1 \\ r &= \sqrt{-2\ln(\xi_2)} \end{aligned}

上記の、2つの一様変数ξ1,ξ2から、2変数の標準正規分布のサンプリングを生成する変換式をBox Muller transformと呼びます。
2016-08-22修正
次回はこの導出を用いて、exp(x)を含むBRDFのImportance Samplingに取り組みたいと思います。
結局次の記事では使いませんでした。ただし、Box-Muller自体は有用で、レンダリングにおいて様々なところで利用されています。