双方向パストレ Part 2：多少の最適化

Part1で一応双方向パストレを動かせるようになったけど、アルゴリズムや数式の理解を中心にしていたためかなり愚直に計算していた。 しかしそれだと無駄が多い。 そこで早くも最適化を適用してみる。

双方向パストレ：シリーズ目次

• 前提内容：MIS、3点形式、経路積分
• Part 1：複数経路でMIS適用
• Part 2：多少の最適化　←本記事

幾何項Gと確率密度の打ち消し

スループットの計算時、拡散反射面の場合に幾何項Gを確率密度で除算するんだけどどちらもほぼ同じ項：距離二乗とコサインを含んでいる （3点形式の式(8.3)、(8.10)を参照）。 しかも重点的サンプリングでコサインに比例した確率密度で反射ベクトルを生成している場合には完全に打ち消される！

これを利用することで、calculateThroughputがだいぶ削れる：

Vec calculateThroughput(const Vertex& vi, const Vertex& vj) {
#if !NO_CANCEL_G_AND_PDF
  // コサイン重み付けサンプリングでは pA = G / pi なので、
  // bsdf * G / pA = bsdf * G / (G / pi) = bsdf * pi となり G が打ち消し合う
  // スペキュラ面では vj.pdf_fwd = vi.rr、拡散反射面でも同様にvi.rrを使う
  double k = vi.refl == DIFF ? M_PI : 1.0;
  return vi.throughput.mult(vi.brdf * (k / vi.rr));
#else
  ...愚直に計算するコード
#endif
}

•  Part 1に比べて、関数内で新たにGを計算する必要も、pdfで除算する必要もなくなった
• Part 1ではvj.pdf_fwdで除算していたのに対し、この修正ではvi.rrでの除算に変更されたことに違和感を覚えるが、pdf_fwdは前の頂点からのロシアンルーレットの確率rrを含めた値なので、rrだけを得たい場合には前の頂点の値を使う必要がある

論文内の解説箇所

Veach論文10.2「数学的公式化」 “重みなしの寄与 \(C_{s,t}^*\)“ 以降に記述がある：

Also, notice that the geometry factors \(G(\bm{x} \leftrightarrow \bm{x^\prime})\) do not appear in the formulas for \(\alpha_i^L\) and \(\alpha_i^E\), because these factors occur in both the numerator and denominator of (10.5) (see the definitions of \(p_i^L\) and \(p_i^E\)).

式(10.5)は：

$$ C_{s,t}^* \equiv \frac{f_j(\bar{x})}{p_{s,t}(\bar{x}_{s,t})}. \tag{10.5} $$

で、\(f_j(\bar{x})\)は式(10.2)：

$$ \begin{align*} f_j(\bar{x}) =& L_e(\bm{x_0} \to \bm{x_1}) G(\bm{x_0} \leftrightarrow \bm{x_1}) W_e^{(j)}(\bm{x_{k-1} \to \bm{x_k}}) \\ & \cdot \prod_{i=1}^{k-1} f_s(\bm{x_{i-1} \to \bm{x_i} \to \bm{x_{i+1}}}) G(\bm{x_i} \leftrightarrow \bm{x_{i+1}}). \tag{10.2} \end{align*} $$

、\(p_{s,t}(\bar{x}_{s,t})\)は式(10.4)：

$$ p_{s,t}(\bar{x}_{s,t}) = p_s^L p_t^E. \tag{10.4} $$

光源部分経路の確率密度 \(p_s^L\)や、視点部分経路の確率密度\(p_t^E\)はその少し前に書かれていて（式番号なし）、

$$ \begin{align*} p_i^L &= P_{\delta^\perp}(\bm{y_{i-2}} \to \bm{y_{i-1}}) G(\bm{y_{i-2}} \leftrightarrow \bm{y_{i-1}}) p_{i-1}^L \quad for \ i \ge 2, \\ p_i^E &= P_{\delta^\perp}(\bm{z_{i-2}} \to \bm{z_{i-1}}) G(\bm{z_{i-2}} \leftrightarrow \bm{z_{i-1}}) p_{i-1}^E \quad for \ i \ge 2. \end{align*} $$

再帰的になっているため全頂点の確率密度が掛け合わされる（細かいので省くが、頂点0,1は光源・視点の開始点と方向成分で別の計算式）。

