Bézier curves

Deep dive

How it works

In: control points P0..P3 (2-D or 3-D) and a parameter t in [0,1], or a tolerance for flattening · Out: the point B(t), the two split halves, a bounding box, or a polyline · State: none · Cost: 6 lerps per cubic evaluation; flattening n = ceil(sqrt(6M / (8 tol))) segments

Control points P0..P3, parameter t in [0,1].

1. Evaluate (de Casteljau). a=lerp(P0,P1,t), b=lerp(P1,P2,t), c=lerp(P2,P3,t); d=lerp(a,b,t), e=lerp(b,c,t); B(t)=lerp(d,e,t).
2. Split at t. The first half has control points P0, a, d, B(t), the second B(t), e, c, P3. No approximation.
3. Bounding box. Per axis, B'(t)/3 is a quadratic a t² + b t + c with a=-p0+3p1-3p2+p3, b=2(p0-2p1+p2), c=p1-p0. Solve it, keep roots in (0,1), evaluate there, and take min/max together with the endpoints.
4. Flatten. A uniform polyline with n segments deviates from the curve by at most max|B''| / (8 n²), and |B''| ≤ 6M where M = max(|P0-2P1+P2|, |P1-2P2+P3|). So n = ceil(sqrt(6M / (8 tol))), called Wang's formula.

level 0:  P0   P1   P2   P3
level 1:    a    b    c        a,b,c = lerp of neighbours
level 2:      d    e           d = lerp(a,b), e = lerp(b,c)
level 3:        m              m = B(t); halves: [P0 a d m] [m e c P3]

Code

#[derive(Clone, Copy, Debug, PartialEq)]
pub struct P { pub x: f64, pub y: f64 }

fn lerp(a: P, b: P, t: f64) -> P { P { x: a.x + (b.x - a.x) * t, y: a.y + (b.y - a.y) * t } }
fn dist(a: P, b: P) -> f64 { (a.x - b.x).hypot(a.y - b.y) }

#[derive(Clone, Copy, Debug)]
pub struct Cubic(pub [P; 4]);

impl Cubic {
    /// de Casteljau: three rounds of lerp; numerically stable for t in [0,1].
    pub fn eval(&self, t: f64) -> P { self.split(t).0 .0[3] }

    /// Split at t into two cubics that together trace exactly the same curve.
    pub fn split(&self, t: f64) -> (Cubic, Cubic) {
        let [p0, p1, p2, p3] = self.0;
        let (a, b, c) = (lerp(p0, p1, t), lerp(p1, p2, t), lerp(p2, p3, t));
        let (d, e) = (lerp(a, b, t), lerp(b, c, t));
        let m = lerp(d, e, t); // the point on the curve
        (Cubic([p0, a, d, m]), Cubic([m, e, c, p3]))
    }

    /// Tight axis-aligned box: endpoints plus roots of B'(t) = 0 in (0,1) per axis.
    pub fn bbox(&self) -> (P, P) {
        let mut lo = P { x: self.0[0].x.min(self.0[3].x), y: self.0[0].y.min(self.0[3].y) };
        let mut hi = P { x: self.0[0].x.max(self.0[3].x), y: self.0[0].y.max(self.0[3].y) };
        for axis in 0..2 {
            let v = |p: P| if axis == 0 { p.x } else { p.y };
            let [p0, p1, p2, p3] = self.0.map(v);
            // B'(t)/3 = a t^2 + b t + c
            let (a, b, c) = (-p0 + 3. * p1 - 3. * p2 + p3, 2. * (p0 - 2. * p1 + p2), p1 - p0);
            let mut roots = [f64::NAN; 2];
            if a.abs() < 1e-12 {
                if b != 0. { roots[0] = -c / b; }
            } else {
                let disc = b * b - 4. * a * c;
                if disc >= 0. {
                    let s = disc.sqrt();
                    roots = [(-b + s) / (2. * a), (-b - s) / (2. * a)];
                }
            }
            for t in roots.into_iter().filter(|t| *t > 0. && *t < 1.) {
                let p = self.eval(t);
                lo = P { x: lo.x.min(p.x), y: lo.y.min(p.y) };
                hi = P { x: hi.x.max(p.x), y: hi.y.max(p.y) };
            }
        }
        (lo, hi)
    }

    /// Wang's formula: segment count n so a uniform-in-t polyline stays within `tol`.
    /// |B''| <= 6*M with M = max second difference of control points; error <= |B''|/(8 n^2).
    pub fn wang_segments(&self, tol: f64) -> usize {
        let [p0, p1, p2, p3] = self.0;
        let dd = |a: P, b: P, c: P| (a.x - 2. * b.x + c.x).hypot(a.y - 2. * b.y + c.y);
        let m = dd(p0, p1, p2).max(dd(p1, p2, p3));
        ((0.75 * m / tol).sqrt().ceil() as usize).max(1)
    }

