Right around the start of 2020, I came across the blog of Fabien Sanglard (fabiensanglard.net). I believe it was because his current series of blog posts (at the time of writing/stumbling-upon) about Another World and its ports were hitting the front page on HN.

But at the time a different post caught my eye and led me down a deep rabbit hole lined with copious amounts of vector math… it’s the one about “Deciphering The Postcard Sized Raytracer”.

While writing this it became clear that I have to split it up into multiple posts:

• The Postcard Pathtracer, Part 1: Deciphering, Continued (you are here)
• The Postcard Pathtracer, Part 2: Of Letters, Coords & Squircles
• The Postcard Pathtracer, Part 3: Of Boxes, Lines & Rays
• The Postcard Pathtracer, Part 4: Hemispheres on vector bases

Here’s the pièce de résistance:

Yep, that’s a full-fledged pathtracer on a postcard, written by Andrew Kensler. Here’s what running it yields:

Fabien’s post is a great stab at making sense of this wholly opaque, inscrutable mess and makes its constituent parts quite a bit more digestible. This was my first time seeing the inner workings of a path tracer, or any ray tracer for that matter, and the first time I heard about CSG (Constructive Solid Geometry) and SDF (Signed Distance Functions/Fields). I also had to take an extensive refresher on vector math and trigonometric functions to get everything figured out and tidied up without breaking anything.

The two big things I really wanted to figure out are:

• how do the letters get the shape that they have, especially with the rounded off edges?
• what in the world are those calculations for the diffuse reflection?

I highly suggest that you read Fabien’s post first, as this one will be taking his findings and deobfuscated code as its basis. And who knows, maybe you’ll be down the rabbit hole alongside me before you know it?

Once you’re ready, let’s dig in further!

— Overview —

For reference, here’s the version of the program we’re starting out with, as Fabien left it:

(And if you have javascript enabled it should hopefully fold itself and the highlighted functions neatly. I hope.)

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196

#include <stdlib.h> // card > pixar.ppm
#include <stdio.h>
#include <math.h>

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 + r.x, y + r.y, z + r.z); }

    Vec operator*(Vec r) { return Vec(x * r.x, y * r.y, z * r.z); }

    float operator%(Vec r) { return x * r.x + y * r.y + z * r.z; }

    // intnv square root
    Vec operator!() {
      return *this * (1 / sqrtf(*this % *this)
      );
    }
};

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));
}

#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)
  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;
}

// 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;
}

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;
}

int main() {
//  int w = 960, h = 540, samplesCount = 16;
  int w = 960, h = 540, samplesCount = 8;
  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 / o.x, color.y / o.y, color.z / o.z) * 255;
      printf("%c%c%c", (int) color.x, (int) color.y, (int) color.z);
    }
}// Andrew Kensler

Looking at the basic structure of the program, we have the following four functions:

• main() – sets up camera vectors, and traces samples for each pixel, averaging them together. Outputs a NetPBM image.
• Trace() – Traces a single ray for a set amount of reflections off objects — specular (letters) or diffuse (walls) — or until it hits the light source. Returns a color value.
• RayMarching() – Performs what Fabien calls “signed sphere marching” or “distance assisted ray marching”: Gets distance to closest object, and moves the ray by that distance. Repeat until close enough to an object, or exit with no hit after set amount of iterations; returns hit position, normal vector and type of surface.
• QueryDatabase() – Calculates distance to either the PIXAR letters (more on how that works later), to the walls and other surfaces of the room, and to the “sun”; returns smallest distance, and type of object hit. Oh, and the “hit position” of that closest point.

I’m using the function names assigned by Mr. Sanglard, listed in reverse, which is basically also the hierarchy they’re called in, each using the output from calls to the following. Think of the list above as a top-down approach, rather than bottom-up.

Besides that there’s also some helper functions and utilities:

• BoxTest() – Helps get the signed distance from a point to the walls of a box
• randomVal() – Returns a pseudo-random float between 0.0f..1.0f
• min() – a min function for floats
• struct Vec – And lastly, a class for a 3d vector of floats

