Parametric Cubic Curves

In computer graphics you'll often need to draw curves, and there are times when a line list just won't cut it.  What you want is something continuous.   Something that is differentiable.  Something that can be rendered with infinite precision.  The cubic curve, my friend, is just the beast you seek.

The parametric function is, of course, our weapon of choice.   Slope-intercept and general form equations have an annoying attribute that make them impossible to work with: namely, the geometric slope of a function may be infinite, creating a mathematical singularity.  Parametric curves, however, never have infinite slopes.  Our curve is described by a piecewise polynomial, that is, by a function x, a function y, and a function z.   In a slope-intercept form function, the slope of a curve is defined by the infamous three words: "rise over run."  When the "run" is zero, a dimensional vortex opens and sucks you away to the outer limits.   Parametric functions, on the other hand, define slope as the sum of the parametric velocities: dx/dt i + dy/dt j + dz/dt k.   Because we are limiting ourselves to parametric polynomials, these derivatives are always defined.  Oh, happy day.

Our cubic curve can exist in any dimension.  A one dimensional curve looks like a line segment; that isn't too interesting.  A parametric curve defined with four or more piecewise polynomials is impossible to draw.  Applications always use two and three dimensional cubic curves.

We define our cubic curve as the function Q(t).  This function is in turn defined as a set of polynomials:

Q(t) = [ f₁(t)   f₂(t)   ...   f_n-1(t)   f_n(t) ]    , n = dimension of curve

f₁(t) = a₁t³ + b₁t² + c₁t + d₁      , 0 £ t £ 1

                    ...

f_n(t) = a_nt³ + b_nt² + c_nt + d_n

We define a time matrix T = [ t³   t²   t   1 ] and a coefficient matrix:

, where n is dimension of the curve.

We can express Q(t) as the product Q(t) = T · C.  We now have the definitions needed to start the cubic curve experience. 

The Hermite Form

The cubic curves discussed here are all controlled by four characteristic points.  The simplest of the cubic curves, the Hermite form, is defined by two endpoints and the tangent vectors at these endpoints.  The following nine curves are defined by a Hermite characteristic matrix.  The points represent the endpoint parameters, and the arrows point in the direction of the tangent vectors.   Note that the arrow lengths do not correspond to the tangent vector magnitudes.

We define the Hermite basis matrix, M_H, which relates the endpoints and tangent vectors on the hermite curve to the polynomial coefficients.  Multiplying the Hermite geometry matrix, G_H, with M_H, effectively transforms G_H from the Hermite geometry to produce the coefficients, C_H, of the parametric function Q_H(t).  Recall that the coefficients of the function Q_H(t) are described by the coefficient matrix, C_H:

C_H = M_H · G_H. 

Our operand matrix, G_H, parameterizes the endpoints and tangent vectors of the Hermite curve.  The matrix takes the form

, where n is dimension of the curve.

The elements P₁ and P₄ are the t = 0 and t = 1 positions of the curve.  R₁ and R₄ are the t = 0 and t = 1 tangent vectors.   This particular notation is used for consistency with later forms.  We will now derive the Hermite basis matrix M_H.

Recall from the previous section that for a generic parametric curve, Q(t) = T · C.  Therefore, Q_H(t) = T · C_H.  We can evaluate the function at t = 0 and t = 1 to find the Hermite basis matrix that transform the P₁ and P₄ rows of the geometry matrix:

Q_H(t) = at³ + bt² + ct + d = T · M_H · G_H = [ t³   t²   t   1 ] · M_H · G_H

Q_H(0) = P₁ = a(0)³ + b(0)² + c(0) + d = T · M_H · G_H = [ 0   0   0   1 ] · M_H · G_H

Q_H(1) = P₄ = a(1)³ + b(1)² + c(1) + d = T · M_H · G_H = [ 1   1   1   1 ] · M_H · G_H

Because the tangent vectors R₁ and R₄ are simply the velocities of the curve at  t = 0 and t = 1, we can differentiate the function Q_H(t) and evaluate at t = 0 and t = 1.

