フォトンマッピングの実装

この記事はフォトンマッピングの理論の続きです。前回はフォトンマッピングの理論的な定式化ついて扱いましたが、今回はその内容を元に実装について解説します。

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

https://github.com/yumcyaWiz/photon_mapping

計算対象の式

まずはフォトンマッピングで計算したい対象の式を確認します。

今、カメラ側からレイを飛ばして $x_0, x_{1}, \cdots, x_t$ という経路を構築したとします。最後の点 $x_t$ においてフォトンマップを用いて放射輝度推定を行うとします。

この時 $j$ 番目の画素への寄与は、前回の式(3)から以下のように計算できます。

$N$ : カメラ側からの経路構築のサンプル数
$N_p$ : 光源から飛ばしたフォトンの数
$x_0^i, x_1^i, \cdots, x_{t(i)}^i$ : $i$ 番目のカメラ側からの経路(経路の長さは $t(i)$ )
$y_0^p, y_1^p, \cdots, y_{s(p)}^p$ : $p$ 番目のフォトンの光源からの経路(経路の長さは $s(p)$ )
$\beta_p$ : $p$ 番目のフォトンのThroughput
$K$ : カーネル関数
$f$ : 点 $x_{t(i)}^i$ におけるBSDF
$T_i$ : $i$ 番目のカメラ側からの経路のThroughput

I_j \approx \frac{1}{N}\sum_{i=1}^N \bigg(\frac{1}{N_p}\sum_{p=1}^{N_p}\beta_p K(y_{s(p)}^p, x_{t(i)}^i)f(y_{s(p)-1}^p \to y_{s(p)}^p, x_{t(i)}^i \to x_{t(i)-1}^i)\bigg)T_i \tag{1}

この式はフォトンのThroughput $\times$ カーネルの値 $\times$ $x_t$ におけるBSDF $\times$ カメラ側からのレイのThroughputの平均になっています。

$\beta_p, T_i$ は以下のようになります。

\beta_p = \frac{L_e(y_0^p \to y_1^p)G(y_0^p, y_1^p)}{p(y_0^p)}\frac{\prod_{k=1}^{s(p)-1} f(y_{k-1}^p \to y_k^p, \to y_{k+1}^p)G(y_{k}^p, y_{k+1}^p)}{\prod_{k=1}^{s(p)} p(y_k^p|y_{k-1}^p)} \tag{2}

T_i = \frac{\prod_{k=1}^{t(i)-1}f(x_{k+1}^i \to x_k^i \to x_{k-1}^i)G(x_{k+1}^i, x_k^i)}{\prod_{k=1}^{t(i)} p(x_k^i|x_{k-1}^i)} \frac{G(x_1^i, x_0^i)W_e^{j}(x_1^i \to x_0^i)}{p(x_0^i)} \tag{3}

複雑に見えますが、基本的には経路の各点(光源とカメラ上の点を除く)でBSDFの値 $\times$ 幾何項 $\div$ その点をサンプリングする確率密度の値を計算して掛け合わせていく形です。 $T_i$ は通常のパストレーシングのThroughputと同じ形になっています。

フォトンマッピングでは式(1)を数値計算していきます。以降ではその実装方法について見ていきます。

フォトンマップの構築

データ構造

式(1)をもう一度見てみます。

I_j \approx \frac{1}{N}\sum_{i=1}^N \bigg(\frac{1}{N_p}\sum_{p=1}^{N_p}\beta_p K(y_{s(p)}^p, x_{t(i)}^i)f(y_{s(p)-1}^p \to y_{s(p)}^p, x_{t(i)}^i \to x_{t(i)-1}^i)\bigg)T_i

これを計算するためにはフォトンのThroughput $\beta_p$ , フォトンが格納された時のフォトンの位置 $y_{s(p)}^p$ , フォトンが格納された点におけるフォトンの入射方向 $y_{s(p)-1}^p \to y_{s(p)}^p$ が必要です。

そのためフォトンマップにはこの3つのデータを保存しておくようにします。フォトンを表現する構造体Photonは以下のように定義すると良いでしょう。

struct Photon {
  Vec3f throughput; // スループット
  Vec3f position; // フォトンの位置
  Vec3f wi; // フォトンの入射方向
};

式(1)の計算では全フォトンに関して総和を取っていますが、前回も述べたように実際にはカーネル関数としては

