thodg/smallpt

diff --git a/smallpt.cpp b/smallpt.cpp
index 61e2eda..da985ab 100644
--- a/smallpt.cpp
+++ b/smallpt.cpp
@@ -1,33 +1,111 @@
-#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);}
+/* smallpt, a Path Tracer by Kevin Beason, 2008
+ * Make : g++ -O3 -fopenmp smallpt.cpp -o smallpt
+ *        Remove "-fopenmp" for g++ version < 4.2
+ * Usage: time ./smallpt 5000 && xv image.ppm
+ */
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+typedef double F;
+inline F Fsqrt (F x) { return sqrt(x); }
+inline F Fabs (F x) { return fabs(x); }
+
+struct Vec {
+  F x; // position, also color (r,g,b)
+  F y;
+  F z;
+
+  Vec (F x_ = 0, F y_ = 0, F 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 * (F 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 & normalize () {
+    return *this = *this * (1 / Fsqrt(x * x + y * y + z * z));
+  }
+
+  F dot (const Vec &b) const
+  {
+    return x * b.x + y * b.y + z * b.z;
+  }
+
+  // cross product
+  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 origin;
+  Vec direction;
+
+  Ray (Vec origin_, Vec direction_)
+    : origin(origin_), direction(direction_)
+  {}
 };
-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()
+
+// material types, used in radiance()
+enum Refl_t { DIFF, SPEC, REFR };
+
 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);
+  F radius;            // radius
+  Vec position;
+  Vec emission;
+  Vec color;
+  Refl_t reflection_type;
+  
+  Sphere (F radius_, Vec position_, Vec emission_, Vec color_,
+          Refl_t reflection_type_)
+    : radius(radius_), p(p_), e(e_), c(c_), refl(refl_)
+  {}
+
+  // returns distance, 0 if nohit
+  F intersect (const Ray &r) const
+  {
+    /* Solve :
+     * t^2 * d . d + 2 * t * (o - p) . d + (o - p) . (o - p) - R^2 = 0
+     */
+    Vec op = p - r.origin;
+    F t;
+    F eps = 1e-4;
+    F b = op.dot(r.direction);
+    F det = b * b - op.dot(op) + radius * radius;
+    if (det < 0)
+      return 0;
+    else
+      det = Fsqrt(det);
+    return (t = b-det) > eps ? t : ((t = b + det) > eps ? t : 0);
   }
 };
-Sphere spheres[] = {//Scene: radius, position, emission, color, material
+
+//Scene: radius, position, emission, color, material
+Sphere g_spheres[] = {
   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
@@ -38,56 +116,116 @@ Sphere spheres[] = {//Scene: radius, position, emission, color, material
   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;
+
+inline F clamp (F x)
+{
+  return x < 0 ? 0 : x > 1 ? 1 : x;
+}
+
+inline int toInt (F x)
+{
+  return int(pow(clamp(x), 1 / 2.2) * 255 + .5);
 }
-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);
+
+inline bool intersect (const Ray &r, F &t, int &id)
+{
+  int i = sizeof(g_spheres) / sizeof(Sphere);
+  F d;
+  F inf = 1e20;
+  t = inf;
+  while (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)
+{
+  F t;                                  // Distance to intersection
+  int id = 0;                           // Id of intersected object
+  if (! intersect(r, t, id))            // If miss, return black
+    return Vec();
+  const Sphere &obj = spheres[id];      // The hit object
+  Vec x = r.origin + r.direction * t;
+  Vec n = (x - obj.position).normalize();
+  Vec nl = n.dot(r.direction) < 0 ? n : n * -1;
+  Vec f = obj.c;
+  F 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.emission; //R.R.
+  if (obj.reflection_type == DIFF) {    // Ideal diffuse reflection
+    F r1 = 2 * M_PI * erand48(Xi);
+    F r2 = erand48(Xi);
+    F r2s = Fsqrt(r2);
+    Vec w = nl;
+    Vec u = ((Fabs(w.x) > .1 ? Vec(0,1) : Vec(1)) % w).normalize();
+    Vec v= w % u;
+    Vec d = (u * cos(r1) * r2s +
+             v * sin(r1) * r2s +
+             w * Fsqrt(1 - r2)).normalize();
+    return obj.emission + f.mult(radiance(Ray(x, d), depth, Xi));
+  }
+  else if (obj.refl == SPEC)            // Ideal specular reflection
+    return obj.emission +
+      f.mult(radiance(Ray(x, r.direction - n * 2 * n.dot(r.direction)),
+                      depth, Xi));
+  else {                                // Ideal dielectric refraction
+  Ray reflRay(x, r.direction - n * 2 * n.dot(r.direction));
+  bool into = n.dot(nl) > 0;            // Ray from outside going in?
+  F nc = 1;
+  F nt = 1.5;
+  F nnt = into ? nc / nt : nt / nc;
+  F ddn = r.direction.dot(nl);
+  F cos2t;
+  cos2t = 1 - nnt * nnt * (1 - ddn * ddn);
+  if (cos2t < 0)                        // Total internal reflection
+    return obj.emission + f.mult(radiance(reflRay, depth, Xi));
+  Vec tdir = (r.direction * nnt - n *
+              ((into ? 1 : -1) *
+               (ddn * nnt + Fsqrt(cos2t)))).normalize();
+  F a = nt - nc;
+  F b = nt + nc;
+  F R0 = a * a / (b * b);
+  F c = 1 - (into ? -ddn : tdir.dot(n));
+  F Re = R0 + (1 - R0) * c * c * c * c * c;
+  F Tr = 1 - Re;
+  F P = .25 + .5 * Re;
+  F RP = Re / P;
+  F TP = Tr / (1 - P);
+  return obj.emission +
+    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];
+
+int main (int argc, char *argv[])
+{
+  int w = 1024;
+  int h = 768;
+  int samples = argc == 2 ? atoi(argv[1]) / 4 : 1;
+  Ray camera(Vec(50,52,295.6), Vec(0,-0.042612,-1).normalize());
+  Vec cx=Vec(w*.5135/h), cy=(cx%camera.direction).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
+    fprintf(stderr,"\rRendering (%d spp) %5.2f%%",samples*4,100.*y/(h-1));
+    for (unsigned short x=0, Xi[3]={0,0,static_cast<unsigned short>(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);
+          for (int s=0; s<samples; s++){
+            F r1=2*erand48(Xi), dx=r1<1 ? Fsqrt(r1)-1: 1-Fsqrt(2-r1);
+            F r2=2*erand48(Xi), dy=r2<1 ? Fsqrt(r2)-1: 1-Fsqrt(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);
+              cy*( ( (sy+.5 + dy)/2 + y)/h - .5) + camera.direction;
+            r = r + radiance(Ray(camera.origin + d*140,d.norm()),0,Xi)*(1./samples);
           } // Camera rays are pushed ^^^^^ forward to start in interior
           c[i] = c[i] + Vec(clamp(r.x),clamp(r.y),clamp(r.z))*.25;
         }