Approximation of minimum energy curves

Applied Mathematics and Computation 108 (2000) 153±166 www.elsevier.nl/locate/amc

Approximation of minimum energy curves Ruibin Qu *, Jieping Ye Department of Mathematics, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260, Singapore

Abstract The problem of interpolating or approximating a given set of data points obtained empirically by measurement frequently arises in a vast number of scienti®c and engineering applications, for example, in the design of airplane bodies, cross sections of ship hull and turbine blades, in signal processing or even in less classical things like ¯ow lines and moving boundaries from chemical processes. All these areas require fast, ecient, stable and ¯exible algorithms for smooth interpolation and approximation to such data. Given a set of empirical data points in a plane, there are quite a few methods to estimate the curve by using only these data points. In this paper, we consider using polynomial least squares approximation, polynomial interpolation, cubic spline interpolation, exponential spline interpolation and interpolatory subdivision algorithms. Through the investigation of a lot of examples, we ®nd a `reasonable good' ®tting curve to the data. Ó 2000 Elsevier Science Inc. All rights reserved. Keywords: Minimal energy curve; Spline; Subdivision algorithms; Approximation

1. Introduction The problem of interpolating or approximating a given set of discrete data points in the plane using smooth ®tting curves is an important area in computer aided geometric design. Sometimes curves which maximize `smoothness' are used. While in engineering, the commonly used and recognized criterion for curve modelling is the least bending energy. The total energy of an elastica of length l is proportional to the integral of the square of the curvature taken along the elastica [1,10,11], i.e.,

*

Corresponding author. E-mail: [email protected]

0096-3003/00/$ - see front matter Ó 2000 Elsevier Science Inc. All rights reserved. PII: S 0 0 9 6 - 3 0 0 3 ( 9 9 ) 0 0 0 1 2 - 0

154

R. Qu, J. Ye / Appl. Math. Comput. 108 (2000) 153±166

Z energy /

l 0

2

j…s† ds;

where j is the curvature and s is the arc length. A curve obtained by minimizing this exact strain energy is called a minimal-energy spline [7±9]. However, the functional expressing the exact energy is complicated, and as a consequence, this method is computationally very demanding and takes a lot of time. In practice, approximation energy functionals are often used. In the subsequent sections we describe two commonly used approximation energy functionals and three other methods to approximate the curves. Eight examples are given to compare these methods. 2. Discrete least squares approximation via orthogonal polynomials and Polynomial interpolation n

Given a set of data points …xi ; yi †0 corresponding to measurements on the test function f 2 C‰a; bŠ, a linear discrete approximation of f minimizes the deviation between f and U with respect to a prescribed error norm. U shall depend on the parameter t 2 ‰a; bŠ and certain coecient set A ˆ fA0 ; . . . ; AN g in R2 , in the form N X Ak /k …t†; U…t; A† ˆ A0 /0 …t† ‡


‡ AN /N …t† ˆ kˆ0

where /k …t† 2 C‰a; bŠ; k ˆ 0 : N are given basis functions. For simplicity, the basis functions /0 , /1 ; . . . ; /N , are orthogonal functions. The uniqueness and existence of the best approximation U is guaranteed by the linearly independence of /0 ; . . . ; /N . The optimum solution U can be obtained by calculating inner products [12]. Polynomial interpolation is the most fundamental of all interpolation conn cepts. Given a set of data points …xi ; yi †0 with the associated parameter values ti 2 ‰a; bŠ, there is a unique polynomial T of degree 6 n ‡ 2 interpolating these points and satisfying the boundary conditions: T 0 …a† ˆ P00 and T 0 …b† ˆ Pn0 . Fast and accurate methods can be found in Ref. [5]. 3. Cubic spline interpolation Let X be a parametric curve in R2 : T

X …t† ˆ …x…t†; y…t†† ; t 2 ‰a; bŠ  R: The curvature of the plane curve is then given by j…t† ˆ

jX 0 …t†  X 00 …t†j jX 0 …t†j3

:

…1†

R. Qu, J. Ye / Appl. Math. Comput. 108 (2000) 153±166

155