K(x, y) = \begin{cases} \frac{1}{\pi r^2} & (\|x-y\| \le r) \\ 0 & (otherwise) \end{cases}

のように距離が $r$ 以上離れたフォトンに対してはカーネルの値が0になるようなものを使います。そのため $x_{t(i)}^i$ の近傍にあるフォトンに対してのみ総和を取れば良いです。そのためフォトンマップには近傍探索が効率良く行えるデータ構造を採用します。

近傍探索が効率良く行えるデータ構造としてはKd-treeやOctreeが挙げられます。ここではKd-treeを実装することにします。

Kd-treeの実装

これについては既に良い記事があるのでここでは省略します。以下を参照してください。

私の実装例は以下で見れます。

https://github.com/yumcyaWiz/photon_mapping/blob/b5d25012ce0817e83c260f2a5b5a034ae0891561/include/photon_map.h#L46

フォトントレーシング

Kd-treeを実装できたら、次は光源からシーンに向けてフォトンを飛ばし、フォトンマップにフォトンを格納していく処理を書いていきます。

フォトンの生成

まずは光源からフォトンを生成し、シーン中にフォトンを飛ばす必要があります。以下のような流れでこれを行っていきます。

シーンから光源を1つ選ぶ
選んだ光源上で点 $y_0$ をサンプリングする
方向 $\vec{\omega_{y_0}}$ をサンプリングする
点 $y_0$ から方向 $\vec{\omega_{y_0}}$ にフォトンを飛ばす

光源の選択

今シーンに光源が $N_l$ 個あるとします。1に関しては $N_l$ 個の光源から等確率でランダムに1つの光源を選ぶような処理が簡単でしょう。この場合 $i$ 番目の光源を選ぶ確率 $P_L$ は

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

となります。

光源上の点サンプリング

光源を選んだら光源上で点 $y_0$ をサンプリングします。光源の表面で一様に点をサンプリングすると仮定すると、光源の表面積を $A$ として点サンプリングの確率密度は

p(y_0) = \frac{1}{A} \tag{5}

となります。

光源上の方向サンプリング

続いて方向サンプリングを行います。拡散光源ならコサイン比例半球サンプリングをし、指向性光源の場合は特定の方向だけをサンプリングするようにすると良いでしょう。今、方向サンプリングとしてコサイン比例半球サンプリングを使ったと仮定すると、方向 $\vec{\omega_{y_0}}$ をサンプリングする確率密度は

p(\vec{\omega_{y_0}}) = \frac{1}{\pi}\cos{\theta_{y_0}} \tag{6}

となります。ここで $\theta_{y_0}$ は方向 $\vec{\omega_{y_0}}$ と点 $y_0$ における法線のなす角度です。

フォトンのThroughputの初期化

最後にフォトンのThroughputを初期化してフォトンを飛ばします。この時点でのフォトンのThroughputは式(2)から

\beta_p = \frac{L_e(y_0 \to y_1)G(y_0, y_1)}{P_L(i) p(y_0)p(y_1 | y_0)} \tag{7}

となります。ここで光源をランダムに1つ選択することによって分母に光源を選ぶ確率 $P_L$ が入っていることに注意します。また、点 $x_0$ から方向 $\vec{\omega_0}$ を選んだことによって暗黙的に次の点 $y_1$ が定まっていることに注意します。

式(7)を $\vec{\omega_0}$ を使って書き直すことを考えます。 $y_1$ における法線と $-\vec{\omega_{y_0}}$ のなす角度を $\theta_1$ とすると、幾何項は定義から

G(y_0, y_1) = \frac{\cos{\theta_{y_0}}\cos{\theta_{y_1}}}{\|y_1 - y_0\|^2}

となります。また、 $p(y_1 | y_0)$ は点サンプリングと方向サンプリングの確率密度の変換公式

p(y_1 | y_0) = \frac{G(y_0, y_1)}{\cos\theta_{y_0}}p(\vec{\omega_{y_0}})

を使うことで、式(7)は以下のように書き換えられます。

\beta_p = \frac{L_e(y_0 \to y_1)\cos\theta_{y_0}}{P_L(i)p(y_0)p(\vec{\omega_{y_0}})} \tag{8}

式(4), 式(5), 式(6)を代入すると最終的に以下のような形になります。

