De Casteljau's algorithm
From Free net encyclopedia
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 Bézier curves.
Although the algorithm is slower for most architectures when compared with the direct approach it is numerically more stable.
Contents |
[edit]
Definition
Given a polynomial B in Bernstein form of degree n
- <math>B(t) = \sum_{i=0}^{n}\beta_{i}b_{i,n}(t),</math>
where b is a Bernstein basis polynomial, the polynomial at point t0 can be evaluated with the recurrence relation
- <math>\beta_i^{(0)} := \beta_i \mbox{ , } i=0,\ldots,n</math>
- <math>\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</math>
with
- <math>B(t_0)=\beta_0^{(n)}</math>.
[edit]
Notes
When doing the calculation by hand it is useful to write down the coefficients in a triangle scheme as
- <math>
\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} </math>
When choosing a point t0 to evaluate a Bernstein polynomial we can use the two diagonals of the triangle scheme to construct a division of the polynomial
- <math>B(t) = \sum_{i=0}^n \beta_i^{(0)} b_{i,n}(t) \qquad \mbox{ , } \in [0,1]</math>
into
- <math>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]</math>
and
- <math>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]</math>
[edit]
Example
We want to evaluate the Bernstein polynomial of degree 2 with the Bernstein coefficients
- <math>\beta_0^{(0)} = \beta_0</math>
- <math>\beta_1^{(0)} = \beta_1</math>
- <math>\beta_2^{(0)} = \beta_2</math>
at the point t0.
We start the recursion with
- <math>\beta_0^{(1)} = \beta_0^{(0)} (1-t_0) + \beta_1^{(0)}t_0 = \beta_0(1-t_0) + \beta_1 t_0</math>
- <math>\beta_1^{(1)} = \beta_1^{(0)} (1-t_0) + \beta_2^{(0)}t_0 = \beta_1(1-t_0) + \beta_2 t_0</math>
and with the second iteration the recursion stops with
- <math>
\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} </math>
which is the expected Bernstein polynomial of degree n.
[edit]
Bézier curve
When evaluating a Bézier curve of degree n in 3 dimensional space with n+1 control points Pi
- <math>\mathbf{B}(t) = \sum_{i=0}^{n} \mathbf{P}_i b_{i,n}(t) \mbox{ , } t \in [0,1]</math>
with
- <math>\mathbf{P}_i :=
\begin{pmatrix} x_i \\ y_i \\ z_i \end{pmatrix} </math>.
we split the Bézier curve into three separate equations
- <math>B_1(t) = \sum_{i=0}^{n} x_i b_{i,n}(t) \mbox{ , } t \in [0,1]</math>
- <math>B_2(t) = \sum_{i=0}^{n} y_i b_{i,n}(t) \mbox{ , } t \in [0,1]</math>
- <math>B_3(t) = \sum_{i=0}^{n} z_i b_{i,n}(t) \mbox{ , } t \in [0,1]</math>
which we evaluate individually using de Casteljau's algorithm.
[edit]
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
[edit]
Sample implementations
[edit]
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:[]) = [linearInterp 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 = [deCas (fromIntegral x / fromIntegral n) b | x<-[0 .. n] ]
[edit]
Python
(This code requires the Python Imaging Library.)
import Image
import ImageDraw
SIZE=128
image = Image.new("RGB", (SIZE, SIZE))
d = ImageDraw.Draw(image)
def midpoint((x1, y1), (x2, y2)):
return ((x1+x2)/2, (y1+y2)/2)
MAX_LEVEL = 5
def draw_curve(P1, P2, P3, P4, level=1):
if level == MAX_LEVEL:
d.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)
draw_curve((10.0,10.0),(100.0,100.0),(100.0,10.0),(100.0,100.0))
image.save("deCasteljau.png", "PNG")
print "ok."
[edit]
Object Pascal/Delphi
{Simple Bezier drawing code using de Casteljau's algorithm, by Roy Smeding.
Take this out or whatever you want, it's not like I care, I just thought it
would be nice for you to know who created this, just for the record.
A possible optimalization would be to use threads, only I don't know anything
about thread programming, so convert it if you want.
Note:
(This code was only tested with the original Delphi IDE/compiler. The form has
a TPaintBox(TImage works, too) called PaintBox and a button called CalcButton.)}
{--Standard header code generated by Delphi, not needed to understand this.--}
var
MainForm: TMainForm;
MaxLevel: Integer = 5; {Maximum iteration level}
implementation
{$R *.dfm}
{Declaration of a more precise TPoint type, to decrease errors because of rounding.}
type TPointDouble = record
X,Y: Double;
end;
function MidPoint(P1, P2: TPointDouble): TPointDouble;
begin
{This calculates the point in between two points. Rather simple.}
MidPoint.X:=(P1.X+P2.X)/2;
MidPoint.Y:=(P1.Y+P2.Y)/2;
end;
procedure DrawCurve(P1,P2,P3,P4: TPointDouble; Level: Integer);
var
L1, L2, L3, L4, R1, R2, R3, R4, H: TPointDouble;
begin
{This draws the actual curve. It very much resembles the pseudocode described above.}
if Level>MaxLevel then begin
MainForm.PaintBox.Canvas.Pen.Color:=clBlack;
{In this wee bit we most definitely need rounding, to convert from the
floating point numbers to integers needed for MoveTo and LineTo.}
MainForm.PaintBox.Canvas.MoveTo(Round(P1.X),Round(P1.Y));
MainForm.PaintBox.Canvas.LineTo(Round(P4.X),Round(P4.Y));
end else begin
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;
DrawCurve(L1,L2,L3,L4,level+1);
DrawCurve(R1,R2,R3,R4,level+1);
Application.ProcessMessages;
end;
end;
procedure TMainForm.CalcButtonClick(Sender: TObject);
begin
{What happens when CalcButton is clicked, this simply starts drawing the curve.}
PaintBox.Refresh;
DrawCurve(Point(3,3),Point(3,PaintBox.Height-3),
Point(PaintBox.Width-3,PaintBox.Height-3),Point(PaintBox.Width-3,3),0);
end;
[edit]
[[C++]] / OpenGL
#include <utility> // std::pair<>
#include <GL/gl.h>
typedef std::pair<double,double> Point;
inline Point midpoint(const Point& p1, const Point& p2)
{
return Point((p1.first+p2.first)/2, (p1.second+p2.second)/2);
}
inline void draw_line(const Point& p1, const Point& p2)
{
glBegin(GL_LINES);
glVertex2d(p1.first, p1.second);
glVertex2d(p2.first, p2.second);
glEnd();
}
void draw_curve(const Point& p1, const Point& p2, const Point& p3, const Point& p4, int level=5)
{
if (level == 0)
draw_line(p1, p4);
else {
const Point& L1 = p1;
const Point& L2 = midpoint(p1, p2);
const Point& H = midpoint(p2, p3);
const Point& R3 = midpoint(p3, p4);
const Point& R4 = p4;
const Point& L3 = midpoint(L2, H);
const Point& R2 = midpoint(R3, H);
const Point& L4 = midpoint(L3, R2);
const Point& R1 = L4;
draw_curve(L1, L2, L3, L4, level-1);
draw_curve(R1, R2, R3, R4, level-1);
}
}
...
draw_curve( Point(10, 10),
Point(100, 100),
Point(100, 10),
Point(100, 100));
[edit]
ARM
AREA bezier, CODE ENTRY Start ldfs F0, X0 stfs F0, [SP, #-4]! ldfs F0, Y0 stfs F0, [SP, #-4]! ldfs F0, X1 stfs F0, [SP, #-4]! ldfs F0, Y1 stfs F0, [SP, #-4]! ldfs F0, X2 stfs F0, [SP, #-4]! ldfs F0, Y2 stfs F0, [SP, #-4]! ldfs F0, X3 stfs F0, [SP, #-4]! ldfs F0, Y3 stfs F0, [SP, #-4]! bl DeCasteljau add SP, SP, #32 swi &11 ;----------------------------------------------------- X0 DCFS 0 Y0 DCFS 0 X1 DCFS 10 Y1 DCFS 10 X2 DCFS 5 Y2 DCFS 2 X3 DCFS 10 Y3 DCFS 1 Thresh DCFS 2 ;----------------------------------------------------- sX0 EQU 36 sY0 EQU 32 sX1 EQU 28 sY1 EQU 24 sX2 EQU 20 sY2 EQU 16 sX3 EQU 12 sY3 EQU 8 DeCasteljau str LR, [SP, #-4]! ; +4 str ip, [SP, #-4]! ; +0 mov ip, SP stfs F0, [SP, #-4]! ; -4 stfs F1, [SP, #-4]! ; -8 stfs F2, [SP, #-4]! ; -12 stfs F3, [SP, #-4]! ; -16 ldfs F0, [ip, #sX0] ldfs F1, [ip, #sY0] ldfs F2, [ip, #sX3] ldfs F3, [ip, #sY3] bl Distance ldfs F1, Thresh cmf F0, F1 bge iterjau ldfs F0, [ip, #sX0] ldfs F1, [ip, #sY0] bl Plot_Line ldfs F3, [SP], #4 ldfs F2, [SP], #4 ldfs F1, [SP], #4 ldfs F0, [SP], #4 ldr ip, [SP], #4 ldr PC, [SP], #4 iterjau stfs F4, [SP, #-4]! ; -20 stfs F5, [SP, #-4]! ; -24 stfs F6, [SP, #-4]! ; -28 stfs F7, [SP, #-4]! ; -32 sub SP, SP, #16 ; -36, -40, -44, -48 ldfs F0, [ip, #sX0] ldfs F1, [ip, #sX1] adfs F0, F0, F1 ; P01x dvfs F7, F0, #2 ; X0 + X1 / 2 ldfs F0, [ip, #sY0] ldfs F1, [ip, #sY1] adfs F0, F0, F1 ; P01y dvfs F6, F0, #2 ; Y0 + Y1 / 2 ldfs F0, [ip, #sX1] ldfs F1, [ip, #sX2] adfs F0, F0, F1 ; P12x dvfs F5, F0, #2 ; X1 + X2 / 2 ldfs F0, [ip, #sY1] ldfs F1, [ip, #sY2] adfs F0, F0, F1 ; P12y dvfs F4, F0, #2 ; Y1 + Y2 / 2 ldfs F0, [ip, #sX2] ldfs F1, [ip, #sX3] adfs F0, F0, F1 ; P23x dvfs F3, F0, #2 ; X2 + X3 / 2 ldfs F0, [ip, #sY2] ldfs F1, [ip, #sY3] adfs F0, F0, F1 ; P23y dvfs F2, F0, #2 ; Y2 + Y3 / 2 adfs F1, F7, F5 ; P012x dvfs F1, F1, #2 ; P01x + P12x / 2 adfs F0, F6, F4 ; P012y dvfs F0, F0, #2 ; P01y + P12y / 2 adfs F5, F5, F3 ; P123x dvfs F5, F5, #2 ; P12x + P23x / 2 adfs F4, F4, F2 ; P123y dvfs F4, F4, #2 ; P12y + P23y / 2 stfs F5, [ip, #-36] stfs F4, [ip, #-40] stfs F3, [ip, #-44] stfs F2, [ip, #-48] adfs F3, F1, F5 ; P0123x dvfs F3, F3, #2 ; P012x + P123x / 2 adfs F2, F0, F4 ; P0123y dvfs F2, F2, #2 ; P012y + P123y / 2 ldfs F5, [ip, #sX0] ldfs F4, [ip, #sY0] stfs F5, [SP, #-4]! stfs F4, [SP, #-4]! stfs F7, [SP, #-4]! stfs F6, [SP, #-4]! stfs F1, [SP, #-4]! stfs F0, [SP, #-4]! stfs F3, [SP, #-4]! stfs F2, [SP, #-4]! bl DeCasteljau ; don't add to SP.. ldfs F5, [ip, #-36] ldfs F4, [ip, #-40] ldfs F7, [ip, #-44] ldfs F6, [ip, #-48] ldfs F1, [ip, #sX3] ldfs F0, [ip, #sY3] stfs F3, [SP, #28] ; ..just change the stfs F2, [SP, #24] ; values on the stack. stfs F5, [SP, #20] stfs F4, [SP, #16] stfs F7, [SP, #12] stfs F6, [SP, #8] stfs F1, [SP, #4] stfs F0, [SP] bl DeCasteljau add SP, SP, #48 ldfs F7, [SP], #4 ldfs F6, [SP], #4 ldfs F5, [SP], #4 ldfs F4, [SP], #4 ldfs F3, [SP], #4 ldfs F2, [SP], #4 ldfs F1, [SP], #4 ldfs F0, [SP], #4 ldr ip, [SP], #4 ldr PC, [SP], #4 ;----------------------------------------------------- Plot_Line str LR, [SP, #-4]! str R0, [SP, #-4]! mov R0, #'<' swi &0 fix R0, F0 str R0, [SP, #-4]! bl Print_Number mov R0, #',' swi &0 fix R0, F1 str R0, [SP] bl Print_Number mov R0, #'>' swi &0 mov R0, #' ' swi &0 mov R0, #'-' swi &0 swi &0 mov R0, #' ' swi &0 mov R0, #'<' swi &0 fix R0, F2 str R0, [SP] bl Print_Number mov R0, #',' swi &0 fix R0, F3 str R0, [SP] bl Print_Number mov R0, #'>' swi &0 mov R0, #10 swi &0 add SP, SP, #4 ldr R0, [SP], #4 ldr PC, [SP], #4 ;----------------------------------------------------- Distance stfs F1, [SP, #-4]! sufs F0, F2, F0 sufs F1, F3, F1 mufs F1, F1, F1 mufs F0, F0, F0 adfs F0, F0, F1 sqts F0, F0 ldfs F1, [SP], #4 mov PC, LR ;----------------------------------------------------- Print_Number str LR, [SP, #-4]! ; +20 str ip, [SP, #-4]! ; +16 str R0, [SP, #-4]! ; +12 str R1, [SP, #-4]! ; +8 str R2, [SP, #-4]! ; +4 str R3, [SP, #-4]! ; +0 mov ip, SP ldr R0, [ip, #24] ; number, via stack mov R3, #0 tst R0, #2,2 beq skipNegLp sub R1, R0, #1 mov R0, #'-' swi &0 mvn R0, R1 skipNegLp add R3, R3, #1 bl UDiv10 str R1, [SP, #-4]! cmp R0, #0 bne skipNegLp pDigitsLp ldr R0, [SP], #4 add R0, R0, #'0' swi &0 sub R3, R3, #1 cmp R3, #0 bne pDigitsLp ldr R3, [SP], #4 ldr R2, [SP], #4 ldr R1, [SP], #4 ldr R0, [SP], #4 ldr ip, [SP], #4 ldr PC, [SP], #4 ;----------------------------------------------------- UDiv10 sub R1, R0, #10 sub R0, R0, R0, LSR #2 add R0, R0, R0, LSR #4 add R0, R0, R0, LSR #8 add R0, R0, R0, LSR #16 mov R0, R0, LSR #3 add R2, R0, R0, LSL #2 subs R1, R1, R2, LSL #1 addpl R0, R0, #1 addmi R1, R1, #10 mov PC, LR
[edit]
References
- Farin, Gerald & Hansford, Dianne (2000). The Essentials of CAGD. Natic, MA: A K Peters, Ltd. ISBN 1-56881-123-3
[edit]
See also
- Horner scheme to evaluate polynomials in Monomial form
- Clenshaw algorithm to evaluate polynomials in Chebyshev formcs:Racionální Algoritmus de Casteljau
de:De Casteljau-Algorithmus fr:Algorithme de De Casteljau zh:De Casteljau算法