thodg/smallpt

diff --git a/Makefile b/Makefile
index 10bd1c3..c632a9d 100644
--- a/Makefile
+++ b/Makefile
@@ -1,118 +1,12 @@
-# smallpt Makefile
-#
-# 
 
+PROG = smallpt
 
-# Select appropriate compiler.
-#CPP=g++
-CPP=g++-4.2
-CPPFLAGS=-O3 -fopenmp       # Use this for gcc >= 4.2
-#CPPFLAGS=-O3               # Use this for gcc <  4.2
+CXX_SOURCES = smallpt.cpp
 
-ifdef DEBUG
-CPPFLAGS=-g
-endif
-
-default: all
-
-all: smallpt smallpt4k
-
-smallpt: smallpt.cpp Makefile
-	$(CPP) $(CPPFLAGS) -o smallpt smallpt.cpp 
+${PROG}: ${CXX_SOURCES}
+	${CXX} ${CPPFLAGS} ${CXXFLAGS} ${LDFLAGS} ${CXX_SOURCES} -o ${PROG}
 
 clean:
-	-rm smallpt smallpt4k
-
-tkdiff:
-	for i in `svn st -q | cut -c 8-`; do tkdiff $$i; done
-
-# Requires sstrip utility, which is part of ELFkickers:
-# http://www.muppetlabs.com/~breadbox/software/elfkickers.html
-# Also requires, 7z, part of p7zip:
-# http://p7zip.sourceforge.net/
-# The compile options were found with GC masher:
-# http://ftp.kameli.net/pub/fit/misc/gcmasher11082005.tar.gz
-# Guides I followed:
-# http://in4k.untergrund.net/index.php?title=Linux
-# http://ftp.kameli.net/pub/fit/misc/presis_asm06.pdf
-# http://www.muppetlabs.com/~breadbox/software/tiny/teensy.html
-#
-# compiles to 4049 bytes on 32-bit kubuntu fiesty
-# compiles to 4054 bytes on 64-bit kubuntu hardy (to a 32 bit binary)
-smallpt4k: smallpt4k.cpp Makefile
-	$(CPP) -m32 -march=native -O1 -nostdlib -fpeephole2 -ftree-vectorize -fopenmp -ffast-math -fomit-frame-pointer -fno-exceptions -c smallpt4k.cpp
-	ld -dynamic-linker /lib/ld-linux.so.2 smallpt4k.o /usr/local/lib/libgomp.so /lib/libm.so.6 -o usmallpt4k            # 32-bit OS
-#	ld -dynamic-linker /lib32/ld-linux.so.2 smallpt4k.o /usr/lib32/libgomp.so /lib32/libm.so.6 -o usmallpt4k -melf_i386 # 64-bit OS
-	strip -R .note -R .comment -R .eh_frame -R .eh_frame_hdr -R .gnu.version usmallpt4k
-	./sstrip usmallpt4k         # 32-bit OS
-#	./sstrip32 usmallpt4k       # 64-bit OS
-	echo 'HOME=/tmp/S;cp $$0 ~;tail -n+2 $$0|zcat>~;~;rm ~;exit' > smallpt4k
-	7z a -tgzip -mx=9 -so dummy usmallpt4k >> smallpt4k
-	#gzip -n --best -c usmallpt4k >> smallpt4k
-	chmod +x smallpt4k
-	ls -l smallpt4k
-
-# This builds a self assembling (SA) binary.
-# File size is 2666 bytes (2.7KB) on 32-bit Ubuntu gutsy with g++ 4.2.
-smallptSA: smallpt.cpp Makefile
-	echo 'a=/tmp/I;tail -n+2 $$0|zcat>$$a.cpp;g++-4.2 -O3 -fopenmp -o $$a $$a.cpp;$$a $$*;rm $$a $$a.cpp;exit' > smallptSA
-	7z a -tgzip -mx=9 -so dummy smallpt.cpp >> smallptSA
-	chmod +x smallptSA
-	ls -l smallptSA
-
-DISTNAME = smallpt
-DISTFILES = README LICENSE.txt smallpt.cpp smallpt4k.cpp Makefile result_25k.png
-
-SRCS = $(CFILES) $(C++FILES) $(C++FILES:%.cxx=%.h)
-
-SCENES=sky nightsky island vista overlap wada wada2 forest
-SAMPLES=8 40 200 1000 5k 25k
-IMAGES=$(SCENES:%=%.png) $(SCENES:%=%_t.jpg) $(SAMPLES:%=result_%.png) $(SAMPLES:%=result_t_%.jpg) $(SAMPLES:%=result_crop_%.jpg)
-
-WEBFILES = $(DISTFILES) smallpt.txt smallpt4k.txt index.php cpp2html.css cpp2html_9.css result640.jpg smallpt.tar.gz top.html bot.html smallpt.html extraScenes.txt $(IMAGES)
-
-results: smallpt smallpt4k
-	-mv result{,4k}_5k.ppm /tmp
-	-(time ./smallpt 5000 && mv image.ppm result_5k.ppm) 2>&1 |tee result.txt
-	-(time ./smallpt4k && mv image.ppm result4k_5k.ppm) 2>&1 |tee result4k.txt
-	-compare -metric AE result5k.ppm result_5k.png /tmp/garbage.ppm
-	-compare -metric AE result4k_5k.ppm result_5k.png /tmp/garbage.ppm
-
-results2: smallpt smallpt4k
-	for i in $(SAMPLES); do \
-	  (time ./smallpt $i && mv image.ppm image_$i.ppm) 2>&1 |tee result_$i.txt; \
-	  compare -metric AE image_$i.ppm result_$i.png /tmp/diff_$i.png; \
-        done
-
-dist: $(DISTFILES)
-	-rm -rf $(DISTNAME)
-	mkdir $(DISTNAME)
-	ln $(DISTFILES) $(DISTNAME)
-	tar cvzf $(DISTNAME).tar.gz $(DISTNAME)
-	chmod 644 $(DISTNAME).tar.gz
-	-rm -rf $(DISTNAME)
-
-thumbs:
-	for i in $(SCENES); do \
-convert "$${i}.ppm" "$${i}.png"; \
-convert -quality 92 -resize 150x150 "$${i}.ppm" "$${i}_t.jpg"; \
-done
-
-errthumbs:
-	-for i in 8 40 200 1000 5k 25k; do \
-convert "image_$${i}.ppm" -quality 92 "result_$${i}.png"; \
-convert "image_$${i}.ppm" -resize 100x100 -quality 92 "result_t_$${i}.jpg"; \
-convert "image_$${i}.ppm" -crop 100x100+205+530 -quality 92 "result_crop_$${i}.jpg"; \
-        done
-	-convert -resize 640x640 -quality 92 result_25k.png result640.jpg
-
-upload: $(WEBFILES)
-	perl -pe 's/\n/ \n/g' smallpt.cpp >/tmp/tmp.smallpt.cpp
-	source-highlight -f html -c cpp2html.css --no-doc -n -i /tmp/tmp.smallpt.cpp -o /tmp/tmp.smallpt.html
-	./source-highlight-finish.pl < /tmp/tmp.smallpt.html > smallpt.html
-	convert -geometry 640x640 -quality 92 result_25k.png result640.jpg
-	rsync -e ssh -avz $(WEBFILES) kbeason@kevinbeason.com:kevinbeason.com/smallpt/
-
-others:
-	$(CPP) $(CPPFLAGS)
+	rm -rf ${PROG}
 
