Final gatheringについて

この記事はフォトンマッピングの実装の続きです。今回はフォトンマッピングの誤差を目立たなくするFinal gatheringという手法について解説します。

最初に書くのもあれですが、綺麗な画像を出すという目的では現在はこの手法を実装するメリットはあまりないです。前回のようなナイーブなフォトンマッピングを実装できた後はFinal gatheringをやるよりも、(Stochastic) progressive photon mappingの実装に移るのがおすすめです。

TL;DR

ある点に入射する光を直接照明 + 集光模様 + 間接照明に分解する
直接照明はNext event estimationで計算
集光模様は集光模様フォトンマップを用いて計算
間接照明はレイを一つ奥に飛ばして、その先で大域フォトンマップを用いて計算

Final gatheringとは

前回の内容を実装してコーネルボックスをレンダリングしてみると、以下のような低周波のノイズが出てしまうことを見ました。

このような低周波ノイズが出るのは、カーネルの半径が有限ゆえに隣接ピクセル間で同じフォトンを放射輝度の推定に用いてしまっていることが原因です。

理論的にはフォトンマップ構築に使うフォトン数を無限大にすれば、カーネルの半径が無限小になるのでこのノイズは消えてくれるはずです。

しかし実際にはメモリ的な制約により1000000フォトン程度が限界でしょう。この画像は1000000フォトンを使ってレンダリングしたものですが、それでもまだノイズが残っています。ナイーブなフォトンマッピングの実装ではこのノイズを消すことはほぼ不可能と言えます。

そこで、Jensenの青い本では Final gathering と呼ばれる手法を用いてこのノイズを目立たなくすることを行っています。

Final gatheringでは入射光を 直接照明(Direct Illumination), 集光模様(Caustics), 間接照明(Indirect Illumination) の3つに分解し、それぞれを別の方法で計算するようにします。以降ではこれについて解説していきます。

入射光の分割

ある点 $x$ に入射した光が視点 $z$ の方向にどれだけ反射するかを考えます。3点形式のレンダリング方程式を用いることで、これは以下のように計算できます。

$x$ : 注目点
$z$ : 視点
$L_o(x \to z)$ : $x$ から $z$ に反射される放射輝度
$f(y \to x \to z)$ : BSDF
$G(y, x)$ : 幾何項
$L_e(x \to z)$ : 点 $x$ での発光による放射輝度
$L_i(y \to x)$ : 点 $y$ から $x$ に入射する放射輝度
$\mathcal{M}$ : シーン上の点全体

L_o(x \to z) = L_e(x \to z) + \int_{\mathcal{M}} f(y \to x \to z)G(y, x)L_i(y \to x)dA(y) \tag{1}

ここで点 $x$ での入射光 $L_i(y \to x)$ を次のような3成分に分解します。

光源から $x$ に直接到達する光(直接照明)
光源からSpecular面を経て $x$ に到達する光(集光模様)
光源からDiffuse面を経て $x$ に到達する光(間接照明)

厳密には2, 3番目は両方間接照明と言えますが、ここでは2を集光模様, 3を間接照明と呼ぶことにします。

光源上の点集合を $\mathcal{M_l}$ , Specular面上の点集合を $\mathcal{M_s}$ , Diffuse面上の点集合を $\mathcal{M_d}$ とします。今 $\mathcal{M}$ が $\mathcal{M} = \mathcal{M_l} \cup \mathcal{M_s} \cup \mathcal{M_d}$ と分解できたとすると、式(1)は以下のように変形できます。

L_o(x \to z) = L_e(x \to z) + L_l(x \to z) + L_s(x \to z) + L_d(x \to z) \tag{2}

\begin{aligned} L_l(x \to z) &= \int_{\mathcal{M_d}} f(y \to x \to z)G(y, x)L_i(y \to x)dA(y) \\ L_s(x \to z) &= \int_{\mathcal{M_s}} f(y \to x \to z)G(y, x)L_i(y \to x)dA(y) \\ L_d(x \to z) &= \int_{\mathcal{M_d}} f(y \to x \to z)G(y, x)L_i(y \to x)dA(y) \end{aligned}

