Time to jump right back into aek’s pathtracer. What’s “5O5_” “5W9W” “5_9_" and how do we get the 3D letter “P” in PIXAR from that? Now we’re starting to get somewhere in our quest to understand every line of code, or rather every symbol printed on that postcard.


This is the second in a series of posts:

— Overview —

In the last post I introduced Andrew Kensler’s postcard-sized pathtracer and the deciphered version by Fabien Sanglard that I’m basing this series on. I still recommend that you read his post first, it covers a lot of the basics of how the program works! And at least skim over my first post, where I went over the process of replacing Andrew’s minimalist Vec class with Urho3D’s Vector3.

— Cleanup, Continued —

Function: QueryDatabase

So now we can dig deeper and take a closer look at the inner workings of our functions, starting with QueryDatabase():

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

The first thing that we’ll unpack are those PIXAR letter coordinates encoded as… letters. Or more accurately, ASCII symbols.

In an earlier version of this I had turned that char array into an Urho String object, because I didn’t understand a warning the compiler was giving me.

It came from a subtle change: Fabien had replaced the hard-coded loop condition of i < 60 with i < sizeof(letters).

We’re iterating in steps of 4, so our i will be 0, 4, 8, …, 56; where we’ll look at the letter at index i and the three after it. Our last bunch should thus be indices 56 to 59. Which is what we get with the original end condition of i < 60, where for the next loop at i = 56+4 we jump right past it and out of the loop we go. That should work out either way, hard-coded value or not, no?

Well, there’s a subtle catch in the substitution: C’s string literals are zero-terminated, and Fabien even made that visible in the declaration letters[15*4+1]. So sizeof(letters) gives us 61, and off we go pretending letters[60...63] exist where they don’t.

To fix that my adjusted loop end condition has become i+3 < sizeof(letters) which will correctly tell us that 63 would be well past the end of our string. (Or for an even easier fix, using C’s strlen() in lieu of sizeof() would’ve given us the “60” we were looking to replace!)

Unpack those coordinates

Alright, so what kind of values can we expect to get out of this anyway? Well, if we look at an ASCII chart for a moment, we’ll find that decimal code 79 is the letter O. Since this is the “baseline” value we subtract from all our chars, the big O somewhat conveniently maps to zero! What is the lowest and highest we’re looking at otherwise? That’d be the character “5” at 53 and the “i” at 105, so they represent coordinate values 53 - 79 = -26 and 105 - 79 = 26 respectively.

But in reality we’re using somewhat different coordinate values, because right after the x/y-pairs end up in Vector2s a *.5 follows right after, i.e. dividing them in half. Is there any reason to store our coordinates this way? In ASCII everything from decimal 32 up to 126 is a printable character. So with the letter O as our center we have 32 - 79 = -57 as our lowest and 126 - 79 = 47 as our highest usable coordinate. That range from -57 to 47 is just about twice as much as we’re using, and in reality we only need to encode values between -13 and 13!

So let’s do that instead:

1
2
5O5_5W9W5_9_COC_AOEOA_E_IOQ_I_QOUOY_Y_]OWW[WaOa_aWeWa_e_cWiO
BOBWBSDSBWDWHOJOIOIWHWJWLOPWLWPOROTWTWVOSSUSXOXWXSZSXWZWYS\O

The top is our original string. The bottom is everything scaled down by that factor of 1/2, still centered at O (notice how it sits at the same spots in both strings). We could almost have everything covered within the alphabet, since we’re using a range across 27 values.

Here’s a few more possible representations:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
>K>S>O@O>S@SDKFKEKESDSFSHKLSHSLKNKPSPSRKOOQOTKTSTOVOTSVSUOXK
?L?T?PAP?TATELGLFLFTETGTILMTITMLOLQTQTSLPPRPULUTUPWPUTWTVPYL
@M@U@QBQ@UBUFMHMGMGUFUHUJMNUJUNMPMRURUTMQQSQVMVUVQXQVUXUWQZM
ANAVARCRAVCVGNINHNHVGVIVKNOVKVONQNSVSVUNRRTRWNWVWRYRWVYVXR[N
BOBWBSDSBWDWHOJOIOIWHWJWLOPWLWPOROTWTWVOSSUSXOXWXSZSXWZWYS\O
CPCXCTETCXEXIPKPJPJXIXKXMPQXMXQPSPUXUXWPTTVTYPYXYT[TYX[XZT]P
^k^s^o`o^s`sdkfkekesdsfshklshslknkpspsrkooqotktstovotsvsuoxk
_l_t_pap_tatelglflftetgtilmtitmlolqtqtslpprpulutupwputwtvpyl
`m`u`qbq`ubufmhmgmgufuhujmnujunmpmrurutmqqsqvmvuvqxqvuxuwqzm
anavarcravcvgninhnhvgvivknovkvonqnsvsvunrrtrwnwvwryrwvyvxr{n
bobwbsdsbwdwhojoioiwhwjwlopwlwporotwtwvossusxoxwxszsxwzwys|o
cpcxctetcxexipkpjpjxixkxmpqxmxqpspuxuxwpttvtypyxyt{tyx{xzt}p

I do kinda like the variants I shall now nickname “Ana” and “Bob”. I never noticed “Bob” is just Ana++. Learn something new every day, eh?

I’m kind of assuming the decision to divide everything in half came after Andrew settled on the nicest of the double range representations. Using “Ana” or “Bob” in the minified original and skipping the division could’ve saved a whole 6 characters. Or maybe the code evolved this way through trial & error on what looks good.

But hey, how about this alternate representation for the coordinates:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
std::vector<int> pixar_points = {
   -13, 0,  -13, 8,         // P stem
   -13, 4,  -11, 4,         // P top bar
   -13, 8,  -11, 8,         // P mid bar

    -7, 0,   -5, 0,         // I bottom bar
    -6, 0,   -6, 8,         // I stem
    -7, 8,   -5, 8,         // I top bar

    -3, 0,    1, 8,         // X slash
    -3, 8,    1, 0,         // X backslash

     3, 0,    5, 8,         // A slash
     5, 8,    7, 0,         // A backslash
     4, 4,    6, 4,         // A bar

     9, 0,    9, 8,         // R stem
     9, 4,   11, 4,         // R mid bar
     9, 8,   11, 8,         // R top bar
    10, 4,   13, 0          // R leg
  };

And here’s what they look like graphed out:

2d graph showing the lines that form the PIXAR letters

Graphed out with python using matplotlib on repl.it

Putting the readable coordinates in place is in this commit.

Distance querying - part 1: Letter lines

So let’s focus on that first for-loop in QueryDatabase for now:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
float distance = 1e9;
Vector2 f(position.x_, position.y_); // Flattened position (z=0)

for (int i = 0; i+3 < pixar_points.size(); i += 4) {
    Vector2 begin = Vector2(pixar_points[i], pixar_points[i + 1]);
    Vector2 e = Vector2(pixar_points[i + 2], pixar_points[i + 3]) - 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.

Even if we don’t know what this does, that line for Vector o is just dying to be pulled apart to help read through what’s going on. Here’s a first step in that direction:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
float distance = 1e9;
Vector2 pos_in_z0(position.x_, position.y_); // Flattened position (z=0)

for (int i = 0; i+3 < pixar_points.size(); i += 4) {
   Vector2 line_begin = Vector2(pixar_points[i], pixar_points[i + 1]);
   Vector2 line_end = Vector2(pixar_points[i + 2], pixar_points[i + 3]);
   Vector2 line_vector = line_end - line_begin;

   float scale_by;
   scale_by = Vector2(line_begin - pos_in_z0).DotProduct(line_vector/line_vector.LengthSquared());
   scale_by = std::min(-std::min(scale_by,  0.0f), 1.0f);

   Vector2 letter_distance = pos_in_z0 - (line_begin + line_vector * scale_factor);

   distance = std::min(distance, letter_distance.LengthSquared()); // compare squared distance.
}
distance = sqrtf(distance); // Get real distance, not square distance.

Let’s take a closer look at this. We have a distance variable that starts out well beyond anything we’ll encounter, because we’ll compare this with whatever distances we find and always keep the smaller one. Right now our focus is on the letters, which are in the z=0 plane — so we project our queried position down onto that and work in 2d here with what I’ve now named pos_in_z0.

Above we had Vector2 o = f - (begin + e * min()) with some factor that we’re scaling our line by, so I pulled that into scale_by; first the bit inside the nested min() calls (line 10) and then applying those min calls (line 11).

What caught my eye in particular in scale_by is that peculiar bit of line_vector / line_vector.LengthSquared()… if we divide by Length we normalize, but why divide by the length squared? Perhaps this “second” division by Length belongs somewhere else?

If either of the vectors in the dot product is multiplied with a scalar, we can pull that out of the dot product (it’s “associative”, but only when it comes to scalar multiplication!) without affecting the result. Since we’re multiplying line_vector by a factor of 1 / line_vector.LengthSquared() or (1 / line_vector.Length())^2, let’s pull one of these lengths further outwards. See lines 4 & 5 in the upcoming code block for this.

It’s also worth noting that our distance comparison involves LengthSquared() instead of Length() and we aren’t taking the square root until after the loop. We know that the Length and LengthSquared of a vector will never be negative. Thus we can assert that sqrt(a) < sqrt(b) <=> a < b holds true for a >= 0 and b >= 0. So we could just as well use letter_distance.Length() in place of letter_distance.LengthSquared() in the last line of the loop, and skip the sqrt() of the very last line here, although that might theoretically sacrifice speed for clarity.

$$\forall a,b\in\N_{0}: \sqrt{a} < \sqrt{b} \Leftrightarrow a < b$$

If we do all this we end up at:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
for (...) {
   ...

   scale_by = Vector2(line_begin - pos_in_z0).DotProduct(line_vector.Normalized());
   scale_by = scale_by / line_vector.Length();
   scale_by = std::min(-std::min(scale_by,  0.0f), 1.0f);

   Vector2 letter_distance = pos_in_z0 - (line_begin + line_vector * scale_by);

   distance = std::min(distance, letter_distance.Length()); // compare distance.
}

// distance = sqrtf(distance);

Now, to make this easier to reason about, here’s a simple substitution that’ll help us, but took me days of brain-wracking to notice:

1
-min(-a,  -b) == max(a,  b)

Yeah. It’s obvious when spelled out like this: If you negate a and b, the less/greater comparisons trade places so to speak, turning a min() into a max(), except we negate that as well to cancel out negating the arguments.

And all this just because Andrew defined his min() and never used a max() function! If we want to make use of this we’ll have to negate scale_by at some point before we plunk it into what shall in a moment be max(). We can pull that factor of -1 all the way “up” into our dot product at the top, and line_begin - pos_in_z0 becomes pos_in_z0 - line_begin. This, alongside the fact that we’re comparing to 0.0f in this instance (which isn’t bothered by negations at all) made this little adjustment camouflage from my view.

1
2
3
scale_by = Vector2(pos_in_z0 - line_begin).DotProduct(line_vector.Normalized());
scale_by = scale_by / line_vector.Length();
scale_by = std::min(std::max(scale_by,  0.0f), 1.0f);

Alright! I promise we’re getting somewhere. Our left-hand side of the dot product is now a vector that points from line_begin to pos_in_z0. We then project it onto our line: The dot product with the normalized vector gives us the factor to scale said normal by to turn it into the projected vector, i.e. the one that points to where we’re at a 90° angle to our position on our line. If that right angle happens to occur at the start of our line, the factor is zero. If it’s at the end of our line, the factor is the length of our line. Anywhere outside of that, we’re also outside this range of factors.

Suddenly the second normalization and the min/max are starting to make sense: After the division 0.0f <= scale_by <= 1.0f means we’re falling within the “top” and “bottom” of our line. And if scale_by falls below 0.0f, the max() will clamp it; if it’s above 1.0f the min() will clamp it. So in the end, scale_by tells us at which point along our line we’re closest to pos_in_z0; but not on an infinite line - only within the bounds of the line segment we’re looking at.

And exactly that’s what we’re calculating here: The vector letter_distance is the smallest distance from our position to our line.

If it helps you visualize what I explained above, we can also write it like this:

1
2
3
4
5
6
7
8
9
scale_by = Vector2(pos_in_z0 - line_begin).DotProduct(line_vector.Normalized());

if(scale_by < 0.0f)
   scale_by = 0.0f;

if(scale_by > line_vector.Length())
   scale_by = line_vector.Length();

Vector2 letter_distance = pos_in_z0 - (line_begin + line_vector.Normalized() * scale_by);

See this and this commit for these steps.

Distance querying - part 2: Curves

Next up: The curves of the P and R.

Once we’ve recovered from all the mathsing and reasoning about (which took me at least a week or two when I first unraveled it) we’re left staring at this piece of code that comes right after:

1
2
3
4
5
6
7
8
// 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 = pos_in_z0 - curves[i];
    distance = std::min(distance,
             o.x_ > 0 ? fabsf(o.Length() - 2) : (o.y_ += o.y_ > 0 ? -2 : 2, o.Length())
                     );
}

Okay, we get two points, one for each curve. If we look them up in our graph of the letters we find them right between the right ends of those little arms where the P and R expect their curve to attach. We have a for-loop checking our distance to each, still checking if they’re closer than what distance is holding so far. But what in the world is happening in that “single” line of code, notably in line 6 in our snippet?

Let’s put the distance = min(...) on its own line again, after calculating whatever we’re comparing to. And it looks like our calculation is split into two cases by a ternary operator based around o.x_ > 0. I’ve gone ahead and turned this into an if-statement and renamed Vector2 o to something slightly more descriptive:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
for (...) {
Vector2 center_to_pos = pos_in_z0 - curves[i];

float comparison_dist;
if (center_to_pos.x_ > 0) {
   comparison_dist = fabsf(center_to_pos.Length() - 2);
} else {
   comparison_dist = (center_to_pos.y_ += center_to_pos.y_ > 0 ? -2 : 2, center_to_pos.Length());
}

distance = std::min(distance, comparison_dist);
}

Okay, let’s deal with the first case first. We made a Vector2 that points from the point we’re checking against towards our position in the plane. If the x coordinate of this vector is > 0 that means it points to the right. So our position is to the right of the center of the curve. For all these points on the right we subtract 2 from the distance to our curve center, and take the absolute value of that. This means if we’re exactly 2 away from that center, we consider ourselves a direct hit, yielding a distance of zero. Anything closer to or further away from the center returns its distance from this radius of 2. Which is perfect, since what we wanted to know here is our point’s distance from a circle with radius 2, but only the right half of that circle.

Now, looking at line 8 above: What are we doing with points that aren’t to the right of our center point? This bit is still rather convoluted, and this is the first time I learned about the comma operator in C/C++. Basically it looks like (<statement>,<statement>) and merely evaluates the first statement and returns the value of the second. So with that in mind we can unpack that and the ternary operator it contains:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
Vector2 center_to_pos = pos_in_z0 - curves[i];

float comparison_dist;
if(center_to_pos.x_ > 0) {
    comparison_dist = std::abs(center_to_pos.Length() - 2.0f);
} else {
    if(center_to_pos.y_ > 0) {
        center_to_pos.y_ -= 2;
    } else {
        center_to_pos.y_ += 2;
    }

    comparison_dist = center_to_pos.Length();
}

distance = std::min(distance, comparison_dist);

(I also replaced the call to fabsf() with std::abs to be more in line with the other uses of the std library)

Look at it this way: Where before we had center_to_pos.y_ = pos_in_z0.y_ - curves[i].y_ now we have center_to_pos.y_ = pos_in_z0.y_ - (curves[i].y_ + 2) if we’re above the point or ... (... - 2) if we’re below/at its height. Or expressed differently yet again: If we’re above (and to the left) we check our distance against the point 2 above the center, otherwise 2 below; these are just the end points of our curve. We actually don’t need this case, because those points are already covered by our straight line segments. The original postcard code could save another 27 characters if we replace the entire second-half return of the outer ternary with a single d for distance.

So with that in mind our thinned out version of this whole section looks like this:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
// 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 center_to_pos = pos_in_z0 - curves[i];

   if(center_to_pos.x_ > 0) {
       float comparison_dist = std::abs(center_to_pos.Length() - 2.0f);
       distance = std::min(distance, comparison_dist);
   }
}

Not too shabby.1

Corresponding commits for this are here and here.

Distance querying - part 3: Squircles

At this point there’s only one step remaining, at least when it comes to letter-distance querying:

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

…Head back into the third dimension! But wait, why are we using \(\cdots^{8}\) and \(\sqrt[8]{\cdots}\)? If we’re gonna take into account our z-axis component, shouldn’t our calculation be 3d distance = sqrt(2d_distance^2 + z^2)? Well, yes and no. If you change the numbers here and run the code, you will see the following:

A render of the PIXAR letters with round edges in all dimensions.

Our letters are perfectly nice tubes with a radius of 0.5f. This bit of distance calculation is very similar to what we saw for the curve, subtracting the radius from our distance. Except in this case we intentionally do not want to remove the sign of the result, to keep it a signed distance: Returning a positive number outside, negative number inside of our letters.

But back the roots and powers at play. What is being used here could also be called a “squircle”.

Preview picture of geogebra link.

Here’s a little interactive squircle visualization I made on geogebra.org. Try setting p=1 or p=0.5! (Warning: Probably only suitable for desktop browser and even there it can lag out quite a bit.)

Close-up of letter edges & line ends.

Close-up of letter edges & line ends.

Also worth noting is the fact that the squircle cross-sections only apply for adding the z-axis back in. The ends of our lines as we calculated them in 2d still use regular euclidean distance. Replacing that with Minkowski distance of p=8 gives all around blockier, interesting letters too, but the extra powf() calls drive up render time on our poor single-threaded implementation. Turning our squircles back to circles (as in the render above) does also give a slight performance boost.


That’s the last bit of code concerning the letters themselves, so it’s time to wrap up this post. Here’s how far we got on revising QueryDatabase()

 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
std::vector<int> pixar_points = {
         -13, 0,  -13, 8, // P stem
         -13, 4,  -11, 4, // P top bar
         -13, 8,  -11, 8, // P mid bar

          -7, 0,   -5, 0, // I bottom bar
          -6, 0,   -6, 8, // I stem
          -7, 8,   -5, 8, // I top bar

          -3, 0,    1, 8, // X slash
          -3, 8,    1, 0, // X backslash

           3, 0,    5, 8, // A slash
           5, 8,    7, 0, // A backslash
           4, 4,    6, 4, // A bar

           9, 0,    9, 8, // R stem
           9, 4,   11, 4, // R mid bar
           9, 8,   11, 8, // R top bar
          10, 4,   13, 0  // R leg
        };

// Sample the world using Signed Distance Fields.
float QueryDatabase(Vector3 position, int &hitType) {
    float distance = 1e9;
    Vector2 pos_in_z0(position.x_, position.y_); // Flattened position (z=0)

    for (int i = 0; i+3 < pixar_points.size(); i += 4) {
        Vector2 line_begin = Vector2(pixar_points[i], pixar_points[i + 1]);
        Vector2 line_end = Vector2(pixar_points[i + 2], pixar_points[i + 3]);
        Vector2 line_vector = line_end - line_begin;

        // factor to scale line_vector by to project position onto line
        // i.e. with this factor alone (line_begin + line_vector.Normalized() * scale_factor)
        // gives us the closest point on an infinite line
        float scale_factor = Vector2(pos_in_z0 - line_begin).DotProduct(line_vector.Normalized());

        // But now we clamp the scaling within 0.0f <= scale_factor <= line length
        scale_factor = std::min(std::max(scale_factor,  0.0f), line_vector.Length());

        Vector2 letter_distance = pos_in_z0 - (line_begin + line_vector.Normalized() * scale_factor);

        distance = std::min(distance, letter_distance.Length()); // compare 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 center_to_pos = pos_in_z0 - curves[i];

        if(center_to_pos.x_ > 0) {
            float comparison_dist = std::abs(center_to_pos.Length() - 2.0f);
            distance = std::min(distance, comparison_dist);
        }
    }

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

There’s still some classic CSG calculations at the end that we haven’t covered yet, plus the quick check for hitting our “sun” light source.

As part of the next post the BoxTest() function will become its own class to pair it up with the hard-coded vectors in its current function calls. And hopefully I will cover RayMarching() too, as I promised somewhere along the way. And part 4 for all the juicy vector reflection maths. Stay tuned!

— Interim result —

And once again for completeness' sake, here’s the entire body of code we got so far, minus the Urho3D scaffolding:

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


std::vector<int> pixar_points = {
         -13, 0,  -13, 8, // P stem
         -13, 4,  -11, 4, // P top bar
         -13, 8,  -11, 8, // P mid bar

          -7, 0,   -5, 0, // I bottom bar
          -6, 0,   -6, 8, // I stem
          -7, 8,   -5, 8, // I top bar

          -3, 0,    1, 8, // X slash
          -3, 8,    1, 0, // X backslash

           3, 0,    5, 8, // A slash
           5, 8,    7, 0, // A backslash
           4, 4,    6, 4, // A bar

           9, 0,    9, 8, // R stem
           9, 4,   11, 4, // R mid bar
           9, 8,   11, 8, // R top bar
          10, 4,   13, 0  // R leg
        };

// Sample the world using Signed Distance Fields.
float QueryDatabase(Vector3 position, int &hitType) {
    float distance = 1e9;
    Vector2 pos_in_z0(position.x_, position.y_); // Flattened position (z=0)

    for (int i = 0; i+3 < pixar_points.size(); i += 4) {
        Vector2 line_begin = Vector2(pixar_points[i], pixar_points[i + 1]);
        Vector2 line_end = Vector2(pixar_points[i + 2], pixar_points[i + 3]);
        Vector2 line_vector = line_end - line_begin;

        // factor to scale line_vector by to project position onto line
        // i.e. with this factor alone (line_begin + line_vector.Normalized() * scale_factor)
        // gives us the closest point on an infinite line
        float scale_factor = Vector2(pos_in_z0 - line_begin).DotProduct(line_vector.Normalized());

        // But now we clamp the scaling within 0.0f <= scale_factor <= line length
        scale_factor = std::min(std::max(scale_factor,  0.0f), line_vector.Length());

        Vector2 letter_distance = pos_in_z0 - (line_begin + line_vector.Normalized() * scale_factor);

        distance = std::min(distance, letter_distance.Length()); // compare 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 center_to_pos = pos_in_z0 - curves[i];

        if(center_to_pos.x_ > 0) {
            float comparison_dist = std::abs(center_to_pos.Length() - 2.0f);
            distance = std::min(distance, comparison_dist);
        }
    }

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

And here’s the relevant part of main() (unchanged from last post):

421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
   //  int w = 960, h = 540, samplesCount = 8;
   Vector3 position(-22, 5, 25);
   Vector3 goal = Vector3(-3, 4, 0) - position;
   goal.Normalize();

   // this vector is actually to the right of goal... flipped signs of x and z for left.
   // Up is down. Hence the subtractions instead of additions for trace directions
   // Note to self: Yes, it was rendering upside-down, but only because originally it
   // wrote values for y = h .. 0 and x = w .. 0, so from bottom right to top left, sequentially
   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)


Thank you again for reading so far. I hope you found this post interesting now that we’re getting into the really fun bits of code. And thanks again to Andrew and Fabien!

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


  1. (I kind of wish we could also use LengthSquared in the remaining first half of the conditional, since then we could delay the sqrt() we did in the original code and save those 6 characters in this bit; but that’d break our nice curves; or at least I haven’t found any way to transform it to get the sqrt() out of it and I strongly suspect it’s mathematically impossible.) ↩︎