Catmull-Rom Spline Interpolation in GLSL

This article explains how to evaluate a cubic Catmull-Rom spline across four control points directly inside an OpenGL Shading Language (GLSL) shader. Catmull-Rom splines provide smooth, \(C^1\)-continuous curves that pass through their internal control points, making them ideal for camera paths, particle trajectories, and smooth texture sampling. Below, you will find the mathematical derivation of the spline and an optimized GLSL implementation using standard vector operations.

Mathematical Formulation

A cubic Catmull-Rom spline evaluates a point along a segment bounded by control points \(P_1\) and \(P_2\), with control points \(P_0\) and \(P_3\) defining the entry and exit tangents. The interpolation factor \(t\) ranges from \(0.0\) to \(1.0\), where \(t = 0.0\) returns \(P_1\) and \(t = 1.0\) returns \(P_2\).

The standard uniform Catmull-Rom spline can be expressed as a matrix multiplication:

\[ P(t) = \frac{1}{2} \begin{bmatrix} 1 & t & t^2 & t^3 \end{bmatrix} \begin{bmatrix} 0 & 2 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 2 & -5 & 4 & -1 \\ -1 & 3 & -3 & 1 \end{bmatrix} \begin{bmatrix} P_0 \\ P_1 \\ P_2 \\ P_3 \end{bmatrix} \]

Expanding this into polynomial coefficients yields:

  • \(a = -P_0 + 3P_1 - 3P_2 + P_3\)
  • \(b = 2P_0 - 5P_1 + 4P_2 - P_3\)
  • \(c = -P_0 + P_2\)
  • \(d = 2P_1\)

The interpolated position is computed as:

\[ P(t) = 0.5 \cdot (d + c \cdot t + b \cdot t^2 + a \cdot t^3) \]

GLSL Implementation

Using Horner's rule, the polynomial evaluation requires fewer multiplications, minimizing instruction count on the GPU.

// Evaluates a cubic Catmull-Rom spline between p1 and p2 using p0 and p3 as tangents.
// Parameter t should be normalized in the range [0.0, 1.0].
vec3 catmullRom(vec3 p0, vec3 p1, vec3 p2, vec3 p3, float t) {
    vec3 a = -p0 + 3.0 * p1 - 3.0 * p2 + p3;
    vec3 b = 2.0 * p0 - 5.0 * p1 + 4.0 * p2 - p3;
    vec3 c = -p0 + p2;
    vec3 d = 2.0 * p1;

    return 0.5 * (((a * t + b) * t + c) * t + d);
}

This function can be overloaded for float, vec2, or vec4 by changing the component types accordingly, as GLSL supports component-wise arithmetic on all vector types.

Boundary Handling and Usage

When using Catmull-Rom splines across a longer sequence of points:

  1. Segment Selection: To interpolate between point \(i\) and point \(i+1\), pass points \(i-1\), \(i\), \(i+1\), and \(i+2\) as \(P_0, P_1, P_2, P_3\).
  2. Endpoints: For the first segment, where \(P_{-1}\) does not exist, approximate \(P_0\) using linear extrapolation: \(P_0 = 2P_1 - P_2\). Similarly, for the final segment, approximate the trailing point as \(P_3 = 2P_2 - P_1\).
  3. Parameter Normalization: Ensure \(t\) is clamped using clamp(t, 0.0, 1.0) to avoid overshooting outside the segment bounds.