Bezier Curves

An N-dimensional Bezier curve of order M is a curve of the form p(s) = åi=0M  MCi si(1-s)M-ipi where the control points p0,...,pM are N-dimensional vectors. Of greatest interest to programmers are Beziers of order 3. These are given by
 p(s) = (1-s)3p0 + 3s(1-s)2p1 + 3s2(1-s)p2 + s3p3 = s3a3 + s2a2 + sa1 + a0
where a0 = p0 ; a1 = -3p0+3p1 ; a2 = 3p0-6p1+3p2 ; a3 = -p0+3p1-3p2+p3 are the coefficient points of the bezier.

Subdividing Beziers
A Bezier curve with control points p0,p1,p2,p3 can be split (at s=½) into two Beziers having control points p00,p01,p02,p03 and  p10,p11,p12,p13 where
p00 = p(0)=p0
p01 = (p0+p1)/2
p02 = (p0+2p1+p2)/4
p03 = p10 = p(½) = (p0+3(p1+p2)+p3)/8
p11 = (p3+2p2+p1)/4
p12 = (p2+p3)/2
p13 = p(1)=p3.

To split at a general s=t, it is easier to manipulate coefficient rather than control points. We can construct a bezier for the [0,t] portion by substituting s=ut to obtain coefficients (a0, a1t, a2t2, a3t3). The [t,1] portion is given by the coefficients of u in a(t + u(1-t)).

Bezier Approximations

To Cubic
The cubic a3s3 + a2s2 + a1s + a0 is precisely equivalent to a Bezier with control points
p3 = (a0+a1+a2+a3)/4 ; p2 = (3a0+2a1+a2)/3 ; p1 = (3a0+a1)/3 ; p0 = a1 .

To Circular Arc
With acknowledgment to the work of David Seal
Suppose we wish to construct a 2D bezier curve   which approximates the unit circular arc from ( cosa,- sina) to ( cosa, sina).
If we impose a second order fit at the arc endpoints then our bezier must have control points
p0 = ( cosa,- sina) ; p1 = p0 + b( sina, cosa) ; p2 = p3 + b( sina,- cosa) ; p3 = ( cosa, sina) ;
whence the curve is given by
 p1(s) = 3b sinas(1-s) + cosa p2(s) = (-4 sina +6b cosa)s3   + (6 sina - 9b cosa)s2 + 3b cosas - sina = (2s-1)((2 sina-3b cosa)s(1-s) + sina) = (2s-1)((1+2s(1-s)) sina -3b cosas(1-s)) Writing z = s(1-s) p1(s) = 3b sinaz + cosa p2(s) = (2s-1)((1+2z) sina -3b cosaz)

Consider the error terms g2(s) = p(s)2 - 1 and g(s) = |p(s)| - 1.
g2(s)= p(s)2-1 = (|p(s)|+1)g(s) = (2+g(s))g(s) whence g(s) » g2(s)/2.
 g2(s) = (3b sinaz + cosa)2 +   (2s-1)2((1+2z) sina -3b cosaz)2 -1 = (9b2sin2az2 + cos2a) -1 + (1-4z)((1+2z)2sin2a -6b cosaz(1+2z) sina +9b2cos2az2) = (9b2sin2az2 + cos2a) -1 + (1-4z)(z2(4sin2a - 12b cosa sina + 9b2cos2a) +z(4sin2a -6b cosa sina) + sin2a) = z3(48b cosa sina - 16sin2a - 36b2cos2a) +z2( 9b2sin2a + (4sin2a - 12b cosa sina + 9b2cos2a) - 4 ( 4 sin2a -6b cosa sina )) = -4z3(3b cosa-2 sina)2 + z2(9b2 + 12b cosa sina - 12sin2a)

 dg2/ds = (1-2s)(-12z2(3b cosa-2 sina)2 +2z(9b2 + 12b cosa sina - 12sin2a)) = (1-2s)z(-12z(3b cosa-2 sina)2 +2(9b2 + 12b cosa sina - 12sin2a)) = 6(1-2s)z(-2z(3b cosa-2 sina)2 + (3b2 + 4b cosa sina - 4sin2a))

