smallpm - Global Illumination in 128 lines of C++

2009-10-15
#cg #gi

Smallpm

Proce55ingでフォトンマッピングを試してみたけど、20万フォトンくらいの放射でメモリの限界になってしまい画質があまりあげられなかったので、C++に置き換えてみる。

Beeさんのsmallppmを元に作る。これなら何か間違ったときに確認できるのでいいね。

…なんとかsmallppmと同じ79文字×128行に詰め込んだ。

Here is the list of features (majority of them are from smallppm and smallpt):

  • Global illumination via Photon Mapping
  • 128 lines of 79-column (or less) open source C++ code
  • Point light source
  • Specular, Diffuse, and Glass BRDFs
  • Ray-sphere intersection
  • Modified Cornell box scene description contains LSDSE path
  • Cosine importance sampling of the hemisphere for diffuse reflection
  • Russian roulette for path termination
  • Russian roulette and splitting for selecting reflection and/or refraction for glass BRDF
  • Quasi Monte Carlo sampling using Halton sequence
  • Antialiasing via 2x2 super-sampling
  • Using kd-tree for radiance estimation

以下、感想：

  • 500万フォトンくらいが限界か？
  • 画像は1024x768で500万フォトン放射、放射輝度推定に500フォトンしたものを512x384に縮小（はてなでさらに縮小されてる…）
  • 500万フォトンでも全然きれいにならない。メモリの制限に引っかかってしまうので、時間かければかけただけきれいにできないのが辛い。
  • フォトンが増えるとkd-treeの構築にすげー時間がかかるようになるのがヤバイ
  • 推定用フォトンの数を増やすと如実に時間がかかるようになるのがヤバイ。断じてO(k + log n)ではない。kd-treeがうまく動いてるのかどうか疑問
  • メモリ節約のため浮動小数型をfloatに変えると、精度不足のためか画像が乱れてしまったので断念
  • 一応ちゃんとsmallppmと同じような絵が出たので満足
    • 自分の勘違いかわからないんですが、smallppmでは最終的にピクセルの放射輝度を求めるところ（c[i]=c[i]+hp->flux*(1.0/(PIhp->r2num_photon*1000.0))）で拡散面のBRDF 1/π がかかってない？ように思うんですがどうでしょう？そのためこちらではシーン全体が暗くなってしまったので光源のパワーをπ倍して帳尻を合わせてみました。
  • OpenMPの#pragmaはsmallppmにあわせて入れてみたけど、動いてるのかどうか未確認
  • smallptからの現象で、ベクトルの破壊的な正規化を行うVec::norm()とそのベクトルを使う箇所があって
r = r + radiance(Ray(cam.o+d*140,d.norm()),0,Xi)*(1./samps);
  • dに対してnorm()呼び出しとスカラー乗算が同時に使われているけど、たぶんC++は関数への引数の評価順序は未定義なので（ググッても確証は得られなかった）、処理系依存になってるんではないかと思う（VC, gccとも正規化が先に実行されてるぽい）。なのでVec::norm()は破壊的操作を行わずに正規化したベクトルを返すように変えてみた。
  • std::vectorだと連続するメモリを要求してしまい厳しいため、ポインタにしている
  • LSDSEが難しい理由というのがよくわかってない。フォトンマッピングの場合光源から放射したフォトンはスペキュラ面以外にぶつかったら格納されてしまうので、L(SD){2,}Eは難しいのかなと思うんだけど。まあ拡散面に２回以上あたったものは明るさ的に重要度は低いだろうから必要はないだろうけど。

以下ソース：smallpm.cpp

#include <math.h>   // smallpm, Photon Mapping
#include <stdlib.h> // originally smallpt, a path tracer by Kevin Beason, 2008
#include <stdio.h>  // derived    smallppm, by T. Hachisuka
#include <vector>   // Usage: ./smallpm 1000000 100 && xv image.ppm
#include <algorithm>                 // #emit photons, #estimate photons
const double PI=3.14159265358979; int estimate=10;
int primes[61]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,
83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,
191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283};
inline int rev(const int i,const int p) {if (i==0) return i; else return p-i;}
double hal(const int b, int j) { // Halton sequence with reverse permutation
  const int p = primes[b]; double h = 0.0, f = 1.0 / (double)p, fct = f;
  while (j > 0) {h += rev(j % p, p) * fct; j /= p; fct *= f;} return h;}
