フォトンマッピングでファイナルギャザリングを適用

2026-01-04
#graphics

相当昔にフォトンマッピングを動かしたことがあるが、放射するフォトン数が少ないと斑点ノイズが発生してしまっていた。 これを低減させるファイナルギャザリングを実装してみる。

フォトンマッピングとは？

光源からフォトンを飛ばして交点を記録したフォトンマップ（以降PM）を作成し、画像生成時の各ピクセルの色を求める際に利用する。

（双方向）パストレーサーに比べての利点は、集光模様（コースティクス）を綺麗に描画できることでしょう。

単純なフォトンマッピングの問題点

ピクセルの色を求める際にカメラからの交点位置に対してPMから輝度推定すると、推定には一定量のフォトンを近傍から探すためフォトン数が少ないと斑点ノイズが発生してしまう：

フォトンマッピングで発生する斑点ノイズ

この問題に対してファイナルギャザリングで対応する。

ファイナルギャザリングとは？

正直「ファイナルギャザリング」というのがなにを指すのかいまいちはっきりしないけど言葉的には「最後にかき集める」ということで、 輝度計算をPMだけで行うのではなく直接照明＋間接照明＋集光模様に分けて計算することにして、 間接照明を周りからかき集めることから命名してるのだと思う。

書籍「フォトンマッピング」での解説

フォトンマッピングの開祖たるHenrik Wann Jensen氏の書かれた書籍の和訳本「フォトンマッピング」を読み返してみた （本にはファイナルギャザリング（以降FG）という単語自体は出てこないようで、索引にもない）。

９章「実践的な２段階アルゴリズム」でレンダリング方程式を、

  1. 直接照明
  2. 鏡面反射
  3. 集光模様
  4. 間接照明

の４つに分けて計算することが提示される。

9.3 第１段階：フォトンの追跡

フォトンの記録を集光模様PMと大域PMに分ける。

  • 集光模様PMには光源から放射され鏡面反射または透過してから拡散面に到達したフォトンを格納　（光伝達の表記法：LS+D）
  • 大域PMには拡散面にぶつかるすべてのフォトンを格納、ということで集光模様PMに格納されるフォトンも格納する　（L(S|D)*D）

9.4 第２段階：描画

４つに分けた各計算法が記述されている（が説明がいまいちわかりづらい…）。

  • 直接照明：影光線を使って直接計算
  • 鏡面反射：反射先で計算
  • 集光模様：集光模様PMで推定
  • 拡散反射：何本もサンプル光線を利用して間接照明を正確に計算、サンプル光線での近似計算には大域PMから推定

この拡散反射を正確に計算するための分散光線追跡法（Distributed Ray Tracing、半球状を多数サンプリング）がFGということになるだろう。

すべての経路が追跡されるかの確認

この方法で正しいかどうか、あらゆる光伝達経路が網羅されるかを確認する必要がある。

拡散面での計算は直接照明＋間接照明＋集光模様PMで、直接照明はLD、集光模様PMはLS+Dとなる。

間接照明はFGでDS*D、FG先のPM推定でL(S|D)*Dということで、つなげるとL(S|D)*DS*Dになる（PM推定最後のDとFG最初のDは同じもののため１つにまとめる）。

これらを視点から辿る逆向きの経路DS*Eとつなげると、

タイプ 光経路 備考
直接照明 LDS*E 拡散面1回
集光模様 LS+DS*E 鏡面反射後に拡散面1回
間接照明 L(S|D)*DS*DS*E 拡散面2回以上

という具合で、任意の経路L(S|D)*Eが必ずどれか１つだけに当てはまるようになっている。

  • 今回は点光源のみのため視点からのトレースで光源にはぶつからないので、光源を直接見る経路LS*Eは存在しないので扱ってない
    • 光源が大きさを持ち拡散面を経ずに到達する場合、光源の色を返すようにする必要がある
    • 同じくFGでS経由で光源にぶつかる場合は黒を返す必要がある（そのような経路は集光模様PMに含まれるため）

実装

てなわけでみっちり実装：

int estimate=10, estimate_caustic=5;
std::vector<Photon*> g_photons, g_caustic_photons; KDTree g_kdtree, g_caustic_kdtree;
Vec estimate_flux(KDTree &tree, const Vec& x, const Vec& albedo, int count) {
  unsigned char* buf = (unsigned char*)(alloca(sizeof(Photon*)*count + sizeof(double)*count));
  Photon** b=(Photon**)(buf); double* rr=(double*)(buf + sizeof(Photon*)*count);
  tree.nearest(b,rr,count,x,1e20); Vec flux; const double fr=1.0/M_PI;
  int j,k; for (j=k=0; j<count; k=j++) { Photon* q=b[j]; if (!q) break;
    flux=flux+q->fl*fr; } return flux.mul(albedo)*(1.0/(M_PI*rr[k])); }
