Marching Tetrahedra in GLSL: Solving Ambiguities

This article explains how to implement the Marching Tetrahedra algorithm in GLSL to eliminate the topological ambiguities and hole artifacts common to Marching Cubes. It details the mathematical cause of these ambiguities, demonstrates how decomposing cubic cells into tetrahedra solves the problem, and provides a structured approach for implementing lookup tables, edge interpolation, and triangle generation in GPU shaders.

The Topological Ambiguity Problem in Marching Cubes

The standard Marching Cubes algorithm evaluates a scalar field at the eight vertices of a grid cell to approximate an isosurface. Because each vertex can be either inside or outside the surface, there are \(2^8 = 256\) possible configurations, typically reduced to 15 canonical cases via symmetry.

Ambiguities occur when diagonally opposite vertices on a shared square face have identical signs while the remaining two vertices have the opposite sign (notably in cases 3, 6, 7, 10, 12, and 13). Because a square face can be legitimately contoured in two conflicting ways, adjacent cubes evaluating the shared face independently can choose mismatched surface topologies. This mismatch introduces visible cracks and non-manifold geometry in the rendered mesh unless complex asymptotic deciders or extended 33-case lookup tables are applied.

How Marching Tetrahedra Resolves Ambiguity

Marching Tetrahedra circumvents this issue by subdividing each cubic cell into smaller simplices (tetrahedra). A tetrahedron has only four vertices, yielding \(2^4 = 16\) possible binary states. When reduced by rotational symmetry, only three distinct geometric topologies exist:

  1. All vertices inside or outside (0 or 4 active): No surface passes through the tetrahedron (0 triangles).
  2. One vertex inside or outside (1 or 3 active): The surface cuts off one corner, producing a single triangle (1 triangle).
  3. Two vertices inside and two outside (2 active): The surface cuts through four edges, forming a quadrilateral that splits into two triangles (2 triangles).

Because every shared face between tetrahedra is a triangle, no configuration allows alternating diagonal signs across a face. A triangular face can only ever have zero, one, or two vertices marked active, each having a mathematically unique intersection line. Consequently, topological ambiguity is impossible, and generated meshes are inherently manifold and watertight.

Cell Decomposition Strategies

A regular cube can be partitioned into tetrahedra in two primary ways:

  • 5-Tetrahedra Decomposition: Requires alternating the orientation of the subdivision across adjacent cells (using a 3D parity check: (x + y + z) % 2) to ensure that internal diagonal cuts align across cube boundaries.
  • 6-Tetrahedra Decomposition: Standardizes edge alignment across all cells along a primary diagonal axis without requiring parity alternating. While it produces more triangles, it is often simpler to vectorize in compute shaders.

For the 6-tetrahedron approach, a cube with local corners labeled 0 to 7 can be split into six tetrahedra using consistent index tuples:

  • Tet 0: (0, 5, 1, 6)
  • Tet 1: (0, 1, 2, 6)
  • Tet 2: (0, 2, 3, 6)
  • Tet 3: (0, 3, 7, 6)
  • Tet 4: (0, 7, 4, 6)
  • Tet 5: (0, 4, 5, 6)

GLSL Implementation Architecture

Because modern GLSL implementations typically process volumetric data in Compute Shaders or Geometry Shaders, Marching Tetrahedra can be implemented cleanly with constant lookup tables.

1. Case and Edge Tables

A 16-element lookup table determines the edge intersections for any 4-bit configuration. Each tetrahedron contains six edges:

  • Edge 0: Vertex 0–1
  • Edge 1: Vertex 1–2
  • Edge 2: Vertex 2–0
  • Edge 3: Vertex 0–3
  • Edge 4: Vertex 1–3
  • Edge 5: Vertex 2–3

The table maps the active vertex bitmask to the edges forming the resulting triangles (terminated by -1):