• \(P_{\delta^\perp}(\bm{y_{i-1}} \to \bm{y_i})\)は\(\bm{y_{i-1}}\)から\(\bm{y_i}\)の方向を選んだ確率密度の投影立体角測度

でこれらを使うと大もとの重みなしの寄与が、光源部分経路と視点部分経路の頂点ごとに事前に計算できて、

$$ C_{s,t}^* = \alpha_s^L c_{s,t} \alpha_t^E, \tag{10.8} $$

となり、\(\alpha_s^L\)や\(\alpha_t^E\)がコードで言うところのスループット、 \(c_{s,t}\)が接続係数（両BRDFと幾何項G）ということですな。

寄与が少ない経路の接続をロシアンルーレットで省く

パストレで一番処理負荷がかかるのは交差判定。 双方向パストレでは部分経路の追跡以外にも接続時にシャドウレイでの判定が必要になる。 「そのシャドウレイで発生する交差判定のコスト、すべての経路接続には正当化できませんよね？」とVeachに詰められる。

「これは仕方なくない？必要経費でしょ」と思うんだけど、そういうところも手を抜かない。 反射回数が多くて寄与が小さい経路ではロシアンルーレットで棄却してやることで、分散は増えてしまうが計算量は減らすことができる。

connectVertices で接続判定をする際にスループットの値を考慮して、閾値以下になったらロシアンルーレットで弾くようにすれば簡単に組み込めるでしょう：

double connectVertices(const Vertex& ys, const Vertex& zt, Vec* out, unsigned short *Xi) {
  Vec d = ys.p - zt.p;
  double dot_ys = -ys.n.dot(d);
  if (dot_ys <= 0)
    return 0;
  double dot_zt = zt.n.dot(d);
  if (dot_zt <= 0)
    return 0;
  double dist2 = d.squareLength();

  // dを正規化してない状態で内積値を計算してるので余分に掛かっているdist^2を除算することで、
  // 幾何因子 G = (cosθo * cosθi') / dist^2 を計算
  double G = (dot_ys * dot_zt) / (dist2 * dist2);

  // 接続ロシアンルーレット
  Vec unoccluded = ys.brdf.mult(zt.brdf).mult(ys.throughput).mult(zt.throughput) * G;
  double lum = unoccluded.max();

  constexpr double THRESHOLD = 0.1;
  double q = lum / THRESHOLD;
  double weight_rr = 1.0;
  if (q < 1) {
    if (erand48(Xi) >= q)
      return 0;
    weight_rr = 1.0 / q;
  }

  // 生き残った場合のみシャドウレイを撃つ
  double dist = sqrt(dist2);
  if (!is_visible(Ray(zt.p, d * (1.0 / dist)), dist))
    return 0;

  *out = unoccluded * weight_rr;
  return G;
}

• 接続の可視判定を行う前にスループットと接続係数で計算した値の大きさから閾値を決めて、ランダムで接続を失敗させる
  • THRESHOLD = 0.1 としてみた
•  Part 1では接続係数を返すだけでスループットは乗じてなかったが今回は内部で計算済み（unoccluded）なのでそれも含めた値を返すように変更、そのため呼出側では乗じる計算を削除する必要あり

論文内の解説箇所

Veach論文10.4「可視テストの数を減らす」以降に記述あり

• 10.4.1「効率性最適化ロシアンルーレット」でどのくらい効果があるのかが指針として解説されるがよく理解できず…

MISの分母計算のループを除去

前回のMISの重みを求める計算で、同じ経路長の各戦略でサンプリングした際の確率密度を求めるには頂点の伸長・縮退の係数を順に掛けていけば求められるというのを使い、それらを足し合わせることで分母を求めた。 ただこの方式だと経路接続時に各頂点に対する計算が必要になってしまっていた。 これを事前に計算しておくことでより効率的にできる。

• この最適化がどこで提唱されたのかググったけどわからなかった…

天下り的に

突然だけど天下り的に次のような数列：

$$ m_i = 1 + m_{i-1} \left(\frac{p(\bm{x_{i-1}} \gets \bm{x_i})}{p(\bm{x_{i-2}} \to \bm{x_{i-1}})}\right)^\beta $$

