Marching Cubes Lookup Table Decoding in GLSL Compute

This article explains the lookup table decoding algorithm used for Marching Cubes surface polygonization within OpenGL Shading Language (GLSL) compute shaders. The process extracts an isosurface from a 3D scalar field by classifying each voxel's eight corners, deriving an 8-bit configuration index, querying precomputed edge and triangle lookup tables, and linearly interpolating vertex coordinates along intersected cell edges. Executing this pipeline entirely in compute shaders allows massive parallelization and direct output into GPU buffers for immediate rendering.

The Lookup Tables

Marching Cubes relies on two standard precomputed lookup tables (originally popularized by Paul Bourke):

  1. edgeTable[256]: A 16-bit integer array indicating which of the 12 voxel edges intersect the isosurface for any given 8-bit corner configuration.
  2. triTable[256][16]: A 2D array mapping each 8-bit configuration to a list of edge indices that form triangles. Each entry contains sequences of three edge indices, terminated by -1 (or 255 in unsigned byte formats), with a maximum of 5 triangles (15 indices) per cell.

In GLSL compute shaders, these tables are commonly bound as read-only Shader Storage Buffer Objects (SSBOs), 1D/2D integer textures, or embedded directly as constant uniform arrays.

Algorithm Pipeline

The compute shader executes one thread per voxel or grid cell, decoding the tables through five discrete steps:

1. Compute the Configuration Index

Sample the scalar field at the 8 corners of the voxel. Compare each corner value against the target threshold (isolevel). If a corner's value is below (or above, depending on convention) the threshold, set the corresponding bit in an 8-bit integer mask:

uint cubeIndex = 0u;
if (cornerValues[0] < isolevel) cubeIndex |= 1u;
if (cornerValues[1] < isolevel) cubeIndex |= 2u;
if (cornerValues[2] < isolevel) cubeIndex |= 4u;
if (cornerValues[3] < isolevel) cubeIndex |= 8u;
if (cornerValues[4] < isolevel) cubeIndex |= 16u;
if (cornerValues[5] < isolevel) cubeIndex |= 32u;
if (cornerValues[6] < isolevel) cubeIndex |= 64u;
if (cornerValues[7] < isolevel) cubeIndex |= 128u;

2. Early Out Check

If cubeIndex is 0u or 255u, the voxel is entirely inside or outside the surface and generates no geometry. The thread terminates immediately:

if (cubeIndex == 0u || cubeIndex == 255u) {
    return;
}

3. Edge Intersection and Vertex Interpolation

The 12 potential edges connect defined pairs of corners:

  • Edge 0: Corner 0 to Corner 1
  • Edge 1: Corner 1 to Corner 2
  • Edge 2: Corner 2 to Corner 3
  • Edge 3: Corner 3 to Corner 0
  • Edge 4: Corner 4 to Corner 5
  • Edge 5: Corner 5 to Corner 6
  • Edge 6: Corner 6 to Corner 7
  • Edge 7: Corner 7 to Corner 4
  • Edge 8: Corner 0 to Corner 4
  • Edge 9: Corner 1 to Corner 5
  • Edge 10: Corner 2 to Corner 6
  • Edge 11: Corner 3 to Corner 7

For every edge marked active in edgeTable[cubeIndex], compute the precise vertex position along the edge using linear interpolation:

vec3 interpolateVertex(vec3 p1, vec3 p2, float val1, float val2, float isolevel) {
    if (abs(isolevel - val1) < 0.00001) return p1;
    if (abs(isolevel - val2) < 0.00001) return p2;
    if (abs(val1 - val2) < 0.00001) return p1;
    float mu = (isolevel - val1) / (val2 - val1);
    return mix(p1, p2, mu);
}

An array of 12 temporary vec3 positions holds the computed vertices for the current voxel.

4. Decoding Triangle Topology

Read the triTable at the row corresponding to cubeIndex in strides of three. Each set of three edge IDs defines a triangle. Loop until the terminator value (-1 or 255) is encountered:

// Assuming triTable is flattened: triTable[cubeIndex * 16 + i]
for (int i = 0; i < 16; i += 3) {
    int edge0 = triTable[cubeIndex * 16 + i];
    if (edge0 == -1) break;

    int edge1 = triTable[cubeIndex * 16 + i + 1];
    int edge2 = triTable[cubeIndex * 16 + i + 2];

    vec3 v0 = edgeVertices[edge0];
    vec3 v1 = edgeVertices[edge1];
    vec3 v2 = edgeVertices[edge2];

    emitTriangle(v0, v1, v2);
}

5. Output Buffer Writing

To write triangles into an output vertex buffer without CPU intervention, use an atomic counter in an SSBO:

layout(binding = 2, std430) buffer CounterBuffer {
    uint vertexCount;
};

layout(binding = 3, std430) buffer OutputVertexBuffer {
    vec3 outPositions[];
};

void emitTriangle(vec3 a, vec3 b, vec3 c) {
    uint baseIndex = atomicAdd(vertexCount, 3u);
    outPositions[baseIndex]     = a;
    outPositions[baseIndex + 1] = b;
    outPositions[baseIndex + 2] = c;
}

Optimizations for GLSL Compute

  • Texture Packing: Pack triTable into an isampler2D (format R8I or RGBA8I) of size 16x256 to leverage texture caches and minimize SSBO memory footprint.
  • Compact Storage: Flatten the 16 elements per table entry to a fixed layout so threads in a warp/subgroup follow homogeneous memory access patterns.
  • Normal Generation: Compute surface normals during the edge interpolation step either by evaluating central differences in the scalar field texture using textureOffset or by taking the cross product of the decoded triangle edges: normalize(cross(v1 - v0, v2 - v0)).