Fair interpolation and approximation of B-splines by energy ...
Computer-A/ded
0010~4485(95)ooo87-9
Design. Vol. 28. No. 9, pp. 7X3-760. 1996 Copyright 0 1994 Elsevier Science Ltd Printed in Great Britain. All rights reserved l3o10-44w96 $15.cQ+o.o0
ELSEVIER
Fair interpolation and approximation of B-splines by energy minimization and points insertion Tzvetomir
lvanov Vassilev
centripetal methods fail to give good results. Several other algorithms for ‘optimal’ parameterization have been described in the literatured-6 but once the optimum has been found, the curve is determined simply by using the interpolation conditions or by a least squares fit. No other shape constraints are applied. One approach for constructing a fair curve is to apply additional mathematical conditions to preserve its shape. For example Carl de Boor2 proposed the ‘taut spline’ that preserves convexity of the data, and Fritsch and Carlson’ constructed a C’ cubic interpolant that preserves the monotonicity. Another possibility is to use an interactive fairing scheme8,9, where the user decides which of the control points should be moved. Sapidis and Farin” suggested an automation of this process showing how to determine where a knot should be removed and a new one inserted using the curvature plots as a fairness criterion. Note that the resulting curve after the fairing is always only an approximation to the original one. An alternative approach exploits minimization of the energy of the curve (surface). This dates back to 1966 when Schweikert suggested the idea of a spline in tension”. Here the standard cubic spline is enriched by adding new exponential terms and degrees of freedom, called tension parameters. The works of Salkauskasl2 and Foley13 use a C’ interpolant that minimizes weighted L2-seminorms of the second, and the first and the second derivatives, respectively. Gelniker and Gossard14 constructed a shape that naturally resists stretching and bending. This leads to a functional expressed as an integral of a weighted sum of the fust and second parametric derivatives of the curve (or surface). In this method C ’ continuous cubic Hermite polynomials are used as a basis for curves and C’ triangular elements for surfaces. Nowacki and Lti” went further and adopted a fairness criterion that contains L2-seminorms of the second and the third order parametric derivatives. Their implementation uses C2 Hermite quintic polynomial curves with area constraints resulting in a non-linear system of equations that is solved numerically. An interpolation scheme, using the fairness norm14 together with Catmull-Clark surfaces, is described by
An efficient method for interpolation and approximation of both curve and surface points using B-splines is described. Automatic faking is presented based on minimizing an energy functional. Additional data points, used as degrees of freedom for the faking, are inserted only where the curve (the surface)
needs them. This reduces the number of the unknowns to a minimum which makes the algorithm very fast and efficient especially when a huge amount of data is concerned. Results of applying the algorithm for about 15,000 face data points, subject to measurement errors due to the digitization, are presented at the end of the paper. Copyright 0 1996 Elsevier Science Ltd Keywords: B-splines, interpolation, smoothing
INTRODUCTION To construct visually acceptable interpolations for sets of data points has been a long-standing research problem. The main difficulty arises from the fact that even when very efficient schemes based on splines are used the resulting curve (or surface) is often not ‘fair’ enough and it has extraneous ‘bumps’ and ‘wiggles’ (for example, see figures in this paper). The objective of the work, described in this paper, is to develop a method for fair spline interpolation and approximation for 3D curves and surfaces that is efficient enough to be applied on large amounts of data. Determining the parameterization of the curve is a fundamental problem to interpolation, Three basic methods have been widely used: uniform’, chord length2 and ‘centripetal’ parameterization3. Uniform knot sequence does not usually give satisfactory results because it does not take into account the geometry of the points. There are also cases where chord length and Department of Informatics, TU of Rousse, 8 Studentska Street, 7017 Rousse, Bulgaria Paper received: 10 March 1995. Revised: 28 November 1995
753
Interpolation
and approximation
of B-splines:
T I Vassilev
Halstead et ~1.‘~ for data sets with arbitrary topology. In order to add degrees of freedom for the fairing, the original control mesh is subdivided twice to decouple the interpolation constraints from each other. This, however, dramatically increases the number of control points. As a result the method is too slow and practically inapplicable for large amounts of data. This paper describes an algorithm for fair interpolation and approximation with C’ continuous cubic B-splines. The integral fairness criterion of Celniker and Gossard14 is applied for curves, but a new one is suggested for surfaces. Its main advantage is that it preserves the twodimensionality of the arrays which increases the efficiency. Additional data points (ADP) are used as degrees of freedom for minimizing the energy of the curve (surface), unlike other methods that use new control points. The ADP are not distributed evenly but inserted only where necessary, where a lack of data is observed. This results in a linear system of equations with a decreased number of unknowns. The parameterization is based on a new method for incorporating the data geometry and the method is presented for both curves and surfaces. The scheme was developed to be applicable for large amounts of data. The rest of the paper is organized in the following way. In the second section the construction of a fair spline curve is described, including evaluation of the energy functional, interpolation, approximation and choice of additional data points. In the third section the method is extended to tensor product surfaces. Finally, the results and conclusions are presented in the fourth and fifth sections, respectively.
FAIR B-SPLINE
CURVES
Evaluating the energy functional Let w(u) = [x(u), y(u), Z(U)] be a space curve parameterized by u and let V be a column vector of its B-spline control points. Celniker and Gossardr4 suggested the following energy functional
formula’
one can easily find that
.I ‘.
K,, =
”
7
-12
-1
7
34
-29
-12
-12
-29
34
7
1-1
-12
7
6 1
1
B,B, du = 120 .T
2
-3
0
1
-3
6
-3
0
o
-3
6
-3
1 ’ h,&!du J0
Kzs =
1
6
r
= A 6
1’
O
-3
1 (4)
21
Now. after we have expressed the local fairness measure for a single segment, the global fairness measure E for a B-spline curve can be found as a sum of the local fairness measures over each of the B-spline polynomial segments E = VTKV.
K = aK, + OK2
(5)
where K is a quadratic symmetric band matrix obtained from the local stiffness matrices. Minimizing Equation 5 for a particular choice of degrees of freedom will ensure the fairing of the curve.
Fair B-spline interpolation with ADP We will modify slightly the original problem of interpolation for the purposes of fairing. The new task is: given a sequence of data points pi, i = 1) . ~n - 1, some of them unknown and used as degrees of freedom, find the vector V of B-spline control points of a curve, which passes through the points pi and minimizes the fairness measure in Equation 5. This is demonstrated in Figuw I, where the original data set is {pI , , ps, pg. . p13} and the unknown data points are (p6, p7, ps}. Using the interpolation constraints with end conditions set as degrees of freedom, the control points can be found by solving AV = P
E =
Ub(aiv2(u) + @ti2(U))du s u‘l
(1)
where A is a coefficient
where W(u) and K(u) are the first and second derivatives in respect to the parameter u. The energy of the curve is represented as a weighted sum of its stretching and bending terms, where Q and p are freely selected coefficients. A single cubic B-spline segment can be written as w,(u) = VTB,, where B, is the B-spline basis functions vector. Evaluating the integral of Equation 1 for the function w, will result in the following quadratic form Es = V;K,V,
(6) matrix
and P is a column
vector
Y
(2)
where K, is a 4 x 4 symmetric matrix, whose entries are calculated by evaluating the integral of Equation 1. Es is called the ‘local fairness measure’ and K, is the ‘local stiffness matrix’ for the segment wS(u). K, can be represented as a weighted sum of its stretching and bending terms K, = oK~s + PK2s For
754
a uniform
X
(3) knot
distribution,
using
de Boor’s
Figure I
Data set for fair interpolation
with additional
data points
interpolation
and approximation
of B-splines:
T I Vassilev
cannot be performed for physical reasons. As shown by the various experiments the author conducted, the curve tends to have extraneous inflection points and undulations in and around these intervals. Therefore the additional data points should be inserted precisely where a lack of data is observed. Figure 2 demonstrates the insertion algorithm for the data set in Figure 1. It is as follows. (1) Calculate the parameter values u,, i = 1, . II, using the chord length or centripetal formula (Figure 2~). (2) Compute the average length between two parameter values Au = u,/(n - 1). (3) For each interval [u,, ul+,] find the number n, = round(Au,/Au) of the new parameter values (data points) and insert them (Figure 2h). (4) Recalculate the sequence u, this time using the uniform formula (Figure 2~). We could use the parameterization from Step 3 to perform the actual interpolation (or approximation). Ln this case, however, we will have to use a numerical method to evaluate the integral in Equation 1, e.g. Gaussian quadrature18. That is why the author suggests that Step 4 should be carried out. Then the stiffness matrix K can be computed in advance and directly incorporated as a constant array into the program. This increases the efficiency especially when a large data set is concerned. The final parameterization (Step 4) is uniform in respect to the ADP-enriched data set but it is non-uniform in respect to the original one (see the sequence {u,. us, uy.. , u13} in Figure 2c and the data set {p, ~. ,ps,p9.. . pli} in Figuw I). In fact, the knot vector (Figure 2c) is only an approximation of the chord length parameterization in Figure 2h. but it is good enough, because the smoothness of the curve is ensured by the minimizing of the fairness measure in Equation S.
FAIR B-SPLINE
-104
4 k
Figure 4 intcrpolant line)
Data set I. Cubic B-spline interpolant (dashed line); fair with ACP (dotted line); fair interpolant with ADP (solid
V be a tzl, x n,. matrix of its control points. The following energy functional is considered in this paper
SURFACES
Fairness criterion Let w(u, U) be a surface parameterized
by u and IS.and let
where the suffixes mean partial derivatives in respect to the parameters u and ?I. A tensor product B-spline surface can be expressed as wju, 11)= BIfVB,. where B,, and functions.
(17)
B,. are
vectors
of the
B-spline
basis
Interpolation and approximation with ADP The interpolation (approximation) the following system:
constraints
result in
A,,VA,, = D or V = A, ‘DA,. ’
(18)
Suppose we have inserted additional data points only in the direction of the parameter u (see Figure 3), then
Table 1
Figure 3 parameter
Surface data set with new data points inscrted 111direction 11
01
.\ )
Data set
0.0
2.0
0.0
0.5
I, Figure4 3.0 6.5
4.0 7.6
5.0 6.5
6.0 1.2
7.0 0.62
8.0 0.6
9.0 0.6
Interpolation
and approximation
of B-splines:
T I Vassilev
Y i
x
7-e54-32I--
X
a
b
2
Figure 6 Face data set: (a) least-squares squares approximation with ADP
IO-
approximation;
(b) fair least-
16.
freedom need to be slightly modified for tensor product surfaces. New parameter values (data points) for the sequence ui will appear in each one of the n, isoparametric curves (see Figure 3), which may be computed as follows:
lo-6-
X 0
Figure 5 Data set 2. Cubic B-spline interpolant interpolant (I = 1, /3 = 0.2 (dotted line); ADP p = 1 (solid line)
(dashed line); ADP interpolant a = 0.2,
Equation 18 becomes: V = A;‘(D,X,
+ Ds)A,*
(19)
where X, is a r, x n, matrix of the unknowns and D, is a II, x n, matrix of the known data points. Substituting Equations 17 and 19 in Equation 16 and finding its minimum leads to a system D;rKAUD,X, + D;rKAUDO= 0
(20)
where KAu and D, are the same as in Equation 11 but computed for the knot sequence Ui. Finding the solution of Equation 20 requires solving rz, linear systems with r, unknowns. Note that the system matrix DzKAUD, is common and it can be computed and decomposed in advance which speeds up the method. A detailed proof of Equation 20 is given in the Appendix. If we insert additional data points in the direction of w we will obtain nearly the same system as Equation 20 but with matrices computed for the sequence vi. Once all the data points have been found the control points are calculated from Equation 18.
Parameterization and choice of degrees of freedom The parameterization
and the choice of degrees of
Table 2
I, 2 and 5
x Y
0.0 0.0
Data set 2, Figures 1.0 0.1
2.0 0.2
3.0 0.4
4.0 1.0
5.0 6.6
6.0 7.2
7.0 7.4
8.0 7.5
9.0 7.6
(1) For each interval [ui, Ui+,], using Steps l-3 from the previous section for each curve, find the number = lElX(nij,~ = l,...,n,). %max (2) For each interval [Ui,Ui+1]insert nimaxnew parameter values. (3) Recalculate the sequence Ui, using the uniform formula Ui = i, thus making the parameterization common for all the isoparametric curves, which is necessary for tensor product surfaces. Controlling the error of approximation can be done in the same way as for curves but this time the values L, and L, for both parameter directions should be varied.
RESULTS The algorithms were implemented on a Silicon Graphics Indy workstation using the open GL library to render the surfaces. The results are shown in Figures 4-8. The original data points are marked with crosses and the new ones (inserted by the algorithm) with squares. The curvature and slope distributions are given below every curve as indicators for the fairness. All the experiments (except Figure 5) were conducted with coefficient values inEquation 1 (~=/3= 1. An alternative method for fair interpolation has been implemented that follows from the idea for inserting Table 3 curves
Performance
of the different
Method
Fair interpolation Fair interpolation Approximation Fair approximation
methods
Figure
with ACP with ADP with ADP
Figure 6a Figure 6b
for approximating
Number of polynomial segments L
Time (s)
354 125 70 70
71.3 4.9 0.05 1.2
757
Interpolation
and approximation
of B-splines:
T I Vassilev
(4
(b) Figure 7 (a) B-spline interpolation
surface
interpolation:
(b) fair B-spline
surface
enough additional control vertices so that the interpolation constraints are decoupled from each other16. This method uses the same energy functional but additional control points (ACP) as degrees of freedom. Figure 4 compares the quality of the schemes with the data set given in Table 1. Type 1 end conditions (given tangent vectors) have been imposed for the cubic B-spline interpolation. The method with ACP obviously improves the smoothness of the curve but it neither cures the violated monotonicity in interval one nor removes the extraneous inflection points in intervals 6, 7. Figure 5 demonstrates the results with another data set (Table 2). It also shows the influence of the parameters N and p, which can be used as shape controls. If N > fl the seminorm of the first derivative takes a bigger part in the
758
Figure 8 (a) Least-squares surface approximation for face data set; (b) fair least-squares surface approximation for face data set
fairness measure and the interpolant (dotted line) is more stretched (as if somebody pulls it at the ends). This effect can be best observed for the limit case (/I = 0) but then the interpolant might not be smooth enough. When (I < ,!3the second derivative seminorm takes the bigger part and the resulting curve (solid line) is more bent. Fair B-spline least squares approximation (Figure 6) was applied on 119 points of a face profile obtained through digitization and subject to substantial noise. The effect of fairing is obvious in the areas of the forehead and the mouth. The performance of the different methods for this data set is given in Table 3. The efficiency of the scheme with ADP improves by a factor of 14 compared with that of ACP. Figure 7 shows B-spline surface interpolation of the
Interpolation Table 4
Performance
and approximation
of B-splines:
T I Vassilev
of the different methods for approximating
surfaces Method
Fair interpolation Approximation Fair approximation
Figure
Figure 8a Figure 8b
Number of segments
Farin, G Curves and Surfaces for Computer Aided Geometric Design Academic Press (1993) de Boor, C A Practical Guide to Splines Springer-Verlag, New
Time (s)
L”
L”
130
130
21.1
56 56
70 70
3.2 6.8 6
data points in Figure 3. The unwanted bumps and wiggles (Figure 7~2)were removed by the fairing (Figure 7b). Least squares B-spline fair approximation has been
applied on a ‘noisy’ data set, of about 15,000 points. The extraneous undulations (Figure 8~) were smoothed by the algorithm (Figure 8b), automatically inserting new data points where necessary. The performance of the different methods for this data set is given in Table 4. It is obvious that the fairing is not too expensive compared to the pure approximation, the quality, however, is much better.
10 11 12
CONCLUSIONS
13 14
A method for fair interpolation and approximation of both curve and surface points using B-splines has been described. The suggested energy functional for surfaces (Equation 16) preserves the two-dimensionality of the arrays of control and data points. Solving this kind of system, as Farin’ shows, leads to 0(n4) computations, which are normally 0(n6) if the points are arranged in a vector. Although the algorithms were implemented for the cubic case, they could be easily extended for splines of any order k. Additional data points (ADP) are inserted as degrees of freedom for the fairing only where necessary. This approach proved to be very suitable for the case of measured (or digitized) data when the points have been originally planned to be evenly distributed, but for some practical reason measurements cannot be performed in all of the intervals. An important feature is that the algorithms for ADP insertion and for the actual fairing are independent. So, if for some practical applications another insertion scheme proves to be more appropriate, it can be used without any limitations with the described fairing procedure. The idea of ADP can be extended to other methods for global fair interpolation (approximation). It will improve the performance, especially when a larger amount of data is concerned.
ACKNOWLEDGEMENTS The author would like to thank the British Council and Microsoft who sponsored his work. Figures 6 and 8 are produced from data obtained from a digitization of a Spitting Image puppet of Margaret Thatcher. Thanks to Spitting Image production for providing the puppet and to the UK Advisory Committee on Computer Graphics for funding the digitization. The work described in this paper was done at Queen Mary and Westfield College, University of London. The author is very grateful to his supervisor Mel Slater and to the referees for their help and useful comments.
York (1978) Lee, E ‘Choosing nodes in parametric curve interpolation’ Comput.-Aided Des. Vol21 No 6 (1989) pp 363-370 Hartley, P and Judd, C ‘Parameterization and shape of B-spline curves’ Comput.-Aided Des. Vol 12 No 5 (1980) pp 235-238 Nielson, G and Foley, T ‘A survey of applications of an affine invariant norm’ in Lyche, T. and Shumaker, L. (Eds.) Mathematical Methods in CAGD Academic Press (1989) pp 445-468 Marin, S P ‘An approach to data parameterization in parametric cubic spline interpolation problems’ Approx. Theory Vol41 (1984) pp 64-86 Fritsch, F N and Carlson, R E ‘Monotone piecewise cubic internolation’ SIAM J. Numer. Anal. Vol 17 No 2 (1980) pp 238-246 Kjellander, J A P ‘Smoothing of cubic parametric splines’ Comput.-Aided Des. Vol 15 No 5 (1983) pp 175-179 Rogers. D F. Satterfield. S G and Rodrieuez, F A ‘Shiphull, B-&ink surfaces and CAD/CAM' in Computer Graphics-Theory and Applications Springer-Verlag (1983) Sapidis, N and Farin, G ‘Automatic fairing algorithm for B-spline curves’ Comput.-Aided Des. Vo122 No 2 (1990) pp 121-129 Schweikert, D G ‘An interpolation curve using a spline in tension’ J. Math. Phys. Vol45 (1966) pp 312-317 Salkauskas, K ‘C’ splines for interpolation of rapidly varying data’ Rocky Mountain J. Math. Vol 14 (1984) pp 239-250 Foley, T A ‘A shape preserving interpolant with tension controls’ Comma. Aided Geom. Des. Vol 5 (1988) DD 105-l 18 Celniker, G and Gossard, D ‘Deformable curve and surface finite elements for free-form shape design’ Proc. SIGGRAPHVol25 No 4 (1991) pp 257-266
15 16 17 18
Nowacki, H and Lii, X ‘Fairing composite polynomial curves with constraints’ Comput. Aided Geom. Des. Vol 11 (1994) pp 1-15 Halstead, M, Kass, M and DeRose, T ‘Efficient, fair interpolation using Catmull-Clark surfaces’ Proc. SIGGRAPH, Annual series (1993) pp 35-44 Forsythe, G and Moler, C Computer Solution of Linear Algebraic Systems Prentice-Hall (1967) pp 114-l 19 Press, W H, Flannery, B P, Teukolsky, S A and Vetterling, W T Numerical Recipes, The Art of Scientt$c Computations Cambridge University Press, New York (1989) pp 121-126
APPENDIX Equation 16 is a sum of four similar terms and here the solution to the first one is given. From Equation 17 w:v = BTVTB 2) uBTVB 11 ‘v
(Al)
where the dots mean first derivative. We have to compute d cl w;,dudv = cl 8% d = Clax,,
f f
B;A;~(D,x,
B=V=i B=VB dudv Cl u I4 21
+ Do)=A;=IiUti;fA,’
(D,X, + Ds)A,%,dudv = 2c,
D~A,=i3,li~A,‘(D,X,
+ D,,)
A,‘&,iS;A,=dudv = 2c,D;A;=
(ji.B:du)A;l(D,,Xu
= 2c,D;fA,TK1uA;1(D,X,
+ D,,)A,’
+ Do)A,‘KI,A,=
(-42)
759
Interpolation
and approximation
of B-splines:
T I Vassilev
Results similar to Equation A2 will be obtained evaluating the remaining three terms of Equation 16. Then the system of equations reads (++ means if and only 4
D:A&,KI,
+
LZKd A,’ (DuX,, + Do)A, ’
((Y,K,~ + ,/&K2,)A,T = 0 # D;A,TK,A,‘(D,X, D:K,,(D,X,
+ Do)A,‘K,,A,.T
+ Do)KA,. = 0
= 0 ++ (A3)
Here cl = (Y,,cY,.,cz = ,&,a,,, c3 = c@,,, c4 = /L&. The matrix KAt, is symmetric, band and positive definite and its inverse KA:, exists. Multiplying both sides of
760
Equation
A3 with K,&! results in
D,TKA~D,X, + D,TKAuDo= 0
(Ad)
Tzvetomir Ivanov Vussilev gruduuted jiom the Technicrrl b’niversi!,~ qf Roussr, Bulgaria in 1990. He i.5 currently a lecturer at the Department of Informatics at the same university. His research interests include computer-aided geometric d~@yn, data processing, and computer graphics
0010~4485(95)ooo87-9
Design. Vol. 28. No. 9, pp. 7X3-760. 1996 Copyright 0 1994 Elsevier Science Ltd Printed in Great Britain. All rights reserved l3o10-44w96 $15.cQ+o.o0
ELSEVIER
Fair interpolation and approximation of B-splines by energy minimization and points insertion Tzvetomir
lvanov Vassilev
centripetal methods fail to give good results. Several other algorithms for ‘optimal’ parameterization have been described in the literatured-6 but once the optimum has been found, the curve is determined simply by using the interpolation conditions or by a least squares fit. No other shape constraints are applied. One approach for constructing a fair curve is to apply additional mathematical conditions to preserve its shape. For example Carl de Boor2 proposed the ‘taut spline’ that preserves convexity of the data, and Fritsch and Carlson’ constructed a C’ cubic interpolant that preserves the monotonicity. Another possibility is to use an interactive fairing scheme8,9, where the user decides which of the control points should be moved. Sapidis and Farin” suggested an automation of this process showing how to determine where a knot should be removed and a new one inserted using the curvature plots as a fairness criterion. Note that the resulting curve after the fairing is always only an approximation to the original one. An alternative approach exploits minimization of the energy of the curve (surface). This dates back to 1966 when Schweikert suggested the idea of a spline in tension”. Here the standard cubic spline is enriched by adding new exponential terms and degrees of freedom, called tension parameters. The works of Salkauskasl2 and Foley13 use a C’ interpolant that minimizes weighted L2-seminorms of the second, and the first and the second derivatives, respectively. Gelniker and Gossard14 constructed a shape that naturally resists stretching and bending. This leads to a functional expressed as an integral of a weighted sum of the fust and second parametric derivatives of the curve (or surface). In this method C ’ continuous cubic Hermite polynomials are used as a basis for curves and C’ triangular elements for surfaces. Nowacki and Lti” went further and adopted a fairness criterion that contains L2-seminorms of the second and the third order parametric derivatives. Their implementation uses C2 Hermite quintic polynomial curves with area constraints resulting in a non-linear system of equations that is solved numerically. An interpolation scheme, using the fairness norm14 together with Catmull-Clark surfaces, is described by
An efficient method for interpolation and approximation of both curve and surface points using B-splines is described. Automatic faking is presented based on minimizing an energy functional. Additional data points, used as degrees of freedom for the faking, are inserted only where the curve (the surface)
needs them. This reduces the number of the unknowns to a minimum which makes the algorithm very fast and efficient especially when a huge amount of data is concerned. Results of applying the algorithm for about 15,000 face data points, subject to measurement errors due to the digitization, are presented at the end of the paper. Copyright 0 1996 Elsevier Science Ltd Keywords: B-splines, interpolation, smoothing
INTRODUCTION To construct visually acceptable interpolations for sets of data points has been a long-standing research problem. The main difficulty arises from the fact that even when very efficient schemes based on splines are used the resulting curve (or surface) is often not ‘fair’ enough and it has extraneous ‘bumps’ and ‘wiggles’ (for example, see figures in this paper). The objective of the work, described in this paper, is to develop a method for fair spline interpolation and approximation for 3D curves and surfaces that is efficient enough to be applied on large amounts of data. Determining the parameterization of the curve is a fundamental problem to interpolation, Three basic methods have been widely used: uniform’, chord length2 and ‘centripetal’ parameterization3. Uniform knot sequence does not usually give satisfactory results because it does not take into account the geometry of the points. There are also cases where chord length and Department of Informatics, TU of Rousse, 8 Studentska Street, 7017 Rousse, Bulgaria Paper received: 10 March 1995. Revised: 28 November 1995
753
Interpolation
and approximation
of B-splines:
T I Vassilev
Halstead et ~1.‘~ for data sets with arbitrary topology. In order to add degrees of freedom for the fairing, the original control mesh is subdivided twice to decouple the interpolation constraints from each other. This, however, dramatically increases the number of control points. As a result the method is too slow and practically inapplicable for large amounts of data. This paper describes an algorithm for fair interpolation and approximation with C’ continuous cubic B-splines. The integral fairness criterion of Celniker and Gossard14 is applied for curves, but a new one is suggested for surfaces. Its main advantage is that it preserves the twodimensionality of the arrays which increases the efficiency. Additional data points (ADP) are used as degrees of freedom for minimizing the energy of the curve (surface), unlike other methods that use new control points. The ADP are not distributed evenly but inserted only where necessary, where a lack of data is observed. This results in a linear system of equations with a decreased number of unknowns. The parameterization is based on a new method for incorporating the data geometry and the method is presented for both curves and surfaces. The scheme was developed to be applicable for large amounts of data. The rest of the paper is organized in the following way. In the second section the construction of a fair spline curve is described, including evaluation of the energy functional, interpolation, approximation and choice of additional data points. In the third section the method is extended to tensor product surfaces. Finally, the results and conclusions are presented in the fourth and fifth sections, respectively.
FAIR B-SPLINE
CURVES
Evaluating the energy functional Let w(u) = [x(u), y(u), Z(U)] be a space curve parameterized by u and let V be a column vector of its B-spline control points. Celniker and Gossardr4 suggested the following energy functional
formula’
one can easily find that
.I ‘.
K,, =
”
7
-12
-1
7
34
-29
-12
-12
-29
34
7
1-1
-12
7
6 1
1
B,B, du = 120 .T
2
-3
0
1
-3
6
-3
0
o
-3
6
-3
1 ’ h,&!du J0
Kzs =
1
6
r
= A 6
1’
O
-3
1 (4)
21
Now. after we have expressed the local fairness measure for a single segment, the global fairness measure E for a B-spline curve can be found as a sum of the local fairness measures over each of the B-spline polynomial segments E = VTKV.
K = aK, + OK2
(5)
where K is a quadratic symmetric band matrix obtained from the local stiffness matrices. Minimizing Equation 5 for a particular choice of degrees of freedom will ensure the fairing of the curve.
Fair B-spline interpolation with ADP We will modify slightly the original problem of interpolation for the purposes of fairing. The new task is: given a sequence of data points pi, i = 1) . ~n - 1, some of them unknown and used as degrees of freedom, find the vector V of B-spline control points of a curve, which passes through the points pi and minimizes the fairness measure in Equation 5. This is demonstrated in Figuw I, where the original data set is {pI , , ps, pg. . p13} and the unknown data points are (p6, p7, ps}. Using the interpolation constraints with end conditions set as degrees of freedom, the control points can be found by solving AV = P
E =
Ub(aiv2(u) + @ti2(U))du s u‘l
(1)
where A is a coefficient
where W(u) and K(u) are the first and second derivatives in respect to the parameter u. The energy of the curve is represented as a weighted sum of its stretching and bending terms, where Q and p are freely selected coefficients. A single cubic B-spline segment can be written as w,(u) = VTB,, where B, is the B-spline basis functions vector. Evaluating the integral of Equation 1 for the function w, will result in the following quadratic form Es = V;K,V,
(6) matrix
and P is a column
vector
Y
(2)
where K, is a 4 x 4 symmetric matrix, whose entries are calculated by evaluating the integral of Equation 1. Es is called the ‘local fairness measure’ and K, is the ‘local stiffness matrix’ for the segment wS(u). K, can be represented as a weighted sum of its stretching and bending terms K, = oK~s + PK2s For
754
a uniform
X
(3) knot
distribution,
using
de Boor’s
Figure I
Data set for fair interpolation
with additional
data points
interpolation
and approximation
of B-splines:
T I Vassilev
cannot be performed for physical reasons. As shown by the various experiments the author conducted, the curve tends to have extraneous inflection points and undulations in and around these intervals. Therefore the additional data points should be inserted precisely where a lack of data is observed. Figure 2 demonstrates the insertion algorithm for the data set in Figure 1. It is as follows. (1) Calculate the parameter values u,, i = 1, . II, using the chord length or centripetal formula (Figure 2~). (2) Compute the average length between two parameter values Au = u,/(n - 1). (3) For each interval [u,, ul+,] find the number n, = round(Au,/Au) of the new parameter values (data points) and insert them (Figure 2h). (4) Recalculate the sequence u, this time using the uniform formula (Figure 2~). We could use the parameterization from Step 3 to perform the actual interpolation (or approximation). Ln this case, however, we will have to use a numerical method to evaluate the integral in Equation 1, e.g. Gaussian quadrature18. That is why the author suggests that Step 4 should be carried out. Then the stiffness matrix K can be computed in advance and directly incorporated as a constant array into the program. This increases the efficiency especially when a large data set is concerned. The final parameterization (Step 4) is uniform in respect to the ADP-enriched data set but it is non-uniform in respect to the original one (see the sequence {u,. us, uy.. , u13} in Figure 2c and the data set {p, ~. ,ps,p9.. . pli} in Figuw I). In fact, the knot vector (Figure 2c) is only an approximation of the chord length parameterization in Figure 2h. but it is good enough, because the smoothness of the curve is ensured by the minimizing of the fairness measure in Equation S.
FAIR B-SPLINE
-104
4 k
Figure 4 intcrpolant line)
Data set I. Cubic B-spline interpolant (dashed line); fair with ACP (dotted line); fair interpolant with ADP (solid
V be a tzl, x n,. matrix of its control points. The following energy functional is considered in this paper
SURFACES
Fairness criterion Let w(u, U) be a surface parameterized
by u and IS.and let
where the suffixes mean partial derivatives in respect to the parameters u and ?I. A tensor product B-spline surface can be expressed as wju, 11)= BIfVB,. where B,, and functions.
(17)
B,. are
vectors
of the
B-spline
basis
Interpolation and approximation with ADP The interpolation (approximation) the following system:
constraints
result in
A,,VA,, = D or V = A, ‘DA,. ’
(18)
Suppose we have inserted additional data points only in the direction of the parameter u (see Figure 3), then
Table 1
Figure 3 parameter
Surface data set with new data points inscrted 111direction 11
01
.\ )
Data set
0.0
2.0
0.0
0.5
I, Figure4 3.0 6.5
4.0 7.6
5.0 6.5
6.0 1.2
7.0 0.62
8.0 0.6
9.0 0.6
Interpolation
and approximation
of B-splines:
T I Vassilev
Y i
x
7-e54-32I--
X
a
b
2
Figure 6 Face data set: (a) least-squares squares approximation with ADP
IO-
approximation;
(b) fair least-
16.
freedom need to be slightly modified for tensor product surfaces. New parameter values (data points) for the sequence ui will appear in each one of the n, isoparametric curves (see Figure 3), which may be computed as follows:
lo-6-
X 0
Figure 5 Data set 2. Cubic B-spline interpolant interpolant (I = 1, /3 = 0.2 (dotted line); ADP p = 1 (solid line)
(dashed line); ADP interpolant a = 0.2,
Equation 18 becomes: V = A;‘(D,X,
+ Ds)A,*
(19)
where X, is a r, x n, matrix of the unknowns and D, is a II, x n, matrix of the known data points. Substituting Equations 17 and 19 in Equation 16 and finding its minimum leads to a system D;rKAUD,X, + D;rKAUDO= 0
(20)
where KAu and D, are the same as in Equation 11 but computed for the knot sequence Ui. Finding the solution of Equation 20 requires solving rz, linear systems with r, unknowns. Note that the system matrix DzKAUD, is common and it can be computed and decomposed in advance which speeds up the method. A detailed proof of Equation 20 is given in the Appendix. If we insert additional data points in the direction of w we will obtain nearly the same system as Equation 20 but with matrices computed for the sequence vi. Once all the data points have been found the control points are calculated from Equation 18.
Parameterization and choice of degrees of freedom The parameterization
and the choice of degrees of
Table 2
I, 2 and 5
x Y
0.0 0.0
Data set 2, Figures 1.0 0.1
2.0 0.2
3.0 0.4
4.0 1.0
5.0 6.6
6.0 7.2
7.0 7.4
8.0 7.5
9.0 7.6
(1) For each interval [ui, Ui+,], using Steps l-3 from the previous section for each curve, find the number = lElX(nij,~ = l,...,n,). %max (2) For each interval [Ui,Ui+1]insert nimaxnew parameter values. (3) Recalculate the sequence Ui, using the uniform formula Ui = i, thus making the parameterization common for all the isoparametric curves, which is necessary for tensor product surfaces. Controlling the error of approximation can be done in the same way as for curves but this time the values L, and L, for both parameter directions should be varied.
RESULTS The algorithms were implemented on a Silicon Graphics Indy workstation using the open GL library to render the surfaces. The results are shown in Figures 4-8. The original data points are marked with crosses and the new ones (inserted by the algorithm) with squares. The curvature and slope distributions are given below every curve as indicators for the fairness. All the experiments (except Figure 5) were conducted with coefficient values inEquation 1 (~=/3= 1. An alternative method for fair interpolation has been implemented that follows from the idea for inserting Table 3 curves
Performance
of the different
Method
Fair interpolation Fair interpolation Approximation Fair approximation
methods
Figure
with ACP with ADP with ADP
Figure 6a Figure 6b
for approximating
Number of polynomial segments L
Time (s)
354 125 70 70
71.3 4.9 0.05 1.2
757
Interpolation
and approximation
of B-splines:
T I Vassilev
(4
(b) Figure 7 (a) B-spline interpolation
surface
interpolation:
(b) fair B-spline
surface
enough additional control vertices so that the interpolation constraints are decoupled from each other16. This method uses the same energy functional but additional control points (ACP) as degrees of freedom. Figure 4 compares the quality of the schemes with the data set given in Table 1. Type 1 end conditions (given tangent vectors) have been imposed for the cubic B-spline interpolation. The method with ACP obviously improves the smoothness of the curve but it neither cures the violated monotonicity in interval one nor removes the extraneous inflection points in intervals 6, 7. Figure 5 demonstrates the results with another data set (Table 2). It also shows the influence of the parameters N and p, which can be used as shape controls. If N > fl the seminorm of the first derivative takes a bigger part in the
758
Figure 8 (a) Least-squares surface approximation for face data set; (b) fair least-squares surface approximation for face data set
fairness measure and the interpolant (dotted line) is more stretched (as if somebody pulls it at the ends). This effect can be best observed for the limit case (/I = 0) but then the interpolant might not be smooth enough. When (I < ,!3the second derivative seminorm takes the bigger part and the resulting curve (solid line) is more bent. Fair B-spline least squares approximation (Figure 6) was applied on 119 points of a face profile obtained through digitization and subject to substantial noise. The effect of fairing is obvious in the areas of the forehead and the mouth. The performance of the different methods for this data set is given in Table 3. The efficiency of the scheme with ADP improves by a factor of 14 compared with that of ACP. Figure 7 shows B-spline surface interpolation of the
Interpolation Table 4
Performance
and approximation
of B-splines:
T I Vassilev
of the different methods for approximating
surfaces Method
Fair interpolation Approximation Fair approximation
Figure
Figure 8a Figure 8b
Number of segments
Farin, G Curves and Surfaces for Computer Aided Geometric Design Academic Press (1993) de Boor, C A Practical Guide to Splines Springer-Verlag, New
Time (s)
L”
L”
130
130
21.1
56 56
70 70
3.2 6.8 6
data points in Figure 3. The unwanted bumps and wiggles (Figure 7~2)were removed by the fairing (Figure 7b). Least squares B-spline fair approximation has been
applied on a ‘noisy’ data set, of about 15,000 points. The extraneous undulations (Figure 8~) were smoothed by the algorithm (Figure 8b), automatically inserting new data points where necessary. The performance of the different methods for this data set is given in Table 4. It is obvious that the fairing is not too expensive compared to the pure approximation, the quality, however, is much better.
10 11 12
CONCLUSIONS
13 14
A method for fair interpolation and approximation of both curve and surface points using B-splines has been described. The suggested energy functional for surfaces (Equation 16) preserves the two-dimensionality of the arrays of control and data points. Solving this kind of system, as Farin’ shows, leads to 0(n4) computations, which are normally 0(n6) if the points are arranged in a vector. Although the algorithms were implemented for the cubic case, they could be easily extended for splines of any order k. Additional data points (ADP) are inserted as degrees of freedom for the fairing only where necessary. This approach proved to be very suitable for the case of measured (or digitized) data when the points have been originally planned to be evenly distributed, but for some practical reason measurements cannot be performed in all of the intervals. An important feature is that the algorithms for ADP insertion and for the actual fairing are independent. So, if for some practical applications another insertion scheme proves to be more appropriate, it can be used without any limitations with the described fairing procedure. The idea of ADP can be extended to other methods for global fair interpolation (approximation). It will improve the performance, especially when a larger amount of data is concerned.
ACKNOWLEDGEMENTS The author would like to thank the British Council and Microsoft who sponsored his work. Figures 6 and 8 are produced from data obtained from a digitization of a Spitting Image puppet of Margaret Thatcher. Thanks to Spitting Image production for providing the puppet and to the UK Advisory Committee on Computer Graphics for funding the digitization. The work described in this paper was done at Queen Mary and Westfield College, University of London. The author is very grateful to his supervisor Mel Slater and to the referees for their help and useful comments.
York (1978) Lee, E ‘Choosing nodes in parametric curve interpolation’ Comput.-Aided Des. Vol21 No 6 (1989) pp 363-370 Hartley, P and Judd, C ‘Parameterization and shape of B-spline curves’ Comput.-Aided Des. Vol 12 No 5 (1980) pp 235-238 Nielson, G and Foley, T ‘A survey of applications of an affine invariant norm’ in Lyche, T. and Shumaker, L. (Eds.) Mathematical Methods in CAGD Academic Press (1989) pp 445-468 Marin, S P ‘An approach to data parameterization in parametric cubic spline interpolation problems’ Approx. Theory Vol41 (1984) pp 64-86 Fritsch, F N and Carlson, R E ‘Monotone piecewise cubic internolation’ SIAM J. Numer. Anal. Vol 17 No 2 (1980) pp 238-246 Kjellander, J A P ‘Smoothing of cubic parametric splines’ Comput.-Aided Des. Vol 15 No 5 (1983) pp 175-179 Rogers. D F. Satterfield. S G and Rodrieuez, F A ‘Shiphull, B-&ink surfaces and CAD/CAM' in Computer Graphics-Theory and Applications Springer-Verlag (1983) Sapidis, N and Farin, G ‘Automatic fairing algorithm for B-spline curves’ Comput.-Aided Des. Vo122 No 2 (1990) pp 121-129 Schweikert, D G ‘An interpolation curve using a spline in tension’ J. Math. Phys. Vol45 (1966) pp 312-317 Salkauskas, K ‘C’ splines for interpolation of rapidly varying data’ Rocky Mountain J. Math. Vol 14 (1984) pp 239-250 Foley, T A ‘A shape preserving interpolant with tension controls’ Comma. Aided Geom. Des. Vol 5 (1988) DD 105-l 18 Celniker, G and Gossard, D ‘Deformable curve and surface finite elements for free-form shape design’ Proc. SIGGRAPHVol25 No 4 (1991) pp 257-266
15 16 17 18
Nowacki, H and Lii, X ‘Fairing composite polynomial curves with constraints’ Comput. Aided Geom. Des. Vol 11 (1994) pp 1-15 Halstead, M, Kass, M and DeRose, T ‘Efficient, fair interpolation using Catmull-Clark surfaces’ Proc. SIGGRAPH, Annual series (1993) pp 35-44 Forsythe, G and Moler, C Computer Solution of Linear Algebraic Systems Prentice-Hall (1967) pp 114-l 19 Press, W H, Flannery, B P, Teukolsky, S A and Vetterling, W T Numerical Recipes, The Art of Scientt$c Computations Cambridge University Press, New York (1989) pp 121-126
APPENDIX Equation 16 is a sum of four similar terms and here the solution to the first one is given. From Equation 17 w:v = BTVTB 2) uBTVB 11 ‘v
(Al)
where the dots mean first derivative. We have to compute d cl w;,dudv = cl 8% d = Clax,,
f f
B;A;~(D,x,
B=V=i B=VB dudv Cl u I4 21
+ Do)=A;=IiUti;fA,’
(D,X, + Ds)A,%,dudv = 2c,
D~A,=i3,li~A,‘(D,X,
+ D,,)
A,‘&,iS;A,=dudv = 2c,D;A;=
(ji.B:du)A;l(D,,Xu
= 2c,D;fA,TK1uA;1(D,X,
+ D,,)A,’
+ Do)A,‘KI,A,=
(-42)
759
Interpolation
and approximation
of B-splines:
T I Vassilev
Results similar to Equation A2 will be obtained evaluating the remaining three terms of Equation 16. Then the system of equations reads (++ means if and only 4
D:A&,KI,
+
LZKd A,’ (DuX,, + Do)A, ’
((Y,K,~ + ,/&K2,)A,T = 0 # D;A,TK,A,‘(D,X, D:K,,(D,X,
+ Do)A,‘K,,A,.T
+ Do)KA,. = 0
= 0 ++ (A3)
Here cl = (Y,,cY,.,cz = ,&,a,,, c3 = c@,,, c4 = /L&. The matrix KAt, is symmetric, band and positive definite and its inverse KA:, exists. Multiplying both sides of
760
Equation
A3 with K,&! results in
D,TKA~D,X, + D,TKAuDo= 0
(Ad)
Tzvetomir Ivanov Vussilev gruduuted jiom the Technicrrl b’niversi!,~ qf Roussr, Bulgaria in 1990. He i.5 currently a lecturer at the Department of Informatics at the same university. His research interests include computer-aided geometric d~@yn, data processing, and computer graphics
