3D Tri-Quadratic B-Spline Interpolation in GLSL

This article provides the mathematical foundation and GLSL implementation details for tri-quadratic B-spline interpolation over 3D voxel grids. Unlike standard trilinear interpolation, which exhibits \(C^0\) continuity and visible derivative discontinuities along voxel boundaries, tri-quadratic B-spline interpolation yields a \(C^1\)-continuous scalar field. This creates smoother gradients and eliminates grid artifacts in volume rendering and implicit surface extraction.

Mathematical Formulation

A uniform quadratic B-spline basis function is defined over three intervals across a support of three units. When evaluating an arbitrary continuous coordinate \(\mathbf{x} = (x, y, z)\) within a discrete voxel grid, the coordinate is first shifted by \(-0.5\) to center the basis functions on grid points.

1. Discrete Indexing and Local Parameterization

Given a continuous coordinate \(x\) in voxel units:

\[i = \lfloor x - 0.5 \rfloor\] \[f = x - 0.5 - i \quad \text{where } f \in [0, 1)\]

The continuous coordinate is reconstructed using the discrete index \(i\) and three consecutive samples along the axis: \(i\), \(i+1\), and \(i+2\).

2. One-Dimensional Basis Functions

The 1D uniform quadratic B-spline weights for the local fractional offset \(f \in [0, 1)\) are defined as:

\[w_0(f) = \frac{1}{2}(1 - f)^2\] \[w_1(f) = -f^2 + f + \frac{1}{2} = \frac{3}{4} - \left(f - \frac{1}{2}\right)^2\] \[w_2(f) = \frac{1}{2}f^2\]

These weights satisfy the partition of unity property:

\[\sum_{k=0}^{2} w_k(f) = w_0(f) + w_1(f) + w_2(f) = 1\]

3. 3D Tensor Product Expansion

Because the 3D B-spline is separable, the tri-quadratic reconstruction \(S(x, y, z)\) is the tensor product of the 1D basis functions across all three spatial dimensions:

\[S(x, y, z) = \sum_{k=0}^{2} \sum_{j=0}^{2} \sum_{i=0}^{2} w_i(f_x) \, w_j(f_y) \, w_k(f_z) \, c_{(i_0 + i, \, j_0 + j, \, k_0 + k)}\]

where:

\(i_0 = \lfloor x - 0.5 \rfloor\), \(j_0 = \lfloor y - 0.5 \rfloor\), \(k_0 = \lfloor z - 0.5 \rfloor\)
\(f_x, f_y, f_z\) are the fractional offsets for each axis
\(c_{(u, v, w)}\) are the voxel grid coefficient values (or raw voxel samples if approximating)

Hardware-Accelerated GLSL Evaluation (8-Tap Filter)

A direct implementation of tri-quadratic B-spline evaluation requires \(3 \times 3 \times 3 = 27\) texture lookups. By leveraging the GPU's native trilinear filtering hardware, pairs of adjacent B-spline weights can be combined into single texture fetches.

For each dimension, the three weights \(w_0, w_1, w_2\) are combined into two evaluation points:

Sample 0 fetches coordinate offset \(\lambda_0 = 0\) with weight \(W_0 = w_0\).
Samples 1 and 2 are combined into an effective single sample with combined weight \(W_1 = w_1 + w_2\) at the interpolated offset \(d_1 = \frac{w_2}{w_1 + w_2}\).

This reduces the required texture fetches from 27 to 8 (\(2 \times 2 \times 2\)).

GLSL Implementation

// Computes 1D quadratic B-spline weights and hardware texture offsets
void getBsplineWeights(float x, out vec2 offsets, out vec2 weights) {
    float i = floor(x - 0.5);
    float f = (x - 0.5) - i;

    float w0 = 0.5 * (1.0 - f) * (1.0 - f);
    float w1 = 0.5 + f * (1.0 - f);
    float w2 = 0.5 * f * f;

    // Combine weights w1 and w2 for hardware linear sampling
    weights.x = w0;
    weights.y = w1 + w2;

    // Coordinate offsets relative to floor(x - 0.5)
    offsets.x = i;
    offsets.y = i + 1.0 + (w2 / weights.y);
}

// Samples a 3D texture using tri-quadratic B-spline interpolation via 8 texture taps
vec4 texture3D_TriQuadratic(sampler3D tex, vec3 coord, vec3 texSize) {
    vec3 voxelCoord = coord * texSize;

    vec2 ox, wx;
    vec2 oy, wy;
    vec2 oz, wz;

    getBsplineWeights(voxelCoord.x, ox, wx);
    getBsplineWeights(voxelCoord.y, oy, wy);
    getBsplineWeights(voxelCoord.z, oz, wz);

    // Convert voxel-space tap locations back to normalized [0, 1] texture coordinates
    vec3 invTexSize = 1.0 / texSize;
    vec3 p000 = vec3(ox.x + 0.5, oy.x + 0.5, oz.x + 0.5) * invTexSize;
    vec3 p100 = vec3(ox.y + 0.5, oy.x + 0.5, oz.x + 0.5) * invTexSize;
    vec3 p010 = vec3(ox.x + 0.5, oy.y + 0.5, oz.x + 0.5) * invTexSize;
    vec3 p110 = vec3(ox.y + 0.5, oy.y + 0.5, oz.x + 0.5) * invTexSize;
    vec3 p001 = vec3(ox.x + 0.5, oy.x + 0.5, oz.y + 0.5) * invTexSize;
    vec3 p101 = vec3(ox.y + 0.5, oy.x + 0.5, oz.y + 0.5) * invTexSize;
    vec3 p011 = vec3(ox.x + 0.5, oy.y + 0.5, oz.y + 0.5) * invTexSize;
    vec3 p111 = vec3(ox.y + 0.5, oy.y + 0.5, oz.y + 0.5) * invTexSize;

    // 8-tap trilinear fetch accumulation
    vec4 result = 
        texture(tex, p000) * (wx.x * wy.x * wz.x) +
        texture(tex, p100) * (wx.y * wy.x * wz.x) +
        texture(tex, p010) * (wx.x * wy.y * wz.x) +
        texture(tex, p110) * (wx.y * wy.y * wz.x) +
        texture(tex, p001) * (wx.x * wy.x * wz.y) +
        texture(tex, p101) * (wx.y * wy.x * wz.y) +
        texture(tex, p011) * (wx.x * wy.y * wz.y) +
        texture(tex, p111) * (wx.y * wy.y * wz.y);

    return result;
}

Interpolation vs. Approximation

Uniform quadratic B-splines are approximating filters by default; evaluating the spline directly on the raw data samples produces low-pass smoothing where \(S(i, j, k) \neq P_{i, j, k}\).

To make the formulation strictly interpolating—such that the reconstructed field matches the original voxel samples exactly at integer coordinates—the original voxel array must undergo an offline 3D recursive prefiltering step. This prefilter applies an inverse B-spline filter along each axis to convert raw data samples into B-spline coefficient controls prior to upload to the GPU 3D texture.