Bézier subdivision

From CGAFaq

A Bézier curve is a parametric polynomial function from the interval [0..1] to a space, usually 2D or 3D. Common Bézier curves use cubic polynomials, so have the form

f(t) = a₃ t³ + a₂ t² + a₁ t + a₀,

where the coefficients are points in 3D. Blossoming converts this polynomial to a more helpful form. Let s = 1-t, and think of tri-linear interpolation:

F([s0,t0],[s1,t1],[s2,t2])	=	s0(s1(s2 c30 + t2 c21) + t1(s2 c21 + t2 c12)) +
		t0(s1(s2 c21 + t2 c12) + t1(s2 c12 + t2 c03))
	=	c30(s0 s1 s2) +
		c21(s0 s1 t2 + s0 t1 s2 + t0 s1 s2) +
		c12(s0 t1 t2 + t0 s1 t2 + t0 t1 s2) +
		c03(t0 t1 t2).

The data points c₃₀, c₂₁, c₁₂, and c₀₃ have been used in such a way as to make the resulting function give the same value if any two arguments, say [s₀,t₀] and [s₂,t₂], are swapped. When [1-t,t] is used for all three arguments, the result is the cubic Bézier curve with those data points controlling it:

f(t)	=	F([1-t,t],[1-t,t],[1-t,t])
	=	(1-t)³ c30 +
		3(1-t)² t c21 +
		3(1-t) t² c12 +
		t³ c03

Notice that

F([1,0],[1,0],[1,0])	=	c30
F([1,0],[1,0],[0,1])	=	c21
F([1,0],[0,1],[0,1])	=	c12
F([0,1],[0,1],[0,1])	=	c03

In other words, c_ij is obtained by giving F argument t's i of which are 0 and j of which are 1. Since F is a linear polynomial in each argument, we can find f(t) using a series of simple steps. Begin with

f₀₀₀ = c₃₀, f₀₀₁ = c₂₁, f₀₁₁ = c₁₂, f₁₁₁ = c₀₃.

Then compute in succession:

f00t = s f000 + t f001	f01t = s f001 + t f011	f11t = s f011 + t f111
f0tt = s f00t + t f01t	f1tt = s f01t + t f11t
fttt = s f0tt + t f1tt

This is the de Casteljau algorithm for computing f(t) = f_ttt.

It also has split the curve for the intervals [0..t] and [t..1]. The control points for the first interval are f₀₀₀, f_00t, f_0tt, f_ttt, while those for the second interval are f_ttt, f_1tt, f_11t, f₁₁₁.

If you evaluate 3⋅F([1-t,t],[1-t,t],[-1,1]) you will get the derivative of f at t. Similarly, 2⋅3⋅F([1-t,t],[-1,1],[-1,1]) gives the second derivative of f at t, and finally using 1⋅2⋅3⋅F([-1,1],[-1,1],[-1,1]) gives the third derivative.

This procedure is easily generalized to different degrees, triangular patches, and B-spline curves.