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の反射分布が正しいかというところまでは確認できていません．ここからテスト項目を増やしていって，確認できることを広げていこうと思います．

ZinBox

2016年8月7日日曜日