Computer Aided Geometric Design 4 (1987) 91-103 North-Holland
91
Fairing cubic B-spline curves G. FAR_IN 1, G. R E I N
2, N.
SAPIDIS, and A.J. WORSEY 3
Department of Mathematics, University of Utah, Salt Lake City, UT 84112, U.S.A. Presented at Wolfenbuettel 24 June 1986 by G. Farin Received 4 September 1986 Abstract. Algorithms are presented for locally fairing a B-spline curve in order to produce a pleasant curvature plot. The methods are based on inversions of the knot insertion algorithm for B-spline curves. Keywords. B-splines, curvature, smoothing, knot insertion, knot removal.
1. Introduction
1.1. Fair curves What does it take to make a curve look 'pleasant', 'sweet', or 'fair'? Probably the most competent people to provide an answer are those who deal with the shape of curves professionally: stylists and designers. Typically, they consider a curve to be 'fair' if it can be drawn using a small number of french curves. These are templates of curves with a particularly simple form, spirals, for example, since their curvature varies linearly with arc length. Mathematically, this notion of 'fairness' of a curve amounts to requiring that its curvature be almost piecewise linear, with only a small number of segments. Continuity of curvature is an obvious additional requirement since discontinuities give rise to aesthetic flaws in a curve that are easily detected by a trained designer. (A similar definition of fairness may be found in [Birkhoff '33, p. 74].) Therefore, it appears that curvature is an important tool in judging the aesthetic quality of a curve. We use it now to develop a concept of fairness which is integral to the curve fairing methods presented in this paper.
1.2. Relative fairness Recall that for a twice differentiable curve x(t), the curvature x is given by ~(t)=
IIx'(t) X x " ( t ) l l II x ' ( t ) II 3
If the curve is 'reasonably' parametrized, II x' II does not vary too much, and so • is mostly influenced by x " . For the sake of simplicity, we therefore restrict our attention to x " , rather 1 Present address: Dept. of Computer Science, Antona State University, Tempe, AZ 85 287, U.S.A. 2 Present address: Mathematisches Institut, Ludwig-Maximilians-Universitiit, Theresienstr. 39, D8000 Miinchen 2, West-Germany. 3 Present address: Dept. of Mathematical Sciences, University of North Carolina at Wilmington, Wilmington, NC 28403-3297, U.S.A. 0167-8396/87/$3.50 © 1987, Elsevier Science Publishers B.V. (North-Holland)
G. Farin et al. / Fairing cubic B-spline curves
92
than r, in order to assess the fairness of a curve, with the advantage that x " depends linearly on x, while K does not. If a fair curve should have regions of (almost) linear curvature, then its second derivative there should also be close to linear. For piecewise cubics, this implies that the third derivative in such a region should be close to constant and suggests the following definition of relative
fairness: Definition. Let x(t) and y(t) be two C 2 piecewise cubic curves defined over the interval [a, b]. We say that x is fairer than y at r ~ (a, b) if
[x'" ( r + ) - x ' " (r_)]2<~ [ y "" ( r + ) - y " " ( r _ ) ] 2, where r÷, r_ denote right and left hand limits, respectively. This definition is strictly local, since it may well be that x is fairer than y in some regions but not in others. It is fundamental to the results of this paper in that it is the basis of the algorithms for fairing C 2 cubic B-spline curves presented in Section 3. The methodology used in our approach derives from an existing fairing method due to [KjeUander '83a], together with the knot insertion algorithm for B-splines; see [Boehm '81] and [Boehm, Farin, K a h m a n n '84]. We review these schemes in the next section and also describe the underlying principle behind our fairing methods. Before doing so however, we conclude our introduction with a discussion which is pertinent to sections 4 and 5 of the paper.
1.3. Curoature plots During the design stage, a curve is usually displayed on a graphical display unit, which eliminates the possibility of judging curve quality directly: the curve is typically scaled down considerably, also screen resolution interferes with a realistic perception of the curve. On the other hand, a precision plot of the curve is expensive and time-consuming, and so a method is desirable that permits a designer to judge the quality of a curve directly on the screen. Curvature plots have turned out to be a reliable tool for this task; see [Irving '70] and [Kjellander '83a]. A curvature plot is the graph of curvature vs. arc length or vs. the parameter that defines the curve. For a three-dimensional curve, curvature is defined to be nonnegative, but for a two-dimensional curve, signed curvature is meaningful: it distinguishes between convex and concave regions and also provides an easier way to detect points of inflection than would unsigned curvature. In a typical design situation, one deals with the projections of a curve into the principal planes, and so can assign a sign to curvature. The curvature plot is a highly sensitive indicator of the shape of a curve: changes in a curve that are only noticeable to highly trained specialists are easily spotted on the curvature plot. If a curve has a 'pleasant' curvature plot, a designer can be sure that a precision plot of that curve will be flawless. Consequently, we shall use curvature plots to assess the quality of our fairing methods.
2. The fairing principle We shall consider curves that need to be faired because their curvature plots exhibit unwanted wiggles. Such curves are typically obtained by interpolation to points obtained from digitizing. This is the process of extracting point locations from either a physical object or from
G. Farin et aL / Fairing cubic B-spline curves
93
a drawing and can be done in many ways. For example, by using mechanical devices which determine the coordinates of points from tactile information or optical devices which employ laser techniques of stereo photography. All digitizing methods suffer from the inherent difficulty that the point locations they generate are only exact to within some device-dependent tolerance. Therefore any curve that is then passed through those points will exhibit an unwanted behavior due to the digitizing errors. Since the data points are known to be in error, it seems reasonable to change them within the known tolerance in order to obtain a curve that is fairer than the interpolant to the original data. This is the idea behind the following curve fairing scheme.
2.1. Kjellander's method We describe briefly an existing approach to fairing a C 2 piecewise cubic curve, due to [Kjellander '83a]. Referring to Fig. 1, let pj = x(tj) be a junction point on the curve x(t), i.e., two cubic segments meet at pj. Suppose the curve is to be faired at pj. The idea is to move pj to a 'better' location, determined by the following two steps: Step 1. Construct the cubic curve c that is uniquely determined by interpolation to the positions x(tj_l), x(tj+l) and to the tangents x ' ( t j _ l ) , x'(tj+l). This cubic c is obtained by cubic Hermite interpolation. Redefine the original curve x to agree with c over [tj_ 1, tj+l] and call the redefined curve x c. Step 2. Compute the C 2 cubic spline x r to the original data p~, but with pj being replaced by c(tj). Using this method, the modified curve x r is usually fairer than the original curve x at tj. Step 2 in the above algorithm is motivated by the following observation: the modified curve x c is C 3 at tj, thus it is fairer than x at tj. On the other hand, x c is only a C 1 curve due to the second derivative discontinuities introduced at tj_ 1 and tj+ 1- The goal, however, was to stay within the class of C 2 cubic splines. We annotate that this fairing method is global, due to the spline interpolation necessary to find the final curve x K. The fairing methods we propose, although based on the ideas behind Kjellander's method, have the desirable quality that they are local schemes. However, before moving on to discuss the details of our methods, it will be helpful to review the process of knot insertion for cubic B-spline curves, since it is the key to our approach to fairing. More details can be found in [Boehm '81] and [Boehm, Farin, K a h m a n n '84].
2.2. Knot insertion for cubic B-splines Let a B-spline curve x ( t ) be defined over the knot sequence { t i },.z'=0, where L is the number of curve segments. The associated B-spline, or de Boor, polygon for x ( t ) is determined by its coefficients, or de Boor points, {d~}~=+2. Throughout, we shall adopt the convention that the knots t o and t L have multiplicity three so that the B-spline polygon is tangent to the curve at the endpoints.
¢(tj)
S
xj+l~
Fig. 1. Kjellander's scheme: The point pj is replaced by e(tj). The solid line denotes the old curve, the thin line denotes the cubic Hermite interpolant to xj_ 1, p
xj-1, Xj+l,
t -1:3"+1.
94
G. Farin et aL / Fairing cubic B-spline curoes
To introduce the scheme of knot insertion, consider a B-spline curve defined over a knot sequence of strictly increasing points {to,..., tj_~, tj+l,..., t L ) and with a B-spline polygon { D i }. The same B-spline curve can also be defined over a refined knot sequence, that is, one with an additional knot tj; tj_l ~< tj < tj+l. This simply means that the cubic segment over [tj-1, tj+l] is split into two segments, defined over [tj_l, t/] and [tj, tj+l]. Although the curve does not change, its de Boor polygon does. The de Boor points d i of the new polygon are related to the original D; via the equations: di=D,; i = 0 , 1 .... , j - 3 ,
(--) ti+ 3
di=
tj
Di-1 +
li+3 - ti_ 1
d,+~=D,;
( tj-li-1)
Di;
i = j - 2, j - 1, j,
(2.1)
ti+3 - t i _ 1
i=j, j + I , . . . , L + I .
The application of KjeUander's method to smoothing a cubic B-spline curve can now be viewed from a different standpoint. The original curve x is defined over the knot sequence {t~ }~=0 where each knot has multiplicity one. The modified curve x c is defined over a somewhat different knot sequence. A particular knot, tj, is removed whilst the knots tj_~ and t/+l now have multiplicity two. Ultimately, the final curve xK is defined over the original knot sequence. The essential part in Kjellander's method is the knot removal, which can be viewed as the inverse process of knot insertion. When a knot tj is inserted into the knot sequence, the new de Boor points d;; i = 0 , . . . , L + 2, are determined from the given D~; i = 0 , . . . , L + 1, by the (L + 2) equations (2.1). On the other hand, if we want to remove tj, given the de Boor points d~, we must solve equations (2.1) for the unknowns D~ which amounts to solving the (L + 2 ) × (L + 1) overdetermined linear system ,:01
t~+3-__--t~I D,_ 1 + t2-3-_-~i_1 Di=di; D~=di+l;
i = j - 2, j - 1, j, i=j, j+I,...,L+I.
(2.2)
In general, this system - and hence the knot removal problem - has only approximate solutions. An exact solution exists if and only if the equations (2.2) are consistent, which is to say that the original curve was not simply C 2, but in fact C 3 at tj. Several methods could be used to obtain solutions for the unknowns D~ in (2.2). Typically, as with some of our methods, these might be approximation schemes based on some minimization criterion. However, regardless of the methodology used to solve (2.2), all the schemes lead to curve fairing algorithms that tackle the problem: Problem. Modify a given C 2 cubic B-spline curve locally so that the resulting curve is fairer than the original at a given knot tj, D3- 1
Di D j-3
uj-3
Fig. 2. Knot insertion:The old polygon is formed by the Di, the new polygon by the di.
G. Farinet aL / FairingcubicB-splinecurves
95
through the following two step procedure: Step I. Remove knot tj from the knot sequence and compute a new de Boor polygon approximating the original. Step 2. Using this polygon, reinsert tj so that the new curve is again defined over the original knot sequence. In the next section we shall present methods to perform Step 1, corresponding to a more geometric (Section 3.1) or to a more algebraic (Section 3.2) interpretation of the problem of knot removal.
3. Knot removal algorithms
3.1. Explicit algorithms Two obvious ways to solve equations (2.2) for the D~'s are to solve from the left with increasing i's, or from the right, with decreasing i's. These two approaches lead to the following two solutions.
D[=d~;
Dti=
i=0,..., j-
~-'t~-I
3,
di-
~--1
Di = di+l;
i=j+l,...,L+l,
Dr=di+l;
i=L+I
D[-1;
i=j-
2 .... , j ,
(3.1 l)
and .... , j ,
Dr 1= ( ti+3- ti-1 ) D:=d~;
tj-- ti-1
i = j . . . . . j - 2,
(3.1 r)
i=j-4,...,O,
see also Fig. 3. In general, since the original curve is n o t C 3 at tj, we have
O;_ 3 ~ dj_ 3 = DJ_ 3 and
DJ ~ dj+ I = D;.
We overcome this inconsistency by combining the intermediate results of (3.1 l) and (3.1 r) to obtain:
Di=d~;
i=0,..., j-
3,
D j _ 3 f DJ_3=dj_3, Dj_2 = aDJ_2 + (1 -
Ol)Dj-2, r
(3.1.2)
Dj_I = flDJ-1 + ( 1 - fl)Df-1, Dj=D/fdj+I, Di=di+l;
i=j+l
.... , L + I ,
where a and/~ are free parameters that must be suitably chosen. Based on these D~, we now reinsert the knot tj and denote the resulting (final) de Boor points by { d i )L=+02.Since of the original de Boor points only d j_ 2, dj_ 1, and d j, are changed by this process, the scheme is local. This method uses extrapolation to obtain the D,.'s from the given initial de Boor points d~.
96
G. Farin et aL / Fairing cubic B-spline curves
D~-2
J-V-
\ d3 ~
D~ D~_a
-
D3_3 o
Fig. 3. Knot removal: Starting with the polygon of the di, one can construct two polygons by starting either from d j_ 3 and proceeding to the right or by starting from d/+l and proceeding to the left. Consequently, it is in general a nontrivial p r o b l e m to choose the p a r a m e t e r s or, fl so that the shape of the original curve is retained as m u c h as possible. As Fig. 4 illustrates, an u n f o r t u n a t e choice for the p a r a m e t e r s m a y induce u n w a n t e d 'wiggles' in the curve. Based on various experimental results, we propose the following three choices for a a n d / 3 . 3.1(a)
a = 1, f l =
0.
As a consequence of this choice, dj_2=dj_2,
dj-~dj,
(3.1.3)
so that only dg_ a is changed. 3.10)
c~c~~+ c~c~ c~c~ + c~c~ + c~c~' C2C~
/~=
(3.1.4)
c?c~ + c?c~ + c7c~'
f
o
Fig. 4. Knot r~noval: unwanted wiggles appear in the new polygon if the polygons starting from the fight and starting from the left are blended together with constant weights.
G. Farinet al. / FairingcubicB-splinecurves
97
where
( tj+l--tJ )dj_3-dj_2+ ( tj-tj-3 )D r Cl = tj+a _ tj_3 tjS?'-..~ tj'~_3 j-2'
C3=
tj+ tj÷3 - tj_,
t j ÷ 3 - tj_
dj+,.
(Note that both a and 13 are automatically in [0,1].) The Cj are chosen so that they bound the difference vectors dj - dj. This choice for a, 13 minimizes
dist_lldj_2_dj_2llZ+l d j _ a _ d j _ , 2+ a~j_dj[12,
(3.1.6)
which follows from observing, after a lengthy calculation, that d i s t = (1 - a)2 C~ + ( a - 13)ZC~ + 132C32. Therefore, given the general approach of algorithm 3.1, with this scheme the new polygon is the least squares approximation to the original. It turns out that even without the restriction imposed on Dj - 2 and Dj_ 1 by (3.1.2) (they can only vary on straight lines), this method yields the least squares approximation to the original polygon. This will become obvious in Section 3.2. 3.1(0 Take a, 13~ [0, 1] so as to minimize: the distance of Dj_ 2 from the polygon leg d j_ 2 , d j_ 1 and also the distance of Dg_ 1 from the polygon leg d j_ 1, dj. In this case, it is more efficient to calculate Dj_2 and Dj_ 1 directly than to calculate the corresponding values of a and 13 first. Method (a) has the advantage that only the de Boor point d j_ 1 is changed and hence the scheme is as local and fast as possible. In many cases (see Example in Section 5) this is sufficient and the shape of the curve is not changed considerably. However, the method may run into the difficulties illustrated in Fig. 4. Method (b) is interesting in case one has a given tolerance for the change in the de Boor points, see section 4.2. But one should note that the new de Boor points may be very close to the original ones and still the new curve may exhibit unwanted 'wiggles'. This problem does not occur in method (c). It has been our experience that this method is the only one presented here which can handle examples like the one illustrated in Fig. 4 without inducing inflection points.
3.2. Matrix algorithms From an algebraic viewpoint, perhaps the most obvious approach to solving the overdetermined linear system (2.2), and hence to the problem of knot removal, is to use the method of least squares. The ( L + 2) × (L + 1) matrix, A say, arising from (2.2) for the unknowns Di, is of full rank. Therefore, there is a unique least squares solution to the equations. Moreover, because A is sparse and has a particularly simple structure, this solution can be easily
G. Farin et al. / Fairing cubic B-spline curves
98
computed. Specifically, we get
3.2(a) Di=d~;
i=O,...,j-4,
(3.2.1)
and Di = d~+l;
(3.2.2)
i=j+l,...,L+l,
whilst Dj_ 3, Dj_ 2, Dj_ 1, and Dj are computed as the least squares solution to the 5 × 4 linear system: 1 u 0 o 0
0 1-u v o 0
0 0 1-v w 0
0
o
{dj-3 idj_ 2
][ I oj_2
=
0 1 --w
1
DJ-1
//9:-
(3.2.3)
dj-a Idj dj+l
where, U=
( j__._+l t --tj ) tj+l __ tj-3 , tj+2-tj tj+ 2 __ tj_ 2
v=
w =
(3.2.4)
and
(3.2.5)
"
A severe drawback of this approach is that five of the initial de Boor points are altered, namely dj-3, ds-2, dj_,, dj, and dj+V Hence the method is not as local as the previous algorithm and in certain cases this may be undesirable. Moreover, one cannot smooth at the first or last interior knots in the knot sequence without modifying the first and last de Boor points. An obvious way to circumvent this possible difficulty, whilst retaining the nature of the least squares approximation, is to reduce the size of the linear system in (3.2.3). That is, we set
(3.2.7)
3.2~) Di=di; Di=di+l;
/lU
i = 0 .... , j - 3 , i=j,...,L
(3.2.8)
+ l,
/
and take Dj_ 2, /):._1 as the least squares solution to the 3 × 2 linear system:
v 0
1-v W
=
I jDj-1 ]
dj-1 Idj_ (l_w)dj+l
(3.2.9)
where u, v, and w are again given by (3.2.4)-(3.2.6). The solution may be obtained explicitly, and it turns out that it is the same as that generated by algorithm 3.1(b). Therefore, perhaps surprisingly, these two methods are exactly the same even though the motivations behind them are somewhat different. (Note that the quantity which is minimized by the least squares solution to (3.2.9) is actually equal to the quantity dist introduced in (3.1.6).) A natural alternative to the least squares approach for solving (3.2.9), is to seek solutions in
G. Farin et al. / Fairing cubic B-spline curves
99
the uniform norm and calculate minimax approximations. Indeed, given the least squares solution to (3.2.9), it is a straightforward matter to then calculate the minimax solution; see [Cheney '82]. As in (3.2.7) and (3.2.8), set 3.2(c)
Di=d~;
i=0,...,j-3,
(3.2.10)
and Di = d~+l;
i=j,...,L+l.
(3.2.11)
Then, given the least squares solution to (3.2.9), take Dj_ 2 and Dj_ 1 to be the minimax solution to the same overdetermined system. This involves solving a 2 × 2 linear system, see [Cheney '82]. Experimental evidence suggests that the question of whether or not this additional computational effort and expense is worthwhile is problem dependent. As with all the other schemes, once the new polygon has been calculated the last step in the algorithm is to insert the knot tj back into the knot sequence.
4. Practical considerations
We now comment on some practical issues that arise in the implementation of the above methods.
4.1. At which knot should we fair the curve? Several strategies are possible to determine at which knots the curve should be faired: (1) Find the largest jump discontinuity in the third derivative and fair there. (2) Find the largest slope discontinuity in the curvature and fair there. (3) Leave the decision as to where to fair to the designer. The first two options seem to be suited for use with batch mode computing, while interactive design really demands option (3). Option (2) should be used with care: a peak in curvature may well be a feature of the curve rather than an indication of misbehavior - in that case it should not be smoothed out.
4.2. Tolerances It frequently happens that the fairing process generates a curve that differs too much from its progenitor. This cannot happen if the difference vectors of the old and new de Boor points are small enough. Therefore, if one of them is larger than a prescribed tolerance, this difference vector should be used to move the corresponding de Boor point in that direction, but by no more than the given tolerance. Since the smoothing process may be repeated several times, one must be careful not to accumulate changes beyond the given tolerance.
4. 3. Limits of the fairing process The above methods were designed to remove 'noise' from data (or from the curvature plot). If this noise is not random, then either the digitizing process or the model itself were flawed, and one should not expect to obtain an acceptable curvature plot. In such cases, the data have to be modified by other methods after the fairing process has revealed a 'wrong trend' in the curvature plot.
lOO
G. Farin et al. / Fairing cubic B-spline curves
4. 4. Which algorithm should we use?
Algorithm 3.1(a) is the fastest of all and should probably be tried out first. If the result is not satisfactory (old and new polygon too far apart) one should try method 3.1(b) (or its matrix equivalent). If the new polygon is still not close enough, method 3.2(a) can be invoked to give the theoretically closest polygon to the original. This hierarchy of methods has decreasing locality, which should be taken into account if smoothing should happen near the curve endpoints. If unwanted wiggles appear, method 3.1(c) may be able to produce a more satisfactory result. These various choices can be implemented by means of a 'do again' option, which seems to be becoming popular in some CAD systems. However, if a systems programmer seeks a method that is easy to program and optimal in execution time, yet does not always yield optimal results, then method 3.1(a) is the one to pick. We have also experimented with algorithms that were more 'geometric' than the two last explicit algorithms. Such algorithms tend to be somewhat ad hoc, and although they perform reasonably well in most cases, we decided to restrict our presentation to methods with a firm theoretical foundation. We finally point out that cases may exist where only one data point is allowed to be changed. For such cases, Kjellander's method is the most appropriate one.
5. Examples We shall illustrate the above methods by considering two specific examples.
Example 1. We simulated digitized data from a curve that is a reasonable approximation to a feature line of a car trunk lid (Fig. 5(a)). The horizontal extension of the curve is 1200 ram, and the data points were in error by approximately 0.2 ram. The smoothing processes (with no tolerances prescribed)produced smoother curves that were within the required tolerance to the
original. Example 2. We simulated free-form design of a tennis racket. The number of curve segments used was relatively small (given the complicated shape of the racket), and deviations from the original could be handled generously. The vertical extension of the racket is 650 ram. The smoothing processes produced deviations of 36 mm. We should note here that those deviations are differences of old and new polygons, and that the deviations in the curves may be much smaller. This is illustrated in Fig. 6(a) and 6(d) by the fact that the de Boor points mainly move parallel to the curve.
(a)
Before smoothing Fig. 5(a). A scaled drawing of an example curve together with its B-spline polygon. The horizontal extension of the curve is approximately1200 ram1
G. Farin et al.
/ Fairing cubic B-spline curves
101
(b) O f
0"-"--'~
O
0
Fig. 5(b). Curvature plot of the curve from Fig. 5(a). The curvature plot reveals that ditizing errors are present in the curve.
(c) 0,
u
0
0
C
Fig. 5(c). Curvature plot after smoothing twice with method 3.1(b). The first time, all knots with an even index were smoothed, the second time all knots with an odd index. The maximum deviation of old and new polygon is approximately 0.2 ram.
(d) O
u
0
C
Fig. 5(d). The curve corresponding to Fig. 5(c) was smoothed only at the first interior knot. Note the 'nicer' curvature plot.
(e) u
&
8
0
F
Fig. 5(e). The original curve was smoothed with method 3.1(a) first at all even knots, then at all odd knots. The maximum deviation is now 0.3 ram.
~vocx
w ~
en I n ~ m l i e a
Am,~rc~cr,
102
G. Farin et al. / Fairing cubic B-spline curves
(a) E3. . . . . . . . . . . . . . . . . . . . . .
El
J
Before s m o o t h i n g
Before s m o o t h i n g
Fig. 6(a). A scaled drawing of a curve and its B-spline polygon. The vertical extension of the curve in the y-direction is approximately 650 ram.
Fig. 6(b). The corresponding curvature plot.
(c) [] . . . . . . . . . . . . . . . . . . . . .
E1
(d)
-%
J
M t e r s m o o t h i n g with m e t h o d 3.1.b) Fig. 6(c). The curve plus polygon corresponding to the curvature plot from Fig. 6(d). The maximal deviation of the old and new de Boor polygons is 36 ram.
Fig. 6(d). The curvature plot after smoothing at knots 2 and 14 with method 3.2(b).
G. Farin et al. / Fairing cubic B-spline curves
103
6. Conclusion and remarks (1) We have discussed methods to fair C 2 cubic B-spline curves. The key idea - remove a knot and then reinsert it - can also be applied to C 1 and even to C O cubic splines, as well as to higher degree B-spline curves. (2) Knot insertion and knot removal correspond respectively to increasing and decreasing the number of de Boor points of the B-spline polygon. Analogous procedures for B6zier curves would be degree elevation and degree reduction. This led us to define an obvious algorithm for fairing B6zier curves. Namely, reduce the polynomial degree by one and then elevate it. Although this is not, of course, a local procedure, results we obtained using this algorithm were similar to those obtained for B-spline curves. (3) KjeUander has generalized his curve smoothing method to tensor product surfaces [Kjellander '83b]. He did not make use, however, of the fact that his scheme is a linear one and thus has a natural tensor product generalization. This generalization is elementary and easy to implement, see [de Boor '78]. Similarly, method 3.1(a) can be generalized to a tensor product scheme. The other methods are not linear and do not permit a straightforward tensor product generalization.
Acknowledgements This research was supported in part by the Department of Energy with Contract no. DE-AC02-85ER12 046 and by NSF grant DCR-8502858 to the University of Utah. N. Sapidis was supported by a scholarship from the BODOSSAKI foundation in Athens, Greece.
References Boeb_m, W. (1981), Inserting new knots into B-spline curves, CAD 12, 199-201. Boehm, W., Farin, G. and Kahmann, J. (1984), A survey of curve and surface methods for CAGD, Computer Aided Geometric Design 1, 1-60. Birkhoff, G.D. (1933), Aesthetic Measure, Harvard University Press, Cambridge, MA. de Boor, C. (1978), A Practical G~uide to Splines, Springer, New York. Cheney, E.W. (1982), Introduction to Approximation Theory, Chelsea, New York. Irving, J.J. (1970), A system for designing and approximating aesthetically smooth curves with interactive graphics controls, Dissertation, Detroit. Kjellander, J.A. (1983a), Smoothing of cubic parametric splines, CAD 15, 175-179. KjeUander, J.A. (1983b), Smoothing of bicubic parametric surfaces, CAD 15, 288-294.