これをさらに点 $x$ におけるBSDFの種類で分割していきます。点 $x$ におけるBSDFが以下のようにDiffuse成分とSpecular成分の和で表されているとします。

$f_d$ : Diffuse部分に対応するBSDF
$f_s$ : Specular成分に対応するBSDF

f(y \to x \to z) = f_d(y \to x \to z) + f_s(y \to x \to z)

これを式(2)に代入してあげると、積分がさらに分割されて次のようになります。

L_o(x \to z) = L_e(x \to z) + L_{l,d}(x \to z) + L_{l,s}(x \to z) + L_{s,d}(x \to z) + L_{s,s}(x \to z) + L_{d,d}(x \to z) + L_{d,s}(x \to z) \tag{3}

\begin{aligned} L_{l,d}(x \to z) &= \int_{\mathcal{M_d}} f_d(y \to x \to z)G(y, x)L_i(y \to x)dA(y) \\ L_{l,s}(x \to z) &= \int_{\mathcal{M_d}} f_s(y \to x \to z)G(y, x)L_i(y \to x)dA(y) \\ L_{s,d}(x \to z) &= \int_{\mathcal{M_s}} f_d(y \to x \to z)G(y, x)L_i(y \to x)dA(y) \\ L_{s,s}(x \to z) &= \int_{\mathcal{M_s}} f_s(y \to x \to z)G(y, x)L_i(y \to x)dA(y) \\ L_{d,d}(x \to z) &= \int_{\mathcal{M_d}} f_d(y \to x \to z)G(y, x)L_i(y \to x)dA(y) \\ L_{d,s}(x \to z) &= \int_{\mathcal{M_d}} f_s(y \to x \to z)G(y, x)L_i(y \to x)dA(y) \end{aligned}

Final gatheringではこの6種類の項を別々の方法で計算していきます。以降では各項の計算方法について見ていきます。

各項の計算方法

$L_{l, d}, L_{l, s}$ の計算(直接照明)

$L_{l, d}, L_{l, s}$ の項(直接照明)はフォトンマップを使わずに、パストレーシングでよく使われるNext event estimationによって計算を行います。

手順は以下のような流れです。

シーン中から光源を1つ選ぶ
選んだ光源上で点 $y$ をサンプリング
現在の点 $x$ から $y$ に向けてレイを飛ばし、光源が見えるか判定
光源が見える場合は寄与を加算

今シーンに $N_l$ 個の光源があるとします。光源を一様にランダムに選ぶことにすると、 $i$ 番目の光源を選ぶ確率は以下で与えられます。

P_L(i) = \frac{1}{N_l}

続いて選んだ光源上で点をサンプリングします。光源上で一様に点サンプリングすることにすると、点サンプリングの確率密度は以下で与えられます。

$A$ : 光源の表面積

p(y) = \frac{1}{A}

あとは $L_{l, d}, L_{l, s}$ をモンテカルロ積分するだけです。以下のような形になります。

$N_{direct}$ : サンプル数
$A_i$ : $i$ 番目のサンプルで選択された光源の表面積
$y_i$ : $i$ 番目のサンプルでサンプリングされた光源上の点
$\theta_i$ : 点 $x$ における法線と $x \to y_i$ がなす角度
$\theta_{l_i}$ : 点 $y_i$ における法線と $y_i \to x$ がなす角度
$r_i$ : 点 $x$ と点 $y_i$ の距離

\begin{aligned} L_d(x \to z) &\approx \frac{1}{N_{direct}}\sum_{i=1}^{N_{direct}} \frac{f(y_i \to x \to z)G(y_i, x)L_i(y_i \to x)}{P_L p(y_i)} \\ &= \frac{N_l}{N_{direct}}\sum_{i=1}^{N_{direct}} \frac{A_i f(y_i \to x \to z)L_i(y_i \to x) \cos \theta_i \cos \theta_{l_i}}{r_i^2} \tag{4} \end{aligned}

