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 finally 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). But here’s the catch: 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 56 - 59. Which is what we get with the original end condition of < 60, where from 56+4 we jump right past it and out of the loop we go.

But maybe you found the catch in the substitution already: C’s string literals are zero-terminated, and fabien even made that explicit by declaring 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 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 useable 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!

 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 everthing 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^oo^ssdkfkekesdsfshklshslknkpspsrkooqotktstovotsvsuoxk _l_t_pap_tatelglflftetgtilmtitmlolqtqtslpprpulutupwputwtvpyl muqbqubufmhmgmgufuhujmnujunmpmrurutmqqsqvmvuvqxqvuxuwqzm 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.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21  std::vector 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:

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 (l. 10) and then applying those min calls (l. 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 ^8 and ^(1/8)? 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:

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

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

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 themself, 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 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 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.) ↩︎