— Plan of action —

This series will roughly consist of the following sections:

Part 1

Methodically switching the Vec class out and making references to that more readable. I’m afraid that’s all I got covered in this first post.

Part 2

Go over QueryDatabase() as mentioned above, which covers the math involved in distance checks, which in turn give the letters their distinct shapes!

Part 3

Quickly go over RayMarching() while adding helpers and variables that give us more flexibility concerning the camera & light positions.

Part 4

Finally addressing the diffuse reflections and everything that happens in Trace().

And Beyond

Tone-mapping, tracing & rendering optimizations, frame-by-frame animations, …?

Whenever I’m addressing specific code changes that I applied, I will be referencing commits in a cleaned up git repository I set up at:

https://git.redflames.space/rf/aek_postcard_pathtracer/commits/branch/master

Before doing anything else I converted the code to an Urho3D project:

1. For one to make it easier to render the output to an image and display it as a sprite;
2. Later on as quality-of-life improvements I added some UI to tweak rendering parameters, the option to iterate over increasing sample sizes, and periodically update the output while rendering.
3. Oh, and to save to PNG with a single function call. :)

For the sake of clarity I’ll omit all that scaffolding code, and focus solely on the original pathtracer code snippets. I will include links to the relevant commits, where applicable.

(Sidenote: I hadn’t really used Urho3D before this, so I based my code off some sample projects. I built Urho from source with MinGW64 on Windows, and set up a CodeBlocks project.)

— Cleanup duty —

First overhaul: Vec class

Let’s start by looking at the Vec class in Fabien’s version:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

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 + r.x, y + r.y, z + r.z); }

    Vec operator*(Vec r) { return Vec(x * r.x, y * r.y, z * r.z); }

    float operator%(Vec r) { return x * r.x + y * r.y + z * r.z; }

    // intnv square root
    Vec operator!() {
      return *this * (1 / sqrtf(*this % *this)
      );
    }
};

It’s a vector composed of three float values. The first constructor takes a single float to initialize all three components to that same value. This one is important, because it is used to implicitly convert float to Vec (and int too, of course!) where necessary. The second constructor takes two or three float arguments, as the last component may default to 0.

Then we have overloaded operators for addition and multiplication with another Vec. But because of the implicit cast mentioned above, as we’ll later see in the rest of the code, these operators also get used when the left-hand side is a Vec and the right-hand side is an int or float. These two are also used when subtraction and division are needed, since you can just do + (x * -1) and * (1 / x) for those respectively.

The other two operators overloading % and ! are the dot product and the normalized version of the Vec: The comment there isn’t entirely accurate since the dot product of a vector with itself is x*x + y*y + z*z — the square root of which is its length!

And then, dividing a vector by its length gives us a the normalized vector, i.e. same direction but length 1.

What’s notably absent here is the cross-product: It is only used once in the main function to construct the normal Vec up. The code there is

V 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);

compared to adding an overloaded operator and calling that instead:

V O&(V r){R V(y*r.z-z*r.y,z*r.x-x*r.z,x*r.y-y*r.x);} … V u=g&l;

…which does use up 7 precious extra characters and is thus obviously rather unfit for the purpose of squeezing a ray tracer onto a postcard. :)¹

Now, for immediate comparison, here’s my first overhaul of the Vec class:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

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 rhs) { return Vec(x + rhs.x, y + rhs.y, z + rhs.z); }
    Vec operator-(Vec rhs) { return Vec(x - rhs.x, y - rhs.y, z - rhs.z); }

    Vec operator*(Vec rhs) { return Vec(x * rhs.x, y * rhs.y, z * rhs.z); }
    Vec operator/(Vec rhs) { return Vec(x / rhs.x, y / rhs.y, z / rhs.z); }

    // dot product
    float operator%(Vec rhs) { return x * rhs.x + y * rhs.y + z * rhs.z; }

    // cross-product
    Vec operator&(Vec rhs) {
        return Vec(y * rhs.z - z * rhs.y,
                   z * rhs.x - x * rhs.z,
                   x * rhs.y - y * rhs.x);
    }

    // actually, this returns a normalized vector
    Vec operator!() {
        return *this / sqrtf(*this % *this);
    }
};