The energy integral can be derived (with respect to parameter t) as Z b Z l 2 …x00 y 0 ÿ y 00 x0 † 2 j ds ˆ dt: …2† Eˆ 2 2 5=2 0 a ……x0 † ‡ …y 0 † † For the case where X ˆ …x; y…x††, the exact strain energy of X can be written as Z …y 00 †2 dx: …3† Eˆ 2 5=2 …1 ‡ …y 0 † † 0 So with the assumption R 00 2 jy j  1, Eq. (3) reduces approximately to the energy functional E ˆ …y † dx, which leads to the cubic interpolatory spline, the most popular spline of the lowest degree with C 2 continuity [5,8,12]. De®nition 3.1. Given a set of n ‡ 1 data points Pi ˆ …xi ; yi †, i ˆ 0 : n, in R2 , with associated parameter values ti 2 ‰a; bŠ satisfying a ˆ t0 <


< tn ˆ b, a function S…t† is cubic interpolatory spline function if the following conditions are satis®ed: · S 2 C 2 ‰a; bŠ; · S is a cubic polynomial Si de®ned on the subinterval ‰ti ; ti‡1 Š, i ˆ 0 : n ÿ 1; · S…ti † ˆ Pi for i ˆ 0 : n, i.e., S interpolates the points Pi ; · one of the following boundary conditions: (a) S 00 …t0 † ˆ S 00 …tn † ˆ 0; (b) S 0 …t0 † ˆ P00 , S 0 …tn † ˆ Pn0 ; (c) S…t0 † ˆ S…tn †, S 0 …t0 † ˆ S 0 …tn †, S 00 …t0 † ˆ S 00 …tn †: With appropriate boundary conditions, a linear system of equations can be formed. The coecient matrix is always diagonally dominant. Thus, the linear system is uniquely solvable and the spline interpolation problem always has a unique solution. We end this section by stating the minimum curvature property [16]. Theorem 3.1 (Minimum Curvature Property). Let f …t† 2 C 2 be any function de®ned on ‰a; bŠ, and moreover, let it satisfy the same interpolating and endpoint conditions as the cubic spline S…t†. Then Z b Z b 2 2 00 ‰S …t†Š dt 6 ‰ f 00 …t†Š dt; a

a

with equality holding if and only if f …t†  S…t† over ‰a; bŠ. 4. Exponential spline interpolation Unwanted infection points in an interpolating curve can be removed by shortening the arc length of the curve. This fact leads Schweikert to investigate the extremes of the functional [6,15]:

156

R. Qu, J. Ye / Appl. Math. Comput. 108 (2000) 153±166

Z

b a

00 2

2

…f † dt ‡ r

Z

b a

2

…f 0 † dt:

…4†

The Schwarz inequality gives the following upper bound for the arc length L of a parametrized curve X …t† ˆ …x…t†; y…t††: 1=2  Z b ……x0 …t††2 ‡ …y 0 …t††2 † dt : L 6 …b ÿ a† a