$L_{s,d}$ の計算(集光模様)

$L_{s,d}$ の項(集光模様)は点 $x$ の近傍にあるフォトンを 集光模様フォトンマップ(Caustics photon map) から持ってきて計算します。

集光模様フォトンマップはSpecular面が続いた後にDiffuse面に当たったフォトンだけを格納したフォトンマップです。正規表現を用いた経路表現を使うと、集光模様フォトンマップはL(S+)Dという経路を持つフォトンを格納しています。

何故このようなフォトンマップを用意する必要があるかについては後ほど説明します。

集光模様フォトンマップを用いた放射輝度推定の方法は前回説明した通常のフォトンマップを用いた放射輝度推定の方法と同じなので、ここでは説明を省略します。

$L_{d,d}$ の計算(間接照明)

計算する対象の式は以下の形でした。

L_{d,d}(x \to z) = \int_{\mathcal{M_d}} f(y \to x \to z)G(y, x)L_i(y \to x)dA(y)

間接照明の計算では $L_i$ を 大域フォトンマップ(Global photon map) を用いて計算します。大域フォトンマップとは、集光模様フォトンマップのような制約を何も入れていない通常のフォトンマップのことです。

点 $y$ は点 $x$ からBSDFを用いた方向サンプリングによって方向 $\vec{\omega}$ を生成し、その方向にレイを飛ばすことで生成します。

この時 $L_{d, d}$ のモンテカルロ積分は以下のように書けます。

$N_{indirect}$ : 方向サンプリングのサンプル数
$\vec{\omega_i}$ : $i$ 番目の方向サンプル
$y_i$ : 点 $x$ から方向 $\vec{\omega_i}$ にレイを飛ばして得られた点
$\theta_i$ : 点 $x$ における法線と方向 $\vec{\omega_i}$ のなす角度
$p(\vec{\omega_i})$ : 方向サンプリングの確率密度

\begin{aligned} L_{d, d}(x \to z) &\approx \frac{1}{N_{indirect}}\sum_{i=1}^{N_{indirect}} \frac{f(y_i \to x \to z)G(y_i, x)L_i(y_i \to x)}{p(y_i)} \\ &= \frac{1}{N_{indirect}}\sum_{i=1}^{N_{indirect}} \frac{f(y_i \to x \to z)L_i(y_i \to x)\cos\theta_i}{p(\vec{\omega_i})} \tag{5} \end{aligned}

$y$ がSpecular面の場合は $y$ の近傍にはフォトンがないので、大域フォトンマップを用いて $L_i$ を推定することができません。この場合は $y$ からさらにレイを再帰的に追跡するようにします。するとどこかでDiffuse面に当たるので、その点で大域フォトンマップを用いた放射輝度推定を行って追跡が終了します。

$L_{s,s}, L_{d,s}$ の計算

最後に残りの $L_{s,s}, L_{s,d}$ の計算方法について見ていきます。この2つの項は点 $x$ がSpecular面の場合に対応します。

まず $L_{s,s}$ の計算について考えます。点 $x$ はSpecular面なので集光模様を計算しようと思っても点 $x$ の近傍にはフォトンがありません。なのでこの場合は $L_{d,d}$ の計算と同様に方向サンプリングを行ってレイを再帰的に追跡するようにします。

$L_{d,s}$ の計算も同様で、こちらも $L_{s,s}$ の計算と同様にレイを再帰的に追跡します。

まとめると点 $x$ がSpecular面の場合はパストレーシング同様にレイを再帰的に追跡します。再帰的に追跡した結果、どこかでDiffuse面と当たります。するとそこで $L_{l,d}, L_{s,d}, L_{d,d}$ が計算されて追跡が終了します。