Q_H'(t) = 3at² + 2bt + c = T · M_H · G_H = [ 3t²   2t   1   0 ] · M_H · G_H

Q_H'(0) = R₁ = 3a(0)² + 2b(0) + c = T · M_H · G_H = [ 0   0   1   0 ] · M_H · G_H

Q_H'(1) = R₄ = 3a(1)² + 2b(1) + c = T · M_H · G_H = [ 3   2   1   0 ] · M_H · G_H

These four equalities can be written in matrix form:

We can now easily solve for the basis matrix M_H by reducing the geometry matrix from both sides of the equation and inverting M_H:

.

We now have found the Hermite basis matrix.  The Hermite parametric cubic curve function is:

.

This efficient representation of the Hermite form allows the parametric function coefficients to be calculated in a single matrix multiplication.  Unfortunately the magnitudes of the tangent vectors are somewhat of a mystery.  Often the magnitudes need to be determined by trial and error, as a very large magnitude may be required to perturb the curve a small amount, or a small magnitude perturbs the curve too much.  It is this problem that leads us to a more useful form of the parametric cubic curve: the Bezier form.

The Bezier Form

The Bezier form of the parametric cubic curve solves the problem that plagued the Hermite form: it provides a mechanism for relating the tangent vector magnitudes to points in space.  The endpoints, P₁ and P₄, and the intermediate points, P₂ and P₃, of the Bezier geometry are resolved into the tangent vectors R₁ and R₄ of the Hermite geometry.  We can transform matrices of the Bezier geometry to the Hermite geometry using the Hermite-Bezier change of basis matrix, and then transform from the Hermite basis to polynomial coefficients with the Hermite basis matrix presented in the previous section. 

To help understand the effects of the intermediate control points, we will walk through the evaluation of the Bezier curve function Q_B(t).   When t = 0, the function evaluates to the point P₁.   As t increases, the curve moves towards the other endpoint, P₄, in a straight line (it is, after all, the shortest path between the two points).   However, the intermediate control points deform this straight line.  When t = .01, the value of the intermediate point P₂ has more influence on the curve than the intermediate point P₃.  As t increases to .4, the intermediate control point P₂ still has more control over the curve than P₃.  At t = .5, the intermediate points have equal influence over the curve.  When t > .5, the control point P₃ has the greater influence.  We will get into the details of the intermediate control points in the next section.  Imagine the intermediate control points deforming the line P₁P₄ as you consider the following curves.

We now define the relevant Bezier form matrices.   The Bezier geometry matrix is similar to the Hermite geometry matrix, but it takes four control points, rather than two control points and two tangent vectors:

, where n is dimension of the curve.

We also define the Hermite-Bezier change of basis: G_H = M_HB · G_B. 

To derive the Hermite-Bezier change of basis matrix, consider a curve in one dimension with constant velocity.  The Bezier curve exposes an interesting property for this scenario: the control points are equally spaced.   This graph shows a one dimensional Bezier curve with constant velocity.

We can use this graph to find the Hermite-Bezier change of basis matrix.  This curve's function is obvious: Q(t) = (P₄-P₁)t + P₁.  The coefficient matrix follows: C = [ 0   0   P₄-P₁    P₁ ]^T.  We can express the coefficient matrix as the product of the Hermite basis matrix and a Hermite geometry matrix.  We can solve for the Hermite geometry matrix to find a relation between the Bezier intermediate control points and the Hermite tangent vectors:

.

Because P₄ - P₁ = 3l and P₂ - P₁ = l, R₁ = 3(P₂ - P₁).   Likewise, since P₄ - P₃ = l, R₄ = 3(P₄ - P₃).  We now have the four equations that relate the Hermite geometry to the Bezier geometry:

P₁ = P₁

P₄ = P₄

R₁ = -3P₁ + 3P₂

R₄ = -3P₃ + 3P₄

