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:

//
// 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.
// Gist here https://gist.github.com/dwilliamson/65a81bf6fcd0e2000039
// Thanks to @rompa for a couple of optimisations.
//

// Generate initial octahedron and subdivide it
var octahedron_geometry = CreateOctahedronGeometry();
var depth = 3;
for (var i = 0; i < depth; i++)
	SubdivideGeometryTriangleList(octahedron_geometry, false);

// Use nlerp to project onto sphere. Would be better to slerp to get less
// distortion but that breaks the indexing scheme. UNCOMMENT to see source
// octahedron.
ProjectVerticesToSphere(octahedron_geometry.Vertices, 1);

var mesh = scene.AddMesh(DrawType.WIREFRAME_TRIS, octahedron_geometry, undefined, OctahedronShader);

//===================================================================
//@buffer(OctahedronShader)
//===================================================================

precision highp float;

uniform vec3 glColour;

varying vec3 ls_Position;

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

int GetTriangleIndex(vec3 P)
{
	// 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;
	face_dir.x = float(side.x) * 2.0 - 1.0;
	face_dir.y = float(side.y) * 1.0 - 1.0;
	face_dir.z = float(side.z) * 0.5 - 1.0;

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

// 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.
	// Note use of unnormalised face direction as length factors out
	float d = dot(v0 - P, face_dir) / dot(P, face_dir);
	vec3 point_on_plane = P + P * d;

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

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

// Normalise plane x position to units of sub triangle half-edges
	float sub_tri_edge_len = inv_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 * inv_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;
}

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

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:

Each subdivision pass splits edges at their midpoint. Project the split onto the sphere after each pass.
Split edges at the midpoint but project onto the sphere only once after all passes complete.
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):

// -----------------------------------------------
// Configuration parameters
// -----------------------------------------------

// How many subdivision passes
var OCTAHEDRON_SUBDIVISION_DEPTH = 3;

// Toggle to see the source Octahedron
var PROJECT_ONTO_SPHERE = true;

// -----------------------------------------------
// Create initial Octahedron Geometry
// -----------------------------------------------

function CreateOctahedronGeometry()
{
	// http://en.wikipedia.org/wiki/Octahedron
	// http://en.wikipedia.org/wiki/Octahedron#mediaviewer/File:Dual_Cube-Octahedron.svg
	// Octahedron is dual to the cube with a point in the centre of each of its faces

// Create vertex positions on the 6 cube faces
	var positions = [
		[  0,  1,  0 ],
		[  0,  0, -1 ],
		[  1,  0,  0 ],
		[  0,  0,  1 ],
		[ -1,  0,  0 ],
		[  0, -1,  0 ],
	];
	var position_array = new Array();
	for (var i in positions)
	{
		var p = vec3.create();
		p[0] = positions[i][0];
		p[1] = positions[i][1];
		p[2] = positions[i][2];
		position_array.push(p);
	}

// Link octahedron triangles
	var indices = [
		0, 1, 2,
		0, 2, 3,
		0, 3, 4,
		0, 4, 1,
		5, 1, 2,
		5, 2, 3,
		5, 3, 4,
		5, 4, 1
	];
	var index_array = new Array();
	for (var i in indices)
		index_array.push(indices[i]);

return new Geometry(IndexType.TRIANGLE_LIST, position_array, index_array);
}

// -----------------------------------------------
// Sphere Subdivision
// -----------------------------------------------

