Marching Cubes on the GPU, count / scan / emit

// status: seedling

The Marching Cubes algorithm is embarrassingly parallel right up until it has to write its answer. A cell can emit anywhere from zero to five triangles, so there is no fixed address for the next cell's output. The useful pattern is to make the variable-length part explicit:

1. Count how much each cell will emit.
2. Scan those counts into exclusive output offsets.
3. Emit triangles into the ranges that the scan reserved.

That is the same shape as GPU compaction, particle spawning, and visibility lists. In Marching Cubes the CPU version can simply append to a vector; on the GPU, count–scan–emit gives every invocation a deterministic write range. It is a good first compute pipeline for Raym, and it exercises almost every important idea in Compute shaders.

The data contract

Assume a scalar field with Nx × Ny × Nz cells. There are (Nx + 1) × (Ny + 1) × (Nz + 1) samples because adjacent cells share corners. A compact implementation keeps these buffers:

The cell's linear index must be shared by every pass. I normally use x-major indexing:

1
2
3
4
5
6
7
8uint cellIndex(uvec3 c, uvec3 cellCount) {
    return c.x + cellCount.x * (c.y + cellCount.y * c.z);
}

uint sampleIndex(uvec3 p, uvec3 cellCount) {
    uvec3 sampleCount = cellCount + uvec3(1);
    return p.x + sampleCount.x * (p.y + sampleCount.y * p.z);
}

Keep the scalar field in one coordinate convention. If a chunk owns [0, N) cells, it still needs to read the positive-side sample at N for the last cell. Sharing that boundary sample is what keeps this pipeline consistent with the seam rules in Marching Cubes and Voxel rendering techniques.

Pass 1: classify and count

Each invocation owns one cell. It reads the eight corners, constructs the case bit mask, and looks up the triangle count. The full triTable is not shown here; it is the usual table from Marching Cubes. Keeping cases is optional — recomputing the case in the emit pass saves a buffer, while storing it saves eight scalar reads per cell.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32layout(local_size_x = 8, local_size_y = 8, local_size_z = 8) in;

layout(std430, binding = 0) readonly buffer Density {
    float value[];
};
layout(std430, binding = 1) writeonly buffer CellCases {
    uint caseIndex[];
};
layout(std430, binding = 2) writeonly buffer CellCounts {
    uint triangleCount[];
};

layout(push_constant) uniform Params {
    uvec3 cellCount;
    float isoLevel;
} params;

void main() {
    uvec3 c = gl_GlobalInvocationID;
    if (any(greaterThanEqual(c, params.cellCount))) return;

    uint cell = cellIndex(c, params.cellCount);
    uint mask = 0u;
    for (uint corner = 0u; corner < 8u; ++corner) {
        uvec3 p = c + cornerOffset[corner];
        float d = value[sampleIndex(p, params.cellCount)];
        if (d < params.isoLevel) mask |= 1u << corner;
    }

    caseIndex[cell] = mask;
    triangleCount[cell] = triCount[mask]; // 0…5, derived from triTable
}

There are two small details here that are easy to miss:

• The dispatch is rounded up to a whole number of workgroups, so the bounds check is mandatory.
• The caseIndex convention has to match the corner and edge order used by triTable. A table from one implementation combined with offsets from another produces plausible-looking but broken meshes.

Pass 2: exclusive prefix scan

The scan turns counts such as [2, 0, 1, 4] into offsets [0, 2, 2, 3]. The total output is offset[last] + count[last]. Exclusive offsets are convenient because the first cell writes at zero and each cell receives its own base address.

For a small array, a workgroup can use a Blelloch scan in shared memory:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22load one count per invocation into shared x[]

// up-sweep: build a reduction tree
for offset = 1; offset < localSize; offset *= 2:
    barrier()
    if localIndex >= offset:
        x[localIndex] += x[localIndex - offset]

barrier()
if localIndex == localSize - 1:
    workgroupTotal = x[localIndex]
    x[localIndex] = 0                 // exclusive scan

// down-sweep: propagate prefixes
for offset = localSize / 2; offset > 0; offset /= 2:
    barrier()
    if localIndex >= offset:
        t = x[localIndex - offset]
        x[localIndex] += t
barrier()

offsets[cell] = x[localIndex] + prefixOfPreviousWorkgroups

The exact in-place indexing differs between a textbook Blelloch tree and a simple Hillis–Steele loop; the invariant matters more than the spelling: every shared-memory write must be followed by the appropriate barrier(), and the operation must be associative. Addition is safe. Floating-point min/max scans are also possible, but beware of NaNs.

For a terrain-sized buffer there are more cells than fit in one workgroup. Use a hierarchical scan:

1. Scan each block in shared memory, writing per-cell offsets and one blockTotal.
2. Recursively scan blockTotal (the array is much smaller each level).
3. Add the scanned block prefix to every offset in that block.

This is three or four dispatches depending on the size. A subgroup inclusive/exclusive add can replace much of the shared-memory tree when the target supports subgroup operations; see the subgroup discussion in Compute shaders. The scan's extra passes are often cheaper than a contended global atomic, especially when a dense terrain produces millions of triangles.

Choosing an output representation

The simplest output is three fresh vertices per triangle. It wastes duplicate edge vertices but makes the emit pass independent: each cell writes triangleCount[cell] × 3 consecutive vertices, with no cross-workgroup coordination. That is usually the right first version for an interactive sculpting prototype.

An indexed mesh needs a second problem: neighbouring cells share edge intersections, but workgroups cannot safely use a normal hash map without a carefully designed GPU table and collision policy. Options are:

• Emit non-indexed triangles, then run a GPU or CPU deduplication pass.
• Give every canonical grid edge a slot and compute that slot directly from the cell coordinate and edge axis. This is deterministic, but the edge buffer is larger than the final vertex set.
• Use a hash table with atomics, accepting contention and a more complicated failure path.

I would start with non-indexed output and measure. Upload bandwidth and vertex shader cost are easier to reason about than a lock-free edge cache.

Pass 3: emit

The emit pass reads the case and offset, interpolates the crossed edges, and writes only inside its reserved slice. No append counter is necessary, so the order of cells does not matter.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22uint baseVertex = offsets[cell] * 3u;
uint writeVertex = 0u;

for (uint t = 0u; triTable[mask][t] != -1; t += 3u) {
    for (uint corner = 0u; corner < 3u; ++corner) {
        uint edge = triTable[mask][t + corner];
        uint a = edgeCorners[edge].x;
        uint b = edgeCorners[edge].y;

        float va = cornerValue[a];
        float vb = cornerValue[b];
        float denominator = vb - va;
        float u = abs(denominator) < 1e-7
                ? 0.5
                : clamp((isoLevel - va) / denominator, 0.0, 1.0);

        Vertex v;
        v.position = mix(cornerPosition[a], cornerPosition[b], u);
        v.normal = normalize(gradientAt(v.position));
        vertices[baseVertex + writeVertex++] = v;
    }
}

gradientAt can sample the scalar field with central differences, as in Marching Cubes, or interpolate precomputed gradients. Central differences cost more reads but make edits and chunk boundaries straightforward. If the field is generated procedurally, evaluate the same function in the gradient pass rather than sampling a stale texture.

The last pass writes the actual draw count. It can read the final offset and count, then populate a DrawArraysIndirectCommand or DrawIndexedIndirectCommand:

1
2
3
4indirect.vertexCount = totalTriangles * 3;
indirect.instanceCount = 1;
indirect.firstVertex = 0;
indirect.firstInstance = 0;

On OpenGL, call glMemoryBarrier before an indirect draw. The relevant bits depend on the buffers' later uses; a typical non-indexed path is:

1
2
3glMemoryBarrier(GL_SHADER_STORAGE_BARRIER_BIT |
                GL_VERTEX_ATTRIB_ARRAY_BARRIER_BIT |
                GL_COMMAND_BARRIER_BIT);

Add GL_ELEMENT_ARRAY_BARRIER_BIT for an index buffer. In Vulkan, insert a barrier from compute shader writes to the next use: COMPUTE_SHADER / SHADER_WRITE to DRAW_INDIRECT / INDIRECT_COMMAND_READ, and to VERTEX_INPUT / VERTEX_ATTRIBUTE_READ for generated vertex data. The exact vkCmdPipelineBarrier2 stages and access masks belong to the resource's actual use, not to a cargo-cult global barrier; compare Vulkan - learning log.

A fast but less scalable alternative: atomic append

One dispatch can classify, interpolate, and reserve output with an atomic counter:

1
2
3
4
5uint n = triCount[mask] * 3u;
uint base = atomicAdd(vertexCount, n);
if (base + n <= vertexCapacity) {
    writeTriangles(base, mask);
}

This is excellent for proving the mesh path and for sparse fields where only a few cells emit. It degrades when many invocations hit the same counter. Reserve one counter per workgroup and combine block totals later if the global atomic becomes the bottleneck. Also put a hard capacity check around the write; an overrun corrupts unrelated resources and can look like a shader or synchronization bug.

Making it interactive

For a Raym-style edit loop, do not remesh the entire world because one brush touched one voxel. Keep a dirty-chunk queue, dispatch the three passes per dirty chunk (or batch several chunks into one dispatch), and include one-cell ghost samples around each chunk. The CPU can update the scalar field, or a brush compute pass can do it in place; after the brush, mark the affected chunks and their boundary neighbours dirty.

Useful measurements are:

• scalar reads and writes per cell;
• scan time versus count and emit time;
• generated triangle count and vertex bandwidth;
• occupancy / register pressure for the interpolation and gradient code;
• time spent waiting for a CPU readback (ideally zero).

Do not read the triangle count back just to call a draw. Use an indirect draw, and keep a conservative maximum allocation per chunk if the API lacks an indirect-count feature. If a CPU-side collision mesh is also required, make that an explicit asynchronous copy; rendering should not stall on collision data.

Things that tripped me up

• Counting triangles is not counting vertices. A cell count of n means 3n non-indexed vertices. Mixing those units shifts every later write.
• The positive boundary sample is real data. A chunk with N³ cells needs (N+1)³ corner samples, even if the last plane is shared with its neighbour.
• An exclusive scan has a different total than an inclusive scan. Test [2, 0, 1] and assert offsets [0, 2, 2], total 3 before trying a terrain.
• A barrier cannot repair a bad dependency. barrier() only synchronises invocations inside one workgroup. Cross-dispatch visibility needs an API memory barrier; cross-workgroup communication needs a different algorithm.
• The case table is a contract. Corner winding, inside/outside sign, edge numbering, and triangle winding must all come from the same convention.
• A fast GPU path can still lose to a CPU path. Small dirty regions pay dispatch and synchronization overhead. Measure the edit size at which GPU generation wins.

References

• Lorensen & Cline, “Marching Cubes: A High Resolution 3D Surface Construction Algorithm” (ACM)
• Blelloch, “Prefix Sums and Their Applications” (CMU)
• CUDA C++ Programming Guide — asynchronous and cooperative execution patterns (NVIDIA)
• OpenGL 4.6 Core Specification — shader memory barriers and indirect commands (Khronos)
• Vulkan API Specification — compute, synchronization, and indirect drawing (Khronos)

Back to Notes · see also Marching Cubes · Compute shaders · Raym · Vulkan - learning log