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

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

変数変換の式からPDFを求める

とある空間:Sに属する2変数:s1,s2が既知の確率密度:ps(s1, s2)で分布するとして、これをSの写像空間:Tへ変換したいとします。写像空間:T内のある領域:Bにt1,t2が内包される確率:P(B)は以下のように表せます。pt()は、t1,t2のTにおける確率密度です。

P(B) = \int_{B} p_t(t_1, t_2) \textrm{d}t_1 \textrm{d}t_2

空間:Tの領域:Bの像に相当する、空間:Sの領域:ABにs1,s2が内包される確率は、上記P(B)と等しくなるはずです。従って下記の等式が成り立つはずです。

\int_{A_B} p_s(s_1, s_2) \textrm{d}s_1 \textrm{d}s_2 = \int_{B} p_t(t_1, t_2) \textrm{d}t_1 \textrm{d}t_2

これを踏まえ、s1,s2からt1,t2への変数変換を行うと以下の式が得られます。

\int_{A_B} p_s(s_1, s_2) \textrm{d}s_1 \textrm{d}s_2 = \int_{B} p_s(s_1(t_1, t_2), s_2(t_1, t_2)) |det(J)| \textrm{d}t_1 \textrm{d}t_2 \\ |det(J)| = | \frac{\partial s_1}{\partial t_1}\frac{\partial s_2}{\partial t_2} - \frac{\partial s_2}{\partial t_1}\frac{\partial s_1}{\partial t_2} |

このことから、以下の式が成り立つことがわかります。

p_t(t_1, t_2) = p_s(s_1(t_1, t_2), s_2(t_1, t_2)) |det(J)|

上記の式を簡単な例に当てはめてみます。まず、空間:Sに属するs1,s2が0~1の一様乱数:u,vで、空間:Tに属するt1,t2が、BRDFをサンプリングする際に使用する半球上の極座標の偏角:θ,φとして考えます。これを半球上の一様分布として変換するケースを考えます。変数は以下のように変換します。この変換は、最初の記事で求めた式にになります。

\begin{aligned} u &= \cos\theta\\ v &= \frac{1}{2\pi}\phi \end{aligned}

ヤコビアンは以下のようになります。

| \frac{\partial u}{\partial \theta}\frac{\partial v}{\partial \phi} - \frac{\partial v}{\partial \theta}\frac{\partial u}{\partial \phi} | = | -\sin\theta\frac{1}{2\pi} - 0 | = \frac{\sin\theta}{2\pi}

ps(u,v)は、u,vの分布は一様で空間の大きさは1なので、1です。 したがって、pt(θ,φ)は以下のようになります。

p_t(\theta,\phi) = \frac{\sin\theta}{2\pi}

変数変換の式からPDFを求めることができました。求められたpt()は、半球全域で積分するとちょうど1になります。値はdθ,dφが形成する半球上の微小面積の大きさに沿っています。また、pt()をベクトルの関数とするならば、微小立体角はdωで表現され、半球全域において一様なので以下のようになります。

p_t(\textbf{t}) = \frac{1}{2\pi}

したがって、pt()を極座標の偏角:θ,φの関数とした場合と、ベクトルの関数とした場合では、以下のような関係となります。

p_t(\theta,\phi)\frac{1}{\sin\theta} = p_t(\textbf{t})

この関係は、 dθ,dφの形成する微小面積と、dωの形成する微小立体角の比によるもので、PDFの関数が連続である限り成り立ちます。

PDFがハーフベクトルの関数の場合について

上記を踏まえ、もう少し実際のケースに当てはめていきたいと思います。多くのBRDFはハーフベクトル:hを関数に取ります。 したがって、 BRDFと比例関係であることが望ましいPDFも、ハーフベクトルの関数であれば、BRDFに近いプロファイルのPDFを作ることができるかもしれないと考えることは自然と言えると思います。ここでは、一様な0~1の変数u,vより、hを生成します。加えて、hのPDF:ph()を算出したいと思います。

p_h(\textbf{h}) = p_{uv}(u(\theta_h, \phi_h), v(\theta_h, \phi_h)) | \frac{\partial u}{\partial \theta_h}\frac{\partial v}{\partial \phi_h} - \frac{\partial v}{\partial \theta_h}\frac{\partial u}{\partial \phi_h} | \frac{1}{\sin\theta_h}

p_uv()は一様分布で、空間の大きさが１なので、以下のようになります。

p_h(\textbf{h}) = | \frac{\partial u}{\partial \theta_h}\frac{\partial v}{\partial \phi_h} - \frac{\partial v}{\partial \theta_h}\frac{\partial u}{\partial \phi_h} | \frac{1}{\sin\theta_h}

これでph()を算出することが可能になるのですが、実際のモンテカルロ積分でサンプリングしたいのは、ハーフベクトル:hの方向ではなく、反射ベクトル:oの方向です。したがって、oのPDF:po()が必要となります。
ハーフベクトル:hと反射ベクトル:oの関係は、入射ベクトル:iを全球の天頂（+z軸）にとる極座標系を用いることで、非常にシンプルに表すことができます。 この座標系の偏角は*をつけて区別します。

