Introduction to Bézier Curves

Bézier curves turn a short list of points into a smooth, editable path. They are the curves behind SVG paths, OpenType outlines, illustration software, animation timing functions, CAD sketches, and many tool paths. Their practical strength is not merely that they look smooth: the same representation supports stable evaluation, subdivision, derivatives, bounds, flattening, hit testing, and intersection algorithms.

Start with the curve below. The black points are endpoints and the outlined points are controls. Drag any point. The curve always begins at \(\mathbf P_0\), ends at \(\mathbf P_3\), and remains inside the convex hull of its control points. Everything developed below follows from one operation: linear interpolation.

Points, Parameters, and Linear Interpolation

A parametric curve gives every coordinate as a function of one shared parameter. For a planar curve,

\[ \mathbf B(t)=\begin{pmatrix}x(t)\\y(t)\end{pmatrix},\qquad 0\leq t\leq1. \]

The parameter orders points along the curve, but it is not generally distance. Halfway in parameter space need not be halfway along the path. The elementary building block is the affine combination

\[ \operatorname{lerp}(\mathbf P,\mathbf Q;t) =(1-t)\mathbf P+t\mathbf Q =\mathbf P+t(\mathbf Q-\mathbf P). \]

At \(t=0\) it returns \(\mathbf P\); at \(t=1\) it returns \(\mathbf Q\). For \(0<t<1\), both coefficients are nonnegative and sum to one, so the result lies on the segment between the two points. That simple coefficient rule will later explain the convex-hull property of an entire Bézier curve.

de Casteljau's Construction

Paul de Casteljau developed this recursive construction at Citroën in 1959. Pierre Bézier's independent work at Renault later made the polynomial curve representation widely known in geometric design. Every evaluation in either formulation is built from repeated linear interpolation.

Given control points \(\mathbf P_0,\ldots,\mathbf P_n\), interpolate every adjacent pair using the same \(t\). Repeat on the shorter list until one point remains. Define

\[ \mathbf P_i^{(0)}=\mathbf P_i,\qquad \mathbf P_i^{(r)}(t)=(1-t)\mathbf P_i^{(r-1)}(t)+t\mathbf P_{i+1}^{(r-1)}(t). \]

After \(n\) rounds, \(\mathbf B(t)=\mathbf P_0^{(n)}(t)\). Change the degree, move the slider, or drag the control points below. Colored segments show the successive interpolation levels; the larger red point is the final result.

This recursive construction is not merely a drawing aid. It is numerically stable, works in any dimension, and produces all control points needed to split the curve. The implementation is short:

function evaluate(points, t) {
  let level = points.map(point => ({ ...point }));
  while (level.length > 1) {
    level = level.slice(0, -1).map((point, i) => ({
      x: point.x + (level[i + 1].x - point.x) * t,
      y: point.y + (level[i + 1].y - point.y) * t
    }));
  }
  return level[0];
}

The Bernstein Form

Expanding the repeated interpolations yields a polynomial basis. A degree-\(n\) Bézier curve is

\[ \mathbf B(t)=\sum_{i=0}^{n} b_{i,n}(t)\mathbf P_i, \qquad b_{i,n}(t)=\binom ni(1-t)^{n-i}t^i. \]

The functions \(b_{i,n}\) are the Bernstein basis polynomials. For every \(t\in[0,1]\), they are nonnegative and form a partition of unity:

\[ \sum_{i=0}^{n}b_{i,n}(t)=((1-t)+t)^n=1. \]

Therefore \(\mathbf B(t)\) is always a convex combination of the control points. The graph shows the four cubic basis weights. At any vertical slice they add to one; the curve point is the control-point average using exactly those weights.

Several additional identities explain why this basis is so useful. Each basis function has an order- \(i\) zero at \(t=0\), an order-\(n-i\) zero at \(t=1\), and the symmetry \(b_{i,n}(t)=b_{n-i,n}(1-t)\). Its maximum on \([0,1]\) occurs at \(t=i/n\). The first two normalized factorial moments are

\[ t=\sum_{i=0}^{n}\frac{i}{n}b_{i,n}(t),\qquad t^2=\sum_{i=0}^{n}\frac{i(i-1)}{n(n-1)}b_{i,n}(t)\quad(n\ge2). \]

Equivalently, if \(I\) is binomially distributed with parameters \(n\) and \(t\), then the Bernstein weights are its probabilities. The first identity says \(\operatorname E[I/n]=t\); the variance \(t(1-t)/n\) measures how tightly those weights concentrate around \(t\). The basis itself can be generated without powers or binomial coefficients by the Pascal-like recursion

\[ b_{i,m}(t)=(1-t)b_{i,m-1}(t)+t b_{i-1,m-1}(t), \]

with \(b_{0,0}=1\) and out-of-range indices interpreted as zero. This is the scalar counterpart of de Casteljau's construction.

Bernstein Approximation and the Weierstrass Theorem

For a continuous scalar function \(f:[0,1]\to\mathbb R\), sample \(f\) at the uniform nodes and use those samples as coefficients:

\[ \mathcal B_n[f](t)=\sum_{i=0}^{n}f\!\left(\frac{i}{n}\right)b_{i,n}(t). \]

