3D Jump Flooding Distance Fields in GLSL Shaders

The 3D Jump Flooding Algorithm (JFA) allows GPU compute shaders to generate high-performance Euclidean Distance Fields (EDF) across three-dimensional voxel grids in logarithmic time. By iteratively propagating nearest-feature coordinates across exponentially decreasing spatial offsets, JFA transforms a computationally prohibitive \(O(N^6)\) global distance check into an efficient \(O(N^3 \log N)\) parallel routine. This article breaks down how the 3D jump flooding scheme is structured, executed, and refined within GLSL compute shaders to produce accurate Euclidean distance volumes.

Initialization and Data Representation

The calculation begins by allocating two 3D textures of identical resolution (\(N \times N \times N\)) to act as ping-pong render targets via GLSL's image3D interface. Unlike simple distance grids that store scalar values, JFA requires storing the spatial coordinates of the closest seed features.

  1. Seed Initialization: A setup compute pass processes the initial boundary geometry or voxelized mesh.
  2. Buffer Encoding: Voxels containing boundary geometry are written with their own 3D grid coordinates ivec3(x, y, z). Voxels with no geometry are tagged with an invalid or infinite coordinate vector, such as ivec3(-1).
  3. Format Selection: A high-precision integer or floating-point format like RGBA16UI, RGBA32UI, or RGBA32F is utilized to hold the \((x, y, z)\) coordinates of the closest identified seed for every voxel.

Logarithmic Step Propagation

The core principle of jump flooding is passing candidate seed positions across large distances first, progressively narrowing the search stride. For a grid with a maximum dimension of \(N\), the algorithm runs for \(M = \lceil \log_2(N) \rceil\) passes.

The CPU manages the dispatch loop, binding one texture as the read source and the other as the write destination, alternating their roles each iteration. The step length \(k\) begins at half the volume dimension:

\[k = \left\lfloor \frac{N}{2} \right\rfloor\]

With each successive dispatch, \(k\) is halved (\(k = \lfloor k / 2 \rfloor\)) until the step size reaches \(1\).

The Compute Shader Kernel Execution

During each pass, a GLSL compute shader evaluates the \(3 \times 3 \times 3\) neighborhood around each voxel. With a work group configuration such as layout(local_size_x = 8, local_size_y = 8, local_size_z = 8), each thread is mapped directly to a voxel at gl_GlobalInvocationID:

layout(local_size_x = 8, local_size_y = 8, local_size_z = 8) in;

layout(binding = 0, rgba32i) uniform readonly iimage3D u_InputGrid;
layout(binding = 1, rgba32i) uniform writeonly iimage3D u_OutputGrid;

uniform int u_StepSize;
uniform ivec3 u_GridResolution;

void main() {
    ivec3 currentPos = ivec3(gl_GlobalInvocationID.xyz);
    if (any(greaterThanEqual(currentPos, u_GridResolution))) {
        return;
    }

    ivec3 bestSeed = imageLoad(u_InputGrid, currentPos).xyz;
    float bestDist = (bestSeed.x >= 0) ? distance(vec3(currentPos), vec3(bestSeed)) : 1e20;

    // Sample the 26 neighbors around the current cell
    for (int z = -1; z <= 1; ++z) {
        for (int y = -1; y <= 1; ++y) {
            for (int x = -1; x <= 1; ++x) {
                if (x == 0 && y == 0 && z == 0) continue;

                ivec3 neighborPos = currentPos + ivec3(x, y, z) * u_StepSize;

                if (all(greaterThanEqual(neighborPos, ivec3(0))) &&
                    all(lessThan(neighborPos, u_GridResolution))) {
                    
                    ivec3 candidateSeed = imageLoad(u_InputGrid, neighborPos).xyz;
                    if (candidateSeed.x >= 0) {
                        float dist = distance(vec3(currentPos), vec3(candidateSeed));
                        if (dist < bestDist) {
                            bestDist = dist;
                            bestSeed = candidateSeed;
                        }
                    }
                }
            }
        }
    }

    imageStore(u_OutputGrid, currentPos, ivec4(bestSeed, 0));
}

By querying coordinates at offset intervals of u_StepSize, voxels incorporate seeds that originated far across the volume without requiring sequential voxel-by-voxel propagation.

Preserving True Euclidean Accuracy

A classic limitation of Jump Flooding is the emergence of propagation errors caused by grid-alignment blind spots, where a correct seed path is obstructed by an intermediate neighbor that favors an alternate seed. To achieve high Euclidean accuracy, two adjustments are applied:

  1. Storing Seed Positions Rather than Distance Accumulations: Instead of computing a distance field through dynamic programming (such as adding the neighbor step length to a scalar value), the shader explicitly propagates the unquantized index of the origin seed. The metric calculated is always the true continuous vector magnitude: \[\text{Distance} = \|\vec{P}_{\text{current}} - \vec{P}_{\text{seed}}\|_2\]
  2. JFA Variants (JFA+1 and JFA+2): Running an additional iteration at step size \(k=1\) or \(k=2\) (repeating passes after the primary sequence finishes) repairs small sampling errors. This reduces geometric inaccuracies across the Voronoi boundaries to near zero.

Final Distance Field Generation

Once the step size reaches \(1\) and any correction passes are complete, the resulting texture contains the exact coordinate of the closest seed for every voxel in the 3D volume. A lightweight final compute pass transforms this coordinate map into a scalar distance field by evaluating the magnitude of the offset vector:

float finalDistance = length(vec3(currentPos) - vec3(closestSeed));

The resulting value can then be normalized, converted to a signed distance field (SDF) by testing inside/outside status against the source geometry, and saved into a single-channel floating-point 3D texture (R16F or R32F) ready for volumetric rendering, collision detection, or ambient occlusion sampling.