Bob is a programmer/analyst for Computer Sciences Corp. and an independent animator. He can be contacted at rkauffma@ slb.isd.csc.com.
Splines let you create smooth curves simply by specifying a few points. The precise shape of the curve is determined by a set of "basis functions" which are selected to provide certain properties. For example, some selections of basis functions guarantee that the spline will actually go through each control point.
The most-common splines use four curves, each bounded more or less within the unit square over the interval [0,1] on the x- and y-axes. As Figure 1 shows, the resulting curve is just the sum of each basis function multiplied by the corresponding control point. By placing the basis functions side by side, you can see how they regulate the smoothness of the spline.
These splines are called "uniform," because a single set of basis functions is used for all curve segments. Nonuniform splines use different and overlapping sets of basis functions for every curve segment. They are more powerful than uniform splines because they offer a greater degree of control over the curve, and remain the method of choice for high-end applications. Uniform splines are still in wide use because they are intuitively easier to understand and implement.
Cubic Splines
Splines are traditionally based on cubic functions because they are simple and easy to control, yet they provide the inflection point needed for changing the curve's direction and curvature. The most popular uniform splines are the BSpline and Catmull-Rom spline. The BSpline is an exterpolating cubic spline, which means that the resulting curve does not pass through its control points, but stays inside the bounding box generated by its control points. The Catmull-Rom spline is an interpolating cubic spline that passes through every control point, but is not as smooth as the BSpline.
Uniform cubic splines are expressed as a set of parametric equations for each segment of the spline. The most compact way to write these equations is as a single matrix equation. For a cubic spline, a segment can be written as r(t)=T
MP. Here, P is a matrix holding the coordinates of the four control points, M holds the coefficients of the four cubic basis functions, and T expresses powers of t; see Figure 2. Generating one curve segment just requires stepping t from 0 to 1, and computing the matrix equation at each point. The result is a succession of X, Y pairs to be plotted.
Uniform cubic splines are easy to compute, but do have some drawbacks. The main trade-off between the Catmull-Rom spline and the BSpline involves continuity. Wherever two curve segments join, BSplines are continuous in the second derivative (curvature). Catmull-Rom splines, on the other hand, are only continuous in the first derivative (slope). However, BSplines do not pass through the control points. Nonuniform splines (such as NURBS) solve this problem, but at the expense of intuitiveness and ease of implementation.
Trigonometric Splines
In this article, I'll present a class of uniform spline curves whose basis functions are trigonometric. These curves take the best features of their popular cubic spline counterparts, but are rarely addressed in the popular literature about graphics. In fact, though I understand work has been done in this area, I did not discover these new splines in literature. I discovered them through trial and error in search of a better spline curve. I'll refer to these trigonometric spline curves as "TSplines."
I examine two varieties of TSplines here -- an interpolating TSpline, which passes through all of its control points, and an exterpolating TSpline, which does not. TSplines can be expressed using a matrix equation similar to the one used for cubic splines. The difference is that the product TM is replaced by a single matrix F, whose entries are functions of t; see Figure 3.
TSplines have two important advantages over cubic splines. Both varieties of TSpline are continuous through the second derivative and beyond. TSplines can also be used to draw exact circles and other conic sections. TSplines share all of the useful and desirable properties of their cubic counterparts; see Table 1. Figure 4 shows exterpolating and interpolating TSplines imposed over their cubic counterparts.
Designing the TSpline
Figure 5 shows the trigonometric basis curves imposed over the cubic spline basis curves. For a set of basis curves to be useful for producing splines, they have to satisfy certain conditions:
- All of the basis functions, when summed together, will equal 1. This ensures that the resulting spline will lie within a bounding box defined by the minimum and maximum coordinates of its control points in each dimension.
- For an exterpolating spline, the basis functions must all be between 0 and 1. This ensures that the spline will be within the bounding box of the control points.
- For interpolating splines, the functions each evaluate to zero where they join, except for where F2 and F3 join, where they evaluate to 1.
- Separate spline segments should have smooth joins. Mathematically, the derivatives should be continuous where the curves join.
- Splines should be invariant under affine transformation. That is, you can scale, rotate, or translate a spline simply by transforming the control points. You get this for free if you write your splines as a matrix product.
In the case of TSplines, continuity in every derivative is achieved by finding a set of basis equations that are pieces of the same sine curve. Each piece is successively phase shifted in 90 degree increments so as to be constrained on the interval x = [0,1]. Continuity between F1 and F4 is achieved by setting the period of the functions that define F1...F4 to 4.
To fulfill the third condition for interpolating TSplines, the basis was constructed as the product of two sine functions with zero crossings at every junction except where F2 and F3 join. Figure 6 shows the basis functions for the interpolating and exterpolating TSplines overlaying one another, accompanied by their equations.
Rendering Conics with TSplines
Interpolating TSplines can draw exact circles and ellipses; just use the four corners of a rhombus as control points. If the rhombus is a square, the result will be a circle. You can also use rectangles, but rhombuses make it easy to plot ellipses by specifying a center point, semi-major and semi-minor axes, and an angle of orientation. From these parameters, the control points (which form the corners of a rhombus) can be generated for an interpolating TSpline to generate all possible ellipses and circles.
Conics may also be generated using the exterpolating TSplines, but are not as useful. The ellipse generated is about 1/2 the dimensions of the bounding box formed by its control points. Figure 7 shows a circle and an ellipse generated using exterpolating TSplines.
Implementation
Listing One is a Pascal implementation of both TSplines and cubic splines. A more complete set of listings is available electronically (see "Availability," page 3), including an interactive demonstration program that lets you build and edit splines with the mouse.
Since TSplines require many sine and cosine calculations, speed is a major concern. Listing One addresses this with the SetBasis procedure, which precomputes the chosen basis equation and stores the results in an array. The PlotCurve procedure can then use the values from the array to render each curve segment.
I tested this approach on a 486DX66, plotting 500 points per curve segment. Without precomputing the basis functions, the cubic splines took 130 milliseconds to render on average, while the TSplines required about 220 milliseconds. When I changed the code to precalculate the basis functions, only 45 milliseconds was required to render each curve segment for both cubic splines and TSplines.
References
Courtney, Mike. "A Cubic Spline Extrema Algorithm," DDJ, April 1996.
Johnson, Stephen P. and Tom McReynolds. "Implementing Curves in C++," DDJ, December 1992.
Lauzzana, Raymond G. and Denise E.M. Penrose, "Implementing Bicubic Splines," DDJ, August 1990.
Rabbitz, Rich. Course Material for "Interactive Graphics," Lockheed Martin, Fall 1993, Moorestown, NJ.
DDJ
Listing One
{--------------------------TYPES------------------------------------}type POINT = record { define a 2-D point}
X, Y : integer;
end;
</p>
PT4 = array[1..5] of POINT; {Array of points controlling 1 seg }
</p>
A4_4 = array[1..4,1..4] of real; {Basis matrices for cubics}
A1000_4 = array[1..1000,1..4] of real; {Lookup table for basis
curve points; up to 1000 points
are pre-calculated in advance
and stored for each of the
four basis equations }
Curve_Type = ( Cat_Rom, {Cubic spline which passes through every
control point continuity in first derivative}
B_Spline, {Cubic spline which does not pass through
control points; continuity in 2nd derivative}
Ext_Tspline, {Triginometric spline which does not pass
through control points; continuity in
all derivatives}
Int_Tspline ); {Cubic spline which passes through all
control points; continuity in all derivatives}
</p>
{--------------------------VARIABLES--------------------------------}
var
Max_Points : real; {Max number of points to interpolate per segment}
Curve : A1000_4; {Current interpolating curve}
Curr_Basis : Curve_Type; {Type of basis function being used}
</p>
{-------------------CUBIC CURVE BASIS MATRICES----------------------}
const
Bspline : A4_4 =
( ( -0.166667, 0.5, -0.5, 0.166667 ),
( 0.5, -1.0, 0.0, 0.666667 ),
( -0.5, 0.5, 0.5, 0.166667 ),
( 0.166667, 0.0, 0.0, 0.0 ) );
Catrom : A4_4 =
( ( -0.5, 1.0, -0.5, 0.0 ),
( 1.5, -2.5, 0.0, 1.0 ),
( -1.5, 2.0, 0.5, 0.0 ),
( 0.5, -0.5, 0.0, 0.0 ) );
{--------------------------------------------------------------------}
function CubicSpline( t : real;
n : integer;
M : A4_4 ) : real;
{ Calculates value of nth cubic basis function specified by M for parameter
t; calculates basis curve one point at a time; called by SetBasis }
var T1, T2, T3 : real;
begin
T1 := t; T2 := t*t;T3 := t*t*t;
CubicSpline := T3*M[n,1] + T2*M[n,2] + T1*M[n,3] + M[n,4];
end; {CubicSpline}
</p>
{-------------------------------------------------------------------}
function ExTSpline( t : real;
n : integer ) : real;
{ Calculates the value of the n-th basis function for the
exterpolating triginometric spline for parameter t; used to
calculate the basis curve one point at a time; Called by
SetBasis }
var Sn, Cs : real;
begin
Sn := sin( 0.5*Pi*t ); Cs := cos( 0.5*Pi*t );
case n of
1 : ExTSpline := 0.25*( 1 - Cs);
2 : ExTspline := 0.25*( 1 + Sn);
3 : ExTspline := 0.25*( 1 + Cs);
4 : ExTspline := 0.25*( 1 - Sn);
end;
end; {ExTspline}
{-------------------------------------------------------------------}
function InTspline( t : real;
n : integer ) : real;
{ Calculates the value of the n-th basis function for the
interpolating triginometric spline for parameter t; used to
calculate the basis curve one point at a time; Called by
SetBasis }
var Sn, Cs : real;
begin
Sn := sin( 0.5*Pi*t ); Cs := cos( 0.5*Pi*t );
case n of
1 : InTspline := 0.5*Cs*(Cs - 1 );
2 : InTspline := 0.5*Sn*(Sn + 1 );
3 : InTspline := 0.5*Cs*(Cs + 1 );
4 : InTspline := 0.5*Sn*(Sn - 1 );
end;
end; {InTspline}
{-------------------------------------------------------------------}
procedure SetBasis( Basis : Curve_Type;
Max : integer);
{Pre-calculate basis function at specified granularity}
var t : real;
i, j : integer;
begin
Curr_Basis := Basis;
if Max > 1000 then Max_Points := 1000 {set granularity <= 1000}
else Max_Points := Max;
for i:= 1 to Max do begin
t:= i/Max;
case Basis of {calculate basis curves on point a time}
Ext_Tspline : {exterpolating trig spline}
for j:= 1 to 4 do Curve[i,j] := ExTspline( t, j );
Int_Tspline : {interpolating trig spline}
for j:= 1 to 4 do Curve[i,j] := InTspline( t, j )
Cat_Rom : {Catmull-Rom cubic spline}
for j:= 1 to 4 do Curve[i,j] := CubicSpline( t, j, CatRom );
B_Spline : {BSpline cubic spline}
for j:= 1 to 4 do Curve[i,j] := CubicSpline( t, j, BSpline );
end;
end;
end; {SetBasis}
{-------------------------------------------------------------------}
procedure PlotCurve ( P1,P2,P3,P4 : POINT;
C : integer );
{Plots a segment of the defined curve in the specified color;
called by Render_Curve }
var Xo, X1, X2, X3, X4, t, T1, T2, T3,
Yo, Y1, Y2, Y3, Y4, X, Y : real;
i, j, Max : integer;
begin
Max := round( Max_Points );
X1 := P1.X; X2 := P2.X; X3 := P3.X; X4 := P4.X;
Y1 := P1.Y; Y2 := P2.Y; Y3 := P3.Y; Y4 := P4.Y;
for i:= 1 to Max do begin { vary thru [0,1] & plot basis functn}
t:= i/Max;
X := X1*Curve[i,1] + X2*Curve[i,2] +
X3*Curve[i,3] + X4*Curve[i,4] ;
Y := Y1*Curve[i,1] + Y2*Curve[i,2] +
Y3*Curve[i,3] + Y4*Curve[i,4] ;
put_pixel( X, Y, C );
end;
end; {PlotCurve}
</p>
Copyright © 1997, Dr. Dobb's Journal
