Fast, Constant Time Sphere Indexing

The problem statement is simple: Assuming a triangle-subdivided sphere, map any 3D point to a triangle on that sphere 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 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.


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 viewport is interactive and code 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 this is an acceptable trade-off for the simplicity of the final implementation.

Initial Mapping Method

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

  1. Identify which octant of the Octahedron the point is on.
  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 ((ONLY WORKS DUE TO TESSELLATION METHOD).
  5. Combine these indices with 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:

int x_side = P.x >= 0.0 ? 1 : 0;
int y_side = P.y >= 0.0 ? 1 : 0;
int z_side = P.z >= 0.0 ? 1 : 0;
int octant_index = x_side + y_side * 2 + z_side * 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;
face_dir.x = float(x_side) * 2.0 - 1.0;
face_dir.y = float(y_side) * 2.0 - 1.0;
face_dir.z = float(z_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
int x_side = P.x >= 0.0 ? 1 : 0;
int y_side = P.y >= 0.0 ? 2 : 0;
int z_side = P.z >= 0.0 ? 4 : 0;
int octant_index = x_side + y_side + z_side;

// Generate face direction from sidedness
vec3 face_dir;
face_dir.x = float(x_side) * 2.0 - 1.0;
face_dir.y = float(y_side) * 1.0 - 1.0;
face_dir.z = float(z_side) * 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 \(\vec{I}\) on the octahedron face then follows from the intersection of a ray with a plane:

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

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

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

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

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

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

Where \(\vec{P}\) is the ray origin, \(\vec{D}\) is the ray direction, \(\vec{X}\) is a point on the plane and \(\vec{N}\) is the plane normal. As the ray direction is toward the octahedron origin, it's equal to the normalised \(-\vec{P}\), so:

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

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

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

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

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

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

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

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

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

The negation of the numerator/denominator also cancels:

\[\vec{I} = \vec{P} + \vec{P} { \frac{(\vec{X} - \vec{P}).\vec{N}}{ \vec{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, requring 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 using the midpoint of the meridian edge as a temporary origin:

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

\[X = |v_1 - O|\]

\[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 hypoteneuse of a right-triangle with side lengths of \(1\).

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

The vector \(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 \(X\). Or just use a multiply of the inverse:

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

The length of \(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:

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