enum RayType { PHOTON_EMITTED, PHOTON_CAUSTIC, PHOTON_DIFFUSED, EYE_GATHER, EYE_ESTIMATE };
Vec trace(const Ray &r, int dpt, RayType raytype, const Vec &fl, int i);
inline Vec direct(const Vec& x, const Vec& n, const Vec& albedo) {
  Vec d = LightPos - x; double dist2 = d.sqlen(), dist = sqrt(dist2); d = d.norm();
  double t; int id; Ray shadowRay(x, d);
  if (intersect(shadowRay, t, id) && t < dist - 1e-4) return Vec();
  Vec brdf = albedo * M_1_PI; double NoD = std::max(0.0, n.dot(d));
  return LightPow.mul(brdf) * (NoD / (4*M_PI*dist2));
}
inline Vec indirect(const Vec& x, const Vec& n, const Vec& albedo, int i) {
  Vec w = n, u = ((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm(), v = w%u;
  double r1=2*M_PI*hal(0, i), r2=hal(1, i), r2s=sqrt(r2);
  Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();
  Vec flux = trace(Ray(x+d*1e-4, d), 1, EYE_ESTIMATE, Vec(), i);
  return flux.mul(albedo);
}
Vec trace(const Ray &r, int dpt, RayType raytype, const Vec &fl, int i) {
  double t;int id,d3=++dpt*3; if((dpt>=20)||!intersect(r,t,id)) return Vec();
  const Sphere&obj=sph[id]; const Vec x=r.o+r.d*t,n=(x-obj.p).norm(),&f=obj.c;
  Vec nl=n.dot(r.d)<0?n:n*-1; if (obj.refl == DIFF) {
    if (raytype==EYE_GATHER) return direct(x,nl,f) + indirect(x,nl,f,i) + estimate_flux(g_caustic_kdtree,x,f,estimate_caustic);
    if (raytype==EYE_ESTIMATE) return estimate_flux(g_kdtree,x,f,estimate);
    Photon* q=new Photon; q->x=x; q->fl=fl; double p=fmax(f.x,fmax(f.y,f.z));
#pragma omp critical
    {
      g_photons.push_back(q);
      if (raytype == PHOTON_CAUSTIC) g_caustic_photons.push_back(q);
    }
    if (hal(d3+1,i)>=p)return Vec();
    double r1=2.*M_PI*hal(d3-1,i), r2=hal(d3+0,i), r2s=sqrt(r2);
    Vec &w=nl, u=((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm(), v=w%u;
    Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2));
    return trace(Ray(x,d),dpt,PHOTON_DIFFUSED,f.mul(fl)*(1./p),i);
  }
  if (raytype == PHOTON_EMITTED) raytype = PHOTON_CAUSTIC;
  if (obj.refl == SPEC) {
    return trace(Ray(x,r.d-n*2*n.dot(r.d)),dpt,raytype,f.mul(fl),i).mul(f);
  } else {Ray lr(x,r.d-n*2.0*n.dot(r.d)); bool into = (n.dot(nl)>0.0);
    double nc = 1.0, nt=1.5, nnt = into?nc/nt:nt/nc, ddn = r.d.dot(nl), cos2t;
    if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0) return trace(lr,dpt,raytype,fl,i);
    Vec td = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
    double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:td.dot(n));
    double Re=R0+(1-R0)*c*c*c*c*c,P=Re;Ray rr(x,td);Vec rfl=f.mul(fl);
    return raytype>=EYE_GATHER ? (trace(lr,dpt,raytype,fl,i)*Re + trace(rr,dpt,raytype,fl,i)*(1-Re)).mul(f)
                                : hal(d3-1,i)<P ?trace(lr,dpt,raytype,rfl,i) :trace(rr,dpt,raytype,rfl,i);}}
