Smallpt: Global Illumination in 99 lines of C++
smallpt: Global Illumination in 99 lines of C++

smallpt is a global illumination renderer.
It is 99 lines of C++, is open source,
and renders the above scene using unbiased Monte Carlo path
tracing (click for full size).


 #include <math.h>   // smallpt, a Path Tracer by Kevin Beason, 2008 
 #include <stdlib.h> // Make : g++ -O3 -fopenmp smallpt.cpp -o smallpt 
 #include <stdio.h>  //        Remove "-fopenmp" for g++ version < 4.2 
 struct Vec {        // Usage: time ./smallpt 5000 && xv image.ppm 
   double x, y, z;                  // position, also color (r,g,b) 
   Vec(double x_=0, double y_=0, double z_=0){ x=x_; y=y_; z=z_; } 
   Vec operator+(const Vec &b) const { return Vec(x+b.x,y+b.y,z+b.z); } 
   Vec operator-(const Vec &b) const { return Vec(x-b.x,y-b.y,z-b.z); } 
   Vec operator*(double b) const { return Vec(x*b,y*b,z*b); } 
   Vec mult(const Vec &b) const { return Vec(x*b.x,y*b.y,z*b.z); } 
   Vec& norm(){ return *this = *this * (1/sqrt(x*x+y*y+z*z)); } 
   double dot(const Vec &b) const { return x*b.x+y*b.y+z*b.z; } // cross: 
   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(Vec o_, Vec d_) : o(o_), d(d_) {} }; 
 enum Refl_t { DIFF, SPEC, REFR };  // material types, used in radiance() 
 struct Sphere { 
   double rad;       // radius 
   Vec p, e, c;      // position, emission, color 
   Refl_t refl;      // reflection type (DIFFuse, SPECular, REFRactive) 
   Sphere(double rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_): 
     rad(rad_), p(p_), e(e_), c(c_), refl(refl_) {} 
   double intersect(const Ray &r) const { // returns distance, 0 if nohit 
     Vec op = p-r.o; // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0 
     double t, eps=1e-4, b=op.dot(r.d), det=b*b-op.dot(op)+rad*rad; 
     if (det<0) return 0; else det=sqrt(det); 
     return (t=b-det)>eps ? t : ((t=b+det)>eps ? t : 0); 
   } 
 }; 
 Sphere spheres[] = {//Scene: radius, position, emission, color, material 
   Sphere(1e5, Vec( 1e5+1,40.8,81.6), Vec(),Vec(.75,.25,.25),DIFF),//Left 
   Sphere(1e5, Vec(-1e5+99,40.8,81.6),Vec(),Vec(.25,.25,.75),DIFF),//Rght 
   Sphere(1e5, Vec(50,40.8, 1e5),     Vec(),Vec(.75,.75,.75),DIFF),//Back 
   Sphere(1e5, Vec(50,40.8,-1e5+170), Vec(),Vec(),           DIFF),//Frnt 
   Sphere(1e5, Vec(50, 1e5, 81.6),    Vec(),Vec(.75,.75,.75),DIFF),//Botm 
   Sphere(1e5, Vec(50,-1e5+81.6,81.6),Vec(),Vec(.75,.75,.75),DIFF),//Top 
   Sphere(16.5,Vec(27,16.5,47),       Vec(),Vec(1,1,1)*.999, SPEC),//Mirr 
   Sphere(16.5,Vec(73,16.5,78),       Vec(),Vec(1,1,1)*.999, REFR),//Glas 
   Sphere(600, Vec(50,681.6-.27,81.6),Vec(12,12,12),  Vec(), DIFF) //Lite 
 }; 
 inline double clamp(double x){ return x<0 ? 0 : x>1 ? 1 : x; } 
 inline int toInt(double x){ return int(pow(clamp(x),1/2.2)*255+.5); } 
 inline bool intersect(const Ray &r, double &t, int &id){ 
   double n=sizeof(spheres)/sizeof(Sphere), d, inf=t=1e20; 
   for(int i=int(n);i--;) if((d=spheres[i].intersect(r))&&d<t){t=d;id=i;} 
   return t<inf; 
 } 
 Vec radiance(const Ray &r, int depth, unsigned short *Xi){ 
   double t;                               // distance to intersection 
   int id=0;                               // id of intersected object 
   if (!intersect(r, t, id)) return Vec(); // if miss, return black 
   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 
     double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), 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)).norm(); 
     return obj.e + f.mult(radiance(Ray(x,d),depth,Xi)); 
   } else if (obj.refl == SPEC)            // Ideal SPECULAR reflection 
     return obj.e + f.mult(radiance(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi)); 
   Ray reflRay(x, r.d-n*2*n.dot(r.d));     // Ideal dielectric REFRACTION 
   bool into = n.dot(nl)>0;                // Ray from outside going in? 
   double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t; 
   if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0)    // Total internal reflection 
     return obj.e + f.mult(radiance(reflRay,depth,Xi)); 
   Vec tdir = (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:tdir.dot(n)); 
   double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P); 
   return obj.e + f.mult(depth>2 ? (erand48(Xi)<P ?   // Russian roulette 
     radiance(reflRay,depth,Xi)*RP:radiance(Ray(x,tdir),depth,Xi)*TP) : 
     radiance(reflRay,depth,Xi)*Re+radiance(Ray(x,tdir),depth,Xi)*Tr); 
 } 
 int main(int argc, char *argv[]){ 
   int w=1024, h=768, samps = argc==2 ? atoi(argv[1])/4 : 1; // # samples 
   Ray cam(Vec(50,52,295.6), Vec(0,-0.042612,-1).norm()); // cam pos, dir 
   Vec cx=Vec(w*.5135/h), cy=(cx%cam.d).norm()*.5135, r, *c=new Vec[w*h]; 
 #pragma omp parallel for schedule(dynamic, 1) private(r)       // OpenMP 
   for (int y=0; y<h; y++){                       // Loop over image rows 
     fprintf(stderr,"\rRendering (%d spp) %5.2f%%",samps*4,100.*y/(h-1)); 
     for (unsigned short x=0, Xi[3]={0,0,y*y*y}; x<w; x++)   // Loop cols 
       for (int sy=0, i=(h-y-1)*w+x; sy<2; sy++)     // 2x2 subpixel rows 
         for (int sx=0; sx<2; sx++, r=Vec()){        // 2x2 subpixel cols 
           for (int s=0; s<samps; s++){ 
             double r1=2*erand48(Xi), dx=r1<1 ? sqrt(r1)-1: 1-sqrt(2-r1); 
             double r2=2*erand48(Xi), dy=r2<1 ? sqrt(r2)-1: 1-sqrt(2-r2); 
             Vec d = cx*( ( (sx+.5 + dx)/2 + x)/w - .5) + 
                     cy*( ( (sy+.5 + dy)/2 + y)/h - .5) + cam.d; 
             r = r + radiance(Ray(cam.o+d*140,d.norm()),0,Xi)*(1./samps); 
           } // Camera rays are pushed ^^^^^ forward to start in interior 
           c[i] = c[i] + Vec(clamp(r.x),clamp(r.y),clamp(r.z))*.25; 
         } 
   } 
   FILE *f = fopen("image.ppm", "w");         // Write image to PPM file. 
   fprintf(f, "P3\n%d %d\n%d\n", w, h, 255); 
   for (int i=0; i<w*h; i++) 
     fprintf(f,"%d %d %d ", toInt(c[i].x), toInt(c[i].y), toInt(c[i].z)); 
 } 
 The argument to smallpt ("5000") is the number of samples per pixel, and must be