// Format: Up to 2 triangles per tetrahedron (6 edge indices max)
const int tetTriTable[16][7] = int[16][7](
    int[7](-1, -1, -1, -1, -1, -1, -1), // 0000: None
    int[7]( 0,  2,  3, -1, -1, -1, -1), // 0001: 1 triangle
    int[7]( 0,  1,  4, -1, -1, -1, -1), // 0010: 1 triangle
    int[7]( 1,  4,  2,  2,  4,  3, -1), // 0011: 2 triangles
    int[7]( 1,  2,  5, -1, -1, -1, -1), // 0100: 1 triangle
    int[7]( 0,  5,  3,  0,  1,  5, -1), // 0101: 2 triangles
    int[7]( 0,  2,  4,  4,  2,  5, -1), // 0110: 2 triangles
    int[7]( 3,  4,  5, -1, -1, -1, -1), // 0111: 1 triangle
    int[7]( 3,  4,  5, -1, -1, -1, -1), // 1000: 1 triangle
    int[7]( 0,  4,  2,  4,  5,  2, -1), // 1001: 2 triangles
    int[7]( 0,  3,  5,  0,  5,  1, -1), // 1010: 2 triangles
    int[7]( 1,  5,  2, -1, -1, -1, -1), // 1011: 1 triangle
    int[7]( 1,  2,  4,  2,  3,  4, -1), // 1100: 2 triangles
    int[7]( 0,  4,  1, -1, -1, -1, -1), // 1101: 1 triangle
    int[7]( 0,  3,  2, -1, -1, -1, -1), // 1110: 1 triangle
    int[7](-1, -1, -1, -1, -1, -1, -1)  // 1111: None
);

2. Linear Edge Interpolation

When an edge crosses the isosurface threshold, the vertex position is linearly interpolated based on the scalar field values at the edge endpoints:

vec3 interpolateEdge(vec3 p1, vec3 p2, float v1, float v2, float isoLevel) {
    if (abs(isoLevel - v1) < 0.00001) return p1;
    if (abs(isoLevel - v2) < 0.00001) return p2;
    if (abs(v1 - v2) < 0.00001) return p1;
    
    float t = (isoLevel - v1) / (v2 - v1);
    return mix(p1, p2, t);
}

3. Tetrahedron Evaluation Loop

In the main shader routine (such as a compute invocation per voxel), sample the eight corners of the current voxel, assemble the tetrahedra, and process each individually:

void processTetrahedron(vec3 p[4], float v[4], float isoLevel) {
    int caseIndex = 0;
    if (v[0] > isoLevel) caseIndex |= 1;
    if (v[1] > isoLevel) caseIndex |= 2;
    if (v[2] > isoLevel) caseIndex |= 4;
    if (v[3] > isoLevel) caseIndex |= 8;

    // Edge-to-vertex index mappings
    const ivec2 edgeVerts[6] = ivec2[6](
        ivec2(0, 1), ivec2(1, 2), ivec2(2, 0),
        ivec2(0, 3), ivec2(1, 3), ivec2(2, 3)
    );

    for (int i = 0; i < 6; i += 3) {
        if (tetTriTable[caseIndex][i] == -1) break;

        // Generate triangle vertices by interpolating along active edges
        for (int j = 0; j < 3; j++) {
            int edge = tetTriTable[caseIndex][i + j];
            int e0 = edgeVerts[edge].x;
            int e1 = edgeVerts[edge].y;
            vec3 vertPos = interpolateEdge(p[e0], p[e1], v[e0], v[e1], isoLevel);
            
            // Output vertex to Shader Storage Buffer Object (SSBO) or append stream
            writeVertex(vertPos);
        }
    }
}

Performance and Trade-Offs

While Marching Tetrahedra guarantees manifold, hole-free meshes without branch-heavy ambiguity deciders, it comes with specific architectural considerations:

  • Triangle Count: Marching Tetrahedra generates approximately 1.5 to 2 times more triangles than Marching Cubes for the same volume resolution.
  • Vertex Redundancy: The 6-fold decomposition creates shared internal diagonals, which can produce skinny or poorly conditioned triangles without post-generation smoothing or vertex welding via compute-based hash tables.
  • Instruction Simplicity: The lookup logic contains only 16 cases and simple branch-free bit operations, which execute with high warp coherence on modern GPUs compared to the 256-entry table lookups and conditional tests of disambiguated Marching Cubes variants.