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

#graphics

光輸送問題・レンダリング方程式を体得するため、3点形式で表して解いてみたり、経路積分で描画してみて多少は理解が深まった。これで一応準備が整った、いよいよ双方向パストレーシングをやってみる。

特に、複数経路の結果を効率的に組み合わせる多重重点的サンプリングを用いる。

640x480, 10,000 spp, 30分20秒（8コア）

双方向パストレーシングの概要

経路積分の確認時にも光源・視点の双方からシーンを追跡したので一応双方向といえないこともない、が単に終端頂点同士を接続して経路を1つ生成しただけなのであまり効率はよくなかった。

いわゆる実際の双方向パストレでは途中の頂点同士も結んで複数の経路を一度に生成する。そうすることで1回のサンプリングでより多くの光の経路の情報を取り込めて、結果的により速くノイズを減らすことができる。

とはいえ多数の経路から得られる光の寄与情報をどのようにまとめるかというのが問題になる。多重重点的サンプリングという、重みをつけて合成する。

部分経路同士を接続して多数の経路を生成

光源と視点から追跡した部分経路の終端頂点同士だけではなく、それぞれの各頂点を結んで複数の経路を生成する。光源部分経路側の頂点数が $n_L$、視点部分経路側の頂点数が $n_E$の場合、最大で $n_L \times n_E$ 本の経路ができる。

両頂点のどちらかでもスペキュラ反射面、または頂点間がさえぎられている場合には接続できず無効となる
今回は視点に直接はつなげてない
部分経路追跡で光源にたどり着いてしまうこともありえて、その場合は探索経路によって対処する：
- 光源部分経路：取り除く
- 視点部分経路：その頂点とは接続は行わないでその部分経路自体を1つの経路として使用する（光源までの頂点は接続としても利用する）

重みの計算：多重重点的サンプリング

複数経路の光の寄与をどうやって組み合わせるかというと、多重重点的サンプリング(Multiple Importance Sampling, 以降MIS)を使う。

ただMISを使うということは概要的に把握してたけど、なにをグループにするのかというのはすぐには理解できてなかった。答えとしては、経路長ごとにグループ化する。同じ経路長の中でも光源・視点側の頂点数に自由度があり、それらの間でMISを適用する。

レンダリング方程式を展開して経路長ごとに分かれた積分の項ごとにMISする、というイメージ
一度のサンプルで同じ経路長のものが複数生成されるわけだけどMISというのはそれらの間で計算するわけではなく、個々で別々の計算となるのでお互いを直接考慮する必要はない

おのおのの接続経路に対して、同じ頂点列で光源・視点の頂点数が異なる選ばれ方との間で重みの計算をする。異なる選ばれ方の確率密度が必要になるが、それらをどうやって計算するか？

今回求まった接続経路の確率密度は計算できて、そこから接続箇所を1つずらす確率密度が順々に計算できるという方法で算出する。さらには比だけが必要なことを利用して、実際の確率密度自体を算出しない（！）という巧妙なトリックを用いる。

内容詳細

MISの適用方法：経路長毎

以前に直接光計算でMISを適用したときには、扱う戦略が(1)光源サンプリングと(2)BRDFサンプリングという2つの固定で、サンプルも戦略ごとに1つのみに限定していたのである意味理解しやすかった。双方向パストレでの、様々な長さの経路が多数サンプリングされる状況ではどうやって適用したらいいのか？場合によっては遮蔽で接続できなかったり、スペキュラ面で戦略が無効だったり、状況が複雑そうでひるんでいた。

双方向パストレでのサンプリング戦略というのはなにかというと、光源・視点部分経路の頂点数(s, t)になる。経路接続で一度に複数の経路がサンプルされるわけだけど幸いなことにそれぞれ戦略が異なり各戦略のサンプル経路はたかだか1つであり、同時にサンプルされた他の経路を考慮することなく計算できる。

また同じ経路長で戦略の異なる経路が複数サンプルされるわけだけどそれら同士の間で計算する必要があるわけではなく、個々のサンプルに対して仮に同じ経路（頂点列）だけど違う戦略（違う位置での頂点接続）でサンプルした場合の確率密度を計算して高いものに高い重みを割り振るようにする、というのがMISなのであった。

別戦略での確率密度の計算方法

MISの重み計算にはサンプル候補を発生した戦略での確率密度と、他の戦略でそのサンプルが発生する確率密度を用いる。

双方向パストレでいえば光源から視点までの経路がサンプルとなる（経路積分）。経路が発生する確率密度は、光源部分経路の確率密度と視点部分経路の確率密度の乗算で求められる。部分経路の確率密度は始点の確率密度から、続く各頂点を生成する確率密度の積で求められる。

例えば具体例として光源部分経路の頂点が$\bm{x_0}, \bm{x_1}$（頂点数2）、視点部分経路の頂点が$\bm{x_2}, \bm{x_3}, \bm{x_4}$（頂点数3）の接続経路$\bar{x}$の確率密度は、