function SubdivideTriangleList(positions, indices, do_slerp)
{
	var out_indices = new Array();

// Use an object as an associative map of edge pairs
	// Not ideal as JS will convert the integer keys to strings, but, it's quick and dirty...
	var MAX_NB_INDICES = 32768;
	var edge_splits = new Object();

function make_edge_key(i0, i1)
	{
		// Ensure lowest edge index is first to ensure consistent key between different edge directions
		if (i1 < i0)
		{
			var temp = i0;
			i0 = i1;
			i1 = temp;
		}

// Check for overflow
		if (i1 >= MAX_NB_INDICES)
			throw new Error("Edge split table not big enough to handle index: " + i1);

return i0 * MAX_NB_INDICES + i1;
	}

function slerp(out, a, b, t)
	{
		// Cosine of angle between two vectors
		var cos_theta = vec3.dot(a, b);
		cos_theta = Math.min(Math.max(cos_theta, -1), 1);

// Calculate interpolated angle
		var theta = Math.acos(cos_theta) * t;

var relative = vec3.create();
		vec3.scale(relative, a, cos_theta);
		vec3.sub(relative, b, relative);
		vec3.normalize(relative, relative);

vec3.scale(relative, relative, Math.sin(theta));
		vec3.scale(out, a, Math.cos(theta));
		vec3.add(out, out, relative);
	}

function split_edge(positions, i0, i1)
	{
		// Check first to see if the edge has been split before
		var edge_key = make_edge_key(i0, i1);
		var edge_split = edge_splits[edge_key];
		if (edge_split !== undefined)
			return edge_split;

// Lookup edge vertex positions
		var p0 = positions[i0];
		var p1 = positions[i1];

// Generate a midpoint
		var p01 = vec3.create();
		if (do_slerp)
			slerp(p01, p0, p1, 0.5);
		else
			vec3.lerp(p01, p0, p1, 0.5);

// Add the new vertex position and index to the edge split table
		positions.push(p01);
		edge_split = positions.length - 1;
		edge_splits[edge_key] = edge_split;

return edge_split;
	}

for (var i = 0; i < indices.length; i += 3)
	{
		// Get vertex indices of the triangle
		var i0 = indices[i + 0];
		var i1 = indices[i + 1];
		var i2 = indices[i + 2];

// Split the edges of the triangle
		var i01 = split_edge(positions, i0, i1);
		var i12 = split_edge(positions, i1, i2);
		var i20 = split_edge(positions, i2, i0);

// Add indices for the 4 new triangles
		var new_indices = [
			i0,  i01, i20,
			i01, i1,  i12,
			i20, i01, i12,
			i20, i12, i2 
		];
		for (var j in new_indices)
			out_indices.push(new_indices[j]);
	}

return out_indices;
}

// -----------------------------------------------
// Scene Setup
// -----------------------------------------------

// Create initial geometries for each of the 3 subdivision methods
var octahedron_geometry_prj = CreateOctahedronGeometry();
var octahedron_geometry = CreateOctahedronGeometry();
var octahedron_geometry_slerp = CreateOctahedronGeometry();

for (var i = 0; i < OCTAHEDRON_SUBDIVISION_DEPTH; i++)
{
    // This geometry projects onto the sphere after each subdivision step
    SubdivideGeometryTriangleList(octahedron_geometry_prj, false);
    if (PROJECT_ONTO_SPHERE)
        ProjectVerticesToSphere(octahedron_geometry_prj.Vertices, 1);

// This geometry waits until all subdivision is complete before projecting onto the sphere
    SubdivideGeometryTriangleList(octahedron_geometry, false);

// This geometry does no sphere projection, relying on slerp to subdivide instead
    SubdivideGeometryTriangleList(octahedron_geometry_slerp, true);
}

// Sphere-project the geometry that only requires it once
if (PROJECT_ONTO_SPHERE)
    ProjectVerticesToSphere(octahedron_geometry.Vertices, 1);

var mesh_prj = scene.AddMesh(DrawType.WIREFRAME_TRIS, octahedron_geometry_prj);
var mesh = scene.AddMesh(DrawType.WIREFRAME_TRIS, octahedron_geometry);
var mesh_slerp = scene.AddMesh(DrawType.WIREFRAME_TRIS, octahedron_geometry_slerp);
mesh_prj.SetPosition(-2, 0, 0);
mesh_slerp.SetPosition(2, 0, 0);

scene.SetCameraPosition(0, 0, 7);

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:

Identify which octant of the Octahedron the point faces.
Cast a ray from the query point to the Octahedron origin and find an intersection point with the Octahedron face.
Construct a basis in the plane of the face and project the intersection point into 2D, relative to bottom-left face point.
Calculate triangle column and row indices the point lies within. This is made simple due to the resultant subdivision containing equilateral triangles.
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:

var octahedron_geometry = CreateOctahedronGeometry();
var mesh = scene.AddMesh(DrawType.WIREFRAME_TRIS, octahedron_geometry, undefined, OctahedronShader);

function AddLabel(x, y, z, text)
{
	var pos = vec3_create(x, y, z);
	scene.AddFloatingText(text, pos, pos);
}

var u = 0.4;
var v = 0.5;
AddLabel(-u, -v, -u, "tri_0");
AddLabel( u,  v,  u, "tri_1");
AddLabel(-u,  v, -u, "tri_2");
AddLabel( u,  v, -u, "tri_3");
AddLabel(-u, -v,  u, "tri_4");
AddLabel( u, -v,  u, "tri_5");
AddLabel(-u,  v,  u, "tri_6");
AddLabel( u, -v, -u, "tri_7");

scene.SetCameraRotation(-0.2, -0.25, 0);

//===================================================================
//@buffer(OctahedronShader)
//===================================================================

precision highp float;

uniform vec3 glColour;

varying vec3 ls_Position;

int GetTriangleIndex(vec3 P)
{
	// 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;
	return x_side + y_side * 2 + z_side * 4;
}

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

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

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:

function CutGeometry(geom, start_index, end_index)
{
	// Create geometry which references the original
	var cut_geom = new Geometry(geom.IndexType, geom.Vertices, geom.Indices);

// Cut out the requested triangles
	cut_geom.Indices = cut_geom.Indices.slice(start_index, end_index);

// Reduce the vertex list to only those used
	var new_vertices = [ ];
	for (var i = 0; i < cut_geom.Indices.length; i++)
	{
		var index = cut_geom.Indices[i];
		new_vertices.push(cut_geom.Vertices[index]);
		cut_geom.Indices[i] = i;
	}
	cut_geom.Vertices = new_vertices;

return cut_geom;
}

function GetDocumentBackgroundColour()
{
	var background_colour = getComputedStyle(document.body).backgroundColor; 
	var m = background_colour.match(/^rgb\s*$\s*(\d+)\s*,\s*(\d+)\s*,\s*(\d+)\s*$$/i);
	return [ m[1] / 255.0, m[2] / 255.0, m[3] / 255.0 ];	
}

function TriangleCenter(geom, triangle)
{
	var a = geom.Vertices[triangle * 3 + 0];
	var b = geom.Vertices[triangle * 3 + 1];
	var c = geom.Vertices[triangle * 3 + 2];

var m = vec3.create();
	vec3.add(m, a, b);
	vec3.add(m, m, c);
	vec3.scale(m, m, 1.0 / 3.0);

return m;
}

// Subdivided octahedron face corner lays down no-overwrite stencil
var geom_face = CreateOctahedronFaceGeometry();
for (var i = 0; i < 3; i++)
    SubdivideGeometryTriangleList(geom_face, false);
var mesh_face = scene.AddMesh(DrawType.WIREFRAME, geom_face);
var gl = mesh_face.gl;
mesh_face.SetPosition(-0.25, -0.25, 0);
mesh_face.Colour = [ 0.8, 0.8, 0.8 ];
mesh_face.StencilRef = 1;
mesh_face.StencilOpZPass = gl.REPLACE;

var geom_tri_a = CutGeometry(geom_face, 123, 126)
var geom_tri_b = CutGeometry(geom_face, 171, 174);
var mesh_tri_a = scene.AddMesh(DrawType.WIREFRAME_TRIS, geom_tri_a);
var mesh_tri_b = scene.AddMesh(DrawType.WIREFRAME_TRIS, geom_tri_b);
mesh_tri_a.SetPosition(-0.25, -0.25, 0);
mesh_tri_b.SetPosition(-0.25, -0.25, 0);
mesh_tri_a.FillColour = [ 0.0, 0.8, 0.4 ];
mesh_tri_b.FillColour = [ 0.0, 0.4, 0.8 ];
mesh_tri_a.StencilRef = 1;
mesh_tri_a.StencilOpZPass = gl.REPLACE;
mesh_tri_b.StencilRef = 1;
mesh_tri_b.StencilOpZPass = gl.REPLACE;

