GLSL 3D LUT Trilinear Interpolation Algorithm

This article provides a technical overview and implementation guide for trilinear interpolation across 3D color lookup tables (LUTs) in GLSL. It explains the underlying mathematical formulation for sampling within a discretized 3D lattice, covers manual 8-point interpolation using direct texel fetches, and demonstrates how to apply this algorithm to both native 3D textures and 2D unwrapped LUT strips.

Mathematical Concept of Trilinear Interpolation

A 3D LUT discretizes an RGB color space into an \(N \times N \times N\) cube of output colors. When an input color does not fall directly on an integer grid point, trilinear interpolation approximates the value by calculating a weighted average of the eight nearest surrounding lattice vertices (corners of the surrounding voxel).

Given an input color \(C = (R, G, B)\) normalized to \([0.0, 1.0]\) and a LUT size \(N\):

  1. Scale Coordinates: Convert normalized coordinates to grid space coordinates: \[\mathbf{P} = C \times (N - 1)\]
  2. Determine Bounds: Compute lower indices \(\mathbf{P}_0 = \lfloor \mathbf{P} \rfloor\) and upper indices \(\mathbf{P}_1 = \min(\mathbf{P}_0 + 1, N - 1)\).
  3. Calculate Fractional Weights: Compute the fractional distance within the cell: \[\mathbf{f} = \mathbf{P} - \mathbf{P}_0\]
  4. Interpolate Across Dimensions:
    • Linearly interpolate along the X-axis (Red) to produce 4 intermediate values.
    • Linearly interpolate those 4 values along the Y-axis (Green) to produce 2 values.
    • Linearly interpolate the remaining 2 values along the Z-axis (Blue) to produce the final color.

GLSL Implementation: Manual 8-Point Interpolation (sampler3D)

While setting texture filtering to GL_LINEAR on a native 3D texture delegates trilinear filtering directly to the GPU hardware, writing the algorithm manually is required when using texelFetch or custom precision routines.

#version 330 core

uniform sampler3D u_LutTexture;
uniform float u_LutSize; // e.g., 32.0, 64.0

vec3 sampleLutTrilinear(vec3 color) {
    // Clamp input to [0.0, 1.0] to prevent sampling outside the lattice
    color = clamp(color, 0.0, 1.0);

    // Scale to the lattice domain [0, LUT_SIZE - 1]
    vec3 coord = color * (u_LutSize - 1.0);

    // Get the base voxel corner and the fractional offset
    ivec3 i0 = ivec3(floor(coord));
    ivec3 i1 = min(i0 + ivec3(1), ivec3(int(u_LutSize) - 1));
    vec3 f = fract(coord);

    // Fetch the 8 surrounding lattice vertices
    vec3 c000 = texelFetch(u_LutTexture, ivec3(i0.x, i0.y, i0.z), 0).rgb;
    vec3 c100 = texelFetch(u_LutTexture, ivec3(i1.x, i0.y, i0.z), 0).rgb;
    vec3 c010 = texelFetch(u_LutTexture, ivec3(i0.x, i1.y, i0.z), 0).rgb;
    vec3 c110 = texelFetch(u_LutTexture, ivec3(i1.x, i1.y, i0.z), 0).rgb;
    vec3 c001 = texelFetch(u_LutTexture, ivec3(i0.x, i0.y, i1.z), 0).rgb;
    vec3 c101 = texelFetch(u_LutTexture, ivec3(i1.x, i0.y, i1.z), 0).rgb;
    vec3 c011 = texelFetch(u_LutTexture, ivec3(i0.x, i1.y, i1.z), 0).rgb;
    vec3 c111 = texelFetch(u_LutTexture, ivec3(i1.x, i1.y, i1.z), 0).rgb;

    // Interpolate along the X-axis (Red)
    vec3 c00 = mix(c000, c100, f.x);
    vec3 c10 = mix(c010, c110, f.x);
    vec3 c01 = mix(c001, c101, f.x);
    vec3 c11 = mix(c011, c111, f.x);

    // Interpolate along the Y-axis (Green)
    vec3 c0 = mix(c00, c10, f.y);
    vec3 c1 = mix(c01, c11, f.y);

    // Interpolate along the Z-axis (Blue)
    return mix(c0, c1, f.z);
}

GLSL Implementation: 2D Texture Atlas (Horizontal Strip)

In environments without native 3D texture support, 3D LUTs are typically stored as a 2D strip of size \((N \times N) \times N\). Trilinear interpolation in this configuration combines hardware bilinear interpolation across the 2D tiles with a manual mix() between the two nearest Z-slices.

#version 330 core

uniform sampler2D u_LutAtlas;
uniform float u_LutSize; // e.g., 32.0 or 64.0

vec3 sample2DStripLutTrilinear(vec3 color) {
    color = clamp(color, 0.0, 1.0);

    float maxCoord = u_LutSize - 1.0;
    float blue = color.b * maxCoord;

    // Determine neighboring blue slices
    float slice0 = floor(blue);
    float slice1 = min(slice0 + 1.0, maxCoord);
    float fz = fract(blue);

    // Standard half-texel offset for accurate UV sampling
    vec2 halfTexel = vec2(0.5 / (u_LutSize * u_LutSize), 0.5 / u_LutSize);

    // Compute base UV for slice 0 and slice 1
    float uOffset0 = slice0 * u_LutSize;
    float uOffset1 = slice1 * u_LutSize;

    vec2 uv0 = vec2(
        (uOffset0 + 0.5 + color.r * maxCoord) / (u_LutSize * u_LutSize),
        (0.5 + color.g * maxCoord) / u_LutSize
    );

    vec2 uv1 = vec2(
        (uOffset1 + 0.5 + color.r * maxCoord) / (u_LutSize * u_LutSize),
        (0.5 + color.g * maxCoord) / u_LutSize
    );

    // Hardware bilinear filtering performs X and Y interpolation
    vec3 sample0 = texture(u_LutAtlas, uv0).rgb;
    vec3 sample1 = texture(u_LutAtlas, uv1).rgb;

    // Linear interpolation along the Z axis completes trilinear filtering
    return mix(sample0, sample1, fz);
}

Hardware Half-Texel Compensation

When converting continuous \([0.0, 1.0]\) coordinates to normalized texture coordinates, sampling exact cell centers requires adding an offset of half a texel (\(\frac{0.5}{N}\)):

\[\text{UV} = \frac{C \times (N - 1) + 0.5}{N}\]

Failing to apply this offset compresses the dynamic range, causing boundary clamping artifacts where the outer edges of the color gamut sample toward adjacent texture boundaries.