Postcard-sized ray tracer decoding
- Transfer
“He did it again!” - this is what first came to my mind when I looked at the back of the Pixar flyer [1] , which was completely filled with code. The accumulation of constructions and expressions was signed in the lower right corner by none other than Andrew Kensler. For those who do not know him, I will say: Andrew is a programmer who invented in 2009 a 1337-byte ray tracer the size of a business card .
This time, Andrew came up with something more voluminous, but with a much more interesting visual result. Since I finished writing my Game Engine Black Books about Wolf3D and DOOMI had time to study the insides of his mysterious code. And almost immediately I was literally fascinated by the techniques found in it. They were very different from Andrew’s previous work, based on the “standard” ray tracer. I was interested in learning about ray marching, the functions of constructive volumetric geometry, Monte Carlo rendering / tracing paths, as well as many other tricks he used to squeeze code into such a small piece of paper.
Source
The front of the flyer is an advertisement for the Pixar recruitment department. Printed on the back side are 3,037 bytes of C ++ code obfuscated to occupy as little surface as possible.
#include<stdlib.h> // card > pixar.ppm#include<stdio.h>#include<math.h>#define R return#define O operatortypedeffloat F;typedefint I;structV{F x,y,z;V(F v=0){x=y=z=v;}V(F a,F b,F
c=0){x=a;y=b;z=c;}V O+(V r){R V(x+r.x,y+r.y,z+r.z);}V O*(V r){R V(x*r.x,y*r.
y,z*r.z);}F O%(V r){R x*r.x+y*r.y+z*r.z;}V O!(){R*this*(1/sqrtf(*this%*this)
);}};F L(F l,F r){R l<r?l:r;}F U(){R(F)rand()/RAND_MAX;}F B(V p,V l,V h){l=p
+l*-1;h=h+p*-1;R-L(L(L(l.x,h.x),L(l.y,h.y)),L(l.z,h.z));}F S(V p,I&m){F d=1\
e9;V f=p;f.z=0;char l[]="5O5_5W9W5_9_COC_AOEOA_E_IOQ_I_QOUOY_Y_]OWW[WaOa_aW\
eWa_e_cWiO";for(I i=0;i<60;i+=4){V b=V(l[i]-79,l[i+1]-79)*.5,e=V(l[i+2]-79,l
[i+3]-79)*.5+b*-1,o=f+(b+e*L(-L((b+f*-1)%e/(e%e),0),1))*-1;d=L(d,o%o);}d=sq\
rtf(d);V a[]={V(-11,6),V(11,6)};for(I i=2;i--;){V o=f+a[i]*-1;d=L(d,o.x>0?f\
absf(sqrtf(o%o)-2):(o.y+=o.y>0?-2:2,sqrtf(o%o)));}d=powf(powf(d,8)+powf(p.z,
8),.125)-.5;m=1;F r=L(-L(B(p,V(-30,-.5,-30),V(30,18,30)),B(p,V(-25,17,-25),V
(25,20,25))),B(V(fmodf(fabsf(p.x),8),p.y,p.z),V(1.5,18.5,-25),V(6.5,20,25)))
;if(r<d)d=r,m=2;F s=19.9-p.y;if(s<d)d=s,m=3;R d;}I M(V o,V d,V&h,V&n){I m,s=
0;F t=0,c;for(;t<100;t+=c)if((c=S(h=o+d*t,m))<.01||++s>99)R n=!V(S(h+V(.01,0
),s)-c,S(h+V(0,.01),s)-c,S(h+V(0,0,.01),s)-c),m;R 0;}V T(V o,V d){V h,n,r,t=
1,l(!V(.6,.6,1));for(I b=3;b--;){I m=M(o,d,h,n);if(!m)break;if(m==1){d=d+n*(
n%d*-2);o=h+d*.1;t=t*.2;}if(m==2){F i=n%l,p=6.283185*U(),c=U(),s=sqrtf(1-c),
g=n.z<0?-1:1,u=-1/(g+n.z),v=n.x*n.y*u;d=V(v,g+n.y*n.y*u,-n.y)*(cosf(p)*s)+V(
1+g*n.x*n.x*u,g*v,-g*n.x)*(sinf(p)*s)+n*sqrtf(c);o=h+d*.1;t=t*.2;if(i>0&&M(h
+n*.1,l,h,n)==3)r=r+t*V(500,400,100)*i;}if(m==3){r=r+t*V(50,80,100);break;}}
R r;}I main(){I w=960,h=540,s=16;V e(-22,5,25),g=!(V(-3,4,0)+e*-1),l=!V(g.z,
0,-g.x)*(1./w),u(g.y*l.z-g.z*l.y,g.z*l.x-g.x*l.z,g.x*l.y-g.y*l.x);printf("P\
6 %d %d 255 ",w,h);for(I y=h;y--;)for(I x=w;x--;){V c;for(I p=s;p--;)c=c+T(e
,!(g+l*(x-w/2+U())+u*(y-h/2+U())));c=c*(1./s)+14./241;V o=c+1;c=V(c.x/o.x,c.
y/o.y,c.z/o.z)*255;printf("%c%c%c",(I)c.x,(I)c.y,(I)c.z);}}// Andrew KenslerDoes he even work?
With the code there is an instruction for its launch. The idea is to redirect standard output to a file. By extension, we can assume that the output format is a text-based image format called NetPBM [2] .
$ clang -o card2 -O3 raytracer.cpp $ time ./card> pixar.ppm real 2m58.524s user 2m57.567s sys 0m0.415s
After two minutes and fifty eight seconds [3] the next image is generated. It's amazing how little code is required for it.
Much can be learned from the above image. Grit is an obvious sign of a path tracer. This type of renderer differs from raytracing in that the rays are not traced back to the light sources. In this method, thousands of rays per pixel are emitted from sources and the program watches them, hoping that they will find the source of illumination. This is an interesting technique that is much better than ray tracing, coping with the rendering of ambient occlusion, soft shadows, caustics and radiosity.
Break the code into parts
Passing input to CLion formats the code (see the output here ) and breaks it into smaller parts / tasks.
#include <stdlib.h> // card> pixar.ppm #include <stdio.h> #include <math.h>
#define R return #define O operator typedef float F; typedef int I;
struct V {F x, y, z; V (F v = 0) {x = y = z = v;} V (F a, F b, F
c = 0) {x = a; y = b; z = c;} V O + (V r) {RV (x + rx, y + ry, z + rz);} VO * (V r) {RV ( x * rx, y * r.
y, z * rz);} FO% (V r) {R x * r.x + y * r.y + z * rz;} VO! () {R * this * (1 / sqrtf (* this% * this)
);}};FL (F l, F r) {R l <r? L: r;} FU () {R (F) rand () / RAND_MAX;} FB (V p, V l, V h) {l = p
+ l * -1; h = h + p * -1; RL (L (L (lx, hx), L (ly, hy)), L (lz, hz));}FS (V p, I & m) {F d = 1 \
e9; V f = p; fz = 0; char l [] = "5O5_5W9W5_9_COC_AOEOA_E_IOQ_I_QOUOY_Y_] OWW [WaOa_aW \
eWa_e_cWiO "; for (I i = 0; i <60; i + = 4) {V b = V (l [i] -79, l [i + 1] -79) *. 5, e = V (l [ i + 2] -79, l
[i + 3] -79) *. 5 + b * -1, o = f + (b + e * L (-L ((b + f * -1)% e / (e% e), 0), 1)) * - 1; d = L (d, o% o);} d = sq \
rtf (d); V a [] = {V (-11,6), V (11,6)}; for (I i = 2; i -;) {V o = f + a [i] * -1; d = L (d, ox> 0? F \
absf (sqrtf (o% o) -2) :( o.y + = oy> 0? -2: 2, sqrtf (o% o)));} d = powf (powf (d, 8) + powf (pz ,
8),. 125) -. 5; m = 1; F r = L (-L (B (p, V (-30, -. 5, -30), V (30,18,30)), B (p, V (-25,17, -25), V
(25,20,25))), B (V (fmodf (fabsf (px), 8), py, pz), V (1.5,18.5, -25), V (6.5,20,25)))
; if (r <d) d = r, m = 2; F s = 19.9-py; if (s <d) d = s, m = 3; R d;}IM (V o, V d, V & h, V & n) {I m, s =
0; F t = 0, c; for (; t <100; t + = c) if ((c = S (h = o + d * t, m)) <. 01 || ++ s> 99) R n =! V (S (h + V (.01,0
), s) -c, S (h + V (0, .01), s) -c, S (h + V (0,0, .01), s) -c), m; R 0;}VT (V o, V d) {V h, n, r, t =
1, l (! V (.6, .6,1)); for (I b = 3; b -;) {I m = M (o, d, h, n); if (! M) break ; if (m == 1) {d = d + n * (
n% d * -2); o = h + d * .1; t = t * .2;} if (m == 2) {F i = n% l, p = 6.283185 * U (), c = U (), s = sqrtf (1-c),
g = nz <0? -1: 1, u = -1 / (g + nz), v = nx * ny * u; d = V (v, g + ny * ny * u, -ny) * (cosf (p) * s) + V (
1 + g * nx * nx * u, g * v, -g * nx) * (sinf (p) * s) + n * sqrtf (c); o = h + d * .1; t = t *. 2; if (i> 0 && M (h
+ n * .1, l, h, n) == 3) r = r + t * V (500,400,100) * i;} if (m == 3) {r = r + t * V (50,80,100) ; break;}}
R r;}I main () {I w = 960, h = 540, s = 16; V e (-22,5,25), g =! (V (-3,4,0) + e * -1), l =! V (gz,
0, -gx) * (1./w), u (gy * lz-gz * ly, gz * lx-gx * lz, gx * ly-gy * lx); printf ("P \
6% d% d 255 ", w, h); for (I y = h; y -;) for (I x = w; x -;) {V c; for (I p = s; p- -;) c = c + T (e
,! (g + l * (xw / 2 + U ()) + u * (yh / 2 + U ()))); c = c * (1./s) + 14. / 241; V o = c + 1; c = V (cx / ox, c.
y / oy, cz / oz) * 255; printf ("% c% c% c", (I) cx, (I) cy, (I) cz);}}// Andrew Kensler
Each of the sections is described in detail in the rest of the article:
■ - simple tricks, ■ is the Vector class, ■ is the auxiliary code, ■ is the database, ■ is Ray marching, ■ is the sampling, ■ is the main code.
Common tricks with #define and typedef
Common tricks are using #define and typedef to significantly reduce the amount of code. Here we denote F = float, I = int, R = return, and O = operator. Reverse engineering is trivial.
Class V
Next comes class V, which I renamed to Vec (even though, as we will see below, it is also used to store RGB channels in float format).
structVec {float x, y, z;
Vec(float v = 0) { x = y = z = v; }
Vec(float a, float b, float c = 0) { x = a; y = b; z = c;}
Vec operator+(Vec r) { return Vec(x + r.x, y + r.y, z + r.z); }
Vec operator*(Vec r) { return Vec(x * r.x, y * r.y, z * r.z); }
// dot productfloatoperator%(Vec r) { return x * r.x + y * r.y + z * r.z; }
// inverse square root
Vec operator!() {return *this * (1 / sqrtf(*this % *this) );}
};Note that there is no subtraction operator (-), so instead of writing "X = A - B", you use "X = A + B * -1". The inverse square root is further useful for normalizing vectors.
Main function
main () is the only character that cannot be obfuscated, because it is called by the _start function of the libc library. Usually it is worth starting with it, because it will be easier to work this way. It took me some time to figure out the meanings of the first letters, but still managed to create something readable.
intmain(){
int w = 960, h = 540, samplesCount = 16;
Vec position(-22, 5, 25);
Vec goal = !(Vec(-3, 4, 0) + position * -1);
Vec left = !Vec(goal.z, 0, -goal.x) * (1. / w);
// Cross-product to get the up vectorVec up(goal.y * left.z - goal.z * left.y,
goal.z * left.x - goal.x * left.z,
goal.x * left.y - goal.y * left.x);
printf("P6 %d %d 255 ", w, h);
for (int y = h; y--;)
for (int x = w; x--;) {
Vec color;
for (int p = samplesCount; p--;)
color = color + Trace(position, !(goal + left * (x - w / 2 + randomVal())+
up * (y - h / 2 + randomVal())));
// Reinhard tone mapping
color = color * (1. / samplesCount) + 14. / 241;
Vec o = color + 1;
color = Vec(color.x / o.x, color.y / o.y, color.z / o.z) * 255;
printf("%c%c%c", (int) color.x, (int) color.y, (int) color.z);
}
}Note that float literals do not contain the letter “f”, and the fractional part is discarded to save space. The same trick is used below, where the integer part is dropped (float x = .5). Also unusual is the “for” construction with an iteration expression inserted into the break condition.
This is a fairly standard main function for the ray / path tracer. Here, camera vectors are defined and rays are emitted for each pixel. The difference between the ray tracer and the path tracer is that several rays are emitted per pixel in the TS, which are slightly shifted randomly. Then the color obtained for each ray in the pixel accumulates in the three float channels R, B, G. At the end, a tonal correction of the result of the Reinhard method is performed.
The most important part is sampleCount, which theoretically you can assign a value of 1 to speed up rendering and iterations. Here are examples of visualizations with sample values from 1 to 2048.
Spoiler header
1
2
4
8
16
32
64
128
256
512
1024
2048
1
2
4
8
16
32
64
128
256
512
1024
2048
Auxiliary code
Another simple code snippet is helper functions. In this case, we have a trivial min () function, a random value generator in the interval [0,1], and a much more interesting boxTest (), which is part of the Constructive Solid Geometry (CSG) system used to cut the world. CSG is covered in the next section.
floatmin(float l, float r){ return l < r ? l : r; }
floatrandomVal(){ return (float) rand() / RAND_MAX; }
// Rectangle CSG equation. Returns minimum signed distance from // space carved by lowerLeft vertex and opposite rectangle// vertex upperRight.floatBoxTest(Vec position, Vec lowerLeft, Vec upperRight){
lowerLeft = position + lowerLeft * -1;
upperRight = upperRight + position * -1;
return -min(
min(
min(lowerLeft.x, upperRight.x),
min(lowerLeft.y, upperRight.y)
),
min(lowerLeft.z, upperRight.z));
}Functions of constructive volumetric geometry
There are no vertices in the code. Everything is done using the functions of the CSG. If you are unfamiliar with them, then simply say that these are functions that describe whether the coordinate is inside or outside the object. If the function returns a positive distance, then the point is inside the object. A negative distance indicates that the point is outside the object. There are many functions for describing different objects, but for the sake of simplification, let's take a sphere and two points, A and B, for example.
// Signed distance point(p) to sphere(c,r)floattestSphere(Vec p, Vec c, float r){
Vec delta = c - p;
float distance = sqrtf(delta%delta);
return radius - distance;
}
Vec A {4, 6};
Vec B {3, 2};
Vec C {4, 2};
float r = 2.;
testSphere(A, C, r); // == -1 (outside)
testSphere(B, C, r); // == 1 (inside)The testSphere () function returns -1 for point A (that is, it is outside) and 1 for B (that is, it is inside). The signs for distances are just a trick, allowing you to get two pieces of information instead of one in the case of a single value. A similar type of function can also be written to describe a parallelogram (this is exactly what is performed in the function BoxTest function).
// Signed distance point(p) to Box(c1,c2)floattestRectangle(Vec p, Vec c1, Vec c2){
c1 = p + c1 * -1;
c2 = c2 + position * -1;
return min(
min(
min(c1.x, c2.x),
min(c1.y, c2.y)),
min(c1.z, c2.z));
}
Vec A {3, 3}; Vec B {4, 6};
Vec C1 {2, 2}; Vec C2 {5, 4};
testRectangle(A, C1, C2); // 1.41 (inside)
testRectangle(B, C1, C2); // -2.23 (outside)Now let's see what happens if we flip the sign of the return value.
// Signed distance point(p) to carved box(c1,c2)floattestCarveBox(Vec p, Vec c1, Vec c2){
c1 = p + c1 * -1;
c2 = c2 + position * -1;
return -min(
min(
min(c1.x, c2.x),
min(c1.y, c2.y)),
min(c1.z, c2.z));
}
Vec A {3, 3}; Vec B {4, 6};
Vec C1 {2, 2}; Vec C2 {5, 4};
testCarveBox(A, C1, C2); // == -1.41 (outside)
testCarveBox(B, C1, C2); // == 2.23 (inside)Now we do not describe a solid object, but declare the whole world solid and cut out the empty space in it. Functions can be used as building blocks, which, when combined with an I can combine, describe more complex shapes. Using the logical addition operator (min function), we can cut a pair of rectangles one above the other and the result will look like this.
// Signed distance point to room floattestRoom(Vec p){
Vec C1 {2, 4}; Vec C2 {5, 2}; // Lower room
Vec C3 {3, 5}; Vec C4 {4, 4}; // Upper room// min() is the union of the two carved volumes.return min(testCarvedBox(p, C1, C2),
testCarvedBox(p, C3, C4));
}
Vec A {3, 3};
Vec B {4, 6};
testRoom(A, C1, C2); // == -1.41 (outside)
testRoom(B, C1, C2); // == 1.00 (inside)If you think about it, it looks like the room we are examining, because this is how the lower room is expressed - with the help of two parallelograms cut out.
Now, having mastered the powerful knowledge of the CSG, we can return to the code and consider the function of the database, which is the hardest to deal with.
#define HIT_NONE 0#define HIT_LETTER 1#define HIT_WALL 2#define HIT_SUN 3// Sample the world using Signed Distance Fields.floatQueryDatabase(Vec position, int &hitType){
float distance = 1e9;
Vec f = position; // Flattened position (z=0)
f.z = 0;
char letters[15*4+1] = // 15 two points lines"5O5_""5W9W""5_9_"// P (without curve)"AOEO""COC_""A_E_"// I"IOQ_""I_QO"// X"UOY_""Y_]O""WW[W"// A"aOa_""aWeW""a_e_""cWiO"; // R (without curve)for (int i = 0; i < sizeof(letters); i += 4) {
Vec begin = Vec(letters[i] - 79, letters[i + 1] - 79) * .5;
Vec e = Vec(letters[i + 2] - 79, letters[i + 3] - 79) * .5 + begin * -1;
Vec o = f + (begin + e * min(-min((begin + f * -1) % e / (e % e),
0),
1)
) * -1;
distance = min(distance, o % o); // compare squared distance.
}
distance = sqrtf(distance); // Get real distance, not square distance.// Two curves (for P and R in PixaR) with hard-coded locations.
Vec curves[] = {Vec(-11, 6), Vec(11, 6)};
for (int i = 2; i--;) {
Vec o = f + curves[i] * -1;
distance = min(distance,
o.x > 0 ? fabsf(sqrtf(o % o) - 2)
: (o.y += o.y > 0 ? -2 : 2, sqrtf(o % o))
);
}
distance = powf(powf(distance, 8) + powf(position.z, 8), .125) - .5;
hitType = HIT_LETTER;
float roomDist ;
roomDist = min(// min(A,B) = Union with Constructive solid geometry//-min carves an empty space
-min(// Lower room
BoxTest(position, Vec(-30, -.5, -30), Vec(30, 18, 30)),
// Upper room
BoxTest(position, Vec(-25, 17, -25), Vec(25, 20, 25))
),
BoxTest( // Ceiling "planks" spaced 8 units apart.
Vec(fmodf(fabsf(position.x), 8),
position.y,
position.z),
Vec(1.5, 18.5, -25),
Vec(6.5, 20, 25)
)
);
if (roomDist < distance) distance = roomDist, hitType = HIT_WALL;
float sun = 19.9 - position.y ; // Everything above 19.9 is light source.if (sun < distance)distance = sun, hitType = HIT_SUN;
return distance;
}You can see here the function of “cutting out” a parallelogram in which only two rectangles are used to build the whole room (our brain does the rest, it represents that these are walls). A horizontal ladder is a slightly more complex CSG function using division with remainder. And finally, the letters of the word PIXAR are composed of 15 lines with a pair of "origin / delta" and two special cases for the curves in the letters P and R.
Ray Marching
Having a database of CSG functions describing the world, it is enough for us to skip all the rays emitted in the main () function. In ray marching, the distance function is used. This means that the sampling position shifts forward to the nearest obstacle.
// Perform signed sphere marching// Returns hitType 0, 1, 2, or 3 and update hit position/normalintRayMarching(Vec origin, Vec direction, Vec &hitPos, Vec &hitNorm){
int hitType = HIT_NONE;
int noHitCount = 0;
float d; // distance from closest object in world.// Signed distance marchingfor (float total_d=0; total_d < 100; total_d += d)
if ((d = QueryDatabase(hitPos = origin + direction * total_d, hitType)) < .01
|| ++noHitCount > 99)
return hitNorm =
!Vec(QueryDatabase(hitPos + Vec(.01, 0), noHitCount) - d,
QueryDatabase(hitPos + Vec(0, .01), noHitCount) - d,
QueryDatabase(hitPos + Vec(0, 0, .01), noHitCount) - d)
, hitType; // Weird return statement where a variable is also updated.return0;
}The idea of ray marching with regard to the distance is to move forward to the distance to the nearest object. In the end, the beam will approach the surface so much that the point can be considered the point of falling.
Note that ray marching does not return a true intersection with a surface, but an approximation. That is why the code marching stops when d <0.01f.
Putting it all together: sampling
Trace path research is almost complete. We lack a bridge that would connect the main () function to ray marcher. This last part, which I renamed “Trace”, is the “brain” in which the rays are reflected or stopped, depending on what they encounter.
Vec Trace(Vec origin, Vec direction){
Vec sampledPosition, normal, color, attenuation = 1;
Vec lightDirection(!Vec(.6, .6, 1)); // Directional lightfor (int bounceCount = 3; bounceCount--;) {
int hitType = RayMarching(origin, direction, sampledPosition, normal);
if (hitType == HIT_NONE) break; // No hit. This is over, return color.if (hitType == HIT_LETTER) { // Specular bounce on a letter. No color acc.
direction = direction + normal * ( normal % direction * -2);
origin = sampledPosition + direction * 0.1;
attenuation = attenuation * 0.2; // Attenuation via distance traveled.
}
if (hitType == HIT_WALL) { // Wall hit uses color yellow?float incidence = normal % lightDirection;
float p = 6.283185 * randomVal();
float c = randomVal();
float s = sqrtf(1 - c);
float g = normal.z < 0 ? -1 : 1;
float u = -1 / (g + normal.z);
float v = normal.x * normal.y * u;
direction = Vec(v,
g + normal.y * normal.y * u,
-normal.y) * (cosf(p) * s)
+
Vec(1 + g * normal.x * normal.x * u,
g * v,
-g * normal.x) * (sinf(p) * s) + normal * sqrtf(c);
origin = sampledPosition + direction * .1;
attenuation = attenuation * 0.2;
if (incidence > 0 &&
RayMarching(sampledPosition + normal * .1,
lightDirection,
sampledPosition,
normal) == HIT_SUN)
color = color + attenuation * Vec(500, 400, 100) * incidence;
}
if (hitType == HIT_SUN) { //
color = color + attenuation * Vec(50, 80, 100); break; // Sun Color
}
}
return color;
}I experimented a bit with this feature to change the maximum number of permissible reflections of the beam. The value “2” gives the letters a surprisingly beautiful Vantablack lacquer paint [4] .
one
2
3
four
Fully cleaned source code
To put everything together, I created a completely cleaned source code .
Reference materials
[1] Source: Post on Twitter lexfrench October 8, 2018
[2] Source: Wikipedia: NetPBM image format
[3] Source: Rendered on the most powerful MacBook Pro, 2017
[4] Source: Wikipedia: Vantab