Postcard-sized ray tracer decoding

“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 DOOM , I had time to study the insides of its mysterious code. And almost immediately I was literally fascinated by the techniques found in it. They were very different from the previous work of Andrew, based on the "standard" ray tracer. I was interested in learning about ray marching, the functions of constructive volumetric geometry, Monte-Carlo rendering / tracing, and the 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 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

Does 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 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.

Very much can be learned from the image shown above. 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 up 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 section is described in detail in the remainder of this article:
■ - common tricks, ■ - Vector class, ■ - auxiliary code, ■ - database, ■ - Ray marching, ■ - sampling, ■ - 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 the 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).

 struct Vec { 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 + rx, y + ry, z + rz); } Vec operator*(Vec r) { return Vec(x * rx, y * ry, z * rz); } // dot product float operator%(Vec r) { return x * rx + y * ry + z * rz; } // 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. It is usually 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.

 int main() { 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 vector Vec 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 / ox, color.y / oy, color.z / oz) * 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 can be assigned the value 1 to speed up rendering and iterations. Here are examples of visualizations with sample values from 1 to 2048.

Spoiler header

one

four

eight

sixteen

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.

 float min(float l, float r) { return l < r ? l : r; } float randomVal() { return (float) rand() / RAND_MAX; } // Rectangle CSG equation. Returns minimum signed distance from // space carved by lowerLeft vertex and opposite rectangle // vertex upperRight. float BoxTest(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 it is enough just to 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 for example the sphere and two points, A and B.

 // Signed distance point(p) to sphere(c,r) float testSphere(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 as a function of BoxTest).

  // Signed distance point(p) to Box(c1,c2) float testRectangle(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) float testCarveBox(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 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 float testRoom(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. float QueryDatabase(Vec position, int &hitType) { float distance = 1e9; Vec f = position; // Flattened position (z=0) fz = 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, ox > 0 ? fabsf(sqrtf(o % o) - 2) : (oy += oy > 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 construct 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. Ray marching uses the distance function. 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/normal int RayMarching(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 marching for (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. return 0; }

The idea of ray marching taking into account 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 light for (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] .