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

前回に引き続きImportanceSamplingの基礎について勉強したいと思います。
今回も、Bruce WalterさんのNotes on the Ward BRDFを読みたいと思います。

Ward BRDFのspecular項

前回の記事では、Ward BRDFをisotropicな場合に限定したBRDFのImportanceSamplingを考えてみましたが、今回はanisotropicなケースを考えてみたいと思います。

f_{r}(i,o) = \frac{\rho_s}{4\pi\alpha_{x}\alpha_{y}\sqrt{\cos\theta_i\cos\theta_o}} e^{-\tan^2\theta_h(\frac{\cos^2\phi_h}{\alpha_{x}^2}+\frac{\sin^2\phi_h}{\alpha_{y}^2})}

ρｓはスペキュラ係数、αx,αyは、接空間におけるx,y方向のroughnessに相当します。αxとαyが同じ値のときはisotropicになり、前回の記事で取り扱った式と同一になります。前回の記事と同様に、一様乱数からの変数変換の式を先に決め、それを偏微分してPDFを導出し、ImportanceSamplingできるようにしたいと思います。

変数変換の式を考える

0~1の一様乱数u,vから、ハーフベクトルの偏角θh,φhへの変換を考えます。変換式は、参照のNotes on the Ward BRDFに記述されているので、以下の立式の理由は後付けです。
まずθhについてですが、前回同様にBRDFのexp()の部分をそのまま使います。理由も前回と同様で、exp()関数を偏微分すれば、そのまま係数項として出てくるはずで、BRDFとPDFを比例関係に近づけることができると思います。

u=e^{-\tan^2\theta_h(\frac{\cos^2\phi_h}{\alpha_{x}^2}+\frac{\sin^2\phi_h}{\alpha_{y}^2})}

次にφhについてですが、BRDFの式を見ると、αx,αyの値に従って、φ周りで値の大きな領域（specular）が楕円状になることがわかります。例えば、αx>αyならば接空間のｘ軸方向に広いspecularが現れるはずです。この場合ならｘ軸に近い角度のφのサンプリングを重点的に行うほうが好ましいと言えます。つまり、φのサンプリングはαx,αyの値にそった楕円状なるのが好ましいといえると思います。以下のように、tan()を使ってx,y軸でαx,αyの比になる楕円状にサンプリングを行います。

\phi_h = \arctan(\frac{\alpha_y}{\alpha_x}\tan(2\pi v))

上記の逆関数は以下のようになります。

v = \frac{1}{2\pi} \arctan(\frac{\alpha_x}{\alpha_y}\tan\phi_h)

偏微分を先にやる

この後の計算の流れは、前回の記事と同様なので、先にu,vの変数変換式の偏微分を計算します。

u = e^{-\tan^2\theta_h(\frac{\cos^2\phi_h}{\alpha_{x}^2}+\frac{\sin^2\phi_h}{\alpha_{y}^2})} \\ \frac{\partial u}{\partial \theta_h} = e^{-\tan^2\theta_h(\frac{\cos^2\phi_h}{\alpha_{x}^2}+\frac{\sin^2\phi_h}{\alpha_{y}^2})} \frac{-2\sin\theta_h}{\cos^3\theta_h} (\frac{\cos^2\phi_h}{\alpha_{x}^2}+\frac{\sin^2\phi_h}{\alpha_{y}^2})

v = \frac{1}{2\pi} \arctan(\frac{\alpha_x}{\alpha_y}\tan\phi_h) \\ \begin{aligned} \frac{\partial v}{\partial \phi_h} &= \frac{\alpha_x}{\alpha_y}\frac{1}{\cos^2\phi_h} \frac{1}{2\pi}\frac{1}{1+(\frac{\alpha_x}{\alpha_y}\tan\phi_h)^2} \\ &= \frac{1}{2\pi} \frac{\alpha_x} {\alpha_y\cos^2\phi_h(1+(\frac{\alpha_x}{\alpha_y}\tan\phi_h)^2)} \\ &= \frac{1}{2\pi} \frac{\alpha_x\alpha_y} {\alpha_y^2\cos^2\phi_h+\alpha_x^2\sin^2\phi_h} \end{aligned}

PDFの計算

次に、上記の計算結果を用いてPDFの計算をします。exp()の項は最後まで保存されるので省略して書きます。