We now present the Hermite-Bezier change of basis matrix: 

.

By definition, M_B = M_H · M_HB.   Therefore,

.

Our derivation of the Bezier form parametric polynomial is now complete.  The function is as follows:

Q_B(t) = T · M_B · G_B.

Splines and Continuity

Four control points is often not enough to define your curve with enough precision.  The "natural cubic spline" accomodates n control points, but calculating the polynomial coefficients requires solving for a n + 1 by n + 1 matrix.  This can become extremely demanding on the computer.  Additionally, moving a single control point may change the shape of the natural cubic spline wildly.  In response to these problems, graphics programmers usually join together multiple cubic curves to approximate a curve fit to an arbitrary number of control points .  To achieve continuity, the curve segments must share control points.  Below are two Bezier curves joined together through control point P₄:

This looks swell.  Look how smooth the curve is.  Now lets move control point P₅ and see what happens:

Curses!  Two curves joined by a single point are not smooth.   Let's define Cⁿ continuity as the equivalence of the n^th derivatives of the left and right curves at their joint.  Two Bezier curves that share P₄ achieve C⁰ continuity - that is, their endpoints have equal positions.  They do not, however, achieve C¹ continuity - the tangent vectors at their endpoints are not necessarily equal.  Is there a way to join two Bezier curves together to achieve C¹ continuity?  Yes.

We define a function L(t) that refers to the left Bezier curve of our two-curve system, and similarly define a function R(t) that refers to the right Bezier curve.  C⁰ continuity is achieved when L(1) = R(0).  In the example above, this is true because L(1) = P₄ and R(0) = P₄.  C¹ continuity relates to the velocity of the curve, and is achieved when L'(1) = R'(0).  Recall the Bezier-Hermite geometry equivalencies from the previous section:

R₁ = -3P₁ + 3P₂

R₄ = -3P₃ + 3P₄

R₁ and R₄ are the velocities of the curves evaluated at t = 0 and t = 1.  Therefore, C¹ continuity of two joined Bezier curves is achieved when 3(P₄ - P₃) = 3(P₅ - P₄), or P₄ - P₃ = P₅ - P₄.  The following two systems are joined Bezier curves with C⁰ and C¹ continuity:

This scheme allows the joining of Bezier curves with first degree continuity.  There is no evident relationship between the Bezier control points and the second derivative of the curve.  C², the second degree of continuity, is achieved when the acceleration of the curves are equal at the joints.   We introduce a new form of the parametric cubic curve that achieves C² continuity: the B-spline. 

The B-spline is defined by a set of control points P₀, P₁, ..., P_m-1, P_m, for m ³ 3 with m - 2 curve segments.   The curve segments are named Q₃(t), Q₄(t), ..., Q_m-1(t), Q_m(t).  We define each curve segment Q_i(t), for 3 £ i £ m in its own time domain (that is, each segment is defined for 0 £ t £ 1), and adjust the spline domain t_i to 0£ t £ 1, for 3 £ i £ m, when evaluating segments.   B-splines come in many flavours: nonrational uniform, nonrational nonunfirom, and rational nonuniform.  In a uniform spline, every curve segment is defined over a domain 0 £ t £ 1.  The nonuniform spline is a generalization of the uniform spline - this restriction has been lifted.  Rational splines are defined in homogenous space (i.e. evaluation of the function requires homogenizing each term by some function W).  The B-spline covered here is uniform and nonrational.

Like the two cubic curves discussed in the previous sections, segments of the nonrational uniform B-spline are defined by four control points.  The function for each segment is defined in the now-familiar way:  Q_i(t) = T_i · M_s · G_si, for 3 £ i £ m, where M_s is the B-spline basis matrix and G_si is the geometry matrix for the i^th segment.  Each geometry matrix is composed of four consecutive control points:

, where n is the dimension and 3 £ i £ m.

The following two systems are B-splines constructed from four cubic curves.  Where each curve segment joins, we can draw a point called a "knot."  The knots are named the spline domain at that point, that is, t_i for 3 £ i £ m.