$$ p_{2,3}(\bar{x}) = \underbrace{p(\bm{x_0}) p(\bm{x_0} \to \bm{x_1})}_{\text{光源側}} \times \underbrace{p(\bm{x_2} \gets \bm{x_3}) p(\bm{x_3} \gets \bm{x_4}) p(\bm{x_4})}_{\text{視点側}} $$

で計算できる。ここで仮にまったく同じ頂点列が隣の戦略：光源側が1頂点、視点側が4頂点でサンプルされた場合の確率密度を考えると、

$$ p_{1,4}(\bar{x}) = \underbrace{p(\bm{x_0})}_{\text{光源側}} \times \underbrace{p(\bm{x_1} \gets \bm{x_2}) p(\bm{x_2} \gets \bm{x_3}) p(\bm{x_3} \gets \bm{x_4}) p(\bm{x_4})}_{\text{視点側}} $$

と計算できる。これは頂点$\bm{x_1}$が光源側の経路の延長で生成されていたものが、視点側からの延長に切り替わったところだけが異なり、他は共通となる。このことを用いて：

$$ \begin{align*} p_{1,4}(\bar{x}) &= p_{1,4}(\bar{x}) \cdot \frac{p_{2,3}(\bar{x})}{p_{2,3}(\bar{x})} \\ &= \frac{p_{1,4}(\bar{x})}{p_{2,3}(\bar{x})} \cdot p_{2,3}(\bar{x}) \\ &= \frac{p(\bm{x_0}) \times p(\bm{x_1} \gets \bm{x_2}) p(\bm{x_2} \gets \bm{x_3}) p(\bm{x_3} \gets \bm{x_4}) p(\bm{x_4})}{p(\bm{x_0}) p(\bm{x_0} \to \bm{x_1}) \times p(\bm{x_2} \gets \bm{x_3}) p(\bm{x_3} \gets \bm{x_4}) p(\bm{x_4})} \cdot p_{2,3}(\bar{x}) \\ &= \frac{p(\bm{x_1} \gets \bm{x_2})}{p(\bm{x_0} \to \bm{x_1})} \cdot p_{2,3}(\bar{x}) \\ \end{align*} $$

大半の項が約分されて残るのは光源側に含まれていた頂点0から頂点1に進む確率密度を取り消す分母と、視点側に追加する頂点2から頂点1に進める確率密度の分子だけで隣の戦略での確率密度を求めることができる。これを順方向・逆方向に繰り返すことで、同じ経路長の異なる戦略での確率密度を次々に計算することができる。そしてそれらを使ってMISの重みが計算できる。

順・逆、接続確率密度

頂点を進める際の確率密度は以前と同様にシーン追跡時に計算できる。上の例でいえば光源側は$p(\bm{x_0} \to \bm{x_1})$、視点側は$p(\bm{x_2} \gets \bm{x_3})$, $p(\bm{x_3} \gets \bm{x_4})$になる。

掛け合わせていく係数にはこの順方向だけじゃなく、 逆向きの確率密度 も必要になる。光源側でいえば$p(\bm{x_0} \gets \bm{x_1})$、視点側でいえば$p(\bm{x_2} \to \bm{x_3})$, $p(\bm{x_3} \to \bm{x_4})$が必要になる。これらはシーン追跡時の向きとは逆なんだけど、仮に逆順に追跡を行う場合にはどういう確率密度になるかを計算して記録しておく。

シーンの追跡方向と接続経路の方向には注意が必要。部分経路を接続する際に光源から視点に向かう方向を順方向だとみなすことにすると、視点から開始してシーンを追跡した部分経路の「順方向」と「逆方向」は反転することになる。

また接続エッジでの確率密度（$p(\bm{x_1} \to \bm{x_2})$と$p(\bm{x_1} \gets \bm{x_2})$）はシーン追跡時には現れるものではないので、経路接続時に計算することになる。

スペキュラ面の扱い

スペキュラ面の場合、確率密度はデルタ関数なので数式上で考えると$\infty$になってしまう。例えば頂点3だけがスペキュラ面だとすると、頂点3から次や前の頂点をたどる確率密度が$\infty$になる。頂点1と2を接続する確率密度から後方のサンプル戦略の確率密度を順に求める場合を考えると、

接続エッジ	確率密度	経路接続
1↔︎2	$\left(p_{0 \to 1}\right) \times \left(\bm{\delta_{2 \gets 3}} \cdot p_{3 \gets 4} \cdot p_{4 \gets 5} \cdot p_{5 \gets 6}\right)$	可
2↔︎3	$\left(p_{0 \to 1} \cdot p_{1 \to 2}\right) \times \left(p_{3 \gets 4} \cdot p_{4 \gets 5} \cdot p_{5 \gets 6}\right)$	不可
3↔︎4	$\left(p_{0 \to 1} \cdot p_{1 \to 2} \cdot p_{2 \to 3}\right) \times \left(p_{4 \gets 5} \cdot p_{5 \gets 6}\right)$	不可
4↔︎5	$\left(p_{0 \to 1} \cdot p_{1 \to 2} \cdot p_{2 \to 3} \cdot \bm{\delta_{3 \to 4}}\right) \times \left(p_{5 \gets 6}\right)$	可

