(L) [2009/07/12] [toxie] [smallppm] Wayback!thanx for this nice piece of code.. i took the freedom to fix compile errors (vc and ic) and fix the sampling (especially with openmp turned on [SMILEY :)])..
(hoping that i didn't introduce some garbage with these changes.. [SMILEY ;)])
so the tiny updated version is here:
edit: seems my openmp 'sampling fix' introduced some qmc related patterns on the walls.. [SMILEY :(] somebody has an idea?!
edit2: hmmm.. somehow i had to trick openmp by using a threadprivate()/copyin() instead of private()?! then it worked.. very strange.. updated source below
edit3: using hammersley seemed a bit suspicous for PPM (as it was used in the code), so i replaced the first sampling component also with halton.. plus saved some more code lines and added some tiny optimzations.. and added simple tonemapping in addition to the already existant gamma correction..
edit4: some additional code squeeze and tiny optimizations..
edit5: 116 lines! plus a bit of speedup
#include <math.h> // smallppm, Progressive Photon Mapping by T. Hachisuka
#include <stdlib.h> // originally smallpt, a Path Tracer by Kevin Beason, 2008
#include <stdio.h> // Usage: time ./smallppm 100000 && xv image.ppm
const double PI=3.14159265358979,inf=1e20,eps=1e-4; // ^^^^^^:emitted photons
const double ALPHA = .7; // the alpha parameter of PPM
const int prim[59] = {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};
inline double hal(unsigned int j,const unsigned int pb) { // Halton sequence
double h=0., f=1./prim[pb], fct=f;
while(j>0) {h+=(j%prim[pb])*fct; j/=prim[pb]; fct*=f;} return h; }
struct Vec { double x, y, z; Vec(double v=0) : x(v),y(v),z(v) {}; // pos/r,g,b
inline Vec(double x_, double y_, double z_) : 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*(const 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 Vec& norm() {return (*this) = (*this) * (1./sqrt(x*x+y*y+z*z));}
Vec operator%(const Vec&b){return Vec(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);}
inline double dot(const Vec &b) const {return x*b.x + y*b.y + z*b.z;} };
struct AABB { Vec min, max; // axis aligned bounding box
inline void fit(const Vec &p) {
if(p.x<min.x)min.x=p.x; if(p.y<min.y)min.y=p.y; if(p.z<min.z)min.z=p.z;
if(p.x>max.x)max.x=p.x; if(p.y>max.y)max.y=p.y; if(p.z>max.z)max.z=p.z; }
inline void reset() {min=Vec(inf); max=Vec(-inf);} };
struct HPoint {Vec f,pos,nrm,flux; double r2; unsigned int n,pix;};
struct List {HPoint *id; List *next;}; unsigned int num_hash,num_photon;
List* ListAdd(HPoint *i,List* h){List* p=new List;p->id=i;p->next=h;return p;}
double hashs; List **hashg; List *hitps = NULL; AABB hpbbox;
inline unsigned int hash(const unsigned int x, const unsigned int y, const
unsigned int z) { return (x*73856093 ^ y*19349663 ^ z*83492791)%num_hash; }
void build_hash_grid(const double wh) { hpbbox.reset(); const List *lst=hitps;
while(lst != NULL) {HPoint *hp=lst->id; lst=lst->next; hpbbox.fit(hp->pos);}
const Vec ssize = hpbbox.max - hpbbox.min; // compute initial radius
const double irad = (ssize.x + ssize.y + ssize.z) * (4.0/3.0) / wh;
hpbbox.reset(); lst = hitps; num_hash = 0; // determine hash size
while(lst != NULL) { HPoint *hp=lst->id; lst=lst->next; hp->r2=irad*irad;
hp->n=0; hp->flux=Vec(); ++num_hash; hpbbox.fit(hp->pos-irad);
hpbbox.fit(hp->pos+irad); } hashs=0.5/irad; hashg=new List*[num_hash];
for (unsigned int i=0; i<num_hash;++i) hashg[i] = NULL;
lst = hitps; while(lst != NULL) { // store hitpoints in the hashed grid
HPoint *hp=lst->id; const Vec BMin=((hp->pos-irad)-hpbbox.min)*hashs,
BMax=((hp->pos+irad)-hpbbox.min)*hashs; lst=lst->next;
for (int iz = abs(int(BMin.z)); iz <= abs(int(BMax.z)); ++iz)
for (int iy = abs(int(BMin.y)); iy <= abs(int(BMax.y)); ++iy)
for (int ix = abs(int(BMin.x)); ix <= abs(int(BMax.x)); ++ix)
{const int hv=hash(ix,iy,iz); hashg[hv]=ListAdd(hp,hashg[hv]);}} }
struct Ray {Vec o, d; Ray(const Vec &o_, const Vec &d_) : o(o_), d(d_) {} };
enum Refl_t {DIFF, SPEC, REFR}; // material types, used in trace()
struct Sphere { const double rads; Vec p, c; Refl_t refl;
Sphere(double r_,Vec p_,Vec c_,Refl_t e_):rads(r_*r_),p(p_),c(c_),refl(e_){}
inline double intersect(const Ray &r) const { // returns distance
const Vec op=p-r.o; double t, b=op.dot(r.d), det=b*b-op.dot(op)+rads;
if(det<0) return inf; det=sqrt(det);
return (t=b-det)>eps ? t : ((t=b+det)>eps ? t : inf); } };
const Sphere sph[9] = { // 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) ,DIFF),//Back
Sphere(1e5, Vec(50,40.8,-1e5+170), Vec(), DIFF),//Frnt
Sphere(1e5, Vec(50, 1e5, 81.6), Vec(.75) ,DIFF),//Botm
Sphere(1e5, Vec(50,-1e5+81.6,81.6),Vec(.75) ,DIFF),//Top
Sphere(16.5,Vec(27,16.5,47), Vec(.999) ,SPEC),//Mirr
Sphere(16.5,Vec(73,16.5,88), Vec(.999) ,REFR),//Glas
Sphere(8.5, Vec(50,8.5,60), Vec(.999) ,DIFF) };//Mid
int toInt(double x){return int(pow((x < 0) ? 0:(x > 1) ? 1:x,1/2.2)*255.+.5);}
inline bool intersect(const Ray &r, double &t, int &id) {// ray-sphere inters.
const int n=sizeof(sph)/sizeof(Sphere); t=inf; for(int i=0;i<n;++i)
{const double d=sph[i].intersect(r); if(d<t){t=d;id=i;} } return t<inf; }
void genp(Ray &pr,unsigned int j) { double p=(2.*PI)*hal(j,0),t=2.*acos(sqrt(
hal(j,1))),st=sin(t);pr=Ray(Vec(50,60,85),Vec(cos(p)*st,cos(t),sin(p)*st));}
void trace(const Ray &r,int dpt,bool m,const Vec &flux,const Vec &adj,int j) {
double t; int id; if (!intersect(r, t, id) || (dpt++ >= 19)) return;
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; double p=f.x>f.y&&f.x>f.z?f.x:f.y>f.z?f.y:f.z;
if (obj.refl == DIFF) { if (m) { HPoint* hp=new HPoint; hp->f=f.mul(adj);
hp->pos=x; hp->nrm=n; hp->pix=j; hitps = ListAdd(hp, hitps);
} else { const Vec u=((fabs(nl.x)>.1?Vec(0,1,0):Vec(1,0,0))%nl).norm(),
v=nl%u, hh=(x-hpbbox.min)*hashs;
#pragma omp critical // critical section in OpenMP
{ List* hp=hashg[hash(abs(int(hh.x)),abs(int(hh.y)),abs(int(hh.z)))];
while(hp!=NULL) { HPoint *hitp=hp->id;hp=hp->next; Vec v=hitp->pos-x;
if ((hitp->nrm.dot(n) > eps) && (v.dot(v) <= hitp->r2)) {
double g = (hitp->n*ALPHA+ALPHA) / (hitp->n*ALPHA+1.);
hitp->r2*=g; hitp->n++;
hitp->flux=(hitp->flux + hitp->f.mul(flux))*g; } } }
if(hal(j,dpt*3-1)<p) { const double r1=(2.*PI)*hal(j,dpt*3),r2=hal(
j,dpt*3+1),r2s=sqrt(r2); trace(Ray(x,(u*(cos(r1)*r2s)+v*(sin(r1)*r2s)+
nl*sqrt(1-r2)).norm()),dpt,0,f.mul(flux)*(1./p),adj,j); } }
} else { const Ray lr(x,r.d-n*(2.*n.dot(r.d))); if (obj.refl == SPEC)
trace(lr, dpt, m, f.mul(flux), f.mul(adj), j);
else { const bool into = (n.dot(nl)>0.); const double nc=1., nt=1.5,
nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t=1-nnt*nnt*(1.-ddn*ddn);
if (cos2t<0.) return trace(lr,dpt,m,flux,adj,j);
Vec td = (r.d*nnt - n*((into?1.:-1.)*(ddn*nnt+sqrt(cos2t)))).norm();
const double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c=1.-(into?-ddn:td.dot(n));
const double Re=R0+(1-R0)*(c*c)*(c*c)*c; Ray rr(x,td);Vec fa=f.mul(adj);
if(m) {trace(lr,dpt,m,flux,fa*Re,j); trace(rr,dpt,m,flux,fa*(1.-Re),j);}
else trace((hal(j,dpt*3-1)<Re)?lr:rr,dpt,m,flux,fa,j); } } }
int main(int argc, char *argv[]) {
int w=1024, h=768; num_photon = (argc==2) ? atoi(argv[1])/1000+1 : 1;
const Vec camd(0,-.04257,-.99909); unsigned int pj=0;
Vec cx=Vec(w*.5135/h,0,0), cy=(cx%camd).norm()*.5135, *c=new Vec[w*h];
for(int y=0;y<h;++y) {fprintf(stderr,"\rHitPointPass %3.2f%%",100.*y/(h-1));
for(int x=0;x<w;++x,++pj) { Vec d=cx*((x+.5)/w-.5)-cy*((y+.5)/h-.5)+camd;
trace(Ray(Vec(50,52,295.6)+d*140.,d.norm()),0,1,Vec(),Vec(1),pj); } }
fprintf(stderr,"\n"); build_hash_grid(w+h);
#pragma omp parallel for schedule(dynamic, 1)
for(int i=0;i<num_photon;++i) { unsigned int pj=i*1000; Vec b,f;Ray r(b,b);
fprintf(stderr,"\rPhotonPass %3.2f%%",100.*(i+1)/num_photon);
for(int j=0;j<1000;++j,++pj){genp(r,pj);trace(r,0,0,Vec(2e4),Vec(1),pj);}}
List* lst=hitps; while(lst != NULL) { HPoint* hp=lst->id;lst=lst->next;
int i=hp->pix; c[i]=c[i]+hp->flux*((1./(PI*1000.))/(hp->r2*num_photon)); }
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){double a=.2125*c[i].x+.7154*c[i].y+.0721*c[i].z + 1.;
fprintf(f,"%d %d %d ",toInt(c[i].x/a),toInt(c[i].y/a),toInt(c[i].z/a));} }