The probabilistic reading gives \(\mathcal B_n[f](t)=\operatorname E[f(I/n)]\). As \(n\) grows, the variance of \(I/n\) tends uniformly to zero. Uniform continuity of \(f\) then makes the contribution near \(t\) arbitrarily accurate, while the probability of the remaining samples tends uniformly to zero. Hence \(\mathcal B_n[f]\to f\) uniformly. This constructive form of the Weierstrass approximation theorem uses the same control-point machinery as a Bézier curve: the graph is the curve with controls \((i/n,f(i/n))\), because its x-coordinate is exactly \(t\) by the first-moment identity.

Linear, Quadratic, and Cubic Curves

The first three degrees are worth knowing explicitly:

\[ \begin{aligned} \mathbf L(t)&=(1-t)\mathbf P_0+t\mathbf P_1,\\ \mathbf Q(t)&=(1-t)^2\mathbf P_0+2(1-t)t\mathbf P_1+t^2\mathbf P_2,\\ \mathbf C(t)&=(1-t)^3\mathbf P_0+3(1-t)^2t\mathbf P_1+3(1-t)t^2\mathbf P_2+t^3\mathbf P_3. \end{aligned} \]

Quadratics have one internal control point and cannot inflect. Cubics have two internal controls and can form an S-shape, cusp, or loop. Graphics systems usually join many low-degree segments instead of using one high-degree curve: local edits remain local, evaluation stays cheap, and numerical behavior is easier to control.

Geometric Properties

The Bernstein representation gives several useful facts immediately:

• Endpoint interpolation: \(\mathbf B(0)=\mathbf P_0\) and \(\mathbf B(1)=\mathbf P_n\).
• Convex hull: the curve lies inside the convex hull of its controls.
• Affine invariance: translating, rotating, scaling, or shearing the controls transforms the curve in exactly the same way.
• Variation diminution: a line intersects the curve no more often than it intersects the control polygon, counting multiplicity appropriately.
• Endpoint tangents: the first and last control-polygon edges determine the endpoint directions.

Affine invariance is especially important in software. Evaluation may happen in object coordinates and the result can then be transformed, or the controls may be transformed first; both routes agree. Perspective projection is not affine, which is one reason rational Bézier curves and homogeneous coordinates matter in 3D pipelines.

To see why, write an affine map as \(A(\mathbf x)=M\mathbf x+\mathbf c\). Because the Bernstein weights sum to one,

\[ \begin{aligned} A(\mathbf B(t)) &=M\sum_{i=0}^{n}b_{i,n}(t)\mathbf P_i+\mathbf c\\ &=\sum_{i=0}^{n}b_{i,n}(t)\bigl(M\mathbf P_i+\mathbf c\bigr) =\sum_{i=0}^{n}b_{i,n}(t)A(\mathbf P_i). \end{aligned} \]

Thus an affine transformation can be applied to the few control points rather than to every evaluated point on the curve. Translation, rotation, reflection, shear, and nonuniform scaling all satisfy this identity.

Rational Bézier Curves

Polynomial Bézier curves are closed under affine transformations, but perspective projection is not affine. Introduce one scalar weight \(w_i\) for every control point and form the rational curve

\[ \mathbf R(t)= \frac{\sum_{i=0}^{n}w_i b_{i,n}(t)\mathbf P_i} {\sum_{i=0}^{n}w_i b_{i,n}(t)}. \]

This is not a separate interpolation trick. Lift each \(d\)-dimensional control point into homogeneous \((d+1)\)-space,

\[ \widehat{\mathbf P}_i=(w_i\mathbf P_i,w_i). \]

Ordinary de Casteljau interpolation of the lifted controls gives \(\widehat{\mathbf R}(t)=(\mathbf X(t),W(t))\). Projecting back by \(\mathbf R(t)=\mathbf X(t)/W(t)\) produces the quotient above. Consequently, the same stable triangular construction works for every degree and any number of spatial coordinates; division happens only once, after the final interpolation. A zero \(W(t)\) is a projective point at infinity and cannot be returned as a finite Euclidean point.

Multiplying every weight by the same nonzero constant leaves the curve unchanged. Equal weights recover the polynomial Bézier curve. Positive weights preserve endpoint interpolation and the convex-hull property because the normalized coefficients \(w_i b_{i,n}(t)/\sum_j w_j b_{j,n}(t)\) remain nonnegative and sum to one. Increasing one positive weight pulls the curve toward its control point. Zero or negative weights are meaningful in projective constructions, but they can remove the convex-hull guarantee or make the denominator vanish.

An Exact Quarter Circle

A nonconstant polynomial parameterization cannot satisfy \(x(t)^2+y(t)^2=1\) identically, but a rational quadratic can. Choose

\[ \mathbf P_0=(1,0),\qquad \mathbf P_1=(1,1),\qquad \mathbf P_2=(0,1), \qquad (w_0,w_1,w_2)=\left(1,\frac{\sqrt 2}{2},1\right). \]

The endpoint derivatives follow the coordinate axes, and symmetry places the midpoint on \(x=y\). Requiring that midpoint to equal \((1/\sqrt 2,1/\sqrt 2)\) determines \(w_1=\sqrt 2/2\). With \(u=1-t\), its homogeneous coordinates are

\[ X=u^2+\sqrt 2ut,\qquad Y=\sqrt 2ut+t^2,\qquad W=u^2+\sqrt 2ut+t^2. \]