Obviously, the second part of the functional Eq. (4) has a `control' e€ect on the arc length of the curve. The Euler equation of the variational problem de®ned by Eq. (4) yields for each subinterval …4†

fi …t† ÿ r2 fi00 …t† ˆ 0;

t 2 ‰tiÿ1 ; ti Š;

and has a general solution of the form: fi …t† ˆ Ai ‡ Bi t ‡ Ci eÿrt ‡ Di ert :

…5†

Therefore, the extremes of Eq. (4) are also called exponential splines. Among all functions of F2 ˆ ff 2 C 2 ‰a; bŠ j f …ti † ˆ yi ; i ˆ 0; . . . ; n and f 0 …t0 † ˆ y00 ; f 0 …tn † ˆ yn0 g, the spline under tension r with the spline segments given by Eq. (5) minimizes the functional (4) [13]. The idea of exponential splines has been generalized by H. Spath and others by introducing di€erent tension parameters ri …i ˆ 1; . . . ; n† for each subinterval. This approach allows a local smoothing of the curve by varying the tension parameters. Using the following representation of the exponential spline segment fi , fi …t† ˆ ai …t ÿ tiÿ1 † ‡ bi …ti ÿ t† ‡ ci wi …t ÿ tiÿ1 † ‡ di wi …ti ÿ t†;

…6†

where wi …t† ˆ

sinh …ri t† ÿ t sinh …ri hi †=hi ; … sinh …ri hi †=hi † ÿ ri

i ˆ 1; . . . ; n;

…7†

and hi ˆ ti ÿ tiÿ1 , one gets from the interpolation conditions: ai ˆ yi =hi ;

bi ˆ yiÿ1 =hi ;

i ˆ 1; . . . ; n:

…8†

Introducing the quantities vi :ˆ w0i …hi †;

0 ‡ yiÿ1 :ˆ fi0 …tiÿ1 †;

yi0 :ˆ fi0 …tiÿ †;

i ˆ 1; . . . ; n;

…9†

yields ci ˆ

0 yiÿ1 ‡ vi yi0 ÿ …vi ‡ 1†…yi ÿ yiÿ1 †=hi ; v2i ÿ 1

…10†

di ˆ

0 yi0 ‡ vi yiÿ1 ÿ …vi ‡ 1†…yi ÿ yiÿ1 †=hi 1 ÿ v2i

…11†

R. Qu, J. Ye / Appl. Math. Comput. 108 (2000) 153±166

157

0 for i ˆ 1; . . . ; n. The unknowns y20 ; . . . ; ynÿ1 can be determined from the tridiagonal system of equations 0 0 ‡ …ci vi ‡ ci‡1 vi‡1 †yi0 ‡ ci‡1 yi‡1 ci yiÿ1

ˆ ci …vi ‡ 1†…yi ÿ yiÿ1 †=hi ‡ ci‡1 …vi‡1 ‡ 1†…yi‡1 ÿ yi †=hi‡1 ;

…12†

for i ˆ 1; . . . ; n ÿ 1, where ci ˆ

…v2i

r2i sinh …ri hi † : ÿ 1†… sinh …ri hi †=hi ÿ ri †

The coecient matrix of this system is strictly diagonal dominant if ri P 0 for i ˆ 1; . . . ; n. For vanishing tension parameters ri the above exponential spline degenerates to the common cubic spline.

5. The 6-point interpolatory subdivision scheme Subdivision schemes constitute a class of techniques in the area of computer-aided geometric design for obtaining satisfactory piecewise linear approximations to smooth curves and surfaces. In this section, a simple interpolatory scheme based on the 6-point recursive subdivision with a tension parameter is presented. The interpolatory algorithm discussed here is unique in combining the four important ingredients: subdivision, locality, interpolation and shape control. Moreover, like any other subdivision algorithms, it is fast and convenient. It is developed closely with the one using 4points [3]. n‡4 Given control points fPi giˆÿ4 , Pi 2 R2 , intermediate points are successively inserted by the scheme:     9 1 ‡ 2h …Pi ‡ Pi‡1 † ÿ ‡ 3h …Piÿ1 ‡ Pi‡2 † ‡ h…Piÿ2 ‡ Pi‡3 †; Pi‡12 ˆ 16 16 …13† 2n‡4

where ÿ2 6 i 6 n ‡ 1. The augmented set of points fPi giˆÿ4 is then considered as the new control points of the polygon in the next stage. By inserting the new points in®nite times, the scheme (13) de®nes an in®nite set points in R2 . Geometrically, the sequences of control polygons generated by this algorithm can be seen to converge to a continuous curve in R2 . Mathematically, this is a scheme which can produce continuous curves with continuous tangent and curvature, depending on the values of the tension parameter h. For h ˆ 0, the scheme (13) corresponds to the 4-point subdivision scheme with the special value w ˆ 1=16 and for which all polynomials of degree 6 3

158

R. Qu, J. Ye / Appl. Math. Comput. 108 (2000) 153±166

can be reproduced [3] under certain assumption. The range of h that guarantees the continuity of the curvature of the curve is 0 < h < 0:02: Remark 5.1. For general 2r-point schemes of the form f2ik‡1 ˆ fik ; ÿ1 6 i 6 1; m X k‡1 k ˆ wj fi‡j ; ÿ1 6 i 6 1; f2i‡1

…14†

jˆÿm‡1

the values generated by Eq. (14) can be extended to a limiting C r function on …ÿ1; 1† only if the process reproduces polynomials up to degree r:  l m X 1 wj jl ˆ ; l ˆ 0; . . . ; r: …15† 2 jˆÿm‡1 This is the necessary condition for C r continuity of the function de®ned by Eq. (14) [2]. Remark 5.2. For general interpolatory algorithm for curves: f2ik‡1 ˆ fik ; n   X k‡1 k‡1 k‡1 ˆ Ln;j fiÿj ‡ fi‡j‡1 ; f2i‡1

…16†

jˆ0

where n is called the degree of the scheme and fLn;j g are the free parameters to be determined. This symmetric scheme produces C n=2 interpolatory curves provided that fLn;j g are given by the formula Ln;j

‰…2n ‡ 1†!!Š2 …ÿ1†j

ˆ 2
4n
…2n ‡ 1†! 2j ‡ 1



 2n ‡ 1 ; nÿj

j ˆ 0; 1; . . . ; n:

…17†

Furthermore, for this choice of the coecients the scheme reproduces all parametric polynomial curves of degree less than or equal to 2n ‡ 1 [14]. Assuming the tension parameter h is ®xed at 0.01 and by increasing the number of level of the subdivision algorithm, numerical examples show that certain level of the subdivision algorithm is suciently stable in determining the energy values accurately. Then by ®xing at this level, the tension parameter h is changed in order to determine the optimal curve with the least energy. Numerical examples show that the optimum tension parameter h which produces the least energy curve by interpolatory subdivision algorithm can always be determined.

R. Qu, J. Ye / Appl. Math. Comput. 108 (2000) 153±166

159

6. Examples Eight examples are used to examine the methods described in the prior sections. The quality of di€erent methods is illustrated with the exact strained energies and lengths of the curves produced. The exact strained energy is obtained through Eq. (2). We use the Composite Simpson's rule [4] with equally spaced parameter values. A similar way is used to compute the lengths. Generally, curves with less strained energy and normal length look `nice'. In the examples, the data point marked with an `o' is sampled from the test functions. In each example, we have summarized the numerical results and rank di€erent methods according to the energy. Tables 1±8 show the results. Table 1 Rank