1↔︎2間から2↔︎3間を求める際に分母に現れる$\delta_{2 \gets 3}$が、3↔︎4間から4↔︎5間を求める際に分子に現れる$\delta_{3 \to 4}$と相殺されることにより、ちゃんと計算できるのでしょう。

2↔︎3間や3↔︎4間の確率密度は有限の値だけど、それらは接続エッジにスペキュラ頂点3を含むため、接続時の係数が0となり影響を与えない。

そして確率密度それ自体を単独で扱うのではなく、BRDFや他の確率密度との比でのみ用いるので、デルタが打ち消されることになるのでコード上では1.0を入れておいて問題ない、ということなんでしょう。

部分経路生成：シーン内の追跡

光源・視点を始点としてシーン内でレイを追跡して部分経路を生成するのは経路積分のときと同じだけど、頂点情報に上記の順・逆確率密度やあとで利用するための要素も記録する。

順方向の確率密度

順方向は経路積分のときと同じで、直前の頂点から新たな頂点をBRDFに即して反射方向を抽出する際の（角度に対する）確率密度を用いて、それと距離・向きによる補正を掛けて面積に対する確率密度を算出する。

直前が拡散反射面の場合：拡散反射時の角度サンプリングにcos重点的サンプリングを用いていて、それと角度のpdfから面積のpdfへの変換（3点形式の式(8.10)）を組み合わせると、幾何項$G/\pi$で計算できる
直前がスペキュラ面の場合：確率密度としては$\infty$となるが、便宜上1.0として記録する（扱う際は事前にスペキュラ面であることを判定した上で利用）

逆方向の確率密度

逆方向も同じ形で計算できる。順方向は直前の頂点の反射特性から今回の頂点を計算したのとは逆に、今回の頂点の反射特性から直前の値を計算する。

ロシアンルーレット

経路積分のときは望みの経路長を得るためにあらかじめ指定の反射回数だけ追跡して途中で抜けるようなことはしてなかったが、今回は先に長さを決めずに・でも無限ループに陥らずに・それでいて真の値に収束させるためにロシアンルーレットを適用する。ロシアンルーレットの確率も記録する。

さらにロシアンルーレットの確率は、順・逆方向の確率密度にもあらかじめ掛け合わせてやる（注意点として、１つ前（逆方向なら次）の頂点からの確率を掛ける）。

Veach論文 10.3.3「部分経路長の決定法」に記述あり

スループット

経路が運ぶ光の計算には、反射面の特性(BRDF)や3点形式の幾何項$G$が関わってくる。部分経路の各頂点で、その頂点までの係数をスループットとして算出する。

経路積分のときは終端同士を結ぶだけなので接続時に計算すればよかったが、双方向パストレでは複数箇所で接続してその各頂点での値は共通なので、何度も同じ計算を避けるために先に計算して頂点情報に記録しておく。

光源にたどり着いた場合

部分経路生成で再度光源にたどり着いた場合には、扱いを注意する必要がある。今回、光源は放射だけで反射はしないということを前提に、追跡を打ち切る。その上で、光源に達した頂点を、

光源部分経路で光源に達した場合：単純に取り除く
視点部分経路で光源に達した場合：接続せずに得られる経路として利用する

MISの重み計算

MISの重み計算で、例えばパワーヒューリスティック（$\beta=2$）で$s$番目の戦略のMIS重みを求める場合、

$$ w_s = \frac{p_s^2}{p_0^2 + p_1^2 + \cdots + p_s^2 + \cdots + p_k^2} $$

で表される（$w_s$の合計が1.0となるような重み配分であることが重要）。

先に述べたように、まず接続した経路の確率密度$p_s$を計算して、それに対して接続位置を１つずらした戦略での確率密度を順に求めて、同じ経路長の全戦略の合計を分母にすることで計算することができる。

そもそもの確率密度の値が不要な件（！）

分母分子を$p_s^2$で除算すれば、

$$ w_s = \frac{1}{(p_0/p_s)^2 + (p_1/p_s)^2 + \cdots + 1 + \cdots + (p_k/p_s)^2} $$

となり、各項は$p_s$との比だけで求められる。各確率密度は$p_s$から順に係数を掛けて求めるので、$p_s$の実際の確率密度を用いる代わりに1.0としてしまっても正しく計算できる、というのが巧妙なトリック。このことによって順・逆方向の確率密度の比だけから計算できる。

スペキュラ面の扱い

接続する経路双方どちらかの頂点がスペキュラ面の場合には接続できない（接続する頂点方向の確率密度は0となるため）。なので分母の累計には蓄積させない。

1つずらした経路の確率密度を計算する際には、頂点に記録する順・逆方向の確率密度を$\infty$ではなく$1.0$としておくことで特別扱いすることなく計算できる。

経路積分のときには頂点がスペキュラ面でつなげない場合の割合を考慮する必要があったことをふまえると、MISでも必要で分母に含めるべきなのではないかと思っちゃうんだけど結局のところその確率密度は0なので追加しないのが正しい（$\infty$の代わりに1.0を入れていることもあり）
Veach論文 10.3.5「スペキュラ表面の扱い」で解説されている