Mirroring the addition we have subtraction, alongside multiplication there’s division, and we’ve gained a handy cross-product, which I chose to overload onto &. The operators for dot- and cross-product don’t help readability very much — quite the opposite — but since my goal in the long run is to replace the whole class with Urho3D’s Vector3, we won’t see them around for much longer. The constructors remain unchanged.

After that I went over the rest of the code and cleaned up instances where subtraction and division were clearly expressed through the other two operators. You can see all the changes in this commit.

I’ll show one instance of cleaner code for comparison. At the start of the original main() we had:

   Vec goal = !(Vec(-3, 4, 0) + position * -1);
   Vec left = !Vec(goal.z, 0, -goal.x) * (1. / w_);

vs.

   Vec goal = !(Vec(-3, 4, 0) - position);
   Vec left = !Vec(goal.z, 0, -goal.x) / w_;

One other change worth mentioning is this:

   color = Vec(color.x / o.x, color.y / o.y, color.z / o.z) * 255;

   color = color / o;

…much nicer, but I did also omit the * 255 since Urho’s Color class expects values from 0.0 to 1.0, and a few lines later I would’ve just done /255.0 anyways.

I am actually a bit surprised the code didn’t have any need for Vec division. We can only really get away with writing division by a scalar via multiplication, since * (1 / v) will not compile without an operator/ (besides the fact that it also will not implicitly cast the int to a Vec on that left-hand side).

Next I threw in these:

   float LengthSquared() {
       return *this % *this;
   }

   float Length() {
       return sqrtf(this->LengthSquared());
   }

As we already noted, the cross-product of a vector with itself is the square of its length. There are a few places where that value actually comes into use by itself, but for the cases where we need the Length() itself, there’s a function for that too now. There aren’t all that many opportunities for substitutions with these in this commit but we’ll be glad for the boost in readability soon enough.

Two other small addition we’ll need for upcoming changes:

   Vec(Vector3 v) {
       x = v.x_;
       y = v.y_;
       z = v.z_;
   }

   Vector3 v3() {
       return Vector3(x, y, z);
   }

Just two helpers in our quest to gradually transition to Urho3D::Vector3.²

I now turned BoxTest’s parameters into Vector3:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

// Rectangle CSG equation. Returns minimum signed distance from
// space carved by lowerLeft vertex and opposite corner vertex upperRight.
float BoxTest(Vector3 position, Vector3 lowerLeft, Vector3 upperRight) {
    Vector3 toLowerLeft = position - lowerLeft;
    Vector3 toUpperRight = upperRight - position;

    return -min( // <- minimum distance out of all, i.e. to closest wall of the box
                min( // <- minimum distance on x- OR y-axis
                    min(toLowerLeft.x_, toUpperRight.x_),  // <- minimum distance on x-axis
                    min(toLowerLeft.y_, toUpperRight.y_)   // <- minimum distance on y-axis
                ),
                min(toLowerLeft.z_, toUpperRight.z_)       // <- minimum distance on z-axis
            );
}

The original code modified the parameters whereas I created two local vectors to make the naming clearer. Fabien already had a solid introduction to the ins and outs of CSG, along with visual explanations using a 2d rectangle analogous to the 3d box, so hopefully that gave you a decent understanding of what this function does. ³

All the call-sites to BoxTest have their arguments altered to construct Vector3s from scratch or create them from existing Vec variables via the v3() function for now. See this commit for the changes.

Yet more substitutions:

• replaced the sole Vec argument (position) of QueryDatabase with a Vector3 but nothing else yet.

• replaced all four arguments of RayMarching with Vector3.

