# ==Marching Cubes on the GPU==, count / scan / emit <p class="doc-sub">// status: seedling</p> 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 - Interactive Terrain Generation with Marching Cubes|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: | Buffer | One element per | Contents | Written by | |---|---|---|---| | `field` | sample | density / signed value | field generation or brush pass | | `cases` | cell | 8-bit Marching Cubes case index | count pass | | `counts` | cell | number of triangles, `0…5` | count pass | | `offsets` | cell | exclusive prefix sum of `counts` | scan pass | | `vertices` | triangle corner | position, normal, material | emit pass | | `indirect` | draw | indexed or non-indexed draw arguments | finalize pass | The cell's linear index must be shared by every pass. I normally use x-major indexing: ```c uint 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. ```glsl layout(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: ```text load 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. ```glsl uint 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`: ```c indirect.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: ```c glMemoryBarrier(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: ```glsl uint 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. | Approach | Strength | Cost / failure mode | |---|---|---| | Atomic append | one pass, minimal bookkeeping | global contention, nondeterministic order, capacity management | | Count / scan / emit | deterministic ranges, easy indirect draw | extra dispatches and scan storage | | Fixed output per cell | simplest indexing | wastes space badly for sparse fields | | Count / scan / indexed edge cache | compact final mesh | cross-workgroup deduplication is complex | ## 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)](https://dl.acm.org/doi/10.1145/37402.37422) - [Blelloch, “Prefix Sums and Their Applications” (CMU)](https://www.cs.cmu.edu/~guyb/papers/Ble93.pdf) - [CUDA C++ Programming Guide — asynchronous and cooperative execution patterns (NVIDIA)](https://docs.nvidia.com/cuda/cuda-c-programming-guide/) - [OpenGL 4.6 Core Specification — shader memory barriers and indirect commands (Khronos)](https://registry.khronos.org/OpenGL/specs/gl/glspec46.core.pdf) - [Vulkan API Specification — compute, synchronization, and indirect drawing (Khronos)](https://registry.khronos.org/vulkan/specs/latest/html/) --- Back to [[Notes/Index|Notes]] · see also [[Marching Cubes]] · [[Compute shaders]] · [[Raym - Interactive Terrain Generation with Marching Cubes|Raym]] · [[Vulkan - learning log]]