\begin{aligned} \theta^{\ast}_{o} &= 2\theta^{\ast}_{h} \\ \phi^{\ast}_{o} &= \phi^{\ast}_{h} \end{aligned}

これを変数変換と考えて、先の式に当てはめれば、po()は以下のように求めることができます。ヤコビアンの後ろの項は、h,oベクトルを偏角で表した場合の微小面積の比になります。

\begin{aligned} p_o(\textbf{o}) &= p_h(\textbf{h}) | \frac{\partial \theta^{\ast}_{h}}{\partial \theta^{\ast}_{o}}\frac{\partial \phi^{\ast}_{h}}{\partial \phi^{\ast}_{o}} - \frac{\partial \phi^{\ast}_{h}}{\partial \theta^{\ast}_{o}}\frac{\partial \theta^{\ast}_{h}}{\partial \phi^{\ast}_{o}} | \frac{\sin\theta^{\ast}_{h}}{\sin\theta^{\ast}_{o}} \\ &= p_h(\textbf{h}) \frac{1}{2} \frac{\sin\theta^{\ast}_{h}}{\sin 2\theta^{\ast}_{h}} \\ &= \frac{p_h(\textbf{h})}{4\cos\theta^{\ast}_{h}} \\ &= \frac{p_h(\textbf{h})}{4(\textbf{h}\cdot\textbf{i})} \end{aligned}

上記の通り、式を整理すると、ph(h)とhとiの内積のみが残り、*を用いた座標系の表現は無くなります。したがって、この変換を考える際に導入した座標系は、モンテカルロ積分する際に考慮する必要はありません。

isotropic Ward BRDFのspecular項を考える

次にWardのspecular項をisotropicなケースに限定したBRDFのImportanceSamplingを考えてみたいと思います。

f^{iso}_{r}(\textbf{i},\textbf{o}) = \frac{\rho_s}{4\pi\alpha^2\sqrt{\cos\theta_i\cos\theta_o}}e^{-\frac{\tan^2\theta_h}{\alpha^2}}

ρsとαは、それぞれspecular係数とroughnessで、定数として考えます。まず、先の通りu,vとθh,φhの関係式を考えます。下記の関係式は参照のNotes on the Ward BRDFに記述されています。 この式は、anisotropic版で使用されているものを簡単にしたものです。以前の記事で説明した通り、PDFとImportanceSamplingの対象になる式は、比例関係にあるのが望ましいので、uとθhの関係式はBRDFのexpの部分をそのまま使っています。 なぜならば、この関係式の偏微分がPDFの係数項として出てくるからです。偏微分可能であればこれが望ましいといえると思います。vとφhの関係式は、BRDF中にφhが無いので、全周2πに線形で変換します。

\begin{aligned} u &= e^{-\frac{\tan^2\theta_h}{\alpha^2}}\\ v &= \frac{\phi_h}{2\pi} \end{aligned}

これを前述の式に当てはめれば、以下のようになります。

\begin{aligned} p^{iso}_{h}(\textbf{h}) &= | (e^{-\frac{\tan^2\theta_h}{\alpha^2}} \frac{2\tan\theta_h}{\alpha^2\cos^2\theta_h}) (\frac{1}{2\pi}) -(0)(0) | \frac{1}{\sin\theta_h} \\ &= \frac{1}{\pi\alpha^2cos^3\theta_h}e^{-\frac{\tan^2\theta_h}{\alpha^2}} \end{aligned}

さらに、反射ベクトルのPDFとして変換すると以下のようになります。

p^{iso}_{o}(\textbf{o}) = \frac{1}{4(\textbf{h}\cdot\textbf{i})\pi\alpha^2cos^3\theta_h} e^{-\frac{\tan^2\theta_h}{\alpha^2}}

これで反射ベクトルをサンプリングする際のPDFを求めることができました。また、実際のサンプリングでは、u,vからθh,φhを求めるので、上記の変換式の逆関数が必要になります。

\begin{aligned} \theta_h &= \arctan(\alpha\sqrt{-ln(u)})\\ \phi_h &= 2\pi v \end{aligned}

したがって、サンプリングの関数は以下の様になります。

vec3 SampleIsoWard(float u1, float u2, float alpha, vec3 V)
{
  float theta = atan(alpha * sqrt(-log(u1)));
  float phi = 2.0f * PI * u2;

  vec3 H = vec3(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta))

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

モンテカルロ積分のWeight関数は、以下のようになります。

W^{iso}(\textbf{i},\textbf{o}) = \frac{f^{iso}(\textbf{i},\textbf{o})}{p^{iso}_o(\textbf{o})} = \frac{\rho_s(\textbf{h}\cdot\textbf{i})\cos^3\theta_h}{\sqrt{\cos\theta_i\cos\theta_o}}