Fast, Constant Time Sphere Indexing, Part 1

The problem statement is: Assuming a triangle-subdivided sphere, map any 3D point on or above it to a triangle index in constant-time, without using recursion, table lookups or complicated branching logic. Even better; make it so simple you can use it in a pixel shader. This is useful for when your playing field is on/above a sphere and you have some lookup tables you want to reference based on player/camera position.

This is typically an offline problem considered in map projections of planets and their skies. It allows you to map the sphere surface to a 2D domain and perform all manner of simulations on it. As such, there are a lot of complicated solutions that give great results with expensive requirements. One example is HEALPix but we need something far simpler.

A more recent development in the realtime world is Spherical Fibonacci Mapping. This will map any point on the sphere to its closest Spherical Fibonacci point in constant-time but may not be fast enough.

The view below shows an interactive result of this post:

The index of each triangle is calculated in the fragment shader using only the 3D world position and then mapped to an arbitrary colour:

void main(void)
    int tri_index = GetTriangleIndex(ls_Position);
    tri_index = int(mod(float(tri_index), 8.0));

    vec3 debug_colour;
    if (tri_index == 0) debug_colour = vec3(0.7, 0.5, 0.2);
    if (tri_index == 1) debug_colour = vec3(0.6, 0.7, 0.2);
    if (tri_index == 2) debug_colour = vec3(0.2, 0.7, 0.2);
    if (tri_index == 3) debug_colour = vec3(0.2, 0.7, 0.6);
    if (tri_index == 4) debug_colour = vec3(0.2, 0.5, 0.7);
    if (tri_index == 5) debug_colour = vec3(0.3, 0.2, 0.7);
    if (tri_index == 6) debug_colour = vec3(0.7, 0.2, 0.7);
    if (tri_index == 7) debug_colour = vec3(0.7, 0.2, 0.3);

    gl_FragColor = vec4(glColour + debug_colour * 1.1, 1);

We will be implementing GetTriangleIndex with a configurable sphere subdivision count.

Subdividing the Sphere

Each triangle on the sphere should be as close to equal area as possible, with minimal angular distortion for a balanced distribution. A good approximation is to start with an Octahedron and subdivide its edges, creating 4 new triangles for each triangle, until the required density is met.

Three of the many ways of performing this subdivision are:

  1. Each subdivision pass splits edges at their midpoint. Project the split onto the sphere after each pass.
  2. Split edges at the midpoint but project onto the sphere only once after all passes complete.
  3. Split edges using a Vector Slerp with no sphere projection required.

The three techniques are implemented below (the code to the right can be edited live):

Notable is that methods 1 and 3 result in the same subdivision, however it's the second that we'll be using. It suffers from slight contraction at the poles but is an acceptable trade-off as it allows the technique to work without any non-uniform direction remapping.

Initial Mapping Method

Assuming a point floating above the unit sphere, the basic method is:

  1. Identify which octant of the Octahedron the point faces.
  2. Cast a ray from the query point to the Octahedron origin and find an intersection point with the Octahedron face.
  3. Construct a basis in the plane of the face and project the intersection point into 2D, relative to bottom-left face point.
  4. Calculate triangle column and row indices the point lies within. This is made simple due to the resultant subdivision containing equilateral triangles.
  5. Combine these indices with the octant index to give a unique index for the final triangle.

The Octahedron looks like this with colour-coded octant indices from 0 to 8:

Finding the octant index from a unit sphere point is a simple case of inspecting the signs of its individual elements, as this point is also a direction from the octahedron:

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

Given that the backward raycast undoes the final projection of any subdivided points on to the sphere, we end up back on a uniformly subdivided plane of the octahedron face. This can be visualised:

In order to intersect the ray with the face plane, the normal and a point on the plane need to be determined. The face normal is derived easily from face sidedness:

vec3 face_dir = vec3(side) * 2.0 - 1.0;
vec3 plane_normal = normalize(face_dir);

The calculation of octant_index can be further simplified by moving the multiplcations to the side select and accounting for scale in the face_dir calculation:

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

A point in the plane can be taken from one of the 3 triangle vertices in this octant. We're going to need all 3 vertices later in the function so we might as well calcuate them up front here. As we're on the unit octahedron, these points come simply from the face direction as they're already in the range \([-1, 1]\):

