パストレーサーで直接光part2：サンプリングを組み合わせてノイズを減らす (Next Event Estimation, Multiple Importance Sampling)

#graphics #probability

パストレーサーで直接光を考慮に入れるようにしてみたが、世に言うNext Event Estimationというのはちょっと違うっぽい。

640x480, 1,000サンプル/ピクセル, 1分09秒, M1MacbookAir (8コア)

解説動画：

Next Event Estimation (NEE)

前回は確率密度関数(pdf)の合成というのを使って、直接光のサンプリングか・反射面の特性によるサンプリングか、のどちらか一方を固定の確率で選ぶという方式だった。どうもググった感じだとNext Event Estimationというのはどちらかを選ぶというのではなく、直接光の計算を行った上で従来のトレースも続けるので両方行うっぽい。

前回調べたsmallpt-explicitが同様の方法を採っていた（重み付けが「BRDFサンプリングで光源だった場合には係数を0にする」という違いはある）

Multiple Importance Sampling（MIS）

じゃあNEEを使用して両方トレースする場合どうやって合成するのかというのが疑問だが、それを行うのが多重重点的サンプリング（Multiple Importance Sampling、以降MIS）ということになる。

合成pdfで選択する確率を先に固定で与えるのが問題なので、トレース時に互いの重要度を比較してその率に応じてどちらか片方をサンプルするのかと勝手に想像していたがどうやら違って、両方のサンプリング方法でサンプルした上で、ウェイトをおのおの計算する。

MIS具体例

解説を読むとサンプリング方式がn通りだとかサンプル点がN個とか数式に二重のシグマやら、一般化されていて正直自分には理解できない。なので具体的に最小のケースを当てはめると、

直接光サンプリングとBRDFサンプリングという2つのサンプリング方式を用いる
サンプリング方式それぞれでサンプル点を1つずつピックアップ（計2つ）
サンプル点それぞれに対して、各サンプリング方式でのpdfを算出しそれに従ってウェイトを計算する
- 直接光サンプリングで選んだサンプル1に対して直接光サンプリングでのpdfとBRDFサンプリングでのpdfを算出してウェイト計算
  - 衝突点の面法線が光源方向とは逆向きだったり他の物体に遮られていたら、そもそも直接光の影響は受けない
- BRDFサンプリングで選んだサンプル2に対しても同様
  - 光源に当たった場合のみ、当たらなかった場合は光源サンプリング側のpdfが0となる
ウェイトを加味した上ですべてを合計する

という手順になる。

pdfをどうやって求めるかは、光源サンプリング方式のpdfは光源を正方形とすれば前回の計算式（光源の直接サンプリング）、 BRDFサンプリング方式のpdfは反射面がランバート拡散反射で入射角とのコサインに比例した重点的サンプリングを用いることと積分値が1.0になるように係数を求めると $\cos \theta / \pi$ となる。

ウェイトの算出方法

各サンプルごとに、各サンプリング方式でのpdfが $p_A$ と $p_B$ で、サンプル点がA方式から選ばれたものだとすると、ウェイトを $\frac{p_A^2}{p_A^2 + p_B^2}$ と算出する（パワーヒューリスティック）。

ファイアフライノイズが発生する原因として、寄与度が高い（光源）けど発生する確率密度が低い（拡散反射面でのサンプリング）場合に、pdfで除算すると結果的に稀に大きな値が発生することになり分散が増えノイズとなってしまう。上記のウェイトを使うと他によりよいサンプリング方式があってそっちのpdfが高ければ（光源サンプリングでは半球状よりも限られた範囲になるので高い）低いpdf側のウェイトは低く抑えられるので、値の分散が減らせると期待できる（理論的には分散の期待値が小さくなることから示される）。

ソース

smallptと合成pdfの適用を元にして改造：

// 従来のintersectを修正し、距離上限を指定できるようにする
inline bool intersect(const Ray &r, double &t, int &id, double max_dist=INFINITY){
  size_t n=sizeof(spheres)/sizeof(Sphere); double d; t=max_dist;
  for(size_t i=n;i>0;) if((d=spheres[--i].intersect(r))&&d<t){t=d;id=i;}
  return t<max_dist;
}