Methods

Length

Energy

1 2 3 4 5

Least square approximation Interpolatory subdivision (h ˆ 0:014) Exponential spline interpolation Cubic spline interpolation Polynomial interpolation Reference (estimated)

6.278373 6.284724 6.279369 6.279287 6.283596 6.283000

6.139000 6.282587 6.290680 6.291333 6.299659 6.283000

Methods

Length

Energy

Interpolatory subdivision (h ˆ 0:027) Exponential spline interpolation Cubic spline interpolation Least square approximation Polynomial interpolation Reference (estimated)

2.746439 2.736654 2.736536 2.780524 3.367553 2.740000

15.577460 16.864245 16.876663 33.848561 157.93847 16.190000

Rank

Methods

Length

Energy

1 2 3 4 5

Least square approximation Interpolatory subdivision (h ˆ 0:015) Exponential spline interpolation Cubic spline interpolation Polynomial interpolation Reference (estimated)

34.977903 35.013768 34.934167 34.930996 34.996828 34.990000

7.524365 7.727407 7.757023 7.758500 8.106722 7.730000

Table 2 Rank 1 2 3 4 5



Table 3

160

R. Qu, J. Ye / Appl. Math. Comput. 108 (2000) 153±166

Table 4 Rank

Methods

Length

Energy

1 2 3 4 5

Exponential spline interpolation Cubic spline interpolation Interpolatory subdivision (h ˆ 0:013) Least square approximation Polynomial interpolation Reference (estimated)

2.261849 2.261812 2.261498 2.278442 2.308188 2.260000

4.061943 4.070289 4.107883 9.791447 19.53611 3.980000

Rank

Methods

Length

Energy

1 2 3 4 5

Cubic spline interpolation Exponential spline interpolation Interpolatory subdivision (h ˆ 0:017) Polynomial interpolation Least square approximation Reference (estimated)

13.755117 13.756338 13.777447 15.897343 17.123659 13.970000

68.818295 68.835072 69.877990 116.20545 249.32671 69.880000

Rank

Methods

Length

Energy

1 2 3 4 5

Exponential spline interpolation Interpolatory subdivision (h ˆ 0:019) Least square approximation Polynomial interpolation Cubic spline interpolation Reference (estimated)

9.440642 9.449382 9.418153 9.428613 9.421462 9.430000

42.469459 44.552279 44.700214 44.860168 45.191475 44.840000

Methods

Length

Energy

Interpolatory Subdivision (h ˆ 0:018) Least Square Approximation Exponential Spline Interpolation Cubic Spline Interpolation Polynomial Interpolation Reference (estimated)

3.199885 3.187952 3.185101 3.184821 3.273062 3.190000

197.395353 197.604170 200.768101 200.899386 244.559839 198.510000

Methods

Length

Energy

4.180119 4.182243 4.187301 5.923992 34.54480 4.180000

8.777637 8.785522 8.908945 293.3368 505.5519 8.600000