アルゴリズムの全体像

Final gatheringのアルゴリズムの全体をまとめると以下のようになります。

$x$ がDiffuse面の場合
- $L_{l,d}$ (直接照明)を式(4)を用いて計算して寄与を加算
- $L_{s,d}$ (集光模様)を集光模様フォトンマップを用いて計算して寄与を加算
- $L_{d,d}$ (間接照明)を式(5)を用いて計算して寄与を加算
$x$ がSpecular面の場合
- Specular面の反射方向にレイを再帰的に追跡

実装

以下ではC++によるFinal gatheringの実装について説明していきます。

コード全体は以下で見ることができるので、あわせて参考にして下さい。

https://github.com/yumcyaWiz/photon_mapping

大域フォトンマップの構築

ここは前回と同じなので省略します。

集光模様フォトンマップの構築

集光模様フォトンマップはL(S+)Dのような経路を持つフォトンを格納しているのでした。Specular面が続いた後にDiffuse面に当たったらそこでフォトンを格納し、追跡を打ち切るようにします。

実装はほとんど大域フォトンマップの構築と同じです。前回いた面がSpecular面かどうかを表すフラグを追加して、そのフラグを用いてフォトンの追加を行うか判定します。

// causticsPhotonMap: 集光模様フォトンマップ

// OpenMPを用いてフォトンを並列に追跡する
#pragma omp parallel for
for (int i = 0; i < nPhotonsCaustics; ++i) {
  auto& sampler_per_thread = *samplers[omp_get_thread_num()];

  // 光源から出るレイをサンプリングする
  Vec3f throughput;
  Ray ray = sampleRayFromLight(scene, sampler_per_thread, throughput);

  // 前回いた面がSpecularかを表すフラグ
  bool prev_specular = false;

  // フォトンの追跡
  for (int k = 0; k < maxDepth; ++k) {
    IntersectInfo info;
    if (scene.intersect(ray, info)) {
      // 交差点のBSDFの種類を取得
      const BxDFType bxdf_type = info.hitPrimitive->getBxDFType();

      // Diffuse面が連続した場合は終了
      if (!prev_specular && bxdf_type == BxDFType::DIFFUSE) {
        break;
      }

      // Specular面の後にDiffuse面に当たったらフォトンを追加する
      if (prev_specular && bxdf_type == BxDFType::DIFFUSE) {
#pragma omp critical
        {
          photons.emplace_back(throughput, info.surfaceInfo.position,
                                -ray.direction);
        }
        break;
      }

      prev_specular = (bxdf_type == BxDFType::SPECULAR);

      // ロシアンルーレット
      if (k > 0) {
        const float russian_roulette_prob =
            std::min(std::max(throughput[0],
                              std::max(throughput[1], throughput[2])),
                      1.0f);
        if (sampler_per_thread.getNext1D() >= russian_roulette_prob) {
          break;
        }
        throughput /= russian_roulette_prob;
      }

      // BSDFを用いた方向サンプリング
      Vec3f dir;
      float pdf_dir;
      const Vec3f f =
          info.hitPrimitive->sampleBxDF(-ray.direction, info.surfaceInfo,
                                        TransportDirection::FROM_LIGHT,
                                        sampler_per_thread, dir, pdf_dir);

      // throughputとレイを更新
      throughput *= f *
                    cosTerm(-ray.direction, dir, info.surfaceInfo,
                            TransportDirection::FROM_LIGHT) /
                    pdf_dir;
      ray = Ray(info.surfaceInfo.position, dir);
    } else {
      // フォトンがシーン外に飛んでいった場合
      break;
    }
  }
}

// 集光模様フォトンマップのビルド(kd-treeの構築)
causticsPhotonMap.setPhotons(photons);
causticsPhotonMap.build();

$L_{l,d}, L_{l,s}$ (直接照明)の計算

式(4)の計算を実装します。ここでは $N_{direct} = 1$ としています。