struct Vec {double x, y, z; // vector: position, also color (r,g,b)
  Vec(double x_ = 0, double y_ = 0, double z_ = 0) {x = x_; y = y_; z = z_;}
  inline Vec operator+(const Vec &b) const {return Vec(x+b.x, y+b.y, z+b.z);}
  inline Vec operator-(const Vec &b) const {return Vec(x-b.x, y-b.y, z-b.z);}
  inline Vec operator*(double b) const {return Vec(x * b, y * b, z * b);}
  inline Vec mul(const Vec &b) const {return Vec(x * b.x, y * b.y , z * b.z);}
  inline double sqlen() const {return dot(*this);}
  inline Vec norm() {return (*this) * (1.0 / sqrt(sqlen()));}
  inline double dot(const Vec &b) const {return x * b.x + y * b.y + z * b.z;}
  Vec operator%(Vec&b) {return Vec(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);}};
struct Ray {Vec o, d; Ray(){}; Ray(Vec o_, Vec d_) : o(o_), d(d_) {}};
enum Refl_t {DIFF, SPEC, REFR};  // material types
struct Sphere {double rad; Vec p, c; Refl_t refl;
  Sphere(double r_,Vec p_,Vec c_,Refl_t re_) : rad(r_),p(p_),c(c_),refl(re_){}
  inline double intersect(const Ray &r) const { // returns distance
    Vec op=p-r.o; double t, b=op.dot(r.d), det=b*b-op.dot(op)+rad*rad;
    if (det < 0) return 1e20; else det = sqrt(det);
    return (t=b-det) > 1e-4 ? t : ((t=b+det)>1e-4 ? t : 1e20);}};
const Vec LightPos(50,60,85), LightPow(PI*10000,PI*10000,PI*10000);
const Sphere sph[] = { // Scene: radius, position, color, material
  Sphere(1e5, Vec( 1e5+1,40.8,81.6), Vec(.75,.25,.25),DIFF),//Left
  Sphere(1e5, Vec(-1e5+99,40.8,81.6),Vec(.25,.25,.75),SPEC),//Rght
  Sphere(1e5, Vec(50,40.8, 1e5),     Vec(.75,.75,.75),DIFF),//Back
  Sphere(1e5, Vec(50,40.8,-1e5+170), Vec(),           DIFF),//Frnt
  Sphere(1e5, Vec(50, 1e5, 81.6),    Vec(.75,.75,.75),DIFF),//Botm
  Sphere(1e5, Vec(50,-1e5+81.6,81.6),Vec(.75,.75,.75),DIFF),//Top
  Sphere(16.5,Vec(27,16.5,47),       Vec(1,1,1)*.999, SPEC),//Mirr
  Sphere(16.5,Vec(73,16.5,88),       Vec(1,1,1)*.999, REFR),//Glas
  Sphere(8.5, Vec(50,8.5,60),        Vec(1,1,1)*.999, DIFF)};//Mid
int toInt(double x){return int(pow(1-exp(-x),1/2.2)*255+.5);} //tone mapping
inline bool intersect(const Ray &r,double &t,int &id){//ray-sphere intrsect.
  int n = sizeof(sph) / sizeof(Sphere); double d, inf = 1e20; t = inf;
  for(int i=0;i<n;++i){d=sph[i].intersect(r);if(d<t){t=d;id=i;}}return t<inf;}
void genp(Ray* pr, Vec* f, int i) { *f = LightPow;
  double p=2*PI*hal(0,i), t=2*acos(sqrt(1.-hal(1,i))), st=sin(t);
  pr->d=Vec(cos(p)*st, cos(t), sin(p)*st); pr->o=LightPos; }
struct Photon { Vec x, fl; };
struct AABB { Vec l, u; AABB():l(1e20,1e20,1e20),u(-1e20,-1e20,-1e20){}
  inline void fit(const Vec &p) {
    l.x=std::min(p.x,l.x); l.y=std::min(p.y,l.y); l.z=std::min(p.z,l.z);
    u.x=std::max(p.x,u.x); u.y=std::max(p.y,u.y); u.z=std::max(p.z,u.z);}};
inline int sepaxis(Photon** photons, unsigned n) { AABB aabb;
  for (unsigned i=0; i<n; ++i) aabb.fit(photons[i]->x); Vec w=aabb.u-aabb.l;
  return w.x>w.y?(w.x>w.z?0:(w.y>w.z?1:2)):(w.y>w.z?1:(w.x>w.z?0:2));}
inline int log_2(int x) {int n; for (n=0; x>1; x>>=1,++n); return n;}
inline int median(int n) {int h=log_2(n),s=1<<h,d=n-s,s2=s/2;
  if (s2>0 && d>=s2) d=s2-1; return s2+d;}
