2016年8月7日日曜日

BRDFのテスト

BRDFのテスト

このページは レイトレ合宿4!? アドベントカレンダー第7週目の記事です. レイトレ合宿まで1ヶ月をきりました. そろそろ提出用のシーン作成やレンダラーのデバッグに集中しているころではないでしょうか.

私はレンダラーの実装に不具合があるかどうかは, レンダリングした画像を見て違和感を持つかどうかで判断しています. この方法は主観的であり, 画像に違和感を持たなくても不具合が含まれている可能性は大いにあります. そこで,自動テストを実装して客観的に不具合の判定をするようにします.

このページでは物理ベースなBRDFの自動テストについて考えます.

物理ベースなBRDF

初めに,物理ベースなBRDFについて確認しておきます.

定義済み記号 説明
\(\omega_{n}\) 法線ベクトル
\(\omega_{i}\) 入射方向の単位ベクトル
\(\omega_{o}\) 反射方向の単位ベクトル
\(f_{r}\) BRDF

物理ベースなBRDFは,反射する放射エネルギーが負になることはないので

\[ \begin{align} \forall \boldsymbol{\omega_{i}} \ \forall \boldsymbol{\omega_{o}}, \ f_{r}(x, \boldsymbol{\omega_{o}} \leftarrow \boldsymbol{\omega_{i}}) \ cos \theta_{\boldsymbol{\omega_{o}}} \ge 0 \end{align} \]

\(cos \theta_{\boldsymbol{\omega_{o}}} \ge 0\) であるため,

\[ \begin{align} \forall \boldsymbol{\omega_{i}} \ \forall \boldsymbol{\omega_{o}}, \ f_{r}(x, \boldsymbol{\omega_{o}} \leftarrow \boldsymbol{\omega_{i}}) \ge 0 \tag{1} \end{align} \]

となります. また,物理ベースなBRDFは以下の2つの特性を持ちます.

ヘルムホルツの相反性 (Helmholtz Reciprocity)

BRDFの値は入射・出射方向が入れ替わっても変化しません.

\[ \begin{align} f_{r}(x, \boldsymbol{\omega_{o}} \leftarrow \boldsymbol{\omega_{i}}) = f_{r}(x, \boldsymbol{\omega_{i}} \leftarrow \boldsymbol{\omega_{o}}) \tag{2} \end{align} \]

エネルギー保存則 (Energy Conservation)

ある点から反射した総出射光束は総入射光束以下となります.

\[ \begin{align} \forall \boldsymbol{\omega_{i}}, \ \int_{\Omega} f_{r}(x, \boldsymbol{\omega_{o}} \leftarrow \boldsymbol{\omega_{i}}) \ cos \theta_{\boldsymbol{\omega_{o}}} d \sigma(\boldsymbol{\omega_{o}}) \le 1 \tag{3} \end{align} \]

モンテカルロ法によるテスト

\((1)\)\((2)\)\((3)\)の自動テストを実装します. 式\((3)\)は解析的に解くことが難しいため,モンテカルロ法を用いてテストを行います. モンテカルロ法を用いた場合,式\((3)\)は以下のように書きかえられます.

\[ \begin{align} \forall \boldsymbol{\omega_{i}}, \frac{1}{N} \sum_{k=1}^{N} \frac{f_{r}(x, \boldsymbol{\omega_{ok}} \leftarrow \boldsymbol{\omega_{i}}) cos \theta_{\boldsymbol{\omega_{ok}}}} {p_{u} \left( \boldsymbol{\omega_{ok}} \right)} \le 1 \tag{4} \end{align} \]

ここで,\({p_{u} \left( \boldsymbol{\omega_{o}} \right)}\)\(\boldsymbol{\omega_{o}}\)の確率密度関数 (PDF) です. このテストでは,半球上の方向を一様にサンプルするようにします \(\left({p_{u} \left( \boldsymbol{\omega_{o}} \right)} = \frac{1}{2 \pi} \right)\)

C++言語で GoogleTest を用いたテストの擬似コードを記述します.

