 fabiensanglard.net Go back Open original

# Raytracing  I recently came across Paul Heckbert's business card raytracer. For those that have never heard of it: It is a very famous challenge in the Computer Graphics field that started on May 4th, 1984 via a post on `comp.graphics` by Paul Heckbert ( More about this in his article "A Minimal Ray Tracer" from the book Graphics Gems IV). The goal was to produce the source code for a raytracer...that would fit on the back of a business card.

Andrew Kensler's version is mesmerizing and one of the greatest hack that I have seen. Since I am curious, I decided to deconstruct it: Here is what I understood.

Edit : Andrew Kensler himself commented on Hacker News and adressed/elaborated on my observations: Thanks you so much Andrew :) !

` #include <stdlib.h> // card > aek.ppm #include <stdio.h> #include <math.h> typedef int i;typedef float f;struct v{ f x,y,z;v operator+(v r){return v(x+r.x ,y+r.y,z+r.z);}v operator*(f r){return v(x*r,y*r,z*r);}f operator%(v r){return x*r.x+y*r.y+z*r.z;}v(){}v operator^(v r ){return v(y*r.z-z*r.y,z*r.x-x*r.z,x*r. y-y*r.x);}v(f a,f b,f c){x=a;y=b;z=c;}v operator!(){return*this*(1/sqrt(*this%* this));}};i G[]={247570,280596,280600, 249748,18578,18577,231184,16,16};f R(){ return(f)rand()/RAND_MAX;}i T(v o,v d,f &t,v&n){t=1e9;i m=0;f p=-o.z/d.z;if(.01 <p)t=p,n=v(0,0,1),m=1;for(i k=19;k--;) for(i j=9;j--;)if(G[j]&1<<k){v p=o+v(-k ,0,-j-4);f b=p%d,c=p%p-1,q=b*b-c;if(q>0 ){f s=-b-sqrt(q);if(s<t&&s>.01)t=s,n=!( p+d*t),m=2;}}return m;}v S(v o,v d){f t ;v n;i m=T(o,d,t,n);if(!m)return v(.7, .6,1)*pow(1-d.z,4);v h=o+d*t,l=!(v(9+R( ),9+R(),16)+h*-1),r=d+n*(n%d*-2);f b=l% n;if(b<0||T(h,l,t,n))b=0;f p=pow(l%r*(b >0),99);if(m&1){h=h*.2;return((i)(ceil( h.x)+ceil(h.y))&1?v(3,1,1):v(3,3,3))*(b *.2+.1);}return v(p,p,p)+S(h,r)*.5;}i main(){printf("P6 512 512 255 ");v g=!v (-6,-16,0),a=!(v(0,0,1)^g)*.002,b=!(g^a )*.002,c=(a+b)*-256+g;for(i y=512;y--;) for(i x=512;x--;){v p(13,13,13);for(i r =64;r--;){v t=a*(R()-.5)*99+b*(R()-.5)* 99;p=S(v(17,16,8)+t,!(t*-1+(a*(R()+x)+b *(y+R())+c)*16))*3.5+p;}printf("%c%c%c" ,(i)p.x,(i)p.y,(i)p.z);}} `
The code looks confusing but it compiles and run flawlessly !! On my MacBook Air, I can save this file on my Desktop as `card.cpp`, open Terminal and type :
` c++ -O3 -o card card.cpp ./card > card.ppm `

Within 27 seconds, the following image is generated: The list of features is impressive :

• The world is made of spheres which are carefully organized.
• There is a floor and it is textured.
• There is a sky and it has a nice gradient.
• OMG, Depth of View blur: Are you kidding me ?

All that fitting on the back of a business card ! Let's see how it works !