// Next event estimationで直接照明を計算する
// wo: 出射方向
Vec3f computeDirectIllumination(const Scene& scene, const Vec3f& wo,
                                const IntersectInfo& info,
                                Sampler& sampler) const {
  Vec3f Ld;

  // 光源のサンプリング
  float pdf_choose_light;
  const std::shared_ptr<Light> light =
      scene.sampleLight(sampler, pdf_choose_light);

  // 光源上の点サンプリング
  float pdf_pos_light;
  const SurfaceInfo light_surf = light->samplePoint(sampler, pdf_pos_light);

  // 点サンプリングの確率密度を方向サンプリングの確率密度に変換
  const Vec3f wi = normalize(light_surf.position - info.surfaceInfo.position);
  const float r = length(light_surf.position - info.surfaceInfo.position);
  const float pdf_dir =
      pdf_pos_light * r * r / std::abs(dot(-wi, light_surf.shadingNormal));

  // 光源に向かうレイ(影光線)の生成
  Ray ray_shadow(info.surfaceInfo.position, wi);
  ray_shadow.tmax = r - RAY_EPS;

  // 光源の可視性判定
  // 光源が見える場合は寄与を加算
  IntersectInfo info_shadow;
  if (!scene.intersect(ray_shadow, info_shadow)) {
    const Vec3f Le = light->Le(light_surf, -wi);
    const Vec3f f = info.hitPrimitive->evaluateBxDF(
        wo, wi, info.surfaceInfo, TransportDirection::FROM_CAMERA);
    const float cos = std::abs(dot(wi, info.surfaceInfo.shadingNormal));
    Ld = f * cos * Le / (pdf_choose_light * pdf_dir);
  }

  return Ld;
}

$L_{s,d}$ (集光模様)の計算

集光模様フォトンマップを用いた放射輝度推定を行います。やっていることは前回説明したフォトンマップを用いた放射輝度推定と全く同じです。

// 集光模様フォトンマップを用いて反射放射輝度の推定を行う
// wo: 出射方向
Vec3f computeCausticsWithPhotonMap(const Vec3f& wo,
                                    const IntersectInfo& info) const {
  // フォトンマップから近傍フォトンのインデックスを取得
  float max_dist2;
  const std::vector<int> photon_indices =
      causticsPhotonMap.queryKNearestPhotons(info.surfaceInfo.position,
                                              nEstimationCaustics, max_dist2);

  // 放射輝度の計算
  Vec3f Lo;
  for (const int photon_idx : photon_indices) {
    const Photon& photon = causticsPhotonMap.getIthPhoton(photon_idx);
    const Vec3f f = info.hitPrimitive->evaluateBxDF(
        wo, photon.wi, info.surfaceInfo, TransportDirection::FROM_CAMERA);
    Lo += f * photon.throughput;
  }
  if (photon_indices.size() > 0) {
    Lo /= (nPhotonsCaustics * PI * max_dist2);
  }

  return Lo;
}

$L_{d,d}$ (間接照明)の計算

式(5)の計算を実装します。ここでは $N_{indirect}=1$ としています。

computeRadianceWithPhotonMapはcomputeCausticsWithPhotonMapの中身を大域フォトンマップを使うように置き換えただけの関数です。

レイを飛ばした先がSpecular面の場合は再帰的に関数を呼び出す必要があるので、ここでは再帰用の関数を別に分けています。

