Given a 3D point on the sphere, map to a unique, sub-divided triangle index in O(1) with no recursive searching or table lookups.

SphereIndexing.glsl

// Some shitty boolean logic code because WebGL...
bool and(int x, int y)
{
	return mod(float(x), exp2(float(y))) != 0.0;
}


int GetTriangleIndex(vec3 position)
{
	// TODO: Add integer arithmetic to simplify the boolean logic.
	// TODO: Move to HLSL... it's just more concise.

	// Normalize input position to unit sphere to simplify required ops
	vec3 P = normalize(position);

	// Get octant index
	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;

	// Generate an un-normalised direction vector for the face from 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;

	// The normalised face direction can be used for the plane normal
	// As we're on the unit sphere, we know the length of face_dir
	float face_dir_len = sqrt(3.0);
	vec3 plane_normal = face_dir / face_dir_len;

	// Get triangle vertices for the current face
	// 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 x0 = vec3(0, 0, face_dir.z);
	vec3 x1 = vec3(face_dir.x, 0, 0);

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

	// Find a point on the octahedron face that maps to the current point on the sphere surface
	// As edge subdivisions were generated using normalized midpoints, tracing a ray back from the
	// sphere to the origin and finding where it intersects the face will yield the original position.
	// There is less distortion generating midpoints with a slerp, but that then makes this back-trace
	// technique a lot trickier to use.
	float d = dot(x0 - P, plane_normal) / dot(-P, plane_normal);
	vec3 point_on_plane = P - P * d;

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

	//
	// Make 2D basis in the face plane using the meridian edge midpoint as the origin:
	//
	//    O = (x0 + x1) * 0.5
	//
	// Since x0 and x1 have all zeroes where there is addition, this becomes:
	//
	//    O = [x0.x, 0, x1.z] * 0.5
	//
	// The 2D basis vectors are then:
	//
	//    X = normalize(x1 - O)
	//    Y = normalize(x2 - O)
	//
	// As we're on the unit sphere, the length of the meridian edge is known:
	//
	//    oct_side_len = 2.0 * sqrt(0.5)
	//
	// Normalisation of X then becomes:
	//
	//    X = (x1 - O) / (oct_side_len * 0.5)
	//
	// The final point to note is that O/X/Y are constructed using elements of face_dir
	// and zeroes. This leads to the final simplification.
	//
	vec3 X = vec3(face_dir.x / oct_side_len, 0, -face_dir.z / oct_side_len);
	vec3 Y = vec3(face_dir.x * -0.5, face_dir.y, face_dir.z * -0.5);

	//
	// What remains is normalisation of Y, which has a length equal to the octahedron triangle's
	// height. Each component of Y is composed of face_dir and has a fixed magnitude, only differing
	// in sign. This leads to:
	//
	//    oct_tri_height = sqrt(0.5 * 0.5 + 1 * 1 + 0.5 * 0.5);
	//
	float oct_tri_height = sqrt(1.5);
	Y /= oct_tri_height;

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

	// Normalise plane x position to units of sub triangle half-edges
	float sub_tri_edge_len = oct_side_len / 8.0;
	uv.x = uv.x / (sub_tri_edge_len * 0.5);

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

	// 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)
			x--;
	}
	else
	{
		if (v - u > 0.0)
			x--;
	}

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