// Winding order relative to other faces is irrelevant so we're free to assume all meridian
// edges start at one of the two points on the z-axis and point toward one of the two
// points on the x-axis...
vec3 v0 = vec3(0, 0, face_dir.z);
vec3 v1 = vec3(face_dir.x, 0, 0);

// ...the last vertex is one of the two poles on the y-axis
vec3 v2 = vec3(0, face_dir.y, 0);

The intersection point \(I\) on the octahedron face can then be calculated. Start with the equations of a ray and plane:

\[R = P + \vec{D}t\]

\[(V - R).\vec{N} = 0\]

where \(R\) is an arbitrary point in 3D, \(P\) is the ray origin, \(\vec{D}\) is the ray direction, \(t\) is distance along the ray, \(V\) is a point on the plane and \(\vec{N}\) is the plane normal. Substitute \(R\) in the plane equation with the ray and make \(t\) the subject:

\[(V - P - \vec{D}t).\vec{N} = 0\]

\[(V - P).\vec{N} = \vec{D}t.\vec{N}\]

\[t =\frac{(V - P).\vec{N}}{\vec{D}.\vec{N}}\]

This can be put back into the ray equation to solve for \(I\):

\[I = P + \vec{D}t\]

As the ray direction is toward the octahedron origin, it's equal to the normalised \(-P\), so:

\[\hat{P} = \frac{P}{|P|}\]

\[t =\frac{(V - P).\vec{N}}{-\hat{P}.\vec{N}}\]

\[I = P - \hat{P}t\]

This can be simplified a little by substituting for \(t\):

\[I = P - \hat{P} { \frac{(V - P).\vec{N}}{ -\hat{P}.\vec{N} } }\]

and substituting for \(\hat{P}\) to see that the normalisation cancels and is not needed:

\[I = P - \frac{P}{|P|} { \frac{(V - P).\vec{N}}{ -\frac{P}{|P|}.\vec{N} } }\]

\[I = P - \frac{P.|P|}{|P|} { \frac{(V - P).\vec{N}}{ -P.\vec{N} } }\]

\[I = P - P { \frac{(V - P).\vec{N}}{ -P.\vec{N} } }\]

The negation of the numerator/denominator also cancels:

\[I = P + P { \frac{(V - P).\vec{N}}{ P.\vec{N} } }\]

leading to the code:

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

As a bonus the calculation of the plane normal using the face direction length similarly factors out. This allows use of the raw face direction instead of having to normalise it as already done above:

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

This point then needs to be projected to 2D on the plane, requiring a completed basis. The red and green unit vectors below are the vectors required, as the blue vector is the plane normal. These vectors must originate from one of the tetrahedron vertices for the 2D projection co-ordinates to remain positive:

Getting the basis vectors can be achieved starting with the midpoint of the meridian edge, \(O\), as a temporary origin and then shifting them to \(v_0\):

\[O = (v_0 + v_1) * 0.5\]

\[\vec{X} = |v_1 - O|\]

\[\vec{Y} = |v_2 - O|\]

// Make 2D basis in the face plane using the meridian edge midpoint as the origin:
vec3 O = (v0 + v1) * 0.5;
vec3 X = normalize(v1 - O);
vec3 Y = normalize(v2 - O);

This requires two normalisations that can be avoided by taking advantage of the fact that the octahedron is circumscribed by the unit sphere. All octahedron side lengths are thus equal to the hypotenuse of a right-triangle with side lengths of \(1\).

\[l_m = \sqrt{1^2 + 1^2}\]

The vector \(\vec{X}\) is half this size since it's calculated using the edge midpoint as its origin:

\[l_x = 0.5 \sqrt{2}\]

So just divide by this length to normalise \(\vec{X}\). Or just use a multiply of the inverse:

\[\frac{1}{l_x} = \frac{2}{\sqrt{2}}\]

The length of \(\vec{Y}\) is then equal to the height of the initial octahedron triangle:

\[l_x^2 = l_y^2 + \sqrt{2}^2\]

\[\sqrt{2}^2 = l_y^2 + (0.5 \sqrt{2})^2\]

\[2 = l_y^2 + 0.5\]

\[l_y = \sqrt{1.5}\]