Direct expansion gives \(X^2+Y^2=W^2\), hence the projected coordinates satisfy \(x(t)^2+y(t)^2=1\) for every \(t\), not merely at sampled parameters. More generally, for a circular sweep \(|\theta|<\pi\), use the arc endpoints as the outer controls, their tangent-line intersection as the middle control, and \(w_1=\cos(\theta/2)\) with endpoint weights one. Split semicircles and longer sweeps into smaller segments to avoid a zero or negative middle weight.

Subdivision

The edges of the de Casteljau triangle already contain two new control polygons. At a split parameter \(s\), the left edge controls the interval \([0,s]\), and the right edge in reverse order controls \([s,1]\). Together the two curves reproduce the original exactly.

For a cubic \(\mathbf P_0,\mathbf P_1,\mathbf P_2,\mathbf P_3\), write one complete de Casteljau triangle as

\[ \begin{aligned} \mathbf A&=(1-s)\mathbf P_0+s\mathbf P_1,& \mathbf B&=(1-s)\mathbf P_1+s\mathbf P_2,& \mathbf C&=(1-s)\mathbf P_2+s\mathbf P_3,\\ \mathbf D&=(1-s)\mathbf A+s\mathbf B,& \mathbf E&=(1-s)\mathbf B+s\mathbf C,& \mathbf F&=(1-s)\mathbf D+s\mathbf E. \end{aligned} \]

The split point is \(\mathbf F=\mathbf B(s)\). The left child has controls \(\mathbf P_0,\mathbf A,\mathbf D,\mathbf F\), and the right child has controls \(\mathbf F,\mathbf E,\mathbf C,\mathbf P_3\). At \(s=1/2\), every operation is an average; no separate midpoint formula is needed.

The controls are read directly from the two boundary chains of the triangle, not fitted to sampled points. Moreover, \(\mathbf D,\mathbf F,\mathbf E\) are collinear because \(\mathbf F=(1-s)\mathbf D+s\mathbf E\). The last edge of the left polygon and the first edge of the right polygon therefore share the original tangent direction at the split. With local child parameters \(u\), their endpoint derivatives are \(s\mathbf B'(s)\) and \((1-s)\mathbf B'(s)\), respectively.

Subdivision is the workhorse behind adaptive rendering and intersections. It avoids deriving a new equation for each operation: repeatedly split, reject pieces that cannot matter, and continue only where more resolution is needed. A subcurve over an arbitrary interval \([a,b]\) follows from splitting at \(b\), retaining the left curve, then splitting that curve at \(a/b\) and retaining the right part.

Repeated subdivision does not change the curve. It replaces its control polygon by smaller child polygons whose maximum distance from the curve tends to zero. Connecting the accepted child endpoints therefore yields a convergent polyline approximation; the flatness test below turns this fact into an adaptive stopping rule.

Degree Elevation and Quadratic-to-Cubic Conversion

Any degree-\(n\) curve can be represented exactly at degree \(n+1\). The elevated controls are