光源・視点自体の頂点の特別処理

追跡の開始点である、光源や視点の頂点は特別な扱いが必要になる。

光源：
- 確率密度は面積の逆数
- スループットは光源が発する放射輝度を確率密度で除算したもので初期化
視点：
- 確率密度を1.0、スループットも1.0とする

またどちらの開始点も反射した光ではないので、他の頂点と異なり接続時に反射係数BRDFを掛けないようにする必要があるが、1.0を入れておくことで場合分けしなくて済むようにする。

コード

経路積分のソースを元に、修正部分を。

Vertex: 頂点情報

struct Vertex {
  Vec p, n;        // 交点、法線
  Refl_t refl;     // 反射タイプ (DIFF, SPEC, REFR)
  Vec brdf;        // BRDF / BSDF 値 (DIFF: アルベド/π, SPEC・REFR: アルベド, 光源・カメラ始点: 1.0)
  double pdf_fwd;  // 前方確率密度（面積PDF、直前頂点のRR確率を含む）
  double pdf_rev;  // 後方確率密度（面積PDF、直後頂点のRR確率を含む）
  Vec throughput = Vec(1);  // 累積スループット
  double rr = 1;   // ロシアンルーレット継続確率
  const Sphere* light = nullptr;  // 衝突した光源オブジェクト
  Vertex(const Vec& p_, const Vec& n_, Refl_t refl_, const Vec& brdf_=Vec(1), double pdf_fwd_=1)
    : p(p_), n(n_), refl(refl_), brdf(brdf_), pdf_fwd(pdf_fwd_), pdf_rev(1) {}
};

brdf: アルベド(Sphereのc)の代わりにBRDFを保持するよう変更
pdf_fwd, pdf_rev: 順方向と逆方向の確率密度
throughput: スループット
rr: ロシアンルーレットの確率
light: 光源の開始点、またはレイの追跡で光源に当たった場合に保持

calculateThroughput: スループットの計算

Vec calculateThroughput(const Vertex& vi, const Vertex& vj) {
  if (vi.refl == DIFF) {
    Vec wi = vj.p - vi.p;
    double dist = wi.length();
    wi = wi * (1 / dist);

    double cos_out = vi.n.dot(wi);
    double cos_in = -vj.n.dot(wi);

    // エッジ (vi, vj) の幾何因子 G(vi, vj)
    double G = (cos_out * cos_in) / (dist * dist);

    // 経路積分形式におけるスループットの更新比: bsdf * (G / pdf)
    return vi.throughput.mult(vi.brdf * (G / vj.pdf_fwd));
  } else {
    // Specular / Refractive (Delta BSDF)
    // デルタ関数の評価により bsdf * (G / pA) = vi.brdf * (1 / pdf)
    return vi.throughput.mult(vi.brdf * (1.0 / vj.pdf_fwd));
  }
}

直前の頂点viから現在の頂点vjへのスループットを計算
- vjの反射特性refl・交点座標p・法線n・順方向確率密度pdf_fwdを参照しているので、あらかじめ設定しておく必要がある
経路積分のupdateThroughputでは確率密度の除算を分離したが、やっぱり計算に入れるようにした
拡散反射面の場合：BRDFと幾何項Gと確率密度で計算
スペキュラの場合：本当ならBRDFも確率密度もデルタ関数で分母分子で打ち消し合う、コードでは$\infty$ではなく1.0を入れていてそのまま計算できる

traceScene: 経路をいったんバッファに記録する

void traceScene(Ray r, unsigned short *Xi, vector<Vertex>* path, int maxDepth=INT_MAX){
  for (int depth = 0; depth < maxDepth; ) {
    double t; int id=0;
    if (!intersect(r, t, id)) break;
    const Sphere &obj = spheres[id];        // the hit object
    Vec x=r.o+r.d*t, n=(x-obj.p).norm(), nl=n.dot(r.d)<0?n:n*-1;
    Vec brdf = obj.c * (obj.refl == DIFF ? 1.0 / M_PI : 1.0);
    auto& v = path->emplace_back(x, nl, obj.refl, brdf);

    // 直前頂点 prev とのエッジ幾何因子 G および面積 PDF の計算
    Vertex& prev = (*path)[path->size() - 2];  // すでに今回の頂点が追加済みなのでさらに-1
    double dist2 = t * t;
    double cos_prev = prev.n.dot(r.d);
    double cos_curr = -nl.dot(r.d);
    double G = dist2 > 1e-8 ? fabs(cos_prev * cos_curr) / dist2 : 0;

    // 前方面積PDF (直前頂点 prev から現在の頂点 v へのサンプリング面積PDF)
    v.pdf_fwd = (prev.refl == DIFF ? G / M_PI : 1.0) * prev.rr;

    // 後方面積PDF (現在の頂点 v から直前頂点 prev への逆サンプリング面積PDF)
    prev.pdf_rev = v.refl == DIFF ? G / M_PI : 1.0;  // 逆方向のロシアンルーレットはまだあとで掛ける

    // スループットの更新
    v.throughput = calculateThroughput(prev, v);

    if (obj.e.x + obj.e.y + obj.e.z > 0) {  // 光源に到達した場合は終了
      v.light = &obj;
      return;
    }
    double maxrefl = obj.c.max();
    v.rr = ++depth > 1 || !maxrefl ? maxrefl : 1.0;
    prev.pdf_rev *= v.rr;
    if (v.rr < 1 && erand48(Xi)>=v.rr) break; //R.R
    if (obj.refl == DIFF) {                     // Ideal DIFFUSE reflection
      Vec d = cosine_weighted_random_hemisphere(nl, Xi);
      r = Ray(x, d);
    } else {
      Vec wi = reflect(r.d, n);
      if (obj.refl == REFR) {
        bool into = n.dot(nl)>0;            // Ray from outside going in?
        double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
        if ((cos2t=1-nnt*nnt*(1-ddn*ddn))>=0) { // Total internal reflection
          Vec tdir = refract(r.d, nl, nnt);
          double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
          double Re=R0+(1-R0)*c*c*c*c*c;
          if (erand48(Xi)>=Re)
            wi = tdir;
        }
      }
      r = Ray(x, wi);
    }
  }
}