## Vector class

 ` #include // card > aek.ppm #include #include typedef int i;typedef float f;struct v{ f x,y,z;v operator+(v r){return v(x+r.x ,y+r.y,z+r.z);}v operator*(f r){return v(x*r,y*r,z*r);}f operator%(v r){return x*r.x+y*r.y+z*r.z;}v(){}v operator^(v r ){return v(y*r.z-z*r.y,z*r.x-x*r.z,x*r. y-y*r.x);}v(f a,f b,f c){x=a;y=b;z=c;}v operator!(){return*this*(1/sqrt(*this%* this));}};i G[]={247570,280596,280600, 249748,18578,18577,231184,16,16};f R(){ return(f)rand()/RAND_MAX;}i T(v o,v d,f &t,v&n){t=1e9;i m=0;f p=-o.z/d.z;if(.01 0 ){f s=-b-sqrt(q);if(s.01)t=s,n=!( p+d*t),m=2;}}return m;}v S(v o,v d){f t ;v n;i m=T(o,d,t,n);if(!m)return v(.7, .6,1)*pow(1-d.z,4);v h=o+d*t,l=!(v(9+R( ),9+R(),16)+h*-1),r=d+n*(n%d*-2);f b=l% n;if(b<0||T(h,l,t,n))b=0;f p=pow(l%r*(b >0),99);if(m&1){h=h*.2;return((i)(ceil( h.x)+ceil(h.y))&1?v(3,1,1):v(3,3,3))*(b *.2+.1);}return v(p,p,p)+S(h,r)*.5;}i main(){printf("P6 512 512 255 ");v g=!v (-6,-16,0),a=!(v(0,0,1)^g)*.002,b=!(g^a )*.002,c=(a+b)*-256+g;for(i y=512;y--;) for(i x=512;x--;){v p(13,13,13);for(i r =64;r--;){v t=a*(R()-.5)*99+b*(R()-.5)* 99;p=S(v(17,16,8)+t,!(t*-1+(a*(R()+x)+b *(y+R())+c)*16))*3.5+p;}printf("%c%c%c" ,(i)p.x,(i)p.y,(i)p.z);}} ` The main trick that helps reducing the code size dramatically it to use two `define` : Define `i` as `int`. Define `f` as `float`. The other big space saver is the 'v' C++ class which is used for vector but also pixel processing. Here is the code formated, highlighted and commented: ` #include #include #include typedef int i; typedef float f; struct v{ f x,y,z; v operator+(v r){return v(x+r.x,y+r.y,z+r.z);} v operator*(f r){return v(x*r,y*r,z*r);} f operator%(v r){return x*r.x+y*r.y+z*r.z;} v(){} v operator^(v r){return v(y*r.z-z*r.y,z*r.x-x*r.z,x*r.y-y*r.x);} v(f a,f b,f c){x=a;y=b;z=c;} v operator!(){return *this*(1 /sqrt(*this%*this));} }; `

## Random() and World database

 ` #include // card > aek.ppm #include #include typedef int i;typedef float f;struct v{ f x,y,z;v operator+(v r){return v(x+r.x ,y+r.y,z+r.z);}v operator*(f r){return v(x*r,y*r,z*r);}f operator%(v r){return x*r.x+y*r.y+z*r.z;}v(){}v operator^(v r ){return v(y*r.z-z*r.y,z*r.x-x*r.z,x*r. y-y*r.x);}v(f a,f b,f c){x=a;y=b;z=c;}v operator!(){return*this*(1/sqrt(*this%* this));}};i G[]={247570,280596,280600, 249748,18578,18577,231184,16,16};f R(){ return(f)rand()/RAND_MAX;}i T(v o,v d,f &t,v&n){t=1e9;i m=0;f p=-o.z/d.z;if(.01 0 ){f s=-b-sqrt(q);if(s.01)t=s,n=!( p+d*t),m=2;}}return m;}v S(v o,v d){f t ;v n;i m=T(o,d,t,n);if(!m)return v(.7, .6,1)*pow(1-d.z,4);v h=o+d*t,l=!(v(9+R( ),9+R(),16)+h*-1),r=d+n*(n%d*-2);f b=l% n;if(b<0||T(h,l,t,n))b=0;f p=pow(l%r*(b >0),99);if(m&1){h=h*.2;return((i)(ceil( h.x)+ceil(h.y))&1?v(3,1,1):v(3,3,3))*(b *.2+.1);}return v(p,p,p)+S(h,r)*.5;}i main(){printf("P6 512 512 255 ");v g=!v (-6,-16,0),a=!(v(0,0,1)^g)*.002,b=!(g^a )*.002,c=(a+b)*-256+g;for(i y=512;y--;) for(i x=512;x--;){v p(13,13,13);for(i r =64;r--;){v t=a*(R()-.5)*99+b*(R()-.5)* 99;p=S(v(17,16,8)+t,!(t*-1+(a*(R()+x)+b *(y+R())+c)*16))*3.5+p;}printf("%c%c%c" ,(i)p.x,(i)p.y,(i)p.z);}} ` The next section also saves a lot of code by defining a `R` function that return a random value within `[0.0f-1.0f]`. It is very useful for the stochastic sampling resulting in blur and soft-shadows. The `G` array encodes the location of the spheres. The integers are interpreted as 9 lines and 19 columns bit vector. Here is the code formated, highlighted and commented : ` i G[]={247570,280596,280600,249748,18578,18577,231184,16,16}; f R(){return(f)rand()/RAND_MAX;} `

