Discrete Comput Geom 1:123-140 (1986) )isa,rete & ('qmqmtatiiMud

Geometrv
© 1986S~;er-Verla8 NewYorkInc. o'

S m o o t h , E a s y to C o m p u t e Interpolating Splines*

John D. Hobby
Computer Science Department, Stanford University, Stanford, CA 94305

Abstract. We present a system of interpolating splines with first-order and

approximate second-order geometric continuity. The curves are easily com-
puted in linear time by solving a diagonally dominant, tridiagonal system of
linear equations. Emphasis is placed on the need to find aesthetically pleasing
curves in a wide range of circumstances; favorable results are obtained even
when the knots are very unequally spaced or widely separated. The curves are
invariant under translation, rotation, and scaling, and the effects of a local
change fall off exponentially as one moves away from the disturbed knot.
Approximate second-order continuity is achieved by using a linear "mock
curvature" function in place of the actual endpoint curvature for each spline
segment and choosing tangent directions at knots so as to equalize these. This
avoids extraneous solutions and other forms of undesirable behavior without
seriously compromising the quality of the results.
The actual spline segments can come from any family of curves whose
endpoint curvatures can be suitably approximated, but we propose a specific
family of parametric cubics. There is freedom to allow tangent directions and
"tension" parameters to be specified at knots, and special "curl" parameters
may be given for additional control near the endpoints of open curves.

1. Introduction

The problem of fitting a smooth curve through a set of points on the plane has
many important applications in computer graphics, computer-aided design, and
typesetting. Often there is no preexisting curve to approximate except possibly a
freehand drawing, and the only requirement is to find an aesthetically pleasing

*This research was supported in part by the National Science Foundation under grants
IST-820-1926 and MCS-83-00984 and by the Systems Development Foundation.
124 J.D. Hobby

curve that the computer can easily manipulate. For some interactive applications
the curves can be controlled by manipulating points that do not lie on the curve,
but many applications require the control points to lie on the curve. For example,
the control points may be obtained by digitizing key points on a drawing, or there
m a y be a priori knowledge that the curve must pass through certain points.
Suppose the curve must pass through points z0, z 1..... z,; and either z 0 and
z n are to be the endpoints of the curve or z 0 = zn and the curve is to be a closed
loop. Optionally, there may be direction vectors wi specifying the slope of the
curve at some z i. For example, some of the zi may have been selected as vertical
extrema so that the curve must pass through them horizontally. It is desirable for
the curve to be invariant under translation, rotation, and scaling in the sense that
if T is such a transformation, then applying T to the computed curve should yield
the same result as computing a new curve through T z i for 0 < i < n with direction
vectors T w i - T(0,0).
The curve should have at least approximate continuity of slope and curvature
at knots where no directions are given, and it would also be desirable to have
some notion of e x t e n s i b i l i t y a n d locality. A system of splines is extensible if the
curve generated from knots z 0, z 1..... z~ is identical to that generated from knots
z~, z~..... z'n+l where z,' = zi for i < k, z~ = zi_ t for i > k, and z~ is on the curve
segment joining z k _ ~ a n d z k. In other words, adding a new knot already on the
curve must not change it. In practice it is extremely difficult to achieve exact
extensibility. The only well-known extensible spline family is the "curve of least
energy" that minimizes the integral of squared curvature with respect to arc
length [5, 9], but this curve is difficult to work with. It is interesting to note that
when the knots are nearly collinear, the curve of least energy approaches the
simple nonparametric cubic spline passing through the given knots with continu-
ous second derivative. The splines that we deal with here will share this property.
The concept of locality is that if one of the knots or direction vectors is
perturbed, the changes should be confined to a few surrounding spline segments.
Specifically, there should be some constant k such that changes to z~ or w~ effect
only the curve segments between z~_ k and zi+ k. As the example of Fig. 1 shows,
it is difficult to have both locality and continuity of curvature even when the
knots are collinear. If the knots zo, z 1, z 2, and z 3 are as shown in the figure but
wo is in the direction of zl - z o, then the desired curve is obviously a straight line.
if wo is changed, then locality demands that the spline segments between z k and
n must remain straight, yet there is no way a cubic curve can join such a straight
line with continuous curvature.
Most local interpolating splines of moderate complexity do not have continu-
ous curvature. Perhaps the best known splines of this type are the cardinal splines
given by a cubic function z ( t ) , where the velocity vector z ' ( i ) at knot zi is
½(zi÷ 1 - z i _ l ) . Figure 2 shows an example of cardinal splines in comparison to
the splines that will be developed in this paper. The price paid for locality is that
Fig. 2(b) has wild discontinuities in curvature while Fig. 2(a) does not.

Fig. 1. The effect of changing wo while preserving

zo zl 2~ z~3 exactlocality.
Smooth, Easy to ComputeInterpolatingSplines 125

(a) (b) (c)

Fig. 2. (a) Cubics with mockcurvatureconstraints.(b) Cardinalsplines. (c) Barsky-DeRosesplines
with B, = IIz,+a- z,ll/llz, - z,-db

Other well-known local cubic interpolating splines can be viewed as generali-

zations of cardinal splines. To see this, it is convenient to parameterize the cubic
spline in terms of the incoming and outgoing velocity vectors D 7 and D~+ at each
knot so that the spline segment between z~ and z~+l has B~zier control points z,,
z i + ½D/+, zi+ 1 13Di+l,
-
-
- and z i + p For cardinal splines as described above,
D Z = D + = ½(zi+ 1 - Z~_l); in the most general form

D Z = a , ( z , - z,_1) + bi(z,+ 1 - z,) and

D? = ¢ , ( z , - z,_l) + d , ( z , + l - z,), (1)

where a;, b i, c i , and di are parameters that may be used to control the shape of
the curve. For instance, Kochanek and Bartels [7] give a scheme where the user
can give continuity, tension, and bias parameters for each knot, and these
determine ai, bi, cj, and d~. The default settings of these parameters result in
a i = b i = c i = d i = ½, giving simple cardinal splines as shown in Fig. 2(b).
Another approach to local interpolating splines can be found in [3], where
DeRose and Barsky give an infinite family of spline curves based on two integer
parameters: One parameter gives the desired order of continuity, and the other
determines whether the splines interpolate the knots or just pass near them.
DeRose and Barsky give explicit equations for cubic interpolating splines of this
form with first-order continuity, but their approach requires polynomials of
degree at least 5 in order to obtain interpolating splines with continuous curva-
ture. In the cubic interpolating splines derived explicitly in [3], each knot zi has a
single shape parameter fl~ associated with it. The splines are exactly those
126 J. D. Hobby

(a~ (b) (c)

Fig. 3. (a) Natural cubics, (b) Cubics with chordal parameterization. (c) Cubics with mock curvature
constraints.

produced by (1) where

t3, 1 B~ and d, = 1
a i = fli+l, bi = f l i ( f l i + l ) , ci = f l i + l , fli+---~"

The default suggested in [3] is to set all fl, = 1 so as to obtain cardinal splines as in
Fig. 2(b), but somewhat better results can be obtained by setting fli equal to the
ratio of Euclidean lengths Ilzi+l - zA[/llzi- z~-ltl as shown in Fig. 2(c).
In this article we shall see how to obtain smother curves by sacrificing locality
and settling for an exponential decline in influence of local changes. The best
known interpolating splines that behave this way are the C 2 continuous paramet-
ric cubic splines, sometimes called "natural cubic splines." If no directions are
given, there is a unique piecewise parametric cubic, closed curve that is C 2
continuous with respect to the parameter and passes through n given points in
order. This curve can easily be computed by solving linear equations.
As Epstein shows in [4], natural cubic splines do not perform well for
unequally spaced knots because the spacing of parameter values at knots does not
reflect the spacing of the knots. Better results can be obtained by setting the
parameter at each knot z i to a value ti, where tj - tj_ 1 = Ilzj - zj-ll[ for 1 _< j < n,
and requiring second-order continuity with respect to the parameter as shown in
Fig. 3(b). This chordal parameterization improves on the uniform parameteriza-
tion of Fig. 3(a), but the splines that we shall develop still have more gentle
curvature in this case as shown in Fig. 3(c).
Figure 3 points out the difference between geometric and parametric continu-
ity. Natural cubic splines have second-order parametric continuity because they
are described by a function z ( t ) that is C 2 continuous. Geometric continuity is a
Smooth, Easy to ComputeInterpolatingSplines 127

relaxation of parametric continuity: A spline function z(t) has k th-order geomet-

ric continuity if there is a continuous, monotonic increasing function f such that
z ( f ( t ) ) is C k continuous. In other words, geometric continuity is independent of
the parameterization.
Requiring first- and second-order continuity with respect to the parameter
uses up four degrees of freedom per knot, enough to completely determine a
parametric cubic spline as shown in Fig. 3(a). In Fig. 3(b), one of these degrees of
freedom is reclaimed and put to better use by altering the parameter spacing, but
another degree of freedom can be made available by requiring only continuity of
slope and curvature.
Other results on geometric continuity in cubic splines can be found in [1] and
[2] where Barsky and Beatty show how two extra degrees of freedom can be
obtained for B-splines by requiring only geometric continuity. We need to obtain
similar degrees of freedom for interpolating splines, but we shall first concentrate
on finding good defaults to work from. The new parameters will be of little value
if it is difficult to set them so as to obtain reasonable results. We cannot afford to
assume arbitrarily that the best defaults are those that yield some form of
parametric C 2 continuity.
In [8], Manning takes an interesting approach to this problem. He defines a
specific family of curves so that there is a unique one for each pair of initial and
final points and directions. He then selects spline directions at each knot so as to
achieve geometric continuity. His approach provides a certain degree of locality in
that effects of local perturbations do not propagate past knots where a direction is
given.
With Manning's approach, both degrees of freedom are available to control
the shape of the curve, and defaults can be selected so as to obtain the most
pleasing curves. Section 2 explains how to select the defaults by choosing two
functions and using them to determine the magnitudes of the velocities at each
knot in such a way as to guarantee that the curves generated will be independent
of scaling, rotation, and reflection. We can then provide two "tension" parame-
ters for each knot by simply dividing them into these functions. Essentially the
same approach would work for other kinds of curves, although there may be more
parameters to choose. We select parametric cubics here because they are essen-
tially the simplest curves that can pass through two arbitrary points in two
arbitrary directions.
On apparent disadvantage to this approach is the difficulty in solving for the
directions that provide continuity of curvature. Manning proposes an iterative
approximation scheme that seems to work well in practice, but he admits that
there is not always a unique solution and there is no guarantee that the iteration
always converges to the desired solution. Cubic splines often have very low
curvature at their endpoints when they have very sharp bends internally, and this
can introduce extraneous solutions as shown in Fig. 4. The three curves shown are
all curvature continuous open curves that have given directions at z o and z 2, but
regardless of the initial conditions, Manning's iteration always converges to one
of the asymmetrical ones with sharp bends. If z 0 is raised and z 2 lowered until
the angle ZoZ~Z2 is about 122 °, the asymmetrical solutions merge with the
symmetrical ones and the rate of convergence for Manning's iteration approaches
zero.
128 J. D. Hobby

Fig. 4. Three splines of the type proposedin [8].

While these kinds of problems do not seem to occur when the angles involved
are not so large, much additional testing would be necessary in order to verify
this. In Section 3, we show how all such problems can be avoided by setting up a
system of linear equations that are easy to solve and guarantee approximate
continuity of curvature. We derive the specific equations appropriate for the
families of curves discussed in Section 2, but similar equations could be derived
for many different classes of curves.

2. Choosing a Family of Parametric Cubics

The subproblem to be solved in this section can be stated as follows: Given two
points z~ = (x~, y~) and Z~+l -- (x~+p Y~+x), and two unit vectors wj and w~+x, find
an aesthetically pleasing parametric cubic z(t) so that z ( 0 ) = zi, z ( 1 ) = z~+ 1,
z'(0) = aw i, and z'(1) = flwi+ x, where a and fl are positive real numbers and z'(t)
is the componentwise derivative (x'(t), y'(t)) of z(t) = (x(t), y(t)). We also wish
to introduce two "tension" parameters ~-~and ~+1 such that pleasing curves will
be obtained when ~'~= ~÷1 =1, and as the tensions approach o0, the curves will
approach the line segment joining z~ to z~+ 1.
In order to guarantee that the results are independent of translation, rotation,
and scaling, we shall begin by finding a function 2(t) such that

~(0) = (0,0), ~(1) = (1,0), ~'(0) -- l # ( o , ¢ ) ( c o s O , s i n O ) , and

2'(1) = _1 o ( O , , ~ ) ( c o s e p , - s i n , ) , (2)
"ri+1

where 0 and o are positive real functions to be determined later, 0 = arg wi -

arg(zi+ 1 - zi), and ~ = arg(zi+ t - z i ) - a r g w i + t . We refer to p and o as velocity
parameter functions because they determine the magnitude of the velocity at t = 0
and at t = 1. (Here arg(x, y ) is the angle ,o such that (x, y ) is a positive multiple
of (cos,0,sin~o).) We then set

z(t)= .,,, +[y,+_y, (3)

It is not hard to see that the parametric cubic satisfying (2) has Btzier control
points (0,0), (p/3~,)(cos0,sin0), ( 1 - ( o / 3 ~ + t ) c o s ~ , ( o / 3 ~ + O s i n ~ ) , and (1,0),
Smooth, Easy to Compute Interpolating Splines 129

so that

~(t) = Pit(1-t)2(cosO,sinO)+ t2(1-t)(3 - ° s i n q ~ o s i n q ) ) + t 3 ( 1 , 0 ) . (4)

• ~ ~/+1' ~+1

It only remains to choose positive functions O(0, qS) and o(0, q0 so that O(0, qs) =
o(q,, 0 ) = p ( - 0, - ,t,).
In [8] Manning chooses
2
o(O,q,) = a+(l-c)cosO+ccose: and
2
o(0, q,) = 1 + ccos qb+ (1 - c)cos 0 (5)

and then empirically selects c = ] to obtain the most pleasing family of curves.
Here we shall attempt to do a systematic analysis of the vast range of possible
functions to determine whether slightly more complicated velocity parameter
functions will yield better results. These functions will have to be evaluated only
once for each segment of the spline curve, and they have a strong influence on the
final shape.

2.1. Mathematical Measures of Smoothness

One common way of evaluating the smoothness of curves is to integrate the
square of curvature with respect to arc length. For 2 = (~, .~)

f k~'~ = fo (Y"~' ~"~")~ tit.

(::+ - ~,,~)'/~ (6)

This can be simplified somewhat but it still proves to be analytically intractable.

This is not surprising considering the complex behavior of numerical solutions.
Equation (6) is exactly the energy function that the curve of least energy
minimizes, but if we restrict 2 to be the cubic spline (4), we can investigate the
velocity parameter functions that minimize (6). Actually we should consider the
smallest local minimum since (6) approaches 0 as p and o approach oo for fixed 0
and q~.
Unfortunately, numerical integration of k 2 ds proves to be slow and impre-
cise, and it would have to be repeated a large number of times in order to get a
good idea what the velocity parameter functions p(0, q~) and o(0, 4) should be.
Instead we shall introduce two other measures of smoothness that behave
similarly:

max Jk(t)l (7)

ONt<l
and

O<t_<l
130 J.D. Hobby

2.5" 2.5-
pl 9,%)
2.o- ¢) ,~,¢) 2.0-

1.5- 1.5-

1.0 = ...... 1.0-

- p(O ~) --"--- /
.5-

.0-
0° i0 o 20 ° 30 ° 40 ° 50 ° 0o 10 ° 20 ° 300 40 ° 50 o

(a) (b)

Fig. 5. The functions p(O, ~b) and 0(0, ~) versus 0 for ff = 90 ° - 0, minimizing (a) the magnitude of
curvature and (b) the magnitude of curvature change.

Using (7) to measure smoothness corresponds to taking the o0-norm of curvature

instead of the 2-norm; using (8) gives roughly similar results but applies a greater
penalty to short periods of relatively high curvature. The velocity parameter
functions that minimize (8) turn out to be somewhat better behaved than those
that minimize (7), but the overall character is the same for both measures of
smoothness.
For fixed 0 and ~, the smoothness measures can have multiple local minima
and the relative smoothness at the local minima can change as 0 and ,~ change.
Therefore, it should not be surprising that the "optimum" velocity parameter
functions have large discontinuities where they catastrophically change from one
local minimum to another. When this happens there tend to be (p, o) values
between the two minima that also generate relatively "smooth" curves, so it is not
really necessary to use discontinuous functions for p and o.
Figure 5 illustrates the most basic catastrophe. Near 0 = ~ = 45 °, the "opti-
m u m " p increases and the o decreases as 0 decreases. This action tends to reduce
the curvature where it is maximum near t = 1 without introducing other points of
high curvature. When 0 ~-0 << ~, the situation is entirely different. The high
curvature near t = 1 is best controlled by making o large, and increasing p
beyond what is needed to control the curvature at t = 0 just makes the problem
worse.
When p and o are chosen to minimize (7) as shown in Fig. 5(a), there are
actually two catastrophes, and the short segment between them is particularly
interesting. Ordinarily, extremely small o values lead to high curvature near t = 1,
but at 0---34.1 °, ~ = 55.9 °, the curvature is actually minimized by choosing
o = 0.261. The choice of p = 2.423 here is extremely critical. As shown in Fig. 6,
this has the effect of making the last three Bdzier control points almost collinear,
so that the endpoint curvature is not too large in spite of the low velocity. (The
bold dots in the figure are the Bdzier control points of (4).)

Fig. 6. The BCzier control polygon for (4) with 0 34.1 ° , ~ 55.9 o, where O a n d
o are chosen to minimize curvature.
S m o o t h , E a s y to C o m p u t e I n t e r p o l a t i n g Splines 131

3.0
. . . . • tQ~/~ 3.0-

2.5- / 2.5

2.0--- / -~ 2.0- - - - -
/~r.t. (7i J
1.5 J 1.5 J
,.r.t.__..A~~
1.0 1.0 ,.l.t. 3)

.5 .5
0o 20 ° 40 ° 60 ° 80 ° 100 ° 0o 20 ° 40 ° 60 ° 80 ° 100 °

(a} (b}

Fig. 7. (a) " ' O p t i m u m " O a n d o versus 0 w h e n (a) 0 = ~ a n d (b) 0 = - ~ , for b o t h s m o o t h n e s s

measures.

The bizarre situation shown in Fig. 6 does not occur when minimizing the
derivative of curvature, but there is still a catastrophe near 25 °. When 0 + ~ = 90 °
and 0 < 25 o the " o p t i m u m " cubic has a point of inflection, but when 0 > 25 ° it
has none.
Figure 7 shows how the "optimal" velocity parameter functions grow as 0
and q~ increase. Minimizing (8) produces a catastrophe at 68 ° where p and o
increase to about 1.58 and the cubic acquires a single point of maximum
curvature at t = 0.5. For # = o ~ 1.34, the maximum curvature occurs at t ~ 0.08
and t --- 0.92. Intermediate values for p and a produce relatively high change in
curvature near t = 0 and t = 1.
Minimizing (7) instead of (8) avoids the catastrophe at O = q~= 68 °, but then
p and a do not approach unique limits as (0, ~) ---, (0,0). Along 0 = - ~ the limit
is ~ - / 2 , while p and a approach 1 as 0 ~ 0 when O = q~. Under the approxima-
tion k = d2y/dx 2, when either (6) or (8) is used as the measure of smoothness, it
can be shown that the optimum curves are cubics where p cos0 = o cosq~ = 1.
Thus, it seems reasonable to let (p, a) ~ (1,1) as (8, ~,) ~ (0,0).

2.2. Velocity Parameter Functions

We are now ready to choose velocity parameter functions p(0, +) and o(0, +) so
as to obtain an aesthetically pleasing family of curves from (4). In order to be
useful as splines these curves should be fairly easy to compute, and they should
respond predictably to changes in 0 and +. The " o p t i m u m " velocity parameter
functions illustrated in Figs. 5 and 7 fail on both counts, but they still provide
useful guidelines. The catastrophes in the " o p t i m u m " velocity parameter func-
tions are not features to be preserved, but rather points where many different
possible values for p and a are feasible. For instance when 0 = 28 ° and + = 62 °,
good results are obtained for almost any values of p and o as long as p + o ~ 2.5.
The actual choice of velocity parameters is necessarily somewhat arbitrary,
and there is a trade-off between "smoothness" and simplicity. Other properties
such as approximate extensibility and predictable response to changes are also
132 J. D, Hobby

important, but empirical studies indicate that these goals also tend to be served
by maximizing "smoothness" and smoothing out the catastrophes.
We have already decided that p ( 0 , 0 ) = o ( 0 , 0 ) = 1. Thus for small angles we
can approximate the behavior of the curve of least energy and achieve very good
approximate extensibility. For 0 = ~, we should approximate the behavior of the
functions shown in Fig. 7(a), but these functions increase too rapidly or large
angles: They seem to approach ~ well before 0 and g, reach 180 °. It is
convenient to let

2
p = o = l+cos0 for0 = ~ (9)

so as to obtain good approximations to circles.

Because of the symmetry requirements p(0, ~) = o(g,, 0) = p ( - 0, - ~), it
suffices to choose p and o for 0 < 10l _< ~ < 180 °. Figure 8(a) shows the p(0, ,p)
and o(0, ~) that minimize (8) for g, = 80 °. Figure 8(b) shows practical functions p
and o that smooth out the catastrophes and are consistent with (9). Similar plots
for smaller g, would have p and o closer to each other and closer to 1. (The slope
discontinuities at 64.7 ° and 69.6 ° are due to changes in the relative sizes of
extrema in d k / d s at different parts of the cubic curve.)
If complexity is of no concern, we might want to choose p(0, g,) and o(0, g,)
as follows for 0 < 101 --<'/' < ~r with angles measured in radians:

p = f(O, ep) + ~,(g})sin ff~, ,