// エネルギー保存則のテスト
TEST(BrdfTest, EnergyConservationTest)
{
  Brdf brdf; // テスト対象であるBRDFの宣言,初期化
  const int N = 1000000;
  // 様々な入射方向に対してテストを行う
  for (int i = 0; i < N; ++i) {
    const float energySum = 0.0;
    const Vector3 vin = sampleDirectionUniformly(); // 半球上の方向を一様にサンプル
    for (int o = 0; o < N; ++i) {
      const Vector3 vout = sampleDirectionUniformly();
      const float f = brdf.f(vin, vout); // BRDFの評価
      ASSERT_LE(0.0, f); // 式(1)のテスト
      const float cosTheta = dot(normal, vout);
      energySum += f * cosTheta * (2.0 * pi); // pi は円周率
    }
    const float energy = energySum / (float)N;
    ASSERT_GE(1.0, energy); // 式(4)のテスト
  }
}
// ヘルムホルツの相反性のテスト
TEST(BrdfTest, HelmHoltzReciprocityTest)
{
  Brdf brdf;
  const int N = 1000000;
  for (int i = 0; i < N; ++i) {
    const Vector3 vin = sampleDirectionUniformly();
    const Vector3 vout = sampleDirectionUniformly();
    const float f1 = brdf.f(vin, vout); // BRDFの評価
    const float f2 = brdf.f(-vout, -vin); // 方向を入れ替えたBRDFの評価
    const float error = 1.0e-5; // 許容できる誤差を定義
    ASSERT_FLOAT_EQ(f1, f2, error); // 式(2)のテスト
  }
}

実際にテストを実装する場合は, モンテカルロ法による\(N\)回試行の和を求める際に情報落ちの危険があるため, 補正加算も考慮したほうがいいと思います.

重点的サンプリング (Importance sampling) のテスト

BRDFによっては,反射方向を重点的サンプリングする式が定義してあることがあります. 重点的サンプリングの際の反射方向のPDFを \(p \left( \boldsymbol{\omega_{o}} \right)\) とします. 確率密度関数の性質から,

\[ \begin{align} \int_{\Omega} p \left( \boldsymbol{\omega_{o}} \right) d \sigma \left( \boldsymbol{\omega_{o}} \right) = 1 \tag{5} \end{align} \]

となります.式\((3)\)と同様に,解析的に解くことが難しいためモンテカルロ法を 用いて以下のように書きかえます.

\[ \begin{align} \frac{1}{N} \sum_{k=1}^{N} \frac{p \left( \boldsymbol{\omega_{ok}} \right)} {p_{u} \left( \boldsymbol{\omega_{ok}} \right)} \approx 1 \tag{6} \end{align} \]

擬似コードを以下に記述します.

// PDFのテスト
TEST(BrdfTest, PdfTest)
{
  Brdf brdf;
  const int N = 1000000;
  // 様々な入射方向に対してテストを行う
  for (int i = 0; i < N; ++i) {
    float pdfSum = 0.0;
    const Vector3 vin = sampleDirectionUniformly();
    for (int o = 0; o < N; ++o) {
      const Vector3 vout = sampleDirectionUniformly();
      const float pdf = brdf.pdf(vin, vout); // 反射方向のPDFを評価
      pdfSum += pdf * (2.0 * pi);
    }
    const float pdf = pdfSum / (float)N;
    const float error = 1.0e-5; // 許容できる誤差を定義
    ASSERT_FLOAT_NEAR(1.0, pdf, error); // 式(6)のテスト
  }
}

終わりに

BRDFのテストについて考えました. このテストで確認できたことは,物理に基づいたBRDFの定義を満たしているか ということであり,BRDFの反射分布が正しいかというところまでは確認できていません. ここからテスト項目を増やしていって,確認できることを広げていこうと思います.

2016年3月5日土曜日

Rough Specular Surface (GGX BRDF)

Rough Specular Surface (GGX BRDF)

本ページでは,モンテカルロレイトレーシング法において, 金属など絶縁体の粗い表面を表現する方法を説明する. 言い換えると,粗い表面モデルのBRDFの実装に必要なパラメータや数式, 重点的サンプリングの方法を説明する. 本ページでは,粗い表面モデルとして[1]のGGXモデルを用いる.

記号について

本ページでは以下の記号は定義済として扱う.

記号 説明
\(\omega_{n}\) Macrosurfaceの単位法線ベクトル
\(\omega_{m}\) Microsurfaceの単位法線ベクトル,本ページではマイクロファセット法線と呼ぶ
\(\omega_{i}\) 入射方向の単位ベクトル
\(\omega_{o}\) 反射方向の単位ベクトル
\(f_{r}\) BRDF
\(n_{i}\) 入射側の媒体の屈折率
\(n_{o}\) 透過側の媒体の屈折率
\(\eta_{o}\) 透過側の媒体の消衰係数
\(n\) \(\frac{n_{o}}{n_{i}}\)
\(\eta\) \(\frac{\eta_{o}}{n_{i}}\)
\(F\) フレネル項
\(\chi^{+} \left( x \right)\) \(\begin{cases} 1 & (0 < x) \\ 0 & (otherwise) \end{cases}\)

