Calculating 3D Morton Codes in GLSL for LBVH

This article explains how to compute 3D Morton codes from integer coordinates using GLSL to facilitate Linear Bounding Volume Hierarchy (LBVH) construction on the GPU. Morton codes, or Z-order curve values, project multi-dimensional data into a single dimension while preserving spatial locality. By quantizing bounding box centers into integer coordinates and interleaving their bits using bitwise operations, compute shaders can generate sorted linear sequences of primitives essential for fast, parallel LBVH generation.

Understanding Morton Codes in LBVHs

Building a Linear BVH requires ordering geometric primitives along a space-filling curve. The 3D Morton code achieves this by interleaving the binary representations of a point's \(x\), \(y\), and \(z\) coordinates.

When applied to 32-bit unsigned integers in GLSL, each axis is typically quantized into a 10-bit integer (representing values from \(0\) to \(1023\)). Interleaving three 10-bit values produces a 30-bit Morton code, which fits cleanly inside a single 32-bit uint. Sorting primitives by this 30-bit key clusters spatially adjacent geometry together, allowing the hierarchy to be constructed in parallel using algorithms such as Karras' LBVH builder.

Quantizing Coordinates to Integers

Before interleaving bits, floating-point positions must be normalized relative to the scene's axis-aligned bounding box (AABB) and scaled to a 10-bit integer range (\([0, 2^{10} - 1]\)):

uvec3 quantizePosition(vec3 point, vec3 sceneMin, vec3 sceneMax) {
    // Normalize coordinates to [0.0, 1.0]
    vec3 normalized = (point - sceneMin) / (sceneMax - sceneMin);
    
    // Clamp to ensure values remain within range
    normalized = clamp(normalized, 0.0, 1.0);
    
    // Scale to 10-bit integer range [0, 1023]
    return uvec3(normalized * 1023.0);
}

Expanding Bits via Dilation

To interleave the bits of \(x\), \(y\), and \(z\), each coordinate's bits must first be spread out so that every bit is separated by two zeros. A 10-bit integer can be dilated into a 30-bit integer using a sequence of bit shifts and masks:

uint expandBits(uint v) {
    v = (v * 0x00010001u) & 0xFF0000FFu;
    v = (v * 0x00000101u) & 0x0F00F00Fu;
    v = (v * 0x00000011u) & 0xC30C30C3u;
    v = (v * 0x00000005u) & 0x49249249u;
    return v;
}

Alternatively, standard shift-or bit manipulation can be used without multiplication:

uint expandBits(uint v) {
    v = (v | (v << 16)) & 0x030000FFu;
    v = (v | (v <<  8)) & 0x0300F00Fu;
    v = (v | (v <<  4)) & 0x030C30C3u;
    v = (v | (v <<  2)) & 0x09249249u;
    return v;
}

Both methods insert two zeros between each of the lowest 10 bits of v.

Interleaving Coordinates into a 3D Morton Code

Once each axis coordinate is dilated, the final 30-bit Morton code is generated by shifting the dilated components and combining them with bitwise OR operations:

uint calculateMorton3D(uvec3 quantizedPos) {
    uint xx = expandBits(quantizedPos.x);
    uint yy = expandBits(quantizedPos.y);
    uint zz = expandBits(quantizedPos.z);
    
    return (xx << 2) | (yy << 1) | zz;
}

In this encoding, bit \(3i + 2\) belongs to the X coordinate, bit \(3i + 1\) belongs to the Y coordinate, and bit \(3i\) belongs to the Z coordinate for any bit index \(i \in [0, 9]\).

Complete GLSL Compute Implementation

The full pipeline integrates coordinate normalization, quantization, dilation, and interleaving within a compute shader stage:

#version 450

layout(local_size_x = 64) in;

struct Primitive {
    vec3 minBound;
    vec3 maxBound;
};

layout(std430, binding = 0) readonly buffer InputPrimitives {
    Primitive primitives[];
};

layout(std430, binding = 1) writeonly buffer OutputMortonCodes {
    uint mortonKeys[];
};

uniform vec3 u_SceneMin;
uniform vec3 u_SceneMax;
uniform uint u_PrimitiveCount;

uint expandBits(uint v) {
    v = (v | (v << 16)) & 0x030000FFu;
    v = (v | (v <<  8)) & 0x0300F00Fu;
    v = (v | (v <<  4)) & 0x030C30C3u;
    v = (v | (v <<  2)) & 0x09249249u;
    return v;
}

uint morton3D(vec3 normalizedPos) {
    normalizedPos = clamp(normalizedPos, 0.0, 1.0);
    uvec3 i = uvec3(normalizedPos * 1023.0);
    return (expandBits(i.x) << 2) | (expandBits(i.y) << 1) | expandBits(i.z);
}

void main() {
    uint index = gl_GlobalInvocationID.x;
    if (index >= u_PrimitiveCount) return;

    vec3 centroid = (primitives[index].minBound + primitives[index].maxBound) * 0.5;
    vec3 normalizedCentroid = (centroid - u_SceneMin) / (u_SceneMax - u_SceneMin);

    mortonKeys[index] = morton3D(normalizedCentroid);
}

The resulting array of Morton codes can then be passed to a GPU sorting pass (such as a parallel Radix Sort) to finalize primitive ordering for LBVH hierarchy construction.