- Email: [email protected]
Knot removal for parametric B-spline curves and surfaces
Computer Aided Geometric Design 4 (1987) 217-230 North-Holland
217
Knot removal for parametric B-spline curves and surfaces Tom L Y C H E and Knut M O R K E N
*
Institutt for informatikk, University of Oslo, Box 1080, Blindern, 0316 Oslo 3, Norway Presented at Mathematisches Forschungsinstitut Oberwolfach February 1987 Received April 1987; revised July 1987
Abstract. This paper is the third in a sequence of papers in which a knot removal strategy for splines, based on certain discrete norms, is developed. In the first paper, approximation methods defined as best approximations in these norms were discussed, while in the second paper a knot removal strategy for spline functions was developed. In this paper the knot removal strategy is extended to parametric spline curves and tensor product surfaces. The method has been implemented and thoroughly tested on a computer. We illustrate with several examples and applications.
Keywords. Splines, B-splines, curves, tensor product surfaces, discrete norms, knot removal.
1. I n t r o d u c t i o n
Piecewise polynomials have become an important tool in many automated design processes and a commonly used class of functions for the approximation of functions and data. As a basis for the representation of piecewise polynomial functions and parametric curves, the B-splines are convenient, and tensor products of B-splines provide a basis for surfaces which is easy to handle both mathematically and in a computer. In two previous papers, ,the authors described a technique for removing knots from a given spline without perturbing the spline more than a given tolerance [Lyche & M~rken '87a, '87b]. We will use the terms 'knot removal' and 'data reduction' interchangeably to denote this process. The purpose of the present paper is to show how this technique can be extended to parametric spline curves and surfaces. To our knowledge, little work has been published on the problem of knot removal. Related problems (dealing mainly with the removal of one knot) are considered in [Kjellander '83a, '83b; Cox & Jones '86; Farin, Rein, Sapidis & Worsey '87]. The paper consists of five sections. In Section 2 we summarize some properties of certain discrete norms which are used in our knot removal processes. In Section 3 we show how the knot removal strategy for spline functions can be extended to parametric spline curves, and illustrate this with two examples. In Section 4 the knot removal strategy is extended first to explicit tensor product spline surfaces, and then to parametric surfaces. Again this is illustrated with two examples. In the final section we list some of the many possible applications of data reduction for splines. The examples in Section 3 and 4 are all similar in that they employ knot removal in a general process to obtain approximations to curves and surfaces. We refer to this technique as * Temporarily at the Center for Industrial Research (SI), Box 124, Blindern, 0314 Oslo 3, Norway from April 1986 to August 1987. 0167-8396/87/$3.50 © 1987, Elsevier Science Publishers B.V. (North-Holland)
T. Lyche, K. Morken / Knot removal
218
approximation by data reduction. The idea behind approximation by data reduction is to start with an initial approximation which is easy to compute, and which is a good approximation because it is given by a large n u m b e r of parameters. T h e n we can perform knot removal to reduce the n u m b e r of parameters, but hopefully retain the quality of the approximation. In the remaining part of this section we introduce some basic notation. Let k and N be positive intergers, and t = (t~)N~ k be a sequence of real numbers with t i < t~+k for i = 1, 2 , . . . , N. We have N B-splines (Bi,k,t)~= N a, associated with the knot vector t, and these we assume to be normalized to sum to one. When no confusion is possible, we will often shorten B~,k,, to B~,~, or B,, or even Bi. The linear space of functions spanned by these B-splines we denote by Sk,,, and we refer to elements of Sk,, written in terms of this basis, as B-spline functions or simply as splines. If we define the vector function b, = ( B L t , . . . , Bu,t) T, a spline of f in Sk,, with coefficients c = ( q , . . . , CN)"r m a y be written as f = F,U=lc,B~,t = bT,e. A B-spline function is a piecewise polynomial of order k (degree k - 1), defined on all of IR, but it is the interval [tk, tu+~] that is of primary interest, since restricted to this interval, the space Sk,, contains all polynomials of order k or less. To avoid degenerate cases we assume in this paper that t k < tk+ ~ and tu < tN+l, a n d that N >/- k. Basic properties of splines and B-splines can be found in [de Boor '78] and [Schumaker '81].
2. Coefficient norms for B-spline curves and surfaces
A parametric B-spline curve f i n R a, of order k, with knots t, is a curve in R a where each c o m p o n e n t is a B-spline function with the same k n o t vector t, i.e., N
N
i=1
i=1
f=(i',...,i<') T: E (cT,...,
N
F ( u , v ) = Y'~ E
C,,jBi,k,,s(u)Bj,k2,t(v)=bs,k,(u)TCbt,k2(V),
i=1 j=l
where the coefficients have been grouped into an M × N matrix C = ( ( c i j ) i ~ l ) f = 1. The function F is called a tensor p r o d u c t B-spline surface. As for curves, we can extend this definition to parametric B-spline surfaces in R a by allowing the coefficients ci,j to be vectors in R a.
In our knOt removal process we try to approximate a given spline f by a spline ~ on a knot vector r = (ri)~+ak which is a subsequence of t, and hence S~,, __ ~k,,- We can therefore always express f and ~ in the same basis, e.g., the B-spline basis of Sk,,. The transformation which maps ~ from Sk, , tO ~k,t is represented by an N × n matrix Ak,,, , relative to the B-spline bases of Sk, , and Sk, ,. It will be convenient to shorten this to A,,, or A, or just A. Thus, if has the B-spline coefficients c = ( c l , . . . , c,) in Sk,,, it has the coefficients A c in Sk, ,. The matrix Ak,~, t is called the B-spline k n o t insertion matrix of order k from ~- to t, and its elements are called discrete B-splines. Discrete B-splines have many properties similar to those of the usual B-splines, and there exist stable and efficient algorithms to c o m p u t e them [Boehm '80; Cohen, Lyche & Riesenfeld '80; Lyche & M o r k e n '86].
219
T. Lyche, K. Msrken / Knot removal
The quality of an approximation is usually measured in a suitable norm. We will use a family of discrete norms which approximate the usual LP-norms well, at least for moderate values of k. Recall that the LP-norm of a function f is defined by l/P
1
(/ 1 lflP
IIf IILP =
9 forl
SUPIf I?
and that the ZP-norm of a vector (LE IR” is defined by
II~IIP= i
(Clailp)‘“, i mala,l, i
forlGp
co.
If f = l$c, the 1P, t-norm of f is defined by
II f II /p,f=
iI
l/P
E
1
lCilP(ri+k-fi)/k
>
i=l ma I ci 1,
forlGpC03, forp=
(2-l)
cc.
i
This may also be written as ]If II,p,,= dimension N with diagonal elements e, ,= I.1
((li+k-fi)/k)l’P,
IIE:/Pc/I,p, where &r/p is a diagonal matrix of
for 1
~0,
i 1, The use of these norms is motivated by the following fundamental inequalities, due to de Boor, which states the equivalence of the Zp, t-norms and the usual LP-norms, see, e.g., [de Boor ‘76, p. 161,
II IP,rQ II f II LP
Dk’llf
G
II f II 1p.t.
(2-4
Here D, is a constant which only depends on k. Since D, is of moderate size for small k, the IF, t-norm provides a good approximation to the LP-norm. If p is a knot vector with t c p, we also have [Lyche & Msrken ‘87b]
(2-3)
II f II fP,p G II f II IP,rFor parametric B-spline curves these inequalities apply to each component. If F = b:,$b~,kz is a tensor product B-spline surface, the inequalities (2.2) imply that
D;‘D$ II F II /p,s,r G II F II LP
G
IIF
Il ~p,s,t,
(2-4
where (IF IIIP,s,tis defined by
II F II 1P.s.r = (1~sl’pcE:‘p(IIP9 and IPII,p is the IP-norm of the M
x N matrix B, viewed as a vector in lWMN(for p = 2 this is the Frobenius norm). For a parametric tensor product B-spline surface the inequalities (2.4) apply in each component similarly to the case of parametric curves. In this paper we shah only use the ZP,t-norms for p = 2 and p = 00. We will use the 12, t-norm to obtain an approximation +* to f = bTc on a knot vector r with T G t. The to f in the 12,t-norm, i.e., as the approximation +* is defined as the best approximation solution of the minimization problem
*yz
.z
Ilf-~IlP,t.
220
T. Lyche, K. Morken / Knot removal
This leads to the least squares problem
rain II E,
TER~
(2.5)
T - ,=)11,=.
If the solution of (2.5) is y*, then
ep* =b~TT.. For m o r e details see [Lyche & Morken '87a,
'87b]. If we try to approximate a parametric B-spline curve with knots t by minimizing the l 2, t-norm of the error in each component, we have to solve problems like (2.5) for the B-spline coefficients of each component, but with E t and A,,t fixed. Similarly, we can approximate a surface F = bT,Cb, by a surface ~ * = boTF . b, with o _ s and , r c t, by minimizing the 12,s,t-norm of the error F - ~ * . This leads to the least squares problem
min
_F'~R rex.
Ey2(Ao,sFA~,t-C)EXt/2
,2,
(2.6)
which can be solved by first finding the solution X * e R m x N of x~Rm×N
leU(aosx-¢)ll,2,
and then finding the solution F * of rain
( F A T - X * ) E , a/2
I~E~ mxn
The matrix /'* is then also the solution of (2.6). Thus, the problem has been reduced to first obtaining best approximations to a set of N curves in the 12,s-norm resulting in a set of m curves which are then best approximated in the 12,t-norm. This is typical for certain types of approximation by tensor products, see [Greville '61, de Boor '79, Cox '86]. Again, the parametric case leads to similar problems for the coefficients of each component, but with A,.~, A,a, E,, and E t fixed. In C A G D applications, the 12,t-norm is n o t very convenient for measuring the error of an approximation since it does n o t provide good pointwise error estimates. We therefore measure the error of an approximation in the l ~, t-norm (l~,s, t-norm for surfaces) which we know to be a reasonable approximation to the L ~ - n o r m by (2.2) ((2.4) for surfaces). The l °~, t-norm is also used for ranking the knots, see the next section.
3. Knot removal for parametric B-spline curves In this section we indicate how the knot removal strategy for B-spline functions introduced in [Lyche & Morken '87b], can be extended to parametric B-spline curves. Suppose f = (fl,..., fd)T is a parametric B-spline curve in R a with k n o t vector t, and let tolerances ~'= ( q , . . . , cd) and a n o r m I1" II on Sk,, be given. The problem we would like to consider is to determine a shortest possible knot vector ~, with ~-c t, and a spline curve q~ = (q~l,..., qfl)T where each q,i is in $k, , such that II fi
_ qy II ~< Cg, for i = 1 , . . . , d.
(3.1)
We do not claim to solve this problem by finding the shortest possible knot vector • with these properties. However, we do believe that our heuristic approach is a reasonable alternative to finding a sphne with as few knots as possible. The first step towards solving the problem is to assign a weight wji to each c o m p o n e n t f i at each interior knot tj of t. These weights should be a measure of the significance of the knot tj in representing the spline ~ A natural approach is to let wji be the distance (as measured by
T. Lyche, K. Morken / Knot removal
221
s o m e n o r m ) between f i a n d the subspace Sj = Sk,,_(,j) o b t a i n e d by r e m o v i n g tj f r o m t. M o r e precisely, if tj is a simple k n o t of t we set wji = d i s t ( f f , S j ) =
rain g~S~
Ilff-glleo,,,
(3.2)
w h e r e II f ' - g II ~,,, is just the max n o r m of the B-spline coefficients of f ' - g written as a spline in Sk,t, cf. (2.1). T h e minimization p r o b l e m (3.2) has a local solution which can be f o u n d in O ( 5 k ) operations for a spline of order k. See [Lyche & M o r k e n '87b] for an algorithm. T h e c o m p u t a t i o n of the weights in the presence of multiple k n o t s can also be f o u n d there. T h e weights are used to determine in which order the k n o t s should be r e m o v e d f r o m t. O n e possibility w o u l d be to set wj = max~w} and remove the k n o t s in order of increasing wj's. This w o r k s well in some cases where all the c o m p o n e n t s of the tolerance vector are equal, b u t has s o m e obvious shortcomings. F o r example, it does n o t w o r k for a circle where all weights are a p p r o x i m a t e l y equal. Also, we believe that it is i m p o r t a n t to allow for different tolerances in the c o m p o n e n t s of the curve. We have instead used the following approach. T o each k n o t tj, we assign a ranking number
vj= maxvj, l
w h e r e vj is such that wji E r i 11). Here I ~i = [ 0 , (//2), 1~i =[(i/2, el), and in general 1]= [2J-2c/, 2 J - l q ) for j = 1, 2, 3, ... Let vjl ~< vj2~< --. ~< vj, be an ordering of the r a n k i n g n u m b e r s of the l interior knots of t in nondecreasing order. If we are to r e m o v e i knots f r o m t, t h e n all k n o t s with r a n k i n g n u m b e r s less than vj, will be removed. T o describe h o w to choose the remaining, say p , k n o t s to be removed, let u a ~< u 2 ~< • • • ~< u r be the k n o t s with r a n k i n g n u m b e r s equal to vs.,, listed in the order in which they occur in t. W e w o u l d then r e m o v e { Udl, Ud:,..., Ud, }, where dj = [(r + 1 ) ( j - ½)/p + ½J for j = 1, 2 , . . . , p (i.e., the n u m b e r ( r = 1 ) ( j - ½)/p is r o u n d e d to the nearest integer). In other words, the last p k n o t s are r e m o v e d a l m o s t u n i f o r m l y on subscripts. T h e k n o t removal process consists of three parts. W e call these Rank, Remove a n d Approximate. T h e c o m p u t a t i o n of weights and ranking n u m b e r s as j u s t described constitute the R a n k phase. In A p p r o x i m a t e , we use a discrete a p p r o x i m a t i o n m e t h o d suggested in [Lyche & M o r k e n '87a]. If the k n o t vector "r is given, with "r c t, the a p p r o x i m a t i o n ~,' t o f i from S , , , is given as the u n i q u e solution of the least squares p r o b l e m
rain Il f~-gllt~,,.
g~Sk.~
I n R e m o v e we use A p p r o x i m a t e on each c o m p o n e n t of f to d e t e r m i n e the m a x i m u m n u m b e r of k n o t s that can be r e m o v e d such that the corresponding spline a p p r o x i m a t i o n yields an error smaller t h a n the given tolerance. A reasonable a s s u m p t i o n is that the error in the a p p r o x i m a tion to f i n c r e a s e s with the n u m b e r of knots that are removed. T h e n u m b e r of k n o t s to r e m o v e c a n therefore be d e t e r m i n e d by a binary search. In other words, we can try to remove half the k n o t s and c o m p u t e the resulting a p p r o x i m a t i o n and its error. If the error is less t h a n the tolerance, we then try to remove 3 / 4 of the knots, if the error is greater t h a n the tolerance, we try to remove 1 / 4 of the knots a n d so on. R a n k , Remove, and A p p r o x i m a t e m a y be used as indicated above to d e t e r m i n e a k n o t vector ,ra with "q c t and an a p p r o x i m a t i o n ~ to f w i t h knots ~5, such that the error is less t h a n ~. If ,q has n o interior knots, we can be very satisfied because we have r e m o v e d all possible knots. I n the other extreme we have ~-~ = t implying that no k n o t s were removed. In general, the result will be s o m e t h i n g in between. We can then try to r e m o v e m o r e k n o t s as follows. R a n k the k n o t s ~'1 of ~1, a n d use R e m o v e to obtain a k n o t vectorz2 with "r2 c ,rI c t, b u t in A p p r o x i m a t e always c o m p u t e a p p r o x i m a t i o n s to the original spline f. In this way we o b t a i n an a p p r o x i m a -
222
T.
Lyche, K Morken / Knot removal
tion ~2 to f o n the k n o t vector ~'2- This process m a y clearly be c o n t i n u e d until n o m o r e knots can be removed. So far we have not said anything a b o u t the n o r m we use. In order to control the point-wise error our programs have used t h e / ° ° , t - n o r m to m e a s u r e the error in (3.1). It can be shown that the Hausdorff distance between f and ~ will then be b o u n d e d by vrd m a x i II f i _ ~, It +=,,- N o t e that the reduced curve q) inherits the parametrization of F In [Lyche & M o r k e n '87b] we indicated a general m e t h o d , based on d a t a reduction, for fitting a spline to univariate functions and data. Here, we show h o w this technique can be extended to parametric curves. Suppose we are given m points ( p--,j ) i,n = a in R d . Typically, we have d = 2 or d = 3. T o fit a cubic spline curve to these p o i n t s , we can proceed as follows. Let f~(t) be the piecewise linear inte~polant to the data, using for example, the a c c u m u l a t e d cord length p a r a m e t r i z a t i o n given bYfl(ti) =.ffi for i = 1, 2 , . . . , m, where
t, =
0,
for i = 1;
t,_l + II p-;-/,-111,=,
for i = 2, 3 , . . . , m.
Explicitly, we have
flt)=((ti--t)ffi1 +(t--ti-1)p-'i)/(ti--ti-1 ) in the interval [ti_ ~, t,]. By i n t r o d u c i n g one extra k n o t at each end, say t o = t 1 a n d t,,+~ = tm, we can represent f~ in terms of linear B-splines. G i v e n a tolerance vector ~ R a, we can then c o m p u t e a parametric, cubic spline a p p r o x i m a t i o n to f~ with error less t h a n ~ at each data point, in two stages: 1. A p p l y d a t a reduction to f~ with k = 2 a n d a tolerance ~ < ~" to obtain a n e w polygonal curve f2. 2. A p p l y data reduction to f~ written as a cubic spline with triple k n o t s at each break point, using a tolerance ~'2 = 2 - ~ (see, e.g., [Lyche & M o r k e n '87b] for details about the conversion of f~ to a cubic spline). This yields a cubic spline a p p r o x i m a t i o n f~ with k n o t vector 'r. N o t e that this m e t h o d m a y be used to obtain spline a p p r o x i m a t i o n s of any o r d e r k >t 2, by writing f~ as a spline of order k instead of as a cubic spline in the second step. Example 3.1. Suppose that we have the rn = 1001 d a t a p o i n t s in R 2 given by
xi= q(t~)
cos ti,
y~=
1
1
q(t~)
sin t~;
-
.
--1
for i = 1, 2 , . . . , rn;
.
.
.
i
Fig. 3.1. The parametric cubic spline approximation to a spiral using the approximation technique outlined in the text (~ = (0.001, 0.001)).
T. Lyche, K. Morken / Knot removal
223
0 .0010
OmOOOe 0-0006 Q o~
0-0004
l: B i
\
°'
\
0-0002 0-0000 --0-0002 --0-000~ --0-000~ --0-0008 --0.0010
I
-0
2 -0
3 -0
~: - 0
5 -0
6
-0
~' - 0
Fig. 3.2. The error in the cubic spline approximation to the spiral. The dashed curve shows the error in the x-component, while the solid curve shows the error in the y-component. The diamonds show the location of the interior knots.
where
5(x) 1(x)2
q(t)=l-~
~
+-~ ~
fort~[0,4~r],
and ti = 4¢r(i - 1 ) / ( m - 1) for i = 1, 2,..., m. The curve defined by (x(t), y(t)) = q(t)(cos t, sin t) is spiral like, and using the multistage technique with 2--(10 -3, 10-3), ~ = (10 -4, 10-4), and 22 = 2 - ~ , we obtained the cubic spline approximation shown in Fig. 3.1. Fig. 3.2 shows the piecewise linear interpolant to the two error components x ( t i ) - x i (dashed) and y ( t i ) - y i (solid), for i = 1 , 2 , . . . , m. The location of the knots is also shown. Note how the density of the knots increase with increasing curvature.
/oo - -. /
/
=
\
"\
• Ai
\
i
J --2
Fig. 3.3. Cubic spline approximation to the offset curve of sin(t) at a distance of 0.4.
T. Lyche, K. Merken / Knot removal
224 0.0010 0.0008 0.0006 0-000~ 0.0002 '
0-0000 --0-0002
t ,t
t t
;
•
:
,
t lt
•
•
.
tt m 0 • *
,
t ,
--O.0OO~
,
,
.
,".,
't
:
t°
t
,
'
,
"
%'
••
~t• • k'
V
1;
t
'
,,
tm
,
•
,,
--0-000~ --0°0008 --0.0010
i
-0
2 -0
3 -0
4 -0
5 -0
6 -0
Fig. 3.4. The error in the cubic spline a p p r o x i m a t i o n to the offset curve. T h e d a s h e d curve shows the error in the x - c o m p o n e n t , while the solid curve shows the error in the y - c o m p o n e n t .
Example 3.2. In this example we use the spline approximation method outlined above to model an offset curve. Recall that if (x(t), y(t)) is a plane curve, then the offset curve (xo(t), yo(t)) at a normal distance a is given by y' X°--X+
(Xt2+)y,2
x' 1/2a'
Y°=Y-
(xt2 +
1/£ a
where x' and y ' denote the derivatives of x(t) and y(t), respectively. Fig. 3.3 shows the curve (x(t), y(t)) = (t, sin ,rt) for t ~ [ - 0 . 5 , 2] (dashed), and a cubic spline approximation to the corresponding offset curve ~solid) using a = 0.4. In the offset curve there is an additional loop around t = 1/2 (not shown). To model this offset curve, we used the same method as above, with 2 = (10 -3, 10 -3) and ~ = (10 -4, 10-4). We sampled 1001 points on the offset curve uniformly in the parameter interval [ - 0 . 5 , 2], removed the loop, and supplemented the remaining set of points with the singular point, which we calculated as the intersection of the two relevant line segments. Fig. 3.4 shows the piecewise linear interpolants to the two error components at the sampled points. Note that the singularity presents no problem. A triple knot is left at the parameter value corresponding to this point. Periodic and closed curves present a problem for our knot removal process, since in general, such a curve would become open. One way out of this is to introduce constraints into the least squares approximation method so that endpoints can be kept fixed. This approach is not completely satisfactory however, since it would distinguish one arbitrary point (the beginning and end point) by keeping it fixed, while the other points are allowed to vary within the tolerance. A more satisfactory solution is obtained if we go back and slightly alter the knot removal strategy for explicit spline functions. The problem is that we only remove interior knots, which makes the beginning and end of a spline special. Suppose that we have k-tuple knots at the beginning and end of the spline (this is no restriction, since on the interval of interest, we can always convert to such a representation). Extending the coefficients and knots periodically, we can then also compute weights for the first k knots (which are equivalent to
T. Lyche, K. Morken / Knot removal
225
the last k knots). The removal of say i knots at the end can then be i m p l e m e n t e d as i constraints in the a p p r o x i m a t i o n routine requiring that the first i - 1 derivatives should be the same at the beginning and end of the spline.
4. Knot removal for tensor product B-spline surfaces One way to do k n o t removal for tensor product B-spline surfaces, is to p e r f o r m a parametric B-spline curve reduction in each p a r a m e t e r direction. Suppose that an explicit surface F with k n o t vectors s and t and coefficient matrix C ~ R M× N is given, together with a tolerance ~. To obtain a d a t a reduction along s, we simply leave out the second variable and consider the parametric curve f ~ = b T c in R N. For this curve we can do knot removal as indicated above, using the tolerance vector ~ = ( c / 2 . . . . , c / 2 ) ~ R N. This results in a new knot vector a, and a parametric curve ~o, with coefficients C ' given by ~,o = b T c '. We then have a surface F'=bTC'b, which deviates less t h a n ~/2 from F in the l°°,s,t-norm. To obtain a data reduction in the second p a r a m e t e r direction, we can proceed similarly. We p e r f o r m data reduction o n the parametric curve f~= C ' b , in R " using the tolerance vector ~'2 = ( c / 2 , . . . , % / 2 ) ~ R " . The result is a parametric curve with knot vector ,r and coefficients F, given by f , = Fb,. W e then have a reduced surface ~ = b T F b , with II F' - • II 1%o,, ~< c/2. For the total error II F - • II L~, we then find
II F -
q~ II ,~ ~ II F -
F ' II z~,s,, + II F ' - ~ II z~,s,,
IIF-F' II z~,~,, + II F ' -
~11 z~o o , , ~ c,
where we have used inequalities (2.2), (2.3) and (2.4). The advantage of this a p p r o a c h to data reduction for surfaces is that it requires very little p r o g r a m m i n g effort when the routines for data reduction for parametric curves are given. On the other hand, there is a lack of s y m m e t r y between the two p a r a m e t e r directions which might n o t be satisfactory. A n alternative m e t h o d which requires m o r e programming, but is aesthetically m o r e pleasing, is the following. Consider the p a r a m e t r i c curves fs = b T c and ft = CbT in R m and R N, respectively. U s i n g the tolerance vectors ~M = (c . . . . , c) ~ R M and ~iv = ( c , . . . , c) ~ R N, we obtain ranking n u m b e r s (/~j) a n d (vj) for the knots of these curves. Suppose that we are to remove i knots, h o w can we pick these from both s and t in a natural and symmetric way? Let rq ~< ~r2 ~< • • • ~< rrt be a n o n d e c r e a s i n g ordering of the ranking n u m b e r s (~ti) u (vj). We then remove all knots with ranking n u m b e r s less than ~'~. T h e question is how to remove the remaining, say p, knots. Let x~ <~ x a <<. • • • <~ Xrl and Yl <~Y2 <~ " " " Yr2 be the knots of s and t respectively, with ranking n u m b e r s equal to rrg, and listed in the order in which they occur in s and t. We can then pick r l p / ( r 1 + r2) knots uniformly on subscripts f r o m ( x j ) a n d r v p / ( r 1 + r2) knots uniformly on subscripts f r o m ( y j ) (see the previous section for a description of ' u n i f o r m l y o n subscripts'). In R e m o v e we can then search for the m a x i m u m n u m b e r of knots to remove as before, while in A p p r o x i m a t e we have to solve least squares problems of the sort indicated at the e n d of Section 2. We have n o t yet i m p l e m e n t e d the second method, so that in our examples we have only used the first method. The data fitting m e t h o d of the previous section can also be used to fit a surface to a set of points / ~ , . j = ( x , , yj, z,.j)
/i=1, ~j=l,
2,...,m; 2 .... ,n.
on the rectangular grid x 1 < x 2 < • • • < x,, and Yl < Y2 < " " " < Y,,- Let F 1 be the piecewise bilinear interpolant to the data, i.e., for all i and j we have F l ( x ~, yj) = z~.j, and F 1 is in each
226
T. Lyche, K. Morken / Knot removal
O
•O
O
Fig. 4.1. Bicubic spline approximation to a polynomial surface of degree 12. The interior knots are shown as diamonds.
rectangle [x,, Xi+l] X [yj, Yj+I] a bilinear function. Given a tolerance E, we can proceed as before. 1. Apply data reduction to F 1 written as a linear combination of bilinear B-splines using a tolerance c I < c. This produces a new tensor product surface F 2. 2. Write F 2 as bicubic tensor product B-spline surface with triple k n o t lines. Perform data reduction on F 2 using a tolerance c 2 = e - q . This results in a third tensor p r o d u c t B-spline surface F 3 which satisfies II Fa - F3 II L~ ~ ~Again we have restricted ourselves to the cubic case, but using the same strategy, we could obviously produce approximations of order k 1 in the first parameter direction and order k 2 in the second parameter direction, where kl and k 2 are arbitrary integers greater than 1. Example 4.1. The method outlined above may be used to c o m p u t e a p p r o x i m a t i o n s to a given surface. In this example we choose the polynomial surface F given by
F(x, y)= (1- 2)6(1-
y ) 6 + 103( 1 _ X ) 3 X 3 ( l _ _ y ) 3 y 3 + yr( 1 __ 2)6
+xr(1-- ~) 6 with (x, y ) ~ [0, 1] 2. To get the required data points we sampled this surface at a uniform 100 x 100 grid in [0, 1] 2. We used c = 0.005 and c 1 = c/10. After the first stage, the grid was reduced to a 72 x 72 grid. The final bicubic approximation is shown in Fig. 4.1, and as can be seen from the figure, it has 5 interior knots in each parameter direction. Fig. 4.2 shows the difference between F and the bicubic approximation sampled at a 30 x 30 grid in [0, 1] 2. This difference is less than 4.7 x 10 -3 in the L~-norm. Example 4.2. In this example we try to approximate the function
F(x, y ) = t a n h ( - 3 g ( x ,
y))+l
for(x, y)~[0,1]
2,
where g(x, y) = b(y - c) 2 - x + a and a = - 10, b = 0.595576, and c = - 3.79762. The plane curve defined by g(x, y ) = 0 is then a parabola which passes through the two points (0, 0.3) and (1, 0.5) with axis parallel to the x-axis. The surface F has a ridge in the n e i g h b o r h o o d of this parabola and is approximately equal to 0 and 2 on the two sides of this ridge. The function
T. Lyche, K. Morken / Knot removal
227
Fig. 4.2. The error in the bicubic approximation in Fig. 4.1.
Fig. 4.3. Bicubic spline approximation to a surface with a ridge along a parabola in the xy-plane.
Fig. 4.4. The error in the bicubic approximation in Fig. 4.3. was sampled at a regular 200 x 100 grid, a n d the a p p r o x i m a t i o n process was run with a tolerance of c = 0.01 a n d c 1 = ~/10. T h e linear data reduction p r o d u c e d a 36 x 59 grid. The final a p p r o x i m a t i o n is s h o w n in Fig. 4.3, and the error surface in Fig. 4.4 (the error is b o u n d e d b y 7.2 x 10-3). Both of the two m e t h o d s outlined at the beginning of this section for p e r f o r m i n g data reduction o n explicit surfaces, generalize naturally to parametric surfaces, the only vague point being h o w to o b t a i n the ranking n u m b e r s for the two knot vectors. If F = b~C'b t is a parametric
228
T. Lyche, K. Morken / Knot removal
surface in R d, we can consider curves of the type f ] = b f C which may be thought of as a collection of d parametric curves ( f ~ , . . . , f~) in R M. Given a tolerance vector g E R d, we can then to each interior knot s i of s compute ranking numbers (/x,)j=l j d for each of the d curves as described above. To obtain a single ranking number for s,, we then set/~i = maxj #,J just as we did for parametric curves. The only other difference from the case of one explicit surface is that in Approximate we have to compute approximations to all the d component surfaces. A parametric surface approximation will be acceptable if the error in each component is smaller than the corresponding component of the tolerance vector ~.
5. Possible applications The knot removal algorithm outlined in this paper has potential applications to data fitting, CAGD, and data compression. In this section we consider briefly some of these applications. Most of these applications are based on the multistage approximation technique outlined for curves just before Example 3.1, and just before Example 4.1 for surfaces. We refer to this technique as approximation by data reduction. The general theme is to use a simple approximation scheme with little restriction on the number of parameters used, to obtain an initial approximation, and then systematically remove knots from this. As an initial approximation we have used piecewise linear interpolation, which has the advantage that the final approximation might reproduce jumps in the first derivative. If smoothness can be assumed, one could also use points on the curve (surface) as control points for say a (bi-) cubic spline.
1. Adaptive curve approximation. Approximation by data reduction can be used to fit any reasonably behaved space curve (isolated singularities in tangent and curvature allowed) by a spline, to within a given tolerance, in any /P-norm. The number and location of knots are chosen automatically by the algorithm. The process is illustrated in Examples 3.1 and 3.2.
2. Data fitting. In [Lyche & Morken '87a, Fig. 5.1; Lyche & Morken '87b, Example 6.2], we have given examples of fitting a spline function to data points in the plane. As an initial spline approximation, we used the piecewise linear interpolant to the data. If the data have large and varying errors, one could instead use a local weighted least squares method of low degree and with many knots, to compute the initial approximation fl. Note that we compute the spline approximation in the knot removal strategy by a least squares method, but that we can control the error to within a given tolerance in any weighted /P-norm.
3. Degree reduction. Suppose one has a piecewise polynomial of fairly high degree and want to approximate it by a spline of lower degree, say three, with few knots, to within a given tolerance. Example 6.5 in [Lyche & Morken '87b] and Example 4.1 in this paper illustrate this process. Note again that the number and location of the knots are chosen automatically. Degree reduction is useful for closely approximating a polynomial geometry by a spline geometry.
4. Conversion of piecewise polynomials to B-spline form. Related to the previous application is the problem of converting a piecewise polynomial from one representation to another. Using data reduction as illustrated in Example 6.5 in [Lyche & Marken '87b] and Example 4.1 in this paper, one can in principle approximate any piecewise polynomial by a B-spline curve or surface with relatively few knots. The problem is indeed a special case of Application 1.
T. Lyche, K. Morken / Knot removal
229
5. Modeling of curves. Data reduction can be used to obtain an accurate spline representation of a digitized curve. This has many applications. We mention the inking problem [Plass & Stone '83]. It could also be useful for designing fonts in automatic type setting systems. 6. Modeling of surfaces. Data reduction can also be used to fit a surface to points given on the surface in a regular fashion [Schmitt & Barsky '86]. We do believe that we can model most curves satisfactory, but in dealing with surfaces there are several problems which have to be dealt with. One problem is related to the 'squaring effect' in storage requirements which one encounters when going from curves to surfaces. Recall that in our multistage approximation scheme one has to sample (scan) the surface at a number of points so that the piecewise linear interpolant approximates the surface to well within the tolerance. For complex shapes it might be difficult to sample the surface at enough points. A possible way out of this problem is to use a more accurate local approximation than linear interpolation. Thus one could use a bicubic Hermite patch or a rectangular Coons patch [Faux & Pratt '79]. Alternatively some local spline approximation method [Lyche & Schumaker '75] could be used for the initial approximation. Another problem is the parametrization of data in space given along curves in space. Any parametrization must apply simultaneously to many curves. We do believe however, that the techniques described in this paper are likely to lead to schemes for rapid modeling of surfaces. 7. Postprocessing a spline approximation. In many applications an approximation or interpolation scheme results in a spline with too many knots. This is the case with for instance the smoothing spline, see, e.g., [de Boor '78, p. 235], which has a knot at every data point, and spline interpolation at a large n u m b e r of points. The knot removal process is well suited for computing an approximation with fewer knots within a given tolerance. One example of this can be found in [Lyche & Morken '87b, Example 6.1], where knots are removed from a smoothing spline. In practice we have often observed that the reduced data seem to be smoother (have fewer 'wiggles') than the original curve (surface). 8. Computing spline approximations to intersection curves. Many algorithms for computing intersection curves between surfaces produce a (large) number of points which somehow have to be represented as a curve. Approximation by data reduction can be applied to this problem. 9. Modeling of offset curves and surfaces. The modeling of offset curves and surfaces has been considered by many, see, e.g. [Tiller & Hanson '84] and [Arnold '86]. Our approach would be to first sample the offset curve or surface at a number of points, and then compute a spline approximation using the techniques described in this paper (cf. Example 3.2.). Using the triple knot technique, singularities parallel to the axes present no problem. 10. Data compression. If a given set of data points can be represented as a linear B-spline curve or bilinear tensor product surface, one can apply the algorithms in this paper with k = 2 to obtain a polygon with fewer points. Note that the points remaining after the reduction will be changed, but only within the given tolerance. In general, it would be recommended to use these changed data values instead of the original values in the reduced set.
Acknowledgement We would like to thank Tor D o k k e n and Richard F. Riesenfeld for useful discussions and comments.
230
T. Lyche, K. Morken / Knot removal
References Arnold, R. (1986), Quadratische und kubische Offset-B6zierkurven, Dissertation, Universit~it Dortmund. de Boor, C. (1976), Splines as linear combinations of B-splines. A survey in: Lorentz, G.G., Chui, C.K., and Schumaker, L.L., eds., Approximation Theory 11, Academic Press, New York, 1-47. de Boor, C. (1978), A Practical Guide to Splines, Springer, New York. Boehm, W. (1980), Inserting new knots into B-spline curves, Computer Aided Design 12, 199-201. Cohen, E., Lyche, T., and Riesenfeld, R. (1980), Discrete B-spline and subdivision techniques in computer-aided geometric design and computer graphics, Computer Graphics and Image Processing 14, 87-111. Cox, M.G. and Jones, H.M. (1987), Shape-preserving spline approximation in the 11 norm, in: Mason, J.C. and Cox, M.G., eds., Algorithms for Approximation, Clarendon Press, Oxford, 115-129. Farin, G., Rein, G., Sapidis, N., Worsey, A.J. (1987), Fairing cubic B-spline curves, Computer Aided Geometric Design 4, 91-103. Faux, I.D. and Pratt, M.J. (1979), Computational Geometry for Design and Manufacture, Ellis Horwood, Chichester, UK. Greville, T.N.E. (1961), Note on fitting of functions of several independent variables, J. Soc. Ind. Appl. Math. 9, 109-115. Kjellander, J.A.P. (1983), Smoothing of cubic parametric splines, Computer Aided Design 15, 175-179. Kjellander, J.A.P. (1983), Smoothing of bicubic parametric surfaces, Computer Aided Design 15, 288-293. Lyche, T. and Morken, K. (1986), Making the Oslo Algorithm more efficient, SIAM J. Numer. Anal. 23, 663-675. Lyche, T. and M~rken, K. (1987a), A discrete approach to knot removal and degree reduction algorithms for splines, in: Mason, J.C., and Cox, M.G., eds., Algorithms for Approximation, Clarendon Press, Oxford, 67-82. Lyche,T. and Morken, K. (1987b), A data reduction strategy for splines, Research Report no. 107 February 1987, Institute of Informatics, University of Oslo, to appear in IMA J. Numer. Anal. Lyche, T. and Schumaker, L.L. (1975), Local spline approximation methods, J. Approximation Theory 15, 294-325. Plass, M. and Stone, M. (1983), Curve-fitting with piecewise parametric cubics, Computer Graphics 17, 229-239. Schmitt, F.J. and Barsky, B.A. (1986), An adaptive subdivision method for surface fitting from sampled data, in: SIGRAPH "86 Conference Proceedings, Computer Graphics, 20, 179-188. Schumaker, L.L. (1981), Spline Functions: Basic Theory, Wiley, New York. Tiller, W. and Hanson, E.G. (1984), Offsets of two-dimensional profiles, IEEE Computer Graphics and Application 4, 36-46.
217
Knot removal for parametric B-spline curves and surfaces Tom L Y C H E and Knut M O R K E N
*
Institutt for informatikk, University of Oslo, Box 1080, Blindern, 0316 Oslo 3, Norway Presented at Mathematisches Forschungsinstitut Oberwolfach February 1987 Received April 1987; revised July 1987
Abstract. This paper is the third in a sequence of papers in which a knot removal strategy for splines, based on certain discrete norms, is developed. In the first paper, approximation methods defined as best approximations in these norms were discussed, while in the second paper a knot removal strategy for spline functions was developed. In this paper the knot removal strategy is extended to parametric spline curves and tensor product surfaces. The method has been implemented and thoroughly tested on a computer. We illustrate with several examples and applications.
Keywords. Splines, B-splines, curves, tensor product surfaces, discrete norms, knot removal.
1. I n t r o d u c t i o n
Piecewise polynomials have become an important tool in many automated design processes and a commonly used class of functions for the approximation of functions and data. As a basis for the representation of piecewise polynomial functions and parametric curves, the B-splines are convenient, and tensor products of B-splines provide a basis for surfaces which is easy to handle both mathematically and in a computer. In two previous papers, ,the authors described a technique for removing knots from a given spline without perturbing the spline more than a given tolerance [Lyche & M~rken '87a, '87b]. We will use the terms 'knot removal' and 'data reduction' interchangeably to denote this process. The purpose of the present paper is to show how this technique can be extended to parametric spline curves and surfaces. To our knowledge, little work has been published on the problem of knot removal. Related problems (dealing mainly with the removal of one knot) are considered in [Kjellander '83a, '83b; Cox & Jones '86; Farin, Rein, Sapidis & Worsey '87]. The paper consists of five sections. In Section 2 we summarize some properties of certain discrete norms which are used in our knot removal processes. In Section 3 we show how the knot removal strategy for spline functions can be extended to parametric spline curves, and illustrate this with two examples. In Section 4 the knot removal strategy is extended first to explicit tensor product spline surfaces, and then to parametric surfaces. Again this is illustrated with two examples. In the final section we list some of the many possible applications of data reduction for splines. The examples in Section 3 and 4 are all similar in that they employ knot removal in a general process to obtain approximations to curves and surfaces. We refer to this technique as * Temporarily at the Center for Industrial Research (SI), Box 124, Blindern, 0314 Oslo 3, Norway from April 1986 to August 1987. 0167-8396/87/$3.50 © 1987, Elsevier Science Publishers B.V. (North-Holland)
T. Lyche, K. Morken / Knot removal
218
approximation by data reduction. The idea behind approximation by data reduction is to start with an initial approximation which is easy to compute, and which is a good approximation because it is given by a large n u m b e r of parameters. T h e n we can perform knot removal to reduce the n u m b e r of parameters, but hopefully retain the quality of the approximation. In the remaining part of this section we introduce some basic notation. Let k and N be positive intergers, and t = (t~)N~ k be a sequence of real numbers with t i < t~+k for i = 1, 2 , . . . , N. We have N B-splines (Bi,k,t)~= N a, associated with the knot vector t, and these we assume to be normalized to sum to one. When no confusion is possible, we will often shorten B~,k,, to B~,~, or B,, or even Bi. The linear space of functions spanned by these B-splines we denote by Sk,,, and we refer to elements of Sk,, written in terms of this basis, as B-spline functions or simply as splines. If we define the vector function b, = ( B L t , . . . , Bu,t) T, a spline of f in Sk,, with coefficients c = ( q , . . . , CN)"r m a y be written as f = F,U=lc,B~,t = bT,e. A B-spline function is a piecewise polynomial of order k (degree k - 1), defined on all of IR, but it is the interval [tk, tu+~] that is of primary interest, since restricted to this interval, the space Sk,, contains all polynomials of order k or less. To avoid degenerate cases we assume in this paper that t k < tk+ ~ and tu < tN+l, a n d that N >/- k. Basic properties of splines and B-splines can be found in [de Boor '78] and [Schumaker '81].
2. Coefficient norms for B-spline curves and surfaces
A parametric B-spline curve f i n R a, of order k, with knots t, is a curve in R a where each c o m p o n e n t is a B-spline function with the same k n o t vector t, i.e., N
N
i=1
i=1
f=(i',...,i<') T: E (cT,...,
N
F ( u , v ) = Y'~ E
C,,jBi,k,,s(u)Bj,k2,t(v)=bs,k,(u)TCbt,k2(V),
i=1 j=l
where the coefficients have been grouped into an M × N matrix C = ( ( c i j ) i ~ l ) f = 1. The function F is called a tensor p r o d u c t B-spline surface. As for curves, we can extend this definition to parametric B-spline surfaces in R a by allowing the coefficients ci,j to be vectors in R a.
In our knOt removal process we try to approximate a given spline f by a spline ~ on a knot vector r = (ri)~+ak which is a subsequence of t, and hence S~,, __ ~k,,- We can therefore always express f and ~ in the same basis, e.g., the B-spline basis of Sk,,. The transformation which maps ~ from Sk, , tO ~k,t is represented by an N × n matrix Ak,,, , relative to the B-spline bases of Sk, , and Sk, ,. It will be convenient to shorten this to A,,, or A, or just A. Thus, if has the B-spline coefficients c = ( c l , . . . , c,) in Sk,,, it has the coefficients A c in Sk, ,. The matrix Ak,~, t is called the B-spline k n o t insertion matrix of order k from ~- to t, and its elements are called discrete B-splines. Discrete B-splines have many properties similar to those of the usual B-splines, and there exist stable and efficient algorithms to c o m p u t e them [Boehm '80; Cohen, Lyche & Riesenfeld '80; Lyche & M o r k e n '86].
219
T. Lyche, K. Msrken / Knot removal
The quality of an approximation is usually measured in a suitable norm. We will use a family of discrete norms which approximate the usual LP-norms well, at least for moderate values of k. Recall that the LP-norm of a function f is defined by l/P
1
(/ 1 lflP
IIf IILP =
9 forl
SUPIf I?
and that the ZP-norm of a vector (LE IR” is defined by
II~IIP= i
(Clailp)‘“, i mala,l, i
forlGp
co.
If f = l$c, the 1P, t-norm of f is defined by
II f II /p,f=
iI
l/P
E
1
lCilP(ri+k-fi)/k
>
i=l ma I ci 1,
forlGpC03, forp=
(2-l)
cc.
i
This may also be written as ]If II,p,,= dimension N with diagonal elements e, ,= I.1
((li+k-fi)/k)l’P,
IIE:/Pc/I,p, where &r/p is a diagonal matrix of
for 1
~0,
i 1, The use of these norms is motivated by the following fundamental inequalities, due to de Boor, which states the equivalence of the Zp, t-norms and the usual LP-norms, see, e.g., [de Boor ‘76, p. 161,
II IP,rQ II f II LP
Dk’llf
G
II f II 1p.t.
(2-4
Here D, is a constant which only depends on k. Since D, is of moderate size for small k, the IF, t-norm provides a good approximation to the LP-norm. If p is a knot vector with t c p, we also have [Lyche & Msrken ‘87b]
(2-3)
II f II fP,p G II f II IP,rFor parametric B-spline curves these inequalities apply to each component. If F = b:,$b~,kz is a tensor product B-spline surface, the inequalities (2.2) imply that
D;‘D$ II F II /p,s,r G II F II LP
G
IIF
Il ~p,s,t,
(2-4
where (IF IIIP,s,tis defined by
II F II 1P.s.r = (1~sl’pcE:‘p(IIP9 and IPII,p is the IP-norm of the M
x N matrix B, viewed as a vector in lWMN(for p = 2 this is the Frobenius norm). For a parametric tensor product B-spline surface the inequalities (2.4) apply in each component similarly to the case of parametric curves. In this paper we shah only use the ZP,t-norms for p = 2 and p = 00. We will use the 12, t-norm to obtain an approximation +* to f = bTc on a knot vector r with T G t. The to f in the 12,t-norm, i.e., as the approximation +* is defined as the best approximation solution of the minimization problem
*yz
.z
Ilf-~IlP,t.
220
T. Lyche, K. Morken / Knot removal
This leads to the least squares problem
rain II E,
TER~
(2.5)
T - ,=)11,=.
If the solution of (2.5) is y*, then
ep* =b~TT.. For m o r e details see [Lyche & Morken '87a,
'87b]. If we try to approximate a parametric B-spline curve with knots t by minimizing the l 2, t-norm of the error in each component, we have to solve problems like (2.5) for the B-spline coefficients of each component, but with E t and A,,t fixed. Similarly, we can approximate a surface F = bT,Cb, by a surface ~ * = boTF . b, with o _ s and , r c t, by minimizing the 12,s,t-norm of the error F - ~ * . This leads to the least squares problem
min
_F'~R rex.
Ey2(Ao,sFA~,t-C)EXt/2
,2,
(2.6)
which can be solved by first finding the solution X * e R m x N of x~Rm×N
leU(aosx-¢)ll,2,
and then finding the solution F * of rain
( F A T - X * ) E , a/2
I~E~ mxn
The matrix /'* is then also the solution of (2.6). Thus, the problem has been reduced to first obtaining best approximations to a set of N curves in the 12,s-norm resulting in a set of m curves which are then best approximated in the 12,t-norm. This is typical for certain types of approximation by tensor products, see [Greville '61, de Boor '79, Cox '86]. Again, the parametric case leads to similar problems for the coefficients of each component, but with A,.~, A,a, E,, and E t fixed. In C A G D applications, the 12,t-norm is n o t very convenient for measuring the error of an approximation since it does n o t provide good pointwise error estimates. We therefore measure the error of an approximation in the l ~, t-norm (l~,s, t-norm for surfaces) which we know to be a reasonable approximation to the L ~ - n o r m by (2.2) ((2.4) for surfaces). The l °~, t-norm is also used for ranking the knots, see the next section.
3. Knot removal for parametric B-spline curves In this section we indicate how the knot removal strategy for B-spline functions introduced in [Lyche & Morken '87b], can be extended to parametric B-spline curves. Suppose f = (fl,..., fd)T is a parametric B-spline curve in R a with k n o t vector t, and let tolerances ~'= ( q , . . . , cd) and a n o r m I1" II on Sk,, be given. The problem we would like to consider is to determine a shortest possible knot vector ~, with ~-c t, and a spline curve q~ = (q~l,..., qfl)T where each q,i is in $k, , such that II fi
_ qy II ~< Cg, for i = 1 , . . . , d.
(3.1)
We do not claim to solve this problem by finding the shortest possible knot vector • with these properties. However, we do believe that our heuristic approach is a reasonable alternative to finding a sphne with as few knots as possible. The first step towards solving the problem is to assign a weight wji to each c o m p o n e n t f i at each interior knot tj of t. These weights should be a measure of the significance of the knot tj in representing the spline ~ A natural approach is to let wji be the distance (as measured by
T. Lyche, K. Morken / Knot removal
221
s o m e n o r m ) between f i a n d the subspace Sj = Sk,,_(,j) o b t a i n e d by r e m o v i n g tj f r o m t. M o r e precisely, if tj is a simple k n o t of t we set wji = d i s t ( f f , S j ) =
rain g~S~
Ilff-glleo,,,
(3.2)
w h e r e II f ' - g II ~,,, is just the max n o r m of the B-spline coefficients of f ' - g written as a spline in Sk,t, cf. (2.1). T h e minimization p r o b l e m (3.2) has a local solution which can be f o u n d in O ( 5 k ) operations for a spline of order k. See [Lyche & M o r k e n '87b] for an algorithm. T h e c o m p u t a t i o n of the weights in the presence of multiple k n o t s can also be f o u n d there. T h e weights are used to determine in which order the k n o t s should be r e m o v e d f r o m t. O n e possibility w o u l d be to set wj = max~w} and remove the k n o t s in order of increasing wj's. This w o r k s well in some cases where all the c o m p o n e n t s of the tolerance vector are equal, b u t has s o m e obvious shortcomings. F o r example, it does n o t w o r k for a circle where all weights are a p p r o x i m a t e l y equal. Also, we believe that it is i m p o r t a n t to allow for different tolerances in the c o m p o n e n t s of the curve. We have instead used the following approach. T o each k n o t tj, we assign a ranking number
vj= maxvj, l
w h e r e vj is such that wji E r i 11). Here I ~i = [ 0 , (//2), 1~i =[(i/2, el), and in general 1]= [2J-2c/, 2 J - l q ) for j = 1, 2, 3, ... Let vjl ~< vj2~< --. ~< vj, be an ordering of the r a n k i n g n u m b e r s of the l interior knots of t in nondecreasing order. If we are to r e m o v e i knots f r o m t, t h e n all k n o t s with r a n k i n g n u m b e r s less than vj, will be removed. T o describe h o w to choose the remaining, say p , k n o t s to be removed, let u a ~< u 2 ~< • • • ~< u r be the k n o t s with r a n k i n g n u m b e r s equal to vs.,, listed in the order in which they occur in t. W e w o u l d then r e m o v e { Udl, Ud:,..., Ud, }, where dj = [(r + 1 ) ( j - ½)/p + ½J for j = 1, 2 , . . . , p (i.e., the n u m b e r ( r = 1 ) ( j - ½)/p is r o u n d e d to the nearest integer). In other words, the last p k n o t s are r e m o v e d a l m o s t u n i f o r m l y on subscripts. T h e k n o t removal process consists of three parts. W e call these Rank, Remove a n d Approximate. T h e c o m p u t a t i o n of weights and ranking n u m b e r s as j u s t described constitute the R a n k phase. In A p p r o x i m a t e , we use a discrete a p p r o x i m a t i o n m e t h o d suggested in [Lyche & M o r k e n '87a]. If the k n o t vector "r is given, with "r c t, the a p p r o x i m a t i o n ~,' t o f i from S , , , is given as the u n i q u e solution of the least squares p r o b l e m
rain Il f~-gllt~,,.
g~Sk.~
I n R e m o v e we use A p p r o x i m a t e on each c o m p o n e n t of f to d e t e r m i n e the m a x i m u m n u m b e r of k n o t s that can be r e m o v e d such that the corresponding spline a p p r o x i m a t i o n yields an error smaller t h a n the given tolerance. A reasonable a s s u m p t i o n is that the error in the a p p r o x i m a tion to f i n c r e a s e s with the n u m b e r of knots that are removed. T h e n u m b e r of k n o t s to r e m o v e c a n therefore be d e t e r m i n e d by a binary search. In other words, we can try to remove half the k n o t s and c o m p u t e the resulting a p p r o x i m a t i o n and its error. If the error is less t h a n the tolerance, we then try to remove 3 / 4 of the knots, if the error is greater t h a n the tolerance, we try to remove 1 / 4 of the knots a n d so on. R a n k , Remove, and A p p r o x i m a t e m a y be used as indicated above to d e t e r m i n e a k n o t vector ,ra with "q c t and an a p p r o x i m a t i o n ~ to f w i t h knots ~5, such that the error is less t h a n ~. If ,q has n o interior knots, we can be very satisfied because we have r e m o v e d all possible knots. I n the other extreme we have ~-~ = t implying that no k n o t s were removed. In general, the result will be s o m e t h i n g in between. We can then try to r e m o v e m o r e k n o t s as follows. R a n k the k n o t s ~'1 of ~1, a n d use R e m o v e to obtain a k n o t vectorz2 with "r2 c ,rI c t, b u t in A p p r o x i m a t e always c o m p u t e a p p r o x i m a t i o n s to the original spline f. In this way we o b t a i n an a p p r o x i m a -
222
T.
Lyche, K Morken / Knot removal
tion ~2 to f o n the k n o t vector ~'2- This process m a y clearly be c o n t i n u e d until n o m o r e knots can be removed. So far we have not said anything a b o u t the n o r m we use. In order to control the point-wise error our programs have used t h e / ° ° , t - n o r m to m e a s u r e the error in (3.1). It can be shown that the Hausdorff distance between f and ~ will then be b o u n d e d by vrd m a x i II f i _ ~, It +=,,- N o t e that the reduced curve q) inherits the parametrization of F In [Lyche & M o r k e n '87b] we indicated a general m e t h o d , based on d a t a reduction, for fitting a spline to univariate functions and data. Here, we show h o w this technique can be extended to parametric curves. Suppose we are given m points ( p--,j ) i,n = a in R d . Typically, we have d = 2 or d = 3. T o fit a cubic spline curve to these p o i n t s , we can proceed as follows. Let f~(t) be the piecewise linear inte~polant to the data, using for example, the a c c u m u l a t e d cord length p a r a m e t r i z a t i o n given bYfl(ti) =.ffi for i = 1, 2 , . . . , m, where
t, =
0,
for i = 1;
t,_l + II p-;-/,-111,=,
for i = 2, 3 , . . . , m.
Explicitly, we have
flt)=((ti--t)ffi1 +(t--ti-1)p-'i)/(ti--ti-1 ) in the interval [ti_ ~, t,]. By i n t r o d u c i n g one extra k n o t at each end, say t o = t 1 a n d t,,+~ = tm, we can represent f~ in terms of linear B-splines. G i v e n a tolerance vector ~ R a, we can then c o m p u t e a parametric, cubic spline a p p r o x i m a t i o n to f~ with error less t h a n ~ at each data point, in two stages: 1. A p p l y d a t a reduction to f~ with k = 2 a n d a tolerance ~ < ~" to obtain a n e w polygonal curve f2. 2. A p p l y data reduction to f~ written as a cubic spline with triple k n o t s at each break point, using a tolerance ~'2 = 2 - ~ (see, e.g., [Lyche & M o r k e n '87b] for details about the conversion of f~ to a cubic spline). This yields a cubic spline a p p r o x i m a t i o n f~ with k n o t vector 'r. N o t e that this m e t h o d m a y be used to obtain spline a p p r o x i m a t i o n s of any o r d e r k >t 2, by writing f~ as a spline of order k instead of as a cubic spline in the second step. Example 3.1. Suppose that we have the rn = 1001 d a t a p o i n t s in R 2 given by
xi= q(t~)
cos ti,
y~=
1
1
q(t~)
sin t~;
-
.
--1
for i = 1, 2 , . . . , rn;
.
.
.
i
Fig. 3.1. The parametric cubic spline approximation to a spiral using the approximation technique outlined in the text (~ = (0.001, 0.001)).
T. Lyche, K. Morken / Knot removal
223
0 .0010
OmOOOe 0-0006 Q o~
0-0004
l: B i
\
°'
\
0-0002 0-0000 --0-0002 --0-000~ --0-000~ --0-0008 --0.0010
I
-0
2 -0
3 -0
~: - 0
5 -0
6
-0
~' - 0
Fig. 3.2. The error in the cubic spline approximation to the spiral. The dashed curve shows the error in the x-component, while the solid curve shows the error in the y-component. The diamonds show the location of the interior knots.
where
5(x) 1(x)2
q(t)=l-~
~
+-~ ~
fort~[0,4~r],
and ti = 4¢r(i - 1 ) / ( m - 1) for i = 1, 2,..., m. The curve defined by (x(t), y(t)) = q(t)(cos t, sin t) is spiral like, and using the multistage technique with 2--(10 -3, 10-3), ~ = (10 -4, 10-4), and 22 = 2 - ~ , we obtained the cubic spline approximation shown in Fig. 3.1. Fig. 3.2 shows the piecewise linear interpolant to the two error components x ( t i ) - x i (dashed) and y ( t i ) - y i (solid), for i = 1 , 2 , . . . , m. The location of the knots is also shown. Note how the density of the knots increase with increasing curvature.
/oo - -. /
/
=
\
"\
• Ai
\
i
J --2
Fig. 3.3. Cubic spline approximation to the offset curve of sin(t) at a distance of 0.4.
T. Lyche, K. Merken / Knot removal
224 0.0010 0.0008 0.0006 0-000~ 0.0002 '
0-0000 --0-0002
t ,t
t t
;
•
:
,
t lt
•
•
.
tt m 0 • *
,
t ,
--O.0OO~
,
,
.
,".,
't
:
t°
t
,
'
,
"
%'
••
~t• • k'
V
1;
t
'
,,
tm
,
•
,,
--0-000~ --0°0008 --0.0010
i
-0
2 -0
3 -0
4 -0
5 -0
6 -0
Fig. 3.4. The error in the cubic spline a p p r o x i m a t i o n to the offset curve. T h e d a s h e d curve shows the error in the x - c o m p o n e n t , while the solid curve shows the error in the y - c o m p o n e n t .
Example 3.2. In this example we use the spline approximation method outlined above to model an offset curve. Recall that if (x(t), y(t)) is a plane curve, then the offset curve (xo(t), yo(t)) at a normal distance a is given by y' X°--X+
(Xt2+)y,2
x' 1/2a'
Y°=Y-
(xt2 +
1/£ a
where x' and y ' denote the derivatives of x(t) and y(t), respectively. Fig. 3.3 shows the curve (x(t), y(t)) = (t, sin ,rt) for t ~ [ - 0 . 5 , 2] (dashed), and a cubic spline approximation to the corresponding offset curve ~solid) using a = 0.4. In the offset curve there is an additional loop around t = 1/2 (not shown). To model this offset curve, we used the same method as above, with 2 = (10 -3, 10 -3) and ~ = (10 -4, 10-4). We sampled 1001 points on the offset curve uniformly in the parameter interval [ - 0 . 5 , 2], removed the loop, and supplemented the remaining set of points with the singular point, which we calculated as the intersection of the two relevant line segments. Fig. 3.4 shows the piecewise linear interpolants to the two error components at the sampled points. Note that the singularity presents no problem. A triple knot is left at the parameter value corresponding to this point. Periodic and closed curves present a problem for our knot removal process, since in general, such a curve would become open. One way out of this is to introduce constraints into the least squares approximation method so that endpoints can be kept fixed. This approach is not completely satisfactory however, since it would distinguish one arbitrary point (the beginning and end point) by keeping it fixed, while the other points are allowed to vary within the tolerance. A more satisfactory solution is obtained if we go back and slightly alter the knot removal strategy for explicit spline functions. The problem is that we only remove interior knots, which makes the beginning and end of a spline special. Suppose that we have k-tuple knots at the beginning and end of the spline (this is no restriction, since on the interval of interest, we can always convert to such a representation). Extending the coefficients and knots periodically, we can then also compute weights for the first k knots (which are equivalent to
T. Lyche, K. Morken / Knot removal
225
the last k knots). The removal of say i knots at the end can then be i m p l e m e n t e d as i constraints in the a p p r o x i m a t i o n routine requiring that the first i - 1 derivatives should be the same at the beginning and end of the spline.
4. Knot removal for tensor product B-spline surfaces One way to do k n o t removal for tensor product B-spline surfaces, is to p e r f o r m a parametric B-spline curve reduction in each p a r a m e t e r direction. Suppose that an explicit surface F with k n o t vectors s and t and coefficient matrix C ~ R M× N is given, together with a tolerance ~. To obtain a d a t a reduction along s, we simply leave out the second variable and consider the parametric curve f ~ = b T c in R N. For this curve we can do knot removal as indicated above, using the tolerance vector ~ = ( c / 2 . . . . , c / 2 ) ~ R N. This results in a new knot vector a, and a parametric curve ~o, with coefficients C ' given by ~,o = b T c '. We then have a surface F'=bTC'b, which deviates less t h a n ~/2 from F in the l°°,s,t-norm. To obtain a data reduction in the second p a r a m e t e r direction, we can proceed similarly. We p e r f o r m data reduction o n the parametric curve f~= C ' b , in R " using the tolerance vector ~'2 = ( c / 2 , . . . , % / 2 ) ~ R " . The result is a parametric curve with knot vector ,r and coefficients F, given by f , = Fb,. W e then have a reduced surface ~ = b T F b , with II F' - • II 1%o,, ~< c/2. For the total error II F - • II L~, we then find
II F -
q~ II ,~ ~ II F -
F ' II z~,s,, + II F ' - ~ II z~,s,,
IIF-F' II z~,~,, + II F ' -
~11 z~o o , , ~ c,
where we have used inequalities (2.2), (2.3) and (2.4). The advantage of this a p p r o a c h to data reduction for surfaces is that it requires very little p r o g r a m m i n g effort when the routines for data reduction for parametric curves are given. On the other hand, there is a lack of s y m m e t r y between the two p a r a m e t e r directions which might n o t be satisfactory. A n alternative m e t h o d which requires m o r e programming, but is aesthetically m o r e pleasing, is the following. Consider the p a r a m e t r i c curves fs = b T c and ft = CbT in R m and R N, respectively. U s i n g the tolerance vectors ~M = (c . . . . , c) ~ R M and ~iv = ( c , . . . , c) ~ R N, we obtain ranking n u m b e r s (/~j) a n d (vj) for the knots of these curves. Suppose that we are to remove i knots, h o w can we pick these from both s and t in a natural and symmetric way? Let rq ~< ~r2 ~< • • • ~< rrt be a n o n d e c r e a s i n g ordering of the ranking n u m b e r s (~ti) u (vj). We then remove all knots with ranking n u m b e r s less than ~'~. T h e question is how to remove the remaining, say p, knots. Let x~ <~ x a <<. • • • <~ Xrl and Yl <~Y2 <~ " " " Yr2 be the knots of s and t respectively, with ranking n u m b e r s equal to rrg, and listed in the order in which they occur in s and t. We can then pick r l p / ( r 1 + r2) knots uniformly on subscripts f r o m ( x j ) a n d r v p / ( r 1 + r2) knots uniformly on subscripts f r o m ( y j ) (see the previous section for a description of ' u n i f o r m l y o n subscripts'). In R e m o v e we can then search for the m a x i m u m n u m b e r of knots to remove as before, while in A p p r o x i m a t e we have to solve least squares problems of the sort indicated at the e n d of Section 2. We have n o t yet i m p l e m e n t e d the second method, so that in our examples we have only used the first method. The data fitting m e t h o d of the previous section can also be used to fit a surface to a set of points / ~ , . j = ( x , , yj, z,.j)
/i=1, ~j=l,
2,...,m; 2 .... ,n.
on the rectangular grid x 1 < x 2 < • • • < x,, and Yl < Y2 < " " " < Y,,- Let F 1 be the piecewise bilinear interpolant to the data, i.e., for all i and j we have F l ( x ~, yj) = z~.j, and F 1 is in each
226
T. Lyche, K. Morken / Knot removal
O
•O
O
Fig. 4.1. Bicubic spline approximation to a polynomial surface of degree 12. The interior knots are shown as diamonds.
rectangle [x,, Xi+l] X [yj, Yj+I] a bilinear function. Given a tolerance E, we can proceed as before. 1. Apply data reduction to F 1 written as a linear combination of bilinear B-splines using a tolerance c I < c. This produces a new tensor product surface F 2. 2. Write F 2 as bicubic tensor product B-spline surface with triple k n o t lines. Perform data reduction on F 2 using a tolerance c 2 = e - q . This results in a third tensor p r o d u c t B-spline surface F 3 which satisfies II Fa - F3 II L~ ~ ~Again we have restricted ourselves to the cubic case, but using the same strategy, we could obviously produce approximations of order k 1 in the first parameter direction and order k 2 in the second parameter direction, where kl and k 2 are arbitrary integers greater than 1. Example 4.1. The method outlined above may be used to c o m p u t e a p p r o x i m a t i o n s to a given surface. In this example we choose the polynomial surface F given by
F(x, y)= (1- 2)6(1-
y ) 6 + 103( 1 _ X ) 3 X 3 ( l _ _ y ) 3 y 3 + yr( 1 __ 2)6
+xr(1-- ~) 6 with (x, y ) ~ [0, 1] 2. To get the required data points we sampled this surface at a uniform 100 x 100 grid in [0, 1] 2. We used c = 0.005 and c 1 = c/10. After the first stage, the grid was reduced to a 72 x 72 grid. The final bicubic approximation is shown in Fig. 4.1, and as can be seen from the figure, it has 5 interior knots in each parameter direction. Fig. 4.2 shows the difference between F and the bicubic approximation sampled at a 30 x 30 grid in [0, 1] 2. This difference is less than 4.7 x 10 -3 in the L~-norm. Example 4.2. In this example we try to approximate the function
F(x, y ) = t a n h ( - 3 g ( x ,
y))+l
for(x, y)~[0,1]
2,
where g(x, y) = b(y - c) 2 - x + a and a = - 10, b = 0.595576, and c = - 3.79762. The plane curve defined by g(x, y ) = 0 is then a parabola which passes through the two points (0, 0.3) and (1, 0.5) with axis parallel to the x-axis. The surface F has a ridge in the n e i g h b o r h o o d of this parabola and is approximately equal to 0 and 2 on the two sides of this ridge. The function
T. Lyche, K. Morken / Knot removal
227
Fig. 4.2. The error in the bicubic approximation in Fig. 4.1.
Fig. 4.3. Bicubic spline approximation to a surface with a ridge along a parabola in the xy-plane.
Fig. 4.4. The error in the bicubic approximation in Fig. 4.3. was sampled at a regular 200 x 100 grid, a n d the a p p r o x i m a t i o n process was run with a tolerance of c = 0.01 a n d c 1 = ~/10. T h e linear data reduction p r o d u c e d a 36 x 59 grid. The final a p p r o x i m a t i o n is s h o w n in Fig. 4.3, and the error surface in Fig. 4.4 (the error is b o u n d e d b y 7.2 x 10-3). Both of the two m e t h o d s outlined at the beginning of this section for p e r f o r m i n g data reduction o n explicit surfaces, generalize naturally to parametric surfaces, the only vague point being h o w to o b t a i n the ranking n u m b e r s for the two knot vectors. If F = b~C'b t is a parametric
228
T. Lyche, K. Morken / Knot removal
surface in R d, we can consider curves of the type f ] = b f C which may be thought of as a collection of d parametric curves ( f ~ , . . . , f~) in R M. Given a tolerance vector g E R d, we can then to each interior knot s i of s compute ranking numbers (/x,)j=l j d for each of the d curves as described above. To obtain a single ranking number for s,, we then set/~i = maxj #,J just as we did for parametric curves. The only other difference from the case of one explicit surface is that in Approximate we have to compute approximations to all the d component surfaces. A parametric surface approximation will be acceptable if the error in each component is smaller than the corresponding component of the tolerance vector ~.
5. Possible applications The knot removal algorithm outlined in this paper has potential applications to data fitting, CAGD, and data compression. In this section we consider briefly some of these applications. Most of these applications are based on the multistage approximation technique outlined for curves just before Example 3.1, and just before Example 4.1 for surfaces. We refer to this technique as approximation by data reduction. The general theme is to use a simple approximation scheme with little restriction on the number of parameters used, to obtain an initial approximation, and then systematically remove knots from this. As an initial approximation we have used piecewise linear interpolation, which has the advantage that the final approximation might reproduce jumps in the first derivative. If smoothness can be assumed, one could also use points on the curve (surface) as control points for say a (bi-) cubic spline.
1. Adaptive curve approximation. Approximation by data reduction can be used to fit any reasonably behaved space curve (isolated singularities in tangent and curvature allowed) by a spline, to within a given tolerance, in any /P-norm. The number and location of knots are chosen automatically by the algorithm. The process is illustrated in Examples 3.1 and 3.2.
2. Data fitting. In [Lyche & Morken '87a, Fig. 5.1; Lyche & Morken '87b, Example 6.2], we have given examples of fitting a spline function to data points in the plane. As an initial spline approximation, we used the piecewise linear interpolant to the data. If the data have large and varying errors, one could instead use a local weighted least squares method of low degree and with many knots, to compute the initial approximation fl. Note that we compute the spline approximation in the knot removal strategy by a least squares method, but that we can control the error to within a given tolerance in any weighted /P-norm.
3. Degree reduction. Suppose one has a piecewise polynomial of fairly high degree and want to approximate it by a spline of lower degree, say three, with few knots, to within a given tolerance. Example 6.5 in [Lyche & Morken '87b] and Example 4.1 in this paper illustrate this process. Note again that the number and location of the knots are chosen automatically. Degree reduction is useful for closely approximating a polynomial geometry by a spline geometry.
4. Conversion of piecewise polynomials to B-spline form. Related to the previous application is the problem of converting a piecewise polynomial from one representation to another. Using data reduction as illustrated in Example 6.5 in [Lyche & Marken '87b] and Example 4.1 in this paper, one can in principle approximate any piecewise polynomial by a B-spline curve or surface with relatively few knots. The problem is indeed a special case of Application 1.
T. Lyche, K. Morken / Knot removal
229
5. Modeling of curves. Data reduction can be used to obtain an accurate spline representation of a digitized curve. This has many applications. We mention the inking problem [Plass & Stone '83]. It could also be useful for designing fonts in automatic type setting systems. 6. Modeling of surfaces. Data reduction can also be used to fit a surface to points given on the surface in a regular fashion [Schmitt & Barsky '86]. We do believe that we can model most curves satisfactory, but in dealing with surfaces there are several problems which have to be dealt with. One problem is related to the 'squaring effect' in storage requirements which one encounters when going from curves to surfaces. Recall that in our multistage approximation scheme one has to sample (scan) the surface at a number of points so that the piecewise linear interpolant approximates the surface to well within the tolerance. For complex shapes it might be difficult to sample the surface at enough points. A possible way out of this problem is to use a more accurate local approximation than linear interpolation. Thus one could use a bicubic Hermite patch or a rectangular Coons patch [Faux & Pratt '79]. Alternatively some local spline approximation method [Lyche & Schumaker '75] could be used for the initial approximation. Another problem is the parametrization of data in space given along curves in space. Any parametrization must apply simultaneously to many curves. We do believe however, that the techniques described in this paper are likely to lead to schemes for rapid modeling of surfaces. 7. Postprocessing a spline approximation. In many applications an approximation or interpolation scheme results in a spline with too many knots. This is the case with for instance the smoothing spline, see, e.g., [de Boor '78, p. 235], which has a knot at every data point, and spline interpolation at a large n u m b e r of points. The knot removal process is well suited for computing an approximation with fewer knots within a given tolerance. One example of this can be found in [Lyche & Morken '87b, Example 6.1], where knots are removed from a smoothing spline. In practice we have often observed that the reduced data seem to be smoother (have fewer 'wiggles') than the original curve (surface). 8. Computing spline approximations to intersection curves. Many algorithms for computing intersection curves between surfaces produce a (large) number of points which somehow have to be represented as a curve. Approximation by data reduction can be applied to this problem. 9. Modeling of offset curves and surfaces. The modeling of offset curves and surfaces has been considered by many, see, e.g. [Tiller & Hanson '84] and [Arnold '86]. Our approach would be to first sample the offset curve or surface at a number of points, and then compute a spline approximation using the techniques described in this paper (cf. Example 3.2.). Using the triple knot technique, singularities parallel to the axes present no problem. 10. Data compression. If a given set of data points can be represented as a linear B-spline curve or bilinear tensor product surface, one can apply the algorithms in this paper with k = 2 to obtain a polygon with fewer points. Note that the points remaining after the reduction will be changed, but only within the given tolerance. In general, it would be recommended to use these changed data values instead of the original values in the reduced set.
Acknowledgement We would like to thank Tor D o k k e n and Richard F. Riesenfeld for useful discussions and comments.
230
T. Lyche, K. Morken / Knot removal
References Arnold, R. (1986), Quadratische und kubische Offset-B6zierkurven, Dissertation, Universit~it Dortmund. de Boor, C. (1976), Splines as linear combinations of B-splines. A survey in: Lorentz, G.G., Chui, C.K., and Schumaker, L.L., eds., Approximation Theory 11, Academic Press, New York, 1-47. de Boor, C. (1978), A Practical Guide to Splines, Springer, New York. Boehm, W. (1980), Inserting new knots into B-spline curves, Computer Aided Design 12, 199-201. Cohen, E., Lyche, T., and Riesenfeld, R. (1980), Discrete B-spline and subdivision techniques in computer-aided geometric design and computer graphics, Computer Graphics and Image Processing 14, 87-111. Cox, M.G. and Jones, H.M. (1987), Shape-preserving spline approximation in the 11 norm, in: Mason, J.C. and Cox, M.G., eds., Algorithms for Approximation, Clarendon Press, Oxford, 115-129. Farin, G., Rein, G., Sapidis, N., Worsey, A.J. (1987), Fairing cubic B-spline curves, Computer Aided Geometric Design 4, 91-103. Faux, I.D. and Pratt, M.J. (1979), Computational Geometry for Design and Manufacture, Ellis Horwood, Chichester, UK. Greville, T.N.E. (1961), Note on fitting of functions of several independent variables, J. Soc. Ind. Appl. Math. 9, 109-115. Kjellander, J.A.P. (1983), Smoothing of cubic parametric splines, Computer Aided Design 15, 175-179. Kjellander, J.A.P. (1983), Smoothing of bicubic parametric surfaces, Computer Aided Design 15, 288-293. Lyche, T. and Morken, K. (1986), Making the Oslo Algorithm more efficient, SIAM J. Numer. Anal. 23, 663-675. Lyche, T. and M~rken, K. (1987a), A discrete approach to knot removal and degree reduction algorithms for splines, in: Mason, J.C., and Cox, M.G., eds., Algorithms for Approximation, Clarendon Press, Oxford, 67-82. Lyche,T. and Morken, K. (1987b), A data reduction strategy for splines, Research Report no. 107 February 1987, Institute of Informatics, University of Oslo, to appear in IMA J. Numer. Anal. Lyche, T. and Schumaker, L.L. (1975), Local spline approximation methods, J. Approximation Theory 15, 294-325. Plass, M. and Stone, M. (1983), Curve-fitting with piecewise parametric cubics, Computer Graphics 17, 229-239. Schmitt, F.J. and Barsky, B.A. (1986), An adaptive subdivision method for surface fitting from sampled data, in: SIGRAPH "86 Conference Proceedings, Computer Graphics, 20, 179-188. Schumaker, L.L. (1981), Spline Functions: Basic Theory, Wiley, New York. Tiller, W. and Hanson, E.G. (1984), Offsets of two-dimensional profiles, IEEE Computer Graphics and Application 4, 36-46.