Bicubic B-Spline Texture Interpolation in GLSL

This article explains how to implement bicubic B-spline texture interpolation manually within a GLSL fragment shader. Standard hardware bilinear filtering only offers \(C^0\) continuity, which introduces noticeable grid artifacts and sharp gradient shifts when textures are magnified. Bicubic B-spline filtering reconstructs continuous \(C^2\) surfaces by taking a weighted average of a \(4 \times 4\) neighborhood of texels using cubic polynomials. While a naive implementation requires 16 texture taps, you can leverage GPU bilinear texture filtering hardware to evaluate the filter using only 4 optimized texture fetches.

The Cubic B-Spline Basis Functions

Uniform cubic B-spline interpolation over an interval \(t \in [0, 1]\) is defined by four polynomial basis weights:

\[\begin{aligned} w_0(t) &= \frac{1}{6}(1 - t)^3 \\ w_1(t) &= \frac{1}{6}(3t^3 - 6t^2 + 4) \\ w_2(t) &= \frac{1}{6}(-3t^3 + 3t^2 + 3t + 1) \\ w_3(t) &= \frac{1}{6}t^3 \end{aligned}\]

Because the partition of unity holds (\(w_0 + w_1 + w_2 + w_3 = 1\)), the reconstructed curve is smooth and normalized across all texels.

Reducing 16 Texture Fetches to 4

A direct 2D cubic filter samples a \(4 \times 4\) grid of texels, requiring 16 dependent texture calls. However, modern GPUs perform hardware bilinear filtering between adjacent texels.

By combining pairs of 1D weights \((w_0, w_1)\) and \((w_2, w_3)\), you can collapse the 4-tap 1D filter into 2 bilinear taps:

\[\begin{aligned} g_0 &= w_0 + w_1, \quad &h_0 &= \frac{w_1}{w_0 + w_1} \\ g_1 &= w_2 + w_3, \quad &h_1 &= \frac{w_3}{w_2 + w_3} \end{aligned}\]

Using \(h_0\) and \(h_1\) as sampling offsets relative to the integer texel grid, the GPU's bilinear sampler will blend the two adjacent texels according to their exact B-spline weight ratios. Extending this across both axes reduces the required samples from 16 to 4 bilinear fetches.

GLSL Implementation

Ensure the texture sampler is configured with GL_LINEAR filtering enabled so that sub-texel coordinates interpolate correctly.

#version 330 core

in vec2 v_TexCoords;
out vec4 FragColor;

uniform sampler2D u_Texture;
uniform vec2 u_TextureResolution;

vec4 sampleBicubicBSpline(sampler2D tex, vec2 uv, vec2 resolution) {
    // Transform coordinates from [0, 1] UV to texel units
    vec2 tc = uv * resolution - 0.5;
    vec2 index = floor(tc);
    vec2 f = fract(tc);

    // Evaluate 1D cubic B-spline weights
    vec2 f2 = f * f;
    vec2 f3 = f2 * f;

    vec2 w0 = (1.0 - 3.0 * f + 3.0 * f2 - f3) / 6.0;
    vec2 w1 = (4.0 - 6.0 * f2 + 3.0 * f3) / 6.0;
    vec2 w2 = (1.0 + 3.0 * f + 3.0 * f2 - 3.0 * f3) / 6.0;
    vec2 w3 = f3 / 6.0;

    // Combine weights to find bilinear sampling offsets
    vec2 g0 = w0 + w1;
    vec2 g1 = w2 + w3;

    vec2 h0 = w1 / g0;
    vec2 h1 = w3 / g1;

    // Compute sample locations in UV space
    vec2 tap0 = (index - 0.5 + h0) / resolution;
    vec2 tap1 = (index + 1.5 + h1) / resolution;

    // Perform the 4 bilinear fetches
    vec4 tex00 = texture(tex, vec2(tap0.x, tap0.y));
    vec4 tex10 = texture(tex, vec2(tap1.x, tap0.y));
    vec4 tex01 = texture(tex, vec2(tap0.x, tap1.y));
    vec4 tex11 = texture(tex, vec2(tap1.x, tap1.y));

    // Combine horizontal bilinear pairs
    vec4 row0 = mix(tex10, tex00, g0.x);
    vec4 row1 = mix(tex11, tex01, g0.x);

    // Combine vertical results
    return mix(row1, row0, g0.y);
}

void main() {
    FragColor = sampleBicubicBSpline(u_Texture, v_TexCoords, u_TextureResolution);
}

How the Coordinate Math Works

  1. Texel Offset: The UV coordinates are multiplied by the texture dimensions and shifted by -0.5 to align UV sample points with the centers of texels.
  2. Fractional Split: floor(tc) locates the top-left texel index of the central \(2 \times 2\) block, while fract(tc) provides the sub-texel interpolation phase \(f \in [0, 1]\).
  3. Hardware Blending: tap0 samples between texels \((i - 1)\) and \(i\), while tap1 samples between texels \((i + 1)\) and \((i + 2)\). The relative offsets \(h_0\) and \(h_1\) shift the hardware tap to reproduce the ratio of \(w_0\) to \(w_1\) and \(w_2\) to \(w_3\).
  4. Final Linear Blend: The mix calls evaluate the outer weights \(g_0\) and \(g_1\). Because \(g_0 + g_1 = 1\), using mix(b, a, g0) executes a * g0 + b * (1.0 - g0), producing the correct composite cubic reconstruction.