bool cmpx(const Photon* a,const Photon* b) {return a->x.x < b->x.x;}
bool cmpy(const Photon* a,const Photon* b) {return a->x.y < b->x.y;}
bool cmpz(const Photon* a,const Photon* b) {return a->x.z < b->x.z;}
struct KDTree { Photon** photons; char* nodes; unsigned total;
  void build(Photon** buf,unsigned total_) {total=total_;
    photons=new Photon*[total]; nodes=new char[total/2]; sep(buf,0,total);}
  void nearest(Photon** buf,double* dist,int n,const Vec& x,double d2) {
    for (int i=0; i<n; ++i) {buf[i]=NULL; dist[i]=d2;} trav(buf,dist,n,x,0);}
  void sep(Photon** buf, unsigned i, unsigned n) {
    if (n==1) photons[i]=*buf; else if (n>1) { int axis=sepaxis(buf, n);
      std::sort(buf, buf+n, (axis==0)?cmpx:(axis==1)?cmpy:cmpz);
      int m=median(n); photons[i]=buf[m]; nodes[i]=axis;
      sep(buf,i*2+1,m); sep(buf+m+1,i*2+2,n-m-1); }}
  void trav(Photon** buf,double* dist,int n,const Vec& x,unsigned i) {
    if (i>=total) return; else if (i*2+1<total) {
      int axis=nodes[i]; Vec& v=photons[i]->x; double* pd2=&dist[n-1];
      double e=(axis==0 ? x.x-v.x : (axis==1 ? x.y-v.y : x.z-v.z));
      if(e<0){trav(buf,dist,n,x,i*2+1);if(e*e<*pd2)trav(buf,dist,n,x,i*2+2);}
      else   {trav(buf,dist,n,x,i*2+2);if(e*e<*pd2)trav(buf,dist,n,x,i*2+1);}}
    double e2=(photons[i]->x-x).sqlen(); if (e2>=dist[n-1]) return;
    int l,r,m; for (l=0,r=n;l<r;) {m=(l+r)/2; (e2<dist[m])?r=m:l=m+1;}
    for (int j=n; --j>l; ) {buf[j]=buf[j-1]; dist[j]=dist[j-1];}
    buf[l]=photons[i]; dist[l]=e2; } };
std::vector<Photon*> g_photons; KDTree g_kdtree;
Vec trace(const Ray &r, int dpt, bool m, 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 (m) { Photon** b=(Photon**)alloca(sizeof(Photon*)*estimate);
      double* rr=(double*)alloca(sizeof(double)*estimate);
      g_kdtree.nearest(b,rr,estimate,x,1e20); Vec flux;const double fr=1.0/PI;
      int j,k; for (j=k=0; j<estimate; k=j++) { Photon* q=b[j]; if (!q) break;
        flux=flux+q->fl*fr; } return flux.mul(f)*(1.0/(PI*rr[k])); }
    Photon* q=new Photon; q->x=x; q->fl=fl; g_photons.push_back(q);
    double p=std::max(f.x,std::max(f.y,f.z)); if (hal(d3+1,i)>=p)return Vec();
    double r1=2.*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,m,f.mul(fl)*(1./p),i);
  } else if (obj.refl == SPEC) {
    return trace(Ray(x,r.d-n#n.dot(r.d)),dpt,m,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,m,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 m ? (trace(lr,dpt,m,fl,i)*Re + trace(rr,dpt,m,fl,i)*(1-Re)).mul(f)
             : hal(d3-1,i)<P ?trace(lr,dpt,m,rfl,i) :trace(rr,dpt,m,rfl,i);}}
int main(int argc, char *argv[]) { int w=1024, h=768; int samps=1, BKT=1000;
  argc>1&&(samps=atoi(argv[1])/BKT); argc>2&&(estimate=atoi(argv[2]));
  for(int i=0; i<samps; ++i) { int m=BKT*i; Ray r; Vec f;
    fprintf(stderr,"\rPhotonPass %5.2f%%",100.*(i+1)/samps);
    for(int j=0;j<BKT;++j) { genp(&r,&f,m+j); trace(r,0,false,f,m+j+1); }}
  fprintf(stderr, "\n#photons:%d\nBuilding kd-tree...\n", g_photons.size());
  if (!g_photons.empty()) g_kdtree.build(&g_photons[0], g_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) {
    fprintf(stderr, "\rHitPointPass %5.2f%%", 100.0*y/(h-1));
    for (int x=0; x<w; ++x) { Vec r, dmy; const double s=1.0/(BKT*samps#2);
      for (int v=0; v<2; ++v) for (int u=0; u<2; ++u) {
        Vec d = cx * ((x+u*.5+.25)/w-.5) + cy * (-(y+v*.5+.25)/h+.5) + cam.d;
        r=r+trace(Ray(cam.o+d*130,d.norm()),0,true,dmy,0);} cc[x+y*w]=r*s;}}
  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));}}

関連記事

新着記事