inline uint32_t hash2d(unsigned int x, unsigned int y) {return (x * 0x3504f333) ^ (y * 0xf1bbcdcb) * 741103597;}
int main(int argc, char *argv[]) { int w=1024, h=768, samps=1000, FG=20;
  argc>1&&(samps=atoi(argv[1])); argc>2&&(FG=atoi(argv[2]));
  argc>3&&(estimate=atoi(argv[3])); argc>4&&(estimate_caustic=atoi(argv[4]));
#pragma omp parallel for schedule(dynamic, 1)
  for(int i=0; i<samps; ++i) { Ray r; Vec f;
    if (omp_get_thread_num() == 0)
      fprintf(stderr,"\rPhotonPass %5.2f%%",100.*i/samps);
    genp(&r,&f,i); trace(r,0,PHOTON_EMITTED,f*(1.0/samps),i+1); }
  fprintf(stderr, "\n#photons:%zu, #caustic:%zu\nBuilding kd-tree...\n", g_photons.size(), g_caustic_photons.size());
  if (!g_photons.empty()) g_kdtree.build(&g_photons[0], g_photons.size());
  if (!g_caustic_photons.empty()) g_caustic_kdtree.build(&g_caustic_photons[0], g_caustic_photons.size());
  Ray cam(Vec(50,52,295.6), Vec(0,-0.042612,-1).norm());
  Vec cx=Vec(w*.5135/h), cy=(cx%cam.d).norm()*.5135, *cc=new Vec[w*h];
#pragma omp parallel for schedule(dynamic, 1)
  for (int y=0; y<h; ++y) {
    if (omp_get_thread_num() == 0)
      fprintf(stderr, "\rHitPointPass %5.2f%%", 100.0*y/h);
    for (int x=0; x<w; ++x) { Vec r; int i=y*w+x;
      for (int fg=0; fg<FG||fg==0; ++fg) { int j=hash2d(x,y)+fg+1;
        double dx=hal(0,j), dy=hal(1,j+1);
        Vec d = (cx * ((x+dx)/w-.5) + cy * (-(y+dy)/h+.5) + cam.d).norm();
        r=r+trace(Ray(cam.o+d*140,d),0,FG>0?EYE_GATHER:EYE_ESTIMATE,Vec(),j);}
      cc[i]=r*(1.0/std::max(FG,1));}}
  FILE* f = fopen("image.ppm","w"); fprintf(f,"P3\n%d %d\n%d\n",w,h,255);
  for (int i=0; i<w*h; ++i) { const Vec& c = cc[i];
    fprintf(f,"%d %d %d ", toInt(c.x), toInt(c.y), toInt(c.z));}}
  • 前回との共通部分は以前のコードから持ってくること
  • FGのメイン部分：trace関数内、拡散反射面の箇所：
    • カメラパスでレイのタイプをEYE_GATHERとEYE_ESTIMATEで判別できるようにしておき、
    • GATHERの場合は直接照明＋間接照明＋集光模様PM
      • 間接照明にESTIMATEレイを撃つ
    • ESTIMATEの場合は大域PMから求める
  • 直接照明directではLightPowに対してBRDFである1.0/PIとコサイン項と、距離減衰1.0/(4*PI*dist2)を乗ずる
  • PMから放射輝度を推定する際には放射したフォトン数で除算する必要がある
    • 以前は最後に一括して除算していたが、今回の計算では直接照明の結果には適用してはいけないので、あらかじめ光源が放射するフォトンの放射輝度を除算してしまっている
  • 間接照明以外にもアンチエイリアス・Glossy光沢反射・半影なども扱うことも考えて、 間接照明計算でループするのではなくメイン側でループさせてみる
  • 元々smallppmを元にしていたので乱数は使わずにhalハルトン数列を利用している
    • フォトンの放射時にはフォトン番号によって偏りなくバラけてくれるが、今回カメラパスで単純にピクセル位置＋FG回数だとパターンが見えてしまった 　→ 2D座標に対する適当なハッシュ関数で回避
    • そのため正の値となるようhal関数の引数jをunsigned intにする
  • その他、改善項目：
    • KDTree構築時に完全にソートする必要はないためnth_elementが使える（O(n)だとか）
    • log_2の計算に__builtin_clzが使える（GCC拡張。C++にはcountl_zeroがあるがC++20が必要なのとunsignedじゃないと怒られる）
  • 集光模様PMのみを可視化した結果：ガラス球の屈折だけじゃなく鏡面球による散らばりでフォトンが散乱したものがマップされるため、推定に用いる個数を少なくすると壁や天井にFG適用前と同様の斑点ノイズが発生するのが盲点だった （ただし直接光に比べれば影響は少ないので合成すれば目立たない）
集光模様マップ

リンク

関連記事

新着記事