順方向の確率密度だけではなく、逆方向の確率密度も計算してやる
- 拡散反射の場合にG / M_PIで計算しているのは、拡散反射でのcos比例の重点的サンプリングによるpdf:$\cos\theta/\pi$と、角度pdfから面積pdfへの変換の係数:$\cos\theta^\prime / ||\bm{x} - \bm{x^\prime}||^2$を組み合わせると、そう計算できるというわけ
スペキュラ面ではBRDFも確率密度も1.0を入れておく（$\infty$ではなく）
1つ前からのロシアンルーレットの確率を掛けることで、この頂点が発生する確率密度となる
calculateThroughput を呼び出してスループットを計算し格納
光源にたどり着いたら抜ける（光源は放射のみで反射はしないという前提）
- 光源に対する逆確率密度は呼び出し側でセットするものとして、ここでは処理しない

部分経路の生成（視点・光源）

void generateEyeSubpath(const Ray& r, unsigned short *Xi, vector<Vertex>* path){
  // ピンホールカメラへの直接接続技法 (S'=k) を除外するため、カメラ頂点は SPEC 扱いとする
  // (S'=k は無効化されるため、カメラ頂点の pdf_fwd/pdf_rev は形式的な初期値 1.0 のままでよい)
  path->emplace_back(r.o, r.d, SPEC);
  traceScene(r, Xi, path);
  auto& last = path->back();
  if (last.light != nullptr) {
    // 光源にたどり着いた場合には、逆向きの確率密度を計算
    last.pdf_rev = 1 / last.light->area();  // 光源上の一様サンプリング面積PDF
  }
}
void generateLightSubpath(unsigned short *Xi, vector<Vertex>* path){
  constexpr int LIGHT_ID = numSpheres - 1;
  const Sphere &light = spheres[LIGHT_ID];
  // 光源の表面をランダムにサンプリング
  Vec y_n = random_unit_direction(Xi);
  Vec y = light.p + y_n * light.rad;
  double pA_light = 1 / light.area(); // 光源上の一様点サンプリング面積PDF

  Vec d = cosine_weighted_random_hemisphere(y_n, Xi);

  // 光源の始点を追加
  auto& v = path->emplace_back(y, y_n, light.refl, Vec(1), pA_light);
  v.throughput = light.e * (1.0 / pA_light);
  v.light = &light;

  traceScene(Ray(y, d), Xi, path);
  auto& last = path->back();
  if (last.light != nullptr) {
    path->pop_back();  // 光源にたどり着いたらそれは取り除く
  }
}

視点部分経路生成：generateEyeSubpath
- 一応開始点も頂点に追加（が接続では利用しない）
  - SPECとすることで経路接続から除外
- 光源にたどり着いた場合には逆向きの確率密度を計算してやる
光源部分経路生成：generateLightSubpath
- 光源である球の表面を一様にサンプリングして開始点とする
  - それゆえ確率密度は面積の逆数
- 初期の追跡方向は法線半球状にcos比例（このことと材質がDIFFであることにより、開始点でも確率密度計算が特別扱いが不要になっている！）
- 初期スループットを光源が発する光light.eと確率密度から計算
- 再度光源にたどり着いた場合には、その最後の頂点は無効

connectVertices: 頂点同士の接続

double connectVertices(const Vertex& ys, const Vertex& zt, Vec* out) {
  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();
  double dist = sqrt(dist2);
  if (!is_visible(Ray(zt.p, d * (1.0 / dist)), dist))
    return 0;

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

  *out = ys.brdf.mult(zt.brdf) * G;
  return G;
}

経路積分のものと計算内容は同じだが、頂点情報に記録したBRDFを利用するよう変更
- 接続するのはランバート面だから方向によらず固定の値だが、Glossyを扱う場合には接続方向と続く頂点方向から求める必要があるだろう
out に接続係数を、戻り値として幾何項Gを返すが0なら接続失敗を表す