(o)
a a 2 + ol + 2 c
f(O,ep)-- a+ccosB+¢'

a=(,-O) 2 to+ '

1.17
-y(~) = ---~#~ - 0.15sin(2~),

+o.,- (10)

A least-squares fit of f to ½(p + o) with p and o chosen to minimize (8) yielded

a = 0.2678306, c = 0.2638750, d =1.402539, and e = 0.7539063. A possible re-
finement is to require p ~ 1.5 sing,/sin(0 + ~) when g, > ½z, and - g , < 0 < 0 so as
to avoid any possibility of generating a curve with a cusp in it. (This only affects
the above functions when ~ > 145°.)
It is desirable to have a simpler approximation that does not use transcenden-
tal functions other than sines and cosines of 8 and ~. One such approximation is
S m o o t h , E a s y to C o m p u t e Interpolating Splines 133

2.5-

22- 2.0

1.5" ~ -gr"~ _,

,..~ ~;o,¢ //
1.s

2X
"E-~\\

°¢) / ._..
1.0- "- -'~ "-~ x 1,0- ~ j , - --
/ T(0,~)
.5- .5 ¸

.0- .0 ]
_80 ° _40 ° 0o 40 ° 80 ° -80 o -40 ° 0o 40 ° 80 °

(a) (b)

