imHo RSSフィード

2009-10-15

smallpm: Global Illumination in 128 lines of C++

f:id:mokehehe:20091016104608j:image:right

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/(PI*hp->r2*num_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<Photon>だと連続するメモリを要求してしまい厳しいため、ポインタにしている
  • LSDSEが難しい理由というのがよくわかってない。フォトンマッピングの場合光源から放射したフォトンはスペキュラ面以外にぶつかったら格納されてしまうので、L(SD){2,}Eは難しいのかなと思うんだけど。まあ拡散面に2回以上あたったものは明るさ的に重要度は低いだろうから必要はないだろうけど。

以下ソース:smallpm.cpp

#include <math.h>   // smallpm, Photon Mapping by mokehehe
#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*2*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*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));}}

BeeBee 2009/10/21 02:20 smallpmいいですね!個人的にはphoton mappingの方が
若干詰め込むのが難しいと思っていたので,凄いと思います.

BRDFのPIの件はご指摘の通りですね.直しておこうと思います.
d.norm()は気がつきませんでした.Intelのコンパイラで挙動が
若干おかしい場合があったので,もしかしたらこれのせいだったのかもしれません.

LSDSEが何故難しいかは,点光源かつピンホールカメラを使った場合に,
どのようにLSDSEというパスを構築するかを考えるとわかりやすいと思います.
例えば,カメラ(ピンホールの位置)からレイを飛ばしてLSDSEを作ろうとすると,
Lの前にSがあるので光源にレイがヒットする必要がありますが,点光源なので不可能です.
光源からレイを飛ばしてもピンホールを通る確率はゼロなので,これも不可能です.

それなら両方から飛ばせばいいのではと考えたとしても,結局は光源と
視点を結ぶDの位置を決めなければならないので,単純な形状以外では
非常に困難になります.smallppmのシーンは簡単に見えますが,unbiasedを謳っている
手法では正しくレンダリングできないと思います.パストレで同じシーンを
レンダリングすることを考えると難しさが実感できるかもしれません.

# 現実には光源とカメラは点ではありませんが,シーンのスケールに対して
非常に小さいことが多いので,ヒットする確率はかなり小さくなります.

mokehehemokehehe 2009/10/21 16:51 ほとんどsmallppmを流用させていただいてますので、なんとか作れた感じです (^^;;

ああ、Photon Mappingでは簡単にLSDSEができるけどbiasedなんですね。
それに対してProgressive Photon Mappingはunbiasedだということですね。
ようやくわかりました、ありがとうございます。

BeeBee 2009/10/22 02:17 > それに対してProgressive Photon Mappingはunbiasedだということですね。

この点,非常に微妙ですがPPMはbiasedです.unbiasedとbiasedの大まかな
違いは,アルゴリズムに適当なサンプル数を入れて何枚も画像を出力して,
平均した結果が沢山サンプル数を取った場合と同じになるかどうか,
というのがわかりやすいと思います.例えばunbiasedなパストレでは
そうなりますが,biasedなフォトンマッピング/PPMでは,
有限のフォトンを使った画像を平均してもぼけた結果が出るだけです.

この「平均」の操作がちょっと特殊なのがPPMで,結果としてbiasを
徐々に減らすことができるという事です.でも定義としてbiasedで
あることは変わりません.まあ,biasedとunbiasedという言葉自体
実用上あまり意味無いと思うので,
「おまえただ(un)biasedって言いたかっただけじゃねえの」
程度に捉えるのがいいかもと思っています.

mokehehemokehehe 2009/10/22 09:45 なるほど、biasって「期待値がある傾向に偏る=全体的に色が明るく/暗くなる」ということかと思ってました。
そうなるとunbiasedはモンテカルロレイトレタイプじゃないと達成できないような。

はてなユーザーのみコメントできます。はてなへログインもしくは新規登録をおこなってください。