To start, let's take the final solution from Part 1. The goal is to reduce this even further so that the technique can be used frequently in shaders.

The code is conceptually simple but is a mess of maths and constants screaming out to be simplified. This could be much faster so let's kick a can down the road of simplification and see if it hits anything.

##### Initial Reduction

Take the original point projection onto the plane:

float d = dot(v0 - P, face_dir) / dot(P, face_dir);
vec3 point_on_plane = P + P * d;


The numerator dot product distributes:

float d = (dot(v0, face_dir) - dot(P, face_dir)) / dot(P, face_dir);


and the second dot(P, face_dir) division cancels:

float d = dot(v0, face_dir) / dot(P, face_dir) - 1.0;


Only v0.z is set here so we can immediately multiply that through:

float d = face_dir.z * face_dir.z / dot(P, face_dir) - 1.0;


Face direction components can only be 1 or -1 so their square will always be 1:

float d = 1.0 / dot(P, face_dir) - 1.0;


Now substitute d back into the projection:

vec3 point_on_plane = P + P * (1.0 / dot(P, face_dir) - 1.0);


Distribute P:

vec3 point_on_plane = P + P  / dot(P, face_dir) - P;


And P disappears!

vec3 point_on_plane = P / dot(P, face_dir);


Next up are the three basis vectors:

vec3 O = (v0 + v1) * 0.5;
vec3 X = (v1 - O) * inv_oct_side_len;
vec3 Y = (v2 - O) * inv_oct_tri_height;


There are a bunch of zero components in here that should allow some simplification. O sums as:

vec3 O = vec3(face_dir.x, 0.0, face_dir.z) * 0.5;


X sums as:

vec3 X = vec3(face_dir.x - face_dir.x * 0.5, 0.0, -face_dir.z * 0.5) * inv_oct_side_len;


Given that k - k * 0.5 is always k * 0.5 this simplifies a little:

vec3 X = vec3(face_dir.x, 0.0, -face_dir.z) * inv_oct_side_len * 0.5;


Y sums as:

vec3 Y = vec3(-face_dir.x * 0.5, face_dir.y, -face_dir.z * 0.5) * inv_oct_tri_height;


If we pair the calculation of X with uv.x it's easy to notice something:

vec3 X = vec3(face_dir.x, 0.0, -face_dir.z) * inv_oct_side_len * 0.5;
float sub_tri_edge_len = inv_oct_side_len / TriHalfBasesPerEdge;
uv.x = uv.x / sub_tri_edge_len;


inv_oct_side_len cancels leaving only:

vec3 X = vec3(face_dir.x, 0.0, -face_dir.z);
uv.x = uv.x * TriHalfBasesPerEdge * 0.5;


inv_oct_tri_height unfortunately doesn't cancel but it distributes, allowing the multiply to be deferred to a single scalar later:

vec3 Y = vec3(-face_dir.x * 0.5, face_dir.y, -face_dir.z * 0.5);
uv.y = (uv.y * inv_oct_tri_height * inv_oct_tri_height) * Rows;


Folding all the immediates into shader constants leaves us with:

// Get octant index
ivec3 side;
side.x = P.x >= 0.0 ? 1 : 0;
side.y = P.y >= 0.0 ? 2 : 0;
side.z = P.z >= 0.0 ? 4 : 0;
int octant_index = side.x + side.y + side.z;

// Generate face direction from sidedness
vec3 face_dir = vec3(side) * vec3(2.0, 1.0, 0.5) - 1.0;

// Projection onto face
vec3 point_on_plane = P / dot(P, face_dir);
vec3 X = face_dir * vec3(1.0, 0.0, -1.0);
vec3 Y = face_dir * vec3(-0.5, 1.0, -0.5);

// 2D projection
vec2 uv;
uv.x = dot(point_on_plane, X) + 1.0;
uv.y = dot(point_on_plane, Y) + 0.5;

// Scale for indexing
uv.x = uv.x * Constants.x;
uv.y = uv.y * Constants.y;


and:

Constants.x = 1 << depth;
Constants.y = (1 << depth) / 1.5;


This is looking great! But we're not done yet...

##### Constraining to an Octant

Something feels redundant here: the octant that a point projects to is used twice in the calculation; first to determine and octant index offset and then to project onto the correct face. Conceptually the projection for all faces can be treated as a projection onto one face only, using the octant index to offset and make unique. Let's try it!

Let's project all vertices onto the face with face_dir={1,1,1}. All points that don't project onto this face are simply flipped around using abs:

// Get octant index
int side_x = P.x >= 0.0 ? 1 : 0;
int side_y = P.y >= 0.0 ? 2 : 0;
int side_z = P.z >= 0.0 ? 4 : 0;
int octant_index = side_x + side_y + side.z;

// Projection onto face
vec3 point_on_plane = abs(P) / dot(abs(P), vec3(1));
vec3 X = vec3(1.0, 0.0, -1.0);
vec3 Y = vec3(-0.5, 1.0, -0.5);

// 2D projection
vec2 uv;
uv.x = dot(point_on_plane, X) + 1.0;
uv.y = dot(point_on_plane, Y) + 0.5;

// Scale for indexing
uv.x = uv.x * Constants.x;
uv.y = uv.y * Constants.y;


There's a couple of very interesting geometric operations in here that are worth investigating in more detail.

##### Visualising the Planes

The last two dot products look like plane/distance equations, except the plane normals aren't unit length. They can be normalised by pushing the the required scale factor into the CPU constants. The length of the first normal is sqrt(2) and the second is sqrt(1.5), so:

// Projection onto face and point/plane distances using 4x vectors
vec4 point_on_plane = vec4(abs(P) / dot(abs(P), vec3(1)), 1.0);
vec4 X = vec4(1.0, 0.0, -1.0, 1.0) / sqrt(2.0);
vec4 Y = vec4(-0.5, 1.0, -0.5, 0.5) / sqrt(1.5);
vec2 uv;
uv.x = dot(point_on_plane, X);
uv.y = dot(point_on_plane, Y);


and the CPU constants are then:

k = 1 << depth;
Constants.x = k * Math.sqrt(2);
Constants.y = k / 1.5 * Math.sqrt(1.5);


Rendering the planes gives:

These two planes are clearly orthogonal and have a line of intersection that passes through the bottom left point, as we worked through in Part 1.

The scale constants also make sense: both distances are initially normalised by the meridian length ($$\frac{2}{\sqrt{2}}$$)/face height ($$\sqrt{1.5}$$) and then multiplied by the the number of rows/columns to get a 2D index.

##### Changing the Frame of Reference

It's important to note that this is an orthogonal projection and that as long as the planes stay orthogonal, you can rotate them and re-run the projection with adjusted constants to get a unique 2D position. Index reconstruction will be difficult, however, as the planes will no longer be aligned with the rows and columns...

...except for one special case that is easier to visualise with an orthographic projection:

Go ahead and rotate the view to get a clear picture of what's going on here. Immediately we can see:

• The planes are much simpler: $$<1,0,0,0>$$ and $$<0, 1, 0, 0>$$.
• The two planes clearly align with rows and columns.
• There is no difference of triangle orientation between rows.
• Rows and columns are now scaled equally with no arbitrary scale factors anywhere.

The initial version of this in code looks like:

// Get octant index
int side_x = P.x >= 0.0 ? 1 : 0;
int side_y = P.y >= 0.0 ? 2 : 0;
int side_z = P.z >= 0.0 ? 4 : 0;
int octant_index = side_x + side_y + side_z;

// Projection onto face
vec3 point_on_plane = abs(P) / dot(abs(P), vec3(1));
vec3 X = vec3(1.0, 0.0, 0.0);
vec3 Y = vec3(0.0, 1.0, 0.0);

// 2D projection
vec2 uv;
uv.x = dot(point_on_plane, X);
uv.y = dot(point_on_plane, Y);

// Scale for indexing
uv.x = uv.x * Constants.x;
uv.y = uv.y * Constants.x;


where:

Constants.x = 1 << depth;


Look at all those zeroes! Multiplying them out, we're left with:

// Get octant index
int side_x = P.x >= 0.0 ? 1 : 0;
int side_y = P.y >= 0.0 ? 2 : 0;
int side_z = P.z >= 0.0 ? 4 : 0;
int octant_index = side_x + side_y + side_z;

// Projection onto face
vec2 point_on_plane = abs(P.xy) / dot(abs(P), vec3(1));

// Projection onto 2D domain
vec2 uv = point_on_plane * Constants.xx;


Even better, there's no need to check for parity anymore as the triangle orientation doesn't change between rows so the nested if statements are gone:

// Get indices and fractionals
ivec2 xy = ivec2(uv);
uv -= vec2(xy);

// Double x to account for two trangles per square
xy.x *= 2;

// Offset by one for the upper triangle
xy.x = (uv.x + uv.y > 1.0) ? xy.x + 1 : xy.x;


Through determined reduction we've managed to derive the core of Octahedral Normal Vectors and use them, maybe, where they weren't expected. The final code is below. While the indices now are different -- something the colour changes visually demonstrate -- they are still unique.