Table 5

Table 6

Table 7 Rank 1 2 3 4 5



Table 8 Rank 1 2 3 4 5



Interpolatory subdivision (h ˆ 0:013) Cubic spline interpolation Exponential spline interpolation Least square approximation Polynomial interpolation Reference (estimated)

R. Qu, J. Ye / Appl. Math. Comput. 108 (2000) 153±166

161

The graphical plots using interpolatory subdivision algorithm are also shown as Figs. 1±8. All the curves are generated by Mathematica. The corresponding lengths and energy of the constructed curve does not vary signi®cantly from the reference length and energy (cf. Remark 6.1).

Fig. 1. Our ®rst example data set is that of 9 points fPi g8iˆ0 sampled from the circle: x ˆ cos…t†, y ˆ sin…t†, for t 2 ‰0; 2pŠ. Here Pi ˆ …x…ti †; y…ti †† …i ˆ 0; . . . ; 8†, where ti ˆ ih and h ˆ p=4.

Fig. 2. 9 points fPi g8iˆ0 from the function y ˆ x2 =…0:16 ‡ x2 † on ‰ÿ1; 1Š are used. Here Pi ˆ …xi ; y…xi †† …i ˆ 0; . . . ; 8†, where xi ˆ ÿ1 ‡ ih and h ˆ 1=4.

162

R. Qu, J. Ye / Appl. Math. Comput. 108 (2000) 153±166

t=10 Fig. 3. 16 points fPi g15 sin …t†, y ˆ et=10 cos…t†, for t 2 ‰0; 15Š are used. iˆ0 from the curve: x ˆ e Here Pi ˆ …x…ti †; y…ti †† …i ˆ 0; . . . ; 15†, where ti ˆ ih and h ˆ 1.

Fig. 4. 8 points fPi g7iˆ0 from the function: y ˆ 1=…e5x ‡ eÿx † on ‰ÿ1; 1Š are used. Here Pi ˆ …xi ; y…xi †† …i ˆ 0; . . . ; 7†, where xi ˆ ÿ1 ‡ ih and h ˆ 2=7.

Remark 6.1. By inserting more points into the original data points, the curves should converge to that of the intended curve (or test function) with the minimal energy. The estimated length and energy of this curve can be computed precisely and used as reference values.

R. Qu, J. Ye / Appl. Math. Comput. 108 (2000) 153±166

163

Fig. 5. 12 points fPi g11 iˆ0 from the function: y ˆ sin…3x† on ‰0; 2pŠ are used. Here Pi ˆ …xi ; y…xi †† …i ˆ 0; . . . ; 11†, where xi ˆ ih and h ˆ 2p=11.

Fig. 6. 15 points fPi g14 iˆ0 from the curve: x ˆ sin …t=2†, y ˆ sin…t†, for t 2 ‰0; 4pŠ are used. Here Pi ˆ …x…ti †; y…ti †† …i ˆ 0; . . . ; 14†, where ti ˆ ih and h ˆ 2p=7.

7. Analysis of experiment results The results are summarized by ranking the curve-®tting methods according to the values of energy for di€erent data sets in Table 9. The following analyses are derived from our investigations:

164

R. Qu, J. Ye / Appl. Math. Comput. 108 (2000) 153±166

Fig. 7. 19 points fPi g18 iˆ0 from the curve: x ˆ …1=4p†t ‡ …1=2p† sin…t†, y ˆ …1=4p† ‡ …1=2p† cos…t†, for t 2 ‰0; 6pŠ are used. Here Pi ˆ …x…ti †; y…ti †† …i ˆ 0; . . . ; 18†, where ti ˆ ih and h ˆ p=3.

5 6 Fig. 8. 16 points fPi g15 iˆ0 from the function: y ˆ x =…0:16 ‡ 4x † on ‰ÿ2; 2Š are used. Here Pi ˆ …xi ; y…xi †† …i ˆ 0; . . . ; 15†, where xi ˆ ih and h ˆ 4=15.

· Table 9 clearly shows that the 6-point interpolatory subdivision scheme and exponential spline interpolation are rated as the best two curve-®tting methods. Moreover, the subdivision scheme is fast and convenient to apply. Surveying the numerical results shows that the energy of the curves generated by using subdivision scheme are often very close to the least energy value when it is not the best. It also generates smooth and well-de®ned curves that ®t the reference curves with a similar shape.