greater than 4. If you don't have ImageMagick's
display
command available, you can use other viewers such as gthumb, gwenview, xv, gimp, etc.

GCC 4.2 or newer is required for multi-threading,
however older versions of gcc will still work without threading.
Remove the "-fopenmp" flag to disable threading support.


smallpt compiles with GCC 4.2 down to a 16 KB
executable and produces a 1024x768 resolution image. Rendering using 5000 paths
per pixel takes 2.1 hours on an Intel Core 2 Quad machine.

  A slightly modified version is provided which compiles to a 4 KB (4049 bytes)
compressed executable. Do "make smallpt4k" to build. 
See smallpt4k.cpp and the Makefile in the full distribution for more information.
Note sstrip
and p7zip are required.
  
If self assembly is used, the binary is only 2.7 KB (2666 bytes).
Do "make smallptSA" to build.

The box scene is constructed out of nine very large spheres which overlap.
The camera is very close to their surface so they appear flat, yet each wall is
actually the side of a sphere. The light is another sphere poking through the ceiling
which is why it's round instead of the normal square.
The black area visible in the mirror ball
is the front of the box colored black to appear like empty space.

The image is computed by solving the
rendering equation  
using
numerical
  integration.
The specific algorithm is Monte Carlo path tracing with
Russian roulette
for path termination. 
I highly recommend Shirley's and Jensen's execellent books (linked below) for
explanations of these ideas.
Due to size constraints and simplicity,
explicit light sampling is not used, nor any ray intersection acceleration data structure.

Parallelism is achieved using an OpenMP #pragma to dynamically
allocate rows of the image to different threads, with a thread for each processor or core.
The variable Xi is used to store the state of the random number
generator erand48(), and is seeded using an arbitrary function of
the row number to decorrelate (at least visually) the sequences from row-to-row.
In this way the sequences are deterministic and consistent from run to run,
and independent of which thread is executing and in what order the rows are executed.

Anti-aliasing is done using supersampling, which removes all the jaggies
except around the light. These are handled by using 2x2 subpixels which
are clamped and then averaged.

Instead of fclose()ing the file, I exploit the
C++ standard which
calls return(0) implicitly, which in turn calls exit(0), which flushes and closes open files.