// 光源（正方形）
const Rectangle light(Vec(50-16,81.5,81.6-16), Vec(32,0,0), Vec(0,0,32), Vec(12,12,12), Vec());

inline bool is_light_visible(const Ray &r, double max_dist) {
  double t; int id;
  if (intersect(r, t, id, max_dist))
    return false;
  // TODO: 光源が１つのみという前提ではない場合、他の光源によって遮られているかの確認が必要になる
  return true;
}

inline Vec cosine_weighted_random_hemisphere(const Vec& n, unsigned short *Xi) {
  double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2);
  Vec w=n, u=((fabs(w.x)>.1?Vec(0,1,0):Vec(1,0,0))%w).norm(), v=w%u;
  return (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();
}

// 多重重点的サンプリングの重み計算
inline double mis_weight(double pdf_a, double pdf_b) {
  // // バランスヒューリスティック
  // double a = pdf_a;
  // double b = pdf_b;
  // パワーヒューリスティック
  double a = pdf_a * pdf_a;
  double b = pdf_b * pdf_b;
  if (a + b == 0) return 0;
  return a / (a + b);
}

inline double calc_light_pdf(const Vec& dir, double distance, const Vec& light_normal, double light_area) {
  double cos_light = -dir.dot(light_normal);
  if (cos_light <= 0)
    return 0;
  cos_light = fmax(cos_light, 1e-6);  // あまり小さい値はクリップして無限大防止
  return (distance * distance) / (cos_light * light_area);
}

inline double calc_lambert_pdf(const Vec& dir, const Vec& surface_normal) {
  double cos_surface = dir.dot(surface_normal);
  return fmax(cos_surface, 0.0) * M_1_PI;
}