\beta_p = N_l A \pi L_e(y_0 \to y_1)

この値でフォトンのThroughputを初期化すれば良いです。

実装

ここまでの処理をコードにすると以下のようになります。

// 光源上から出るレイをサンプリングして返す。この時throughputも初期化する。
Ray sampleRayFromLight(const Scene& scene, Sampler& sampler,
                        Vec3f& throughput) {
  // 光源をシーンからランダムに1つ選ぶ
  float light_choose_pdf;
  const std::shared_ptr<Light> light =
      scene.sampleLight(sampler, light_choose_pdf);

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

  // 方向をサンプリングする
  float light_dir_pdf;
  const Vec3f dir =
      light->sampleDirection(light_surf, sampler, light_dir_pdf);

  // サンプリングした点, 方向を使ってレイを生成
  Ray ray(light_surf.position, dir);

  // Throughputを初期化する
  throughput = light->Le(light_surf, dir) /
                (light_choose_pdf * light_pos_pdf * light_dir_pdf) *
                std::abs(dot(dir, light_surf.shadingNormal));

  return ray;
}

この実装ではThroughputの初期化を式(8)で行っています。

光源側からの経路構築とフォトンの格納

次にフォトンをシーン中に飛ばし、光源側からの経路構築を行っていきます。経路構築の途中でDiffuse面に当たる度にフォトンをフォトンマップに追加していきます。

実装は以下のようになります。ここではフォトンの配列を最初に用意して、そこにフォトンを追加していき、最後にKd-treeを構築するような流れになっています。

また、高速化のためにフォトンの追跡はOpenMPを用いて並列に行っています。throughputは $\beta_p$ に対応します。

std::vector<Photon> photons;

// スレッドごとにSamplerを用意する
std::vector<std::unique_ptr<Sampler>> samplers(omp_get_max_threads());
for (int i = 0; i < samplers.size(); ++i) {
  samplers[i] = sampler.clone();
  samplers[i]->setSeed(samplers[i]->getSeed() * (i + 1));
}