calculatePathWeight: 経路の重み計算

ここが今回の目玉！

// 経路の重み計算（MIS: パワーヒューリスティック、β = 2）
// s: lightSubpath 配列のインデックス（頂点数は s + 1。視点直接ヒット時は s = -1 で頂点数 0）
// t: eyeSubpath 配列のインデックス（頂点数は t + 1。カメラ始点は t = 0）
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;

  // 各頂点の前方・後方面積PDFの計算
  struct Node {
    double pdf_fwd;  // 頂点 i-1 から頂点 i へのサンプリング面積PDF（ロシアンルーレット含む）
    double pdf_rev;  // 頂点 i+1 から頂点 i へのサンプリング面積PDF（ロシアンルーレット含む）
    Refl_t refl;     // 反射タイプ
  };
#if __GNUC__
  // ヒープ確保を避けるためスタック上に可変長配列で確保
  // 経路長があまりに長い場合にはスタック溢れの危険性がある。
  Node nodes[k + 1];
#else
  vector<Node> nodes(k + 1);  // vectorでヒープに確保（安全）
#endif

  // 光源部分経路側 (s >= 0 の場合のみ実行)
  for (int i = 0; i <= s; ++i) {
    Node* node = &nodes[i];
    const Vertex* lp = &lightSubpath[i];
    node->pdf_fwd = lp->pdf_fwd;
    node->pdf_rev = lp->pdf_rev;
    node->refl = lp->refl;
  }

  // 視点部分経路側 (s = -1 の場合は 0..t、s >= 0 の場合は s+1..k)
  for (int j = t; j >= 0; --j) {
    Node* node = &nodes[k - j];
    const Vertex* ep = &eyeSubpath[j];
    node->pdf_fwd = ep->pdf_rev;
    node->pdf_rev = ep->pdf_fwd;
    node->refl = ep->refl;
  }

  if (s >= 0) {
    // 接続エッジ (y_s <-> z_t) の幾何項から確率密度を算出
    nodes[s + 1].pdf_fwd = (G / M_PI) * lightSubpath[s].rr;
    nodes[s    ].pdf_rev = (G / M_PI) * eyeSubpath[t].rr;
  }

  // パワーヒューリスティック。R は各技法の PDF 比なので、分母には R^2 を足す。
  double denom = 1.0;  // 本来なら(s, t)経路の確率密度を初期値とするんだけど、比で扱うため実の値を計算する必要がない。

  // 左（逆）方向 (S' <= s)
  double R = 1.0;  // denomの初期値と同じく、比のため1.0で。
  for (int i = s; i >= 0; --i) {
    Node& node = nodes[i];
    R *= node.pdf_rev / node.pdf_fwd;
    // nodes[i] と nodes[i-1] がともに DIFF の場合のみ有効（S' = 0 は光源自身なので常に有効）。
    if (node.refl == DIFF && (i == 0 || nodes[i - 1].refl == DIFF))
      denom += R * R;
  }

  // 右（順）方向 (S' > s)
  R = 1.0;  // 同じく。
  for (int i = s; ++i < k; ) {
    Node& node = nodes[i];
    R *= node.pdf_fwd / node.pdf_rev;
    // nodes[i] と nodes[i+1] がともに DIFF の場合のみ有効（カメラ頂点 nodes[k] は SPEC なので S' = k は自動的に無効）。
    if (node.refl == DIFF && nodes[i + 1].refl == DIFF)
      denom += R * R;
  }

  return 1.0 / denom;
}

光源部分経路のs番目と、視点部分経路のt番目を接続する際の重みを、MISで計算
- sやtは接続頂点のインデックスを表す
- Veach論文ではsやtは頂点数を表し、頂点が0オリジンで$\bm{y_0}~~\bm{y_{s-1}}$と$\bm{z_0}~~\bm{z_{t-1}}$であり、違いがあるため注意！
視点部分経路で光源にたどり着いた場合には s==-1 として呼び出され、経路の接続は行わない（lightSubpathは使われずにeyeSubpathのみ）
2つの部分経路を組み合わせた1つの経路をインデックスでアクセスできるようnodesに必要な内容をコピー
- 光源から視点に進む方向を順方向として、視点部分経路で探索した方向は逆にする
- ヒープ確保を避けるために可変長配列を用いてみた（※C++標準ではない）
  - 手元の環境ではトータルで2~3割高速化したので相当効果あり
接続エッジの位置(s)をまたぐ際の確率密度を幾何項Gから算出
- connectVerticesでGを計算するので、それを引数で与えてもらい再計算を省略
- 反映させるのは遷移先なのでfwdがs+1、revがsとなる
現在の戦略(s,t)に対して、他の戦略の確率密度の比を計算して、合計の逆数が重みとなる

radiance: ピクセルの色計算