## Main method

 ` #include // card > aek.ppm #include #include typedef int i;typedef float f;struct v{ f x,y,z;v operator+(v r){return v(x+r.x ,y+r.y,z+r.z);}v operator*(f r){return v(x*r,y*r,z*r);}f operator%(v r){return x*r.x+y*r.y+z*r.z;}v(){}v operator^(v r ){return v(y*r.z-z*r.y,z*r.x-x*r.z,x*r. y-y*r.x);}v(f a,f b,f c){x=a;y=b;z=c;}v operator!(){return*this*(1/sqrt(*this%* this));}};i G[]={247570,280596,280600, 249748,18578,18577,231184,16,16};f R(){ return(f)rand()/RAND_MAX;}i T(v o,v d,f &t,v&n){t=1e9;i m=0;f p=-o.z/d.z;if(.01 0 ){f s=-b-sqrt(q);if(s.01)t=s,n=!( p+d*t),m=2;}}return m;}v S(v o,v d){f t ;v n;i m=T(o,d,t,n);if(!m)return v(.7, .6,1)*pow(1-d.z,4);v h=o+d*t,l=!(v(9+R( ),9+R(),16)+h*-1),r=d+n*(n%d*-2);f b=l% n;if(b<0||T(h,l,t,n))b=0;f p=pow(l%r*(b >0),99);if(m&1){h=h*.2;return((i)(ceil( h.x)+ceil(h.y))&1?v(3,1,1):v(3,3,3))*(b *.2+.1);}return v(p,p,p)+S(h,r)*.5;}i main(){printf("P6 512 512 255 ");v g=!v (-6,-16,0),a=!(v(0,0,1)^g)*.002,b=!(g^a )*.002,c=(a+b)*-256+g;for(i y=512;y--;) for(i x=512;x--;){v p(13,13,13);for(i r =64;r--;){v t=a*(R()-.5)*99+b*(R()-.5)* 99;p=S(v(17,16,8)+t,!(t*-1+(a*(R()+x)+b *(y+R())+c)*16))*3.5+p;}printf("%c%c%c" ,(i)p.x,(i)p.y,(i)p.z);}} ` The main method uses a simple image file format known as PPM. It is text based with a ridiculously simple header: `P6 [WIDTH] [HEIGHT] [MAX_VALUE]` followed by the RGB values for each pixels. For each 512x512 pixels in the image, the program samples (S) the color of 64 rays, accumulates the result and write it to the output (stdout). It applies a bit of random delta on each ray origin and ray direction in order to produce a nice Depth of View effect. Here is the code formated, highlighted and commented : ` i main(){ printf("P6 512 512 255 "); v g=!v(-6,-16,0), a=!(v(0,0,1)^g)*.002, b=!(g^a)*.002, c=(a+b)*-256+g; for(i y=512;y--;) for(i x=512;x--;){ v p(13,13,13); // Default pixel color is almost pitch black for(i r=64;r--;){ v t=a*(R()-.5)*99+b*(R()-.5)*99; p=S(v(17,16,8)+t, !(t*-1+(a*(R()+x)+b*(y+R())+c)*16) )*3.5+p; } printf("%c%c%c",(i)p.x,(i)p.y,(i)p.z); } } `