\begin{aligned} p_h(h) &= | (e^{-\tan^2\theta_h(\frac{\cos^2\phi_h}{\alpha_{x}^2}+\frac{\sin^2\phi_h}{\alpha_{y}^2})} \frac{-2\sin\theta_h}{\cos^3\theta_h} (\frac{\cos^2\phi_h}{\alpha_{x}^2}+\frac{\sin^2\phi_h}{\alpha_{y}^2})) (\frac{1}{2\pi} \frac{\alpha_x\alpha_y} {\alpha_y^2\cos^2\phi_h+\alpha_x^2\sin^2\phi_h}) -(0)(0) | \frac{1}{\sin\theta_h} \\ &= | exp(...) \frac{-1}{\pi} \frac{1}{\cos^3\theta_h} (\frac{\cos^2\phi_h}{\alpha_{x}^2}+\frac{\sin^2\phi_h}{\alpha_{y}^2}) \frac{\alpha_x\alpha_y}{\alpha_y^2\cos^2\phi_h+\alpha_x^2\sin^2\phi_h} |\\ &= | exp(...) \frac{-1}{\pi} \frac{1}{\cos^3\theta_h} \frac{\alpha_{y}^2\cos^2\phi_h+\alpha_{x}^2\sin^2\phi_h}{\alpha_x\alpha_y} \frac{1}{\alpha_y^2\cos^2\phi_h+\alpha_x^2\sin^2\phi_h} |\\ &= \frac{1}{\pi\cos^3\theta_h\alpha_x\alpha_y} e^{-\tan^2\theta_h(\frac{\cos^2\phi_h}{\alpha_{x}^2})+\frac{\sin^2\phi_h}{\alpha_{y}^2})} \end{aligned}

上記の通り、ハーフベクトルのPDFが求まりました。前回同様に反射ベクトルのPDFとして変換すると以下のようになります。

\begin{aligned} p_o(o) &= \frac{1}{4(h \cdot i)} \frac{1}{\pi\cos^3\theta_h\alpha_x\alpha_y} e^{-\tan^2\theta_h(\frac{\cos^2\phi_h}{\alpha_{x}^2})+\frac{\sin^2\phi_h}{\alpha_{y}^2})} \end{aligned}

サンプリングの計算

実際のサンプリングでは、u,vからθh,φhを求めるので、先に記した変換式の逆関数が必要になります。今回は、θhはφhの関数となっているのでφhを先に求める必要があります。
\begin{aligned} \phi_h &= \arctan(\frac{\alpha_y}{\alpha_x}\tan(2\pi v)) \\ \theta_h &= \arctan( \sqrt{ \frac {-\ln(u)} {\frac{\cos^2\phi_h}{\alpha_{x}^2}+\frac{\sin^2\phi_h}{\alpha_{y}^2}} } ) \end{aligned}

φhを求める際に、arctan()を使っています。この関数は非連続な関数なので、実際にプログラム内で使用する場合は注意が必要です。直接φhを求めずに、sin(φh),cos(φh)を求めて、そのまま使う方法も考えられます。その場合は、以下の様に連続的に求めることができます。

\begin{aligned} \tan\phi_h &= \frac{\alpha_y}{\alpha_x}\tan(2\pi v) \\ \frac{\sin\phi_h}{\cos\phi_h} &= \frac{\alpha_y\sin(2\pi v)}{\alpha_x\cos(2\pi v)}\\ \end{aligned} \\ \begin{aligned} x &= \alpha_x\cos(2\pi v)\\ y &= \alpha_y\sin(2\pi v)\\ \end{aligned} \\ \begin{aligned} \sin\phi_h &= \frac{y}{\sqrt{x^2+y^2}}\\ \cos\phi_h &= \frac{x}{\sqrt{x^2+y^2}}\\ \end{aligned}

同様の計算が、θhでも可能です。したがって、サンプリングの関数は以下の様になります。

vec3 SampleWard(float u1, float u2, vec3 V)
{
 float sinPhi, cosPhi;
 {
  float x = ax * cos(2.0f * PI *u2);
  float y = ay * sin(2.0f * PI *u2);
  float len = sqrt(x*x+y*y);
  sinPhi = y / len;
  cosPhi = x / len;
 }

 float sinTheta, cosTheta;
 {
  flaot x = sqrt((cosPhi*cosPhi)/(ax*ax)+(sinPhi*sinPhi)/(ay*ay));
  float y= sqrt(-log(u1));
  float len = sqrt(x*x+y*y);
  sinTheta = y / len;
  cosTheta = x / len;
 }

  vec3 H = vec3(sinTheta * cosPhi, sinTheta * sinPhi, cosTheta)

  return -V + 2.0f * H * dot(H, V)
}

一方でモンテカルロ積分のweight関数は以下のようになります。
\begin{aligned} W(i,o) &= \frac{f_r(i,o)}{p_o(o)} \\ &= \frac {\rho_s(h\cdot i)\cos^3\theta_h} {\sqrt{\cos\theta_i\cos\theta_o}} &= \frac {\rho_s(h\cdot i)(n \cdot h)^3} {\sqrt{(n \cdot i)(n \cdot o)}} \end{aligned}