• replaced both arguments to Trace, as well as sampledPosition and normal with Vector3. (Note the different constructors: Vec(1) initializes as (1,1,1) while Vector3 explicitly needs all three arguments!)

• in the rendering code (formerly main()) the second argument to Trace was a complex computation that included the !operator to get a normalized Vec. I split this off into a helper variable, since Vector3’s *.Normalized() call is more verbose. I even remembered to turn the dangerous single-line for-loop into a braced one. Woo!

You can find all the details of these changes in this commit.

And, before we move onto the meat of the actual code and what it does, we have one more round of vector-swapping in this commit.

A couple of notes:

• I took the liberty to convert the vectors used for calculations related to the lines and curves of the letters into Vector2s since they all have their z-component set to zero at all times.

• Split off some !operators into Vector3.Normalize() calls again.

• I made a mistake in the Trace code, where I assumed color starts as ones, but it’s zeroes… when I initially played around with the code I fixed it by adding 3 instead of one to the divisor during tone-mapping… whoops ( https://stackoverflow.com/questions/6838408/how-can-i-declare-and-define-multiple-variables-in-one-line-using-c/6838480 )

Finally we can remove Vec, and I removed the min() function alongside it, replacing it with std::min. Corresponding commit here.

— Interim result —

To wrap up this first post, let’s look at the entire body of code we got so far, minus the Urho3D scaffolding:

 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162

float randomVal() { return (float) rand() / RAND_MAX; }

// Rectangle CSG equation. Returns minimum signed distance from
// space carved by lowerLeft vertex and opposite corner vertex upperRight.
float BoxTest(Vector3 position, Vector3 lowerLeft, Vector3 upperRight) {
    Vector3 toLowerLeft = position - lowerLeft;
    Vector3 toUpperRight = upperRight - position;

    return -std::min( // <- minimum distance out of all, i.e. to closest wall of the box
                std::min( // <- minimum distance on x- OR y-axis
                    std::min(toLowerLeft.x_, toUpperRight.x_),  // <- minimum distance on x-axis
                    std::min(toLowerLeft.y_, toUpperRight.y_)   // <- minimum distance on y-axis
                ),
                std::min(toLowerLeft.z_, toUpperRight.z_)       // <- minimum distance on z-axis
            );
}

#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(Vector3 position, int &hitType) {
    float distance = 1e9;
    Vector2 f(position.x_, position.y_); // Flattened position (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) {
        Vector2 begin = Vector2(letters[i] - 79, letters[i + 1] - 79) * .5;
        Vector2 e = Vector2(letters[i + 2] - 79, letters[i + 3] - 79) * .5 - begin;
        Vector2 o = f - (begin + e * min(-min(Vector2(begin - f).DotProduct(e / e.LengthSquared()),
                                          0.0f),
                                     1.0f)
                    );
        distance = std::min(distance, o.LengthSquared()); // compare squared distance.
    }
    distance = sqrtf(distance); // Get real distance, not square distance.

    // Two curves (for P and R in PixaR) with hard-coded locations.
    Vector2 curves[] = {Vector2(-11, 6), Vector2(11, 6)};
    for (int i = 2; i--;) {
        Vector2 o = f - curves[i];
        distance = std::min(distance,
                       o.x_ > 0 ? fabsf(o.Length() - 2)
                               : (o.y_ += o.y_ > 0 ? -2 : 2, o.Length())
                       );
    }

    distance = powf(powf(distance, 8) + powf(position.z_, 8), .125) - .5;
    hitType = HIT_LETTER;

    float roomDist ;
    roomDist = std::min(// min(A,B) = Union with Constructive solid geometry
                    //-min carves an empty space
                    -std::min(// Lower room
                         BoxTest(position, Vector3(-30, -.5, -30), Vector3(30, 18, 30)),
                         // Upper room
                         BoxTest(position, Vector3(-25, 17, -25), Vector3(25, 20, 25))
                    ),
                    BoxTest( // Ceiling "planks" spaced 8 units apart.
                      Vector3(fmodf(fabsf(position.x_), 8),
                          position.y_,
                          position.z_),
                      Vector3(1.5, 18.5, -25),
                      Vector3(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;
}

// Perform signed sphere marching
// Returns hitType 0, 1, 2, or 3 and update hit position/normal
int RayMarching(Vector3 origin, Vector3 direction, Vector3 &hitPos, Vector3 &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 =
                 Vector3(QueryDatabase(hitPos + Vector3(.01, 0, 0), noHitCount) - d,
                      QueryDatabase(hitPos + Vector3(0, .01, 0), noHitCount) - d,
                      QueryDatabase(hitPos + Vector3(0, 0, .01), noHitCount) - d).Normalized()
             , hitType; // Weird return statement where a variable is also updated.
    return 0;
}

Vector3 Trace(Vector3 origin, Vector3 direction) {
    Vector3 sampledPosition(1, 1, 1);
    Vector3 normal(1, 1, 1);
    Vector3 color(0, 0, 0);
    Vector3 attenuation(1, 1, 1);
    Vector3 lightDirection(.6, .6, 1); // Directional light
    lightDirection.Normalize();

    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.DotProduct(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.DotProduct(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 = Vector3(v,
                          g + normal.y_ * normal.y_ * u,
                          -normal.y_) * (cosf(p) * s)
                      +
                      Vector3(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 * Vector3(500, 400, 100) * incidence;
        }
        if (hitType == HIT_SUN) { //
          color = color + attenuation * Vector3(50, 80, 100); break; // Sun Color
        }
    }
    return color;
}

398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427

   //  int w = 960, h = 540, samplesCount = 8;
   Vector3 position(-22, 5, 25);
   Vector3 goal = Vector3(-3, 4, 0) - position;
   goal.Normalize();
   Vector3 left = Vector3(goal.z_, 0, -goal.x_).Normalized() / w_;

   // Cross-product to get the up vector
   Vector3 up = goal.CrossProduct(left);

   // this used to be the argument to Trace
   Vector3 trace_dir{};

   for (int x = 0; x < w_; x++) {
       Vector3 color;
       for (int p = samplesCount_; p--;) {
           // this was the second Trace argument as both lines in one:
           // !(goal - left * (x - w_ / 2 + randomVal()) - up * (y_ - h_ / 2 + randomVal()))
           trace_dir = goal - left * (x - w_ / 2 + randomVal()) - up * (y_ - h_ / 2 + randomVal());
           trace_dir.Normalize();

           color = color + Trace(position, trace_dir);
       }

       // Reinhard tone mapping
       color = color / samplesCount_ + Vector3(1.0f, 1.0f, 1.0f) * 14. / 241.0;
       Vector3 o = color + Vector3(1.0, 1.0, 1.0);
       color = color / o;

       [...]
   }

(full main.cpp here)

To be continued…

In the next post we will:

• convert that jumble of letters in QueryDatabase() into readable coordinates.
• pick apart the maths involved in checking the distance to the letter lines and how the curve is generated, too.
• explain how the letters get their blocky-but-rounded edges.

If it’s online by the time you read this, I will be linking to Part 2 here, as well as in the list of posts near the start of this one.

Beyond that I am planning future posts covering the remaining code, with topics such as:

• how the rays are reflected once they hit an object, because as I mentioned in the beginning I had quite a hard time working out how the diffuse reflections bounce off the walls.
• how the colors are calculated based off the light source and attenuation.
• and I might cover some basics of tone mapping (think HDR), assuming I actually get a better understanding of it by then.

Thank you for reading so far. I hope you found this post interesting despite how much of a deep-dive into the minutiae of Andrew Kensler’s code it’s turning out to be. Thanks again to him for making this pathtracer and to Fabien for his analysis that kicked off this investigation of my own!

If you have any corrections, constructive criticism or comments, please post them in the comments section below, or write to me via e-mail or Twitter! Feedback highly welcome.

~rf