Vec radiance(const Ray &primaryRay, unsigned short *Xi, int pixelindex){
  vector<Vertex> eyeSubpath;
  generateEyeSubpath(primaryRay, Xi, &eyeSubpath);
  vector<Vertex> lightSubpath;
  generateLightSubpath(Xi, &lightSubpath);

  int nl = lightSubpath.size();
  int ne = eyeSubpath.size();
  Vec color(0);

  // --- 1. 視点部分経路で光源に到達した場合 (s = -1 相当) ---
  if (ne > 1) {
    const Vertex& zt = eyeSubpath.back();
    if (zt.light != nullptr) {
      int t = ne - 1;
      Vec c = zt.throughput.mult(zt.light->e);
      double weight = calculatePathWeight(lightSubpath, -1, eyeSubpath, t);
      color = color + c * weight;
#if REGISTER_CONTRIBUTION
      registerContribution(c, 0, t - 1, pixelindex, weight);
#endif
      // 後続の接続ループ (s >= 0) で光源頂点が接続対象にならないように取り除く
      eyeSubpath.pop_back();
      --ne;
    }
  }

  // --- 2. 部分経路同士の頂点をつなげる(拡散反射面 DIFF 同士の場合のみ) ---
  for (int s = 0; s < nl; ++s) {
    const Vertex& ys = lightSubpath[s];
    if (ys.refl != DIFF)
      continue;

    for (int t = 1; t < ne; ++t) {  // 視点の開始点[0]とは接続させない
      const Vertex& zt = eyeSubpath[t];
      if (zt.refl != DIFF)
        continue;

      Vec connect;
      if (double G = connectVertices(ys, zt, &connect); G > 0) {
        double weight = calculatePathWeight(lightSubpath, s, eyeSubpath, t, G);
        Vec c = connect.mult(ys.throughput).mult(zt.throughput);
        color = color + c * weight;
#if REGISTER_CONTRIBUTION
        registerContribution(c, s + 1, t - 1, pixelindex, weight);
#endif
      }
    }
  }

  return color;
}

光源・視点部分経路を生成し、それぞれを接続・MIS重み計算して組み合わせる
sが光源側頂点のインデックス、tが視点側頂点のインデックス
- 視点部分経路のループはtの開始が1からで、視点自体(t = 0)にはつなげてない
視点部分経路で光源に達した場合にはループ前に判定・処理する
- この経路だけが、経路接続しない＝拡散反射が連続する箇所がない経路を扱える唯一の手法
  - そのような経路は接続では発生しないので、MIS重みは1.0となり100%の影響が加えられる
- 光源に達した頂点自体は、接続対象からは取り除く（光源は放射のみで、反射しない（アルベド=0）という前提）

結果

左：パストレ＋直接光（smallpt_explicit）、右：双方向パストレ

解像度320x240で、ある程度同じ実行時間となるようにピクセルあたりのサンプル数を指定
- パストレ＋直接光：1,300spp、26.4秒
- 双方向パストレ：500spp、24.3秒

この比較を見ると直接光だけでも十分…ということになってしまう。もっと光源がたどり辛い間接照明のようなシーンだと差が出るだろうということで、同条件で光源を左下隅に移動させたシーン：

ファイアフライノイズが目立つとはいえ双方向からの追跡が効果を発揮している、と言えるでしょう💦

サンプリング戦略ごとの寄与（重みあり）

光源・視点部分経路の頂点数に応じて寄与を記録することで、どのサンプリング戦略がどういう画像となるかを確認できる
- Veach論文の図10.4相当
左上が(s=-1,t=1)（光源部分経路の頂点数0, 視点部分経路の頂点数2）で、右方向が光源側sの増加、下方向が視点側tの増加
- s=-1（１列目）は経路接続がない、ナイーブなパストレを表す（ただしMIS重みを掛けているのでスペキュラ経路以外は暗い）
- t=-1というのはピンホールカメラなので発生しない（視点に対して交差判定を行っていないので光源部分経路追跡で視点に達することがない）
- 視点自体と経路接続すればt=0というのはありえる、ので本来ならさらに上に1行列存在する（が、本コードでは行なってないため存在しない）
斜め（／）の組が同じ経路長で、同じグループ内でも寄与は全然違う（MIS重みが掛かっていることもあり）

後書き、考察

シーンの構成・光源のレイアウトによらずに、間接光やスペキュラ反射を含むシーンでも効果的に寄与が得られて素晴らしい
理論も数式もコード化もムズイのなんの…
- コード上の計算式だけじゃなく、分岐や計算順とかも関わってくるので１行たりとも気が抜けない
- 数式的にもデルタが関わってたり打ち消し合うトリックでコード上からは読み解けなかったりしてムズイ
Veach論文の結果画像は500x500x25sppとか、30年弱経た今から思うと慎ましくて涙が出そう（計算能力の飛躍的向上に感謝！）

経路接続で同経路長の全戦略が列挙されるわけではない件

光源・視点部分経路同士を接続することで経路長が1から $k (= n_L + n_E - 1)$ までの経路が複数生成される。しかし経路長が$n_L$または$n_E$より長いものには列挙されないサンプリング戦略が出てしまう。そのような状態でMISを計算するにはどうするんだろうか？同じ長さの全戦略を用いていいのか、それとも列挙できた戦略だけに絞る必要があるんだろうか？というのが疑問だった。