var a = TriangleCenter(geom_tri_a, 0);
var b = vec3.create();
vec3.copy(b, a);
vec3.normalize(b, b);
vec3.scale(b, b, 1.05);
var m = scene.AddLineMesh(b, a, 0.0025, [ 0.8, 0.2, 0.4 ]);
m.SetPosition(-0.25, -0.25, 0);
m.StencilRef = 1;
m.StencilOpZPass = gl.REPLACE;

var smpos = vec3_create(-0.25, -0.25, 0);
vec3.add(smpos, smpos, b);
var sm = scene.AddSphereMesh(smpos, 0.01, 1, [ 1.0, 0.4, 0.8 ]);
sm.StencilRef = 1;
sm.StencilOpZPass = gl.REPLACE;

var a = TriangleCenter(geom_tri_b, 0);
var b = vec3.create();
vec3.copy(b, a);
vec3.normalize(b, b);
vec3.scale(b, b, 1.2);
var m = scene.AddLineMesh(b, a, 0.0025, [ 0.8, 0.2, 0.4 ]);
m.SetPosition(-0.25, -0.25, 0);m.Colour = [ 0.8, 0.2, 0.4 ];
m.StencilRef = 1;
m.StencilOpZPass = gl.REPLACE;

var smpos = vec3_create(-0.25, -0.25, 0);
vec3.add(smpos, smpos, b);
var sm = scene.AddSphereMesh(smpos, 0.01, 1, [ 1.0, 0.4, 0.8 ]);
sm.StencilRef = 1;
sm.StencilOpZPass = gl.REPLACE;

// Subdivided sphere corner
var geom_sphere = CreateOctahedronFaceGeometry();
for (var i = 0; i < 3; i++)
    SubdivideGeometryTriangleList(geom_sphere, false);
ProjectVerticesToSphere(geom_sphere.Vertices, 1);
var mesh1 = scene.AddMesh(DrawType.WIREFRAME_TRIS, geom_sphere);
mesh1.SetPosition(-0.25, -0.25, 0);
mesh1.FillColour = GetDocumentBackgroundColour();
mesh1.Colour = [ 0.5, 0.5, 0.5 ];
mesh1.StencilRef = 1;
mesh1.StencilFunc = gl.NOTEQUAL;

var geom2 = CutGeometry(geom_sphere, 123, 126)
var geom3 = CutGeometry(geom_sphere, 171, 174);
var mesh2 = scene.AddMesh(DrawType.WIREFRAME_TRIS, geom2)
var mesh3 = scene.AddMesh(DrawType.WIREFRAME_TRIS, geom3)
mesh2.SetPosition(-0.25, -0.25, 0);
mesh3.SetPosition(-0.25, -0.25, 0);
mesh2.FillColour = [ 0.0, 0.8, 0.4 ];
mesh3.FillColour = [ 0.0, 0.4, 0.8 ];

scene.SetCameraPosition(0, 0, 2.2);
scene.SetCameraRotation(-0.2, -0.25, 0);

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:

var geom_face = CreateOctahedronFaceGeometry();
var mesh_face = scene.AddMesh(DrawType.WIREFRAME, geom_face);
mesh_face.Colour = [ 0.8, 0.8, 0.8 ];

// Octahedron face points
var v0 = vec3_create(0, 0, 1);
var v1 = vec3_create(1, 0, 0);
var v2 = vec3_create(0, 1, 0);

// Meridian edge center point
var O = vec3.create();
vec3.add(O, O, v0);
vec3.add(O, O, v1);
vec3.scale(O, O, 0.5);

// Construct basis
var X = vec3.create();
vec3.sub(X, v1, v0);
vec3.normalize(X, X);
var Y = vec3.create();
vec3.sub(Y, v2, O);
vec3.normalize(Y, Y);
var Z = vec3.create();
vec3.cross(Z, X, Y);

// Place the basis on the origin point
vec3.add(X, X, v0);
vec3.add(Y, Y, v0);
vec3.add(Z, Z, v0);

