De Casteljau Cubic Bezier Algorithm in GLSL

This article provides a concise guide to implementing De Casteljau's algorithm for evaluating cubic Bézier curves using the OpenGL Shading Language (GLSL). You will explore the mathematical foundation of recursive linear interpolation, examine an optimized GLSL function utilizing the native mix operation, and understand the practical benefits of this approach for GPU-based rendering and procedural geometry.

Understanding De Casteljau's Algorithm

De Casteljau's algorithm evaluates a Bézier curve at parameter \(t \in [0, 1]\) through repeated linear interpolation between control points. A cubic Bézier curve is defined by four control points: \(P_0, P_1, P_2,\) and \(P_3\).

The evaluation proceeds through three sequential reduction steps:

  1. First Reduction (Order 3 to 2): Interpolate between adjacent control points to produce three intermediate points. \[Q_0 = (1 - t)P_0 + tP_1\] \[Q_1 = (1 - t)P_1 + tP_2\] \[Q_2 = (1 - t)P_2 + tP_3\]

  2. Second Reduction (Order 2 to 1): Interpolate between the resulting three points to produce two intermediate points. \[R_0 = (1 - t)Q_0 + tQ_1\] \[R_1 = (1 - t)Q_1 + tQ_2\]

  3. Final Reduction (Order 1 to 0): Interpolate between the final two points to obtain the position on the curve. \[C(t) = (1 - t)R_0 + tR_1\]

GLSL Implementation

In GLSL, linear interpolation is computed efficiently using the built-in mix(a, b, t) function, which maps directly to hardware-level Multiply-Add (MAD) or Fused Multiply-Add (FMA) instructions.

Here is the complete implementation for 3D coordinates (vec3). This function can be adapted to vec2 or vec4 by changing the vector types:

// Evaluates a cubic Bezier curve at parameter t using De Casteljau's algorithm
vec3 evaluateCubicBezier(vec3 p0, vec3 p1, vec3 p2, vec3 p3, float t) {
    // First stage of interpolation (4 points -> 3 points)
    vec3 q0 = mix(p0, p1, t);
    vec3 q1 = mix(p1, p2, t);
    vec3 q2 = mix(p2, p3, t);

    // Second stage of interpolation (3 points -> 2 points)
    vec3 r0 = mix(q0, q1, t);
    vec3 r1 = mix(q1, q2, t);

    // Final interpolation (2 points -> 1 point)
    return mix(r0, r1, t);
}

First Derivative (Tangent Calculation)

To compute the tangent vector along the curve for normal calculations or mesh extrusion, evaluate the difference between the two points from the second-to-last reduction step (\(R_0\) and \(R_1\)):

// Computes the position and tangent vector at parameter t
void evaluateCubicBezierWithTangent(
    vec3 p0, vec3 p1, vec3 p2, vec3 p3, float t, 
    out vec3 position, out vec3 tangent
) {
    vec3 q0 = mix(p0, p1, t);
    vec3 q1 = mix(p1, p2, t);
    vec3 q2 = mix(p2, p3, t);

    vec3 r0 = mix(q0, q1, t);
    vec3 r1 = mix(q1, q2, t);

    position = mix(r0, r1, t);
    tangent = normalize(3.0 * (r1 - r0));
}

Numerical Stability vs. Direct Polynomial Form

Cubic Bézier curves can also be calculated via the standard Bernstein polynomial formula:

\[B(t) = (1-t)^3 P_0 + 3(1-t)^2 t P_1 + 3(1-t) t^2 P_2 + t^3 P_3\]

While the direct polynomial method requires fewer operations overall, De Casteljau's algorithm offers distinct advantages in shader environments:

  • Numerical Precision: Repeated interpolation keeps intermediate values strictly within the convex hull of the control points, minimizing floating-point precision loss when \(t\) approaches \(0.0\) or \(1.0\).
  • Intermediate Subdivision: The intermediate variables (q0, q1, q2, r0, r1) represent the exact split points needed to subdivide the curve into two smaller Bézier segments directly in a geometry or compute shader.