\[ \mathbf P'_0=\mathbf P_0,\qquad \mathbf P'_i=\frac{i}{n+1}\mathbf P_{i-1}+\left(1-\frac{i}{n+1}\right)\mathbf P_i, \qquad \mathbf P'_{n+1}=\mathbf P_n. \]

For an SVG quadratic command with start \(\mathbf Q_0\), control \(\mathbf Q_1\), and end \(\mathbf Q_2\), the equivalent cubic controls are

\[ \mathbf C_0=\mathbf Q_0,\quad \mathbf C_1=\mathbf Q_0+\frac23(\mathbf Q_1-\mathbf Q_0),\quad \mathbf C_2=\mathbf Q_2+\frac23(\mathbf Q_1-\mathbf Q_2),\quad \mathbf C_3=\mathbf Q_2. \]

function quadraticToCubic(q0, q1, q2) {
  return [
    q0,
    add(q0, scale(sub(q1, q0), 2 / 3)),
    add(q2, scale(sub(q1, q2), 2 / 3)),
    q2
  ];
}

Degree reduction is different: a generic cubic cannot be represented exactly as a quadratic. Reduction is an approximation problem, so its result is meaningless until an error measure and any endpoint or tangent constraints are specified. For a degree-\(n\) curve with controls \(\mathbf P_i\), one useful degree-\(m\) approximation minimizes the parameter-space squared error

\[ E=\int_0^1\left\| \sum_{i=0}^{n}b_{i,n}(t)\mathbf P_i- \sum_{j=0}^{m}b_{j,m}(t)\mathbf Q_j \right\|^2\,dt. \]

Differentiating with respect to each unknown \(\mathbf Q_k\) gives the linear system

\[ \sum_{j=0}^{m}G_{kj}\mathbf Q_j= \sum_{i=0}^{n}H_{ki}\mathbf P_i, \qquad G_{kj}=I(m,k,m,j),\quad H_{ki}=I(m,k,n,i), \]

where products of Bernstein basis functions integrate exactly:

\[ I(a,r,b,s)=\int_0^1 b_{r,a}(t)b_{s,b}(t)\,dt =\frac{\binom{a}{r}\binom{b}{s}} {(a+b+1)\binom{a+b}{r+s}}. \]

Endpoint interpolation is imposed by fixing \(\mathbf Q_0=\mathbf P_0\) and \(\mathbf Q_m=\mathbf P_n\), moving their known contributions to the right-hand side, and solving only for the interior controls. A stable linear solver should operate on these small matrices directly; explicit factorial products overflow quickly, while a recurrence specialized to one matrix size is harder to verify. This least-squares projection minimizes error at equal parameter values. It does not minimize Hausdorff distance or compensate for a different parameterization, so geometric reduction may additionally reparameterize, split, and test the result against a spatial tolerance.

const quadratic = cubic.reduce(2);
const unconstrained = cubic.reduce(2, { preserveEndpoints: false });

reduce uses the integral above and preserves both endpoints by default. Setting preserveEndpoints to false computes the unconstrained projection, which can lower the integrated error while moving the ends. Neither mode silently claims a geometric minimax result.

Font Outlines

Font outlines are a concrete reason to distinguish quadratic and cubic segments. Traditional TrueType glyf outlines use straight segments and quadratic Bézier curves. A stored off-curve point is a quadratic control; two consecutive off-curve points imply an on-curve point at their midpoint, which keeps many contours compact. OpenType is a container rather than a single curve format: outlines in glyf remain quadratic, while CFF and CFF2 outlines use cubic Bézier curves.

Rendering and editing need no special curve mathematics beyond the operations already developed. Apply the font and layout transformations to the controls by affine invariance, split contours with de Casteljau when local processing is needed, and convert a quadratic segment to a cubic exactly with the controls above when a cubic-only path system is the destination. Conversion in the opposite direction is generally approximate.

Derivatives, Tangents, and Normals

Differentiating the Bernstein form produces another Bézier curve, one degree lower:

\[ \mathbf B'(t)=n\sum_{i=0}^{n-1}b_{i,n-1}(t)(\mathbf P_{i+1}-\mathbf P_i). \]

Its control points are \(n(\mathbf P_{i+1}-\mathbf P_i)\). In particular, \(\mathbf B'(0)=n(\mathbf P_1-\mathbf P_0)\) and \(\mathbf B'(1)=n(\mathbf P_n-\mathbf P_{n-1})\). When \(\mathbf B'(t)\neq\mathbf0\), the planar unit tangent and one consistently oriented unit normal are

\[ \mathbf T(t)=\frac{\mathbf B'(t)}{\|\mathbf B'(t)\|}, \qquad \mathbf N(t)=\begin{pmatrix}-T_y(t)\\T_x(t)\end{pmatrix}. \]

Curvature and Inflections

For a regular planar parametric curve, signed curvature is

\[ \kappa(t)= \frac{x'(t)y''(t)-y'(t)x''(t)}{\left(x'(t)^2+y'(t)^2\right)^{3/2}}. \]

Its magnitude measures turning per unit arc length; its sign records the turning side. Where \(\kappa\neq0\), the osculating circle has radius \(R=1/|\kappa|\) and center \(\mathbf B+(1/\kappa)\mathbf N\) for the chosen normal convention. An ordinary inflection occurs where \(x'y''-y'x''=0\) and the sign changes. Quadratics do not have interior inflections; cubics can have up to two.

A zero derivative is a singular point: tangent normalization and the curvature formula both fail there. Robust code must detect the small denominator rather than allowing infinities to contaminate later geometry.

Normals in Three Dimensions

In the plane, rotating a tangent by 90 degrees chooses a normal. In 3D, every vector in a whole plane is perpendicular to the tangent, so extra structure is required. For a regular spatial curve, curvature is

\[ \kappa(t)=\frac{\|\mathbf B'(t)\times\mathbf B''(t)\|}{\|\mathbf B'(t)\|^3}. \]

Where the cross product is nonzero, exact polynomial derivatives give the Frenet frame without sampling a nearby parameter:

\[ \mathbf T=\frac{\mathbf B'}{\|\mathbf B'\|},\qquad \mathbf B_f=\frac{\mathbf B'\times\mathbf B''}{\|\mathbf B'\times\mathbf B''\|},\qquad \mathbf N=\mathbf B_f\times\mathbf T. \]

The frame becomes undefined where curvature vanishes and may flip near inflections. For camera rails, swept surfaces, and ribbons, a rotation-minimizing frame transported along sampled points is usually more stable and visually quieter. The library therefore reports null instead of inventing a direction at a straight point.

Tensor-Product Bézier Surfaces

A rectangular control net \(\mathbf P_{ij}\), with \(0\le i\le n\) and \(0\le j\le m\), extends the curve construction to two parameters:

\[ \mathbf S(u,v)=\sum_{i=0}^{n}\sum_{j=0}^{m} b_{i,n}(u)b_{j,m}(v)\mathbf P_{ij},\qquad (u,v)\in[0,1]^2. \]

The product weights are nonnegative and sum to one, so the entire patch remains in the convex hull of its control net. Evaluation requires no new algorithm: first run de Casteljau along every row at \(v\), then run it once more through the resulting points at \(u\). Reversing the two directions gives the same point. Fixing one parameter produces an ordinary Bézier iso-parameter curve.

Differentiating in either parameter turns neighboring control points into a derivative control net:

\[ \begin{aligned} \mathbf S_u(u,v)&=n\sum_{i=0}^{n-1}\sum_{j=0}^{m} (\mathbf P_{i+1,j}-\mathbf P_{ij})b_{i,n-1}(u)b_{j,m}(v),\\ \mathbf S_v(u,v)&=m\sum_{i=0}^{n}\sum_{j=0}^{m-1} (\mathbf P_{i,j+1}-\mathbf P_{ij})b_{i,n}(u)b_{j,m-1}(v). \end{aligned} \]

Where \(\mathbf S_u\times\mathbf S_v\) is nonzero, the unit normal is \(\mathbf N=(\mathbf S_u\times\mathbf S_v)/\|\mathbf S_u\times\mathbf S_v\|\). In the projected patch below, blue and green curves hold \(v\) and \(u\) constant respectively. The short blue and green segments are the two tangent directions; red is the surface normal.

Why B-Splines Behave More Locally

A single degree-\(n\) Bézier segment uses every control point throughout \([0,1]\), so moving one control generally changes the whole segment. B-spline basis functions instead have compact support on a knot vector; moving one de Boor control point changes only the parameter intervals covered by the corresponding basis function. They retain nonnegativity, partition of unity, convex-hull bounds, and recursive evaluation, but use the de Boor algorithm rather than de Casteljau. Piecewise Bézier segments remain preferable when explicit segment boundaries and simple interchange are the main concerns; B-splines are preferable when local editing across a longer smooth curve or surface matters.

Joining Segments Smoothly

Suppose one segment ends where the next begins. Several notions of continuity must be distinguished:

• \(C^0\): positions agree. There is no gap.
• \(G^1\): tangent directions agree up to a positive scale. The join looks smooth but speed may jump.
• \(C^1\): first derivatives agree exactly. Direction and parameter speed agree.
• \(G^2\): signed curvature agrees after allowing reparameterization.
• \(C^2\): first and second derivatives agree exactly.

For adjacent cubics \(\mathbf P_0,\ldots,\mathbf P_3\) and \(\mathbf Q_0,\ldots,\mathbf Q_3\), \(C^0\) requires \(\mathbf P_3=\mathbf Q_0\), while \(C^1\) additionally requires

\[ \mathbf P_3-\mathbf P_2=\mathbf Q_1-\mathbf Q_0. \]

Collinearity with unequal handle lengths gives \(G^1\), not \(C^1\). Editors often expose this useful geometric condition because it preserves a smooth appearance while allowing independent timing on both segments.

Adaptive Flattening and Flatness

Displays and many fabrication formats ultimately consume line segments. Uniformly stepping through \(t\) is easy but wasteful on straight regions and insufficient in sharply curved regions. Adaptive flattening instead asks whether a segment is flat enough; if not, it splits and tests both halves.

A useful flatness measure is the largest perpendicular distance from any internal control point to the baseline through the endpoints. For baseline endpoints \(\mathbf A\), \(\mathbf B\) and control \(\mathbf P_i\),

\[ f=\max_{1\leq i<n} \frac{|(\mathbf B-\mathbf A)\times(\mathbf P_i-\mathbf A)|}{\|\mathbf B-\mathbf A\|}. \]

function flatten(curve, tolerance, output) {
  if (curve.flatness() <= tolerance) {
    output.push(curve.end);
    return;
  }
  const { left, right } = curve.split(0.5);
  flatten(left, tolerance, output);
  flatten(right, tolerance, output);
}

The tolerance belongs in output space. If one device pixel is the smallest meaningful visual error, a tolerance near half a pixel is a sensible starting point. Zooming or applying a nonuniform transform changes the relevant scale, so flattening solely in object coordinates can violate a screen-space error guarantee. Flatness is a conservative practical criterion, not an equality for the maximum curve-to-chord deviation.

Extrema and Bounding Boxes

The control-point box is cheap and always contains the curve, but it is generally loose. A tight axis-aligned box must include endpoints and every interior parameter where an axis derivative vanishes:

\[ x'(t)=0\quad\text{or}\quad y'(t)=0,\qquad 0<t<1. \]

A cubic's derivatives are quadratic, so their roots follow from the quadratic formula. Evaluate the original curve at all accepted roots and take coordinate-wise minima and maxima. The pale rectangle below bounds the control polygon; the red rectangle is the exact axis-aligned curve box.

Bounding boxes are rejection tests, not proof of intersection. Nonoverlapping boxes prove that two pieces cannot intersect; overlapping boxes say only that further work is required. For a tighter oriented box, translate the first endpoint to the origin, rotate the endpoint chord onto the x-axis, compute a box there, and transform its corners back.

Arc Length and Constant-Speed Travel

The arc length from \(0\) to \(u\) is the speed integral

\[ s(u)=\int_0^u\|\mathbf B'(t)\|\,dt =\int_0^u\sqrt{x'(t)^2+y'(t)^2}\,dt. \]

Generic cubic Bézier arc length has no elementary closed form. Numerical quadrature is therefore normal, not a fallback of last resort. Adaptive Simpson integration and Gauss-Legendre quadrature both work well when paired with explicit error control. A flattened polyline gives a fast lower estimate whose accuracy follows its flattening tolerance.

Red markers use equal increments of \(t\); black markers use equal arc-length increments. To move at constant speed, solve \(s(t)=d\) for the desired distance \(d\). Since \(s\) is monotone on a regular curve, bisection is reliable. Repeated animation queries are cheaper with a lookup table of \((t,s)\) pairs followed by local interpolation or Newton refinement.

Separating the Path from Its Timing

A path answers where an object may travel; a timing law answers when it reaches each place. Keeping those questions separate avoids changing the geometry merely to change the speed. Let \(L=s(1)\) and define normalized arc length

\[ \sigma(u)=\frac{s(u)}{L},\qquad U=\sigma^{-1}. \]

A monotone time-distance law \(S:[0,1]\to[0,1]\) then produces the space-time curve

\[ \boxed{\mathbf p(\tau)=\mathbf B\!\left(U(S(\tau))\right)}. \]

The three maps have distinct jobs: \(\mathbf B(u)\) supplies the spatial path, \(U\) converts normalized distance back to the curve parameter, and \(S(\tau)\) controls motion over normalized time. With physical duration \(T\), the speed is

\[ v(\tau)=\frac{L}{T}S'(\tau). \]

Thus \(S(\tau)=\tau\) gives constant speed even when equal parameter increments do not. The cubic smoothstep \(S(\tau)=3\tau^2-2\tau^3\) starts and ends at rest, with \(S'(\tau)=6\tau(1-\tau)\). Its acceleration is not zero at the endpoints; when endpoint acceleration must also vanish, the quintic smootherstep \(6\tau^5-15\tau^4+10\tau^3\) is the usual replacement.

The small markers represent equal time increments. Constant speed spaces the black markers evenly by distance; the red smoothstep markers cluster near both endpoints and spread out near the middle. Acceleration has both a tangential part, controlled by \(S''\), and a normal part caused by curvature: \(\mathbf a=\dot v\,\mathbf T+\kappa v^2\mathbf N\). Consequently, a timing law can smooth changes in speed, but it cannot remove the centripetal acceleration required by a sharply curved path.

CSS timing functions use the same normalized time-distance idea. Their horizontal coordinate is time and their vertical coordinate is progress; the CSS Cubic Bézier Editor evaluates that inversion explicitly. For motion through points assigned to absolute times, apply one monotone timing law per path interval and match endpoint velocities when a speed-continuous join is required.

Projecting a Point onto the Curve

Hit testing, dragging a point on a path, and measuring distance all require the nearest curve point to a query \(\mathbf X\). It minimizes

\[ D(t)=\|\mathbf B(t)-\mathbf X\|^2, \]

whose stationary points satisfy

\[ D'(t)=2(\mathbf B(t)-\mathbf X)\cdot\mathbf B'(t)=0. \]

For a cubic, this is generally a degree-five equation. Isolate every real root in \([0,1]\) by recursively partitioning at the roots of its derivative, refine the sign-changing intervals with bisection, and evaluate \(D\) at every resulting stationary point. Repeated roots occur at derivative boundaries and must be tested explicitly. The endpoints \(t=0\) and \(t=1\) remain candidates because a constrained minimum need not be stationary. Comparing this finite candidate set yields the global nearest point rather than merely the nearest point in a sampled neighborhood. Drag the black query point below.

Molding a Curve

Once projection supplies a parameter \(t\), an editor can let the user grab the curve itself rather than a handle. For a quadratic, solving the Bernstein equation for its only control point is direct:

\[ \mathbf P_1= \frac{\mathbf X-(1-t)^2\mathbf P_0-t^2\mathbf P_2}{2(1-t)t}. \]

A cubic has two unknown controls but only one point constraint, so infinitely many solutions exist. A molding tool must add a policy: preserve the local tangent, minimize control-point displacement, maintain symmetry, or solve a constrained least-squares problem. There is no single geometrically mandatory cubic answer.

const molded = Bezier.quadraticFromPoints(p0, draggedPoint, p2, t);

Intersections

Curve and Line

Apply an affine line predicate to the curve, or translate and rotate both objects until the line lies on the x-axis. Intersections then occur where the resulting scalar polynomial is zero. A degree-\(n\) curve produces a polynomial of degree at most \(n\): quadratic curves generally require quadratic root finding and cubic curves require cubic root finding. Retain roots in \([0,1]\), then verify that the corresponding line parameter lies inside the line segment rather than merely on its infinite extension. The complete line-segment intersection follows directly from this predicate.

Keeping the scalar polynomial in Bernstein form adds a useful rejection test: its value lies between its smallest and largest coefficients. If that interval excludes zero, no root exists. On retained intervals, de Casteljau subdivision produces exact child coefficients, while Bézier clipping intersects the coefficient control hull with the zero axis to contract the parameter range more aggressively. These interval methods extend naturally to higher degrees and preserve tangential candidates that endpoint sign-change tests can miss.

Another global construction recursively isolates the roots of the derivative, whose Bernstein coefficients are simply \(n(c_{i+1}-c_i)\) for scalar controls \(c_i\). Those critical parameters partition \([0,1]\) into monotone intervals. Test every critical point itself to retain even-multiplicity roots, and use bisection only on intervals whose endpoint values have opposite signs. Linear, quadratic, and cubic recursion bases can be solved analytically. This gives an explicit parameter tolerance without relying on a fixed sampling density.

Curve and Curve

A robust general strategy uses subdivision and bounding boxes:

1. Reject a pair whose bounding boxes do not overlap.
2. Split the larger or less flat curve.
3. Test the resulting pairs recursively.
4. When both boxes are below a geometric tolerance, record a parameter pair.
5. Cluster neighboring candidates and refine them with a two-variable solver if higher precision is required.

Tangential contacts are harder than crossings because no sign change occurs. Pure polyline intersection can miss them unless flattening error is included in the acceptance test. Overlapping or coincident curve intervals also require a result model richer than a list of isolated points.

What Subdivision Actually Proves

Disjoint convex hulls are a certificate of non-intersection because each Bézier segment lies in its hull. Overlapping hulls prove only that the pair remains a candidate. Stopping when both boxes are smaller than an arbitrary \(\varepsilon\) can report a near miss as a hit, merge distinct nearby intersections, or return several candidates for one contact. The library method intersects() deliberately has this practical contract: it returns clustered approximate candidates at the requested geometric tolerance, not a topological certificate.

A complete method must distinguish transversal crossings, tangential crossings of odd multiplicity, and tangential non-crossing contacts of even multiplicity. Yap's exact subdivision framework replaces a display-scale cutoff by input-dependent separation bounds. Once a candidate pair is smaller than such a bound, it can contain at most one isolated intersection or relevant critical feature. Elementary convex or concave graph pieces can then be coupled through their normal fields; the signs of the endpoint normal-angle differences decide whether a non-crossing tangency exists. Stationary points, extrema, and inflections must first be isolated so that this elementary-curve criterion applies.

That completeness result assumes exact bigfloat input and requires bounds derived from the degrees and coefficient heights of the original curves. It is substantially more than recursive box overlap, and those guarantees cannot be retrofitted by choosing a smaller floating-point tolerance. Curves sharing a component are a separate output case: their intersection is an interval rather than an isolated parameter pair.

Circle Intersections

For center \(\mathbf C\) and radius \(r\), solve \(\|\mathbf B(t)-\mathbf C\|^2-r^2=0\). A cubic produces a degree-six polynomial, so numerical isolation is the practical route. Partition at derivative roots, retain critical points that satisfy the equation for tangencies, bisect sign-changing monotone intervals, and deduplicate roots within a parameter tolerance.

Offset Curves

An exact constant-distance offset would be

\[ \mathbf O(t)=\mathbf B(t)+d\,\mathbf N(t). \]

The normalized normal contains \(1/\|\mathbf B'(t)\|\), which is generally not polynomial. Consequently, a nontrivial polynomial Bézier curve does not generally have an exact polynomial Bézier offset. Software must approximate it, preserve it as a different curve representation, or flatten it.

The visualization samples exact normal-displaced points; connecting them is only a polyline approximation. A vector workflow usually subdivides at extrema and inflections, approximates simple pieces with new cubics, and checks the geometric error. Large inward offsets can form cusps and self-intersections. Building a stroked outline also requires join and cap policies; offsetting the centerline alone is not enough.

Graduated offsets replace the constant \(d\) by \(d(s)\), preferably as a function of arc length rather than \(t\). This keeps width changes visually uniform even when the curve's parameter speed varies.

The polygonal outline API accepts independent left and right endpoint widths. Each pair is interpolated by the normalized arc length \(s(t)/s(1)\), not directly by \(t\):

const tapered = curve.outline({
  left: [0.4, 0.1],
  right: [0.2, 0.05]
}, 64);

A single number still produces a constant symmetric outline. The profile form permits asymmetric pressure or calligraphic effects while retaining explicit semantics. Both sides are sampled at matching parameters and joined by straight end caps; joins between multiple curve segments remain the caller's responsibility. Pass parameterization: 'parameter' only when variation in \(t\), rather than physical distance, is intentional.

Circular Arcs and Bézier Curves

Polynomial Bézier curves cannot represent a nondegenerate circular arc exactly. A polynomial approximation of a circular arc depends on the degree, endpoint contact order, and geometric error being minimized. The standard cubic construction for an arc of radius \(r\), sweep angle \(\theta\), and endpoint tangents places each internal control a distance

\[ k r,\qquad k=\frac43\tan\frac{\theta}{4} \]

from its endpoint along the tangent. For a quarter circle, \(k\approx0.55228475\). Change the angle below; the reported radial deviation makes the approximation error explicit rather than hiding it behind a smooth stroke.

Split sweeps larger than 90 degrees into multiple segments. If exact conics are required, use rational quadratic Bézier curves: homogeneous weights extend the polynomial basis and can represent circles exactly. If the target format natively supports arcs, retain arcs and convert only at an export boundary.

Approximating a Bézier Curve with Circular Arcs

The reverse conversion is useful for CNC controllers and geometry kernels that prefer lines and arcs. An incenter biarc approximation preserves both endpoints and their tangent directions, then recursively subdivides where radial deviation is largest. For a cubic Bézier spiral, all relevant extrema reduce to quartic equations. General cubics must first be split at singularities, inflections, and curvature extrema before the spiral guarantees can be applied.

Fitting Curves to Data

Control points are not samples the curve generally passes through. Constructing a curve from measured points therefore requires both parameter values and a fitting criterion. A common workflow is:

1. Assign parameters by cumulative chord length or centripetal distance.
2. Fix endpoints and estimate endpoint tangents.
3. Solve a linear least-squares problem for the internal cubic controls.
4. Project each sample back onto the curve and update its parameter.
5. Repeat until the error stabilizes; split where the tolerance cannot be met.

Equal parameter spacing is simple but performs poorly when samples are uneven. Chord-length parameterization is a better default. For long freehand paths, fitting many cubic segments under a maximum error is more stable than fitting one high-degree curve.

When the target is a known function rather than scattered data, fitting a Bézier curve to its graph can preserve the exact endpoint values and tangent directions. The remaining handle lengths then form a small nonlinear optimization problem whose solution depends on the selected graph-error metric.

Alignment and Numerical Conditioning

Many calculations simplify after translating \(\mathbf P_0\) to the origin and rotating \(\mathbf P_n-\mathbf P_0\) onto the x-axis. Distances, incidences, parameter values, and curve classification are preserved by this rigid transform. Alignment removes constants from component equations and makes baseline distances direct y-coordinates.

Production implementations should also observe the following:

• Compare floating-point values with scale-aware tolerances, not exact equality.
• Detect coincident endpoints before dividing by baseline length.
• Keep parameter tolerances separate from geometric tolerances.
• Use de Casteljau evaluation near difficult configurations instead of expanded high-degree power polynomials.
• Retain parameter intervals during subdivision; geometry without its parameter mapping is insufficient for later refinement.
• Cluster numerical roots and intersections, but do not merge genuinely distinct nearby features merely because the drawing resolution is low.

A Compact JavaScript API

The visualizations use the standalone library at /js/lib/bezier.min.js. It has no DOM dependency and works in a browser or through CommonJS:

const curve = new Bezier([
  { x: 0, y: 0 },
  { x: 1, y: 3 },
  { x: 4, y: 3 },
  { x: 5, y: 0 }
]);

const point = curve.get(0.4);
const tangent = curve.tangent(0.4);
const { left, right } = curve.split(0.4);
const reduced = curve.reduce(2);
const bounds = curve.bbox();
const polyline = curve.flatten(0.25);
const length = curve.length();
const nearest = curve.project({ x: 2, y: 1 });
const molded = Bezier.quadraticFromPoints(p0, pointOnCurve, p2, 0.4);
const svgPath = curve.toSVG();
const powerCoefficients = curve.toPoly();
const strokePolygon = curve.outline(0.2, 64);
const taperedPolygon = curve.outline({ left: [0.2, 0.05], right: [0.1, 0.05] }, 64);

const approximation = Bezier.bernsteinApproximation(x => Math.abs(x - 0.5), 12);

const surface = new Bezier.Surface(controlGrid);
const surfacePoint = surface.get(0.4, 0.6);
const tangentU = surface.partialU(0.4, 0.6);
const tangentV = surface.partialV(0.4, 0.6);
const surfaceNormal = surface.normal(0.4, 0.6);

const spatialCurve = new Bezier(spatialControls);
const spatialCurvature = spatialCurve.curvature(0.4);
const frame = spatialCurve.frenet(0.4);

const quarterCircle = new Bezier.Rational([
  { x: 1, y: 0 },
  { x: 1, y: 1 },
  { x: 0, y: 1 }
], [1, Math.SQRT1_2, 1]);

const exactPoint = quarterCircle.get(0.5);

The same API evaluates arbitrary-degree rational curves through homogeneous de Casteljau interpolation. For polynomial curves it provides degree elevation, parameter-space least-squares degree reduction, exact quadratic-to-cubic conversion, construction of a quadratic through a specified point at a specified parameter, curvature, inflection parameters, line and curve intersections, arc-length inversion, sampled offsets, cubic circular-arc approximations, complete four-segment circle approximations, and SVG elliptical-arc conversion with radius correction and flag handling. It also constructs classical Bernstein function approximations and evaluates tensor-product surfaces, their iso-parameter curves, partial derivatives, and unit normals. The geometry code remains separate from Plot.js; rendering is one consumer, not a responsibility of the curve model.

Curves may also be constructed from coordinate pairs, as in new Bezier(0, 0, 1, 3, 4, 3, 5, 0). The aliases eval, getLUT, linearize, subdivide, cut, hodograph, and intersection support familiar terminology without duplicating the underlying algorithms. An outline is a closed polygon assembled from sampled normal offsets; constant symmetric widths and arc-length-interpolated asymmetric profiles use the same representation. It is not presented as an exact polynomial offset curve.

Practical Checklist

• Evaluate with de Casteljau when stability and arbitrary degree matter.
• Evaluate rational curves in homogeneous coordinates and divide only after the final interpolation.
• Use quadratic or cubic segments for editable paths, not one very high-degree curve.
• Evaluate tensor-product surfaces by applying de Casteljau successively in both parameter directions.
• Split adaptively and measure error in the coordinate space that matters to the output.
• Find exact axis-aligned bounds from derivative roots.
• Use numerical quadrature for length and a lookup table for repeated constant-speed traversal.
• Treat projection and intersections as root-isolation problems with explicit tolerances.
• Remember that offsets and ordinary circular arcs are not polynomial Bézier curves.
• Distinguish geometric continuity from parametric continuity at joins.

References

• [deCasteljau1959]Paul de Casteljau, Outillages méthodes calcul, André Citroën internal report, 1959.
• [Bezier1972]Pierre Bézier, Numerical Control: Mathematics and Applications, John Wiley & Sons, 1972.
• [Farin2002]Gerald Farin, Curves and Surfaces for CAGD: A Practical Guide, fifth edition, Morgan Kaufmann, 2002.
• [SederbergNishita1990]Thomas W. Sederberg and Tomoyuki Nishita, Curve Intersection Using Bézier Clipping, Computer-Aided Design 22(9), 1990.
• [Yap2006]Chee K. Yap, Complete Subdivision Algorithms, I: Intersection of Bézier Curves, Proceedings of the 22nd Annual Symposium on Computational Geometry, 2006.
• [StoneDeRose1989]Maureen C. Stone and Tony D. DeRose, A Geometric Characterization of Parametric Cubic Curves, ACM Transactions on Graphics 8(3), 1989.
• [WangEtAl2008]Wenping Wang, Bert Jüttler, Dayue Zheng, and Yang Liu, Computation of Rotation Minimizing Frames, ACM Transactions on Graphics 27(1), 2008.