// OpenMPを用いて並列にフォトンを飛ばす
// nPhotonsGlobal: フォトン数
// globalPhotonMap: フォトンマップ
#pragma omp parallel for
for (int i = 0; i < nPhotonsGlobal; ++i) {
  auto& sampler_per_thread = *samplers[omp_get_thread_num()];

  // 光源からレイをサンプリングし、Throughputを初期化する
  Vec3f throughput;
  Ray ray = sampleRayFromLight(scene, sampler_per_thread, throughput);

  // 光源からの経路構築
  for (int k = 0; k < maxDepth; ++k) {
    IntersectInfo info;
    if (scene.intersect(ray, info)) {
      // レイが当たった点のBSDFの種類を取得
      const BxDFType bxdf_type = info.hitPrimitive->getBxDFType();

      // Diffuse面ならフォトンを配列に追加
      if (bxdf_type == BxDFType::DIFFUSE) {
#pragma omp critical
        {
          photons.emplace_back(throughput, info.surfaceInfo.position,
                                -ray.direction);
        }
      }

      // ロシアンルーレットを行って経路構築を続けるか決める
      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のビルド)
globalPhotonMap.setPhotons(photons);
globalPhotonMap.build();

経路構築の部分はパストレーシングと同様のアルゴリズムになっています。ここではロシアンルーレットを行って経路構築を確率的に打ち切ることを行っています。

Throughputの更新は式(2)の後ろの項

\frac{\prod_{k=1}^{s-1} f(y_{k-1}^p \to y_k^p, \to y_{k+1}^p)G(y_{k}^p, y_{k+1}^p)}{\prod_{k=2}^{s} p(y_k^p|y_{k-1}^p)} \tag{9}

に対応しています。元々の式は分母が $k=1$ から始まっていますが、 $k=1$ の場合は既に式(8)で計算されているので、この式ではスキップされています。

ここでは点 $y_k^p$ は点 $y_{k-1}^p$ からBSDFを用いた方向サンプリングによって暗黙的に取っているので、Throughputの初期化の所で述べたのと同様に、点サンプリングと方向サンプリングの確率密度の変換公式

p(y_k^p | y_{k-1}^p) = \frac{G(y_k^p, y_{k-1}^p)}{\cos_{\theta_{y_{k-1}}}}p(\vec{\omega_{y_{k-1}^p}})

を用いることで、式(9)は以下のように書き換えられます。

\frac{\prod_{k=1}^{s-1} f(y_{k-1}^p \to y_k^p, \to y_{k+1}^p)\cos{\theta_{y_k^p}}}{\prod_{k=2}^{s} p(\vec{\omega_{y_{k-1}^p}})}

Throughputの更新の所のコードはこの式を計算しています。

放射輝度の計算

フォトンマップを構築することが出来たので、次はカメラ側からレイを飛ばしてフォトンマップを用いた放射輝度推定を行います。

式(1)をもう一度見てみます。

I_j \approx \frac{1}{N}\sum_{i=1}^N \bigg(\frac{1}{N_p}\sum_{p=1}^{N_p}\beta_p K(y_{s(p)}^p, x_{t(i)}^i)f(y_{s(p)-1}^p \to y_{s(p)}^p, x_{t(i)}^i \to x_{t(i)-1}^i)\bigg)T_i

フォトンマップの構築で $\beta_p$ , $y_{s(p)}^p$ , $y_{s(p)-1}^p \to y_{s(p)}^p$ は既に求まっています。ここで計算する必要があるのは $K(y_{s(p)}^p, x_{t(i)}^i)f(y_{s(p)-1}^p \to y_{s(p)}^p, x_{t(i)}^i \to x_{t(i)-1}^i)T_i$ の部分です。

カメラ側からの経路構築

今、カメラ側からパストレーシングによって経路構築をして経路 $x_0^i, x_1^i, \cdots, x_{t(i)}^i$ を得て、 $x_{t(i)}^i$ においてフォトンマップを用いて放射輝度推定を行うとします。

式(1)から $x_{t(i)}^i$ がDiffuse面の場合は近傍にあるフォトンを求めることで、 $i$ 番目のサンプルの $j$ 番目の画素への寄与は以下で計算できます。

I_j^i = \frac{1}{N_p}\sum_{p=1}^{N_{x_{t(i)}^i}}\beta_p K(y_{s(p)}^p, x_{t(i)}^i)f(y_{s(p)-1}^p \to y_{s(p)}^p, x_{t(i)}^i \to x_{t(i)-1}^i)T_i \tag{10}

ここで $N_{x_{t(i)}^i}$ は $x_{t(i)}^i$ の近傍にある(カーネルの値が0にならないような)フォトンの数です。ここで総和を近傍にあるフォトン数 $N_{x_{t(i)}^i}$ ではなく、フォトンマップ構築時に飛ばしたフォトン数 $N_p$ で割っていることに注意してください。

$x_{t(i)}^i$ がSpecular面の場合はその近くにフォトンはないので、次の方向をサンプリングしてレイを再帰的に追跡し、どこかでDiffuse面に当たったら式(10)を計算するようにします。

実装

まずは式(10)の計算をコードに落とし込みます。 $T_i$ の項は後で掛け合わせてあげることにして、それ以外の項を計算するようなコードを書くと以下のようになります。

// ある点から反射する放射輝度をフォトンマップを用いて計算する
// nPhotonsGlobal: 光源から飛ばしたフォトンの総数(N_p)
// nEstimationGlobal: 何個の近傍フォトンを推定に使うか(K)
// globalPhotonMap: フォトンマップ
Vec3f computeRadianceWithPhotonMap(const Vec3f& wo,
                                    const IntersectInfo& info) const {
  // 近傍にあるK個フォトンのインデックスを取得する(K近傍探索)
  float max_dist2;
  const std::vector<int> photon_indices =
      globalPhotonMap.queryKNearestPhotons(info.surfaceInfo.position,
                                            nEstimationGlobal, max_dist2);

  // 式(10)の計算(T_iの項はまだない)
  Vec3f Lo;
  for (const int photon_idx : photon_indices) {
    const Photon& photon = globalPhotonMap.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 /= (nPhotonsGlobal * PI * max_dist2);
  }
  return Lo;
}

この実装ではカーネル関数として

K(x, y) = \begin{cases} \frac{1}{\pi r^2} & (\|x-y\| \le r) \\ 0 & (otherwise) \end{cases}

を使用し、半径 $r$ は $K$ 個の近傍フォトンの最大距離になるようにしています。つまり半径 $r$ が点 $x_{t(i)}^i$ 周辺のフォトンの集まり具体によって変わる形になっています。Jensenの青い本ではそのような実装になっていたので、それに合わせました。

もちろん半径 $r$ を固定するカーネル関数を使っても良いです。ただしこの場合は飛ばしたフォトン数が少ないと半径 $r$ 内に近傍フォトンがなく、推定した値が0になってしまうという問題が起こります。

一方で半径 $r$ を $K$ 個の近傍フォトンの最大距離にするような実装では、必ず $K$ 個のフォトンの総和を取ることになるのでそういった問題は起こりません。

カメラ側から経路構築をし、式(10)を計算するコードは以下のようになります。以下のコードではthroughputが $T_i$ に対応します。

// ray_inの方向から来る放射輝度をフォトンマップを用いて計算する
Vec3f integrate(const Ray& ray_in, const Scene& scene,
                Sampler& sampler) const {
  Ray ray = ray_in;
  Vec3f throughput(1, 1, 1);

  // カメラ側からの経路構築
  for (int k = 0; k < maxDepth; ++k) {
    IntersectInfo info;
    if (scene.intersect(ray, info)) {
      // 光源に直接当たった場合は光源の放射輝度を直接返す
      if (info.hitPrimitive->hasAreaLight()) {
        return throughput *
                info.hitPrimitive->Le(info.surfaceInfo, -ray.direction);
      }

      // レイが当たった点のBSDFの種類を取得
      const BxDFType bxdf_type = info.hitPrimitive->getBxDFType();

      // Diffuse面に当たった場合にはフォトンマップを用いて式(10)を計算
      if (bxdf_type == BxDFType::DIFFUSE) {
        return throughput *
                computeRadianceWithPhotonMap(-ray.direction, info);
      }
      // Specular面に当たった場合には次の方向をサンプリングして、再帰的にレイを追跡する
      else if (bxdf_type == BxDFType::SPECULAR) {
        // 方向サンプリング
        Vec3f dir;
        float pdf_dir;
        Vec3f f = info.hitPrimitive->sampleBxDF(
            -ray.direction, info.surfaceInfo, TransportDirection::FROM_CAMERA,
            sampler, dir, pdf_dir);

        // Throughputとレイを更新する
        throughput *= f *
                      cosTerm(-ray.direction, dir, info.surfaceInfo,
                              TransportDirection::FROM_CAMERA) /
                      pdf_dir;
        ray = Ray(info.surfaceInfo.position, dir);
      }
    } else {
      // レイがシーンに当たらなかった場合
      break;
    }
  }

  return Vec3f(0);
}

引数のray_inはカメラから生成したレイです。

フォトンマップを用いた放射輝度推定では光源から直接カメラに入ってくるような経路や、光源からSpecular面が続いてカメラに至るような経路を含めることができないので、レイが光源に直接当たった場合には光源の放射輝度をThroughputをかけて返すようにします。

Throughputの更新では光源側からの経路構築の時と同様の議論から

\frac{\prod_{k=1}^{t(i)-1} f(x_{k-1}^i \to x_k^i, \to x_{k+1}^i)\cos{\theta_{x_k^i}}}{\prod_{k=2}^{t(i)} p(\vec{\omega_{x_{k-1}^i}})}

を計算しています。

以上で式(10)を計算できるようになりました。これで $I_j^i$ が求まるので、あとは $i=1,\cdots,N$ に関して平均を取ってあげれば $j$ 番目の画素の値が求まります。

I_j = \frac{1}{N}\sum_{i=1}^N I_j^i

画像を出してみる

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

パラメーター	値
フォトン数( $N_p$ )	1000000
近傍フォトンの数( $K$ )	100
サンプル数( $N$ )	100

フォトン数はかなり多いにも関わらず、謎の模様が出てしまっていることが分かります。これはカーネル関数近似の誤差によるものです。

この模様を取り除くためにはさらに多くのフォトンを使ってカーネルの半径を縮小させていく必要があります。しかしこれ以上多くのフォトンを使うのはメモリ的に厳しくなってきます。

そこでJensenの青い本ではFinal Gatheringという手法を使ってこの模様を取り除くことを行っています。

長くなってきたので今回はここで終わりにしようと思います。次回はFinal Gatheringについて説明していきます。

次回

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.