var r = 0.015;
var c = r * 3;
scene.AddLineMesh(v0, X, r, [ 0.8, 0.2, 0.4 ], c);
scene.AddLineMesh(v0, Y, r, [ 0.2, 0.8, 0.4 ], c);
scene.AddLineMesh(v0, Z, r, [ 0.2, 0.4, 0.8 ], c);

scene.AddSphereMesh(v0, 0.02, 1, [ 1, 1, 1 ]);
scene.AddSphereMesh(v1, 0.02, 1, [ 1, 1, 1 ]);
scene.AddSphereMesh(v2, 0.02, 1, [ 1, 1, 1 ]);

vec3.scale(v0, v0, 1.1);
vec3.scale(v1, v1, 1.1);
vec3.scale(v2, v2, 1.2);
scene.AddFloatingText("v_0", v0);
scene.AddFloatingText("v_1", v1);
scene.AddFloatingText("v_2", v2);

scene.AddSphereMesh(O, 0.04, 1, [ 1, 1, 1 ]);
vec3.scale(O, O, 1.1);
scene.AddFloatingText("O", O);

//vec3.scale(X, X, 1.2);
//vec3.scale(Y, Y, 1.2);
X[0] = X[0] * 1.1;
Y[1] = Y[1] * 1.3;
scene.AddFloatingText("\\vec{X}", X);
scene.AddFloatingText("\\vec{Y}", Y);

scene.SetCameraPosition(0.25, 0.25, 3);
scene.SetCameraRotation(-0.2, 0.2, 0);

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

var thick = 0.005;
var col = [ 0.8, 0.8, 0.8 ];
scene.AddCircleLineMesh(64, 1, thick, col);

var p0 = vec3_create(0, 0, 1);
var p1 = vec3_create(0, -1, 0);
var p2 = vec3_create(1, 0, 0);
var p3 = vec3_create(0, 1, 0);
var p4 = vec3_create(-1, 0, 0);
scene.AddLineMesh(p2, p3, thick, col);
scene.AddLineMesh(p3, p4, thick, col);
scene.AddLineMesh(p4, p1, thick, col);
scene.AddLineMesh(p1, p2, thick, col);

scene.AddLineMesh(p1, p0, thick, col);
scene.AddLineMesh(p2, p0, thick, col);
scene.AddLineMesh(p3, p0, thick, col);
scene.AddLineMesh(p4, p0, thick, col);

var p5 = vec3_create(p1[0], p1[1], p1[2]);
vec3.add(p5, p5, p2);
vec3.scale(p5, p5, 0.5);
scene.AddLineMesh(p0, p5, thick, col, 0, 0.1);

var o = vec3_create(0, 0, 0);
scene.AddLineMesh(o, p1, thick, col, 0, 0.1);
scene.AddLineMesh(o, p2, thick, col, 0, 0.1);

scene.AddMeasure(p1, p2, col, 0.1, "l_m", -0.1, 0.05);
scene.AddMeasure(p1, p5, col, 0.12, "l_x", 0.1, 0.05);
scene.AddMeasure(p0, p5, col, 0.15, "l_y", 0.1, -0.05);

scene.SetCameraPosition(0, 0, 3);

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

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.

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

var mesh = scene.AddMesh(DrawType.WIREFRAME_TRIS, octahedron_geometry, undefined, OctahedronShader);
mesh.FloatUniforms["Rows"] = (1 << depth);
mesh.FloatUniforms["TriHalfBasesPerEdge"] = (1 << depth) * 2;
mesh.FloatUniforms["GridLUT"] = (1 << depth) * 2 - 1;

//===================================================================
//@buffer(OctahedronShader)
//===================================================================

precision highp float;

uniform vec3 glColour;
uniform float Rows;
uniform float TriHalfBasesPerEdge;
uniform float GridLUT;

varying vec3 ls_Position;

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

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

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

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

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

// Normalise plane x position to units of sub triangle half-edges
	float sub_tri_edge_len = inv_oct_side_len / TriHalfBasesPerEdge;
	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) * Rows;

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

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

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

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.

Gazoo.vrv