Optimal multi-degree reduction of Bézier curves with G ... - ScienceDirect
Computer Aided Geometric Design 23 (2006) 673–683 www.elsevier.com/locate/cagd
Optimal multi-degree reduction of Bézier curves with G2 -continuity Lizheng Lu ∗ , Guozhao Wang Institute of Computer Graphics and Image Processing, Department of Mathematics, Zhejiang University, Hangzhou 310027, China Received 22 November 2005; received in revised form 7 September 2006; accepted 16 September 2006 Available online 30 October 2006
Abstract In this paper we present a novel approach to consider the multi-degree reduction of Bézier curves with G2 -continuity in L2 norm. The optimal approximation is obtained by minimizing the objective function based on the L2 -error between the two curves. In contrast to traditional methods, which typically consider the components of the curve separately, we use geometric information on the curve to generate the degree reduction. So positions, tangents and curvatures are preserved at the two endpoints. For avoiding the singularities at the endpoints, regularization terms are added to the objective function. Finally, numerical examples demonstrate the effectiveness of our algorithms. © 2006 Elsevier B.V. All rights reserved. Keywords: Degree reduction; Bézier curves; Optimal approximation; G2 -continuity; L2 -norm
1. Introduction Degree reduction of polynomial curves and surfaces is a common process in CAGD. It consists of approximating a polynomial by another one of lower degree. This process is of great importance in geometric modelling, such as data exchange, data compression and data comparison. For example, degree reduction is needed when data are transferred from one modelling system to another and these systems have different limitations on the maximum degree of polynomials. Furthermore, it can also be used to generate a piecewise continuous lower degree approximation to a given curve or surface so as to simplify some geometric or graphical algorithms like intersection calculation or rendering. 1.1. Related work There have been many methods developed for degree reduction. Forrest (1972) and Farin (1983) considered the inverse of degree elevation and obtained two sets of control points. This approach is easy to compute, but not optimal in the usual Lp -norm for any p = 1, 2, . . . , +∞. Since degree reduction is an approximation problem in nature, methods in the classical approximation theory (Szegö, 1975) can be employed. In particular, the optimal approximations with respect to the L∞ or L2 metric are of interest. Watkins and Worsey (1988) used Chebyshev economization to produce * Corresponding author.
E-mail addresses: [email protected] (L. Lu), [email protected] (G. Wang). 0167-8396/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cagd.2006.09.002
674
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
the best L∞ -approximation of degree n − 1 to a given degree n polynomial. This best approximation, however, does not interpolate the given curve at the endpoints. The endpoint constraints are frequently required in many applications, especially when degree reduction is combined with subdivision to generate piecewise continuous approximations. Recently, C k -constrained optimal degree reduction has been presented. Eck (1995) used constrained Legendre polynomials to minimize the L2 -norm between the two curves. Ahn (2003) presented one degree reduction in L∞ norm. And multi-degree reduction at one time avoiding stepwise computing was investigated in (Ahn et al., 2004; Chen and Wang, 2002; Zheng and Wang, 2003), which showed that the optimal approximation in L2 -norm can be obtained by different methods and the results are equivalent. For every Lp -norm there exists an optimal approximation, but which one is the best is an open problem. In the above mentioned papers the approximation error is measured between the parameterized curves. Pottmann et al. (2002) applied active curves to degree reduction, thus minimizing the geometric distance between the curves. The key idea is to iteratively change the control points of the active curve so that the active curve deforms towards the given curve. And their method shows the better approximation obtained by the geometric distance. 1.2. Motivation and outline of the present paper Traditional methods for degree reduction are applied to parametric curves simply by matching the position and derivatives at the same parameter values, in general f (r) (t) = g (r) (t),
r = 0, 1, . . . , α.
Since parametric representations of curves are not unique, they will produce different approximations for a given curve. For example, Fig. 1 shows a quartic polynomial (solid): F(x) = (x, 2x − 6x 2 + 5x 3 − x 4 ),
x ∈ [0, 1],
and its cubic best L2 -approximation with C 1 -continuity (dotted) by the method in (Ahn et al., 2004). However, if one reparameterizes F(x) by taking x = φ(t) = (t + t 2 )/2, the curve F(φ(t)) remains unchanged but with a different approximating curve (dashed). The corresponding approximations with G1 -continuity (with “+”) by our method are also shown in Fig. 1. Our main motivation is to consider degree reduction with geometric continuity, which is less dependent on parameterizations. And such method gives better approximation for degree reduction. Geometric continuity was introduced by de Boor et al. (1987) to curve interpolation, which is called Geometric Hermite Interpolation (GHI). They showed that if the curvature at one endpoint is not vanished, a planar curve can be interpolated by cubic spline with G2 -continuity and the approximation order is 6. For more details about GHI, see the recent survey by Degen (2005).
Fig. 1. A comparison of C 1 - and G1 -constrained degree reduction.
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
675
In this paper, we consider multi-degree reduction of Bézier curves with geometric continuity by minimizing the objective function based on the L2 -error between the two curves. For G1 -continuity, it provides two additional parameters compared to C 1 approximation. And for G2 -continuity, the additional parameters are four. By using these parameters, we can optimize the approximation and obtain the optimal approximating curve in L2 -norm. The outline of this paper is as follows. We briefly characterize the problem of Gα -constrained multi-degree reduction in Section 2. In Section 3, we consider G1 -constrained degree reduction and propose the algorithm. We then outline G2 -constrained degree reduction in Section 4. In Section 5 we provide some examples and compare our method with traditional methods. Finally, Section 6 concludes this paper. 2. Preliminaries 2.1. Definitions and notations the space of parametric polynomials in Rd with degree at most n and · denotes the In this paper, Πnd denotes √ Euclidean norm v = v, v. A Bézier curve of degree n is defined by the control points pi ∈ Rd in the form P(t) =
n
Bin (t)pi ,
0 t 1,
i=0
where Bin (t) are the Bernstein polynomials given by Bin (t) = operator as follows, pi = pi+1 − pi ,
r pi = (r−1 pi ),
n n−i t i . And we define the forward difference i (1 − t)
r 2.
The product of two Bernstein polynomials is given by mn i
j
m+n (t) Bim (t)Bjn (t) = m+n Bi+j
(1)
i+j
and the integral is 1 Bin (t) dt = 0
1 . n+1
(2)
We introduce the following lemma, where mij can be easily obtained from (1) and (2). Lemma 1. Let Mm,n = (mij ) be an (m + 1) × (n + 1) matrix with the elements given by 1 mij =
Bim (t)Bjn (t) dt, 0
then Mm,m is a real symmetric positive definite matrix. Proof. For any nonzero column vector ξ with m + 1 elements, we have 2 1 m m Bi (t)ξi dt = ξ T Mm,m ξ > 0. 2 0
i=0
2.2. Problem statement We start with the problem of Gα -constrained degree reduction. Given a degree n Bézier curve P(t) =
n i=0
Bin (t)pi ,
(3)
676
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
find a degree m (m < n) Bézier curve Q(t) =
m
Bim (t)qi
i=0
such that the following two conditions are satisfied: (i) P(t) and Q(t) are Gα -continuous (2α m − 1) at t = 0, 1, i.e., Q(r) (v) = P(r) φ(v) , v = 0, 1, r = 0, 1, . . . , α,
(4)
where φ : [0, 1] → [0, 1] is an strictly increasing function with φ(0) = 0 and φ(1) = 1; 1 (ii) Q(t) minimizes the squared L2 -error ε = 0 P(t) − Q(t)2 dt for all possible polynomials in Πmd . Note that traditional methods consider degree reduction of Bézier curves with the constraints of C α -continuity, which fix the first α + 1 and the last α + 1 control points of the approximating curve. In contrast, Gα -continuity will bring additional parameters to optimize the approximation. By using these parameters, we can obtain the approximating curve with a smaller approximation error. However, it turns out to be a nonlinear problem and brings more difficulties to degree reduction when α 1. 3. G1 -constrained degree reduction 3.1. G1 condition Clearly, for G0 -continuity, the endpoints of Q(t) should coincide with the endpoints of P(t). And for G1 -continuity, the coincidence of the oriented tangents is additionally needed. From (4), we get Q (0) = φ (0)P (0) Denoting
φ (0) = δ02
and
and Q (1) = φ (1)P (1). φ (1) = δ12 ,
(5) G1
we can easily relate condition (5) with the control points. More precisely, n 2 q0 = p0 , q1 = p0 + δ0 (p1 − p0 ), m (6) n qm−1 = pn − δ12 (pn − pn−1 ). qm = pn , m Noting that geometric continuity does not depend on the chosen parameterization, the control point q1 can thus move along the direction − p−0−p→1 without violating G1 condition (see Fig. 2), and the same holds for qm−1 . Thus it provides two additional parameters to optimize the shape of the degree reduced curve. When replacing G1 -continuity with C 1 -continuity, q1 and qm−1 are uniquely determined by δv = 1. Therefore, 1 C approximation is a special case of G1 approximation. We can also imagine that G1 approximation will lead to a smaller error between the original curve and the approximating one. We will propose the algorithm in Section 3.2. 3.2. Algorithm for G1 -constrained degree reduction For a given degree n Bézier curve P(t), the problem of G1 -constrained degree reduction can be solved through two stages. In the first stage, we construct a degree m Bézier curve Q(t) interpolating P(t) according to (6). More precisely, it can be expressed as Q(t) = B0m (t)q0 + B1m (t)q1 +
m−2
m m Bim (t)qi + Bm−1 (t)qm−1 + Bm (t)qm ,
(7)
i=2
where q1 and qm−1 contain the unknown parameters δ0 and δ1 , respectively. We assume that the free parameters δ0 and δ1 are temporarily fixed, then solve the interior control points qi , i = 2, . . . , m − 2, by minimizing 1 ε= 0
P(t) − Q(t) 2 dt.
(8)
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
677
Fig. 2. G1 condition for degree reduction.
Denoting Bn = (B0n (t), . . . , Bnn (t)) and Pn = (p0 , . . . , pn )T , we have 1 ε= 0
1 =
P(t) − Q(t) 2 dt =
1 Bn Pn − Bm Qm 2 dt 0
Bn Pn − Bc Qc − Bfm Qfm 2 dt, m m
(9)
0 f
where Qcm = (q0 , q1 , qm−1 , qm )T and Qm denotes other control points of Qm . f For a minimum of (9), it is necessary that the derivatives of ε with respect to the elements of Qm are zero, which gives 0=
1 ∂ε f f f = Mm,n Pn − Mcm,m Qcm − Mm,m Qm , 2 ∂Qfm
(10)
where f
Mm,n := Mm,n (2, . . . , m − 2; 0, . . . , n), Mcm,m := Mm,m (2, . . . , m − 2; 0, 1, m − 1, m), f
Mm,m := Mm,m (2, . . . , m − 2; 2, . . . , m − 2). The notation A(. . . ; . . .) denotes the submatrix of the matrix A obtained by extracting the specific rows and columns. f Since Mm,m is a real symmetric positive definite matrix (see Lemma 1), Mm,m is invertible. Therefore, ε is minimized by choosing f −1 f f f Qm = Qm (δ0 , δ1 ) = Mm,m Mm,n Pn − Mcm,m Qcm . (11) Note that q1 and qm−1 are quadratic functions of the parameters δ0 and δ1 , respectively. Therefore, from (11), the interior control points qi , i = 2, . . . , m − 2, are quadratic functions of the parameters δ0 and δ1 . As shown by Fig. 3, we can obtain different G1 -constrained approximating curves to the curve F(x) for different parameter values δ0 and δ1 . The second stage is then to determine the two free parameters δv such that the approximation error (8) is minimized. We will propose how to optimize the two parameters so as to obtain the optimal L2 -approximation as follows.
678
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
Fig. 3. Influence of the parameters δ0 and δ1 .
Recall that 1 ε=
Bn Pn − Bm Qm 2 dt
(12)
0 f
with Qm = Qcm ∪ Qm . By (6) and (11), we can express the control points qi , i = 1, . . . , m − 1, in quadratic functions of the parameters δv , rather than constants. After replacing them into (12), the objective function forms a quartic polynomial with two parameters. Since δv2 0, the coincidence of the oriented tangents is always satisfied. However, the approximating curve will be singular at the endpoint when δv = 0 or δv is nearly equal to 0, which is obviously undesirable in shape designing. For this reason, we adjust (12) by adding the regularization terms:
1
ε =
2 2 Bn Pn − Bm Qm 2 dt + λ 1 − δ02 + 1 − δ12 ,
(13)
0
with λ = σc AEL(Pn ), where AEL(Pn ) denotes the average edge length of the control polygon Pn and σc is a small positive number. In (13), the first term represents the L2 -error and the last two regulate δ0 and δ1 , respectively. The regularization terms will reach the minimum when δ02 = δ12 = 1. The parameter σc is a balance factor between the L2 -error and the regularization terms. So σc is adjustable and different values lead to different minima. If one decreases the value of σc , the approximating curve will have a smaller L2 -error; on the contrary, the L2 -error will become bigger since δv2 are forced to converge to 1. If σc is large enough, it will just lead to the C 1 approximation. Therefore, we can solve the singularity problem by minimizing (13). Since our goal is to achieve the optimal approximation, we set it to be 10−4 , which is based on many experiments. Due to the nonlinear minimization, it is impossible to obtain an explicit solution for the approximation. This minimization can be accomplished by many numerical methods (e.g., (Press et al., 2002, Chapter 10)), which usually involve a number of iterations. Since the gradient information can be easily calculated from (12) and (13), we use the conjugate gradient method to solve it. And a “good” initial guess which suffices in most cases is to set δv = 1, which is exactly the condition of C 1 -continuity. After replacing all the parameters δv in Qm with the values solved above, we obtain the multiple degree reduced approximating curve Q(t) which preserves G1 -continuity at the endpoints. We now summarize the algorithm for degree reduction as follows.
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
679
Fig. 4. Degree reduction by Algorithm 1 (from degree 7 to degree 4). Solid: without regularization terms; dashed: with regularization terms.
Algorithm 1. Input: pi Output: qi , ε 1. 2. 3. 4.
Set δ0 and δ1 with the initial values 1; Express q1 and qm−1 by (6) and qi , i = 2, . . . , m − 2, by (11); Use the conjugate gradient method to obtain δ0 and δ1 by minimizing the objective function (12) or (13); Compute qi , i = 1, . . . , m − 1, by (6) and (11) and the approximation error ε by (12).
Fig. 4 illustrates the improvement of the singularities by the regularization terms. The approximating curve is nearly singular at the left endpoint by using the objective function (12), which is eliminated by adding the regularization terms to the objective function. And the approximation effect still remains very good. Remark 1. We use the balance factor σc in the objective function (13) to avoid the singularities. If one obtains a regular curve by the objective function (12), σc can be set to 0. However, it may be a burden to regulate it manually. So we set σc = 10−4 . Remark 2. When δv are solved by the objective function (13), ε is exactly the squared L2 -error of C 1 approximation at the beginning of the iterations. Therefore, the curve obtained by Algorithm 1 will have a smaller L2 -error than traditional methods, no matter what value the parameter σc is. Remark 3. If one subdivides the given curve of degree n at an arbitrary interior point ts ∈ (0, 1)\{0.5} and transforms the parametric domains [0, ts ] of the left part and [ts , 1] of the right part both to [0, 1], the resulting curves will be Gn continuous at the breakpoint (the terminal of curve segments). The reason is that the derivatives of the composite function g = f ◦ φ are different for different functions φ and that the first derivatives of the above two linear transformations are two distinct constants. Obviously, C α -constrained degree reduction will generate curve segments with Gα -continuity at the breakpoints and with C α -continuity in the interior since degree reduction is dealt with over the domain [0, 1]. This result can also be fulfilled by Gα -constrained degree reduction, but with the advantage of having a smaller approximation error. Therefore, better approximation can be obtained if the curve is subdivided at the feature points, e.g., inflection points, or points of extreme curvature. However, if the curve is subdivided at ts = 0.5, the continuity at the breakpoints will be C k and Gk for the above two kinds of degree reduction, respectively.
680
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
4. G2 -constrained degree reduction 4.1. G2 condition Since G2 -continuity implies G1 -continuity, the condition (5) is necessary to be satisfied. Additionally, it has to satisfy (14) Q (v) = P φ(v) , v = 0, 1. Applying the chain rule to (14) leads to 2 Q (v) = φ (v) P (v) + φ (v)P (v),
v = 0, 1.
(15)
After relating (15) with the control points, we get n(n − 1) n (δ0 )4 2 p0 + η0 p0 , m(m − 1) m(m − 1) n(n − 1) n qm−2 = 2qm−1 − qm + (δ1 )4 2 pm−2 + η1 pm−1 , m(m − 1) m(m − 1)
q2 = 2q1 − q0 +
(16)
where ηv = φ (v), v = 0, 1. From now on, we assume that G2 condition is equivalent to the combination of the conditions (6) and (16). From (6), q1 (resp. qm−1 ) is a quadratic function of the parameter δ0 (resp. δ1 ). From (16), q2 (resp. qm−2 ) is quartic of the parameter δ0 (resp. δ1 ) and linear of the parameter η0 (resp. η1 ). Thus, it provides four additional parameters to optimize the approximation. 4.2. Algorithm for G2 -constrained degree reduction By G2 condition, we consider the G2 -constrained degree reduction, which is very similar to that described in Section 3.2. So we only outline the main framework and omit the details. f Let Qcm = (q0 , q1 , q2 , qm−2 , qm−1 , qm )T and Qm denotes other control points of Qm . Similar to (11), we obtain f −1 f f f Qm = Qm (δ0 , η0 , δ1 , η1 ) = Mm,m Mm,n Pn − Mcm,m Qcm . (17) f
Then, we replace Qcm and Qm into (12) for the general case and (13) based on the regularization terms, respectively. Therefore, the objective function (12) or (13), which contains four unknown parameters δv and ηv , forms a polynomial of degree 8. Again, we use the conjugate gradient method to solve the four parameters. A “good” initial guess which suffices in most cases is to set δv = 1 and ηv = 0, which is exactly the condition of C 2 -continuity. After replacing all the four parameters in Qm with the values solved above, we obtain the optimal approximating curve Q(t) which preserves G2 -continuity at the endpoints. The algorithm for G2 -constrained degree reduction is very similar to Algorithm 1 described in Section 3.2, so we omit it for simplicity. 5. Examples We now show some examples for the algorithms discussed in Sections 3 and 4 and make a comparison with traditional methods. Example 1. (Example 2 in (Ahn et al., 2004)) We consider a planar quintic Bézier curve with control points given by (0, 0), (0.2, 1), (0.4, 4), (0.6, 2), (0.8, 5), (1, 0). We want to find quartic and cubic Bézier curves to approximate it with G1 -continuity. Fig. 5 compares G1 -constrained degree reduction (by Algorithm 1) with C 1 -constrained one (Ahn et al., 2004), with the quartic and cubic approximations shown in Figs. 5(a) and 5(b), respectively. It is clearly seen that our method
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
(a)
681
(b)
Fig. 5. Degree reduction of the quintic Bézier curve (solid). Dashed: our method; dotted: the method in (Ahn et al., 2004). (a) Degree 5 to degree 4. (b) Degree 5 to degree 3.
(a)
(b)
Fig. 6. Degree reduction of the Bézier curve (solid) from degree 15 to degree 5. Dashed: our method; dotted: the method in (Ahn et al., 2004). (a) Without subdivision. (b) With one subdivision.
Fig. 7. Reduction from degree 10 to degree 6 with G2 -continuity. Solid: the original curve; dashed: our method; dotted: the method in (Ahn et al., 2004).
682
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
approximates the whole curve better, due to the G1 -continuity nature. Especially, the cubic approximation by our method is nearly the same as the quartic approximation. Example 2. We give a planar Bézier curve of degree fifteen which is the out-lines of the font “S”. The control points are (0, 0), (1.5, −2), (4.5, −1), (9, 0), (4.5, 1.5), (2.5, 3), (0, 5), (−4, 8.5), (3, 9.5), (4.4, 10.5), (6, 12), (8, 11), (9, 10), (9.5, 5), (7, 6), (5, 7). Fig. 6 illustrates multi-degree reduction of the complex curve. In Fig. 6(a), our method represents more features of the original curve. To improve the approximation effect, we subdivide the curve at ts = 0.5 and show the result in Fig. 6(b). Better approximation can be obtained through curve subdivisions. And we find that after one or two subdivisions, the given curve can be well approximated by piecewise continuous curves. Example 3. We consider multi-degree reduction with G2 -continuity, see Fig. 7. By the algorithm in Section 4.2, we obtain the optimal approximating curve in L2 -norm. It is clear that our method has better approximation effect. 6. Conclusion and future work In this paper we have introduced a new framework for multi-degree reduction of Bézier curves with G2 -continuity and have obtained the optimal approximation in L2 -norm. So the position, tangent direction and curvature are preserved at each endpoint. Due to the geometric continuity, our method has a smaller approximation error than traditional methods. We can also investigate Gα -constrained degree reduction (α 3) in the same way as described in Section 4. It is obvious that if a curve of degree n is degree elevated from degree m, the degree m approximating curve is exactly itself whatever method one uses. We assume that this case is not concerned in degree reduction. From our results it is easy to make a general conjecture: if one considers the problem of degree reduction with Gα -continuity, the approximation error will be strictly smaller than that obtained by considering it with C α -continuity. However, to prove such a conjecture seems to be a difficult task. Note that we use the additional parameters provided by geometric continuity. However, the solutions are obtained by numerical methods. So the theoretical analysis on the properties of Gα -constrained degree reduction is currently unavailable. In future work, we wish to tackle this problem with other methods. Another research direction is to develop degree reduction method which is independent on the parameterizations of the original curve. Perhaps one needs to consider the geometric distance between the curves. A method of this kind has been presented by Pottmann et al. (2002). Acknowledgements The authors thank the anonymous referees for their helpful suggestions and comments. This work is supported by the Natural Science Foundation of China (No. 60473130) and Foundation of State Key Basic Research 973 Development Programming Item of China (No. G2004CB318000). References Ahn, Y.J., 2003. Degree reduction of Bézier curves with C k -continuity using Jacobi polynomials. Computer Aided Geometric Design 20, 423–434. Ahn, Y.J., Lee, B.G., Park, Y., Yoo, J., 2004. Constrained polynomial degree reduction in the L2 -norm equals best weighted Euclidean approximation of Bézier coefficients. Computer Aided Geometric Design 21, 181–191. de Boor, C., Höllig, K., Sabin, M., 1987. High accuracy geometric Hermite interpolation. Computer Aided Geometric Design 4, 269–278. Chen, G., Wang, G., 2002. Optimal multi-degree reduction of Bézier curves with constraints of endpoints continuity. Computer Aided Geometric Design 19, 365–377. Degen, W.L.F., 2005. Geometric Hermite interpolation – In memoriam Josef Hoschek. Computer Aided Geometric Design 22, 573–592. Eck, M., 1995. Least squares degree reduction of Bézier curves. Computer-Aided Design 27 (11), 845–851. Farin, G., 1983. Algorithms for rational Bézier curves. Computer-Aided Design 15 (2), 73–79. Forrest, A.R., 1972. Interactive interpolation and approximation by Bézier curve. The Computer Journal 15 (1), 71–79. Pottmann, H., Leopoldseder, S., Hofer, M., 2002. Approximation with active B-spline curves and surfaces. In: Coquillart, S., Shum, H., Hu, S.M. (Eds.), Proc. of Pacific Graphics 2002. IEEE Press, Los Alamitos, CA, pp. 8–25.
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
683
Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 2002. Numerical Recipes in C++: The Art of Scientific Computing, second ed. Cambridge University Press, Cambridge, pp. 398–460. Szegö, G., 1975. Orthogonal Polynomials, fourth ed. American Mathematical Society, Providence, RI. Watkins, M.A., Worsey, A.J., 1988. Degree reduction of Bézier curves. Computer-Aided Design 20 (7), 398–405. Zheng, J., Wang, G., 2003. Perturbing Bézier coefficients for best constrained degree reduction in the L2 -norm. Graphical Models 65, 351–368.
Optimal multi-degree reduction of Bézier curves with G2 -continuity Lizheng Lu ∗ , Guozhao Wang Institute of Computer Graphics and Image Processing, Department of Mathematics, Zhejiang University, Hangzhou 310027, China Received 22 November 2005; received in revised form 7 September 2006; accepted 16 September 2006 Available online 30 October 2006
Abstract In this paper we present a novel approach to consider the multi-degree reduction of Bézier curves with G2 -continuity in L2 norm. The optimal approximation is obtained by minimizing the objective function based on the L2 -error between the two curves. In contrast to traditional methods, which typically consider the components of the curve separately, we use geometric information on the curve to generate the degree reduction. So positions, tangents and curvatures are preserved at the two endpoints. For avoiding the singularities at the endpoints, regularization terms are added to the objective function. Finally, numerical examples demonstrate the effectiveness of our algorithms. © 2006 Elsevier B.V. All rights reserved. Keywords: Degree reduction; Bézier curves; Optimal approximation; G2 -continuity; L2 -norm
1. Introduction Degree reduction of polynomial curves and surfaces is a common process in CAGD. It consists of approximating a polynomial by another one of lower degree. This process is of great importance in geometric modelling, such as data exchange, data compression and data comparison. For example, degree reduction is needed when data are transferred from one modelling system to another and these systems have different limitations on the maximum degree of polynomials. Furthermore, it can also be used to generate a piecewise continuous lower degree approximation to a given curve or surface so as to simplify some geometric or graphical algorithms like intersection calculation or rendering. 1.1. Related work There have been many methods developed for degree reduction. Forrest (1972) and Farin (1983) considered the inverse of degree elevation and obtained two sets of control points. This approach is easy to compute, but not optimal in the usual Lp -norm for any p = 1, 2, . . . , +∞. Since degree reduction is an approximation problem in nature, methods in the classical approximation theory (Szegö, 1975) can be employed. In particular, the optimal approximations with respect to the L∞ or L2 metric are of interest. Watkins and Worsey (1988) used Chebyshev economization to produce * Corresponding author.
E-mail addresses: [email protected] (L. Lu), [email protected] (G. Wang). 0167-8396/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cagd.2006.09.002
674
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
the best L∞ -approximation of degree n − 1 to a given degree n polynomial. This best approximation, however, does not interpolate the given curve at the endpoints. The endpoint constraints are frequently required in many applications, especially when degree reduction is combined with subdivision to generate piecewise continuous approximations. Recently, C k -constrained optimal degree reduction has been presented. Eck (1995) used constrained Legendre polynomials to minimize the L2 -norm between the two curves. Ahn (2003) presented one degree reduction in L∞ norm. And multi-degree reduction at one time avoiding stepwise computing was investigated in (Ahn et al., 2004; Chen and Wang, 2002; Zheng and Wang, 2003), which showed that the optimal approximation in L2 -norm can be obtained by different methods and the results are equivalent. For every Lp -norm there exists an optimal approximation, but which one is the best is an open problem. In the above mentioned papers the approximation error is measured between the parameterized curves. Pottmann et al. (2002) applied active curves to degree reduction, thus minimizing the geometric distance between the curves. The key idea is to iteratively change the control points of the active curve so that the active curve deforms towards the given curve. And their method shows the better approximation obtained by the geometric distance. 1.2. Motivation and outline of the present paper Traditional methods for degree reduction are applied to parametric curves simply by matching the position and derivatives at the same parameter values, in general f (r) (t) = g (r) (t),
r = 0, 1, . . . , α.
Since parametric representations of curves are not unique, they will produce different approximations for a given curve. For example, Fig. 1 shows a quartic polynomial (solid): F(x) = (x, 2x − 6x 2 + 5x 3 − x 4 ),
x ∈ [0, 1],
and its cubic best L2 -approximation with C 1 -continuity (dotted) by the method in (Ahn et al., 2004). However, if one reparameterizes F(x) by taking x = φ(t) = (t + t 2 )/2, the curve F(φ(t)) remains unchanged but with a different approximating curve (dashed). The corresponding approximations with G1 -continuity (with “+”) by our method are also shown in Fig. 1. Our main motivation is to consider degree reduction with geometric continuity, which is less dependent on parameterizations. And such method gives better approximation for degree reduction. Geometric continuity was introduced by de Boor et al. (1987) to curve interpolation, which is called Geometric Hermite Interpolation (GHI). They showed that if the curvature at one endpoint is not vanished, a planar curve can be interpolated by cubic spline with G2 -continuity and the approximation order is 6. For more details about GHI, see the recent survey by Degen (2005).
Fig. 1. A comparison of C 1 - and G1 -constrained degree reduction.
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
675
In this paper, we consider multi-degree reduction of Bézier curves with geometric continuity by minimizing the objective function based on the L2 -error between the two curves. For G1 -continuity, it provides two additional parameters compared to C 1 approximation. And for G2 -continuity, the additional parameters are four. By using these parameters, we can optimize the approximation and obtain the optimal approximating curve in L2 -norm. The outline of this paper is as follows. We briefly characterize the problem of Gα -constrained multi-degree reduction in Section 2. In Section 3, we consider G1 -constrained degree reduction and propose the algorithm. We then outline G2 -constrained degree reduction in Section 4. In Section 5 we provide some examples and compare our method with traditional methods. Finally, Section 6 concludes this paper. 2. Preliminaries 2.1. Definitions and notations the space of parametric polynomials in Rd with degree at most n and · denotes the In this paper, Πnd denotes √ Euclidean norm v = v, v. A Bézier curve of degree n is defined by the control points pi ∈ Rd in the form P(t) =
n
Bin (t)pi ,
0 t 1,
i=0
where Bin (t) are the Bernstein polynomials given by Bin (t) = operator as follows, pi = pi+1 − pi ,
r pi = (r−1 pi ),
n n−i t i . And we define the forward difference i (1 − t)
r 2.
The product of two Bernstein polynomials is given by mn i
j
m+n (t) Bim (t)Bjn (t) = m+n Bi+j
(1)
i+j
and the integral is 1 Bin (t) dt = 0
1 . n+1
(2)
We introduce the following lemma, where mij can be easily obtained from (1) and (2). Lemma 1. Let Mm,n = (mij ) be an (m + 1) × (n + 1) matrix with the elements given by 1 mij =
Bim (t)Bjn (t) dt, 0
then Mm,m is a real symmetric positive definite matrix. Proof. For any nonzero column vector ξ with m + 1 elements, we have 2 1 m m Bi (t)ξi dt = ξ T Mm,m ξ > 0. 2 0
i=0
2.2. Problem statement We start with the problem of Gα -constrained degree reduction. Given a degree n Bézier curve P(t) =
n i=0
Bin (t)pi ,
(3)
676
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
find a degree m (m < n) Bézier curve Q(t) =
m
Bim (t)qi
i=0
such that the following two conditions are satisfied: (i) P(t) and Q(t) are Gα -continuous (2α m − 1) at t = 0, 1, i.e., Q(r) (v) = P(r) φ(v) , v = 0, 1, r = 0, 1, . . . , α,
(4)
where φ : [0, 1] → [0, 1] is an strictly increasing function with φ(0) = 0 and φ(1) = 1; 1 (ii) Q(t) minimizes the squared L2 -error ε = 0 P(t) − Q(t)2 dt for all possible polynomials in Πmd . Note that traditional methods consider degree reduction of Bézier curves with the constraints of C α -continuity, which fix the first α + 1 and the last α + 1 control points of the approximating curve. In contrast, Gα -continuity will bring additional parameters to optimize the approximation. By using these parameters, we can obtain the approximating curve with a smaller approximation error. However, it turns out to be a nonlinear problem and brings more difficulties to degree reduction when α 1. 3. G1 -constrained degree reduction 3.1. G1 condition Clearly, for G0 -continuity, the endpoints of Q(t) should coincide with the endpoints of P(t). And for G1 -continuity, the coincidence of the oriented tangents is additionally needed. From (4), we get Q (0) = φ (0)P (0) Denoting
φ (0) = δ02
and
and Q (1) = φ (1)P (1). φ (1) = δ12 ,
(5) G1
we can easily relate condition (5) with the control points. More precisely, n 2 q0 = p0 , q1 = p0 + δ0 (p1 − p0 ), m (6) n qm−1 = pn − δ12 (pn − pn−1 ). qm = pn , m Noting that geometric continuity does not depend on the chosen parameterization, the control point q1 can thus move along the direction − p−0−p→1 without violating G1 condition (see Fig. 2), and the same holds for qm−1 . Thus it provides two additional parameters to optimize the shape of the degree reduced curve. When replacing G1 -continuity with C 1 -continuity, q1 and qm−1 are uniquely determined by δv = 1. Therefore, 1 C approximation is a special case of G1 approximation. We can also imagine that G1 approximation will lead to a smaller error between the original curve and the approximating one. We will propose the algorithm in Section 3.2. 3.2. Algorithm for G1 -constrained degree reduction For a given degree n Bézier curve P(t), the problem of G1 -constrained degree reduction can be solved through two stages. In the first stage, we construct a degree m Bézier curve Q(t) interpolating P(t) according to (6). More precisely, it can be expressed as Q(t) = B0m (t)q0 + B1m (t)q1 +
m−2
m m Bim (t)qi + Bm−1 (t)qm−1 + Bm (t)qm ,
(7)
i=2
where q1 and qm−1 contain the unknown parameters δ0 and δ1 , respectively. We assume that the free parameters δ0 and δ1 are temporarily fixed, then solve the interior control points qi , i = 2, . . . , m − 2, by minimizing 1 ε= 0
P(t) − Q(t) 2 dt.
(8)
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
677
Fig. 2. G1 condition for degree reduction.
Denoting Bn = (B0n (t), . . . , Bnn (t)) and Pn = (p0 , . . . , pn )T , we have 1 ε= 0
1 =
P(t) − Q(t) 2 dt =
1 Bn Pn − Bm Qm 2 dt 0
Bn Pn − Bc Qc − Bfm Qfm 2 dt, m m
(9)
0 f
where Qcm = (q0 , q1 , qm−1 , qm )T and Qm denotes other control points of Qm . f For a minimum of (9), it is necessary that the derivatives of ε with respect to the elements of Qm are zero, which gives 0=
1 ∂ε f f f = Mm,n Pn − Mcm,m Qcm − Mm,m Qm , 2 ∂Qfm
(10)
where f
Mm,n := Mm,n (2, . . . , m − 2; 0, . . . , n), Mcm,m := Mm,m (2, . . . , m − 2; 0, 1, m − 1, m), f
Mm,m := Mm,m (2, . . . , m − 2; 2, . . . , m − 2). The notation A(. . . ; . . .) denotes the submatrix of the matrix A obtained by extracting the specific rows and columns. f Since Mm,m is a real symmetric positive definite matrix (see Lemma 1), Mm,m is invertible. Therefore, ε is minimized by choosing f −1 f f f Qm = Qm (δ0 , δ1 ) = Mm,m Mm,n Pn − Mcm,m Qcm . (11) Note that q1 and qm−1 are quadratic functions of the parameters δ0 and δ1 , respectively. Therefore, from (11), the interior control points qi , i = 2, . . . , m − 2, are quadratic functions of the parameters δ0 and δ1 . As shown by Fig. 3, we can obtain different G1 -constrained approximating curves to the curve F(x) for different parameter values δ0 and δ1 . The second stage is then to determine the two free parameters δv such that the approximation error (8) is minimized. We will propose how to optimize the two parameters so as to obtain the optimal L2 -approximation as follows.
678
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
Fig. 3. Influence of the parameters δ0 and δ1 .
Recall that 1 ε=
Bn Pn − Bm Qm 2 dt
(12)
0 f
with Qm = Qcm ∪ Qm . By (6) and (11), we can express the control points qi , i = 1, . . . , m − 1, in quadratic functions of the parameters δv , rather than constants. After replacing them into (12), the objective function forms a quartic polynomial with two parameters. Since δv2 0, the coincidence of the oriented tangents is always satisfied. However, the approximating curve will be singular at the endpoint when δv = 0 or δv is nearly equal to 0, which is obviously undesirable in shape designing. For this reason, we adjust (12) by adding the regularization terms:
1
ε =
2 2 Bn Pn − Bm Qm 2 dt + λ 1 − δ02 + 1 − δ12 ,
(13)
0
with λ = σc AEL(Pn ), where AEL(Pn ) denotes the average edge length of the control polygon Pn and σc is a small positive number. In (13), the first term represents the L2 -error and the last two regulate δ0 and δ1 , respectively. The regularization terms will reach the minimum when δ02 = δ12 = 1. The parameter σc is a balance factor between the L2 -error and the regularization terms. So σc is adjustable and different values lead to different minima. If one decreases the value of σc , the approximating curve will have a smaller L2 -error; on the contrary, the L2 -error will become bigger since δv2 are forced to converge to 1. If σc is large enough, it will just lead to the C 1 approximation. Therefore, we can solve the singularity problem by minimizing (13). Since our goal is to achieve the optimal approximation, we set it to be 10−4 , which is based on many experiments. Due to the nonlinear minimization, it is impossible to obtain an explicit solution for the approximation. This minimization can be accomplished by many numerical methods (e.g., (Press et al., 2002, Chapter 10)), which usually involve a number of iterations. Since the gradient information can be easily calculated from (12) and (13), we use the conjugate gradient method to solve it. And a “good” initial guess which suffices in most cases is to set δv = 1, which is exactly the condition of C 1 -continuity. After replacing all the parameters δv in Qm with the values solved above, we obtain the multiple degree reduced approximating curve Q(t) which preserves G1 -continuity at the endpoints. We now summarize the algorithm for degree reduction as follows.
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
679
Fig. 4. Degree reduction by Algorithm 1 (from degree 7 to degree 4). Solid: without regularization terms; dashed: with regularization terms.
Algorithm 1. Input: pi Output: qi , ε 1. 2. 3. 4.
Set δ0 and δ1 with the initial values 1; Express q1 and qm−1 by (6) and qi , i = 2, . . . , m − 2, by (11); Use the conjugate gradient method to obtain δ0 and δ1 by minimizing the objective function (12) or (13); Compute qi , i = 1, . . . , m − 1, by (6) and (11) and the approximation error ε by (12).
Fig. 4 illustrates the improvement of the singularities by the regularization terms. The approximating curve is nearly singular at the left endpoint by using the objective function (12), which is eliminated by adding the regularization terms to the objective function. And the approximation effect still remains very good. Remark 1. We use the balance factor σc in the objective function (13) to avoid the singularities. If one obtains a regular curve by the objective function (12), σc can be set to 0. However, it may be a burden to regulate it manually. So we set σc = 10−4 . Remark 2. When δv are solved by the objective function (13), ε is exactly the squared L2 -error of C 1 approximation at the beginning of the iterations. Therefore, the curve obtained by Algorithm 1 will have a smaller L2 -error than traditional methods, no matter what value the parameter σc is. Remark 3. If one subdivides the given curve of degree n at an arbitrary interior point ts ∈ (0, 1)\{0.5} and transforms the parametric domains [0, ts ] of the left part and [ts , 1] of the right part both to [0, 1], the resulting curves will be Gn continuous at the breakpoint (the terminal of curve segments). The reason is that the derivatives of the composite function g = f ◦ φ are different for different functions φ and that the first derivatives of the above two linear transformations are two distinct constants. Obviously, C α -constrained degree reduction will generate curve segments with Gα -continuity at the breakpoints and with C α -continuity in the interior since degree reduction is dealt with over the domain [0, 1]. This result can also be fulfilled by Gα -constrained degree reduction, but with the advantage of having a smaller approximation error. Therefore, better approximation can be obtained if the curve is subdivided at the feature points, e.g., inflection points, or points of extreme curvature. However, if the curve is subdivided at ts = 0.5, the continuity at the breakpoints will be C k and Gk for the above two kinds of degree reduction, respectively.
680
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
4. G2 -constrained degree reduction 4.1. G2 condition Since G2 -continuity implies G1 -continuity, the condition (5) is necessary to be satisfied. Additionally, it has to satisfy (14) Q (v) = P φ(v) , v = 0, 1. Applying the chain rule to (14) leads to 2 Q (v) = φ (v) P (v) + φ (v)P (v),
v = 0, 1.
(15)
After relating (15) with the control points, we get n(n − 1) n (δ0 )4 2 p0 + η0 p0 , m(m − 1) m(m − 1) n(n − 1) n qm−2 = 2qm−1 − qm + (δ1 )4 2 pm−2 + η1 pm−1 , m(m − 1) m(m − 1)
q2 = 2q1 − q0 +
(16)
where ηv = φ (v), v = 0, 1. From now on, we assume that G2 condition is equivalent to the combination of the conditions (6) and (16). From (6), q1 (resp. qm−1 ) is a quadratic function of the parameter δ0 (resp. δ1 ). From (16), q2 (resp. qm−2 ) is quartic of the parameter δ0 (resp. δ1 ) and linear of the parameter η0 (resp. η1 ). Thus, it provides four additional parameters to optimize the approximation. 4.2. Algorithm for G2 -constrained degree reduction By G2 condition, we consider the G2 -constrained degree reduction, which is very similar to that described in Section 3.2. So we only outline the main framework and omit the details. f Let Qcm = (q0 , q1 , q2 , qm−2 , qm−1 , qm )T and Qm denotes other control points of Qm . Similar to (11), we obtain f −1 f f f Qm = Qm (δ0 , η0 , δ1 , η1 ) = Mm,m Mm,n Pn − Mcm,m Qcm . (17) f
Then, we replace Qcm and Qm into (12) for the general case and (13) based on the regularization terms, respectively. Therefore, the objective function (12) or (13), which contains four unknown parameters δv and ηv , forms a polynomial of degree 8. Again, we use the conjugate gradient method to solve the four parameters. A “good” initial guess which suffices in most cases is to set δv = 1 and ηv = 0, which is exactly the condition of C 2 -continuity. After replacing all the four parameters in Qm with the values solved above, we obtain the optimal approximating curve Q(t) which preserves G2 -continuity at the endpoints. The algorithm for G2 -constrained degree reduction is very similar to Algorithm 1 described in Section 3.2, so we omit it for simplicity. 5. Examples We now show some examples for the algorithms discussed in Sections 3 and 4 and make a comparison with traditional methods. Example 1. (Example 2 in (Ahn et al., 2004)) We consider a planar quintic Bézier curve with control points given by (0, 0), (0.2, 1), (0.4, 4), (0.6, 2), (0.8, 5), (1, 0). We want to find quartic and cubic Bézier curves to approximate it with G1 -continuity. Fig. 5 compares G1 -constrained degree reduction (by Algorithm 1) with C 1 -constrained one (Ahn et al., 2004), with the quartic and cubic approximations shown in Figs. 5(a) and 5(b), respectively. It is clearly seen that our method
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
(a)
681
(b)
Fig. 5. Degree reduction of the quintic Bézier curve (solid). Dashed: our method; dotted: the method in (Ahn et al., 2004). (a) Degree 5 to degree 4. (b) Degree 5 to degree 3.
(a)
(b)
Fig. 6. Degree reduction of the Bézier curve (solid) from degree 15 to degree 5. Dashed: our method; dotted: the method in (Ahn et al., 2004). (a) Without subdivision. (b) With one subdivision.
Fig. 7. Reduction from degree 10 to degree 6 with G2 -continuity. Solid: the original curve; dashed: our method; dotted: the method in (Ahn et al., 2004).
682
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
approximates the whole curve better, due to the G1 -continuity nature. Especially, the cubic approximation by our method is nearly the same as the quartic approximation. Example 2. We give a planar Bézier curve of degree fifteen which is the out-lines of the font “S”. The control points are (0, 0), (1.5, −2), (4.5, −1), (9, 0), (4.5, 1.5), (2.5, 3), (0, 5), (−4, 8.5), (3, 9.5), (4.4, 10.5), (6, 12), (8, 11), (9, 10), (9.5, 5), (7, 6), (5, 7). Fig. 6 illustrates multi-degree reduction of the complex curve. In Fig. 6(a), our method represents more features of the original curve. To improve the approximation effect, we subdivide the curve at ts = 0.5 and show the result in Fig. 6(b). Better approximation can be obtained through curve subdivisions. And we find that after one or two subdivisions, the given curve can be well approximated by piecewise continuous curves. Example 3. We consider multi-degree reduction with G2 -continuity, see Fig. 7. By the algorithm in Section 4.2, we obtain the optimal approximating curve in L2 -norm. It is clear that our method has better approximation effect. 6. Conclusion and future work In this paper we have introduced a new framework for multi-degree reduction of Bézier curves with G2 -continuity and have obtained the optimal approximation in L2 -norm. So the position, tangent direction and curvature are preserved at each endpoint. Due to the geometric continuity, our method has a smaller approximation error than traditional methods. We can also investigate Gα -constrained degree reduction (α 3) in the same way as described in Section 4. It is obvious that if a curve of degree n is degree elevated from degree m, the degree m approximating curve is exactly itself whatever method one uses. We assume that this case is not concerned in degree reduction. From our results it is easy to make a general conjecture: if one considers the problem of degree reduction with Gα -continuity, the approximation error will be strictly smaller than that obtained by considering it with C α -continuity. However, to prove such a conjecture seems to be a difficult task. Note that we use the additional parameters provided by geometric continuity. However, the solutions are obtained by numerical methods. So the theoretical analysis on the properties of Gα -constrained degree reduction is currently unavailable. In future work, we wish to tackle this problem with other methods. Another research direction is to develop degree reduction method which is independent on the parameterizations of the original curve. Perhaps one needs to consider the geometric distance between the curves. A method of this kind has been presented by Pottmann et al. (2002). Acknowledgements The authors thank the anonymous referees for their helpful suggestions and comments. This work is supported by the Natural Science Foundation of China (No. 60473130) and Foundation of State Key Basic Research 973 Development Programming Item of China (No. G2004CB318000). References Ahn, Y.J., 2003. Degree reduction of Bézier curves with C k -continuity using Jacobi polynomials. Computer Aided Geometric Design 20, 423–434. Ahn, Y.J., Lee, B.G., Park, Y., Yoo, J., 2004. Constrained polynomial degree reduction in the L2 -norm equals best weighted Euclidean approximation of Bézier coefficients. Computer Aided Geometric Design 21, 181–191. de Boor, C., Höllig, K., Sabin, M., 1987. High accuracy geometric Hermite interpolation. Computer Aided Geometric Design 4, 269–278. Chen, G., Wang, G., 2002. Optimal multi-degree reduction of Bézier curves with constraints of endpoints continuity. Computer Aided Geometric Design 19, 365–377. Degen, W.L.F., 2005. Geometric Hermite interpolation – In memoriam Josef Hoschek. Computer Aided Geometric Design 22, 573–592. Eck, M., 1995. Least squares degree reduction of Bézier curves. Computer-Aided Design 27 (11), 845–851. Farin, G., 1983. Algorithms for rational Bézier curves. Computer-Aided Design 15 (2), 73–79. Forrest, A.R., 1972. Interactive interpolation and approximation by Bézier curve. The Computer Journal 15 (1), 71–79. Pottmann, H., Leopoldseder, S., Hofer, M., 2002. Approximation with active B-spline curves and surfaces. In: Coquillart, S., Shum, H., Hu, S.M. (Eds.), Proc. of Pacific Graphics 2002. IEEE Press, Los Alamitos, CA, pp. 8–25.
L. Lu, G. Wang / Computer Aided Geometric Design 23 (2006) 673–683
683
Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 2002. Numerical Recipes in C++: The Art of Scientific Computing, second ed. Cambridge University Press, Cambridge, pp. 398–460. Szegö, G., 1975. Orthogonal Polynomials, fourth ed. American Mathematical Society, Providence, RI. Watkins, M.A., Worsey, A.J., 1988. Degree reduction of Bézier curves. Computer-Aided Design 20 (7), 398–405. Zheng, J., Wang, G., 2003. Perturbing Bézier coefficients for best constrained degree reduction in the L2 -norm. Graphical Models 65, 351–368.