inline Vec radiance(const Ray &r, int depth, unsigned short *Xi) {
  Vec radiance(const Ray &r, int depth, unsigned short *Xi, double &t_out, bool &hit_light_out);
  double t_dummy; bool hit_light_dummy;
  return radiance(r, depth, Xi, t_dummy, hit_light_dummy);
}
Vec radiance(const Ray &r, int depth, unsigned short *Xi, double &t_out, bool &hit_light_out){
  hit_light_out = false;

  double t;                               // distance to intersection
  int id=-1;                              // id of intersected object
  intersect(r, t, id);

  if (double t_light = light.intersect(r); t_light < t) {
    t_out = t_light;
    hit_light_out = r.d.dot(light.n) < 0;
    return hit_light_out ? light.e : Vec(0,0,0);
  }

  t_out = t;
  if (!isfinite(t)) { // if miss, return black
    return Vec();
  }

  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, f=obj.c;
  double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z; // max refl
  if (++depth>5) { if (erand48(Xi)<p) f=f*(1/p); else return obj.e; } //R.R.

  if (obj.refl == DIFF){                  // Ideal DIFFUSE reflection with MIS
    Vec direct;
    // 直接光サンプリング(Next Event Estimation)
    {
      Vec l_dir = light.random_sample(Xi) - x;
      double t_direct = l_dir.length();
      l_dir = l_dir * (1.0 / t_direct);
      double pdf_light1 = calc_light_pdf(l_dir, t_direct, light.n, light.area);
      if (pdf_light1 > 0 && is_light_visible(Ray(x, l_dir), t_direct - 1e-4)) {
        double pdf_brdf1 = calc_lambert_pdf(l_dir, nl);
        double weight1 = mis_weight(pdf_light1, pdf_brdf1);
        double cos_surface = fmax(l_dir.dot(nl), 0.0);
        direct = f.mult(light.e) * (M_1_PI * cos_surface * weight1 / pdf_light1);  // f*M_1_PIがbrdf
      }
    }

    // BRDFサンプリング（ランバート）
    Vec indirect;
    {
      Vec i_dir = cosine_weighted_random_hemisphere(nl, Xi);
      Ray diff_ray(x, i_dir);
      double t_indirect; bool hit_light_indirect;
      Vec incoming_rad = radiance(diff_ray, depth, Xi, t_indirect, hit_light_indirect);
      double pdf_brdf2 = calc_lambert_pdf(i_dir, nl);
      double pdf_light2 = 0;
      if (hit_light_indirect) {  // 光源にヒットした場合にはpdfを算出
        pdf_light2 = calc_light_pdf(i_dir, t_indirect, light.n, light.area);
      } else {
        // 光源にヒットしなかった場合には方向が外れている、
        // またはなにかに遮られていて寄与0なので、光源サンプリングのpdfは0。
        // その場合ウェイトは1になり、そのまま加える形になる。
      }
      double weight2 = mis_weight(pdf_brdf2, pdf_light2);
      indirect = f.mult(incoming_rad) * weight2;  // `(コサイン項/(pdf_brdf2*π))` が打ち消し合うので省略
    }

    return obj.e + direct + indirect;
  }
  // 以降、鏡面反射と透過材質の処理（元ソースに合流）

拡散反射の場合にのみNEE+MISを適用（鏡面反射や透過材質の場合は従来通り）
mis_weight関数で2つのpdfからウェイトを計算
- calc_light_pdf、calc_lambert_pdfでサンプル方向から確率密度を逆算
- mis_weight呼び出しはサンプルごとに行う必要があり、上記では2回の呼び出しをまとめたりはできない
直接光のサンプリングではradianceの再帰呼出ではなく、is_light_visibleでシャドウレイを撃って遮られていなかったら、という具合に簡略化
- 光源は自己発光のみで反射光がない（アルベドが0）という想定
BRDFサンプリングで光源にヒットした場合に判別できるようにするために、radiance関数からの戻り値としてカラー値以外にも交点との距離と光源に当たったかを返せるよう引数に追加したオーバーロード関数に変更
合成pdfのRectangleやVecの変更あたりも持ってくること

比較画像

左：合成pdf、右：NEE+MIS

ノイズ軽減が効いた上で、ファイアフライノイズが相当軽減される！
ピクセルあたりのサンプル数　左：合成pdfは2500spp、右：NEE+MISは1000spp
- NEE+MISでは２方向をトレースする分時間が余計にかかるので、同じくらいの時間になるよう適当にサンプル数を調整

考察・あれこれ

MIS、数値トリック巧妙過ぎ…！
- 実際に理解するにはまず期待値を理解する必要がありそう
光源を小さくするとまだ結構ファイアフライノイズが発生する
- 鏡面反射やコースティクスなど、低pdfでも高寄与なパスがあるため？
- ↑すべて拡散反射面にすると光源が小さくてもノイズは少なくてすむ
ウェイトの計算法は理論的にはバランスヒューリスティックが最適らしい、がパストレみたいに小さな範囲が大きな影響を及ぼす場合にはパワーヒューリスティックは外れ値を抑えるのでそっちの方がよい、とのこと

参考

Demystifying multiple importance sampling | lisyarus blog わかりやすい解説。光源サンプリング、コサインウェイトサンプリングの他に、メタリック反射用のVNDFという3つ使っている。 Veach’97にも言及。
多重重点的サンプリング - Computer Graphics - memoRANDOM
NEEとMISの実装 (MIS編) - MochiMochi3D
パストレーシング / Path Tracing - Speaker Deck
smallpt 毎度ベースのパストレーサー
論文：Robust Monte Carlo Methods for Light Transport Simulation Veach97, MISの論文（長くてムズイ、読めてない）
- Optimally Combining Sampling Techniques for Monte Carlo Rendering Veach95

付録

複数光源の場合

光源が複数あるシーンを扱う場合には方法がいくつか考えられる：

光源の1つのみをランダムに選んでサンプリング
すべての光源をサンプリング

光源1つのみを選ぶ方法では、BRDFサンプリングでその光源に当たった場合にはウェイトを乗算し、別の光源やその他に当たった場合にはウェイトを乗算せずにふつうに影響を加える、といった区別が必要になる。散乱点から見たときに２つの光源が重なっているような配置で奥の光源が選択された場合には、必ず手前の光源にヒットしてしまうので直接光としては影響を受けず、なのでウェイトも下げないことになる。

またはウェイトを構わず適用して、直接光の影響を光源選択の割合でスケール1/p倍してやればよい？その方がノイズが少なくなるように見える。

すべての光源をサンプリングして影響を考慮する方法では、反射面でのサンプリングでどれかの光源に当たった場合にはその対象の光源サンプリングのpdfとの関係から算出したウェイトを乗ずればよい（この場合、光源の個数分の光源サンプリング＋BRDFサンプリングでMISしてることになる）。他の光源サンプリングのpdfを考慮しなくていいのは、光源領域外であるかシャドウレイで遮られるかで必ず対象外(pdf=0)になるからだと思う。

ノイズとしてはすべての光源をサンプリングした場合と光源をランダムに選んだ場合とで大差ないように見えたので、計算負荷としては光源ランダム選択方式でスケール1/pがいいのかも。

MISはちゃんと真の値に収束するのか？

MISがアンバイアスである理由を数式では理解できてないのだけど概念的には、例えば両サンプリング方式で「たまたま」全く同じサンプル点が選ばれたとすれば、その場合ウェイトの合計がぴったり1.0になるようにしているため、真の値に等しくなる。

１回のトレース時に両方のサンプリング方式でまったく同じサンプル点が選ばれることはまずないだろうが、多数の試行を繰り返した場合に全体で見るとマッチするケースが対になるため、平均で均すと真の値に収束する。実際にはpdfが異なるため対応するサンプル同士の出現回数が異なるのだけど、お互いMISウェイト以外の自身のpdfによる除算（単体の重点的サンプリング）によって同じレベルに調整されるのでこちらでも釣り合いが取れている。

光源方向ではない領域では光源サンプリング側のpdfが0となり、間接光サンプリングのウェイトが1でそのままになるので、その領域もバイアスは発生しない。

１サンプルモデル

MISだからといって必ず複数サンプルする必要があるわけではないらしい。複数のサンプリング方式のうちの１つをランダムに選び、それだけを調べるということが可能（Veach97, 9.2.4 The one-sample model）。ランダムに選ぶ代わりにその確率で割ることで結果を持ち上げて辻褄を合わせる。

合成pdfと感覚的には似ていて、MISのウェイトがかかるところだけが異なる。合成pdfよりはノイズは低いが普通のMISよりは劣る。あと確率の割合の決め方がやはり難しい、なので普通のMISを使うのがいいんじゃないかという感想。

疑似コード：

光源サンプリングの割合 = 0.2  // この値はあらかじめ適切に決める必要がある
if (乱数 < 光源サンプリングの割合) {
  direct = 黒
  if (シャドウレイが遮られない) {
    direct = 光源サンプリングの計算（MISウェイト込み）
  }
  return 自己発光 + direct * (1.0 / 光源サンプリングの割合)
}

indirect = BRDFサンプリングの計算（MISウェイト込み）
BRDFサンプリングの割合 = 1.0 - 光源サンプリングの割合
return 自己発光 + indirect * (1.0 / BRDFサンプリングの割合)

シャドウレイが遮られたら光源サンプリングの割合を0にしてBRDFサンプリングに流した方が無駄がなくていいかと思ったが、割合が変わってしまうため結果が異なってしまう
- 確率で選ぶより前にシャドウレイで調べるようにもできるが、それだと計算が半分にはなってない
上記のような用途ではなく、例えばBSDFで拡散・光沢・鏡面の各手法のランダム選択に使えるという話

Next Event Estimation (NEE)

Multiple Importance Sampling（MIS）

MIS具体例

ウェイトの算出方法

ソース

比較画像

考察・あれこれ

参考

付録

複数光源の場合

MISはちゃんと真の値に収束するのか？

１サンプルモデル

関連記事

パストレーサーで直接光計算（重点的サンプリング、合成pdf）

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

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

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

CUDAでAO bench

新着記事