- Email: [email protected]
Local control of interval tension using weighted splines
Computer Aided Geometric Design 3 (1986) 281-294 North-Holland
281
Local control of interval tension using weighted splines Thomas A. FOLEY * Computer Science Dept., Arizona State University, Tempe, A Z 85287, U.S.A. Received 4 February 1986 Revised 9 March 1987 Abstract. Cubic spline interpolation and B-spline sums are useful and powerful tools in computer aided design. These are extended by weighted cubic splines which have tension controls that allow the user to tighten or loosen the curve on intervals between interpolation points. The weighted spline is a C 1 piecewise cubic that minimizes a variational problem similar to one that a C 2 cubic spline minimizes. A B-spline like basis is constructed for weighted splines where each basis function is nonnegative and nonzero only on four intervals. The basis functions sum up identically to one, thus curves generated by summing control points multiplied by the basis functions have the convex hull property. Different weights are built into the basis functions so that the control point curves are piecewise cubics with local control of interval tension. If all weights axe equal, then the weighted spline is the C 2 cubic spline a n d the basis functions are B-splines. Keywords. Interpolation, approximation, curves, computer aided design, B-splines.
1. Introduction
The theory and application of C 2 cubic spline interpolation to computer aided design are well documented in [Schumaker '81, Boehm et al. '84, de Boor '78, Faux et al. '79]. Piecewise cubic tension splines in [Nielson '74, Kochanek '84] provide controls for making a curve tighter or looser at the interpolation points. The weighted splines in [Salkauskas '84] and the exponential-based splines under tension in [Barsky '84] provide tension controls that make the curve tighter or looser on intervals between interpolation points, rather than at the discrete interpolation points. The weighted splines are piecewise cubics, while the splines under tension are piecewise exponential functions. Weighted splines are discussed in Section 2 and useful end conditions are developed that parallel those for C 2 cubic spline interpolants. B-splines, described in [Boehm et al. '84, de Boor '78, Faux et al. '79, Cohen et al. '80, Gordon et al. '74], are a useful and powerful tool in computer aided design. They form a basis for the space of C 2 cubic splines, and each B-spline is a nonnegative cubic spline that is nonzero only on four intervals. The B-splines form a partition of unity, that is, they sum up to one. Curves generated by summing control points multiplied by the B-splines have some very desirable shape properties, including the local convex hull property. The Beta-splines in [Barsky '81] and [Barsky et al. '83] are similar to B-splines in that they are nonnegative piecewise cubics, they are nonzero only on four intervals and they sum to one. Two parameters are incorporated into the Beta-splines that control the tension at a point and the bias of a curve generated by summing control points multiplied by the Beta-splines. If * A portion of this work was supported by the U.S. Department of Energy at Lawrence Livermore National Laboratory under contract No. W-7405-ENG-48 and No. DE-FG02-87 ER 25041. 0167-8396/86/$3.50 © 1986, Elsevier Science Publishers B.V. (North-Holland)
T.A. Foley / Local control of interval tension
282
different tension values are used at various points, then an interpolation between parameters is performed that yields a higher degree piecewise polynomial. In Section 3, a B-spline like basis is constructed for weighted splines where each basis function is nonnegative and nonzero on only four intervals. This basis is normalized so that it forms a partition of unity. Section 4 uses this basis to form a control point sum will yield a curve that has the convex hull property. The different weights are built into the basis functions so that the resulting control point curve is a piecewise cubic with local control of interval tension.
2. Weighted spline interpolation This section describes weighted spline interpolation with different end conditions. The weighted spline developed in [Salkauskas '84] did not involve end conditions, although his approach is equivalent to using Type II or natural end conditions as described below. Given a = t 1 < • • • < t n = b, Yl,-.-, Yn, and a positive integrable weight function w(t), the approach used in [Salkauskas '84] finds the C a interpolant to the given data that minimizes the quantity
fabW(t)[f"(t)] 2 dt. When w(t) is a piecewise constant, the solution is a piecewise cubic function and the minimization problem reduces to solving a tridiagonal system of equations whose unknowns are the first derivatives at the data points. The minimizing function is the piecewise cubic Hermite interpolant to the data using the solution of the tridiagonal system for the derivatives. Some motivation for the weighted cubic spline is that the C 2 natural cubic spline is the unique interpolating function that minimizes lab[ f " ( t ) ] 2 dt subject to certain continuity conditions. Minimizing this global measure of curvature is effective if the data set is well behaved because it places equal weight on all intervals. When using the weighted spline, if w(t) is large on one interval relative to the other intervals, then f " ( t ) is forced to be small in magnitude on that interval, hence f ( t ) will be more linear on that interval. Similarly, smaller relative weights allow f " ( t ) to take on larger values. The term relative is used because if all weights are multiplied by a positive constant, then the resulting weighted spline would be the same. As noted in [Salkauskas '84], if the weight function w(t) is constant, then the weighted spline is identical to the C 2 natural cubic spline. Given a = t 1 < t 2 < • • • < t n = b, Yl,..., Yn, and w i > 0 for i = 1 . . . . , n - 1, define w(t) = w i if t i ~< t < ti+ 1. Using a piecewise constant weight function, the weighted spline interpolant is the C 1 piecewise cubic function f ( t ) that minimizes
V= f bw(t)[f"(t)]2 dt
(1)
subject to f(ti) =Yi for i = 1 , . . . , n, and one of the following end conditions: (I) f ' ( a ) = r n l (II) f " ( a + ) = 0 (III) f ' ( a ) = f ' ( b )
and
f'(b)=m~,
and
f"(b-)=0,
and
or
Wxf"(a+)=w~_lf"(b-).
Type I are called first derivative end conditions, Type II are called natural end conditions, and Type III are called periodic end conditions. First derivative and natural end conditions are
T.A. Foley / Local control of interval tension
283
identical with those used by a standard C 2 cubic spline. The rationale for periodic end conditions is discussed later in this section. A constructive approach is used to verify that a unique weighted spline exists, subject to the gix/en end conditions. This approach will yield a method for computing the weighted spline using cardinal piecewise cubic Hermite basis functions. In [de Boor '78], it is shown that given Yi and m i, for i = 1 , . . . , n, then there exists a unique C 1 piecewise cubic that satisfies f ( t i ) =y~ and f ' ( t i ) = m i for i = 1 , . . . , n. Furthermore, f ( t ) can be represented by n
f(t)=
(2)
~,yiPi(t)+miDi(t), i=1
where Pi(t) and Di(t) are the cardinal piecewise cubic Hermite functions that satisfy Pi(tj)=Si,j,
Pi'(tj)=O,
Di(tj)=O,
D;(tj)=~ij,
for i, j = 1 . . . . , n. Since the functions Pi(t) and D~(t) are zero outside of [t~_ 1, t;+l], equation (2) simplifies computationally to f ( t ) = y i P i ( t ) + m i D i ( t ) + yi+lPi+l(t) + m i + a D i + l ( t )
for t i <<.t <~t~+ 1. The unknowns for the weighted spline are the derivative values m~. Necessary conditions for minimiT.ing V in (1) are OV Om k
----~0
for k = 2 , . . . , n - 1. It follows that m k _ l W k _ l h k l _ l + m k 2 ( W k _ l h k l _ l + Wkh-~ 1) + mk+lWkh-~ 1 = 3Wk_l(Yk--Yk_a)h-k21 + 3W/,(yk+ 1 --yl,)h-k 2,
(3)
for k = 2 , . . . , n - 1, where h i = ti+ 1 - t r The n - 2 equations in (3) involve the n unknowns m l , . . . , m,,. The remaining two equations are generated from the given end conditions. For Type I or first derivative end conditions, the other equations are m 1 = value of f ' ( a )
and
m , = value of f ' (b).
The two additional equations for Type II or natural end conditions are m 1 2 h l 1 + m E h l 1 = 3 ( y 2-y1)h1-2,
and --1
-1
--2
m n _ l h n _ 1 + m . 2 h . _ 1 = 3(y. - y n _ l ) h ~ _ l .
For Type III or periodic end conditions, the equations are m l ( 2 W a h f I + 2 W n _ l h n- _1 l ) + m 2 w l h f I + m . . l W . _ l h - n l _ l
= 3(y 2
--
-2
ya)hl 2+ 3(y,,-y~_l)hn_l,
and m 1 = m no
The linear systems of equations that occur when using (3) and either Type I or Type II end conditions are tridiagonal and diagonally dominant. N o t only do they have unique solutions,
284
T.A. Foley / Local control of interval tension
v
w
Fig. 1. This is a parametric weighted spline interpolant with natural end conditions using weights w i = 10, 1, 20, 1, 10.
but they can be solved efficiently using a standard tridiagonal system solver, such as the one in [Conte et al. '80]. If the two equations for Type III end conditions are placed in the first and last rows of the linear system, then with a slight modification of the tridiagonal system solver, the unique solution can be efficiently computed. The previous linear system of equations are also sufficient conditions for minimizing V. An argument similar to those in [Salkauskas '84] and [Foley '87a] can be used to verify this. The weighted splines are C ~ functions in general, but as shown in [Salkauskas '84], if wi_ 1 = wi, then the weighted spline is C 2 at t~. This result follows from the fact that wif
w
+
(t i ) = W i _ l f
.,v
(t;).
(4)
This behavior of the second derivative is the rationale for the Type III end conditions that are termed periodic. In fact, if all the weights are the same, then the weighted spline with specified end conditions is identical to the C 2 standard cubic spline with the same end conditions. For parametric interpolation, assume that we are given (x;, yg) for i = 1 , . . . , n. Let t 1 < --< tn, let X ( t ) be the weighted spline interpolant to the data (tg, x~) and let Y(t) be the weighted spline interpolant to the data ( ti, Y i ) , i = 1 . . . . . n. The parametric curve ( X ( t ), Y( t )), t 1 ~ t ~< t, passes through the given data. For 3-D space curves, the z-component Z ( t ) is computed similarly. Fig. 1 is the parametric weighted spline using natural end conditions and weights w~ = 10, 1, 20, 1, 10 using a uniform parametrization with t i = i. This example shows how easy it is for the curve to be nearly linear on one interval and have high curvature on bordering intervals. Some visual guidelines on the effects of the weights are given in Fig. 2 using the uniform parametrization. All of the curves interpolate the six points (0, - 1 ) , (0, 0), (0, 1), (1, 1), (1, 0) and (1, - 1) using w i = 1 for i ~ 3. F r o m top to bottom on the middle interval, the weighted spline interpolant uses w 3 = 1, 5, 20 and 100, respectively. The v-splines in [Nielson '74] and the Beta-splines in [Barsky et al. '83] are C 1 as functions, but they have continuous curvature as parametric curves. The parametric weighted splines are not G 2, that is, the curvature is not continuous at the interpolation points. If they were G 2 continuous, then Y"(t~-)-
Y"(t~-)=vY'(ti)
for some value of v, as shown in [Farin '86]. Consider the example where t i = 1, 2, 3, 4, 5, y; = 0, 1, 14, 4, 0, and w t = 11, 1, 2, 1. The natural weighted spline then has derivative values of m i = 0, 3, 0, - 1 2 , 0. Since Y ' ( 3 ) = 0 and Y ' ( 3 ÷) :~ Y " ( 3 - ) , the weighted spline is not G 2 in general. As noted earlier, if wi -- wi_l, then each component of the curve is C 2, and thus if the tangent vector is nonzero, the parametric curve has continuous curvature at t i. Although the weighted spline is not generally G 2, it is generally 'smoother' than many C 1 curves in that it has a continuous osculating plane everywhere. A discussion of the osculating plane at a point on a
T.A. Foley / Local control of interval tension
285
• 2 ~ ~'.~._.-:...-_ ..........................
:.3
_-:._-_.__-:. -..
S
',,'i1
.6
:'I
l
J i~ .2
1 .2
8.
Iq
i
.2
.
I
I
I
.6
.8
1.g
..
Fig. 2. T h e s e are w e i g h t e d s p l i n e i n t e r p o l a n t s u s i n g w i = 1 for i ~ 3, a n d f r o m t o p t o b o t t o m o n t h e m i d d l e i n t e r v a l , t h e v a l u e s o f w 3 are 1, 5, 20 a n d 100, respectively.
3-D space curve is given in [Mortenson '85, p. 271] where an e q u a t i o n for the p l a n e at
(X(t), Y(t), Z(t)) is given by [ x - X ( t ) ] [ Y ' ( t ) Z " ( t ) - Y"(t)Z'(t)] -[y- Y(t)][X'(t)Z"(t)-X"(t)Z'(t)] + [z- Z(/)][X'(t)Y"(t)X"(t)Y'(t)]
=0.
U s i n g (4), it follows that the osculating plane at t = t~ is continuous.
3. A B-spline like basis for weighted splines W e present a basis No(t),..., Nn+l(t) such that if f(t) is any weighted spline, that is, a C 1 piecewise cubic that minimizes (1), then there exists u n i q u e constants c; such that
f(t)
is
n+l
f(t)=
E
ciNi(t)
(5)
i=0
for all t I ~< t ~< t,. Each basis function N~(t) is a weighted spline, N ~ ( t ) > / 0 for all t, a n d N~(t) = 0 for all t outside of [ti_ 2, ti+2]. T h e basis functions will be n o r m a l i z e d or f o r m a partition of unity, that is, n+l
N/(t) = 1
(6)
i=0
for all t I ~< t ~< t,. F o r notational convenience, we will first c o n s t r u c t the basis using a u n i f o r m p a r a m e t r i z a t i o n where t, = i a n d then we state the results for an arbitrary parametrization. W e first construct basis functions Bg(t) and later define the n o r m a l i z e d functions N,.(t) in terms of them. G i v e n weights Wa,... , Wn_l, define the weights w 0 = w_ 1 = w_z = Wl a n d w , = w , + 1 =
286
T.A. Foley / Local control of interval tenswn
w~+ 2 = wn_ 1. F o r Bi(t ) to be a C 1 piecewise cubic a n d have c o m p a c t s u p p o r t on [i - 2, i + 2], if follows that
0 = n i ( i - 2) = B/' ( i - 2) = B i ( i + 2) = n/' (i W 2). Since B~'(i - 2 - ) = 0, b y (4) it follows t h a t B f ' ( i - 2) = 0. Similarly, B~'(i + 2) = 0. Since B~(t) is a C 1 piecewise cubic, it can be represented in terms of its B6zier points for j <<.t <~j + l by 3
3 ) u (1 _ u)3-k,
Bi(t) = E bj,k( k k=O
where u = t - j . Based o n the zero properties of B~(t) discussed in the previous paragraph, we have that bi_2, k = bi+l, k = 0 for k = 0, 1, 2. F o r B~(t) to be C ~ at t = i - 1, i, a n d i + 1, we m u s t have that bi-2,3 = bi-l,O ~ ( b i - 2 , 2 q- b / - 1 , 1 ) / 2 b i _ l , 3 -~ bi, 0 = ( b i _ l , 2 q- b i , 1 ) / 2 ,
-~- b i - l , 1 / 2 ,
and
hi,3 = bi /~alo = ( bi,2 + bi+1,1)/2 = bi,2/2.
T h u s the only u n k n o w n B6zier points are bi-l,1, 'b'i2a,2, bi,1 a n d bi, 2. If we arbitrarily set bi-1,~ = 1, then we can c o m p u t e the three u n k n o w n B6zier p o i n t s by requiring that (4) holds at the points i - 1, i a n d i + 1. I n terms of the B6zier points, these three equations are
6 w j ( b j , o - 2bja + bj,2) = 6 w j - a ( b j - l a -
2bj_x, 2 + bj_l,3),
for j = i - 1, i, i + 1. In terms of the three u n k n o w n B6zier points, these three equations simplify to
2wi_lbi_l, 2 = 3wi_ 1 + wi_ 2, (w, + 3wi_l)bi_l, 2 + ( - 3 w i - wi_l)bi, 1 + 2wibi,2 = 2 w i _ l , 2wibi, 1 + ( - 3 w i - wi+l)bi, 2 = O. T h e u n i q u e solution of the above system is hi_l,2 = ( 3 w i - 1 q- w i _ 2 ) / ( 2 w i _ l )
,
bi,1 = (wi+ 1 + 3 w i ) d i _ l / ( 2 d i ) ,
and
bi, 2 = w i d i _ l / (
W i_ldi) ,
where
(7)
d i = 3 w i w i + 1 + W i _ l W i + 1 q- 5w2i --t- 3 W i _ l W i.
Observe that Bi(t ) >1 0 because all of its B6zier points are nonnegative. Before we show that B o ( t ) , . . . , B , + l ( t ) is the basis for w e i g h t e d splines, we first determine a i > 0 and define normalized functions Ni(t ) - - a i B i ( t ) so that (6) holds. A straightforward a p p r o a c h is to set u p a linear system of equations involving a 0 , . . . , an+ 1. However, the following a p p r o a c h yields a simpler solution. Consider the line L l ( t ) that passes t h r o u g h the two B6zier p o i n t s (i - 1 + 1 / 3 , bi_l,1) a n d (i - 1 + 2 / 3 , b;-1,2), a n d the line L2(t ) t h r o u g h the points (i + 1 / 3 , bi,1) and (i + 2 / 3 , b~,2). T h e p o i n t where the two lines intersect is (?, ~), where
-t=i-(wi-Wi_x)/(3wi+
3wi_l)
and
y=d,_a/(2w~_l(w~+Wi_a)
),
TA. Foley / Local control of interval tension
287
where di_ 1 is defined in (7). We claim that if a i =~-1 and N,.(t) = aiBi(t ), then (6) holds. The Btzier points of N~(t) are bi-2,k = bi+l,k+l = 0 for k = 0, 1, 2,
bi - 1,1 = 2 W i _ l ( W i _ 1 + w i ) / d i _ 1, bi- 1,2 = (3w,_ 1 + wi_2)(w,_a + wi)/di_ 1, bi,1 = (3wi+wi+l)(wi-1 + w i ) / d , , bi,2 = 2 W i ( W i _ 1 + wi)/di,
(8)
bi_l, 0 = bi_2, 3 = b i _ 1 , 1 / 2 , bi_ l, 3 = bio = ( b i _ l , 2 + bi,1) / 2
and
bi,3 = bi + 1,0 = b i , 2 / 2 ,
where di is given in (7). Observe that if w i is changed in an interactive application, then the only functions that change are N/_ 1(t), N/(t), N/+ 1(t) and N/+ 2(t). To verify that (6) holds, consider the interior B&ier points for the function g(t) = Ej=0Nj(t) ,+1 on the interval (i, i + 1). These B4zier points are only affected by N~(t) and N,.+l(t ) because the inner B4zier points of N~_l(t ) and N,+2(t ) are zero on that interval. Using (8), the B4zier point for g(t) at i + 1 / 3 is
(3wi +
W i + l ) ( W i _ 1 + w i ) / d i + 2 w i ( w i + w i + l ) / d i = 1.
Similarly, all of the other inner Brzier points for g(t) are equal to one, thus g(t) = 1 and (6) holds. It remains to be shown that N o , . . . , N n÷ 1 is a basis for weighted splines. Consider the n + 2 by n + 2 matrix A defined by -N~(1) No0) 0 A= 0 0 0
NI'(1) Nz'(1) N~(1) N2(1) Na(2) N2(2) 0 : 0 0 0
0 0 0
0 0 N3(2)
0 0 0 :
0 0 0
0 0 0 :
0 0 0 •
0 0 0
Nn_2(n - 1 ) 0 0
Nn_l(n - 1 ) Nn_l(n ) U~'_,(n)
Nn(n - 1 ) N,(n) U'(n)
0 N,,+ l(n )
N'+a(n )
Note that
N,.(i
1)=b,_l,a/2,
N~(i)=(b,_l,2+bia)/2,
N/(i-1)=3N/(i-1),
N,.(i+ l)=bi,2/2.
N / ( i ) = 3(bi,l-b,_l,2)/2
and
N / ( i + 1 ) = - 3 N i ( i + 1). To show that A is nonsingular, consider the matrix E ~that is identical to A except that the first and last rows of E are one third those of A. E is seen to be column diagonally dominant and strictly column diagonally dominant in all but the first and last columns using (8) and the note above. The matrix E is also irreducible, as defined in [Varga '62], thus the irreducibly diagonally dominant matrix E is nonsingular by a result in [Varga '62]. Therefore the matrix A is also nonsingular. Let f ( t ) be any weighted spline. Solve the linear system A C = F for C, where
C = (Co, Cl, ¢2,-.-, on, Cn~-l) t and F= (f'(1), f(1), f(2),..., f ( n -
1), f ( n ) , f ' ( n ) ) t.
T.A. Foley / Local control of interval tension
288
Define the function g ( t ) by n+l g ( t ) = E ciNi(t)" i=0 Since each function N~(t) is a weighted spline, it follows from the linearity of the weighted spline interpolant operator that g ( t ) is also a weighted spline. Since A C = F, it follows that f ( i ) = g(i) for i = 1 . . . . , n, f ' ( 1 ) = g'(1), and f ' ( n ) = g ' ( n ) , h e n c e we have that f ( t ) = g ( t ) for all t I ~< t ~< t n because the weighted spline i n t e r p o l a n t to this d a t a is unique. Therefore, if f ( t ) is a weighted spline, then there exist unique constants q, such that (5) holds. The functions N~(t) are linearly independent, thus No, N a , . . . , N,+a is a basis of weighted splines. F o r an arbitrary partition t 1 < - . . < t,, a similar construction leads to the following definition of N/(t), for i = 0, 1 , . . . , n + 1. For t n o t in the interval [ti_2, //+2], define Ni(t ) = 0, a n d for tj ~< t ~< tj + 1, let u = (t - t j ) / h j and define 3 N i ( t ) = ~_, b j , k ( 3 1 U k ( 1 - U) 3-k, k=O ~'~ where b i_ 2,k = bi+ 1,k + 1 = 0 for k = O, 1, 2,
bi - 1,1
=
(hi_ 2 q- h i _ l ) V i _ l ( V , _ 1 q- v i ) / e i _ l ,
bi _ 1,2 =(/3/_1 + v i ) ( 2 h i _ l V i _ l
+ hi_2vi_ 1 q - h i _ l V i _ 2 ) / e , _ 1,
bi, 1 = (Vi_a + v i ) ( 2 h i v i + hivi+a + h i + a v i ) / e i , bi, 2 = ( h i + h i + l ) v i ( v i - 1
+vi)/ei,
bi-2,3 = bi_l, 0 = h i _ 2 b i _ l a / ( h i _ 2
(9)
+ hi_l),
bi- 1,3 = bi, 0 = (hibi_l, 2 + h i _ l b i , a ) l ( h , _ l
+ hi),
bi, 3 = bi+l, 0 = hi+lbi,2/(hi + hi+l).
1.0
a
.8
,6" ,4" .2" O. .
I.
2.
5.
4.
1.0
b
8
o! O.
-
I.
2.
3.
4.
Fig. 3. These are the normalized weighted spline basis function N2(t ) with weights (a) w i = l for all i; (b) w0=w2=5 and Wl=W3=l.
T.A. Foley / Local control of interval tension
289
where v i = w i / h i and el= [(2hi + h i _ l ) V i + hil)i_l]l.)i+ 1 -k-(hi+ 1 q- 3h i + h,_l)v2i + (2h i + hi+l)Vi_lVi.
(10)
Fig. 3 contains two plots of the normalized weighted spline basis function N2(t ) with different weights using the uniform partition t i = i. In Fig. 3(a), the weights all equal one, thus N2(t ) is siniply a C 2 B-spline. In Fig. 3(b), the weights are w 0 = w 2 = 5, and w 1 = w 3 = 1. The B6zier points are denoted by circles.
4. Control point sums
In computer aided design applications, the normalized weighted spline basis can be used in the same way that B-splines are used. Given control points (xi, Yi) for i = 0, 1 , . . . , n + 1, define the parametric curve ( X ( t ) , Y(t)) for t 1 <~t <~t, by n+l
X(t)=
~ xiN,(t) i=0
n+l
and
Y(t)=
E yiN,(t) •
(11)
i=0
For 3-D space curves, the z-component Z ( t ) is defined similarly. Since N~(t) = 0 for t outside of [ti_ 2, ti+2], if t i <~t <~t~+1, then X ( t ) = X i _ l N i _ l ( t ) "~- x i N i ( t ) ~- x i + l N / + l ( t ) --~ xi+ 2Ni + 2 ( t ), and similarly for Y(t). The curve in (11) is locally controlled in that if the point (x~, y~) is changed, then the only portion of the curve that changes is when t;_ 2 < t < ti+2. If one weight is changed, say wi, then only the basis functions N/_l(t ), N,(t), N/+x(t ) and N/+z(t ) are changed, and as we will show later, the control point curve will change only on the portions where ti_ 2 < t < ti+ 3. Fig. 4(a) shows the parametric curve (11) with all weights equal to one, thus it is a B-spline approximation, and Fig. 4(b) is the parametric curve in (11) with the weights on segments A, B, C and D equal to 8, 40, 40 and 8, respectively, and all other weights equal to one. The uniform parametrization, where t r = i , is used in Fig. 4 through 6. The n = 18 control points (Xl, Yl),.--,(xas, Yas) are denoted by dots and are connected by line segments. There is one control point at (2.0, 0.0) that is obscured by a triangle symbol, and the first and last point are in the lower left at (0.0, 0.0). The points (x0, Y0) = (0.0, - 1 . 0 ) and ( x , + l , Y,,+I) = ( - 0 . 5 , 0.0), which are the linear extensions of the first and last line segments, are not displayed. The triangles represent where one cubic segment of the parametric curve ends and the next begins. More precisely, they represent points on the curve ( X ( t ) , Y(t)) where t = i, for i = 1 , . . . , 1 8 . The triangle points are also interpolation points, in that the curve ( X ( t ) , Y(t)) is the parametric weighted spline interpolant to the triangle data using appropriate first derivative end conditions at (x 1, Yl) and ( x , , y,). The figures that follow adhere to the same notational conventions. Since the functions N,(t) are nonnegative and they sum to one for all tl <~ t <~t,, the functions X ( t ) and Y ( t ) are convex sums and hence the parametric curve ( X ( t ) , Y(t)) satisfies the local convex hull property as described in [Faux et al. '79, Gordon et al. '74]. This is demonstrated in Fig. 4 where the curve defined on t~ ~< t ~< t~+ a does not exceed the convex hull of the four points (xj, yj), j = i - 1, i, i + 1, i + 2. The parametric curve ( X ( t ) , Y(t)) in (11) has C 1 continuity as long as four consecutive control points (xg, yi) are not all equal and the tangent vector is nonzero. Fig. 5 exhibits the behavior when there are two control points at E. The weights and control points in Fig. 5 are the same as in Fig. 4 with the exception that point E is counted twice. The curve comes closer
290
T.A. Foley / Local control o f interval tension
a 2 . 0 -
1.51.0
.5 O.
0'.
I'.
2'.
3'.
4'.
5'.
B 2.01.5-
~
b
C
1.0-
,
i
i
,
~
i
1. 2. 3. 4. 5. Fig. 4. These are the weighted spline basis approximations in (11) to the control points denoted by dots using weights (a) w i =1 for all i (same as using B-splines);'and (b) wi = 8, 40, 40, and 8 on segments A, B, C and D respectively, and 1 otherwise. O.
2 , 0 -
1.51.0-
.~O.
0'.
1'.
2'.
3'.
4'.
5'.
Fig. 5. This is the same as Fig. 4(b) except that the control point E is counted twice. B
a
b
Fig. 6. These are periodic weighted spline basis approximations to the control points using weights (a) wi = 1 for all i (same as using B-splines); and (b) w i ~ 5, 20, 5, and 20 on segments A, B, C and D respectively, and 1 otherwise.
T.A. Foley / Local control of interoal tension
291
to the point E because the convex hull of the neighboring 4 points is narrower near E because of the repeated point. As with B-splines, if three consecutive control points are equal, then the convex hull of the four nearest points is the point itself. In this case the curve must pass through the control point with C O continuity. A common approach in forming the control point curve in (11) is to define the first and last points with a multiplicity of three, thus the curve will interpolate the first and last points as Brzier curves do. If four consecutive control points are collinear, then the curve in (11) is linear on the middle segment and equal to the line segment joining the two interior control points. This follows from the local convex hull property. To generate periodic curves, one approach is to define Wn+k ~ W k + l ,
Xn+k ~ X k + l ,
W-k = Wn-l-k,
Yn+k = Y k + l ,
X-k = Xn-l-k~
Y-k =Yn-l-k
for k = 0, 1, 2. For the curve defined by (11), we have X ( t l ) = X ( t , ) and X ' ( t l ) = X ' ( t n ) and similarly for Y ( t ) , thus the curve is C 1 where the beginning meets the end. If w 1 = w , - 1 , then by (4), the curve is C 2 where the beginning meets the end. Fig. 6 exhibits this behavior using periodic control points and weights. The weights in Fig. 6(a) are all equal to one, hence this is a C 2 B-spline approximation. The weights in Fig. 6(b) are 20.0 on the top and bottom segments B and D, 5 on the neck segments A and C, and one on all other segments. The computation of (11) can be simplified by computing the Brzier points for the parametric curve. The interior Brzier points for (11) on the interval (ti, /i+1) are only dependent on the contributions from N,(t) and N/+l(t ) because the inner Brzier points of N/_l(t ) and N/+.2(t) are zero on this interval. If we denote the control points by Ci = (xi, Yi), or (xi, Yi, zi) in the 3-D case, then for t i ~< t ~< ti+l, the control point curve satisfies n+l
E j=0
3
~---
k=0
_
)3-k
' ,k
'
where u = (t - t i ) / h i and the Brzier points are bi, 1 = C i "}- ( Ci+ 1 - C i ) ( h i q- h i _ l ) l ) i ( o i -I- oi+ 1 ) / e i ,
bi,2 = Ci+ l + ( Ci - Ci+ a)( hi + hi+ a)vi( vi + u i - 1 ) / e i ,
b,,o = b,-1,3 = (
+
(12)
hi +
a n d bi, 3 = bi+l, o, where v i = w i / h i and e i is given in (10). Using (9), (10) and (11), it is straightforward to show that bi, 1 given above is equal to C i ( u i _ 1 + t3i)(2hil3 i q- hil3i+ 1 -t- h i + i v i ) / e i
+Ci+l(hi-1
+hi)ui(ui+oi+l)/ei.
Similarly, bi,2 in (12) can be verified and the remaining Brzier points in (12) follow because each component of the parametric curve is a C 1 function in t. If w} is changed then the Brzier points that change are bi_L0 = bi_2,3, bi_l,1, bi_l, 3 = bi,o, bi,1, bi,2, bi, 3 = bi+l,0 , bi+l,2 and bi+L3 = bi+2,0. Thus if w; is changed, the only portion of the control point curve that changes is when ti_ 2 < t < ti+ 3. Fig. 7 illustrates the Brzier representation of the control point curve using the uniform parametrization with w 3 = 10, w 6 = 100 and w~ = 1, otherwise. The first control point is (1, 1), the control points are denoted by triangles and the Brzier points are denoted by circles. The first and last control points can be considered to be multiple points by simply requiring that bl,0 = C~ and bn-l,3 = Cn. Observe that the inner Brzier points fall on the line segments connecting the control points. Also observe that the curve ( X ( t ) , Y ( t ) ) approximates the middle 60% of the line segment connecting C6 and C 7 because w 6 is large. In general, if w;
T.A. Foley / Local control of interoal tension
292
t
1.
I
I
2.
I
3.
q.
Fig. 7. This is the control point curve with w 3 = 1 0 , w 6 = 1 0 0 and w i = l , otherwise. The control points are denoted by triangles, the B6zier points by circles, and the first point is in the lower left.
I
S.
approaches oo and w~_ 1 and w~+l remain bounded, then the control point curve approaches the line segment connecting bi,0 and bi, 3 for tg ~< t ~< ti+l, where bi,o = Ci + (Ci. 1 - - C i ) h i / ( hi_ 1 -.[-3h i + bi,3=Ci+l + ( C i - C i + l ) h i + l / ( h i - 1
hi.l),
+ 3hi+hi+l).
These follow because, using (12), as w~ approaches oo, biA ~ C~+ (C~+1- C~)(h~-i + h ~ ) / ( h ~ - i + 3h~+ h~+]) , +
bi_l, 2 ~ C i and bi+l, 1 ~ Ci+ 1.
Everything presented so far in this section has dealt with parametric curves. In regards to the limiting results in the previous paragraph, it is interesting to note that as a function on [t~, t~+l], the component Y(t) in (11) has a different visual behavior. If the uniform parametrization is used where t~ = i, then as w i ~ oo, we have that Y ( t ) ~ 0.8y i + 0.2yi+ 1 + 0.6(yi+ 1 - y i ) ( t -
i),
for i ~< t ~< i + 1. Thus Y(t) approaches a line segment with slope 0.6(yg+ 1 -y~), whereas the parametric control point curve approaches a line segment with slope (Y~+I - y i ) / ( X , ÷ l - x ~ ) . For weighted spline interpolation of functional data, [Salkauskas '84] used default weights defined by w
h,
"
(13)
with m = 3, so that V, defined in (1), approximates the L 2 norm of the curvature of f ( t ) . Fig. 8 contains three functional plots of Y(t) in (11) using the weights defined in (13) with m = 0 in part (a), m = 1 in part (b), and m = 3 in part (c). Using m = 3 in (c) yields the same weights as those in [Salkauskas '84], and using m - - 0 in (a) yields a C 2 B-spline approximation because w~ = 1 for all i. There are two additional control points not shown, one at each end, that are linear extensions of the first and last line segments. The points (q, y~) are denoted by dots and
T.A. Foley / Local control of interval tension
293
a .
6. 4. 2. O.
2'.
4'.
6'.
8'. b
.
6. 4. 2. O.
0'.
2'.
4'.
6'.
8'. C
.
6. 4. 2. O.
Fig. 8. These are functional plots of the curve Y(t) in (11) using the weights de-
2'.
4'.
6'.
8'.
fined in (13) with (a) m = 0, (b) m = l and (c) m = 3.
are connected by line segments. The triangle points denote where one cubic segment begins and the other ends.
5. Concluding remarks For the control point sums in Section 4, it is suggested to use the B6zier representation in (12) because it is computationally more efficient than actually computing the B-spline like basis. In fact, we could have directly generated the formulas in (12) without ever constructing the B-spline like basis by using an approach similar to those used in [Boehrn '85] and [Farin '82]. However, there are many other applications in approximation theory that require a B-spline like basis.
294
T.A. Foley / Local control of interoal tension
For three-dimensional surfaces, the weighted bicubic spline in [Foley '87b] is a generalization of the weighted spline interpolant given here. The weights control the tension on the rectangular patches defined by the data. A B-spline like basis for this bivariate interpolant is currently being researched. Unfortunately, the solution is not the tensor product basis N,(s)Nj(t). Surface approximation by
EEz,jNi(s) (t) i
j
will have weights controlling entire rows and columns of the surface. Weighted spline interpolation is a useful and efficient tool in computer aided design when control of tension on intervals connecting interpolation points is needed. The B-spline like basis is well suited for interactive generation of freeform curves, in the same way that B-splines are. The resulting curves satisfy the convex hull property, they are piecewise cubics and the curves can be locally controlled with interval tension in a computationally efficient manner.
Acknowledgements Special thanks are due to the referees for their constructive comments, and to Ralph Carlson and Fred Fritsch of Lawrence Livermore National Laboratory, where a portion of this work was done. This work was supported by the U.S. Department of Energy under contract No. W-7405-ENG-48 and No. DE-FG02-87ER25041.
References Barsky, B.A. (1981), The beta-spline: A local representation based on shape parameters and fundamental geometric measures, Ph.D. dissertation, Dept. of Computer Science, Univ. of Utah. Available in revised form as Computer Graphics and Geometric Modeling using Beta-splines, Springer, Tokyo, 1987. Barsky, B.A. (1984), Exponential and polynomial methods for applying tension to an interpolating spline curve, Computer Vision, Graphics, and Image Processing 27, 1-18. Barsky, B.A. and Beatty, J.C. (1983), Local control of bias and tension in Beta-splines, ACM Trans. Graphics 2, 109-134. Boehm, W. (1985), Curvature continuous curves and surfaces, Computer Aided Geometric Design 2, 313-323. Boehm, W., Farin, G. and Kahmann, J. (1984), A survey of curve and surface methods in CAGD, Computer Aided Geometric Design 1, 1-60. de Boor, C. (1978), A Practical Guide to Splines, Springer, New York. Cohen, E., Lyche, T. and Riesenfeld, R.F. (1980), Discrete B-splines and subdivision techniques in computer aided geometric design and computer graphics, Computer Graphics and Image Processing 14, 87-111. Conte, S.D. and de Boor, C. (1980), Elementary Numerical Analysis, McGraw-Hill, New York. Farin, G.E. (1982), Visually C 2 cubic splines, Computer-aided Design 14, 137-139. Farin, G.E. (1985), Some remarks on V 2 splines, Computer Aided Geometric Design 2, 325-328. Faux, I.E. and Pratt, M.J. (1979), Computational Geometry for Design and Manufacture, Ellis Horwood, Chichester. Foley, T.A. (1987a), Interpolation with interval and point tension controls using cubic weighted p-splines, ACM Trans. Math. Softw. 13, to appear. Foley, T.A. (1987b), Weighted bicubic spline interpolation to rapidly varying data, ACM Trans. Graphics 6, to appear. Gordon, W.J. and Riesenfeld, R.F. (1974), B-spline curves and surfaces, in: Barnhill, R.E. and Riesenfeld, 1LF., eds., Computer Aided Geometric Design, Academic Press, New York. Kochanek, D.H. and Bartels, R.H. (1984), Interpolating splines with local tension, continuity, and bias control, Computer Graphics 18, 33-41. Mortenson, M.E. (1985), Geometric Modeling, Wiley, New York. Nielson, G.M. (1974), Some pi¢cewise polynomial alternatives to splines under tension, in: Barnhin~ R.E. and Riescnfeld, R.F., eds., Computer Aided Geometric Design, Academic Press, New York. Salkauskas, K. (1984), C 1. splines for interpolation of rapidly varying data, Rocky Mountain J. Math. 14, 239-250. Schumaker, L.L. (1981), Spline Functions: Basic Theory, Wiley, New York. Varga, R.S. (1962), Matrix Iterative Analysis, Prentice-Hall, Englewood Cliffs.
281
Local control of interval tension using weighted splines Thomas A. FOLEY * Computer Science Dept., Arizona State University, Tempe, A Z 85287, U.S.A. Received 4 February 1986 Revised 9 March 1987 Abstract. Cubic spline interpolation and B-spline sums are useful and powerful tools in computer aided design. These are extended by weighted cubic splines which have tension controls that allow the user to tighten or loosen the curve on intervals between interpolation points. The weighted spline is a C 1 piecewise cubic that minimizes a variational problem similar to one that a C 2 cubic spline minimizes. A B-spline like basis is constructed for weighted splines where each basis function is nonnegative and nonzero only on four intervals. The basis functions sum up identically to one, thus curves generated by summing control points multiplied by the basis functions have the convex hull property. Different weights are built into the basis functions so that the control point curves are piecewise cubics with local control of interval tension. If all weights axe equal, then the weighted spline is the C 2 cubic spline a n d the basis functions are B-splines. Keywords. Interpolation, approximation, curves, computer aided design, B-splines.
1. Introduction
The theory and application of C 2 cubic spline interpolation to computer aided design are well documented in [Schumaker '81, Boehm et al. '84, de Boor '78, Faux et al. '79]. Piecewise cubic tension splines in [Nielson '74, Kochanek '84] provide controls for making a curve tighter or looser at the interpolation points. The weighted splines in [Salkauskas '84] and the exponential-based splines under tension in [Barsky '84] provide tension controls that make the curve tighter or looser on intervals between interpolation points, rather than at the discrete interpolation points. The weighted splines are piecewise cubics, while the splines under tension are piecewise exponential functions. Weighted splines are discussed in Section 2 and useful end conditions are developed that parallel those for C 2 cubic spline interpolants. B-splines, described in [Boehm et al. '84, de Boor '78, Faux et al. '79, Cohen et al. '80, Gordon et al. '74], are a useful and powerful tool in computer aided design. They form a basis for the space of C 2 cubic splines, and each B-spline is a nonnegative cubic spline that is nonzero only on four intervals. The B-splines form a partition of unity, that is, they sum up to one. Curves generated by summing control points multiplied by the B-splines have some very desirable shape properties, including the local convex hull property. The Beta-splines in [Barsky '81] and [Barsky et al. '83] are similar to B-splines in that they are nonnegative piecewise cubics, they are nonzero only on four intervals and they sum to one. Two parameters are incorporated into the Beta-splines that control the tension at a point and the bias of a curve generated by summing control points multiplied by the Beta-splines. If * A portion of this work was supported by the U.S. Department of Energy at Lawrence Livermore National Laboratory under contract No. W-7405-ENG-48 and No. DE-FG02-87 ER 25041. 0167-8396/86/$3.50 © 1986, Elsevier Science Publishers B.V. (North-Holland)
T.A. Foley / Local control of interval tension
282
different tension values are used at various points, then an interpolation between parameters is performed that yields a higher degree piecewise polynomial. In Section 3, a B-spline like basis is constructed for weighted splines where each basis function is nonnegative and nonzero on only four intervals. This basis is normalized so that it forms a partition of unity. Section 4 uses this basis to form a control point sum will yield a curve that has the convex hull property. The different weights are built into the basis functions so that the resulting control point curve is a piecewise cubic with local control of interval tension.
2. Weighted spline interpolation This section describes weighted spline interpolation with different end conditions. The weighted spline developed in [Salkauskas '84] did not involve end conditions, although his approach is equivalent to using Type II or natural end conditions as described below. Given a = t 1 < • • • < t n = b, Yl,-.-, Yn, and a positive integrable weight function w(t), the approach used in [Salkauskas '84] finds the C a interpolant to the given data that minimizes the quantity
fabW(t)[f"(t)] 2 dt. When w(t) is a piecewise constant, the solution is a piecewise cubic function and the minimization problem reduces to solving a tridiagonal system of equations whose unknowns are the first derivatives at the data points. The minimizing function is the piecewise cubic Hermite interpolant to the data using the solution of the tridiagonal system for the derivatives. Some motivation for the weighted cubic spline is that the C 2 natural cubic spline is the unique interpolating function that minimizes lab[ f " ( t ) ] 2 dt subject to certain continuity conditions. Minimizing this global measure of curvature is effective if the data set is well behaved because it places equal weight on all intervals. When using the weighted spline, if w(t) is large on one interval relative to the other intervals, then f " ( t ) is forced to be small in magnitude on that interval, hence f ( t ) will be more linear on that interval. Similarly, smaller relative weights allow f " ( t ) to take on larger values. The term relative is used because if all weights are multiplied by a positive constant, then the resulting weighted spline would be the same. As noted in [Salkauskas '84], if the weight function w(t) is constant, then the weighted spline is identical to the C 2 natural cubic spline. Given a = t 1 < t 2 < • • • < t n = b, Yl,..., Yn, and w i > 0 for i = 1 . . . . , n - 1, define w(t) = w i if t i ~< t < ti+ 1. Using a piecewise constant weight function, the weighted spline interpolant is the C 1 piecewise cubic function f ( t ) that minimizes
V= f bw(t)[f"(t)]2 dt
(1)
subject to f(ti) =Yi for i = 1 , . . . , n, and one of the following end conditions: (I) f ' ( a ) = r n l (II) f " ( a + ) = 0 (III) f ' ( a ) = f ' ( b )
and
f'(b)=m~,
and
f"(b-)=0,
and
or
Wxf"(a+)=w~_lf"(b-).
Type I are called first derivative end conditions, Type II are called natural end conditions, and Type III are called periodic end conditions. First derivative and natural end conditions are
T.A. Foley / Local control of interval tension
283
identical with those used by a standard C 2 cubic spline. The rationale for periodic end conditions is discussed later in this section. A constructive approach is used to verify that a unique weighted spline exists, subject to the gix/en end conditions. This approach will yield a method for computing the weighted spline using cardinal piecewise cubic Hermite basis functions. In [de Boor '78], it is shown that given Yi and m i, for i = 1 , . . . , n, then there exists a unique C 1 piecewise cubic that satisfies f ( t i ) =y~ and f ' ( t i ) = m i for i = 1 , . . . , n. Furthermore, f ( t ) can be represented by n
f(t)=
(2)
~,yiPi(t)+miDi(t), i=1
where Pi(t) and Di(t) are the cardinal piecewise cubic Hermite functions that satisfy Pi(tj)=Si,j,
Pi'(tj)=O,
Di(tj)=O,
D;(tj)=~ij,
for i, j = 1 . . . . , n. Since the functions Pi(t) and D~(t) are zero outside of [t~_ 1, t;+l], equation (2) simplifies computationally to f ( t ) = y i P i ( t ) + m i D i ( t ) + yi+lPi+l(t) + m i + a D i + l ( t )
for t i <<.t <~t~+ 1. The unknowns for the weighted spline are the derivative values m~. Necessary conditions for minimiT.ing V in (1) are OV Om k
----~0
for k = 2 , . . . , n - 1. It follows that m k _ l W k _ l h k l _ l + m k 2 ( W k _ l h k l _ l + Wkh-~ 1) + mk+lWkh-~ 1 = 3Wk_l(Yk--Yk_a)h-k21 + 3W/,(yk+ 1 --yl,)h-k 2,
(3)
for k = 2 , . . . , n - 1, where h i = ti+ 1 - t r The n - 2 equations in (3) involve the n unknowns m l , . . . , m,,. The remaining two equations are generated from the given end conditions. For Type I or first derivative end conditions, the other equations are m 1 = value of f ' ( a )
and
m , = value of f ' (b).
The two additional equations for Type II or natural end conditions are m 1 2 h l 1 + m E h l 1 = 3 ( y 2-y1)h1-2,
and --1
-1
--2
m n _ l h n _ 1 + m . 2 h . _ 1 = 3(y. - y n _ l ) h ~ _ l .
For Type III or periodic end conditions, the equations are m l ( 2 W a h f I + 2 W n _ l h n- _1 l ) + m 2 w l h f I + m . . l W . _ l h - n l _ l
= 3(y 2
--
-2
ya)hl 2+ 3(y,,-y~_l)hn_l,
and m 1 = m no
The linear systems of equations that occur when using (3) and either Type I or Type II end conditions are tridiagonal and diagonally dominant. N o t only do they have unique solutions,
284
T.A. Foley / Local control of interval tension
v
w
Fig. 1. This is a parametric weighted spline interpolant with natural end conditions using weights w i = 10, 1, 20, 1, 10.
but they can be solved efficiently using a standard tridiagonal system solver, such as the one in [Conte et al. '80]. If the two equations for Type III end conditions are placed in the first and last rows of the linear system, then with a slight modification of the tridiagonal system solver, the unique solution can be efficiently computed. The previous linear system of equations are also sufficient conditions for minimizing V. An argument similar to those in [Salkauskas '84] and [Foley '87a] can be used to verify this. The weighted splines are C ~ functions in general, but as shown in [Salkauskas '84], if wi_ 1 = wi, then the weighted spline is C 2 at t~. This result follows from the fact that wif
w
+
(t i ) = W i _ l f
.,v
(t;).
(4)
This behavior of the second derivative is the rationale for the Type III end conditions that are termed periodic. In fact, if all the weights are the same, then the weighted spline with specified end conditions is identical to the C 2 standard cubic spline with the same end conditions. For parametric interpolation, assume that we are given (x;, yg) for i = 1 , . . . , n. Let t 1 < --< tn, let X ( t ) be the weighted spline interpolant to the data (tg, x~) and let Y(t) be the weighted spline interpolant to the data ( ti, Y i ) , i = 1 . . . . . n. The parametric curve ( X ( t ), Y( t )), t 1 ~ t ~< t, passes through the given data. For 3-D space curves, the z-component Z ( t ) is computed similarly. Fig. 1 is the parametric weighted spline using natural end conditions and weights w~ = 10, 1, 20, 1, 10 using a uniform parametrization with t i = i. This example shows how easy it is for the curve to be nearly linear on one interval and have high curvature on bordering intervals. Some visual guidelines on the effects of the weights are given in Fig. 2 using the uniform parametrization. All of the curves interpolate the six points (0, - 1 ) , (0, 0), (0, 1), (1, 1), (1, 0) and (1, - 1) using w i = 1 for i ~ 3. F r o m top to bottom on the middle interval, the weighted spline interpolant uses w 3 = 1, 5, 20 and 100, respectively. The v-splines in [Nielson '74] and the Beta-splines in [Barsky et al. '83] are C 1 as functions, but they have continuous curvature as parametric curves. The parametric weighted splines are not G 2, that is, the curvature is not continuous at the interpolation points. If they were G 2 continuous, then Y"(t~-)-
Y"(t~-)=vY'(ti)
for some value of v, as shown in [Farin '86]. Consider the example where t i = 1, 2, 3, 4, 5, y; = 0, 1, 14, 4, 0, and w t = 11, 1, 2, 1. The natural weighted spline then has derivative values of m i = 0, 3, 0, - 1 2 , 0. Since Y ' ( 3 ) = 0 and Y ' ( 3 ÷) :~ Y " ( 3 - ) , the weighted spline is not G 2 in general. As noted earlier, if wi -- wi_l, then each component of the curve is C 2, and thus if the tangent vector is nonzero, the parametric curve has continuous curvature at t i. Although the weighted spline is not generally G 2, it is generally 'smoother' than many C 1 curves in that it has a continuous osculating plane everywhere. A discussion of the osculating plane at a point on a
T.A. Foley / Local control of interval tension
285
• 2 ~ ~'.~._.-:...-_ ..........................
:.3
_-:._-_.__-:. -..
S
',,'i1
.6
:'I
l
J i~ .2
1 .2
8.
Iq
i
.2
.
I
I
I
.6
.8
1.g
..
Fig. 2. T h e s e are w e i g h t e d s p l i n e i n t e r p o l a n t s u s i n g w i = 1 for i ~ 3, a n d f r o m t o p t o b o t t o m o n t h e m i d d l e i n t e r v a l , t h e v a l u e s o f w 3 are 1, 5, 20 a n d 100, respectively.
3-D space curve is given in [Mortenson '85, p. 271] where an e q u a t i o n for the p l a n e at
(X(t), Y(t), Z(t)) is given by [ x - X ( t ) ] [ Y ' ( t ) Z " ( t ) - Y"(t)Z'(t)] -[y- Y(t)][X'(t)Z"(t)-X"(t)Z'(t)] + [z- Z(/)][X'(t)Y"(t)X"(t)Y'(t)]
=0.
U s i n g (4), it follows that the osculating plane at t = t~ is continuous.
3. A B-spline like basis for weighted splines W e present a basis No(t),..., Nn+l(t) such that if f(t) is any weighted spline, that is, a C 1 piecewise cubic that minimizes (1), then there exists u n i q u e constants c; such that
f(t)
is
n+l
f(t)=
E
ciNi(t)
(5)
i=0
for all t I ~< t ~< t,. Each basis function N~(t) is a weighted spline, N ~ ( t ) > / 0 for all t, a n d N~(t) = 0 for all t outside of [ti_ 2, ti+2]. T h e basis functions will be n o r m a l i z e d or f o r m a partition of unity, that is, n+l
N/(t) = 1
(6)
i=0
for all t I ~< t ~< t,. F o r notational convenience, we will first c o n s t r u c t the basis using a u n i f o r m p a r a m e t r i z a t i o n where t, = i a n d then we state the results for an arbitrary parametrization. W e first construct basis functions Bg(t) and later define the n o r m a l i z e d functions N,.(t) in terms of them. G i v e n weights Wa,... , Wn_l, define the weights w 0 = w_ 1 = w_z = Wl a n d w , = w , + 1 =
286
T.A. Foley / Local control of interval tenswn
w~+ 2 = wn_ 1. F o r Bi(t ) to be a C 1 piecewise cubic a n d have c o m p a c t s u p p o r t on [i - 2, i + 2], if follows that
0 = n i ( i - 2) = B/' ( i - 2) = B i ( i + 2) = n/' (i W 2). Since B~'(i - 2 - ) = 0, b y (4) it follows t h a t B f ' ( i - 2) = 0. Similarly, B~'(i + 2) = 0. Since B~(t) is a C 1 piecewise cubic, it can be represented in terms of its B6zier points for j <<.t <~j + l by 3
3 ) u (1 _ u)3-k,
Bi(t) = E bj,k( k k=O
where u = t - j . Based o n the zero properties of B~(t) discussed in the previous paragraph, we have that bi_2, k = bi+l, k = 0 for k = 0, 1, 2. F o r B~(t) to be C ~ at t = i - 1, i, a n d i + 1, we m u s t have that bi-2,3 = bi-l,O ~ ( b i - 2 , 2 q- b / - 1 , 1 ) / 2 b i _ l , 3 -~ bi, 0 = ( b i _ l , 2 q- b i , 1 ) / 2 ,
-~- b i - l , 1 / 2 ,
and
hi,3 = bi /~alo = ( bi,2 + bi+1,1)/2 = bi,2/2.
T h u s the only u n k n o w n B6zier points are bi-l,1, 'b'i2a,2, bi,1 a n d bi, 2. If we arbitrarily set bi-1,~ = 1, then we can c o m p u t e the three u n k n o w n B6zier p o i n t s by requiring that (4) holds at the points i - 1, i a n d i + 1. I n terms of the B6zier points, these three equations are
6 w j ( b j , o - 2bja + bj,2) = 6 w j - a ( b j - l a -
2bj_x, 2 + bj_l,3),
for j = i - 1, i, i + 1. In terms of the three u n k n o w n B6zier points, these three equations simplify to
2wi_lbi_l, 2 = 3wi_ 1 + wi_ 2, (w, + 3wi_l)bi_l, 2 + ( - 3 w i - wi_l)bi, 1 + 2wibi,2 = 2 w i _ l , 2wibi, 1 + ( - 3 w i - wi+l)bi, 2 = O. T h e u n i q u e solution of the above system is hi_l,2 = ( 3 w i - 1 q- w i _ 2 ) / ( 2 w i _ l )
,
bi,1 = (wi+ 1 + 3 w i ) d i _ l / ( 2 d i ) ,
and
bi, 2 = w i d i _ l / (
W i_ldi) ,
where
(7)
d i = 3 w i w i + 1 + W i _ l W i + 1 q- 5w2i --t- 3 W i _ l W i.
Observe that Bi(t ) >1 0 because all of its B6zier points are nonnegative. Before we show that B o ( t ) , . . . , B , + l ( t ) is the basis for w e i g h t e d splines, we first determine a i > 0 and define normalized functions Ni(t ) - - a i B i ( t ) so that (6) holds. A straightforward a p p r o a c h is to set u p a linear system of equations involving a 0 , . . . , an+ 1. However, the following a p p r o a c h yields a simpler solution. Consider the line L l ( t ) that passes t h r o u g h the two B6zier p o i n t s (i - 1 + 1 / 3 , bi_l,1) a n d (i - 1 + 2 / 3 , b;-1,2), a n d the line L2(t ) t h r o u g h the points (i + 1 / 3 , bi,1) and (i + 2 / 3 , b~,2). T h e p o i n t where the two lines intersect is (?, ~), where
-t=i-(wi-Wi_x)/(3wi+
3wi_l)
and
y=d,_a/(2w~_l(w~+Wi_a)
),
TA. Foley / Local control of interval tension
287
where di_ 1 is defined in (7). We claim that if a i =~-1 and N,.(t) = aiBi(t ), then (6) holds. The Btzier points of N~(t) are bi-2,k = bi+l,k+l = 0 for k = 0, 1, 2,
bi - 1,1 = 2 W i _ l ( W i _ 1 + w i ) / d i _ 1, bi- 1,2 = (3w,_ 1 + wi_2)(w,_a + wi)/di_ 1, bi,1 = (3wi+wi+l)(wi-1 + w i ) / d , , bi,2 = 2 W i ( W i _ 1 + wi)/di,
(8)
bi_l, 0 = bi_2, 3 = b i _ 1 , 1 / 2 , bi_ l, 3 = bio = ( b i _ l , 2 + bi,1) / 2
and
bi,3 = bi + 1,0 = b i , 2 / 2 ,
where di is given in (7). Observe that if w i is changed in an interactive application, then the only functions that change are N/_ 1(t), N/(t), N/+ 1(t) and N/+ 2(t). To verify that (6) holds, consider the interior B&ier points for the function g(t) = Ej=0Nj(t) ,+1 on the interval (i, i + 1). These B4zier points are only affected by N~(t) and N,.+l(t ) because the inner B4zier points of N~_l(t ) and N,+2(t ) are zero on that interval. Using (8), the B4zier point for g(t) at i + 1 / 3 is
(3wi +
W i + l ) ( W i _ 1 + w i ) / d i + 2 w i ( w i + w i + l ) / d i = 1.
Similarly, all of the other inner Brzier points for g(t) are equal to one, thus g(t) = 1 and (6) holds. It remains to be shown that N o , . . . , N n÷ 1 is a basis for weighted splines. Consider the n + 2 by n + 2 matrix A defined by -N~(1) No0) 0 A= 0 0 0
NI'(1) Nz'(1) N~(1) N2(1) Na(2) N2(2) 0 : 0 0 0
0 0 0
0 0 N3(2)
0 0 0 :
0 0 0
0 0 0 :
0 0 0 •
0 0 0
Nn_2(n - 1 ) 0 0
Nn_l(n - 1 ) Nn_l(n ) U~'_,(n)
Nn(n - 1 ) N,(n) U'(n)
0 N,,+ l(n )
N'+a(n )
Note that
N,.(i
1)=b,_l,a/2,
N~(i)=(b,_l,2+bia)/2,
N/(i-1)=3N/(i-1),
N,.(i+ l)=bi,2/2.
N / ( i ) = 3(bi,l-b,_l,2)/2
and
N / ( i + 1 ) = - 3 N i ( i + 1). To show that A is nonsingular, consider the matrix E ~that is identical to A except that the first and last rows of E are one third those of A. E is seen to be column diagonally dominant and strictly column diagonally dominant in all but the first and last columns using (8) and the note above. The matrix E is also irreducible, as defined in [Varga '62], thus the irreducibly diagonally dominant matrix E is nonsingular by a result in [Varga '62]. Therefore the matrix A is also nonsingular. Let f ( t ) be any weighted spline. Solve the linear system A C = F for C, where
C = (Co, Cl, ¢2,-.-, on, Cn~-l) t and F= (f'(1), f(1), f(2),..., f ( n -
1), f ( n ) , f ' ( n ) ) t.
T.A. Foley / Local control of interval tension
288
Define the function g ( t ) by n+l g ( t ) = E ciNi(t)" i=0 Since each function N~(t) is a weighted spline, it follows from the linearity of the weighted spline interpolant operator that g ( t ) is also a weighted spline. Since A C = F, it follows that f ( i ) = g(i) for i = 1 . . . . , n, f ' ( 1 ) = g'(1), and f ' ( n ) = g ' ( n ) , h e n c e we have that f ( t ) = g ( t ) for all t I ~< t ~< t n because the weighted spline i n t e r p o l a n t to this d a t a is unique. Therefore, if f ( t ) is a weighted spline, then there exist unique constants q, such that (5) holds. The functions N~(t) are linearly independent, thus No, N a , . . . , N,+a is a basis of weighted splines. F o r an arbitrary partition t 1 < - . . < t,, a similar construction leads to the following definition of N/(t), for i = 0, 1 , . . . , n + 1. For t n o t in the interval [ti_2, //+2], define Ni(t ) = 0, a n d for tj ~< t ~< tj + 1, let u = (t - t j ) / h j and define 3 N i ( t ) = ~_, b j , k ( 3 1 U k ( 1 - U) 3-k, k=O ~'~ where b i_ 2,k = bi+ 1,k + 1 = 0 for k = O, 1, 2,
bi - 1,1
=
(hi_ 2 q- h i _ l ) V i _ l ( V , _ 1 q- v i ) / e i _ l ,
bi _ 1,2 =(/3/_1 + v i ) ( 2 h i _ l V i _ l
+ hi_2vi_ 1 q - h i _ l V i _ 2 ) / e , _ 1,
bi, 1 = (Vi_a + v i ) ( 2 h i v i + hivi+a + h i + a v i ) / e i , bi, 2 = ( h i + h i + l ) v i ( v i - 1
+vi)/ei,
bi-2,3 = bi_l, 0 = h i _ 2 b i _ l a / ( h i _ 2
(9)
+ hi_l),
bi- 1,3 = bi, 0 = (hibi_l, 2 + h i _ l b i , a ) l ( h , _ l
+ hi),
bi, 3 = bi+l, 0 = hi+lbi,2/(hi + hi+l).
1.0
a
.8
,6" ,4" .2" O. .
I.
2.
5.
4.
1.0
b
8
o! O.
-
I.
2.
3.
4.
Fig. 3. These are the normalized weighted spline basis function N2(t ) with weights (a) w i = l for all i; (b) w0=w2=5 and Wl=W3=l.
T.A. Foley / Local control of interval tension
289
where v i = w i / h i and el= [(2hi + h i _ l ) V i + hil)i_l]l.)i+ 1 -k-(hi+ 1 q- 3h i + h,_l)v2i + (2h i + hi+l)Vi_lVi.
(10)
Fig. 3 contains two plots of the normalized weighted spline basis function N2(t ) with different weights using the uniform partition t i = i. In Fig. 3(a), the weights all equal one, thus N2(t ) is siniply a C 2 B-spline. In Fig. 3(b), the weights are w 0 = w 2 = 5, and w 1 = w 3 = 1. The B6zier points are denoted by circles.
4. Control point sums
In computer aided design applications, the normalized weighted spline basis can be used in the same way that B-splines are used. Given control points (xi, Yi) for i = 0, 1 , . . . , n + 1, define the parametric curve ( X ( t ) , Y(t)) for t 1 <~t <~t, by n+l
X(t)=
~ xiN,(t) i=0
n+l
and
Y(t)=
E yiN,(t) •
(11)
i=0
For 3-D space curves, the z-component Z ( t ) is defined similarly. Since N~(t) = 0 for t outside of [ti_ 2, ti+2], if t i <~t <~t~+1, then X ( t ) = X i _ l N i _ l ( t ) "~- x i N i ( t ) ~- x i + l N / + l ( t ) --~ xi+ 2Ni + 2 ( t ), and similarly for Y(t). The curve in (11) is locally controlled in that if the point (x~, y~) is changed, then the only portion of the curve that changes is when t;_ 2 < t < ti+2. If one weight is changed, say wi, then only the basis functions N/_l(t ), N,(t), N/+x(t ) and N/+z(t ) are changed, and as we will show later, the control point curve will change only on the portions where ti_ 2 < t < ti+ 3. Fig. 4(a) shows the parametric curve (11) with all weights equal to one, thus it is a B-spline approximation, and Fig. 4(b) is the parametric curve in (11) with the weights on segments A, B, C and D equal to 8, 40, 40 and 8, respectively, and all other weights equal to one. The uniform parametrization, where t r = i , is used in Fig. 4 through 6. The n = 18 control points (Xl, Yl),.--,(xas, Yas) are denoted by dots and are connected by line segments. There is one control point at (2.0, 0.0) that is obscured by a triangle symbol, and the first and last point are in the lower left at (0.0, 0.0). The points (x0, Y0) = (0.0, - 1 . 0 ) and ( x , + l , Y,,+I) = ( - 0 . 5 , 0.0), which are the linear extensions of the first and last line segments, are not displayed. The triangles represent where one cubic segment of the parametric curve ends and the next begins. More precisely, they represent points on the curve ( X ( t ) , Y(t)) where t = i, for i = 1 , . . . , 1 8 . The triangle points are also interpolation points, in that the curve ( X ( t ) , Y(t)) is the parametric weighted spline interpolant to the triangle data using appropriate first derivative end conditions at (x 1, Yl) and ( x , , y,). The figures that follow adhere to the same notational conventions. Since the functions N,(t) are nonnegative and they sum to one for all tl <~ t <~t,, the functions X ( t ) and Y ( t ) are convex sums and hence the parametric curve ( X ( t ) , Y(t)) satisfies the local convex hull property as described in [Faux et al. '79, Gordon et al. '74]. This is demonstrated in Fig. 4 where the curve defined on t~ ~< t ~< t~+ a does not exceed the convex hull of the four points (xj, yj), j = i - 1, i, i + 1, i + 2. The parametric curve ( X ( t ) , Y(t)) in (11) has C 1 continuity as long as four consecutive control points (xg, yi) are not all equal and the tangent vector is nonzero. Fig. 5 exhibits the behavior when there are two control points at E. The weights and control points in Fig. 5 are the same as in Fig. 4 with the exception that point E is counted twice. The curve comes closer
290
T.A. Foley / Local control o f interval tension
a 2 . 0 -
1.51.0
.5 O.
0'.
I'.
2'.
3'.
4'.
5'.
B 2.01.5-
~
b
C
1.0-
,
i
i
,
~
i
1. 2. 3. 4. 5. Fig. 4. These are the weighted spline basis approximations in (11) to the control points denoted by dots using weights (a) w i =1 for all i (same as using B-splines);'and (b) wi = 8, 40, 40, and 8 on segments A, B, C and D respectively, and 1 otherwise. O.
2 , 0 -
1.51.0-
.~O.
0'.
1'.
2'.
3'.
4'.
5'.
Fig. 5. This is the same as Fig. 4(b) except that the control point E is counted twice. B
a
b
Fig. 6. These are periodic weighted spline basis approximations to the control points using weights (a) wi = 1 for all i (same as using B-splines); and (b) w i ~ 5, 20, 5, and 20 on segments A, B, C and D respectively, and 1 otherwise.
T.A. Foley / Local control of interoal tension
291
to the point E because the convex hull of the neighboring 4 points is narrower near E because of the repeated point. As with B-splines, if three consecutive control points are equal, then the convex hull of the four nearest points is the point itself. In this case the curve must pass through the control point with C O continuity. A common approach in forming the control point curve in (11) is to define the first and last points with a multiplicity of three, thus the curve will interpolate the first and last points as Brzier curves do. If four consecutive control points are collinear, then the curve in (11) is linear on the middle segment and equal to the line segment joining the two interior control points. This follows from the local convex hull property. To generate periodic curves, one approach is to define Wn+k ~ W k + l ,
Xn+k ~ X k + l ,
W-k = Wn-l-k,
Yn+k = Y k + l ,
X-k = Xn-l-k~
Y-k =Yn-l-k
for k = 0, 1, 2. For the curve defined by (11), we have X ( t l ) = X ( t , ) and X ' ( t l ) = X ' ( t n ) and similarly for Y ( t ) , thus the curve is C 1 where the beginning meets the end. If w 1 = w , - 1 , then by (4), the curve is C 2 where the beginning meets the end. Fig. 6 exhibits this behavior using periodic control points and weights. The weights in Fig. 6(a) are all equal to one, hence this is a C 2 B-spline approximation. The weights in Fig. 6(b) are 20.0 on the top and bottom segments B and D, 5 on the neck segments A and C, and one on all other segments. The computation of (11) can be simplified by computing the Brzier points for the parametric curve. The interior Brzier points for (11) on the interval (ti, /i+1) are only dependent on the contributions from N,(t) and N/+l(t ) because the inner Brzier points of N/_l(t ) and N/+.2(t) are zero on this interval. If we denote the control points by Ci = (xi, Yi), or (xi, Yi, zi) in the 3-D case, then for t i ~< t ~< ti+l, the control point curve satisfies n+l
E j=0
3
~---
k=0
_
)3-k
' ,k
'
where u = (t - t i ) / h i and the Brzier points are bi, 1 = C i "}- ( Ci+ 1 - C i ) ( h i q- h i _ l ) l ) i ( o i -I- oi+ 1 ) / e i ,
bi,2 = Ci+ l + ( Ci - Ci+ a)( hi + hi+ a)vi( vi + u i - 1 ) / e i ,
b,,o = b,-1,3 = (
+
(12)
hi +
a n d bi, 3 = bi+l, o, where v i = w i / h i and e i is given in (10). Using (9), (10) and (11), it is straightforward to show that bi, 1 given above is equal to C i ( u i _ 1 + t3i)(2hil3 i q- hil3i+ 1 -t- h i + i v i ) / e i
+Ci+l(hi-1
+hi)ui(ui+oi+l)/ei.
Similarly, bi,2 in (12) can be verified and the remaining Brzier points in (12) follow because each component of the parametric curve is a C 1 function in t. If w} is changed then the Brzier points that change are bi_L0 = bi_2,3, bi_l,1, bi_l, 3 = bi,o, bi,1, bi,2, bi, 3 = bi+l,0 , bi+l,2 and bi+L3 = bi+2,0. Thus if w; is changed, the only portion of the control point curve that changes is when ti_ 2 < t < ti+ 3. Fig. 7 illustrates the Brzier representation of the control point curve using the uniform parametrization with w 3 = 10, w 6 = 100 and w~ = 1, otherwise. The first control point is (1, 1), the control points are denoted by triangles and the Brzier points are denoted by circles. The first and last control points can be considered to be multiple points by simply requiring that bl,0 = C~ and bn-l,3 = Cn. Observe that the inner Brzier points fall on the line segments connecting the control points. Also observe that the curve ( X ( t ) , Y ( t ) ) approximates the middle 60% of the line segment connecting C6 and C 7 because w 6 is large. In general, if w;
T.A. Foley / Local control of interoal tension
292
t
1.
I
I
2.
I
3.
q.
Fig. 7. This is the control point curve with w 3 = 1 0 , w 6 = 1 0 0 and w i = l , otherwise. The control points are denoted by triangles, the B6zier points by circles, and the first point is in the lower left.
I
S.
approaches oo and w~_ 1 and w~+l remain bounded, then the control point curve approaches the line segment connecting bi,0 and bi, 3 for tg ~< t ~< ti+l, where bi,o = Ci + (Ci. 1 - - C i ) h i / ( hi_ 1 -.[-3h i + bi,3=Ci+l + ( C i - C i + l ) h i + l / ( h i - 1
hi.l),
+ 3hi+hi+l).
These follow because, using (12), as w~ approaches oo, biA ~ C~+ (C~+1- C~)(h~-i + h ~ ) / ( h ~ - i + 3h~+ h~+]) , +
bi_l, 2 ~ C i and bi+l, 1 ~ Ci+ 1.
Everything presented so far in this section has dealt with parametric curves. In regards to the limiting results in the previous paragraph, it is interesting to note that as a function on [t~, t~+l], the component Y(t) in (11) has a different visual behavior. If the uniform parametrization is used where t~ = i, then as w i ~ oo, we have that Y ( t ) ~ 0.8y i + 0.2yi+ 1 + 0.6(yi+ 1 - y i ) ( t -
i),
for i ~< t ~< i + 1. Thus Y(t) approaches a line segment with slope 0.6(yg+ 1 -y~), whereas the parametric control point curve approaches a line segment with slope (Y~+I - y i ) / ( X , ÷ l - x ~ ) . For weighted spline interpolation of functional data, [Salkauskas '84] used default weights defined by w
h,
"
(13)
with m = 3, so that V, defined in (1), approximates the L 2 norm of the curvature of f ( t ) . Fig. 8 contains three functional plots of Y(t) in (11) using the weights defined in (13) with m = 0 in part (a), m = 1 in part (b), and m = 3 in part (c). Using m = 3 in (c) yields the same weights as those in [Salkauskas '84], and using m - - 0 in (a) yields a C 2 B-spline approximation because w~ = 1 for all i. There are two additional control points not shown, one at each end, that are linear extensions of the first and last line segments. The points (q, y~) are denoted by dots and
T.A. Foley / Local control of interval tension
293
a .
6. 4. 2. O.
2'.
4'.
6'.
8'. b
.
6. 4. 2. O.
0'.
2'.
4'.
6'.
8'. C
.
6. 4. 2. O.
Fig. 8. These are functional plots of the curve Y(t) in (11) using the weights de-
2'.
4'.
6'.
8'.
fined in (13) with (a) m = 0, (b) m = l and (c) m = 3.
are connected by line segments. The triangle points denote where one cubic segment begins and the other ends.
5. Concluding remarks For the control point sums in Section 4, it is suggested to use the B6zier representation in (12) because it is computationally more efficient than actually computing the B-spline like basis. In fact, we could have directly generated the formulas in (12) without ever constructing the B-spline like basis by using an approach similar to those used in [Boehrn '85] and [Farin '82]. However, there are many other applications in approximation theory that require a B-spline like basis.
294
T.A. Foley / Local control of interoal tension
For three-dimensional surfaces, the weighted bicubic spline in [Foley '87b] is a generalization of the weighted spline interpolant given here. The weights control the tension on the rectangular patches defined by the data. A B-spline like basis for this bivariate interpolant is currently being researched. Unfortunately, the solution is not the tensor product basis N,(s)Nj(t). Surface approximation by
EEz,jNi(s) (t) i
j
will have weights controlling entire rows and columns of the surface. Weighted spline interpolation is a useful and efficient tool in computer aided design when control of tension on intervals connecting interpolation points is needed. The B-spline like basis is well suited for interactive generation of freeform curves, in the same way that B-splines are. The resulting curves satisfy the convex hull property, they are piecewise cubics and the curves can be locally controlled with interval tension in a computationally efficient manner.
Acknowledgements Special thanks are due to the referees for their constructive comments, and to Ralph Carlson and Fred Fritsch of Lawrence Livermore National Laboratory, where a portion of this work was done. This work was supported by the U.S. Department of Energy under contract No. W-7405-ENG-48 and No. DE-FG02-87ER25041.
References Barsky, B.A. (1981), The beta-spline: A local representation based on shape parameters and fundamental geometric measures, Ph.D. dissertation, Dept. of Computer Science, Univ. of Utah. Available in revised form as Computer Graphics and Geometric Modeling using Beta-splines, Springer, Tokyo, 1987. Barsky, B.A. (1984), Exponential and polynomial methods for applying tension to an interpolating spline curve, Computer Vision, Graphics, and Image Processing 27, 1-18. Barsky, B.A. and Beatty, J.C. (1983), Local control of bias and tension in Beta-splines, ACM Trans. Graphics 2, 109-134. Boehm, W. (1985), Curvature continuous curves and surfaces, Computer Aided Geometric Design 2, 313-323. Boehm, W., Farin, G. and Kahmann, J. (1984), A survey of curve and surface methods in CAGD, Computer Aided Geometric Design 1, 1-60. de Boor, C. (1978), A Practical Guide to Splines, Springer, New York. Cohen, E., Lyche, T. and Riesenfeld, R.F. (1980), Discrete B-splines and subdivision techniques in computer aided geometric design and computer graphics, Computer Graphics and Image Processing 14, 87-111. Conte, S.D. and de Boor, C. (1980), Elementary Numerical Analysis, McGraw-Hill, New York. Farin, G.E. (1982), Visually C 2 cubic splines, Computer-aided Design 14, 137-139. Farin, G.E. (1985), Some remarks on V 2 splines, Computer Aided Geometric Design 2, 325-328. Faux, I.E. and Pratt, M.J. (1979), Computational Geometry for Design and Manufacture, Ellis Horwood, Chichester. Foley, T.A. (1987a), Interpolation with interval and point tension controls using cubic weighted p-splines, ACM Trans. Math. Softw. 13, to appear. Foley, T.A. (1987b), Weighted bicubic spline interpolation to rapidly varying data, ACM Trans. Graphics 6, to appear. Gordon, W.J. and Riesenfeld, R.F. (1974), B-spline curves and surfaces, in: Barnhill, R.E. and Riesenfeld, 1LF., eds., Computer Aided Geometric Design, Academic Press, New York. Kochanek, D.H. and Bartels, R.H. (1984), Interpolating splines with local tension, continuity, and bias control, Computer Graphics 18, 33-41. Mortenson, M.E. (1985), Geometric Modeling, Wiley, New York. Nielson, G.M. (1974), Some pi¢cewise polynomial alternatives to splines under tension, in: Barnhin~ R.E. and Riescnfeld, R.F., eds., Computer Aided Geometric Design, Academic Press, New York. Salkauskas, K. (1984), C 1. splines for interpolation of rapidly varying data, Rocky Mountain J. Math. 14, 239-250. Schumaker, L.L. (1981), Spline Functions: Basic Theory, Wiley, New York. Varga, R.S. (1962), Matrix Iterative Analysis, Prentice-Hall, Englewood Cliffs.