を考えてみる。 この\(s\)番目の項\(m_s\)と\(p_{s-1}^\beta\)を乗じてみると：

$$ \begin{align*} m_s \cdot p_{s-1}^\beta &= \left(1 + m_{s-1} \left(\frac{p(\bm{x_{s-1}} \gets \bm{x_s})}{p(\bm{x_{s-2}} \to \bm{x_{s-1}})}\right)^\beta \right) p_{s-1}^\beta \\ &= p_{s-1}^\beta + m_{s-1} \left(\frac{p(\bm{x_{s-1}} \gets \bm{x_s})}{p(\bm{x_{s-2}} \to \bm{x_{s-1}})} p_{s-1}\right)^\beta \\ &= p_{s-1}^\beta + m_{s-1} \cdot p_{s-2}^\beta \end{align*} $$

インデックスが1つ減った後ろの項の展開を続けると最終的に：

$$ \begin{align*} m_s \cdot p_{s-1}^\beta &= p_{s-1}^\beta + p_{s-2}^\beta + \cdots + p_1^\beta + p_0^\beta \\ &= \sum_{i=0}^{s-1} p_i^\beta \end{align*} $$

となって\(s\)未満の全戦略での確率密度の\(\beta\)乗の和が得られる。 これを光源・視点部分経路の両方で行うことでMIS重み計算の分母が計算できる。

シーンの追跡時に前の頂点情報を使ってこの\(m_i\)を頂点情報として格納しておき、経路接続時にループせず定数時間で計算できるようになる。

• i番目の項がi-1番目の頂点の伸縮係数を用いる点に注意！

初項の扱い

i=0の初項は特別扱いが必要で、初期値として0または1を入れる： 光源部分経路の場合は1、視点部分経路の場合は0。

この違いはシーンの追跡で到達できるかどうかからきていて、smallpt相当のシーンでは光源は大きさを持っているため追跡できるのに対して、視点はピンホールカメラでありレイの反射で到達しないので、到達の確率密度を0とする必要がある。

スペキュラ面の扱い

分母にはすべてのサンプリング戦略の確率密度を入れていいわけではなく、どちらかでもスペキュラ面の場合には省く必要がある。

その場合は、1の加算をなくすことで取り除くことができる。

接続時の計算

部分経路同士の接続時にどういう計算が必要になるか？ 現在注目している接続エッジ\((s,t)\)から1つずらしたサンプリング戦略\((s-1,t+1)\)と\((s+1,t-1)\)への係数の計算が必要となる。

戦略\((s-1,t+1)\)への係数は、

$$ \frac{p(\bm{x_s} \gets \bm{x_t})}{p(\bm{x_{s-1}} \to \bm{x_s})} $$

戦略\((s+1,t-1)\)への係数は、

$$ \frac{p(\bm{x_s} \to \bm{x_t})}{p(\bm{x_{t}} \gets \bm{x_{t-1}})} $$

で求められる。 双方の分子：\(p(\bm{x_s} \gets \bm{x_t})\)も\(p(\bm{x_s} \to \bm{x_t})\)も現在の接続エッジをまたぐ確率密度：

$$ \frac{\cos\theta_s \cdot \cos\theta_t}{\pi \cdot ||\bm{x_t} - \bm{x_s}||^2} = \frac{G}{\pi} $$

で計算できる（\(\cos\theta_s/\pi\)が拡散反射面でコサイン重み付け重点的サンプリングを行うことに基づいている）。

• 細かくは、ロシアンルーレットの確率も掛ける必要がある

分母はシーンの追跡時に頂点情報として格納済み。

視点の逆向き確率密度を正しい値に修正

経路の伸縮に順・逆方向の確率密度を用いるが、Part 1のコードでは視点の逆向き確率密度に異なる値が入っていたのを修正する必要がある。 というのも現状はピンホールカメラなのでシーンの追跡によって視点に到達することがないため、逆向きの確率密度は0となる。 MISの重み計算時に頂点の属性が DIFF ではないため分母には含まれず、先行の頂点もないため前回の計算結果に影響はなかったが、ここで正しい値である必要がある。

視点の開始頂点の属性を SPEC としていたところを、ピンホールカメラであることを示す項目（PINHOLE_CAM）を追加する。 でシーンの追跡時、1つ前の頂点の逆向き確率密度 pdf_rev を設定する際に PINHOLE_CAM だったら 0 にしてやる。