マイクロファセットモデルとGGX

マイクロファセットモデルは様々な 方向を向いた微小平面で構成される表面(Microsurface)を 単純な平面(Macrosurface)で近似したモデルである. このモデルでは,単純化のため波動光学や表面での2回以上の反射は無視する. マイクロファセットモデルは2つの統計量, 法線分布関数 \(D\) と Shadowing-Masking Function \(G\) を用いて表される.

表面の粗さについて

GGXは微小平面で見れば鏡面反射しかしていないが, 様々な法線を持った微小平面から構成されているため 表面が粗いほど拡散反射をしているように見える. 表面の粗さは \(roughness\) パラメータを用いて制御する. \(roughness\) がとる値の範囲は \(0 < roughness \le 1\) である. 値が \(0\) に近いほど理想的な鏡面反射に近くなり, \(1\) に近いほど理想的な拡散反射に近くなる.

また,計算時に用いる表面粗さの値は, \(roughness\) を直接用いるのではなく, [3]で解説してあるように

\[ \begin{align} \alpha = \left( roughness \right)^{2} \tag{1} \end{align} \]

を用いるほうがより視覚的に線形比例した変化となる. 本ページでも, \(\alpha\) を計算時の表面粗さの値として用いる.

法線分布関数 (Normal Distribution Function) \(D\)

法線分布関数 \(D\) はMicrosurface上にできる法線の統計分布を表したものである. GGXでは \(D\) は以下の式で表される.

\[ \begin{align} D \left( \omega_{m} \right) = \chi^{+} \left( \omega_{n} \cdot \omega_{m} \right) \frac{\alpha^{2}} {\pi \left( \left( \alpha^{2} - 1 \right) \left( \omega_{n} \cdot \omega_{m} \right)^{2} + 1 \right)^{2}} \tag{2} \end{align} \]

ここでは簡易表記のため, \(D \left( \omega_{n}, \omega_{m}, \alpha \right)\)\(D \left( \omega_{m} \right)\) と表記している.

Shadowing-Masking Function \(G\)

Shadowing-Masking Function は[2]の論文では V-cavityモデルとSmithモデルの2通りのモデルを載せているが, ここではSmithモデルを説明する(V-cavityモデルについては[2]を参照).

Shadowing-Masking Function \(G\) はある法線を持つMicrosurfaceが可視かどうかの比率 を表したものである. \(G_{2}\)は入射方向,反射方向の双方向から可視かどうかの比率を表しており, \(G_{1}\)は入射方向もしくは反射方向のどちらか単方向から可視かどうかの比率を表している.

\[ \begin{align} G_{2} \left( \omega_{i}, \omega_{o}, \omega_{m} \right) = & G_{1} \left( \omega_{i}, \omega_{m} \right) G_{1} \left( \omega_{o}, \omega_{m} \right) \tag{3} \\ G_{1} \left( \omega, \omega_{m} \right) = & \chi^{+} \left( \left( \omega \cdot \omega_{n} \right) \left( \omega \cdot \omega_{m} \right) \right) \frac{2 \left( \omega \cdot \omega_{n} \right)} { \left( \omega \cdot \omega_{n} \right) \sqrt{\alpha^{2} + \left( 1 - \alpha^{2} \right) \left( \omega \cdot \omega_{n} \right)^{2} }} \tag{4} \end{align} \]

ここでは簡易表記のため, \(G_{1} \left( \omega, \omega_{n}, \omega_{m}, \alpha \right)\)\(G_{1} \left( \omega, \omega_{m} \right)\) と表記している.

GGX BRDF

GGXのBRDFは以下の式になる.

\[ \begin{align} f_{r} \left( x, \omega_{i}, \omega_{o} \right) = \frac{1}{4 \left( \omega_{i} \cdot \omega_{n} \right) \left( \omega_{o} \cdot \omega_{n} \right)} F \left( \omega_{i}, \omega_{m} \right) G_{2} \left( \omega_{i}, \omega_{o}, \omega_{m} \right) D \left( \omega_{m} \right) \tag{5} \end{align} \]

ここでは簡易表現のため, \(F \left( \omega, \omega_{m}, n, \eta \right)\)\(F \left( \omega, \omega_{m} \right)\) と表記している.

ExplicitConnection のように, マイクロファセット法線 \(\omega_{m}\) がわからない状態で \(f_{r} \left( x, \omega_{i}, \omega_{o} \right)\) を評価する場合は, \(\omega_{m}\) は 反射ハーフベクトル \(\omega_{h_{r}}\) を用いて,