+.PHONY: clean
diff --git a/smallpt.cpp b/smallpt.cpp
index aee3c5c..c5aee6b 100644
--- a/smallpt.cpp
+++ b/smallpt.cpp
@@ -9,7 +9,7 @@
 
 typedef double F;
 inline F Fsqrt (F x) { return sqrt(x); }
-inline F Fabs (F x) { return fabs(x); }
+inline F Fabs  (F x) { return fabs(x); }
 
 struct Vec {
   F x; // position, also color (r,g,b)
@@ -71,18 +71,23 @@ struct Ray {
 };
 
 // material types, used in radiance()
-enum Refl_t { DIFF, SPEC, REFR };
+enum ReflectionType {
+  DIFFUSE,
+  SPECULAR,
+  REFRACTION
+};
 
 struct Sphere {
   F radius;            // radius
   Vec position;
   Vec emission;
   Vec color;
-  Refl_t reflection_type;
+  ReflectionType 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_)
+          ReflectionType reflection_type_)
+    : radius(radius_), position(position_), emission(emission_),
+      color(color_), reflection_type(reflection_type_)
   {}
 
   // returns distance, 0 if nohit
@@ -91,7 +96,7 @@ struct Sphere {
     /* Solve :
      * t^2 * d . d + 2 * t * (o - p) . d + (o - p) . (o - p) - R^2 = 0
      */
-    Vec op = p - r.origin;
+    Vec op = position - r.origin;
     F t;
     F eps = 1e-4;
     F b = op.dot(r.direction);
@@ -100,21 +105,63 @@ struct Sphere {
       return 0;
     else
       det = Fsqrt(det);
-    return (t = b-det) > eps ? t : ((t = b + det) > eps ? t : 0);
+    t = b - det;
+    if (t > eps)
+      return t;
+    t = b + det;
+    if (t > eps)
+      return t;
+    return 0;
   }
 };
 
 //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