R. Qu, J. Ye / Appl. Math. Comput. 108 (2000) 153±166

165

Table 9 Method\example\rank

1

2

3

4

5

6

7

8

1st 2nd 3rd 4th 5th

Least squares approximation Polynomial interpolation Cubic spline interpolation Exponential spline interpolation 6-point interpolatory subdivision

1 5 4 3 2

4 5 3 2 1

1 5 4 3 2

4 5 2 1 3

5 4 1 2 3

3 4 5 1 2

2 5 4 3 1

4 5 2 3 1

2 0 1 2 3

1 0 2 2 3

1 0 1 4 2

3 2 3 0 0

1 6 1 0 0

· The cubic spline interpolation is rated the best only once. But one observes from experimental results that even at its worst ranking, the energy values are often reasonably close to the least energy value. It suces to say that this method typically produces satisfactory data ®tting curves. · For the least squares approximation, on occasions when it is rated fourth or last, the energy and length of the curves constructed often deviate substantially from that of the minimal-energy curve. It highlights the fact that the curves have many unnecessary in¯ections and bending. Thus, it is potentially bad as a curve ®tting method for constructing the smoothest curve in the least-energy sense. · Evidently, from the numerical results the polynomial interpolation is the worst curve ®tting method. 8. Conclusion Five methods have been presented for constructing smooth curves interpolating given empirical data points. From numerical results, it shows that interpolatory subdivision scheme and exponential spline interpolation are wellsuited to this modelling problem. They yield curves of high quality in a reasonable amount of time. But practical application shows that it is not an easy work to manipulate the tension parameters ri , especially when the number of data points increases. In conclusion, using the interpolatory subdivision scheme, one makes a good compromise between curve quality and algorithmic eciency. References [1] J.P. Den Hartog, Strength of Materials. Dover, New York, 1960. [2] N. Dyn, D. Levin, Smooth Interpolation by Bisection Algorithms, in: Chui, Schumaker, Wards (Eds.), Approximation Theory 5, 1986, pp. 335±337. [3] N. Dyn, J. Gregory, D. Levin, A 4-point interpolatory subdivision scheme for curve design, Comput. Aided Geom. Des. 4 (1987) 257±268. [4] G. Engeln-mullges, F. Uhlig, Numerical algorithms with C. Springer, Berlin, 1996.

166

R. Qu, J. Ye / Appl. Math. Comput. 108 (2000) 153±166

[5] G. Farin, Curves and Surfaces for Computer Aided Geometric Design: A practical Guide, 3rd ed., Academic Press, New York, 1993. [6] H. Hagen, G. Schulze, Variational principals in Curve and Surface Design, in: H. Hagen, D. Roller (Eds.), Geometric Modelling. Methods and Applications, Springer, Berlin, 1991, pp. 161±184. [7] B.K.P. Horn, The curve of least energy, ACM Trans. Math. Softw. 9 (1983) 441±460. [8] J. Hoschek, D. Lasser, Fundamentals of Computer Aided Geometric Design, A.K. Peters, Wellesley, MA, 1993. [9] E. Jou, W. Han, Minimal-energy splines: I. Plane curves with angle constrains, Math. Meth. Appl. Sci. 13 (1990) 351±372. [10] E. Jou, W. Han, Elastic and Minimal-energy Splines in Curves and Surfaces (chamonix± Mont±Blanc, 1990), Academic Press, New York, 1991, pp. 247±250. [11] E.H. Lee, G.E. Forsythe, Variational study of nonlinear spline curves, SIAM Rev. 15 (1973) 161±184. [12] W.C. Lim, Investigation of optimum curve ®tting with respect to approximation energy norms, Honours Project in Computational Science and Mathematics, 1997/1998. National University of Singapore. [13] D.T. Pilcher, Smooth parametric surfaces, in: Barnhill-Riesenfeld (Ed.), Computer Aided Geometric Design, Academic Press, New York, 1974. [14] R. Qu, R.P. Agarwal, A cross di€erence approach to the analysis of subdivision algorithms, Neural, Parallel & Scienti®c Computations 3 (3) (1995) 393±416. [15] D.G. Schweikert, An interpolation curve using a spline in tension, J. Math. Phys. 45 (1966) 312±317. [16] A.F. Thomas, Interpolation with interval and point tensior controls using cubic weighted v-splines, ACM Trans. Math. Softw. 13 (1987) 68±96.