// Final gatheringによって間接照明を計算する
// 再帰を表現するために別関数に分けている
// wo: 出射方向
// depth: 再帰の深さ
Vec3f computeIndirectIlluminationRecursive(const Scene& scene,
                                            const Vec3f& wo,
                                            const IntersectInfo& info,
                                            Sampler& sampler,
                                            int depth) const {
  if (depth >= maxDepth) return Vec3f(0);

  Vec3f Li;

  // BSDFを用いた方向サンプリング
  Vec3f dir;
  float pdf_dir;
  const Vec3f f = info.hitPrimitive->sampleBxDF(
      wo, info.surfaceInfo, TransportDirection::FROM_CAMERA, sampler, dir,
      pdf_dir);
  const float cos = std::abs(dot(info.surfaceInfo.shadingNormal, dir));

  // サンプリングした方向にレイを飛ばす
  Ray ray_fg(info.surfaceInfo.position, dir);
  IntersectInfo info_fg;
  if (scene.intersect(ray_fg, info_fg)) {
    // 交差点のBSDFの種類を取得
    const BxDFType bxdf_type = info_fg.hitPrimitive->getBxDFType();

    // Diffuse面の場合は大域フォトンマップを用いて放射輝度推定を行う
    if (bxdf_type == BxDFType::DIFFUSE) {
      Li += f * cos *
            computeRadianceWithPhotonMap(-ray_fg.direction, info_fg) /
            pdf_dir;
    }
    // Specular面の場合はレイを再帰的に追跡する
    else if (bxdf_type == BxDFType::SPECULAR) {
      Li += f * cos *
            computeIndirectIlluminationRecursive(
                scene, -ray_fg.direction, info_fg, sampler, depth + 1) /
            pdf_dir;
    }
  }

  return Li;
}

// Final gatheringによって間接照明を計算する
// wo: 出射方向
Vec3f computeIndirectIllumination(const Scene& scene, const Vec3f& wo,
                                  const IntersectInfo& info,
                                  Sampler& sampler) const {
  // 深さ0から再帰を始める
  return computeIndirectIlluminationRecursive(scene, wo, info, sampler, 0);
}

放射輝度の計算

最後にこれまで定義してきた関数を用いてカメラに到達する放射輝度の計算を実装します。前回の式(10)を計算します。

Diffuse面に当たった場合は直接照明, 集光模様, 間接照明の計算を上で定義してきた関数を用いて行うようにします。

Specular面に当たった場合は再帰的にレイを追跡します。ここでは再帰用に関数を分けています。

// rayの方向から来る放射輝度を計算する
// 再帰用に関数を分けている
// depth: 再帰の深さ
Vec3f integrateRecursive(const Ray& ray, const Scene& scene, Sampler& sampler,
                          int depth) const {
  // 再帰の深さが最大値を超えたら終了
  if (depth >= maxDepth) return Vec3f(0);

  // カメラ側からの経路構築
  IntersectInfo info;
  if (scene.intersect(ray, info)) {
    // 光源に直接当たった場合は寄与をそのまま返す
    if (info.hitPrimitive->hasAreaLight()) {
      return info.hitPrimitive->Le(info.surfaceInfo, -ray.direction);
    }

    // 交差点のBSDFの種類を取得
    const BxDFType bxdf_type = info.hitPrimitive->getBxDFType();

    // Diffuse面に当たった場合はL_{l,d}, L_{s,d}, L_{d,d}を計算
    if (bxdf_type == BxDFType::DIFFUSE) {
      // 直接照明の計算(L_{l,d})
      const Vec3f Ld =
          computeDirectIllumination(scene, -ray.direction, info, sampler);

      // 集光模様の計算(L_{s,d})
      const Vec3f Lc = computeCausticsWithPhotonMap(-ray.direction, info);

      // 間接照明の計算(L_{i,d})
      const Vec3f Li =
          computeIndirectIllumination(scene, -ray.direction, info, sampler);

      // 3つを足し合わせて終了
      return (Ld + Lc + Li);
    }
    // Specular面に当たった場合は再帰的にレイを追跡する
    else if (bxdf_type == BxDFType::SPECULAR) {
      // BSDFを用いた方向サンプリング
      Vec3f dir;
      float pdf_dir;
      const Vec3f f = info.hitPrimitive->sampleBxDF(
          -ray.direction, info.surfaceInfo, TransportDirection::FROM_CAMERA,
          sampler, dir, pdf_dir);

      // サンプリングした方向に向かうレイを生成
      const Ray next_ray(info.surfaceInfo.position, dir);

      // Throughputの計算
      const Vec3f throughput =
          f *
          cosTerm(-ray.direction, dir, info.surfaceInfo,
                  TransportDirection::FROM_CAMERA) /
          pdf_dir;

      // 再帰的にレイを追跡する
      return throughput *
              integrateRecursive(next_ray, scene, sampler, depth + 1);
    } else {
      spdlog::error("[PhotonMapping] invalid BxDF type");
      return Vec3f(0);
    }
  } else {
    // レイがシーンに当たらなかった場合
    return Vec3f(0);
  }
  return Vec3f(0);
}