Ftg. 8. (a) "Optimum" and (b) practical functions p(O,~) a n d e ( 0 , ~ ) v e r s u s 0 f o r q, = 80 ° .

the following functions which were developed for the new METAFONTsystem [6]:
2+a
P = 1 + ( 1 - c)cos0 + ccosq~'
2-a
1 + (1 - c)cos ff + ccos0 '
where

a = a (sin 0 - b sin q~)(sin q~- b sin 0 ) (cos 0 - cos q~). (11)

The constants a, b, and c were chosen to minimize an error function based on the
value of (8) for curves generated by (11) for 116 different (0, q~) pairs. This
suggested a =1.597, b = 0.0700, and c = 0.370, but since empirical evidence
indicated that large values for JO - o] were causing problems, METAFONTuses the
slightly perturbed values a = v~-, b = ~ , and c = ½(3- v~).
Figures 9(a) and (b) show how velocity parameter functions compare in terms
of the maximum magnitude of change in curvature. The vertical axis in each
graph gives the ratio of the value of (8) for the indicated velocity parameter
functions to that for the "optimum" velocity parameter functions. In other words
for each value of 0 and q,, the graphs give the ratio of (8) to its minimum possible
value. Since (10) and (11) were derived in order to minimize (8), it is not
surprising that they perform better than (5) in Figs. 9(a) and (b). We would also
expect the more complicated functions defined in (10) to perform better than the
simpler versions of (11), and this is indeed the case except near 0 - - - 2 0 °,
q~= 80 °, where riO, qJ) is a little too small.
Figures 10(a) and (b) correspond to Figs. 9(a) and (b), except that the ratios
are based on the maximum magnitude of curvature instead of the maximum
magnitude of change in curvature. The graphs show that the curves generated by
(10) and (11) do have less curvature than Manning's curves, especially when q~ is
!airly large and 0 < 0. The figures also show that (7) and (8) do behave similarly
m that they both produce the same overall ranking: Curves from (10) are
smoother than those from (11) which are smoother than Manning's curves.
134 J.D. Hobby