\[ \begin{align} \omega_{m} = & \frac{\omega_{h_{r}}} {\left| \left| \omega_{h_{r}} \right| \right|} \tag{6} \\ \omega_{h_{r}} = & \omega_{i} + \omega_{o} \tag{7} \end{align} \]

から求める.

重点的サンプリング (Importance Sampling)

反射方向 \(\omega_{o}\) の重点的サンプリングを行う. GGX BRDF では, マイクロファセット法線 \(\omega_{m}\) を重点的サンプリングし, 入射方向 \(\omega_{i}\)\(\omega_{m}\) から正反射方向 \(\omega_{o}\) を求める.

\(\omega_{m}\) は 投影法線分布関数(Projected Normal Distribution Function) \(D_{\omega_{i}} \left( \omega_{m} \right)\) に従ってサンプリングする. \(D_{\omega_{i}} \left( \omega_{m} \right)\) は, \(D \left( \omega_{m} \right)\) の中で 入射方向 \(\omega_{i}\) から見える法線のみの分布を示したもので,

\[ \begin{align} D_{\omega_{i}} \left( \omega_{m} \right) = \frac{ \left( \omega_{i} \cdot \omega_{m} \right) } { \left( \omega_{i} \cdot \omega_{n} \right) } G_{1} \left( \omega_{i}, \omega_{m} \right) D \left( \omega_{m} \right) \tag{8} \end{align} \]

で定義される (ここでは簡易表記のため, \(D_{\omega_{i}} \left( \omega_{n}, \omega_{m} \right)\)\(D_{\omega_{i}} \left( \omega_{m} \right)\) と表記している). マイクロファセット法線 \(\omega_{m}\) の pdf \(p \left( \omega_{m} \right) = D_{\omega_{i}} \left(\omega_{m} \right)\) から, 反射方向 \(\omega_{o}\) とそのpdf \(p \left( \omega_{o} \right)\)

\[ \begin{align} \omega_{o} = & 2 \left( \omega_{i} \cdot \omega_{m} \right) \omega_{m} - \omega_{i} \tag{9} \\ p \left( \omega_{o} \right) = & \left| \left| \frac{\delta \omega_{m}}{\delta \omega_{i}} \right| \right| p \left( \omega_{m} \right) \\ = & \frac{1}{4 \left( \omega_{i} \cdot \omega_{m} \right)} D_{\omega_{i}} \left( \omega_{m} \right) \tag{10} \end{align} \]

で求められる. ここで,\(\left| \left| \frac{\delta \omega_{m}}{\delta \omega_{i}} \right| \right| = \frac{1}{4 \left( \omega_{m} \cdot \omega_{i} \right)}\) はヤコビアンである.

重点的サンプリングの際の反射レイのWeightは

\[ \begin{align} \frac{f_{r} \left( x, \omega_{i}, \omega_{o} \right) \left( \omega_{o} \cdot \omega_{n} \right)} {p \left( \omega_{o} \right)} = F \left( \omega_{i}, \omega_{m} \right) G_{1} \left( \omega_{o}, \omega_{m} \right) \tag{11} \end{align} \]

となる.

\(\omega_{m}\) のサンプリング

\(D_{\omega_{i}}(\omega_{m})\) に従ってマイクロファセット法線 \(\omega_{m}\) をサンプリングする. サンプリング方法はV-cavity法とSmith法の2通りの方法があるが, ここではSmith法を載せる(V-cavity法については[2]を参照).

注意として,このサンプリング方法は macrosurfaceの法線が \(\omega_{n} = (0, 0, 1)\) の時の 入射角 \(\omega_{i} = (x_{i}, y_{i}, z_{i})\) を用いる. \(\omega_{n}\) がそれ以外の場合は \(\omega_{n} = (0, 0, 1)\) となるように \(\omega_{i}\) の基底変換を行う. \(\omega_{m}\) のサンプリング後は \(\omega_{i}\)\(\omega_{m}\) を 元の基底に変換する.

Smith法による\(\omega_{m}\)のサンプリング方法は 大きく分けて以下の5つのプロセスからなる.

  1. Stretch \(\omega_{i}\)
  2. Sample \(P_{\omega_{i}}^{22}(x_{\tilde{m}}, y_{\tilde{m}}, 1)\)
  3. Rotate \(\omega_{m}\)
  4. Unstretch \(\omega_{m}\)
  5. Normalize \(\omega_{m}\)

Stretch \(\omega_{i}\)