    /// Append the polyline (excluding the start point) with deviation <= tol.
    pub fn flatten(&self, tol: f64, out: &mut Vec<P>) {
        let n = self.wang_segments(tol);
        for i in 1..=n { out.push(self.eval(i as f64 / n as f64)); }
    }
}

Checked with tests: split halves match the original curve to 1e-9; the flattened polyline stays within tol for tol 1, 0.25, 0.05 (measured worst deviation 0.52, 0.14, 0.03 on a 400-unit S-curve, so the bound is about 2x conservative); bbox matches a 100,000-sample brute force.

Complexity and constants

• Evaluation: 3 + 2 + 1 = 6 lerps per cubic, each 2 mul-adds per axis. Power-basis Horner form is cheaper (3 mul-adds per axis) but gives no split.
• Flatten count grows as sqrt(M/tol): a 100 px quarter circle needs 12 segments at tol 0.25 px (measured, formula above), 19 at 0.1 px. Halving tol multiplies n by about 1.41.
• Data: 4 points x 2 x f32 = 32 bytes per cubic; the flattened output dominates memory, not the input.
• Evaluation is compute-bound and branch-free; the 4 points stay in registers.

Efficiency

• Cost: a cubic evaluation is 6 lerps (2 mul-adds each per axis), compute-bound and branch-free with the 4 points in registers. A 100 px quarter circle needs 12 segments at tolerance 0.25 px and 19 at 0.1 px (measured with the formula above). Memory: 32 B per cubic as 4 x 2 x f32; the flattened output dominates. Bound: flattening compute, then fill rate and memory bandwidth in the rasterizer.
• Levers:
  • Wang's formula for n — one count from four points, no recursion or branching; Skia states it as the minimum evenly spaced count within tolerance — GPU tessellation
  • tolerance in device pixels (order of 0.1-0.25 px at 1x) — zoomed-out paths are not over-tessellated — any renderer
  • cull on the control-point hull before flattening — 4 points, no root solve — clipped or dirty-rectangle drawing
  • forward differencing — 3 additions per axis per point, but accumulates float error; use f64 or n up to a few dozen — fixed n
• Naive → standard → state of the art:
  • fixed 100 segments per curve · 100 points however flat the curve
  • Wang's formula per curve · sqrt(M/tol) segments, 12 for the example above
  • adaptive subdivision on uneven curves (Levien's quadratic flattening) · fewer segments, but branching

Implementation notes

• Use f32 on the GPU; keep f64 on the CPU for geometry that is subdivided repeatedly or has large coordinates.
• Store paths as flat [f32] of points plus a verb array; do not allocate per curve. Reuse the out vector across curves.
• Flatten in device space after applying the transform, or scale tol by the transform's largest singular value.
• Arc length has no closed form: sum the flattened segments, or use 5-point Gauss-Legendre per piece. To get equal spacing along the curve, build a table of (t, cumulative length) and binary-search it.
• Degenerate control points (P1 = P0) make the tangent zero at t=0; stroking and normals need a fallback direction.

Pitfalls

• Treating t as arc length: uniform t clusters points where control points are close. Animation easing is exactly this effect.
• Using the control polygon as the bounding box: it is too large (the curve touches the hull only at the ends); the derivative roots give the tight box.
• Using the single endpoint check for extrema: curves overshoot endpoints (see the test's S-curve).
• Approximating a circle: the standard quarter-arc cubic uses handle length k = 4/3·tan(π/8) ≈ 0.5522847498 times the radius; the radial error is about 2.7e-4 of the radius (measured), fine for display but not a circle.
• Splitting a cubic at its inflection points is not needed for flattening but is needed for offsetting and for cubic-to-quad conversion.

Variants and alternatives

In production

• Typography: OpenType CFF2 charstrings describe each contour with line segments and cubic Bézier curves (spec); TrueType glyf outlines mark points on and off the curve (spec); FreeType and every font rasterizer consume them (font rasterization, OpenType parsing).
• Web and UI: SVG path data has cubic and quadratic Bézier commands (spec); CSS cubic-bezier() is a cubic with fixed endpoints (0,0) and (1,1) (spec, animation easing).
• GPU 2-D rendering: Skia's GPU path tessellation uses Wang's formula (Sources).
• Print and vector formats: PostScript and PDF path operators are Bézier.

Implementations

Sources

• Mike Kamermans (Pomax), A Primer on Bézier Curves, web. online
• Ron Goldman, Pyramid Algorithms: A Dynamic Programming Approach to Curves and Surfaces for Geometric Modeling, Morgan Kaufmann 2003 (section 5.6.3, Wang's formula).
• Thomas W. Sederberg, Computer Aided Geometric Design (course notes). BYU ScholarsArchive
• Raph Levien, Flattening quadratic Béziers, 2019. blog
• Skia, src/gpu/tessellate/WangsFormula.h. code