3.5.

18 3.0 It
1.0- - - . (5)

~.o ~7" ~" ~ 1.o ~ /

-40 ° --20 ° 0o 20 ° 40 ° -80 ° -40 ° 0o 40 ° 80 °
(a) (b)
F i g . 9. (a), ( b ) Smoothness ratios based on (8), comparing optimum velocity parameter functions to
those of (5), (10), a n d (11).

1.25" 1°

1.2 - - - ? --
11)~ 1.5

1.1s-// / / 1.4

1.1" /
lo//~ \
1.05-

1.0-
/ \
1.0
-40 ° -20* 0° 20 ° 40 ° -80 ° -40 ° 0o 40 a 80 °
(a) (b)

F i g . 10, (a), ( b ) Smoothness ratios based on (7), comparing optimum velocity parameter functions to
those of (5), (10), a n d (11).

Figure 11 shows some of the curves generated by (10) and (11). They are
similar for moderate angles, but the simpler equations set p too small and o too
large when q~ = - 9 0 °. Equation (11) does not perform well in such extreme cases
because it does not allow p - o to be large enough when 0 << 0/q, < 1 without
making o - p too large when - 1 < O/q, < 0 or moving the crossover point where
p = o too close to O/q~ = O.

(a) (b)

Fig. II. (a) Curves from (10) and (b) curves from (11), both with 0 = 45%
Smooth, Easy to Compute Interpolating Splines 135