How does the B-spline form of the parametric cubic curve achieve C² continuity?  The answer, of course, is in the B-spline basis matrix M_s.   To derive M_s we need to look back at the Bezier form.

Recall that the Bezier function, Q_B(t), is defined as the product of the time matrix, Bezier basis, and the Bezier geometry matrix:

Q_B(t) = T · M_B · G_B.

Expanding this product we find that:

Q_B(t) = (-t³ + 3t² - 3t + 1)P₁ + (3t³ - 6t² + 3t)P₂ + (-3t³ + 3t²)P₃ + t³P₄

We can use factoring to clean up this function:

Q_B(t) = (-t³ + 3t² - 3t + 1)P₁ + 3t(t² - 2t + 1)P₂ + 3t²(-t + 1)P₃ + t³P₄

Upon seeing this we say, "Holy Cow."  The coefficient on each control point is a row from Pascal's triangle.  This means we can factor each coefficient into a binomial expression:

Q_B(t) = (1 - t)³P₁ + 3t(1 - t)²P₂ + 3t²(1 - t)P₃ + t³P₄

This shows very clearly how each control point affects the curve at time t.  We call each coefficient a "blending function" of the Bezier form:

B₁(t) = (1 - t)³

B₂(t) = 3t(1 - t)²

B₃(t) = 3t²(1 - t)

B₄(t) = t³.

We can think of each blending function as a weight on its respective control point.  By examining the graph of the blending functions we can understand how each control point affects the Bezier curve at some time t:

The first impression is symmetry.  At t = .25, the B₁ blend function has a value equal to the B₄ blend function at time t = .75.  The same is true for B₂ and B₃.   We can also graph the sum of the blend functions:

This graphs claims that S B_i(t) = 1, for 0£ t £ 1.  Can it really be? 

B₁(t) + B₂(t) + B₃(t) + B₄(t)

= (1 - t)³ + 3t(1 - t)² + 3t²(1 - t) + t³

= (-t³ + 3t² - 3t + 1) + 3t(t² - 2t + 1) + 3t²(-t + 1) + t³

= (-t³ + 3t² - 3t + 1) + (3t³ - 6t² + 3t) + (-3t³ + 3t²) + t³

= 1.

Indeed, the sum of the Bezier blend functions is one for all values of t.  This means that Q_B(t) is a weighted average of the control points with respect to t.  This is a significant property, and it will help in the derivation of the B-spline basis matrix.

To avoid subscript proliferation, we redefine B₁(t), B₂(t), B₃(t), and B₄(t) to be the blending functions of the B-spline.  The values of these functions at time t are the weights of the control points on the curve at t.  Like in the Bezier form, we want the sum of the blending functions to equal one at any time.  This means that each B-spline segment function is the weighted average of its control points at t, for 0£ t £ 1:

We construct the B-spline basis matrix such that segment Q_i(t) evaluated at t = 1 will have the same position, derivative, and second derivative as Q_i+1(0) evaluated at t = 0.   That is, we need assurance that C⁰, C¹, and C² continuity is achieved. 

Two consecutive B-spline segments, Q_i(t) and Q_i+1(t) are defined by five control points: P_k for i - 3 £ k £ i + 1.  Using this powerful notation, we define:

Q_i(t) = B₁(t) · P_i-3 + B₂(t) · P_i-2 + B₃(t) · P_i-1 + B₄(t) · P_i

Q_i+1(t) = B₁(t) · P_i-2 + B₂(t) · P_i-1 + B₃(t) · P_i + B₄(t) · P_i+1.

Additionally, to assure continuity:

Q_i(1) = Q_i+1(0)

Q_i'(1) = Q_i+1'(0)

Q_i''(1) = Q_i+1''(0).

To setup the derivation, we simply apply our new definition of the curve segment into the equalities from above:

B₁(1) · P_i-3 + B₂(1) · P_i-2 + B₃(1) · P_i-1 + B₄(1) · P_i = B₁(0) · P_i-2 + B₂(0) · P_i-1 + B₃(0) · P_i + B₄(0) · P_i+1

B₁'(1) · P_i-3 + B₂'(1) · P_i-2 + B₃'(1) · P_i-1 + B₄'(1) · P_i = B₁'(0) · P_i-2 + B₂'(0) · P_i-1 + B₃'(0) · P_i + B₄'(0) · P_i+1

B₁''(1) · P_i-3 + B₂''(1) · P_i-2 + B₃''(1) · P_i-1 + B₄''(1) · P_i = B₁''(0) · P_i-2 + B₂''(0) · P_i-1 + B₃''(0) · P_i + B₄''(0) · P_i+1.

This is not a system you'd want to meet in a dark alley.  Recall that the blending functions are cubic polynomials:

B_i(t) = a_it³ + b_it² + c_it + d_i.

In our discussion of the Bezier blending functions, it was shown that the blending function coefficients are actually the elements from the basis matrix.  The same is true here.  By solving for the 16 blending function coefficients (four functions each with four coefficients) we derive the B-spline basis matrix.

We can solve the three continuity equations above by setting each B_k(1) · P_i+k-1 = B_k-1(0) · P_i+k-1, for 1 £ k £ 4:

We can reduce the control point term from both sides of all equations.  This brings us a step closer to the solution:

Because each blending function is a cubic polynomial, we can rewrite the fifteen equations from above in terms of the blending function coefficients.  For example, B₁(1) = 0 is equivalent to:

(a₁t³ + b₁t² + c₁t + d₁)_t=1 = 0

a₁(1)³ + b₁(1)² + c₁(1) + d₁ = 0

a₁ + b₁ + c₁ + d₁ = 0.

The ten first and second differential equations can be rewritten in a similar way.  For example, B₁'(1) = 0 is equivalent to:

B₁'(t) = 3a₁t² + 2b₁t + c₁

(3a₁t² + 2b₁t + c₁)_t=1 = 0

3a₁(1)² + 2b₁(1) + c₁ = 0

3a₁ + 2b₁ + c₁ = 0

The previous fifteen equations are equivalent to:

We are trying to solve for a matrix with 16 elements, so we'll need 16 equations.  Recall that the sum of the blending functions at any time is one.  Each blending function can be evaluated at any time, but we'll choose to evaluate it at t = 0 because the coefficients on twelve variables in the equation become zero, so those variables drop out entirely:

B₁(0) + B₂(0) + B₃(0) + B₄(0) = 1

d₁ + d₂ + d₃ + d₄ = 1

Finally, it's time to solve for the B-spline basis matrix.  If you are sufficiently crazy, you can load the coefficients on each of the 16 equations into a 16x16 matrix and Gaussian eliminate to solve.  I'll just use Mathematica:

Oh, happy day!  We can factor out a 1/6 to make the basis matrix easier to write:

.

In conclusion, a B-spline segment is defined as:

.

You now have an understanding of the Hermite, Bezier, and B-spline parameterizations of parametric cubic curves.  These cubic curves have   applications that extend beyond the realm of graphics.  Because they are continuous and piecewise differentiatable (the B-spline even has a defined second derivative at segment joints), cubic curves can be used to parameterize object transformations.  Unlike Euler angles, object rotation fit to a cubic curve is unambigious.  Additionally, no expensive trigonometry operations are required for evaluation.  The bicubic surface, an elegant application of parametric cubic curves, can replace the discrete and storage-intensive triangle mesh in many scenarios.  The uses of parametric cubic curves are numerous, and I feel deeply honored to have helped create a spline-friendly world.

References:

James Foley, Andries van Dam, et al., Computer Graphics : Principles and Practice, (Reading, Massachusetts; Addison-Wesley, 1996).

Alan Watt and Mark Watt, Advanced Animation and Rendering Techniques, (New York: Addison-Wesley, 1992).