Gk,l-constrained multi-degree reduction of B´ezier curves Przemysław Gospodarczyk∗, Stanisław Lewanowicz, Paweł Woźny Institute of Computer Science, University of Wrocław, ul. Joliot-Curie 15, 50-383 Wrocław, Poland Abstract We present a new approach to the problem of Gk,l-constrained (k,l 3) multi-degree reduction of B´ezier ≤ curves with respect to the least squares norm. First, to minimize the least squares error, we consider two methods of determining the values of geometric continuity parameters. One of them is based on quadratic andnonlinearprogramming,while the otheruses somesimplifying assumptionsandsolvesa systemoflinear equations. Next, for prescribed values of these parameters, we obtain control points of the multi-degree 5 1 reduced curve, using the properties of constrained dual Bernstein basis polynomials. Assuming that the 0 input and output curves are of degree n and m, respectively, we determine these points with the complexity 2 O(mn), which is significantly less than the cost of other known methods. Finally, we give several examples r to demonstrate the effectiveness of our algorithms. a M Keywords: Constrained dual Bernstein basis, B´ezier curves, Multi-degree reduction, Geometric continuity, Quadratic programming,Nonlinear programming. 9 ] R 1. Introduction G Let Πd denote the space of all parametric polynomials in Rd of degree at most n; Π1 :=Π . . n n n s A B´ezier curve P Πd of degree n N is the following parametric curve: c n ∈ n ∈ [ n 3 Pn(t):= piBin(t) (0≤t≤1), (1.1) v i=0 X 2 3 where p ,p ,...,p Rd are so-called control points, and Bn,Bn,...,Bn are the Bernstein polynomials of 0 degree n0giv1en by n ∈ 0 1 n 3 n 0 Bn(t):= ti(1 t)n−i (0 i n). (1.2) i i − ≤ ≤ 1. (cid:18) (cid:19) 0 In this paper, we consider the following problem. 5 1 Problem 1.1. [Gk,l-constrained multi-degree reduction] : For a given B´ezier curve P of degree n, find a B´ezier curve R of lower degree m, v n m i X m R (t):= r Bm(t) (0 t 1), (1.3) r m i i ≤ ≤ a i=0 X so that the following conditions are satisfied: (i) P and R are Gk,l-continuous ( 1 k,l 3 and k+l <m 1) at the endpoints, i.e., n m − ≤ ≤ − di di R (t)= P (ϕ(t)) (t=0; i=0,1,...,k), dti m dti n (1.4) dj dj  dtjRm(t)= dtjPn(ϕ(t)) (t=1; j =0,1,...,l),  where ϕ:[0,1] [0,1] is a strictly increasing function with ϕ(0)=0 and ϕ(1)=1; → ∗Correspondingauthor. Fax+48713757801 Email addresses: [email protected] (PrzemysławGospodarczyk), [email protected] (Stanisław Lewanowicz), [email protected] (PawełWoźny) (ii) value of the squared L -error 2 1 P R 2 := (1 t)αtβ P (t) R (t) 2dt (α, β > 1) || n− m||L2 − || n − m || − Z0 is minimized in the space Πd , where is the Euclidean vector norm. m ||·|| Problemsoftheabovetypehavebeenrecentlydiscussedinseveralpapers[6,9,10,11,12,16,17],usuallyun- dersimplifying assumptionsϕ′(0)=ϕ′(1)=1,whichimplied,forexample,the hybridC1,1/G2,2-constrained degree reduction, meaning that we impose constraints of C1,1-continuity, followed by G2,2-continuity, at the endpoints. Most of the known algorithms solve a system of normal equations to get control points of the multi-degree reduced curve (1.3). Consequently, solution depends on the inverse of a certain matrix, so the obtainedformulasarenottrulyexplicitandthecostofthemethodishigh(see,e.g.,[6,10,12]). Forextensive lists of references, see the recent papers of Lu [6], or Rababah and Mann [12]. The conventional problem of degree reduction differs from Problem 1.1 in considering, instead of condition (i), the Ck,l-continuity at the endpoints of curves, i.e., R(i)(0)=P(i)(0) (i=0,1,...,k), m n (1.5) (j) (j) Rm (1)=Pn (1) (j =0,1,...,l). ) Inthe past30years,manypapersdealingwiththis problemhavebeenpublished(see,e.g.,[2,3,13,14,15]). In particular, in [15], two of us have proposed a method based on the use of the so-called dual Bernstein polynomials, which has complexity O(mn), the least among the existing algorithms. In the present paper, we apply an extended version of this method as an essential part of the algorithms of solving Problem 1.1. Suchanapproachallowsusto avoidmatrixinversion. Assumingthat 1 k,l 3andincluding thehybrid − ≤ ≤ cases, there are 37 continuity cases which require computation of the continuity parameters. Those variants of the problem differ, and we have not proven that in each case a unique solution exists. The outline of the paper is as follows. Section 2 contains a preliminary material. In Section 3, we relate the Gk,l-continuity conditions with the control points of the curves P and R . Section 4 brings complete n m solutions of Problem 1.1, with and without the simplifying assumptions. Section 5 deals with algorithmic implementation of the proposed methods. In Section 6, we give some examples showing efficiency of our methods. Conclusions are given in Section 7. 2. Preliminaries In this section, we introduce necessary definitions and notation. We define the inner product , by α,β h· ·i 1 f,g := (1 t)αtβf(t)g(t)dt (α, β > 1). (2.1) α,β h i − − Z0 There is a unique dual Bernstein polynomial basis of degree n Dn,Dn,...,Dn Π , 0 1 n ∈ n associated with the basis (1.2), so that Dn,Bn =δ (i,j =0,1,...,n), i j α,β ij (cid:10) (cid:11) where δ equals 1 if i=j, and 0 otherwise. ij Given the integers k,l such that k,l 1 and k+l<n 1, let Π(k,l) be the space of all polynomials of n ≥− − degree at most n, whose derivatives of orders 0,1,...,k at t=0, as well as derivatives of orders 0,1,...,l at t = 1, vanish. We use the convention that derivative of order 0 of a function is the function itself. Clearly, dim Π(k,l) = n k l 1, and the Bernstein polynomials Bn ,Bn ,...,Bn form a basis of this n − − − k+1 k+2 n−l−1 space. There is a unique dual constrained Bernstein polynomial basis of degree n (cid:8) (cid:9) D(n,k,l),D(n,k,l),...,D(n,k,l) Π(k,l) k+1 k+2 n−l−1 ⊂ n n o satisfying the relation D(n,k,l),Bn = δ (i,j = k + 1,k + 2,...,n l 1). Obviously, we have i j α,β ij − − D(n,−1,−1) = Dn, whichDcorrespondsEto the case without any constraints. For properties of the polynomials i i Bn and Dn, see [15, Appendix A]. i i 2 Forward difference operator is given by ∆0q :=q , ∆kq :=∆k−1q ∆k−1q (k =1,2,...). i i i i+1 i − We use Cp,q/Gk,l notation to describe the hybrid constraints, where p,q ,1 and (k 2 or l 2). In the case of k 2 and p = 1, we set ϕ′(0) := 1. Similarly, for l 2 and q =∈ {1,−we}set ϕ′(1)≥:= 1. Se≥tting p:= , q := m≥eans that we do not fix ϕ′(0), ϕ′(1), respectively.≥Clearly, C−,−/Gk,l denotes Gk,l. − − 3. Geometric continuity In this section, we relate the Gk,l-continuity conditions (1.4) with the control points. We limit ourselves to k,l 3 cases, which are the most important from a practical point of view. ≤ Remark that the control points r , r ,..., r depend on the parameters 1 2 k λ :=ϕ(j)(0) (j =1,2,...,k), j while the points rm−1,rm−2,...,rm−l depend on µ :=ϕ(j)(1) (j =1,2,...,l). j Now, let us recall the well known formulas (see [7], also [6, 11]). When k =3, we have: n r =p , r =p + λ ∆p , (3.1) 0 0 1 0 1 0 m n 1 (n 1) r =p + 2λ + λ ∆p + − 2λ2∆2p , (3.2) 2 0 m 1 m 1 2 0 (m 1) 1 0 (cid:20) − (cid:21) − 2 n 3 1 r =p + 3λ + λ + λ ∆p 3 0 1 2 3 0 m m 1 (m 2) (cid:20) − − 2 (cid:21) (n 1) 1 (n 2) +3 − 2 λ2+ λ λ ∆2p + − 3 λ3∆3p . (3.3) (m 1) 1 m 2 1 2 0 (m 2) 1 0 − 2 (cid:20) − (cid:21) − 3 In the case of k = 2, we use (3.1) and (3.2). For k = 1, formulas (3.1) hold. Analogously, when l = 3, we have: n rm =pn, rm−1 =pn µ1∆pn−1, (3.4) − m n 1 (n 1) rm−2 =pn− m 2µ1− m 1µ2 ∆pn−1+ (m− 1)2µ21∆2pn−2, (3.5) (cid:20) − (cid:21) − 2 n 3 1 rm−3 =pn 3µ1 µ2+ µ3 ∆pn−1 − m − m 1 (m 2) (cid:20) − − 2 (cid:21) (n 1) 1 (n 2) +3(m− 1)2 µ21− m 2µ1µ2 ∆2pn−2− (m− 2)3 µ31∆3pn−3. (3.6) − 2 (cid:20) − (cid:21) − 3 In the case of l=2, we use (3.4) and (3.5). For l =1, formulas (3.4) hold. 4. Gk,l-constrained multi-degree reduction problem 4.1. Multi-degree reduction of B´ezier curves with prescribed boundary control points First, we discuss the following model problem of constrained multi-degree reduction: Problem 4.1. [Multi-degree reduction with prescribed boundary control points] Given a B´ezier curve P Πd, n ∈ n n P (t):= p Bn(t), n i i i=0 X 3 we look for a B´ezier curve R Πd (m<n), m ∈ m m R (t):= r Bm(t), (4.1) m i i i=0 X having the prescribed control points r0,r1,...,rk and rm−l,rm−l+1,...,rm, that gives minimum value of the error 1 E(α,β) := P R 2 = (1 t)αtβ P (t) R (t) 2dt (α, β > 1). (4.2) || n− m||L2 − || n − m || − Z0 Giventhepointsp :=(p ,p ...,p ) Rd(i=0,1,...,n)andr :=(r ,r ...,r ) Rd(i=0,1,...,m), i i1 i2 id i i1 i2 id ∈ ∈ we use notation ph, rh for the vectors of hth coordinates of the points p ,p ,...,p and r ,r ,...,r , 0 1 n 0 1 m respectively: ph :=[p ,p ,...,p ], rh :=[r ,r ,...,r ] (h=1,2,...,d). 0h 1h nh 0h 1h mh As an extension of the result given in [15] (see also [5]), we obtain the following theorem. Theorem 4.2. The inner control points r = (r ,r ,...,r ) (k+1 i m l 1) of the curve (4.1), i i1 i2 id ≤ ≤ − − being the solution of the Problem 4.1, are given by n r = υ φ (i=k+1,k+2,...,m l 1), (4.3) i j ij − − j=0 X where φ := Bn,D(m,k,l) , (4.4) ij j i α,β D E and −1 k m n n m m υ :=p + − r (j =0,1,...,n). (4.5) j j h −(cid:18)j(cid:19) h=0 h=m−l!(cid:18)j−h(cid:19)(cid:18)h(cid:19) X X The squared L -error (4.2) is given by 2 d E(α,β) = I (ph,ph)+I (rh,rh) 2I (ph,rh) , (4.6) nn mm nm − h=1 X(cid:8) (cid:9) where for a:=[a ,a ,...,a ] and b:=[b ,b ,...,b ], we define 0 1 N 0 1 M N M B(α+1,β+1) N M INM(a,b):= (α+1)N+M−i−j(β+1)i+jaibj, (α+β+2) i j N+M i=0j=0(cid:18) (cid:19)(cid:18) (cid:19) XX where B(α,β):= Γ(α)Γ(β) is the beta function. Γ(α+β) Proof. Let us write R (t)=S (t)+T (t), m m m where m−l−1 k m S (t):= r Bm(t), T (t):= + r Bm(t). m i i m i i ! i=k+1 i=0 i=m−l X X X Using the degree elevation formula (see, e.g., [4, 6.10]; we adopt the usual convention that u =0 if v <0 § v or v >u) n −1 (cid:0) (cid:1) m n m n Bm(t)= − Bn(t), i i h i h h (cid:18) (cid:19)h=0(cid:18) − (cid:19)(cid:18) (cid:19) X we write n T (t)= d Bn(t), m j j j=0 X 4 where −1 k m n n m m d := + − r . j h j j h h (cid:18) (cid:19) h=0 h=m−l!(cid:18) − (cid:19)(cid:18) (cid:19) X X Now, we observe that d 1 P R 2 = W S 2 = (1 t)αtβ Wh(t) Sh(t) 2dt, || n− m||L2 || n− m||L2 − n − m h=1Z0 X (cid:2) (cid:3) where n W (t):= W1(t),W2(t),...,Wd(t) =P (t) T (t)= υ Bn(t), n n n n n − m i i i=0 (cid:2) (cid:3) X S (t):= S1(t),S2 (t),...,Sd(t) , m m m m with (cid:2) (cid:3) υ :=p d . i i i − Thus, we are looking for the best least squares approximation for Wh (h = 1,2,...,d) in the space Π(k,l). n m Remembering that Bm and D(m,k,l) (k+1 i m l 1) are the dual bases in the space Π(k,l), we obtain i i ≤ ≤ − − m n n r = υ Bn,D(m,k,l) = υ φ (i=k+1,k+2,...,m l 1), i j j i α,β j ij − − Xj=0 D E Xj=0 which is the formula (4.3). Proof of (4.6) uses an argument similar to the one given in [15]. (cid:3) Remark 4.3. Letus define the quantities ψ (i=k+1,k+2,...,m l 1; j =0,1,...,n), relatedto the ij − − coefficients φ (cf. (4.4)) by the following formula: ij m k l 2 m −1 n (α+l+2)n−j(β+k+2)j φ := − − − ψ . (4.7) ij ij i k 1 i j (α+l+2) (β+k+2) (cid:18) − − (cid:19)(cid:18) (cid:19) (cid:18) (cid:19) l+1 k+1 Observe that the quantities ψ can be put in a rectangular table and the entries of this ψ-table can be ij computed using [5, Algorithm 4.2], assuming that c := k+1, c := l+1, α := α and α := β. Note that 1 2 1 2 the complexity of this algorithm is O(mn). 4.2. Gk,l-constrained multi-degree reduction Coming back to the problem of Gk,l-constrained multi-degree reduction (see Problem 1.1), let us notice that the formulas (3.1)–(3.6) with fixed parameters λ and µ (cf. 3) constitute constraints of the form i j § demanded in Problem 4.1. As a result, the control points (4.3) depend on these parameters. Now, the optimum values of the parameters can be obtained by minimizing the error function (4.6), E(α,β) E(α,β)(λ ,λ ,...,λ ,µ ,µ ,...,µ ), (4.8) 1 2 k 1 2 l ≡ depending on λ and µ via formulas (3.1)–(3.6) and (4.3). i j { } { } For a minimum of function (4.8), it is necessarythat the derivativesof E(α,β) with respect to the param- eters are zero, which yields the system d m−l−1 ∂r F (rh) F (ph) jh =0 (u=1,2,...,k), mj nj − ∂λ u h=1 j=u X X (cid:2) (cid:3) (4.9) d m−v ∂r F (rh) F (ph) jh =0 (v =1,2,...,l), mj nj − ∂µ v h=1j=k+1 X X (cid:2) (cid:3) where we use notation t 1 m t Ftj(q):= (α+1)t+m−i−j(β+1)i+jqi (α+β+m+2) j i t(cid:18) (cid:19)i=0(cid:18) (cid:19) X 5 with q=[q ,q ,...,q ]. 0 1 t Inthecaseofk=l =3,wecomputethepartialderivativesofhthcoordinatesofthe controlpoints(3.1)– (3.6). We obtain: n∆p (i=1), m 0h  2n∆p +2λ (n−1)2∆2p (i=2), ∂∂rλi1h = 3mmn∆p00hh+ 2λ11(m+−1m)21−2λ20h3((mn−−11))22∆2p0h+3λ21((mn−−22))33∆3p0h (i=3), h i  0 (i=0; m−3≤i≤m)(,4.10) n ∆p (i=2), (m−1)2 0h ∂r ih = 3 n ∆p +3λ (n−1)2∆2p (i=3), (4.11) ∂λ2  (m−1)2 0h 1(m−2)3 0h 0 (i=0,1; m 3 i m), − ≤ ≤ ∂rih  (m−n2)3∆p0h (i=3), = (4.12) ∂λ 3 ( 0 (i=0,1,2; m 3 i m), − ≤ ≤ −mn∆pn−1,h (i=m−1), ∂∂rµi1h = −−23mmnn∆∆ppnn−−11,,hh++22µµ11((mn−−−11m))212−∆2µ2p2n−32((m,nh−−11))22∆2pn−2,h−3µ21((mn−−22))33∆3pn−3,h ((ii==mm−−32)),, h i  0 (0≤i≤3; i=(m4.)1,3) (m−n1)2∆pn−1,h (i=m−2), ∂r ∂µi2h = 3(m−n1)2∆pn−1,h−3µ1((mn−−12))23∆2pn−2,h (i=m−3), (4.14) 0 (0 i 3; i=m 1,m), ≤ ≤ − ∂rih  −(m−n2)3∆pn−1,h (i=m−3), = (4.15) ∂µ 3 ( 0 (0 i 3; m 2 i m). ≤ ≤ − ≤ ≤ Notice that the partial derivatives of hth coordinates of control points (4.3) depend on (4.10)–(4.15) in the following way: n −1 k ∂r n n m m ∂r ih gh = − φ , (4.16) ij ∂λ − j j g g ∂λ u j=0(cid:18) (cid:19) g=u(cid:18) − (cid:19)(cid:18) (cid:19) u X X n −1 m−v ∂r n n m m ∂r ih gh = − φ . (4.17) ij ∂µ − j j g g ∂µ v j=0(cid:18) (cid:19) g=m−l(cid:18) − (cid:19)(cid:18) (cid:19) v X X One can easily see, that when k,l 3, we compute ∂rih, ∂rih by (4.16), (4.17) if k < i < m l, and ≤ ∂λu ∂µv − by (4.10)–(4.15) otherwise. Finally, we put the expressions (4.10)–(4.17) into the equations of system (4.9). Observe that for k 2 or l 2, system (4.9) is nonlinear, which makes it quite difficult to solve. ≥ ≥ Furthermore, from a practical point of view, we additionally require that λ ,µ > 0, which results in the 1 1 same directions of tangent vectors at the endpoints of curves (1.1) and (1.3). Therefore, to guarantee that these conditions will be satisfied, it is not enough just to solve the system (4.9). Now, let us discuss two possible ways of determining the values of geometric continuity parameters. 4.2.1. Determining the Gk,l parameters using optimization methods It is easy to check that if (k = 1 and l 1) or (l = 1 and k 1), then the error (4.6) is a quadratic ≤ ≤ function of continuity parameters. 6 In the case of (k = 2 and l 2) or (l = 2 and k 2), the error (4.6) is a fourth-degree polynomial ≤ ≤ function of continuity parameters. For (k = 3 and l 3) or (l = 3 and k 3), the error (4.6) is a sixth-degree polynomial function of ≤ ≤ continuity parameters. To find the optimum values of parameters λ , µ in the case of G1,1-constrained multi-degree reduction 1 1 problem, assuming that α,β = 0, Lu and Wang [10] solve the quadratic programming problem, subject to the constraints λ d , µ d , (4.18) 1 0 1 1 ≥ ≥ where d and d are positive lower bounds, usually prescribed to small values (they set 10−4 for both lower 0 1 bounds in the examples section). Such approachcan be used in the cases which result in the quadratic error function (4.6). One can solve the quadratic programmingproblem using, e.g., an iterative active-set method, which is implemented in many software libraries. The active-set mechanism used by standard quadratic solvers is described in [1, 6.5]. § Analogously, one can observe that for k = 2,3 or l = 2,3, the problem of minimizing the error (4.6), subject to the constraints (4.18) is a nonlinear programming problem. To solve it, one can use, for instance, a sequential quadratic programming (SQP) method (see, e.g., [1, 15.1]), which is available in many software § libraries. 4.2.2. Determining the Cp,q/Gk,l parameters by solving a system of linear equations In the case of G2,2, Rababahand Mann [11] simplified the problem by consideringC1,1-continuity at the endpoints, i.e., they set λ = µ := 1. Later, this approach was also used by Lu [6]. In [12], the same idea 1 1 was used to simplify the G3,3 case, and the authors noted that such approach leads to a system of linear equations. Now, we generalize the above-described approach for any k, l such that 1 k,l 3. If k 2, we − ≤ ≤ ≥ set λ := 1, which implies C1-continuity at t = 0 and consequently, Gk,l constraints become C1,q/Gk,l 1 constraints, where q ,1 . Similarly, when l 2, we set µ := 1, which implies C1-continuity at t = 1 1 ∈ {− } ≥ and consequently, Gk,l constraints become Cp,1/Gk,l constraints, where p ,1 . ∈{− } Noticethatinthecasesofk =2,3orl=2,3,theabove-describedmethodleadstothelinearsystem(4.9) and the error (4.6) is a quadratic function of the continuity parameters. However, in the cases of k = 1 or l =1, there is no guarantee that the solution satisfies λ >0 or µ >0, respectively. In the case of solution 1 1 withnonpositivevaluesoftheseparameters,ourchoiceistosolveaquadraticprogrammingproblem,subject to the constraints with prescribed positive lower bounds for the parameters (see (4.18)). Observe that this approach uses no simplifying assumptions for k,l 1. ≤ Letusdenotetheabove-describedapproachtoProblem1.1asCp,q/Gk,l-constrainedmulti-degreereduction of B´ezier curves. 5. Algorithms In this section, we show the details of implementation of the proposedmethod of Gk,l-constrainedmulti- degreereductionofB´eziercurves. Moreover,wegiveashortdescriptionofCp,q/Gk,l-constrainedmulti-degree reduction algorithm. 5.1. Gk,l-constrained multi-degree reduction algorithm Now, we give the method of solving Problem 1.1, summarized in the following two-phase algorithm. PhaseAofthealgorithmconsistsinfindingvaluesoftheparametersλ andµ tominimizetheerror(4.6), i j which—bytheresultsgiveninTheorem4.2—dependsonlyontheseparameters. Theideaisbasedonsolving thequadraticornonlinearprogrammingproblem(see 4.2.1). Noticethatwhenk,l<1,wecancomputethe § φ coefficients (cf. (4.7)) and omit the remaining steps of Phase A, since there are no continuity parameters ij to determine. During Phase B, we use the results of Theorem 4.2 and the obtained values of continuity parametersto compute controlpoints r ,r ,...,r . Most of the knownalgorithms solvea systemof normal 0 1 m equations, to get the inner control points of multi-degree reduced curve (1.3). Such approach makes these points dependent on the inverse of a certain matrix. Our formulas do not require matrix inversion. What is more, the complexity of Phase B is O(mn), which is significantly less than the cost of other known methods for this phase. The algorithm works for any k and l not exceeding 3. Algorithm 5.1. [Gk,l-constrained multi-degree reduction] Data: α, β – parameters of the inner product (2.1); 7 n, p ,p ,...,p – degree and the control points of the B´ezier curve (1.1); 0 1 n m – degree of the reduced B´ezier curve (1.3); k, l – orders of the G-continuity at the endpoints of the curve (1.3); d , d – lower bounds for the parameters λ and µ , respectively (cf. 4.2.1). 0 1 1 1 § Assumptions: n>m>0; 1 k,l 3; k+l <m 1; d ,d >0; α,β > 1. 0 1 − ≤ ≤ − − Result: control points r ,r ,...,r of the Gk,l-constrained multi-degree reduced B´ezier curve (1.3). 0 1 m Phase A Step I Compute φ (i=k+1,k+2,...,m l 1; j =0,1,...,n)by[5, Algorithm 4.2] andformula (4.7) ij { } − − (see Remark 4.3). Step II Check if the remaining steps of Phase A can be omitted: If (k <1) and (l <1) then go to Step VI. Step III Compute E(α,β)(λ ,λ ,...,λ ,µ ,µ ,...,µ ) by (4.6). 1 2 k 1 2 l Step IV Determine set c of constraints: c:= λ d , µ d ; 1 0 1 1 { ≥ ≥ } If (k <1) then c:=c λ d ; 1 0 \{ ≥ } If (l <1) then c:=c µ d . 1 1 \{ ≥ } Step V If (k >1 or l >1) then obtain λ ,λ ,...,λ , and µ ,µ ,...,µ by solving the nonlinear programming problem of minimiz- 1 2 k 1 2 l ing the error (4.6), subject to the constraints c; else obtain λ ,λ ,...,λ , and µ ,µ ,...,µ by solving the quadratic programming problem of minimiz- 1 2 k 1 2 l ing the error (4.6), subject to the constraints c. Phase B Step VI Compute 1. r ,r ,...,r by (3.1)–(3.3); 0 1 k 2. rm−l,rm−l+1,...,rm by (3.4)–(3.6). Step VII Compute υ ,υ ,...,υ by (4.5). 0 1 n Step VIII Compute rk+1,rk+2,...,rm−l−1 by (4.3). Step IX Return the solution, i.e., the control points r ,r ,...,r of the reduced B´ezier curve (1.3). 0 1 m 5.2. Cp,q/Gk,l-constrained multi-degree reduction algorithm Now, let us give the outline of the two-phase Cp,q/Gk,l-constrained multi-degree reduction algorithm. PhaseAofthealgorithmimplementstheideasdiscussedin 4.2.2,therefore,itsolvesthesystemoflinear § equations (4.9) to determine values of the continuity parameters. In the case of solution with nonpositive values of λ or µ , which can happen when k = 1 or l = 1, the algorithm solves a quadratic programming 1 1 problem, subject to the constraints with prescribed positive lower bounds for the parameters (see (4.18)). An example of a resulting B´ezier curve that does not satisfy the positive condition can be found in [6, Fig. 1(a)]. We performed more than 40 different tests (results of some of them are available in the next section). None of them caused such problem. Phase B is the same as for Algorithm 5.1. The algorithm works for any k and l not exceeding 3. For details, see our implementation in MapleTM13 available on the website http://www.ii.uni.wroc.pl/~pgo/GDegRed.mws. Obviously, Algorithm 5.1 costs more, but also produces more accurate results, since for the Cp,q/Gk,l- constrained approach we additionally assume that λ =1 when k >1, and µ =1 when l >1. 1 1 8 6. Examples This section provides of the application of our Gk,l-constrained and Cp,q/Gk,l-constrained multi-degree reductionalgorithms. Ineachcase,wegivetheleastsquareserrorE(α,β) :=√E(α,β) andthemaximumerror 2 E∞ := max Pn(t) Rm(t) max Pn(t) Rm(t) , t∈DN|| − ||≈t∈[0,1]|| − || where D := 0,1/N,2/N,...,1 with N :=500. N { } In our experiments, we consider the “natural” choices for the values of parameters α, β, i.e., (α,β) (0,0), 1,1 , 1,1 , 1, 1 , 1, 1 , and set the lower bounds d ,d of λ ,µ to 10−4 (see (4.18)).∈ 2 2 −2 2 2 −2 −2 −2 0 1 1 1 Taking into account the different types of continuity constraints, we compare the following cases: (cid:8) (cid:0) (cid:1) (cid:0) (cid:1) (cid:0) (cid:1) (cid:0) (cid:1)(cid:9) (i) Ck,l-constrained case (see (1.5)), which can be solved by using Theorem 4.2; (ii) Cp,q/Gk,l-constrained case, solved by the algorithm described in 5.2; § (iii) Gk,l-constrained case, solved by Algorithm 5.1. Results of the experiments have been obtained on a computer with Intel Core i5-3337U 1.8GHz pro- cessor and 8GB of RAM, using 32-digit arithmetic. MapleTM13 worksheet containing implementation of the algorithms and tests is available on the website http://www.ii.uni.wroc.pl/~pgo/GDegRed.mws. We use MapleTM fsolve procedure, in the Cp,q/Gk,l case, to solve the system of linear equations, and QPSolve, NLPSolve procedures, to solve the quadratic and nonlinear programming problems, respectively. QPSolve uses the iterative active-set method, and for NLPSolve we select sqp method. Initial points for both proce- dures correspond to the values of continuity parameters in the Ck,l case. Example 6.1. First, let us consider degree eleven B´ezier curve which is an outline of the font “alpha” (for thecontrolpoints,see[15,Example6.1]). Theresultsofmulti-degreereductionaregiveninTable1. Figs.1a and 1b illustrate two of the considered cases. One can see, that when it comes to minimizing E∞ error, usually a good choice is α = β = 1. As expected, solution to the Gk,l case is the most accurate, while −2 Cp,q/Gk,l approach gives less precise results. Ck,l conditions seem to be too restrictive, especially for k or l exceeding 2. (a) (b) Figure1: Multi-degreereductionofdegreeelevencurve(bluesolidline)todegreesevenwithCk,l(blackdottedline),Cp,q/Gk,l(green dash-dotted line) and Gk,l (red dashed line) continuity constraints; parameters: (a) α = β = −21, p = q = 1, k = l = 2, and (b) α=β=−21,p=1,q=−,k=3andl=1. 9 Parameters Ck,l solution Cp,q/Gk,l solution Gk,l solution m k l p q α β E2(α,β) E∞ E2(α,β) E∞ E2(α,β) E∞ 7 2 2 1 1 0 0 3.73 5.97 2.83 5.27 0.93 2.26 1 1 5.75 5.83 4.40 5.14 1.67 1.91 −2 −2 1 1 3.83 7.53 2.62 6.42 0.79 2.91 −2 2 1 1 3.83 7.69 3.18 5.18 1.26 2.19 2 −2 1 1 2.43 6.10 1.83 5.40 0.53 2.54 2 2 7 3 1 1 0 0 9.13 16.24 2.51 5.11 1.02 2.41 − 1 1 13.97 16.72 4.07 4.95 1.81 2.07 −2 −2 1 1 9.41 19.51 2.18 6.38 0.75 3.11 −2 2 1 1 9.17 18.79 2.98 4.91 1.45 1.99 2 −2 1 1 6.00 15.82 1.56 5.25 0.59 2.69 2 2 Table1: Leastsquareserrorandmaximumerrorinmulti-degreereductionofdegreeelevenB´ezier“alpha”curve. Example 6.2. Let us apply the algorithms to degree thirteen B´ezier “heart” curve (for the control points, see [12, Appendix B]) and consider the case of k = l = 2. The results of experiments are given in Table 2. Notice that the case of α = β = 0 was also considered in [12, 5.2] and [17, Example 4]. As in [17], we can § clearlysee that the solutionto the G2,2 case, in this paper obtainedby Algorithm5.1, is more accuratethan the result given by the approach proposed in [12], which leads to the C1,1/G2,2 case (the same as for the algorithm discussed in 5.2). As our approach considers different weight functions, it can be seen that the § best choice to minimize E∞ error is α=β = 1. Fig. 2 presents α=β =0 case. −2 Now, we focus on the running times. We have implemented G1,1, C1,1/G2,2 and C1,1/G3,3-constrained methods from [12], G1,1 and C1,1/G2,2-constrained methods from [6] as well as G2,1 and G2,2-constrained methods from [17]. The methods of Rababah and Mann and of Lu solve the same problem and give the same results as our Cp,q/Gk,l-constrained method (see 5.2). In Table 3, we compare the running times § of these algorithms. Clearly, our approach is the fastest one. For the comparison of the Gk,l-constrained algorithms, see Table 4. Notice that, in some cases, our Gk,l-constrained approach is slightly faster than the methods from [17]. We use MapleTM fsolve procedure to solve the cubic equation [17, (23)] associated with the G2,1-constrained case. The implementation of G2,2-constrained method from [17] requires the unconstrained nonlinear programming solver. According to our experiments, the nonlinear simplex method (NLPSolvecommandwithoptionmethod = nonlinearsimplexandtheinitialpointλ=η =1)isthefastest solver available in MapleTM13. Therefore, we use this solver for the purpose of the comparison. It is worth mentioning that Zhou et al. have omitted the constraints (4.18). Consequently, in some rare cases, the resulting curve may not preserve the original tangent directions at the endpoints. To avoid this issue, one can implement the improvements proposed by Lu [8]. Parameters C2,2 solution C1,1/G2,2 solution G2,2 solution m α β E2(α,β) E∞ E2(α,β) E∞ E2(α,β) E∞ 8 0 0 1.52 2.52 0.64 1.12 0.36 0.71 1 1 2.37 2.39 1.05 1.00 0.62 0.53 −2 −2 1 1 1.58 3.34 0.64 1.55 0.42 0.94 −2 2 1 1 1.52 3.48 0.74 1.31 0.37 0.89 2 −2 1 1 0.98 2.64 0.39 1.24 0.21 0.90 2 2 Table2: Leastsquareserrorandmaximumerrorinmulti-degreereductionofdegreethirteenB´ezier“heart”curve. 10