3. M o c k Curvature Constraints

Here we need to extend the notation of Section 2 so that 0j = arg wj - arg(zj+ 1 - zj)
for 0 _< j < n, q~y= arg(zj - Z j _ l ) - a r g w j for 0 < j < n, and dj is the Euclidean
length of the vector z j + ~ _ zj for 0 < j < n. If the problem is to find a closed
curve with no directions gwen, it will be convenient to sometimes use alternative
names z,+ 1, On, ~,+1, ~',, and ?,+1 for z 1, 00, ~1, 'r0' and 71, respectively. We can
then define q9 = arg(zJ+l - z j ) - a r g ( z j - zj_l) for 1 < j < n', where n ' = n for
closed curve problems with no directions given and n' = n - 1 otherwise. Unless
stated otherwise, all q9 are at most 180 ° in absolute value.
Where w i has been given in advance, it simply determines q~i and 0i; other w,
need to be determined by solving for O, and q~,. Since the problem of finding
direction angles can be broken into independent subproblems separated by knots
where directions are given, we can assume that no directions other than w0 and w,
are given. F o r closed curve problems we can assume that no directions at all are
given, otherwise the problem could be reduced to one or more open curve
problems.
The requirement that the curvature be continuous at some knot z t, 0 < i < n,
is

kl(Zt_l, Zt,Wt_l,Wt,q'i_l,?i) ~---ko(zi, zi+l,Wt,Wt+l,Tl,~+l),

where k 0 and k 1 are functions that give the curvature at t = 0 and t = 1 in terms
of the endpoints, terminal directions, and tension parameters for the family of
curves being used. Because of the requirement for invariance under translation,
rotation, and scaling, there exists a function k such that

,o(Zj, =

and

kl(Zj, Zj+l,Wj,Wj+l,~,~j+l) = k(d~j+l, Oj,Tj+l, Tj)/d j. (12)

A n y particular family of curves determines a specific function k that satisfies

(12). The corresponding mock curvature function Ic consists of the linear terms in
the Taylor series for k(O, ep, r, ?), expanded about ( 0 , ~ ) = (0,0). For the curves
determined by (3) and (4) with p and o determined by (10) or (11),

k(O, ¢, r, ?) = 20(0, ¢)sin(0 + ¢ ) / ? - 6 s i n 0

(0(0,*)/*) 2

and

~:(O,~,~',?) = r2(2(0;q~)-60), (13)

where the angles are measured in radians. Since the tension parameters are always
known in advance, they are treated like constants in this expansion.
136 J.D. Hobby

Instead of requiring continuity of curvature, we will obtain a system of linear

equations by requiring the mock curvatures to match at each knot. When the
angles 0~ and ¢i are small, this guarantees that the resulting splines will have
approximate second-order continuity. For instance

l k ( 0 , q , , 1 , 1 ) - k(O,q,,1,1) I
max(Ik(O,,l,,1,1)l,lk(q,,O,l,1)l) < 0.11

when 1Ol, Iq'l --< 22.5° and to and o are determined by (11) with a = v~-, b = ~,
and c = ½(3- v~-).
In many applications the angles 0~ and q~ seldom get much larger than 22.5 °
but experience has shown that the splines usually appear very smooth even when
some 0r and ~ are much larger than 22.5 °. An iterative algorithm such as
Manning's can be used to obtain true curvature continuity if desired, but this
necessitates giving up the simplicity and guaranteed behavior of a system of linear
equations.
Continuity of mock curvature requires

k( Oi, ~i+l,'Ti, ~i+l) = 0 forl_< i < n'. (14)

di-1 di
Combining this with the first-order continuity equations

0r+~i= -~k, f o r l < i < n' (15)

gives enough equations to determine all 0 i and q~i for closed curve problems. For
open curve problems when directions w0 and w. are given in advance, these
provide the necessary additional equations by fixing 0o and ep., but otherwise
additional constraints are needed.
The additional constraints are controlled by special "curl" parameters X0 and
X.. There should be one such constraint for each endpoint where no direction is
specified. They have the form

?,(00, ¢,, = Xo?,(¢,, 0o, (16a)

and

= (16b)

The curl parameters give the ratio of the mock curvature at the endpoints to that
at the adjacent knots. They should probably have default values of 1 so that the
first and last spline segments will usually be good approximations to circular arcs
as in Fig. 12(b).
We now have a system of equations consisting of (14), (15), and possibly
(16a) a n d / o r (16b). If 00 or ~n have been given in advance then they may be
Smooth, Easy to Compute Interpolating Splines 137

2o 222 20 22 20 2,2

(a) (b) (c)

Fig. 12. (a) Xo=x2=O.(b) Xo=X2=l.(c) Xo=X2=~e-

regarded as constants. The first step is to rewrite (16a) and (16b) as

0o = ~'o3+Xoq~(3~o - 1 ) q~+X.~'3_l(3q,,-1)
%3(3~1_l)+Xo~3q, 1 and q,. = .7,3(3%_l_l)+x.r~_lO._1
(17)

so that 0o and Cn can be eliminated. Then (15) may be used to eliminate all ¢~ so
that the remaining variables are 0x, 02..... On,, and the remaining equations are
those given by (14) with appropriate substitutions. This system has some im-
portant properties that may be summarized as follows.

Theorem 1. I f n > 2, if all tension parameters satisfy the bound ~',, ~ >_ rmin> 2,
and if any curl parameters satisfy X o, X . > O, then after the aforementioned substitu-
tions, all coefficients of 01, 0 z ..... 0., in (14) are nonnegative, and for each i, the
coefficient of 0 i is at least 3Tmi. - 1 times the sum of all the other coefficients in that
equation.

Proof The bounds on the tension parameters guarantee that the coefficient of 0
in (13) will be negative, the coefficient of ~ will be positive, and the magnitude of
the former will be at least 3~-min - 1 times the latter. When 1 < i < n' in (14), the
only relevant substitutions are ~j = - ~kj - q~j for I J - i[ < 1, so the coefficients of
0,_ 1, 0r, and 0r+ ~ clearly have the required properties. For closed curve problems,
the same holds for i = 1 and i = n', otherwise additional substitutions eliminate 00
and q~. so that ~:(qh,00, ql,~0) depends only on q~l and k(0._l,q~.,%_l,~.)
depends only on 0._ 1. We need only show that both of these variables have
nonpositive coefficients. This is clearly true for given directions, and it also holds
when curl parameters apply since the coefficients in (17) are at most 3~"0 - 1 and
3q. - 1, respectively. []
Theorem 1 shows that subject to certain reasonable limitations on the tension
and curl parameters, the system of equations is diagonally dominant, and hence it
has a unique solution. Actually the solution is unique only up to the choice of the
angles q~i- Ordinarily all q'i should be chosen so that they are at most 180 ° in
absolute value, but it is possible to add multiples of 360 ° to them. The effect of
such a change is usually to add a loop to the curve as in Fig. 13.
Theorem 1 also shows that the splines have approximate locality in the sense
that changes in direction angles fall off exponentially as one moves away from a
disturbance. Specifically, suppose a given direction 0o is displaced by an angle 8
and let A be the matrix of coefficients of 01, 02.... 0., from (14) after the
substitutions. The change in 01, 02..... 0., due to this displacement is given by the
solution vector O to AO = 8e x where e I = (1,0,0 ..... 0) r.
138 J, D, Hobby

Z2

(a) (b)
Fig. 13. A spline computed (a) with qJ2 = - 9 0 ° and (b) with +z = 270%

Fig. 14. Exponential decline in the effect

ZI) Zl 7.2 Z3 z4 of a 45 ° change in direction.

We know that A is tridiagonal with nonnegative entries, and within each row
the diagonal element dominates the sum of the other two elements by at least a
factor of 3 ' r m i n - - 1. It is not hard to see that for any two adjacent components of
O, either O k = 0 or O k _ l / O k _<1 - 3¢min. This is trivial for k = n', and it may be
extended inductively to smaller k using the fact that Akk >_ (3%i n --1)(Ak,k_ 1 +
Ak,k+l). Thus, j knots away from where a given direction is changed; the effect
of the change is reduced by at least a factor of ( 3 % i n - 1 ) L In practice the
reduction is often by a somewhat greater amount as in Fig. 14 where %in = 1 and
0o/01 = ,1
11"
When a knot z~ is displaced, three mock curvature constraints are directly
affected due to changes in d,_ a, d i, ff,-x, ~P~,and qJg+l. The adjustment will cause
some change in q),-i and 0,+ 1, and the effect on 0i+j and 0~ j for j > l is
equivalent to what would happen if directions w,+ 1 and w,_ 1 were given in
advance. The change in 0t+j will be at most 1 / ( 3 % i . - 1 ) j-1 as great as the
change in 0i+ 1, and the change in 0i_j will be at most 1 / ( 3 ~ - ~ . - 1) j-~ as great as
the change in q~_ r If the original problem were to find a closed curve with no
directions given, then these two effects will add together so that the change in 0i+j
will be at most 1/(3%i n - 1 ) j-1 of the change in 0,+ 1 plus 1/(3%i . - 1) " - t - j of
the change in q,i_ ~.

4. Conclusion

We have developed a tridiagonal system of linear equations that can be solved in

linear time to determine the spline direction at each knot so as to match mock
curvatures. It is necessary to use arctangents to set up the system of equations,
and to use sines and cosines to recover the resulting spline directions; but this
work can be reduced to one arctangent, one sine, and one cosine per knot on the
spline. When all unit direction vectors w~ have been determined, the ith spline
segment is the cubic whose B6zier control points are

p(o,,,,+O
zi, zi + ~ I I z i + l - zAIw. z~+l 3~+1 IIz~+l- z~llwi+l,

and Zi+l°
Smooth, Easy to Compute Interpolating Splines 139

0 '3~'~'~4 0

(a) (b) (c)

(d) (e} (t3

Fig. 15. (a) "r,=?,=0.75. (b) "r,=~=0.85. (c) "r,='~=l. (d) ~',=?,=1.2. (e) r , = ~ = l . 5 . (f)
~', = ~, = 2.

Thus, the total work required in order to find a computationally convenient

representation of the desired spline is linear in the number of knots. The constant
factor obviously depends on the implementation, but it is usually not much higher
than the similar factor for natural cubic splines.
The splines do not have locality, but we have shown that changes in direction
angles fall off exponentially. The rate of decline depends on how small the
tension parameters are allowed to be, but at least a factor of 2 per knot is
guaranteed for the default tensions, and decay rates as high as a factor of 4 per
knot are typical. It should be noted that an exponential decline in angular change
does not guarantee that curve displacements decline similarly because it is
technically feasible for dj to be exponential in j.
The curve families discussed in Section 2 and defined by (10) and (11) are
somewhat arbitrary, and the concept of mock curvature could be applied to other
families of curves. It would be desirable to find p and o functions simpler than
(10) that perform better than (11), although even the simplified functions of (11)
produce very good results for problems such as that shown in Figs. 2(a) and 3(c).
In fact the splines produced by (11) have performed outstandingly in extensive
tests with Knuth's new METAFONTsystem [6]. It is easy to compute the splines by
solving the equations given in Section 3, but readers wishing a detailed descrip-
tion of one implementation can refer to [6].
The tension parameters have proved very useful in METAFONT,because it is
often much easier to control a shape by specifying tensions than by specifying
directions or giving more knots. Figure 15 shows the effect of adjusting all the
tension parameters simultaneously, and Fig. 16 shows the effect of adjusting only
one tension parameter and leaving the rest of the tension parameters at the
default values (~'i= ~ = 1). Note how the tension effects the direction chosen at
knot 2 in Fig. 16.
140 J.D. Hobby

(a) (b) (c)

Fig. 16. (a) ;r2 =1.2. (b) ~2=1.5.(c) 72=2.

Acknowledgments
The author is indebted to Professor Donald Knuth of Stanford University for doing the first complete
implementation and helping to refine some of the ideas. Professor Kevin Karplus of Cornell
University helped with the numerical investigation of the "optimal" p and a functions.

References
t. Brian A. Barsky, The Beta Sptine: A Local Representation Based on Shape Parameters and
Fundamental Geometric Measures, Ph.D. thesis, Univ. of Utah, December, 1981.
2. Brian A. Barsky, John C. Beatty, Varying the betas in beta-splines, Univ. of California, Berkeley
TR CSD 82/112, December 1982.
3. Tony D. DeRose and Brian A. Barsky, Geometric continuity and shape parameters for Catmull-
Rom Splines, Graphics Interface '84, 57-64.
4. M.P. Epstein, Parametric interpolation, SIAM Journal of Numerical Analysis 13 (1976), 261-268.
5. B . K . P . Horn, The curve of least energy, Massachusetts Institute of Technology A.I. Memo No.
612, January 1981.
6. Donald E. Knuth, Computers and Typesetting, Vot. 4: METAFONT the Program, Addison
Wesley, to appear.
7. Doris H. U. Kochanek and Richard H. Bartels, Interpolating splines with local tension, continu-
ity, and bias control, Computer Graphics 18 (1984), 33-41.
8. J.R. Manning, Continuity conditions for spline curves, Computer Joumat t7 (1974), 181-186.
9. Even Mehlum, Curve and Surface Fitting Based on Variational Criteriae for Smoothness, Oslo,
1969.

Received Janua~ 31, 1985, and in revised form September 4, 1985.