-  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
+  Sphere(1e5,                           // Left
+         Vec(1e5 + 1, 40.8, 81.6),
+         Vec(),
+         Vec(0.75, 0.25, 0.25),
+         DIFFUSE),
+  Sphere(1e5,                           // Right
+         Vec(-1e5 + 99, 40.8, 81.6),
+         Vec(),
+         Vec(0.25, 0.25, 0.75),
+         DIFFUSE),
+  Sphere(1e5,                           // Back
+         Vec(50, 40.8, 1e5),
+         Vec(),
+         Vec(0.75, 0.75, 0.75),
+         DIFFUSE),
+  Sphere(1e5,                           // Front
+         Vec(50, 40.8, -1e5 + 170),
+         Vec(),
+         Vec(),
+         DIFFUSE),
+  Sphere(1e5,                           // Bottom
+         Vec(50, 1e5, 81.6),
+         Vec(),
+         Vec(0.75, 0.75, 0.75),
+         DIFFUSE),
+  Sphere(1e5,                           // Top
+         Vec(50, -1e5 + 81.6, 81.6),
+         Vec(),
+         Vec(0.75, 0.75, 0.75),
+         DIFFUSE),
+  Sphere(16.5,                          // Mirror ball
+         Vec(27, 16.5, 47),
+         Vec(),
+         Vec(1, 1, 1) * 0.999,
+         SPECULAR),
+  Sphere(16.5,                          // Glass
+         Vec(73, 16.5, 78),
+         Vec(),
+         Vec(1, 1, 1) * 0.999,
+         REFRACTION),
+  Sphere(600,                           // Light
+         Vec(50, 681.6 - 0.27, 81.6),
+         Vec(12, 12, 12),
+         Vec(),
+         DIFFUSE)
 };
 
 inline F clamp (F x)
@@ -134,7 +181,7 @@ inline bool intersect (const Ray &r, F &t, int &id)
   F inf = 1e20;
   t = inf;
   while (i--) {
-    if ((d = spheres[i].intersect(r)) && d < t) {
+    if ((d = g_spheres[i].intersect(r)) && d < t) {
       t = d;
       id = i;
     }
@@ -148,18 +195,19 @@ Vec radiance (const Ray &r, int depth, unsigned short *Xi)
   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
+  const Sphere &obj = g_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;
+  Vec f = obj.color;
   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 (++depth > 5) {
     if (erand48(Xi) < p)
       f = f * (1 / p);
     else
       return obj.emission; //R.R.
-  if (obj.reflection_type == DIFF) {    // Ideal diffuse reflection
+  }
+  if (obj.reflection_type == DIFFUSE) {
     F r1 = 2 * M_PI * erand48(Xi);
     F r2 = erand48(Xi);
     F r2s = Fsqrt(r2);
@@ -171,40 +219,43 @@ Vec radiance (const Ray &r, int depth, unsigned short *Xi)
              w * Fsqrt(1 - r2)).normalize();
     return obj.emission + f.mult(radiance(Ray(x, d), depth, Xi));
   }
-  else if (obj.refl == SPEC)            // Ideal specular reflection
+  else if (obj.reflection_type == SPECULAR) {
     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);
+  }
+  else if (obj.reflection_type == 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);
+  }
+  return Vec();
 }
 
 int main (int argc, char *argv[])
@@ -226,8 +277,9 @@ int main (int argc, char *argv[])
             samples * 4,
             100.0 * y / (h - 1));
     x = 0;
-    *Xi = (unsigned short [])
-      {0, 0, static_cast<unsigned short>(y * y * y)};
+    Xi[0] = 0;
+    Xi[1] = 0;
+    Xi[2] = static_cast<unsigned short>(y * y * y);
     while (x < w) {                     // Loop over image columns
       int sy = 0;
       int i = (h - y - 1) * w + x;