de casteljau's algorithm

In the mathematical subfield of numerical analysis the de Casteljau's algorithm, named after its inventor Paul de Casteljau, is a recursive method to evaluate polynomials in Bernstein form or Bzier curves. Although the algorithm is slower for most architectures when compared with the direct approach it is numerically more stable.

Definition

Given a polynomial B in Bernstein form of degree n

B(t) = \sum_{i=0}^{n}\beta_{i}b_{i,n}(t)

the polynomial at point t₀ can be evaluated with recurrence relation

\beta_i^{(0)} := \beta_i \mbox{ , } i=0,\ldots,n

\beta_i^{(j)} := \beta_i^{(j-1)} (1-t_0) + \beta_{i+1}^{(j-1)} t_0 \mbox{ , } i = 0,\ldots,n-j \mbox{ , } j= 1,\ldots n

with

B(t_0)=\beta_0^{(n)}

Notes

When doing the calculation by hand it is useful to write down the coefficients in a triangle scheme as

\begin{matrix} \beta_0 & = \beta_0^{(0)} & & & \\

             &                     & \beta_0^{(1)}     &         &               \\

\beta_1 & = \beta_1^{(0)} & & & \\

             &                     &                   & \ddots  &               \\

\vdots & & \vdots & & \beta_0^{(n)} \\

             &                     &                   &         &               \\

\beta_{n-1} & = \beta_{n-1}^{(0)} & & & \\

             &                     & \beta_{n-1}^{(1)} &         &               \\

\beta_n & = \beta_n^{(0)} & & & \\ \end{matrix} When choosing a point t₀ to evaluate a Bernstein polynomial we can use the two diagonals of the triangle scheme to construct a division of the polynomial

B(t) = \sum_{i=0}^n \beta_i^{(0)} b_{i,n}(t) \qquad \mbox{ , } \in 0,1

into

B_1(t) = \sum_{i=0}^n \beta_0^{(i)} b_{i,n}\left(\frac{t}{t_0}\right) \qquad \mbox{ , } \in 0,t_0

and

B_2(t) = \sum_{i=0}^n \beta_{n-i}^{(i)} b_{i,n}\left(\frac{t-t_0}{1-t_0}\right) \qquad \mbox{ , } \in t_0,1

Example

We want to evaluate the Bernstein polynomial of degree 2 with the Bernstein coefficients

\beta_0^{(0)} = \beta_0

\beta_1^{(0)} = \beta_1

\beta_2^{(0)} = \beta_2

at the point t₀. We start the recursion with

\beta_0^{(1)} = \beta_0^{(0)} (1-t_0) + \beta_1^{(0)}t = \beta_0(1-t_0) + \beta_1 t_0

\beta_1^{(1)} = \beta_1^{(0)} (1-t_0) + \beta_2^{(0)}t = \beta_1(1-t_0) + \beta_2 t_0

and with the second iteration the recursion stops with

\begin{matrix} \beta_0^{(2)} & = & \beta_0^{(1)} (1-t_0) + \beta_1^{(1)} t_0 \\ \ & = & \beta_0(1-t_0) (1-t_0) + \beta_1 t_0 (1-t_0) + \beta_1(1-t_0)t_0 + \beta_2 t_0 t_0 \\ \ & = & \beta_0 (1-t_0)^2 + \beta_1 2t_0(1-t_0) + \beta_2 t_0^2 \\ \end{matrix} which is the expected Bernstein polynomial of degree n.

Bzier curve

When evaluating a Bzier curve of degree n in 3 dimensional space with n+1 control points P_i

\mathbf{B}(t) = \sum_{i=0}^{n} \mathbf{P}_i b_{i,n}(t) \mbox{ , } t \in 0,1

with

\mathbf{P}_i :=

\begin{pmatrix} x_i \\ y_i \\ z_i \end{pmatrix} . we split the Bzier curve into three separate equations

B_1(t) = \sum_{i=0}^{n} x_i b_{i,n}(t) \mbox{ , } t \in 0,1

B_2(t) = \sum_{i=0}^{n} y_i b_{i,n}(t) \mbox{ , } t \in 0,1

B_3(t) = \sum_{i=0}^{n} z_i b_{i,n}(t) \mbox{ , } t \in 0,1

which we evaluate individually using de Casteljau's algorithm.

Pseudocode example

Here is a pseudo-code example that recursively draws a curve from point P1 to P4, tending towards points P2 and P3. The level parameter is the limiting factor of the recursion. The procedure calls itself recursively with an incremented level parameter. When level reaches the max_level global value, a line is drawn between P1 and P4 at that level. The function midpoint takes two points, and returns the middle point of a line segment between those two points.

      global max_level = 5     procedure draw_curve(P1, P2, P3, P4, level)         if (level > max_level)             draw_line(P1, P4)         else             L1 = P1             L2 = midpoint(P1, P2)             H  = midpoint(P2, P3)             R3 = midpoint(P3, P4)             R4 = P4             L3 = midpoint(L2, H)             R2 = midpoint(R3, H)             L4 = midpoint(L3, R2)             R1 = L4             draw_curve(L1, L2, L3, L4, level + 1)             draw_curve(R1, R2, R3, R4, level + 1)     end procedure draw_curve

Sample implementations

Haskell

 Compute a point R by performing a linear interpolation of points P and Q with parameter t. Usage: linearInterp P Q t            P = list representing a point           Q = list representing a point           t = parametric value for the linear interpolation, t<-0..1 
 Returns: List representing the point (1-t)P + tQ  >	linearInterp :: Float->Float->Float->Float >	linearInterp [] [] _ = [] >	linearInterp (p:ps) (q:qs) t = (1-t)*p + t*q : linearInterp ps qs t   Compute one level of intermediate results by performing a linear interpolation between each pair of control points. Usage: eval t b            t = parametric value for the linear interpolation, t<-0..1           b = list of control points 
 Returns:  For n control points, returns a list of n-1 interpolated points  >	eval :: Float->Float->Float >	eval t (bi:bj:[]) = bi bj t >	eval t (bi:bj:bs) = (linearInterp bi bj t) : eval t (bj:bs)   Compute a point on a Bezier curve using the de Casteljau algorithm. Usage: deCas t b            t = parametric value for the linear interpolation, t<-0..1           b = list of control points 
 Returns:  List representing a single point on the Bezier curve  >	deCas :: Float->Float->Float >	deCas t (bi:[]) = bi >	deCas t bs = deCas t (eval t bs)   Compute the points along a Bezier curve using the de Casteljau algorithm.  Points are returned in a list. Usage: bezierCurve n b           n = number of points along the curve to compute          b = list of Bezier control points 
 Returns: A list of n+1 points on the Bezier curve Example: bezierCurve 50 [[0,0],[1,1],[0,1],[1,0]]  >	bezierCurve :: Int->Float->Float >	bezierCurve n b = (fromIntegral x / fromIntegral n) b | x<-[0 .. n ]

References

Farin, Gerald & Hansford, Dianne (2000). The Essentials of CAGD. Natic, MA: A K Peters, Ltd. ISBN 1-56881-123-3