\[ \begin{align} \omega'_{i} \equiv \frac{(\alpha x_{i}, \alpha y_{i}, z_{i})} {\sqrt{(\alpha x_{i})^{2} + (\alpha y_{i})^{2} + z_{i}^2}} \end{align} \]

Sample \(P_{\omega_{i}}^{22}(x_{\tilde{m}}, y_{\tilde{m}}, 1)\)

\(x_{\tilde{m}}\) のサンプリング

\[ \begin{align} A \equiv & u_{1} \left( 1 + \sqrt{1 + tan^{2} \theta_{\omega'_{i}}} \right) - 1 \\ B \equiv & tan \theta_{\omega'_{i}} \\ x_{\tilde{m}1} \equiv & \frac{B}{A^{2} - 1} - \sqrt{\left(\frac{B}{A^{2} - 1}\right)^{2} - \frac{A^{2} - B^{2}}{A^{2} - 1}} \\ x_{\tilde{m}2} \equiv & \frac{B}{A^{2} - 1} + \sqrt{\left(\frac{B}{A^{2} - 1}\right)^{2} - \frac{A^{2} - B^{2}}{A^{2} - 1}} \\ x_{\tilde{m}} \equiv & \begin{cases} x_{\tilde{m}1} & (A < 0 \ \ or \ \ cot \theta_{\omega'_{i}} < x_{\tilde{m}2}) \\ x_{\tilde{m}2} & (otherwise) \end{cases} \end{align} \]

ここで,\(u_{1}\)\([0, 1)\) の乱数である.

\(y_{\tilde{m}}\) のサンプリング

\[ \begin{align} s \equiv & \begin{cases} 1 & (0.5 \lt u_{2}) \\ -1 & (otherwise) \end{cases} \\ t \equiv & 2s(u_{2} - 0.5) \\ z \equiv & s \frac{0.46341 t − 0.73369 t^{2} + 0.27385 t^{3}} {0.597999 − t + 0.309420 t^{2} + 0.093073 t^{3}} \\ y_{\tilde{m}} \equiv & z \sqrt{1 + x_{\tilde{m}}^{2}} \end{align} \]

ここで,\(u_{2}\)\([0, 1)\) の乱数である.

Rotate \(\omega_{m}\)

\[ \begin{align} \left( \begin{array}{c} x_{\tilde{m}} \\ y_{\tilde{m}} \end{array} \right) = \left( \begin{array}{c} cos \phi_{\omega_{i}'} & -sin \phi_{\omega_{i}'} \\ sin \phi_{\omega_{i}'} & cos \phi_{\omega_{i}'} \end{array} \right) \left( \begin{array}{c} x_{\tilde{m}} \\ y_{\tilde{m}} \end{array} \right) \end{align} \]

Unstretch \(\omega_{m}\)

\[ \begin{align} \left( \begin{array}{c} x_{\tilde{m}} \\ y_{\tilde{m}} \end{array} \right) = \left( \begin{array}{c} \alpha x_{\tilde{m}} \\ \alpha y_{\tilde{m}} \end{array} \right) \end{align} \]

Normalize \(\omega_{m}\)

\[ \begin{align} \omega_{m} \equiv \frac{(-x_{\tilde{m}}, -y_{\tilde{m}}, 1)} {\sqrt{x_{\tilde{m}}^{2} + y_{\tilde{m}}^{2} + 1}} \end{align} \]

まとめ

パラメータ

名称 記号
表面粗さ \(roughness\) \(0 < roughness \le 1\)

レンダリング

名称 記号
BRDF \(f_{r} \left( x, \omega_{i}, \omega_{o} \right)\) \(\frac{1}{4 \left( \omega_{i} \cdot \omega_{n} \right) \left( \omega_{o} \cdot \omega_{n} \right)} F \left( \omega_{i}, \omega_{m} \right) G_{2} \left( \omega_{i}, \omega_{o}, \omega_{m} \right) D \left( \omega_{m} \right)\)
pdf \(p \left( \omega_{o} \right)\) \(\frac{1}{4 \left( \omega_{i} \cdot \omega_{m} \right)} D_{\omega_{i}} \left( \omega_{m} \right)\)
反射レイのWeight \(\frac{f_{r} \left( x, \omega_{i}, \omega_{o} \right) \left( \omega_{o} \cdot \omega_{n} \right)} {p \left( \omega_{o} \right)}\) \(F \left( \omega_{i}, \omega_{m} \right) G_{1} \left( \omega_{o}, \omega_{m} \right)\)

以上,GGX BRDFについてまとめた.

参考文献

  1. B. Walter, S. R. Marschner, H. Li, K. E. Torrance: Microfacet models for refraction through rough surfaces (2007)
  2. E. Heitz, E. D'Eon: Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals (2014)
  3. B. Burley: Physically-Based Shading at Disney (2012)