## Sampler

 ` #include // card > aek.ppm #include #include typedef int i;typedef float f;struct v{ f x,y,z;v operator+(v r){return v(x+r.x ,y+r.y,z+r.z);}v operator*(f r){return v(x*r,y*r,z*r);}f operator%(v r){return x*r.x+y*r.y+z*r.z;}v(){}v operator^(v r ){return v(y*r.z-z*r.y,z*r.x-x*r.z,x*r. y-y*r.x);}v(f a,f b,f c){x=a;y=b;z=c;}v operator!(){return*this*(1/sqrt(*this%* this));}};i G[]={247570,280596,280600, 249748,18578,18577,231184,16,16};f R(){ return(f)rand()/RAND_MAX;}i T(v o,v d,f &t,v&n){t=1e9;i m=0;f p=-o.z/d.z;if(.01 0 ){f s=-b-sqrt(q);if(s.01)t=s,n=!( p+d*t),m=2;}}return m;}v S(v o,v d){f t ;v n;i m=T(o,d,t,n);if(!m)return v(.7, .6,1)*pow(1-d.z,4);v h=o+d*t,l=!(v(9+R( ),9+R(),16)+h*-1),r=d+n*(n%d*-2);f b=l% n;if(b<0||T(h,l,t,n))b=0;f p=pow(l%r*(b >0),99);if(m&1){h=h*.2;return((i)(ceil( h.x)+ceil(h.y))&1?v(3,1,1):v(3,3,3))*(b *.2+.1);}return v(p,p,p)+S(h,r)*.5;}i main(){printf("P6 512 512 255 ");v g=!v (-6,-16,0),a=!(v(0,0,1)^g)*.002,b=!(g^a )*.002,c=(a+b)*-256+g;for(i y=512;y--;) for(i x=512;x--;){v p(13,13,13);for(i r =64;r--;){v t=a*(R()-.5)*99+b*(R()-.5)* 99;p=S(v(17,16,8)+t,!(t*-1+(a*(R()+x)+b *(y+R())+c)*16))*3.5+p;}printf("%c%c%c" ,(i)p.x,(i)p.y,(i)p.z);}} ` The Sampler (S) is a function that returns a pixel color for a given Point Origin (o) and Vector Direction (d). It is recursive if it hits a sphere. If no sphere is hit by the ray (and a ray always ends up missing) then the function either returns a sky color gradient or a floor color based on a checkerboard texture. Remark the calls to `R` when computing the light direction: This is how the soft-shadows are generated. Here is the code formated, highlighted and commented : ` v S(v o,v d){ f t; v n; i m=T(o,d,t,n); if(!m) return v(.7,.6,1)*pow(1-d.z,4); v h=o+d*t, l=!(v(9+R(),9+R(),16)+h*-1), r=d+n*(n%d*-2); f b=l%n; if(b<0<>||T(h,l,t,n)) b=0; f p=pow(l%r*(b>0),99); if(m&1){ h=h*.2; return((i)(ceil(h.x)+ceil(h.y))&1?v(3,1,1):v(3,3,3))*(b*.2+.1); } return v(p,p,p)+S(h,r)*.5; } `