// ray_inの方向から来る放射輝度を計算する
Vec3f integrate(const Ray& ray_in, const Scene& scene,
                Sampler& sampler) const override {
  // 深さ0から再帰を始める
  return integrateRecursive(ray_in, scene, sampler, 0);
}

レンダリング

以上の処理を実装してコーネルボックスをレンダリングしてみると以下のようになります。

パラメーター	値
フォトン数	1000000
近傍フォトン数	100
カメラ側からの経路構築のサンプル数	100

誤差が目立たなくなり、とても綺麗にレンダリングできるようになりました。

以降では補足的な説明をしていきます。

集光模様フォトンマップを別に用意する理由

理由は幾つか考えられるのですが、最も分かりやすいのは以下のような実装上の問題です。

理由1: 実装上の問題

今、集光模様も大域フォトンマップを用いて計算することを考えてみます。

放射輝度推定を行っている点を $x$ とします。大域フォトンマップには光源から直接 $x$ に到達するようなフォトンも含まれています。

集光模様の計算では点 $x$ の近傍にあるフォトンを用いて放射輝度の推定を行います。すると光源から直接 $x$ に到達するようなフォトンもこの時計算に入れてしまっていることが分かります。つまり集光模様だけを計算するはずが、直接照明まで計算してしまっている状態になっています。

直接照明の計算は既にNext event estimationで行っているので、この時実質的に 直接照明計算のサンプル数が増えた状態 になっています。そのため最後にモンテカルロ積分の平均を取る時に直接照明だけ異なるサンプル数で割ってあげないと明るすぎる画像が得られてしまいます。

以下の画像は左はパストレーシングで計算したリファレンス、右は集光模様を大域フォトンマップで計算した画像です。リファレンスと比べて明るくなってしまっていることが分かります。

直接照明だけ異なるサンプル数で割らないといけないのは実装上面倒です。そのため集光模様フォトンマップを用意して直接照明の計算を全てNext event estimationに分離し、この面倒な問題を回避していると考えられます。

理由2: レンダリング品質の調整

大域照明フォトンマップと集光模様フォトンマップの構築に使うフォトン数を調整することで、集光模様や間接照明をどれくらい綺麗にレンダリングするかをユーザーが調整することができます。例えば集光模様が支配的なシーンでは集光模様フォトンマップの構築に使うフォトン数を増やすといった使い方が考えられます。

以下は集光模様が支配的なシーンで集光模様フォトンマップの構築に使うフォトン数を10000, 100000, 1000000と変えた場合の比較画像です。フォトン数が増えるほど集光模様の解像度が上がっていく様子が見えます。

このようにユーザーがレンダリング品質を調整しやすくなるという意味でも、集光模様フォトンマップを別に用意するメリットがあると思います。

なぜFinal gatheringで低周波ノイズが減るのか

最後にFinal gatheringで低周波ノイズが減る理由について考えてみたいと思います。

水入りコーネルボックスのシーンを直接照明、集光模様、間接照明のそれぞれを1サンプルで計算して画像にしてみます。

直接照明はNext event estimationで計算しているため、誤差はパストレーシング同様に高周波ノイズとして画像に現れます。このような高周波ノイズはサンプル数を増やしてあげれば目立たなくなってくれます。