\[\frac{1}{l_y} = \frac{1}{\sqrt{1.5}}\]

This neatly reduces to:

// Assume octahedron is circumscribed by the unit sphere, length of its meridian edges is known
float inv_oct_side_len = 2.0 / sqrt(2.0);

// What remains is normalisation of Y, which has a length equal to the octahedron triangle's height.
float inv_oct_tri_height = 1.0 / sqrt(1.5);

// Make 2D basis in the face plane using the meridian edge midpoint as the origin
vec3 O = (v0 + v1) * 0.5;
vec3 X = (v1 - O) * inv_oct_side_len;
vec3 Y = (v2 - O) * inv_oct_tri_height;

The constructed basis can now be used to project onto the plane, relative to the intended origin \(v_0\):

// Project the intersection point onto this plane with tetrahedron point origin
vec2 uv;
uv.x = dot(point_on_plane - v0, X);
uv.y = dot(point_on_plane - v0, Y);

uv gives us an effective distance of the point along each axis.

The penultimate step is to map this 2D position to a unique triangle index on the face:

There are 8 rows in the above face and 8 triangles with their base touching the bottom edge. Given a subdivision depth \(d\), the row index is simply:

\[row = \left \lfloor{\frac{l_y}{2^d}}\right \rfloor\]

where \(d=0\) represents no subdivision. Calculating the triangle index within a row is slightly more involved. The method works by splitting each triangle in half and mapping rectangles over the top:

Treat the rectangles in pairs where the first triangle has a rising diagonal (left-to-right) and the second triangle has a falling diagonal. Find where the point is within the rectangle it falls upon (\(\{u,v \mid 0 \le u,v < 1 \}\)) and the following rules determine which half of the sub-triangle the point is in:

  • For even rectangles, if \(v - u > 0\) the point is in the top sub-triangle.
  • For odd triangles, if \(u + v < 1\) the point is in the bottom sub-triangle.

You can then use that decision to shift the rectangle index down by one to correct to the triangle index within the row:

// Normalise plane x position to units of sub triangle half-edges
float sub_tri_edge_len = inv_oct_side_len / 16;
uv.x = uv.x / sub_tri_edge_len;

// Normalise plane y position to units of sub triangle heights (the easy bit)
uv.y = (uv.y * inv_oct_tri_height) * 8;

// Get integer indices, y is already in its final form
int x = int(uv.x);
int y = int(uv.y);

// Get fractionals
float u = uv.x - float(x);
float v = uv.y - float(y);

// Assuming a grid of equilateral triangles (guaranteed with octahedron midpoint/normalize subd)
// Shift x index 1 to the left for each rise above the diagonals. Need to alternate diagonal
// direction based on parity of y index.
if (And(x, 1) != And(y, 1))
    if (u + v < 1.0)
    if (v - u > 0.0)

The and function here is a terrible little workaround function for GLSL in WebGL 1 that can't do bitwise integer arithmetic. It's available in the final example below.

The triangle row and index now need to be turned into a linear triangle index within the octant and that's trivially achieved by treating the octant as a grid and ignoring indices outside the octant triangle. The grid width is a function of the subdivision depth, noting that one less triangle than rectangles can be stored:

\[gridwidth = 2^d * 2 - 1\]

Which can then be combined with the octant index to give the final linear index:

// This is a nice way of making sure that triangles on each row start off at index 0
//    x -= y;
// However, can't easily turn that into a linear index as row starting indices form the following sequence
//    0, 15, 28, 39, 48, 55, 60, 63
// Instead, juse assume a square grid and only store data for those triangles within the octahedron face bounds.
return octant_index * 15 * 15 + y * 15 + x;

The resultant code can be seen and edited below. Setting depth to something high like 6 highlights the difference in projected area between big triangles in the middle of an octant and and small triangles on the edge.

While this code is pretty fast (pixel-shader fast, even) it can be optimised much more. These optimisations will be covered in a further post.

Thanks to Mark Wayland, Mykhailo Parfeniuk, Randall Rauwendaal, Mike Ducker and Kai Jouran for optimisations, fixes, suggestions and proof-reading.

The code for this website and sandbox is open-source on Github. I will be hacking on the sandbox code over the next few weeks to try and simplify things a little; the result of which may make its way to its own repo.