## Tracer

 ` #include // card > aek.ppm #include #include typedef int i;typedef float f;struct v{ f x,y,z;v operator+(v r){return v(x+r.x ,y+r.y,z+r.z);}v operator*(f r){return v(x*r,y*r,z*r);}f operator%(v r){return x*r.x+y*r.y+z*r.z;}v(){}v operator^(v r ){return v(y*r.z-z*r.y,z*r.x-x*r.z,x*r. y-y*r.x);}v(f a,f b,f c){x=a;y=b;z=c;}v operator!(){return*this*(1/sqrt(*this%* this));}};i G[]={247570,280596,280600, 249748,18578,18577,231184,16,16};f R(){ return(f)rand()/RAND_MAX;}i T(v o,v d,f &t,v&n){t=1e9;i m=0;f p=-o.z/d.z;if(.01 0 ){f s=-b-sqrt(q);if(s.01)t=s,n=!( p+d*t),m=2;}}return m;}v S(v o,v d){f t ;v n;i m=T(o,d,t,n);if(!m)return v(.7, .6,1)*pow(1-d.z,4);v h=o+d*t,l=!(v(9+R( ),9+R(),16)+h*-1),r=d+n*(n%d*-2);f b=l% n;if(b<0||T(h,l,t,n))b=0;f p=pow(l%r*(b >0),99);if(m&1){h=h*.2;return((i)(ceil( h.x)+ceil(h.y))&1?v(3,1,1):v(3,3,3))*(b *.2+.1);}return v(p,p,p)+S(h,r)*.5;}i main(){printf("P6 512 512 255 ");v g=!v (-6,-16,0),a=!(v(0,0,1)^g)*.002,b=!(g^a )*.002,c=(a+b)*-256+g;for(i y=512;y--;) for(i x=512;x--;){v p(13,13,13);for(i r =64;r--;){v t=a*(R()-.5)*99+b*(R()-.5)* 99;p=S(v(17,16,8)+t,!(t*-1+(a*(R()+x)+b *(y+R())+c)*16))*3.5+p;}printf("%c%c%c" ,(i)p.x,(i)p.y,(i)p.z);}} ` The Tracer is in charge of casting a Ray (Origin o ,Direction d). It return an integer that is a code for the intersection test with the spheres in the world (0=miss toward sky,1=miss toward floor, 2=sphere hit). If a sphere is hit, it updates the references t (the parametric value to compute the distance of intersection) and n (the half-vector where the ray bounce). Here is the code formated, highlighted and commented : ` i T(v o,v d,f& t,v& n){ t=1e9; i m=0; f p=-o.z/d.z; if(.010){ f s=-b-sqrt(q); if(s.01) t=s, n=!(p+d*t), m=2; } } return m; } `

## The Leet Number

Many readers have tried to shorten the code further with more or less success. The original author stopped with the current version for a very interesting reason:

` fabien\$ wc card.cpp 35 95 1337 card.cpp `

The total size of the code is 1337 bytes: The Leet number. Beat that !

## Script kiddie

And now that the business card is fully understood, we can generate our own :

` #include <stdlib.h> // card > aek.ppm #include <stdio.h> #include <math.h> typedef int i;typedef float f;struct v{ f x,y,z;v operator+(v r){return v(x+r.x ,y+r.y,z+r.z);}v operator*(f r){ return v(x*r,y*r,z*r);}f operator%(v r){return x*r.x+y*r.y+z*r.z;}v(){}v operator^(v r ){return v(y*r.z-z*r.y,z*r.x-x*r.z,x*r. y-y*r.x);}v(f a,f b,f c){x=a;y=b;z=c;}v operator!(){return*this*(1/sqrt(*this%* this));}};i G[]={133022, 133266,133266, 133022, 254096, 131216, 131984, 131072, 258048,};f R(){return(f)rand()/RAND_MAX ;}i T(v o,v d,f&t,v&n){t=1e9;i m=0;f p= -o.z/d.z; if(.01<p)t=p, n=v(0,0,1),m=1; for(i k=19;k--;)for(i j=9;j--;)if(G[j]& 1<<k){v p=o+v(-k,0,-j-4);f b=p%d,c=p%p- 1,q=b*b-c;if(q>0){f s=-b-sqrt(q);if(s<t &&s>.01)t=s,n=!(p+d*t),m=2;}}return m;} v S(v o,v d){f t;v n;i m=T(o,d,t,n);if( !m)return v(.7,.6,1)*pow(1-d.z,4);v h=o +d*t,l=!(v(9+R(),9+R(),16)+h*-1),r=d+n* (n%d*-2);f b=l%n;if(b<0||T(h,l,t,n))b=0 ;f p=pow(l%r*(b>0),99);if(m&1){h=h*.2; return((i)(ceil(h.x)+ceil(h.y))&1?v(3,1 ,1):v(3,3,3))*(b*.2+.1);}return v(p,p,p )+S(h,r)*.5;}i main(){printf("P6 512 " "512 255 ");v g=!v(-6,-16,0),a=!(v(0,0, 1)^g)*.002,b=!(g^a)*.002,c=(a+b)*-256+g ;for(i y=512;y--;)for(i x=512;x--;){v p (9,9,9);for(i r=64;r--;){v t=a*(R()-.5) *99+b*(R()-.5)*99;p=S(v(17,16,8)+t,!(t* -1+(a*(R()+x)+b*(y+R())+c)*16))*3.5+p;} printf("%c%c%c",(i)p.x,(i)p.y,(i)p.z);}} ` 