Hence dg2/ds = 0 when s=0,1, ½, or when z = (3b2 + 4b cosa sina - 4sin2a))/(-2z(3b cosa-2 sina)2.     These give g2(s) = 0, 0, (-(3b cosa-2 sina)2+(9b2+12b cosa sina-12sin2a))/16 = (9b2sin2a+24b cosa sina-16sin2a)/16, and
(3b2+4b cosa sina-4sin2a)3/(4(3b cosa-2 sina)4) repsectively.

We can ensure g2(½)=0 by choosing b0 such that p1(½)=1 Þ 3b0 sina½(1-½) + cosa = 1 Þ b0 = 4(1- cosa)/(3 sina).
The maximal error is then
 g2max = (3b02+4b0 cosa sina-4sin2a)3/ (4(3b0 cosa-2 sina)4) = (16(1- cosa)2/(3sin2a)+16(1- cosa) cosa/3-4sin2a)3 / (4(4(1- cosa) cosa/ sina-2 sina)4) = ... = (1- cosa)3/(27(1+ cosa)) Þ gmax = (1- cosa)3/(54(1+ cosa)) » (1- cosa)3/108   when a small

 a gmax Pixel radius for > ½ pixel error p/2 0.019 27 p/3 .0015 324 p/4 2.8 x 10-4 1831 p/8 4.24 x 10-6 117770 p/16 6.6 x 10-8 7538744 d d6/864 432 d-6

To General 1-D Curve
Suppose we wish to construct a 1D bezier curve q(s)=s3q3 + s2q2 + sq1 + q0 with given values q(0), q(1) and gradients q'(0), q'(1) .
 æ 1 1 1 1 ö æ q0 ö = æ q(1) ö ç 0 1 2 3 ÷ ç q1 ÷ ç q'(1) ÷ ç 0 1 0 0 ÷ ç q2 ÷ ç q'(0) ÷ è 1 0 0 0 ø è q3 ø è q(0) ø
is easily solved by q0 = q(0) ; q1 = q'(0) ; q2 = 1/3(q'(1)-2q1-q0) = 1/3(q'(1)-2q'(0)-q(0)) = ; q3 = q(0)-q0-q1-q2 = -(q1+q2) = -1/3(q'(1)+q'(0)-q(0)) .

To General N-D Curve
Suppose we wish to construct a 2D bezier curve q(s)=s3q3 + s2q2 + sq1 + q0 through given points q(0), q(t) and q(1).
Suppose we also   require q'(0) = m0 g0 , q'(1)= m1g1 where g0 and g1 are given. If g0 and g1 are non parallel this may be achieved. If g0 and g1 are (close to) parallel then a (viable) fit is only possible if q(0), q(t) and q(1) are colinear and in accordance with the desired gradient.
 æ 1 t t2 t3 ö æ q0 ö = æ q(t) ö ç 1 1 1 1 ÷ ç q1 ÷ ç q(1) ÷ ç 0 1 2 3 ÷ ç q2 ÷ ç m1g1 ÷ ç 0 1 0 0 ÷ è q3 ø ç m0g0 ÷ è 1 0 0 0 ø è q(0) ø

Substituting q0 = q(0), then q1 = m0g0, and finally q2 = q(1)-q0-m0g0 - q3 we obtain:
 æ t t2 t3 ö æ q1 ö = æ q(t) - q0 ö ç 1 1 1 ÷ ç q2 ÷ ç q(1)-q0 ÷ ç 1 2 3 ÷ è q3 ø ç m1g1 ÷ è 1 0 0 ø è m0g0 ø

Substituting q1 = m0g0 :
 æ t2 t3 ö æ q2 ö = æ q(t) - q0 - tm0g0 ö ç 1 1 ÷ è q3 ø ç q(1)-q0-m0g0 ÷ è 2 3 ø è m1g1-m0g0 ø

Substituting q2 = q(1)-q0-m0g0 - q3 and then eliminating q3:
 t2(q(1)-q0-m0g0 - q3) + t3q3 = q(t) - q0 - tm0g0 2(q(1)-q0-m0g0 - q3) + 3q3 = m1g1-m0g0 Þ m0g0(1-t)/t - m1g1 = q(t)/t2(1-t) + q(1)(2t-3)/(1-t) - q0(2t3 - 3t2 +1)/t2(1-t))

If t is known and g0 and g1 are not parallel this may be solved for m0,m1 whereupon we can set

• q0 = q(0)
• q1 = m0g0
• q3 = m1g1 + q1 -2q(1) +2q0
• q2 = q(1) - q3 - q1 - q0

What is a sound choice for t ?
t =½ is an obvious candidate, the equations for m0 and m1 becoming
m0g0 - m1g1 = 8q(½) - 4q(1) - 4q0
The point of least curvature is given by minimising q"(t)2 = (6q3t + 2q2)2 which is minimised by t = -q2.q3/3 so if this lies comfortably inside (0,1) it may be a good choice. Please let me know if progress is made in this direction.

Glossary   Contents   Author
Copyright (c) Ian C G Bell 1998
Web Source: www.iancgbell.clara.net/maths or www.bigfoot.com/~iancgbell/maths
Latest Edit: 09 Jun 2007.