しかも直接照明は大抵のシーンで最も見た目への影響が大きいです。以下は元の画像から直接照明、集光模様、間接照明のそれぞれを抜き出して画像にしたものですが、直接照明が一番全体的に明るい画像であることが分かります。

視覚的な影響が最も大きい部分にフォトンマッピング特有の低周波ノイズが入らなくなっているので、最終的な合成画像も綺麗に見えます。

集光模様は一般的に光が狭い領域に集まるため、フォトンマップを用いた推定を行ってもカーネルの半径が縮小して誤差が小さくなる傾向にあります(カーネルの半径を $K$ 近傍で変えている場合)。そのため誤差があまり目立たない画像が得られます。

間接照明は方向をランダムに選んでその方向にレイを飛ばし、レイを飛ばした先でフォトンマップを用いた放射輝度推定を行っているのでした。フォトンマップによる推定を行う前にランダムな方向サンプリングが入ることで、隣接ピクセル間で同じフォトンが参照されることがほぼなくなります。その結果、誤差はパストレーシングのような高周波ノイズとして画像に現れます。このノイズは直接照明同様にサンプル数を増やしてあげれば目立たなくなってくれます。

まとめるとFinal gatheringではNext event estimationやランダムな方向サンプリングによって、フォトンマッピング特有の低周波ノイズを高周波ノイズに変換 してしまうことで画像が綺麗に見えています。

Final gatheringがうまくいかない場合

直接照明と間接照明は大抵のシーンで誤差が目立たない画像が得られます。一方で集光模様は光が一箇所に集まるような場合でないと誤差が目立つ可能性があります。

以下は完全鏡面のあるコーネルボックスのシーンをレンダリングしたものです。左がパストレーシングで計算したリファレンス、右がFinal gatheringで計算したものです。

このシーンだと鏡によって光が拡散されてしまうのでフォトンはシーン全体に散らばります。その結果カーネルの半径が広がり、集光模様フォトンマップを用いた推定の誤差が大きくなります。レンダリング結果にはその誤差が目立つ形で現れています。天井の明るくなっている部分が分かりやすいです。

Final gatheringは万能ではないと言えるでしょう。

発展手法

これまでの3回で紹介してきたフォトンマッピングは90年代の終わりに提案されたものです。現在は以下のような発展手法があります。

Final gathering付きのフォトンマッピングは実装が煩雑であり、良いレンダリング結果を得るためのパラメーター調整も難しいです。また、最も大きいデメリットとして、メモリ的な制約により非常に大きなフォトン数を取れない ということが挙げられます。その結果、大量のフォトンを用意して時間をかけて画像を綺麗にしようとしても、そもそもメモリ的な制約によりフォトンマップが構築できないという問題が起こります。

(Stochastic) progressive photon mappingはこの問題を解決し、メモリ制約の中で大量のフォトンを用いてレンダリングすることが可能になっています。しかも実装はナイーブなフォトンマッピングの実装を少しいじるだけなので簡単です。そのため現在では最初から(Stochastic) progressive photon mappingを実装するのがおすすめです。

これらの発展手法についてもそのうち解説記事を書きます。

References

Jensen, Henrik Wann. Realistic image synthesis using photon mapping. AK Peters/crc Press, 2001.
https://pbr-book.org/3ed-2018/Light_Transport_III_Bidirectional_Methods/Stochastic_Progressive_Photon_Mapping#
Jensen, Henrik Wann. "Global illumination using photon maps." Eurographics workshop on Rendering techniques. Springer, Vienna, 1996.
Veach, Eric. Robust Monte Carlo methods for light transport simulation. Stanford University, 1998.
Hachisuka, Toshiya, Jacopo Pantaleoni, and Henrik Wann Jensen. "A path space extension for robust light transport simulation." ACM Transactions on Graphics (TOG) 31.6 (2012): 1-10.