コード

頂点情報に mis_factor を追加して：

struct Vertex {
  ...
#if !NO_PRECALC_MIS_FACTOR
  double mis_factor = 0;  // MIS 再帰累積係数
#endif
  ...
};

シーンの追跡時に値を更新：

void traceScene(Ray r, unsigned short *Xi, vector<Vertex>* path, int maxDepth=INT_MAX){
  for (int depth = 0; depth < maxDepth; ) {
    ...
    // ここまでにロシアンルーレットの確率を設定しておく（光源だった場合は1.0）
    // 光源にあたった場合でも以下の更新を行なった後に抜けること

    prev.pdf_rev = prev.refl == PINHOLE_CAM ? 0.0 :  // ピンホールカメラには到達しないため、逆向き確率密度は0とする必要がある
                   (v.refl == DIFF ? G / M_PI : 1.0) * v.rr;

    // prev.pdf_rev を設定後にmis_factorを計算する
#if !NO_PRECALC_MIS_FACTOR
    // MIS 累積係数の再帰更新
    double m = prev.pdf_rev / prev.pdf_fwd;
    double delta = (prev.refl == DIFF && v.refl == DIFF) ? 1 : 0;
    v.mis_factor = delta + (m * m) * prev.mis_factor;
#endif
    ...

• 両頂点がDIFFで経路接続できる場合のみdelta=1とする
• m * mで\(\beta=2\)（パワーヒューリスティック）

MISの重み計算で利用：

double calculatePathWeight(const vector<Vertex>& lightSubpath, int s,
                           const vector<Vertex>& eyeSubpath, int t, double G = 0) {
  int k = s + 1 + t; // 経路長（エッジ数）
  if (k <= 1)
    return 1.0;

#if !NO_PRECALC_MIS_FACTOR
  const Vertex& zt = eyeSubpath[t];
  double denom;  // MISの分母：各戦略での接続確率密度（の比）のパワーヒューリスティック値合計
  if (s == -1) {  // 視点部分経路で光源にたどり着いた
    // 経路接続がなく新たなエッジが発生してないため、頂点に計算済みのpdf比で計算できる
    double p0 = zt.pdf_rev / zt.pdf_fwd;  // 光源サンプリングをして残りの頂点と接続した場合の確率密度の増減係数
    denom = (p0 * p0) * zt.mis_factor;
  } else {        // 接続エッジをまたぐ確率密度を計算
    const Vertex& ys = lightSubpath[s];
    double ys_pdf_rev = (G / M_PI) * zt.rr;  // 頂点ysがztからサンプリングされた場合の確率密度
    double ps_1 = ys_pdf_rev / ys.pdf_fwd;   // 現在の接続を光源側に1つにずらす際の確率密度の増減係数
    double zt_pdf_rev = (G / M_PI) * ys.rr;  // 頂点ztがysからサンプリングされた場合の確率密度
    double pt_1 = zt_pdf_rev / zt.pdf_fwd;   // 現在の接続を視点側に1つにずらす際の確率密度の増減係数：pdf_fwd/pdf_rev じゃなく逆になっているのは視点部分経路が逆向きだから
    denom = (ps_1 * ps_1) * ys.mis_factor + (pt_1 * pt_1) * zt.mis_factor;
  }
  return 1.0 / (1 + denom);  // `1+` は現在の接続経路の確率密度比の分
#else
  ...愚直に計算するコード
#endif
}

• 経路接続の場合はG/M_PIにrrを掛けた値を用いる
• 視点部分経路の追跡で光源にたどり着いた場合には s==-1 で呼び出され、その場合には経路の接続ではなく単独の経路で、その場合には頂点ztのpdfから求める
  • 光源に達した場合はRRしてないので省かれている
• どちらも部分経路をループでたどることなくO(1)

他、注意点

• 光源部分経路の始点は mis_factor = 1 で初期化する必要がある （視点部分経路の始点は mis_factor = 0）

結果

各変更を順に施した実行時間（640x480, 100spp）：

ロシアンルーレット以外は微々たる改善でしかないけども、仕組みの理解も進んでまあよしとしよう。

• 前：双方向パストレ Part 1：複数経路でMIS適用