結果としては気にせず全戦略を対象にしてしまっていい（むしろしないとダメ）っぽい。数学的な理論はわからないが、ロシアンルーレットで打ち切る確率で補正されるので、それでなんかうまく辻褄が合うんだろうと勝手に推測。直感では光源・視点の矩形で組み合わせるので中央付近の戦略が多く出てしまうんじゃないかと懸念するが、光源・視点のどちらかが極端に短い、長方形のような分布が選ばれれば端の戦略が登場して中央の戦略は登場しないというようなケースもあり得て結果的にUnbiasedになるんちゃう、知らんけど。

Veach論文 10.1「概要」 効率的なサンプルの生成 の「There is an important detail that we have not mentioned yet.」以降に書かれている寄与0の 空経路 $\epsilon$ と再解釈して無限にサンプルしたことにする、という記述が対応するっぽい（数式は9章）。

経路接続がランバート面を前提としている件

現状は材質として完全な拡散反射、鏡面反射、屈折のみ扱っている。そして経路接続できるのは拡散反射のみで、それを前提に確率密度を計算している。また拡散反射面での反射方向のサンプリングがcos重み付きで、そのことも前提になっている。

Glossyな材質も組み込む場合にはどういう修正が必要になるのか、汎用的にできるのか、まだ想像がつかない。

視点接続していないデメリットは小さくない

双方向パストレはコースティクスが得意というのがメリットでもあるんだけど、このコードでは視点に接続してないためその利点は得られてない。光源から発してガラス球で集光されて拡散面である地面に高い確率密度・スループットで得られた頂点が、視点と直接接続できれば効率的に影響が得られるところを、視点側に拡散面を挟まなければいけないため無駄になってしまっている。

またガラス球だけじゃなく鏡の球体によって反射された光がシーン全体に散らばる経路も拾えておらず、それがファイアフライノイズの原因となってしまっている。

サポートするには今までナイーブなパストレと同じでピクセルごとに計算できたところを、ライトトレーシングと同じようにイメージセンサー全体への影響を考慮する必要が出てしまってまたさらに難しい計算になりそう。今後の課題。

リンク

Eric Veach、”Robust Monte Carlo Methods for Light Transport Simulation“、10章双方向パストレーシング至高の聖典
双方向パストレーシング - Computer Graphics - memoRANDOM、実装　日本語で貴重な解説（理解してから読むと理解できる、してない段階ではできないというジレンマ）
edubpt　貴重な解説スライドとソースあり
PBRT: Bidirectional Path Tracing 丁寧に解説されているんだとは思うが、この章までの項目を追えてない自分にとってはソースを見てもチンプンカンプンなところが辛い

接続エッジ	確率密度	経路接続
1↔︎2	\(\left(p_{0 \to 1}\right) \times \left(\bm{\delta_{2 \gets 3}} \cdot p_{3 \gets 4} \cdot p_{4 \gets 5} \cdot p_{5 \gets 6}\right)\)	可
2↔︎3	\(\left(p_{0 \to 1} \cdot p_{1 \to 2}\right) \times \left(p_{3 \gets 4} \cdot p_{4 \gets 5} \cdot p_{5 \gets 6}\right)\)	不可
3↔︎4	\(\left(p_{0 \to 1} \cdot p_{1 \to 2} \cdot p_{2 \to 3}\right) \times \left(p_{4 \gets 5} \cdot p_{5 \gets 6}\right)\)	不可
4↔︎5	\(\left(p_{0 \to 1} \cdot p_{1 \to 2} \cdot p_{2 \to 3} \cdot \bm{\delta_{3 \to 4}}\right) \times \left(p_{5 \gets 6}\right)\)	可

双方向パストレーシングの概要

部分経路同士を接続して多数の経路を生成

重みの計算：多重重点的サンプリング

内容詳細

MISの適用方法：経路長毎

別戦略での確率密度の計算方法

順・逆、接続確率密度

スペキュラ面の扱い

部分経路生成：シーン内の追跡

順方向の確率密度

逆方向の確率密度

ロシアンルーレット

スループット

光源にたどり着いた場合

MISの重み計算

そもそもの確率密度の値が不要な件（！）

スペキュラ面の扱い

光源・視点自体の頂点の特別処理

コード

Vertex: 頂点情報

calculateThroughput: スループットの計算

traceScene: 経路をいったんバッファに記録する

部分経路の生成（視点・光源）

connectVertices: 頂点同士の接続

calculatePathWeight: 経路の重み計算

radiance: ピクセルの色計算

結果

サンプリング戦略ごとの寄与（重みあり）

後書き、考察

経路接続で同経路長の全戦略が列挙されるわけではない件

経路接続がランバート面を前提としている件

視点接続していないデメリットは小さくない

リンク

関連記事

Processing(Java)で.plyファイル読み込み

Processing(Java)で.objファイル読み込み

AABBTreeを組み込んで、レイトレでポリゴン描画

PBRマテリアルの異方性反射

CUDAでAO bench

新着記事