Catmull-Rom Bicubic Filtering with 4 Samples in GLSL

Catmull-Rom bicubic texture filtering produces sharp, artifact-free image magnification by fitting cubic splines over neighboring texels. While a naive bicubic evaluation requires 16 individual texel lookups on a 4x4 grid, you can leverage the GPU's native bilinear sampling hardware to reduce the operation to just 4 bilinear samples. This article explains the mathematical reduction and provides a practical GLSL implementation for real-time rendering pipelines.

The Mathematical Principle

Standard 2D bicubic filtering samples a \(4 \times 4\) neighborhood of discrete texels around a continuous coordinate. The 1D Catmull-Rom spline evaluated at fractional offset \(f \in [0, 1)\) between texel \(1\) and \(2\) uses four weights (\(w_0, w_1, w_2, w_3\)) applied to texels at positions \(-1, 0, 1,\) and \(2\):

\[w_0(f) = -0.5f^3 + f^2 - 0.5f\] \[w_1(f) = 1.5f^3 - 2.5f^2 + 1.0\] \[w_2(f) = -1.5f^3 + 2.0f^2 + 0.5f\] \[w_3(f) = 0.5f^3 - 0.5f^2\]

These weights satisfy the partition of unity (\(\sum w_i = 1\)).

Hardware bilinear filtering computes an interpolated value between two adjacent texels \(T_a\) and \(T_b\) at a continuous sub-texel coordinate \(a + \alpha\):

\[\text{sample}(a + \alpha) = (1 - \alpha)T_a + \alpha T_b\]

By pairing the Catmull-Rom weights into groups of two—\((w_0, w_1)\) and \((w_2, w_3)\)—we can solve for combined weights and sub-texel sampling offsets:

  1. First Pair (\(T_0\) and \(T_1\)):
    • Combined Weight: \(W_{01} = w_0 + w_1\)
    • Normalized Offset: \(\alpha = \frac{w_1}{w_0 + w_1}\)
    • Tap Position: \(P_{01} = \text{texelIndex}_0 - 0.5 + \alpha = \text{texelIndex}_0 - 0.5 + \frac{w_1}{W_{01}}\)
  2. Second Pair (\(T_2\) and \(T_3\)):
    • Combined Weight: \(W_{23} = w_2 + w_3\)
    • Normalized Offset: \(\beta = \frac{w_3}{w_2 + w_3}\)
    • Tap Position: \(P_{23} = \text{texelIndex}_1 + 0.5 + \beta = \text{texelIndex}_1 + 0.5 + \frac{w_3}{W_{23}}\)

Applying this reduction to both the X and Y dimensions yields a 2D tensor product: 2 bilinear taps along X multiplied by 2 bilinear taps along Y, resulting in 4 bilinear texture samples instead of 16 nearest-neighbor samples.

Algorithm Steps

  1. Transform normalized UV coordinates into continuous texel coordinates offset by \(0.5\) texels so that texel centers fall on integer values.
  2. Separate the coordinates into integer indices (floor) and fractional offsets (fract).
  3. Compute the 1D Catmull-Rom polynomial weights for both \(x\) and \(y\).
  4. Calculate the combined weights and bilinear sample UV coordinates for each axis.
  5. Sample the texture at the four coordinate combinations.
  6. Compute the weighted sum of the four samples to produce the filtered output.

GLSL Implementation

vec4 SampleTextureCatmullRom(sampler2D tex, vec2 uv, vec2 texSize) {
    vec2 invTexSize = 1.0 / texSize;
    
    // Scale UV to texel coordinates centered on texel pixels
    vec2 tc = uv * texSize - 0.5;
    vec2 f = fract(tc);
    vec2 index = floor(tc);

    // Evaluate 1D Catmull-Rom cubic weights
    vec2 f2 = f * f;
    vec2 f3 = f2 * f;

    vec2 w0 = -0.5 * f3 + f2 - 0.5 * f;
    vec2 w1 =  1.5 * f3 - 2.5 * f2 + 1.0;
    vec2 w2 = -1.5 * f3 + 2.0 * f2 + 0.5 * f;
    vec2 w3 =  0.5 * f3 - 0.5 * f2;

    // Combine adjacent pairs into bilinear weights
    vec2 w01 = w0 + w1;
    vec2 w23 = w2 + w3;

    // Compute bilinear sample locations in normalized texture coordinates
    vec2 pos01 = (index - 0.5 + (w1 / w01)) * invTexSize;
    vec2 pos23 = (index + 1.5 + (w3 / w23)) * invTexSize;

    // Perform 4 bilinear texture fetches
    vec4 s0 = texture(tex, vec2(pos01.x, pos01.y));
    vec4 s1 = texture(tex, vec2(pos23.x, pos01.y));
    vec4 s2 = texture(tex, vec2(pos01.x, pos23.y));
    vec4 s3 = texture(tex, vec2(pos23.x, pos23.y));

    // Combine horizontal pairs using X weights
    vec4 row0 = mix(s0, s1, w23.x / (w01.x + w23.x));
    vec4 row1 = mix(s2, s3, w23.x / (w01.x + w23.x));

    // Combine vertical results using Y weights
    return mix(row0, row1, w23.y / (w01.y + w23.y));
}

Considerations for Implementation

  • Texture Filtering Mode: The texture sampler bound in OpenGL/Vulkan must be set to GL_LINEAR (minFilter and magFilter) for hardware bilinear filtering to operate correctly.
  • Negative Weights: Catmull-Rom splines contain negative lobes that cause minor undershoot or overshoot near high-contrast edges. If your target requires values within the strict \([0, 1]\) range (e.g., standard LDR output), clamp the final result using clamp(color, 0.0, 1.0).
  • Texel Boundary Clamping: At texture boundaries, UV calculations can sample beyond texel limits. Use GL_CLAMP_TO_EDGE on the sampler to avoid edge bleeding.