- Email: [email protected]
An iterative algorithm for polynomial approximation of rational triangular Bézier surfaces
Applied Mathematics and Computation 219 (2013) 9308–9316
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
An iterative algorithm for polynomial approximation of rational triangular Bézier surfaces Qianqian Hu Department of Mathematics, Zhejiang Gongshang University, Hangzhou 310018, PR China
a r t i c l e
i n f o
Keywords: Rational triangular Bézier surfaces Polynomial approximation Weighted progressive iteration Triangular Bernstein basis Gauss Legendre quadrature rule
a b s t r a c t This paper presents the weighted progressive iteration approximation (WPIA) property for the triangular Bernstein basis over a triangle domain with uniform parameters, which is extended from the PIA property for triangular Bernstein basis proposed by Chen and Wang in [J. Chen, G.J. Wang, Progressive-iterative approximation for triangular Bézier surfaces, Computer-Aided Design 43 (2011) 889–895]. We also provide how to choose an optimal value of the weight to own the fastest convergence rate for triangular Bernstein basis. Furthermore, a new and efficient iterative method is proposed for polynomial approximation of rational triangular Bézier surfaces. The algorithm is reiterated until a halting condition about approximation error is satisfied. And the approximation error in Lp-norm (p = 1, 2, 1) is calculated by the symmetric Gauss Legendre quadrature rule for composite numerical integration over a triangular surface. Finally, several numerical examples are presented to validate the effectiveness of this method. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction Approximating rational curves and surfaces by polynomial curves and surfaces is a significant problem in Computer Aided Geometric Design (CAGD). There are two advantages for polynomial approximation of rational forms. One is to simplify the integral and differential calculations of rational curves and surfaces; another is to exchange data between different curve and surface design systems using polynomial or rational forms. In recent twenty years, many works on polynomial approximation of rational Bézier curves have been developed. In 1991, Sederberg and Kakimoto [1] presented the hybrid algorithm for polynomial approximation of rational curves for the first time. Then the recursive algorithm and convergence condition for the hybrid approximation were put forward in [2]. However, the hybrid approximation involves complicated calculation and its convergence is not guaranteed. Floater [3] showed that a k-degree rational parametric curve can be interpolated by an m-degree polynomial matching 2m 2k + 4 data. Huang et al. [4] constructed a Bézier curve sequence with control points obtained by degree-elevate the control points of a rational Bézier curve. This method is simple, and can guarantee that the approximation curve sequence converges to the rational curve, but cannot approximate the rational curve by polynomial curves with lower degree. The progressive iteration approximation (PIA) property for polynomial curves was first proposed in [5,6], and attracted a good deal of attention [7–10]. Curves generated by normalized totally positive (NTP) basis, such as Bernstein basis and B-spline basis, satisfy the PIA property [11]. Recently, Lu [12] presented an iterative method for polynomial approximation of rational Bézier curves based on the idea of weighted progressive iteration approximation proposed in [13]. As we all know, it is easy to generalize the abovementioned methods to polynomial approximation of rational tensor product surfaces, because a rational tensor product Bézier surface can be regarded as the space moving trace of a Bézier E-mail addresses: [email protected], [email protected] 0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.03.053
Q. Hu / Applied Mathematics and Computation 219 (2013) 9308–9316
9309
curve with moving control points. Compared with surfaces in rectangular domain, Bernstein–Bézier surfaces over a triangle domain possess a series of beautiful geometric properties, such as be flexible for blending and the non-degenerate representation of segments of a sphere [14]. Also triangular surfaces are adapted to surface interpolation from scattered data and surface design on arbitrary topology mesh [15]. Thus study on the polynomial approximation of rational triangular Bézier surfaces is an urgent need for engineering. Up until now, relatively little has been published about this topic, because the expression form of rational triangular Bézier surfaces are more complex than rational Bézier curves. Enlightened by the degree-elevated vertices approximation [4], Xu and Wang [16] constructed a polynomial triangular Bézier approximation surface, which has the same control points and degree of the rational triangular Bézier surface after degree elevation. Zhang and Wang [17] extended the hybrid approximation to rational triangular Bézier surfaces. Inspired by [18], we present a new and efficient iterative method, the weighted progressive iteration approximation, for polynomial approximation of rational triangular Bézier surfaces. Our method differs from the previous technique [18] in that we multiply all the adjusting control points by some weight to accelerate the iteration convergence rate. And then we apply the WPIA method to approximate rational triangular Bézier surfaces by polynomial triangular surfaces. First of all, we sample some points from the original rational triangular Bézier surface with uniform parameters. Second, we construct a triangular Bézier surface with these sample points as control points and take it as an initial iteration surface. Finally, we use the WPIA method to reiterate the initial surface to approximate the original rational surface. The iterative process is repeated until a satisfactory approximation error is reached. The approximation error in the Lp (p = 1, 2, 1)-norm is calculated with high accuracy by the symmetric Gauss Legendre quadrature rule for composite numerical integration over a triangular surface [19]. The paper is laid out as follows. We give a brief introduction to rational triangular Bézier surfaces in Section 2. In Section 3, we propose the weighted progressive iteration approximation for triangular Bézier surfaces and proof its convergence property. Then we present a numerical calculation of the approximation error between the original rational surface and the approximation surface in different norms, and provide an algorithm for polynomial approximation of rational triangular Bézier surfaces using the weighted progressive iteration approximation in Section 4. Finally, numerical examples are presented in Section 5 to confirm the effectiveness of this algorithm. 2. Preliminary For the following developments, we need first some notations: i ¼ ði; j; kÞ; jij ¼ i þ j þ k. When we say jij ¼ m, we mean i þ j þ k ¼ m, always assuming i; j; k P 0. Definition 1. A rational triangular Bézier surface of degree m over a triangle domain T :¼ fðu; v ; wÞ : u; v ; w P 0; u þ v þ w ¼ 1g is represented by [15]
P
m
xi bi Bmi ðuÞ ; m jij¼m xi Bi ðuÞ
jij¼m
b ðuÞ ¼ P
ð1Þ
where
Bm i ðuÞ ¼
m
i
ui ¼
m! i j k uv w ; i!j!k!
jij ¼ m;
juj ¼ 1;
are the triangular Bernstein basis functions of degree m, bi ¼ bi;j;k are the control points, and xi ¼ xi;j;k are the associated weights. Here are ðmþ1Þðmþ2Þ basis functions. 2 Definition 2. Suppose T; bi ; xi are shown as in Definition 1. Then a rational triangular Bézier surface equivalent to the surface (1) is represented by m
b ðu; v Þ ¼
Pm
m iþj¼0 i;j;mij bi;j;mij Bi;j;mij ðu; Pm m iþj¼0 i;j;mij Bi;j;mij ðu; ; 1
x
x
v
v; 1 u vÞ ; u vÞ
ð2Þ
and is called an m degree rational triangular Bézier surface over a triangle domain D :¼ fðu; v Þ : 0 6 u; v 6 1; u þ v 6 1g Definition 3. Lexicographic order [20]: Given two vectors a; b in dimension d P 1. The vector a is arranged before the vector b, denoted as a b, if the first nonzero entry in the difference a b ¼ ða1 b1 ; . . . ; ad bd Þ is positive. Definition 4. The index set of any positive integer m is defined by Km ¼ fðm; 0; 0Þ; ðm 1; 1; 0Þ; ðm 1; 0; 1Þ; ðm 2; 2; 0Þ; ðm 2; 1; 1Þ; ðm 2; 0; 2Þ; . . . ; ð0; m; 0Þ; ð0; m 1; 1Þ; ð0; m 2; 2Þ; . . . ; ð0; 2; m 2Þ; ð0; 1; m 1Þ; ð0; 0; mÞg. The index set has
ðmþ1Þðmþ2Þ 2
elements sorted in lexicographic order.
9310
Q. Hu / Applied Mathematics and Computation 219 (2013) 9308–9316
Example 1. The control points of a rational triangular Bézier surface of degree 5 are arranged according to the subscripts in 5 lexicographic order as b ¼ fbi gi2K5 ¼ fb500 ; b410 ; b401 ; b320 ; b311 ; b302 ; b230 ; b221 ; b212 ; b203 ; b140 ; b131 ; b122 ; b113 ; b104 ; b050 ; b041 ; b032 ; b023 ; b014 ; b005 g. 3. Weighted progressive iteration approximation and convergence analysis At first, we proof that triangular Bernstein basis over a triangle domain with uniform parameters satisfy the weighted progressive iteration approximation property. Given a sequence of points fpi gi2Kn arranged by the subscripts in lexicographic order in R3 and the corresponding uniform parameter sequence in lexicographic order fui gi2Kn ; ui ¼ ni ; jij ¼ n, then the initial approximation triangular surface of degree n over the triangle domain T is represented by
p0 ðuÞ ¼
X p0i Bni ðuÞ; jij¼n
where p0i ¼ pi : For each point pi , we calculate the adjusting vector
D0i ¼ pi p0 ðui Þ;
jij ¼ n;
and let
p1i ¼ p0i þ kD0i ;
jij ¼ n;
where k is a weight. When k is chosen to some special value, the iteration approximation has the fastest convergence rate. Theorem 1 in the following will show readers how to choose the optimal value of k. Then the first iteration triangular surface over T is expressed as
p1 ðuÞ ¼
X p1i Bni ðuÞ: jij¼n
After the (s + 1)-th iteration, the triangular surface over T can be constructed as follows
psþ1 ðuÞ ¼
X psþ1 Bni ðuÞ; i jij¼n
where
pisþ1 ¼ psi þ kDsi ;
Dsi ¼ pi ps ðui Þ;
jij ¼ n:
ð3Þ
Then we construct a triangular Bézier surface sequence ps ðuÞ; s ¼ 0; 1; . . . According to (3), we obtain
Dsi ¼ Ds1 k i
X s1 n Dj Bj ðui Þ;
jij ¼ n:
jjj¼n
In matrix–vector notation, the abovementioned equation can be rewritten as
h iT s T Di i2Kn ¼ ðI kBÞ Ds1 i
n
i2K
h iT ¼ ðI kBÞs D0i
i2Kn
;
ð4Þ
T where the elements of the column vector Dsi i2Kn are arranged according to the subscript i in lexicographic order, I is the identity matrix of order ðnþ1Þðnþ2Þ , and B ¼ Bni nj is the collocation matrix of the triangular Bernstein basis functions 2 i2Kn ;j2Kn n
i Bi i2Kn at the uniform parametric values n i2Kn , expressed as
1 Bnn;0;0 nn ; 0n ; 0n Bnn1;1;0 nn ; 0n ; 0n Bn0;0;n nn ; 0n ; 0n
C B n n1 1 0 B Bn;0;0 n ; n ; n Bnn1;1;0 n1 ; 1n ; 0n Bn0;0;n n1 ; 1n ; 0n C n n C B ¼B C: .. .. .. .. C B . . . . A @
Bnn1;1;0 0n ; 0n ; nn Bn0;0;n 0n ; 0n ; nn Bnn;0;0 0n ; 0n ; nn 0
B :¼
n Bi n
i i2K n i2Kn
!
ð5Þ
Thus, for the triangular Bézier surface sequence ps ðuÞ; s ¼ 0; 1; . . ., if
limps ðui Þ ¼ pi ;
s!1
ui ¼
i ; n
jij ¼ n;
then we call the initial triangular surface p0 ðuÞ has the weighted progressive iteration approximation property (ab. WPIA).
9311
Q. Hu / Applied Mathematics and Computation 219 (2013) 9308–9316
According to [18,21], we know that all the eigenvalues ni of the collocation matrix B (5) satisfy that 0 < ni 6 1; jij ¼ n. We denote by n ¼ minfni g its minimum, and denote by qðBÞ its spectral radius, which is the maximal absolute value of its eigenjij¼n
values. Then the iterative process (4) converges when qðI kBÞ < 1, and its speed depends on the value of qðI kBÞ. The smaller the value is, the faster the iterative algorithm converges. Obviously, when k ¼ 1, the WPIA property turns out to be the progressive iteration approximation for triangular Bézier surfaces [18]. Next, the following theorem will show us how to choose the optimal value of k so that the iteration approximation has the fastest convergence rate. Theorem 1. Let Bni i2Kn be the triangular Bernstein basis functions, then weighted progressive iteration approximation
based on the uniform parameters ni i2Kn has the fastest convergence rate when
k¼
2nn ; nn þ n!
ð6Þ
and
qðI kBÞ ¼
nn n! nn þ n!
Proof.n According
to
[18,20],o the
eigenvalues
n! Xn ¼ nni ¼ ni ðniÞ! j i ¼ 0; 1; . . . ; n .
ni
of
the
collocation
matrix
B
(5)
belong
to
the
set
It is easy to verify that
nniþ1 n i ¼ < 1; n nni
iP1
and the number sequence nni i¼1;...;n is strictly decreasing, and nni 2 ð0; 1. So the maximum eigenvalue of the matrix B is 1, and the minimum eigenvalue of the matrix B is
n¼
n! nn
ð7Þ
Note that
nni ¼
n! n n1 niþ1 ¼ 6 1; ni ðn iÞ! n n n
i.e., 0 < nni 6 1; i ¼ 0; 1; . . . ; n, so it holds 0 < ni 6 1 i 2 Kn , and the eigenvalues of the matrix I B are equal to n 1 ni 2 ½0; 1Þ i 2 K . This result implies
qðI BÞ ¼ 1 n < 1;
ð8Þ
where n is the minimum eigenvalue of the matrix B defined by (7). First, according to the definition of spectral radius, we have
qðI kBÞ ¼ maxfj1 knj; j1 kjg
ð9Þ
Obviously, qðI kBÞ P j1 kj. Due to the convergence of the iterative process, qðI kBÞ < 1 must hold, i.e., 0 < k < 2. Second, for any k 2 ð0; 1Þ, according to (8) and (9), it holds
qðI kBÞ ¼ j1 knj ¼ 1 kn > 1 n ¼ qðI BÞ: Therefore, when k ¼ 1, the convergence rate is always faster than that of k 2 ð0; 1Þ. Finally, for any k 2 ð1; 2Þ, in the case of kn > 1, according to (8) and (9), it holds
qðI kBÞ ¼ k 1 > n1 1 > 1 n ¼ qðI BÞ: It means that when k ¼ 1, the convergence rate is always faster than that in such case. On the other hand, in the case of kn 6 1, we have
(
qðI kBÞ ¼ maxf1 kn; k 1g ¼
1 n; 1 < k < 2ð1 þ nÞ1 k 1; 2ð1 þ nÞ1 6 k < 2
:
Clearly, qðI kBÞ reaches the minimum when k ¼ 2ð1 þ nÞ1 . In other words, the weighted progressive iteration approximation based on the uniform parameters has the fastest convergence rate when
k¼
2 ; 1þn
and in such case, the spectral radius of the matrix ðI kBÞ is
9312
Q. Hu / Applied Mathematics and Computation 219 (2013) 9308–9316
qðI kBÞ ¼
1n : 1þn
Substituting (7) into the two abovementioned equations, the theorem is proved. h
4. Error measurement and algorithm description m
We consider the problem of approximating the rational triangular Bézier surface b ðuÞ (1) by an n degree triangular Bézier surface defined by
pðuÞ ¼
X n pi Bi ðuÞ:
ð10Þ
jij¼n
using the weighted progressive iteration approximation property. Since we have only proved that triangular Bézier surfaces with uniform parameters have the WPIA property, we sample m ðnþ1Þðnþ2Þ points p0i i 2 Kn from the original surface b ðuÞ with the uniform parameters ui ¼ ni ¼ ni ; nj ; nk ; ði 2 Kn Þ, i.e., 2 m
p0i ¼ b ðui Þ;
ui ¼
i ; n
jij ¼ n:
ð11Þ
for the iterative process. It should be emphasized that we are interested in a good approximation to the rational surface m m b ðuÞ rather than the sampled points b ðui Þði 2 Kn Þ. Next, we provide the error measurement for the approximation surface. The approximation error in the Lp-norm is defined by
Fig. 1. (a) The rational triangular Bézier surface of degree 6; (b) The approximation triangular surface after 6 iterations; (c) The original surface, the approximation surfaces after 0 and 6 iterations superimposed; (d) The error plot for the approximation surfaces after different iterations.
Q. Hu / Applied Mathematics and Computation 219 (2013) 9308–9316
ep ¼
Z Z
m jb ðuÞ pðuÞ jp dA
1p ;
1 6 p < 1;
9313
ð12Þ
T m
where b ðuÞ is the original rational surface (1), and A is the area of the triangle domain T, and when p = 1, the error is defined by
e1 ¼ maxkbm ðuÞ pðuÞk: u2T
ð13Þ
For the double integral over the triangle domain T in (12), it is very hard to integrate it explicitly. So we consider using m m numerical integration methods for good results. Since the surface b ðuÞ (1) is equivalent to the surface b ðu; v Þ (2), we can change the triangle domain T :¼ fðu; v ; wÞ : u; v ; w P 0; u þ v þ w ¼ 1g as D :¼ fðu; v Þ : 0 6 u; v 6 1; u þ v 6 1g. Then we denote by pðu; v Þ the approximation surface pðuÞ defined over D. Clearly, the integral region D is the standard triangular surface in the Cartesian two dimensional (u, v) space, and (12) can be revised as
ep ¼
Z
1
dv
Z
0
1v
1p f ðu; v Þdu ;
ð14Þ
0 m
where f ðu; v Þ ¼ kb ðu; v Þ pðu; v Þkp . To yield an accurate and reliable evaluation of the integral (14), a symmetric Gauss Legendre quadrature rule is provided for composite numerical integration over the standard triangular surface D [19]. We discretize the triangle domain D in (u, v) space into N N = N2 right isosceles triangle Di each of area 1/(2N2), and divide each triangle Di into three new triangles of equal area 1/(6N2) by joining the centroidal point to the three vertices of Di. Readers can use the q-order of Gauss Legendre quadrature rule for numerical integration, and refer to [19,22] for more details. As Rathod et al. said in [19], the numerical integration algorithm converges to the exact value of the integral, for sufficiently large value of N and the convergence is
Fig. 2. (a) The rational triangular Bézier surface of degree 5; (b) The original rational surface and the approximation surface using our method superimposed; (c) The original surface and the approximation surface using Xu&Wang’s method superimposed; (d) The error plot for approximation surfaces obtained by two different methods.
9314
Q. Hu / Applied Mathematics and Computation 219 (2013) 9308–9316
Table 1 The Lp-error (p = 1, 2, 1) between the k-th iteration triangular surface and the original rational surface. k
L1
L2
L1
k
L1
L2
L1
0
2.652368e-001 1.627374e-001 1.171966e-001 8.241156e-002 7.545272e-002 6.758371e-002 6.564027e-002 6.215079e-002 6.017676e-002 5.781102e-002 5.593812e-002 4.905877e-002 4.595881e-002 4.535136e-002 4.559194e-002
3.090208e-001 1.836061e-001 1.302147e-001 8.674335e-002 8.375550e-002 7.258315e-002 7.180955e-002 6.739415e-002 6.546666e-002 6.271153e-002 6.071003e-002 5.352006e-002 5.097849e-002 5.064423e-002 5.109564e-002
7.765044e-001 4.574360e-001 3.384068e-001 2.755662e-001 2.686409e-001 2.774455e-001 2.745880e-001 2.763090e-001 2.745905e-001 2.745392e-001 2.736165e-001 2.735843e-001 2.762940e-001 2.795807e-001 2.824704e-001
35 40 45 50 55 60 65 70 75 80 85 90 95 100 105
4.611488e-002 4.664660e-002 4.706643e-002 4.738014e-002 4.760596e-002 4.776406e-002 4.787312e-002 4.794862e-002 4.800051e-002 4.803631e-002 4.806085e-002 4.807760e-002 4.808906e-002 4.809688e-002 4.810221e-002
5.171207e-002 5.227033e-002 5.271267e-002 5.304198e-002 5.327903e-002 5.344630e-002 5.356289e-002 5.364353e-002 5.369901e-002 5.373706e-002 5.376310e-002 5.378088e-002 5.379302e-002 5.380130e-002 5.380694e-002
2.847315e-001 2.864040e-001 2.876034e-001 2.884479e-001 2.890356e-001 2.894418e-001 2.897210e-001 2.899125e-001 2.900434e-001 2.901329e-001 2.901939e-001 2.902355e-001 2.902639e-001 2.902832e-001 2.902964e-001
1 2 3 4 5 6 7 8 9 10 15 20 25 30
much faster for higher order (q-order) rules. To achieve high accuracy of numerical integration, it is desirable that N and q are as large as possible. However, since the computation time increases as N and q increase, we need to seek a balance between the speed and the accuracy. Experimentally, we found that q = 2, N = 10 is suitable for good results. When p = 1, (13) can be estimated by
i N
m e1 ¼ max kb ðui Þ pðui Þk; ui ¼ : N i2K
ð15Þ
After our testing, we found that the estimation works well when 50 6 N 6 100. After the error measurement, we present the iterative algorithm. First, we sample ðnþ1Þðnþ2Þ points p0i ði 2 Kn Þ (11) from the 2 m rational triangular surface b ðuÞ for the iterative process. The initial approximation surface is expressed by P p0 ðuÞ ¼ jij¼n p0i Bni ðuÞ; and clearly has the WPIA property. Thus we can construct a triangular Bézier surface sequence m s p ðuÞ; s ¼ 0; 1; :, and the limit surface interpolates the ðnþ1Þðnþ2Þ points b ðui Þði 2 Kn Þ 2 m Clearly, the approximation surface pðuÞ (10) interpolates three corner points of the original rational surface b ðuÞ (1), since we have un;0;0 ¼ ð1; 0; 0Þ; u0;n;0 ¼ ð0; 1; 0Þ; and u0;0;n ¼ ð0; 0; 1Þ. n Second, according to Theorem 1, we choose the optimal value of weight as k ¼ n2n n þn! for the iteration approximation with the fastest convergence rate. It should be pointed out that the iterative method converges slowly, even for the best choice of the parameter k, compared with Xu&Wang’s method, i.e., the degree-elevated method. However, it has a better approximation result than Xu&Wang’s method. Finally, in order to find a good approximation surface and improve the iteration speed, the iterative algorithm can be reiterated until a halting condition about the approximation error is satisfied. The termination condition iskþ1 the relative differje ek j ence between the two consecutive iteration errors ek and ekþ1 is less than a given threshold value f, i.e., ek 6 f; f 2 ð0; 1Þ; or the latter iteration error ekþ1 is larger than the former error ek as ekþ1 > ek . So the termination condition is ekþ1 P dek ; d 2 ð0; 1Þ: We summarize the iterative algorithm as follows. Algorithm 1 (Polynomial approximation of a rational triangular Bézier surface). m Input: An m degree rational triangular Bézier surface b ðuÞ with the control points fbi gi2Km and weights fxi gi2Km , a userspecified error thresholding d, and p in Lp-error. Output: An n degree triangular Bézier surface pðuÞ with the control points fpi gi2Kn , the iteration number s, and the iteration error e. Step Step Step Step Step
1. 2. 3. 4. 5.
Let k ¼ 0. When p ¼ 1; 2, let q ¼ 2; N ¼ 10; when p ¼ 1, let N ¼ 100: Compute the weight k by (6). Compute p0i i2Kn of the initial iteration by (11), and the initial iteration error e0p by (14) or (15). Computer psþ1 of the ðs þ 1Þ-th iteration by (3) and the iteration error esþ1 by (14) or (15). p i i2Kn If esþ1 6 des ; then set s ¼ s þ 1 and go to Step 4. Else, set pi ¼ psi ; i 2 Kn and e ¼ esp . Then output fpi gi2Kn ; e; s; and stop.
5. Numerical examples
Example 2. Given a rational triangular Bézier surface of degree 6, its control points in lexicographic order are (2, 5, 1), (2, 4, 2), (3, 4.2, 2), (1.5, 3.5, 3), (2.5, 3, 5), (3.5, 3.5, 4), (1, 2, 3), (1.5, 2.8, 4), (2.5, 2.5, 5), (4, 3, 6), (0.8, 1.5, 1.5), (2, 1.2, 2), (3, 2, 3), (4, 2.2, 4), (4.5, 2, 5), (0.4, 0.8, 0.6), (1.2, 1, 2), (2, 0.9, 3), (3, 0.6, 4), (4.5, 1, 3), (5.2, 1, 3), (0, 0, 0), (0.4, 0.4, 1), (1.5, 0.5, 2), (3, 0.2,
9315
Q. Hu / Applied Mathematics and Computation 219 (2013) 9308–9316 Table 2 The L1-errors of two different approximations by triangular Bézier surfaces of degree n. n
Our algorithm (d ¼ 0:9)
3 4 5 6 7 8 9 10 11 12 13 14 15
Our algorithm (d ¼ 0:95)
Xu&Wang’s method
Iteration number
L1-error
Iteration number
L1-error
L1-error
2 4 4 4 6 6 6 6 8 8 8 8 8
2.036672e-001 9.938691e-002 6.342811e-002 5.342232e-002 3.272189e-002 2.834954e-002 2.605764e-002 2.463970e-002 1.676902e-002 1.665214e-002 1.665931e-002 1.672104e-002 1.677523e-002
2 8 6 8 8 8 8 8 10 12 14 16 18
2.036672e-001 7.709331e-002 5.400863e-002 3.539708e-002 2.509612e-002 2.038410e-002 1.819450e-002 1.714959e-002 1.206609e-002 8.892393e-003 6.841984e-003 5.441692e-003 4.436991e-003
N/A N/A 1.794590e-001 1.499794e-001 1.181078e-001 9.779428e-002 8.395377e-002 7.383622e-002 6.603440e-002 5.979490e-002 5.467085e-002 5.037954e-002 4.672831e-002
4), (4, 0.5, 3.5), (5, 0, 3), (6, 0, 2), and the associated weights in lexicographic order are 1, 0.5, 0.8, 1.2, 0.6, 1.5, 0.9, 1, 2, 0.3, 0.8, 1.8, 1, 1, 0.5, 0.9, 1.3, 2.4, 0.8, 0.4, 1, 1.6, 0.2, 0.8, 1.2, 1.8, 0.3, 0.8. Here we approximate the rational triangular Bézier surface in Example 2 by a triangular Bézier surface of degree n ¼ 5. According to Theorem 1, we set the weight k ¼ 1:9260. And we use the L1-error and the user-specified error thresholding d ¼ 0:95 in Algorithm 1. Then the iterative algorithm is terminated after the sixth iteration. The original rational triangular surface and the approximation triangular surface are shown in Fig. 1(a–c), and the error plot for the approximation surfaces after the s-th iteration (s = 0,1,. . .,6) is shown in Fig. 1(d). The control points of the approximation triangular Bézier surface after sixth iterations in lexicographic order are (2, 5, 1), (1.8276, 3.9034, 2.2539), (3.4538, 3.6312, 2.6960), (1.2471, 2.8050, 3.6341), (2.5222, 2.0658, 5.8359), (3.0521, 4.1889, 4.2855), (1.0175, 1.9988, 2.4328), (2.2017, 1.1800, 2.4782), (2.3559, 2.5360, 4.8208), (4.2715, 2.4020, 5.2695), (0.2377, 0.4242, 0.0947), (2.5582, 1.7435, 3.2672), (1.5905, 1.0222, 2.2327), (3.2889, 1.4604, 3.4818), (5.3613, 0.7226, 2.9446), (0, 0, 0), (0.3136, 0.3569, 0.0217), (4.2032, 0.0056, 4.9185), (3.3464, 0.5181, 3.5930), (3.8690, 0.4021, 3.5145), (6, 0, 2). Obviously, the approximation surface interpolates the original rational surface at three corner points. In addition, we list the approximation error in the Lp-norm (p = 1, 2, 1) between the original rational surface and the approximation surface after the s-th iteration in Table 1. It shows how the error changes as we increase the iteration number in our method. And it is easy to notice that the Lp-error (p = 1, 2) decreases at first, and then reaches a minimum, finally increases slightly. The reason is that the iterative method only guarantees the points of the approximation surfaces at the uniform parametric values converge to the corresponding points on the original rational surface, but the approximation surfaces do not converge to the rational surface. In case of n ¼ 5, the L2-error reaches a minimum, e2 ¼ 5:061727e 002, at the 24th iteration. In the case of the L2-norm, it is not very difficult to find the best polynomial approximation of the original rational surface which interpolates the corners using the least squares method. It is enough to solve the appropriate system of normal equations. The approximation error of the best polynomial approximation in the L2-norm is e2 ¼ 2:557000e 002. Clearly, the least squares method has a better approximation result than the weighted iterative algorithm. However, the number of rational functions which should be numerical integrated is 504 for the least squares method, but only 24 for the weighted iterative algorithm. The amounts of computation of the two methods both mostly center on the numerical integration of rational functions. Then the iterative algorithm has less calculation if the degrees of the original and approximation surface are large, such as, n,m > 4. And it still works for the case of the Lp (p = 1, 1)-norm. Example 3. Given a rational triangular Bézier surface of degree 5, its control points in lexicographic order are (0, 6, 1), (0.5, 5, 1), (0.5, 5, 1), (1, 4, 2), (0.1, 3.8, 1.9), (1, 4, 2), (2, 2, 2), (0.7, 2.2, 0.5), (0.7, 1.8, 2.5), (2, 2, 2), (3.5, 1, 1), (1.7, 0.7, 2.1), (0, 1.1, 2), (1.8, 1.2, 2.3), (3.5, 1, 1), (4, 0, 0), (2, 0, 2), (0.5, 0.5, 2), (1, 0.3, 1.5), (2, 0.1, 2), (4, 0, 0), and the associated weights in lexicographic order are 1, 1.5, 0.2, 0.9, 1.6, 1.5, 0.3, 1.1, 2, 0.3, 0.7, 1.7, 0.4, 1, 0.5, 1.6, 1.2, 2, 0.4, 1.3, 1. Here we approximate the rational triangular Bézier surface in Example 3 by a triangular Bézier surface of degree n ¼ 7. According to Theorem 1, we set the weight k ¼ 1:9878. And we use the L1-error and the user-specified error thresholding d ¼ 0:9 in Algorithm 1. Then the iterative algorithm is terminated after six iterations. The original rational triangular surface and the approximation triangular surface are shown in Fig. 2(a–b), and the control points of the approximation triangular surface after six iterations in lexicographic order are (0, 6, 1), (0.5193, 4.9835, 0.8832), (0.1046, 5.7768, 0.4973), (0.6701, 4.7567, 0.8538), (0.3055, 3.7301, 2.0277), (1.2048, 3.6849, 2.1471), (0.3154, 5.0342, 1.5518), (0.3032, 3.3630, 1.4127), (0.4245, 2.6486, 2.1289), (1.0435, 3.8753, 2.0139), (1.4059, 3.3541, 1.9785), (0.4402, 2.2063, 1.4765), (0.3180, 2.7162, 1.7041), (0.4931, 2.3851, 2.2251), (0.9764, 3.9326, 1.9728), (3.9288, 0.1503, 0.7662), (1.0727, 0.8124, 1.7323),
9316
Q. Hu / Applied Mathematics and Computation 219 (2013) 9308–9316
(0.4913, 1.5332, 1.4113), (0.5505, 2.0491, 2.0036), (0.6842, 1.7669, 2.7353), (2.8146, 1.5780, 1.2800), (3.8823, 0.2997, 0.1674), (1.7597, 0.5518, 1.8265), (1.4195, 0.5142, 0.5902), (0.6272, 0.6753, 0.2963), (1.1758, 0.8284, 1.3738), (1.9256, 0.5468, 2.3706), (4.1615, 0.0603, 0.1623), (4, 0, 0), (2.7184, 0.0146, 1.0799), (1.1203, 0.2745, 0.3051), (0.9045, 0.4270, 1.4798), (0.1271, 0.3802, 1.3518), (1.8964, 0.0282, 0.8264), (2.5424, 0.1168, 1.8312), (4, 0, 0). Xu and Wang [16] presented a simple method for approximate an m-degree rational triangular Bézier surface by an n(n P m)-degree Bézier surface with control points obtained by degree-elevate the control points of the rational Bézier surface. The original rational surface in Example 3 and the approximation surface of degree n ¼ 7 using Xu&Wang’s method superimposed are shown in Fig 2(c). The error plot for the approximation surfaces obtained by the iterative method and Xu&Wang’s method is shown in Fig. 2(d), where the horizontal axis refers to the degree of the approximation surface. Furthermore, we list the iteration number required and the L1-errors obtained by Algorithm 1 and Xu&Wang’s method, respectively for different approximation degree n(n = 3,4,. . .,15) in Table 2. Obviously, the iterative algorithm is terminated at different iterations with different error thresholdings. Table 2 shows that the larger the error thresholding d is, the larger the iteration number is and the smaller the L1-error is. Compared with Xu&Wang’s method, our iterative algorithm has better approximation results for n = 3,4,. . .,15. Fig 2(d) illuminates the error plot for our algorithm with the error thresholdings d ¼ 0:9 and d ¼ 0:95, and Xu&Wang’s method. However, it is worth mentioning that the approximation triangular surface sequence obtained by Xu&Wang’s method, whose control points are those of degree-elevated rational triangular Bézier surfaces, converges uniformly to the original rational surface. It means that the Lp-error (p = 1, 2, 1) by Xu&Wang’s method will be arbitrarily small if the approximation degree n is large enough. In Example 3, our algorithm has better results than Xu&Wang’s method when the degree n of the approximation surface is not more than 52. Acknowledgements The author thanks the anonymous referees for their valuable suggestions and comments. This work is supported by the National Natural Science Foundation of China (Grants Nos. 61202201, 61272307, and 11201423). The Open Project Program of the State Key Lab of CAD&CD (Grant No. A1305), Zhejiang University. References [1] T.W. Sederberg, M. Kakimoto, Approximating rational curves using polynomial curves, in: G. Farin (Ed.), NURBS for Curve and Surface Design, SIAM, Philadelphia, 1991, pp. 144–158. [2] G.J. Wang, T.W. Sederberg, F.L. Chen, On the convergence of polynomial approximation of rational functions, Journal of Approximation Theory 89 (3) (1997) 267–288. [3] M.S. Floater, High order approximation of rational curves by polynomial curves, Computer Aided Geometric Design 23 (8) (2006) 621–628. [4] Y.D. Huang, H.M. Su, H.W. Lin, A simple method for approximating rational Bézier curve using Bézier curves, Computer Aided Geometric Design 25 (8) (2008) 697–699. [5] C. de Boor, How does Agee’s smoothing method work? In: Proceedings of the 1979 Army Numerical Analysis and Computers Conference, ARO Report 79–3, Army Research, Office, pp. 299–302. [6] D. Qi, Z. Tian, Y. Zhang, J. Feng, The method of numeric polish in curve fitting, Acta Mathematica Sinica 18 (3) (1975) 173–184 (in Chinese). ˇ a, Progressive iteration approximation and the geometric algorithm, Computer-Aided Design 44 (2) (2012) 143–145. [7] J.M. Carnicer, J. Delgado, J.M. Pen [8] J. Delgado, J.M. Penˇ a, Progressive iterative approximation and bases with the fastest convergence rates, Computer Aided Geometric Design 24 (1) (2007) 10–18. [9] H.W. Lin, Local progressive-iterative approximation format for blending curves and patches, Computer Aided Geometric Design 27 (4) (2010) 322–339. [10] H.W. Lin, G.J. Wang, C.S. Dong, Constructing iterative non-uniform B-spline curve and surface to fit data points, Science in China, Series F 47 (3) (2004) 315–331. [11] H.W. Lin, H.J. Bao, G.J. Wang, Totally positive bases and progressive iteration approximation, Computer & Mathematics with Applications 50 (3–4) (2005) 575–586. [12] L.Z. Lu, Sample-based polynomial approximation of rational Bézier curves, Journal of Computational and Applied Mathematics 235 (2011) 1557–1563. [13] L.Z. Lu, Weighted progressive iteration approximation and convergence analysis, Computer Aided Geometric Design 27 (2010) 129–137. [14] G. Farin, B. Piper, A.J. Worsey, The octant of a sphere as a non-degenerate triangular Bézier patch, Computer Aided Geometric Design 4 (4) (1987) 329– 332. [15] G. Farin, Curves and Surfaces for CAGD, 5th ed., Morgan Kaufmann, San Francisco, 2001. [16] H.X. Xu, G.J. Wang, Approximating rational triangular Bézier surfaces by polynomial triangular Bézier surfaces, Journal of Computational and Applied Mathematics 228 (2009) 287–295. [17] L. Zhang, G.J. Wang, An effective algorithm to approximate rational triangular B-B surfaces using polynomial forms, Chinese Journal of Computers 29 (12) (2006) 2151–2162 (in Chinese). [18] J. Chen, G.J. Wang, Progressive-iterative approximation for triangular Bézier surfaces, Computer-Aided Design 43 (8) (2011) 889–895. [19] H.T. Rathod, K.V. Nagaraja, B. Venkatesudu, On the application of two symmetric Gauss Legendre quadrature rules for composite numerical integration over a triangular surface, Applied Mathematics and Computation 190 (2007) 21–39. [20] C. Dunll, Y. Xu, Orthogonal polynomial of several variables, Cambridge University Press, Cambridge, 2001. [21] S. Cooper, S. Waldron, The eigenstructure of the Bernstein operator, Journal of Approximation Theory 105 (1) (2000) 133–165. [22] P.J. Davis, P. Rabinowitz, Methods of Numerical Integration, Dover Publications, 2007.
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
An iterative algorithm for polynomial approximation of rational triangular Bézier surfaces Qianqian Hu Department of Mathematics, Zhejiang Gongshang University, Hangzhou 310018, PR China
a r t i c l e
i n f o
Keywords: Rational triangular Bézier surfaces Polynomial approximation Weighted progressive iteration Triangular Bernstein basis Gauss Legendre quadrature rule
a b s t r a c t This paper presents the weighted progressive iteration approximation (WPIA) property for the triangular Bernstein basis over a triangle domain with uniform parameters, which is extended from the PIA property for triangular Bernstein basis proposed by Chen and Wang in [J. Chen, G.J. Wang, Progressive-iterative approximation for triangular Bézier surfaces, Computer-Aided Design 43 (2011) 889–895]. We also provide how to choose an optimal value of the weight to own the fastest convergence rate for triangular Bernstein basis. Furthermore, a new and efficient iterative method is proposed for polynomial approximation of rational triangular Bézier surfaces. The algorithm is reiterated until a halting condition about approximation error is satisfied. And the approximation error in Lp-norm (p = 1, 2, 1) is calculated by the symmetric Gauss Legendre quadrature rule for composite numerical integration over a triangular surface. Finally, several numerical examples are presented to validate the effectiveness of this method. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction Approximating rational curves and surfaces by polynomial curves and surfaces is a significant problem in Computer Aided Geometric Design (CAGD). There are two advantages for polynomial approximation of rational forms. One is to simplify the integral and differential calculations of rational curves and surfaces; another is to exchange data between different curve and surface design systems using polynomial or rational forms. In recent twenty years, many works on polynomial approximation of rational Bézier curves have been developed. In 1991, Sederberg and Kakimoto [1] presented the hybrid algorithm for polynomial approximation of rational curves for the first time. Then the recursive algorithm and convergence condition for the hybrid approximation were put forward in [2]. However, the hybrid approximation involves complicated calculation and its convergence is not guaranteed. Floater [3] showed that a k-degree rational parametric curve can be interpolated by an m-degree polynomial matching 2m 2k + 4 data. Huang et al. [4] constructed a Bézier curve sequence with control points obtained by degree-elevate the control points of a rational Bézier curve. This method is simple, and can guarantee that the approximation curve sequence converges to the rational curve, but cannot approximate the rational curve by polynomial curves with lower degree. The progressive iteration approximation (PIA) property for polynomial curves was first proposed in [5,6], and attracted a good deal of attention [7–10]. Curves generated by normalized totally positive (NTP) basis, such as Bernstein basis and B-spline basis, satisfy the PIA property [11]. Recently, Lu [12] presented an iterative method for polynomial approximation of rational Bézier curves based on the idea of weighted progressive iteration approximation proposed in [13]. As we all know, it is easy to generalize the abovementioned methods to polynomial approximation of rational tensor product surfaces, because a rational tensor product Bézier surface can be regarded as the space moving trace of a Bézier E-mail addresses: [email protected], [email protected] 0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.03.053
Q. Hu / Applied Mathematics and Computation 219 (2013) 9308–9316
9309
curve with moving control points. Compared with surfaces in rectangular domain, Bernstein–Bézier surfaces over a triangle domain possess a series of beautiful geometric properties, such as be flexible for blending and the non-degenerate representation of segments of a sphere [14]. Also triangular surfaces are adapted to surface interpolation from scattered data and surface design on arbitrary topology mesh [15]. Thus study on the polynomial approximation of rational triangular Bézier surfaces is an urgent need for engineering. Up until now, relatively little has been published about this topic, because the expression form of rational triangular Bézier surfaces are more complex than rational Bézier curves. Enlightened by the degree-elevated vertices approximation [4], Xu and Wang [16] constructed a polynomial triangular Bézier approximation surface, which has the same control points and degree of the rational triangular Bézier surface after degree elevation. Zhang and Wang [17] extended the hybrid approximation to rational triangular Bézier surfaces. Inspired by [18], we present a new and efficient iterative method, the weighted progressive iteration approximation, for polynomial approximation of rational triangular Bézier surfaces. Our method differs from the previous technique [18] in that we multiply all the adjusting control points by some weight to accelerate the iteration convergence rate. And then we apply the WPIA method to approximate rational triangular Bézier surfaces by polynomial triangular surfaces. First of all, we sample some points from the original rational triangular Bézier surface with uniform parameters. Second, we construct a triangular Bézier surface with these sample points as control points and take it as an initial iteration surface. Finally, we use the WPIA method to reiterate the initial surface to approximate the original rational surface. The iterative process is repeated until a satisfactory approximation error is reached. The approximation error in the Lp (p = 1, 2, 1)-norm is calculated with high accuracy by the symmetric Gauss Legendre quadrature rule for composite numerical integration over a triangular surface [19]. The paper is laid out as follows. We give a brief introduction to rational triangular Bézier surfaces in Section 2. In Section 3, we propose the weighted progressive iteration approximation for triangular Bézier surfaces and proof its convergence property. Then we present a numerical calculation of the approximation error between the original rational surface and the approximation surface in different norms, and provide an algorithm for polynomial approximation of rational triangular Bézier surfaces using the weighted progressive iteration approximation in Section 4. Finally, numerical examples are presented in Section 5 to confirm the effectiveness of this algorithm. 2. Preliminary For the following developments, we need first some notations: i ¼ ði; j; kÞ; jij ¼ i þ j þ k. When we say jij ¼ m, we mean i þ j þ k ¼ m, always assuming i; j; k P 0. Definition 1. A rational triangular Bézier surface of degree m over a triangle domain T :¼ fðu; v ; wÞ : u; v ; w P 0; u þ v þ w ¼ 1g is represented by [15]
P
m
xi bi Bmi ðuÞ ; m jij¼m xi Bi ðuÞ
jij¼m
b ðuÞ ¼ P
ð1Þ
where
Bm i ðuÞ ¼
m
i
ui ¼
m! i j k uv w ; i!j!k!
jij ¼ m;
juj ¼ 1;
are the triangular Bernstein basis functions of degree m, bi ¼ bi;j;k are the control points, and xi ¼ xi;j;k are the associated weights. Here are ðmþ1Þðmþ2Þ basis functions. 2 Definition 2. Suppose T; bi ; xi are shown as in Definition 1. Then a rational triangular Bézier surface equivalent to the surface (1) is represented by m
b ðu; v Þ ¼
Pm
m iþj¼0 i;j;mij bi;j;mij Bi;j;mij ðu; Pm m iþj¼0 i;j;mij Bi;j;mij ðu; ; 1
x
x
v
v; 1 u vÞ ; u vÞ
ð2Þ
and is called an m degree rational triangular Bézier surface over a triangle domain D :¼ fðu; v Þ : 0 6 u; v 6 1; u þ v 6 1g Definition 3. Lexicographic order [20]: Given two vectors a; b in dimension d P 1. The vector a is arranged before the vector b, denoted as a b, if the first nonzero entry in the difference a b ¼ ða1 b1 ; . . . ; ad bd Þ is positive. Definition 4. The index set of any positive integer m is defined by Km ¼ fðm; 0; 0Þ; ðm 1; 1; 0Þ; ðm 1; 0; 1Þ; ðm 2; 2; 0Þ; ðm 2; 1; 1Þ; ðm 2; 0; 2Þ; . . . ; ð0; m; 0Þ; ð0; m 1; 1Þ; ð0; m 2; 2Þ; . . . ; ð0; 2; m 2Þ; ð0; 1; m 1Þ; ð0; 0; mÞg. The index set has
ðmþ1Þðmþ2Þ 2
elements sorted in lexicographic order.
9310
Q. Hu / Applied Mathematics and Computation 219 (2013) 9308–9316
Example 1. The control points of a rational triangular Bézier surface of degree 5 are arranged according to the subscripts in 5 lexicographic order as b ¼ fbi gi2K5 ¼ fb500 ; b410 ; b401 ; b320 ; b311 ; b302 ; b230 ; b221 ; b212 ; b203 ; b140 ; b131 ; b122 ; b113 ; b104 ; b050 ; b041 ; b032 ; b023 ; b014 ; b005 g. 3. Weighted progressive iteration approximation and convergence analysis At first, we proof that triangular Bernstein basis over a triangle domain with uniform parameters satisfy the weighted progressive iteration approximation property. Given a sequence of points fpi gi2Kn arranged by the subscripts in lexicographic order in R3 and the corresponding uniform parameter sequence in lexicographic order fui gi2Kn ; ui ¼ ni ; jij ¼ n, then the initial approximation triangular surface of degree n over the triangle domain T is represented by
p0 ðuÞ ¼
X p0i Bni ðuÞ; jij¼n
where p0i ¼ pi : For each point pi , we calculate the adjusting vector
D0i ¼ pi p0 ðui Þ;
jij ¼ n;
and let
p1i ¼ p0i þ kD0i ;
jij ¼ n;
where k is a weight. When k is chosen to some special value, the iteration approximation has the fastest convergence rate. Theorem 1 in the following will show readers how to choose the optimal value of k. Then the first iteration triangular surface over T is expressed as
p1 ðuÞ ¼
X p1i Bni ðuÞ: jij¼n
After the (s + 1)-th iteration, the triangular surface over T can be constructed as follows
psþ1 ðuÞ ¼
X psþ1 Bni ðuÞ; i jij¼n
where
pisþ1 ¼ psi þ kDsi ;
Dsi ¼ pi ps ðui Þ;
jij ¼ n:
ð3Þ
Then we construct a triangular Bézier surface sequence ps ðuÞ; s ¼ 0; 1; . . . According to (3), we obtain
Dsi ¼ Ds1 k i
X s1 n Dj Bj ðui Þ;
jij ¼ n:
jjj¼n
In matrix–vector notation, the abovementioned equation can be rewritten as
h iT s T Di i2Kn ¼ ðI kBÞ Ds1 i
n
i2K
h iT ¼ ðI kBÞs D0i
i2Kn
;
ð4Þ
T where the elements of the column vector Dsi i2Kn are arranged according to the subscript i in lexicographic order, I is the identity matrix of order ðnþ1Þðnþ2Þ , and B ¼ Bni nj is the collocation matrix of the triangular Bernstein basis functions 2 i2Kn ;j2Kn n
i Bi i2Kn at the uniform parametric values n i2Kn , expressed as
1 Bnn;0;0 nn ; 0n ; 0n Bnn1;1;0 nn ; 0n ; 0n Bn0;0;n nn ; 0n ; 0n
C B n n1 1 0 B Bn;0;0 n ; n ; n Bnn1;1;0 n1 ; 1n ; 0n Bn0;0;n n1 ; 1n ; 0n C n n C B ¼B C: .. .. .. .. C B . . . . A @
Bnn1;1;0 0n ; 0n ; nn Bn0;0;n 0n ; 0n ; nn Bnn;0;0 0n ; 0n ; nn 0
B :¼
n Bi n
i i2K n i2Kn
!
ð5Þ
Thus, for the triangular Bézier surface sequence ps ðuÞ; s ¼ 0; 1; . . ., if
limps ðui Þ ¼ pi ;
s!1
ui ¼
i ; n
jij ¼ n;
then we call the initial triangular surface p0 ðuÞ has the weighted progressive iteration approximation property (ab. WPIA).
9311
Q. Hu / Applied Mathematics and Computation 219 (2013) 9308–9316
According to [18,21], we know that all the eigenvalues ni of the collocation matrix B (5) satisfy that 0 < ni 6 1; jij ¼ n. We denote by n ¼ minfni g its minimum, and denote by qðBÞ its spectral radius, which is the maximal absolute value of its eigenjij¼n
values. Then the iterative process (4) converges when qðI kBÞ < 1, and its speed depends on the value of qðI kBÞ. The smaller the value is, the faster the iterative algorithm converges. Obviously, when k ¼ 1, the WPIA property turns out to be the progressive iteration approximation for triangular Bézier surfaces [18]. Next, the following theorem will show us how to choose the optimal value of k so that the iteration approximation has the fastest convergence rate. Theorem 1. Let Bni i2Kn be the triangular Bernstein basis functions, then weighted progressive iteration approximation
based on the uniform parameters ni i2Kn has the fastest convergence rate when
k¼
2nn ; nn þ n!
ð6Þ
and
qðI kBÞ ¼
nn n! nn þ n!
Proof.n According
to
[18,20],o the
eigenvalues
n! Xn ¼ nni ¼ ni ðniÞ! j i ¼ 0; 1; . . . ; n .
ni
of
the
collocation
matrix
B
(5)
belong
to
the
set
It is easy to verify that
nniþ1 n i ¼ < 1; n nni
iP1
and the number sequence nni i¼1;...;n is strictly decreasing, and nni 2 ð0; 1. So the maximum eigenvalue of the matrix B is 1, and the minimum eigenvalue of the matrix B is
n¼
n! nn
ð7Þ
Note that
nni ¼
n! n n1 niþ1 ¼ 6 1; ni ðn iÞ! n n n
i.e., 0 < nni 6 1; i ¼ 0; 1; . . . ; n, so it holds 0 < ni 6 1 i 2 Kn , and the eigenvalues of the matrix I B are equal to n 1 ni 2 ½0; 1Þ i 2 K . This result implies
qðI BÞ ¼ 1 n < 1;
ð8Þ
where n is the minimum eigenvalue of the matrix B defined by (7). First, according to the definition of spectral radius, we have
qðI kBÞ ¼ maxfj1 knj; j1 kjg
ð9Þ
Obviously, qðI kBÞ P j1 kj. Due to the convergence of the iterative process, qðI kBÞ < 1 must hold, i.e., 0 < k < 2. Second, for any k 2 ð0; 1Þ, according to (8) and (9), it holds
qðI kBÞ ¼ j1 knj ¼ 1 kn > 1 n ¼ qðI BÞ: Therefore, when k ¼ 1, the convergence rate is always faster than that of k 2 ð0; 1Þ. Finally, for any k 2 ð1; 2Þ, in the case of kn > 1, according to (8) and (9), it holds
qðI kBÞ ¼ k 1 > n1 1 > 1 n ¼ qðI BÞ: It means that when k ¼ 1, the convergence rate is always faster than that in such case. On the other hand, in the case of kn 6 1, we have
(
qðI kBÞ ¼ maxf1 kn; k 1g ¼
1 n; 1 < k < 2ð1 þ nÞ1 k 1; 2ð1 þ nÞ1 6 k < 2
:
Clearly, qðI kBÞ reaches the minimum when k ¼ 2ð1 þ nÞ1 . In other words, the weighted progressive iteration approximation based on the uniform parameters has the fastest convergence rate when
k¼
2 ; 1þn
and in such case, the spectral radius of the matrix ðI kBÞ is
9312
Q. Hu / Applied Mathematics and Computation 219 (2013) 9308–9316
qðI kBÞ ¼
1n : 1þn
Substituting (7) into the two abovementioned equations, the theorem is proved. h
4. Error measurement and algorithm description m
We consider the problem of approximating the rational triangular Bézier surface b ðuÞ (1) by an n degree triangular Bézier surface defined by
pðuÞ ¼
X n pi Bi ðuÞ:
ð10Þ
jij¼n
using the weighted progressive iteration approximation property. Since we have only proved that triangular Bézier surfaces with uniform parameters have the WPIA property, we sample m ðnþ1Þðnþ2Þ points p0i i 2 Kn from the original surface b ðuÞ with the uniform parameters ui ¼ ni ¼ ni ; nj ; nk ; ði 2 Kn Þ, i.e., 2 m
p0i ¼ b ðui Þ;
ui ¼
i ; n
jij ¼ n:
ð11Þ
for the iterative process. It should be emphasized that we are interested in a good approximation to the rational surface m m b ðuÞ rather than the sampled points b ðui Þði 2 Kn Þ. Next, we provide the error measurement for the approximation surface. The approximation error in the Lp-norm is defined by
Fig. 1. (a) The rational triangular Bézier surface of degree 6; (b) The approximation triangular surface after 6 iterations; (c) The original surface, the approximation surfaces after 0 and 6 iterations superimposed; (d) The error plot for the approximation surfaces after different iterations.
Q. Hu / Applied Mathematics and Computation 219 (2013) 9308–9316
ep ¼
Z Z
m jb ðuÞ pðuÞ jp dA
1p ;
1 6 p < 1;
9313
ð12Þ
T m
where b ðuÞ is the original rational surface (1), and A is the area of the triangle domain T, and when p = 1, the error is defined by
e1 ¼ maxkbm ðuÞ pðuÞk: u2T
ð13Þ
For the double integral over the triangle domain T in (12), it is very hard to integrate it explicitly. So we consider using m m numerical integration methods for good results. Since the surface b ðuÞ (1) is equivalent to the surface b ðu; v Þ (2), we can change the triangle domain T :¼ fðu; v ; wÞ : u; v ; w P 0; u þ v þ w ¼ 1g as D :¼ fðu; v Þ : 0 6 u; v 6 1; u þ v 6 1g. Then we denote by pðu; v Þ the approximation surface pðuÞ defined over D. Clearly, the integral region D is the standard triangular surface in the Cartesian two dimensional (u, v) space, and (12) can be revised as
ep ¼
Z
1
dv
Z
0
1v
1p f ðu; v Þdu ;
ð14Þ
0 m
where f ðu; v Þ ¼ kb ðu; v Þ pðu; v Þkp . To yield an accurate and reliable evaluation of the integral (14), a symmetric Gauss Legendre quadrature rule is provided for composite numerical integration over the standard triangular surface D [19]. We discretize the triangle domain D in (u, v) space into N N = N2 right isosceles triangle Di each of area 1/(2N2), and divide each triangle Di into three new triangles of equal area 1/(6N2) by joining the centroidal point to the three vertices of Di. Readers can use the q-order of Gauss Legendre quadrature rule for numerical integration, and refer to [19,22] for more details. As Rathod et al. said in [19], the numerical integration algorithm converges to the exact value of the integral, for sufficiently large value of N and the convergence is
Fig. 2. (a) The rational triangular Bézier surface of degree 5; (b) The original rational surface and the approximation surface using our method superimposed; (c) The original surface and the approximation surface using Xu&Wang’s method superimposed; (d) The error plot for approximation surfaces obtained by two different methods.
9314
Q. Hu / Applied Mathematics and Computation 219 (2013) 9308–9316
Table 1 The Lp-error (p = 1, 2, 1) between the k-th iteration triangular surface and the original rational surface. k
L1
L2
L1
k
L1
L2
L1
0
2.652368e-001 1.627374e-001 1.171966e-001 8.241156e-002 7.545272e-002 6.758371e-002 6.564027e-002 6.215079e-002 6.017676e-002 5.781102e-002 5.593812e-002 4.905877e-002 4.595881e-002 4.535136e-002 4.559194e-002
3.090208e-001 1.836061e-001 1.302147e-001 8.674335e-002 8.375550e-002 7.258315e-002 7.180955e-002 6.739415e-002 6.546666e-002 6.271153e-002 6.071003e-002 5.352006e-002 5.097849e-002 5.064423e-002 5.109564e-002
7.765044e-001 4.574360e-001 3.384068e-001 2.755662e-001 2.686409e-001 2.774455e-001 2.745880e-001 2.763090e-001 2.745905e-001 2.745392e-001 2.736165e-001 2.735843e-001 2.762940e-001 2.795807e-001 2.824704e-001
35 40 45 50 55 60 65 70 75 80 85 90 95 100 105
4.611488e-002 4.664660e-002 4.706643e-002 4.738014e-002 4.760596e-002 4.776406e-002 4.787312e-002 4.794862e-002 4.800051e-002 4.803631e-002 4.806085e-002 4.807760e-002 4.808906e-002 4.809688e-002 4.810221e-002
5.171207e-002 5.227033e-002 5.271267e-002 5.304198e-002 5.327903e-002 5.344630e-002 5.356289e-002 5.364353e-002 5.369901e-002 5.373706e-002 5.376310e-002 5.378088e-002 5.379302e-002 5.380130e-002 5.380694e-002
2.847315e-001 2.864040e-001 2.876034e-001 2.884479e-001 2.890356e-001 2.894418e-001 2.897210e-001 2.899125e-001 2.900434e-001 2.901329e-001 2.901939e-001 2.902355e-001 2.902639e-001 2.902832e-001 2.902964e-001
1 2 3 4 5 6 7 8 9 10 15 20 25 30
much faster for higher order (q-order) rules. To achieve high accuracy of numerical integration, it is desirable that N and q are as large as possible. However, since the computation time increases as N and q increase, we need to seek a balance between the speed and the accuracy. Experimentally, we found that q = 2, N = 10 is suitable for good results. When p = 1, (13) can be estimated by
i N
m e1 ¼ max kb ðui Þ pðui Þk; ui ¼ : N i2K
ð15Þ
After our testing, we found that the estimation works well when 50 6 N 6 100. After the error measurement, we present the iterative algorithm. First, we sample ðnþ1Þðnþ2Þ points p0i ði 2 Kn Þ (11) from the 2 m rational triangular surface b ðuÞ for the iterative process. The initial approximation surface is expressed by P p0 ðuÞ ¼ jij¼n p0i Bni ðuÞ; and clearly has the WPIA property. Thus we can construct a triangular Bézier surface sequence m s p ðuÞ; s ¼ 0; 1; :, and the limit surface interpolates the ðnþ1Þðnþ2Þ points b ðui Þði 2 Kn Þ 2 m Clearly, the approximation surface pðuÞ (10) interpolates three corner points of the original rational surface b ðuÞ (1), since we have un;0;0 ¼ ð1; 0; 0Þ; u0;n;0 ¼ ð0; 1; 0Þ; and u0;0;n ¼ ð0; 0; 1Þ. n Second, according to Theorem 1, we choose the optimal value of weight as k ¼ n2n n þn! for the iteration approximation with the fastest convergence rate. It should be pointed out that the iterative method converges slowly, even for the best choice of the parameter k, compared with Xu&Wang’s method, i.e., the degree-elevated method. However, it has a better approximation result than Xu&Wang’s method. Finally, in order to find a good approximation surface and improve the iteration speed, the iterative algorithm can be reiterated until a halting condition about the approximation error is satisfied. The termination condition iskþ1 the relative differje ek j ence between the two consecutive iteration errors ek and ekþ1 is less than a given threshold value f, i.e., ek 6 f; f 2 ð0; 1Þ; or the latter iteration error ekþ1 is larger than the former error ek as ekþ1 > ek . So the termination condition is ekþ1 P dek ; d 2 ð0; 1Þ: We summarize the iterative algorithm as follows. Algorithm 1 (Polynomial approximation of a rational triangular Bézier surface). m Input: An m degree rational triangular Bézier surface b ðuÞ with the control points fbi gi2Km and weights fxi gi2Km , a userspecified error thresholding d, and p in Lp-error. Output: An n degree triangular Bézier surface pðuÞ with the control points fpi gi2Kn , the iteration number s, and the iteration error e. Step Step Step Step Step
1. 2. 3. 4. 5.
Let k ¼ 0. When p ¼ 1; 2, let q ¼ 2; N ¼ 10; when p ¼ 1, let N ¼ 100: Compute the weight k by (6). Compute p0i i2Kn of the initial iteration by (11), and the initial iteration error e0p by (14) or (15). Computer psþ1 of the ðs þ 1Þ-th iteration by (3) and the iteration error esþ1 by (14) or (15). p i i2Kn If esþ1 6 des ; then set s ¼ s þ 1 and go to Step 4. Else, set pi ¼ psi ; i 2 Kn and e ¼ esp . Then output fpi gi2Kn ; e; s; and stop.
5. Numerical examples
Example 2. Given a rational triangular Bézier surface of degree 6, its control points in lexicographic order are (2, 5, 1), (2, 4, 2), (3, 4.2, 2), (1.5, 3.5, 3), (2.5, 3, 5), (3.5, 3.5, 4), (1, 2, 3), (1.5, 2.8, 4), (2.5, 2.5, 5), (4, 3, 6), (0.8, 1.5, 1.5), (2, 1.2, 2), (3, 2, 3), (4, 2.2, 4), (4.5, 2, 5), (0.4, 0.8, 0.6), (1.2, 1, 2), (2, 0.9, 3), (3, 0.6, 4), (4.5, 1, 3), (5.2, 1, 3), (0, 0, 0), (0.4, 0.4, 1), (1.5, 0.5, 2), (3, 0.2,
9315
Q. Hu / Applied Mathematics and Computation 219 (2013) 9308–9316 Table 2 The L1-errors of two different approximations by triangular Bézier surfaces of degree n. n
Our algorithm (d ¼ 0:9)
3 4 5 6 7 8 9 10 11 12 13 14 15
Our algorithm (d ¼ 0:95)
Xu&Wang’s method
Iteration number
L1-error
Iteration number
L1-error
L1-error
2 4 4 4 6 6 6 6 8 8 8 8 8
2.036672e-001 9.938691e-002 6.342811e-002 5.342232e-002 3.272189e-002 2.834954e-002 2.605764e-002 2.463970e-002 1.676902e-002 1.665214e-002 1.665931e-002 1.672104e-002 1.677523e-002
2 8 6 8 8 8 8 8 10 12 14 16 18
2.036672e-001 7.709331e-002 5.400863e-002 3.539708e-002 2.509612e-002 2.038410e-002 1.819450e-002 1.714959e-002 1.206609e-002 8.892393e-003 6.841984e-003 5.441692e-003 4.436991e-003
N/A N/A 1.794590e-001 1.499794e-001 1.181078e-001 9.779428e-002 8.395377e-002 7.383622e-002 6.603440e-002 5.979490e-002 5.467085e-002 5.037954e-002 4.672831e-002
4), (4, 0.5, 3.5), (5, 0, 3), (6, 0, 2), and the associated weights in lexicographic order are 1, 0.5, 0.8, 1.2, 0.6, 1.5, 0.9, 1, 2, 0.3, 0.8, 1.8, 1, 1, 0.5, 0.9, 1.3, 2.4, 0.8, 0.4, 1, 1.6, 0.2, 0.8, 1.2, 1.8, 0.3, 0.8. Here we approximate the rational triangular Bézier surface in Example 2 by a triangular Bézier surface of degree n ¼ 5. According to Theorem 1, we set the weight k ¼ 1:9260. And we use the L1-error and the user-specified error thresholding d ¼ 0:95 in Algorithm 1. Then the iterative algorithm is terminated after the sixth iteration. The original rational triangular surface and the approximation triangular surface are shown in Fig. 1(a–c), and the error plot for the approximation surfaces after the s-th iteration (s = 0,1,. . .,6) is shown in Fig. 1(d). The control points of the approximation triangular Bézier surface after sixth iterations in lexicographic order are (2, 5, 1), (1.8276, 3.9034, 2.2539), (3.4538, 3.6312, 2.6960), (1.2471, 2.8050, 3.6341), (2.5222, 2.0658, 5.8359), (3.0521, 4.1889, 4.2855), (1.0175, 1.9988, 2.4328), (2.2017, 1.1800, 2.4782), (2.3559, 2.5360, 4.8208), (4.2715, 2.4020, 5.2695), (0.2377, 0.4242, 0.0947), (2.5582, 1.7435, 3.2672), (1.5905, 1.0222, 2.2327), (3.2889, 1.4604, 3.4818), (5.3613, 0.7226, 2.9446), (0, 0, 0), (0.3136, 0.3569, 0.0217), (4.2032, 0.0056, 4.9185), (3.3464, 0.5181, 3.5930), (3.8690, 0.4021, 3.5145), (6, 0, 2). Obviously, the approximation surface interpolates the original rational surface at three corner points. In addition, we list the approximation error in the Lp-norm (p = 1, 2, 1) between the original rational surface and the approximation surface after the s-th iteration in Table 1. It shows how the error changes as we increase the iteration number in our method. And it is easy to notice that the Lp-error (p = 1, 2) decreases at first, and then reaches a minimum, finally increases slightly. The reason is that the iterative method only guarantees the points of the approximation surfaces at the uniform parametric values converge to the corresponding points on the original rational surface, but the approximation surfaces do not converge to the rational surface. In case of n ¼ 5, the L2-error reaches a minimum, e2 ¼ 5:061727e 002, at the 24th iteration. In the case of the L2-norm, it is not very difficult to find the best polynomial approximation of the original rational surface which interpolates the corners using the least squares method. It is enough to solve the appropriate system of normal equations. The approximation error of the best polynomial approximation in the L2-norm is e2 ¼ 2:557000e 002. Clearly, the least squares method has a better approximation result than the weighted iterative algorithm. However, the number of rational functions which should be numerical integrated is 504 for the least squares method, but only 24 for the weighted iterative algorithm. The amounts of computation of the two methods both mostly center on the numerical integration of rational functions. Then the iterative algorithm has less calculation if the degrees of the original and approximation surface are large, such as, n,m > 4. And it still works for the case of the Lp (p = 1, 1)-norm. Example 3. Given a rational triangular Bézier surface of degree 5, its control points in lexicographic order are (0, 6, 1), (0.5, 5, 1), (0.5, 5, 1), (1, 4, 2), (0.1, 3.8, 1.9), (1, 4, 2), (2, 2, 2), (0.7, 2.2, 0.5), (0.7, 1.8, 2.5), (2, 2, 2), (3.5, 1, 1), (1.7, 0.7, 2.1), (0, 1.1, 2), (1.8, 1.2, 2.3), (3.5, 1, 1), (4, 0, 0), (2, 0, 2), (0.5, 0.5, 2), (1, 0.3, 1.5), (2, 0.1, 2), (4, 0, 0), and the associated weights in lexicographic order are 1, 1.5, 0.2, 0.9, 1.6, 1.5, 0.3, 1.1, 2, 0.3, 0.7, 1.7, 0.4, 1, 0.5, 1.6, 1.2, 2, 0.4, 1.3, 1. Here we approximate the rational triangular Bézier surface in Example 3 by a triangular Bézier surface of degree n ¼ 7. According to Theorem 1, we set the weight k ¼ 1:9878. And we use the L1-error and the user-specified error thresholding d ¼ 0:9 in Algorithm 1. Then the iterative algorithm is terminated after six iterations. The original rational triangular surface and the approximation triangular surface are shown in Fig. 2(a–b), and the control points of the approximation triangular surface after six iterations in lexicographic order are (0, 6, 1), (0.5193, 4.9835, 0.8832), (0.1046, 5.7768, 0.4973), (0.6701, 4.7567, 0.8538), (0.3055, 3.7301, 2.0277), (1.2048, 3.6849, 2.1471), (0.3154, 5.0342, 1.5518), (0.3032, 3.3630, 1.4127), (0.4245, 2.6486, 2.1289), (1.0435, 3.8753, 2.0139), (1.4059, 3.3541, 1.9785), (0.4402, 2.2063, 1.4765), (0.3180, 2.7162, 1.7041), (0.4931, 2.3851, 2.2251), (0.9764, 3.9326, 1.9728), (3.9288, 0.1503, 0.7662), (1.0727, 0.8124, 1.7323),
9316
Q. Hu / Applied Mathematics and Computation 219 (2013) 9308–9316
(0.4913, 1.5332, 1.4113), (0.5505, 2.0491, 2.0036), (0.6842, 1.7669, 2.7353), (2.8146, 1.5780, 1.2800), (3.8823, 0.2997, 0.1674), (1.7597, 0.5518, 1.8265), (1.4195, 0.5142, 0.5902), (0.6272, 0.6753, 0.2963), (1.1758, 0.8284, 1.3738), (1.9256, 0.5468, 2.3706), (4.1615, 0.0603, 0.1623), (4, 0, 0), (2.7184, 0.0146, 1.0799), (1.1203, 0.2745, 0.3051), (0.9045, 0.4270, 1.4798), (0.1271, 0.3802, 1.3518), (1.8964, 0.0282, 0.8264), (2.5424, 0.1168, 1.8312), (4, 0, 0). Xu and Wang [16] presented a simple method for approximate an m-degree rational triangular Bézier surface by an n(n P m)-degree Bézier surface with control points obtained by degree-elevate the control points of the rational Bézier surface. The original rational surface in Example 3 and the approximation surface of degree n ¼ 7 using Xu&Wang’s method superimposed are shown in Fig 2(c). The error plot for the approximation surfaces obtained by the iterative method and Xu&Wang’s method is shown in Fig. 2(d), where the horizontal axis refers to the degree of the approximation surface. Furthermore, we list the iteration number required and the L1-errors obtained by Algorithm 1 and Xu&Wang’s method, respectively for different approximation degree n(n = 3,4,. . .,15) in Table 2. Obviously, the iterative algorithm is terminated at different iterations with different error thresholdings. Table 2 shows that the larger the error thresholding d is, the larger the iteration number is and the smaller the L1-error is. Compared with Xu&Wang’s method, our iterative algorithm has better approximation results for n = 3,4,. . .,15. Fig 2(d) illuminates the error plot for our algorithm with the error thresholdings d ¼ 0:9 and d ¼ 0:95, and Xu&Wang’s method. However, it is worth mentioning that the approximation triangular surface sequence obtained by Xu&Wang’s method, whose control points are those of degree-elevated rational triangular Bézier surfaces, converges uniformly to the original rational surface. It means that the Lp-error (p = 1, 2, 1) by Xu&Wang’s method will be arbitrarily small if the approximation degree n is large enough. In Example 3, our algorithm has better results than Xu&Wang’s method when the degree n of the approximation surface is not more than 52. Acknowledgements The author thanks the anonymous referees for their valuable suggestions and comments. This work is supported by the National Natural Science Foundation of China (Grants Nos. 61202201, 61272307, and 11201423). The Open Project Program of the State Key Lab of CAD&CD (Grant No. A1305), Zhejiang University. References [1] T.W. Sederberg, M. Kakimoto, Approximating rational curves using polynomial curves, in: G. Farin (Ed.), NURBS for Curve and Surface Design, SIAM, Philadelphia, 1991, pp. 144–158. [2] G.J. Wang, T.W. Sederberg, F.L. Chen, On the convergence of polynomial approximation of rational functions, Journal of Approximation Theory 89 (3) (1997) 267–288. [3] M.S. Floater, High order approximation of rational curves by polynomial curves, Computer Aided Geometric Design 23 (8) (2006) 621–628. [4] Y.D. Huang, H.M. Su, H.W. Lin, A simple method for approximating rational Bézier curve using Bézier curves, Computer Aided Geometric Design 25 (8) (2008) 697–699. [5] C. de Boor, How does Agee’s smoothing method work? In: Proceedings of the 1979 Army Numerical Analysis and Computers Conference, ARO Report 79–3, Army Research, Office, pp. 299–302. [6] D. Qi, Z. Tian, Y. Zhang, J. Feng, The method of numeric polish in curve fitting, Acta Mathematica Sinica 18 (3) (1975) 173–184 (in Chinese). ˇ a, Progressive iteration approximation and the geometric algorithm, Computer-Aided Design 44 (2) (2012) 143–145. [7] J.M. Carnicer, J. Delgado, J.M. Pen [8] J. Delgado, J.M. Penˇ a, Progressive iterative approximation and bases with the fastest convergence rates, Computer Aided Geometric Design 24 (1) (2007) 10–18. [9] H.W. Lin, Local progressive-iterative approximation format for blending curves and patches, Computer Aided Geometric Design 27 (4) (2010) 322–339. [10] H.W. Lin, G.J. Wang, C.S. Dong, Constructing iterative non-uniform B-spline curve and surface to fit data points, Science in China, Series F 47 (3) (2004) 315–331. [11] H.W. Lin, H.J. Bao, G.J. Wang, Totally positive bases and progressive iteration approximation, Computer & Mathematics with Applications 50 (3–4) (2005) 575–586. [12] L.Z. Lu, Sample-based polynomial approximation of rational Bézier curves, Journal of Computational and Applied Mathematics 235 (2011) 1557–1563. [13] L.Z. Lu, Weighted progressive iteration approximation and convergence analysis, Computer Aided Geometric Design 27 (2010) 129–137. [14] G. Farin, B. Piper, A.J. Worsey, The octant of a sphere as a non-degenerate triangular Bézier patch, Computer Aided Geometric Design 4 (4) (1987) 329– 332. [15] G. Farin, Curves and Surfaces for CAGD, 5th ed., Morgan Kaufmann, San Francisco, 2001. [16] H.X. Xu, G.J. Wang, Approximating rational triangular Bézier surfaces by polynomial triangular Bézier surfaces, Journal of Computational and Applied Mathematics 228 (2009) 287–295. [17] L. Zhang, G.J. Wang, An effective algorithm to approximate rational triangular B-B surfaces using polynomial forms, Chinese Journal of Computers 29 (12) (2006) 2151–2162 (in Chinese). [18] J. Chen, G.J. Wang, Progressive-iterative approximation for triangular Bézier surfaces, Computer-Aided Design 43 (8) (2011) 889–895. [19] H.T. Rathod, K.V. Nagaraja, B. Venkatesudu, On the application of two symmetric Gauss Legendre quadrature rules for composite numerical integration over a triangular surface, Applied Mathematics and Computation 190 (2007) 21–39. [20] C. Dunll, Y. Xu, Orthogonal polynomial of several variables, Cambridge University Press, Cambridge, 2001. [21] S. Cooper, S. Waldron, The eigenstructure of the Bernstein operator, Journal of Approximation Theory 105 (1) (2000) 133–165. [22] P.J. Davis, P. Rabinowitz, Methods of Numerical Integration, Dover Publications, 2007.