- Email: [email protected]
An iterative algorithm for G2 multiwise merging of Bézier curves
Journal of Computational and Applied Mathematics 296 (2016) 352–361
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
An iterative algorithm for G2 multiwise merging of Bézier curves Lizheng Lu ∗ , Chengkai Jiang School of Statistics and Mathematics, Zhejiang Gongshang University, Hangzhou 310018, China
article
abstract
info
This paper presents an iterative algorithm for G2 multiwise merging of Bézier curves. By using the G2 constraint, the L2 distance is represented after simplification as a quartic polynomial in two parameters relating to the magnitudes of end tangents of the merged curve. These two parameters are restricted in a feasible region, in order for the merged curve to preserve the specified directions of end tangents. Then G2 multiwise merging is formulated as a constrained minimization problem, and the classic projected Newton method is applied to find the minimizer. Some extensions of multiwise merging using G3 constraints, other energy functionals and curve representations are also outlined. Several comparative examples are provided to demonstrate the effectiveness of the proposed method. © 2015 Elsevier B.V. All rights reserved.
Article history: Received 4 May 2015 Received in revised form 4 October 2015 Keywords: Bézier curve Merging Geometric continuity L2 distance Iterative method
1. Introduction Assume that the given are K (≥ 2) segments of adjacent Bézier curves of degrees nk , k = 1, . . . , K , each of which is defined in terms of the control points Pknk := (pk0 , . . . , pknk )T as Pk (t ) =
nk
n
Bi k (t )pki = Bnk Pknk ,
t ∈ [0, 1].
(1)
i=0
Throughout, the notation Bn := Bn0 (t ), . . . , Bnn (t ) is used to denote the row vector of the Bernstein basis of degree n, with
Bni (t ) = i t i (1 − t )n−i . The problem of multiwise merging [1,2] aims to find a single degree-n Bézier curve (called merged curve) with control points Rn := (r0 , . . . , rn )T ,
n
R(t ) =
n
Bni (t )ri = Bn Rn ,
t ∈ [0, 1],
(2)
i=0
such that the L2 distance E=
K k=1
∗
tk+1 tk
2 K R(t ) − Pk t − tk dt = ∆ t k ∆t k
Corresponding author. E-mail address: [email protected] (L. Lu).
http://dx.doi.org/10.1016/j.cam.2015.10.007 0377-0427/© 2015 Elsevier B.V. All rights reserved.
k=1
1 0
∥Bn Sk Rn − Bnk Pknk ∥2 dt
(3)
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
353
is minimized. Here, 0 = t1 < t2 < · · · < tK +1 = 1 is a partition of the interval [0, 1], ∆tk = tk+1 − tk , and Sk Rn represents the control points of the sub-curve R(t ), t ∈ [tk , tk+1 ] after it is reparametrized to the standard interval [0, 1]. See [1] for the expression of the (n + 1) × (n + 1) matrix Sk . In geometric modeling, a complex shape is generally constructed by using multiple segments, with the advantage of producing more desirable shapes through local adjustment. It is often required to merge multiple adjacent segments into a single one, while preserving as faithful as possible to the original shape. For this goal, the merging problem is studied as a minimization problem using some suitable distances such as (3). To meet certain smoothness requirement in various applications, it is important to impose continuity conditions on the two endpoints of the merged curve R(t ), which correspond to the initial endpoint of the first segment P1 (t ) and the terminal endpoint of the last segment PK (t ). As multiple segments are merged at a time, multiwise merging does not need to consider continuity conditions at the interior endpoints. It is well known [3,4] that in contrast to the classical parametric continuity (C k ), geometric continuity (Gk ) shares additional parameters that are available for further shape manipulation by direct assignment or optimization techniques. Noting the fact that geometric continuity is a relaxation of parametrization (but not a relaxation of smoothness), it is more suitable to use geometric continuity for the merging problem since the parametrization of unknown merged curve should not be simply taken as that of the first (or last) segment. For instance, G2 continuity (mostly used in applications) is able to retain the positions, tangent directions and curvatures. 1.1. Related work Most of the existing merging approaches can be categorized into matrix-based approaches. The merging problem was firstly introduced in 1987 by Hoschek [5], where many segments of curves are approximately converted into a spline curve through degree reduction and splitting. Although degree reduction was widely investigated by many researchers, merging has not received adequate attention until Hu et al. [6] proposed a pioneering work of pairwise merging where a pair of Bézier curves is merged into one segment instead of applying degree reduction twice. Then, merging has become a special research topic. In terms of matrix representations and operations, many approaches have arisen, including pairwise merging of B-spline curves [7] and Bézier curves [8,9], and multiwise merging of Bézier curves [10,1] and B-spline surfaces [11]. In particular, C r ,s multiwise merging in [1] can preserve C r continuity and C s continuity at the two endpoints respectively. A novel dual approach was proposed by Woźny et al. [2] for C r ,s multiwise merging, which makes use of the advantages of constrained dual Bernstein basis [12,13]. Remarkably, this approach has lower complexity due to the fast evaluation schemes of certain connections between Bernstein and dual Bernstein polynomials. 1.2. Contributions In this paper, we revisit the multiwise merging problem with the constraint of G2 or G3 continuity, and propose an iterative algorithm based on the projected Newton method. According to the G2 constraint, there are four parameters (i.e., α1 , α2 , β1 , β2 in (5)–(8)), and two more parameters by the G3 constraint. Similar to the degree reduction work of Lu [14], we also find that for G2 multiwise merging, α2 and β2 are both expressed as quadratic polynomials of α1 and β1 when the minimum of the L2 distance (3) is reached; therefore, the number of independent parameters is reduced to 2, and the L2 distance is a quartic polynomial in two unknowns. It is important to note that many geometric modeling applications often require the two parameters α1 and β1 to be positive, so that the end tangents of the merged curve coincide with the specified directions, referring to (5) and (7). On the other hand, the values of α1 and β1 indicate the magnitudes of the end tangents. For these reasons, we impose a feasible region on (α1 , β1 ):
D = (α1 , β1 ) ∈ R2 : 0 < l0 ≤ α1 ≤ u0 , 0 < l1 ≤ β1 ≤ u1 .
(4)
By restricting (α1 , β1 ) in the feasible region, the merging problem is changed to a constrained minimization problem which may be efficiently solved by the projected Newton method, as done in [14] for degree reduction. Our contributions are summarized as follows:
• An iterative algorithm is proposed for G2 multiwise merging with the feasible region. • An intuitive interaction tool is provided for shape-preserving merging, in terms of the feasible region. Compared to the previous C r ,s multiwise merging [1,2], our method produces better results due to the used G2 and G3 constraints. Our method can be regarded as a generalization of the previous pairwise merging [8] which uses C 1 G2 and C 1 G3 constraints. Very recently, Gospodarczyk and Woźny [15] proposed a new method for C r ,s multiwise merging by restricting all the interior control points of the merged curve in a specified region; whereas our method restricts (α1 , β1 ) in a specified region. 2. Iterative algorithm for G 2 multiwise merging Given K segments of adjacent Bézier curves Pk (t ) with control points Pknk , k = 1, . . . , K , G2 multiwise merging amounts
to identifying the control points Rn of the merged curve R(t ), such that the L2 distance is minimized and the G2 constraint at the endpoints is satisfied.
354
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
By the G2 constraint, the first three and last three control points in Rn are determined by [3,4]: r0 = p10 ,
n1
r1 = p10 +
r2 = p10 +
2n1 n
rn = pKnK ,
α1 ∆p10 +
α1 ∆p10 , n n1 (n1 − 1) n(n − 1) nK
rn−1 = pKnK − 2nK
rn−2 = pKnK −
n
(5)
α12 ∆2 p10 +
n1 n(n − 1)
α2 ∆p10 ,
(6)
β1 ∆pKnK −1 ,
n
β1 ∆pKnK −1 +
(7)
nK (nK − 1) n( n − 1 )
nK
β12 ∆2 pKnK −2 +
n(n − 1)
β2 ∆pKnK −1 ,
(8)
where α1 , α2 , β1 , β2 ∈ R, and the forward difference operator ∆ is defined by ∆pki = pki+1 − pki and by ∆r pki = ∆r −1 pki+1 − ∆r −1 pki for r > 1. The merged curve is rewritten as R(t ) = Bn0 (t )r0 + Bn1 (t )r1 + Bn2 (t )r2 + Bnn−2 (t )rn−2 + Bnn−1 (t )rn−1 + Bnn (t )rn +
n −3
Bni (t )ri .
i =3
Let Rcn := (r0 , r1 , r2 , rn−2 , rn−1 , rn )T and the L2 distance in (3) is expanded by E=
K
∆tk
:= (r3 , . . . , rn−3 )T be the constrained and free control points respectively. Then
1
0
k=1
f Rn
∥Bn Skc Rcn + Bn Skf Rfn − Bnk Pknk ∥2 dt
=(Rcn )T Hcc Rcn + 2(Rfn )T Hfc Rcn + (Rfn )T Hff Rfn − 2(Rcn )T
K
∆tk (Skc )T Gn,nk Pknk
k=1
− 2(Rfn )T
K
f
∆tk (Sk )T Gn,nk Pknk +
k=1
K
∆tk (Pknk )T Gnk ,nk Pknk
(9)
k=1
n ,m
where Gn,m = (gij )i=0,j=0 ∈ R(n+1)×(m+1) with entries 1
Bni (t )Bm j (t )dt =
gij = 0
1 n+m+1
n
m
n+m
i
j
i+j
and Hcc =
K
∆tk (Skc )T Gn,n Skc ,
k=1
Hfc =
K
f
∆tk (Sk )T Gn,n Skc ,
Hff =
k=1
K
f
f
∆tk (Sk )T Gn,n Sk ,
k=1 f
with Skc := Sk (· · · ; 0, 1, 2, n − 2, n − 1, n) and Sk := Sk (· · · ; 3, . . . , n − 3) denoting the submatrices of Sk formed by selecting the specified columns. f
At first glance, it is natural to determine α1 , β1 , α2 , β2 and Rn using (9) via min
f α1 ,β1 ,α2 ,β2 ,Rn
{E := E (α1 , β1 , α2 , β2 , Rfn )}.
(10)
However, as will be shown in the following subsection, a further analysis reveals that we are able to decrease the complexity of the problem.
2.1. Simplification The expression (9) can be simplified through three steps: 1. Express ri = ri (α1 , β1 , α2 , β2 ), i = 3, . . . , m − 3 as quadratic polynomials, which are quadratic in α1 , β1 and linear in α2 , β2 . 2. Express α2 = α2 (α1 , β1 ) and β2 = β2 (α1 , β1 ) as quadratic polynomials of α1 , β1 . 3. Express E = E (α1 , β1 ) as a quartic polynomial of α1 , β1 .
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
355
f
Firstly, assuming that Rcn and Rn are variables, it follows from (9) that
K ∂E ∂E ∂E ∂E ∂E ∂E ∂E T c T f := , , , , , = 2H R + 2H R − 2 ∆tk (Skc )T Gn,nk Pknk , cc n fc n ∂ Rcn ∂ r0 ∂ r1 ∂ r2 ∂ rn−2 ∂ rn−1 ∂ rn k =1 K ∂E ∂E ∂E T f f c := , . . . , = 2H R + 2H R − 2 ∆tk (Sk )T Gn,nk Pknk . ff fc n n ∂ r3 ∂ rn−3 ∂ Rfn k=1
(11)
(12)
By (12), Rfn = −Hff−1 Hfc Rcn + Hff−1
K
f
∆tk (Sk )T Gn,nk Pknk
(13)
k=1
which implies that ri , i = 3, . . . , m − 3 are quadratic polynomials of α1 , β1 , α2 , β2 . Secondly, note the following facts: r1 is linear in α1 ; r2 is quadratic in α1 and linear in α2 ; rn−1 is linear in β1 ; rn−2 is quadratic in β1 and linear in β2 . After substituting (13) into (9), E is a quartic polynomial of α1 , β1 , α2 , β2 and takes the form E (α1 , β1 , α2 , β2 ) =
1 2
α2
β2
a00
a01 a11
a10
α2 + α2 b0 (α1 , β1 ) + β2 b1 (α1 , β1 ) + c (α1 , β1 ) β2
(14)
where the constants aij and quadratic polynomials bi (α1 , β1 ) will be derived below, and c (α1 , β1 ) = E (α1 , β1 , 0, 0) by (9) and (13). From (5)–(8) and (13), it is easy to show that
∂ Rcn = 0 ∂α2 ∂ Rcn = 0 ∂β2
T
,
T
,
0
v0
0
0
0
0
0
v1
0
0
n
∂ Rfn = −Hff−1 Hfc 0 ∂α2 ∂ Rfn = −Hff−1 Hfc 0 ∂β2
T
,
T
,
0
v0
0
0
0
0
0
v1
0
0
n
K 1 ∆p10 and v1 = n(n− ∆pKnK −1 . Differentiating by the chain rule gives where v0 = n(n− 1) 1)
∂E = 2A · 0 ∂α2
0
v0
0
0
T
0
,
∂E = 2A · 0 ∂β2
0
0
v1
0
T
0
,
(15)
where A = (Hcc − HfcT Hff−1 Hfc )Rcn −
K
∆tk (Skc )T Gn,nk Pknk + HfcT Hff−1
K
f
∆tk (Sk )T Gn,nk Pknk .
k=1
k=1
Using (14) and (15), we obtain b0 (α1 , β1 ) =
∂ E (α1 , β1 , 0, 0) , ∂α2
b1 (α1 , β1 ) =
∂ E (α1 , β1 , 0, 0) ∂β2
(16)
and
∂ 2E ∂α 2 2 = ∂ 2E ∂β2 ∂α2
a00 a10
a01 a11
∂ 2E 2w22 v0 · v0 ∂α2 ∂β2 = 2w32 v1 · v0 ∂ 2E
2w23 v0 · v1 2w33 v1 · v1
(17)
∂β22
5,5 ∂E ∂E where wij are the entries of the 6 × 6 matrix W = wij i=0,j=0 := Hcc − HfcT Hff−1 Hfc . By applying ∂α = ∂β = 0 to (14), α2 2 2 and β2 are solved as α2 =
−a11 b0 (α1 , β1 ) + a01 b1 (α1 , β1 ) , a00 a11 − a201
β2 =
a10 b0 (α1 , β1 ) − a00 b1 (α1 , β1 ) a00 a11 − a201
,
(18)
both of which are linear combinations of quadratic polynomials b0 , b1 and are thus quadratic polynomials of α1 , β1 . Thirdly, after substituting (18) into (14), E is a quartic polynomial of α1 , β1 . Based on the above three steps, we finally arrive at the constrained minimization problem: min {Eˆ (α1 , β1 ) := E (α1 , β1 , α2 (α1 , β1 ), β2 (α1 , β1 ))}.
(α1 ,β1 )∈D
(19)
356
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
2.2. Numerical solution To solve the constrained minimization problem (19), we resort to some numerical methods. Since the objective function is quartic in two unknowns, it is trivial to compute its gradient and Hessian matrix. The projected Newton method [16] is employed, which is an iterative method and achieves a quadratic convergence rate. For numerical optimization, it is advisable to start from an initial point that is not far away from a local minimum. If (1, 1) ∈ D , (α1 , β1 ) = (1, 1) is initially assigned, which corresponds to the C 1 continuity conditions; otherwise, the initial value of (α1 , β1 ) can be taken as the central point of D . 3. Extensions G2 multiwise merging described detailedly in Section 2 can be extended in several ways, with similar derivations. Three directions are outlined in the following. 3.1. G3 constraint Other constraints of geometry continuity can be used, e.g., G1,2 constraint, G2,3 constraint, G3 constraint, etc. We do not suggest to use higher-order geometric continuity of order greater than 3, because G2 continuity and G3 continuity have sufficient smoothness properties for practical applications and more degrees of freedom should be left to the interior control points of merged curves for generating more satisfactory results. For example, G3,2 multiwise merging is briefly described below. In this case, R(t ) and P1 (t ) are G3 continuous at t = 0, while R(t ) and PK (t ) are G2 continuous at t = 1; hence, besides (5)–(8), the control point r3 is expressed using a new parameter α3 ∈ R as 3n1
n 3n1
α1 ∆p10 +
3n1 (n1 − 1)
α12 ∆2 p10 +
n1 (n1 − 1)(n1 − 2)
α13 ∆3 p10 n(n − 1)(n − 2) n1 α2 ∆p10 + α1 α2 ∆2 p10 + α3 ∆p10 . + n( n − 1 ) n(n − 1)(n − 2) n(n − 1)(n − 2)
r3 = p10 +
n(n − 1) 3n1 (n1 − 1)
(20)
In a similar way, the L2 distance is simplified through three steps: 1. Express ri = ri (α1 , α2 , β1 , α3 , β2 ), i = 4, . . . , m − 3 as cubic polynomials, which are cubic in α1 , quadratic in β1 and linear in α2 , α3 , β2 . 2. Express α3 = α3 (α1 , α2 , β1 ) and β2 = β2 (α1 , α2 , β1 ) as cubic polynomials of α1 , α2 , β1 . 3. Express E = E (α1 , α2 , β1 ) as a sextic polynomial of α1 , α2 , β1 . Then, a similar constrained minimization problem like (19) is built, but now the objective function is sextic in three unknowns. 3.2. Energy functional Obviously, the l2 distance in [10,6,8,1] can be used as the energy functional. If geometric information of higher-order derivatives of the given curves is considered, then the L2 distance in (3) is further modified as E˜ =E +
2
ωℓ
ℓ=1
K k=1
tk+1 tk
ℓ 2 d dℓ k t − tk dt R ( t ) − P dt ℓ dt ℓ ∆t k
2 1 K d d 1 k Bn Sk Rn − Bnk Pnk dt = ∆tk ∥Bn Sk Rn − ∥ + ω1 ∆tk 0 dt dt 0 k=1 k=1 K 2 1 2 2 1 d Bn Sk Rn − d Bn Pk dt + ω2 n k k (∆tk )3 0 dt 2 dt 2 k=1 K
1
Bnk Pknk 2 dt
(21)
where the weights ω1 , ω2 ≥ 0 control relative weighting of the three involved terms. Note that E˜ is identical to E when ω1 = ω2 = 0 and that the second derivatives are not included in E˜ when ω1 > 0 and ω2 = 0. We do not suggest to include ˜ since more weights would be introduced that unavoidably increase the difficulty of third (or higher-order) derivatives in E, their adjustments.
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
357
The whole process of Section 2 can be kept to minimize the new energy functional (21), but with the modification that n ,m n ,m for every matrix Gn,m = (gij )i=0,j=0 ∈ R(n+1)×(m+1) , two corresponding matrices Gℓn,m = (gijℓ )i=0,j=0 ∈ R(n+1)×(m+1) , ℓ = 1, 2 are additionally required, where gijℓ =
ℓ m dℓ Bni (t ) d Bj (t )
1
dt ℓ
0
dt ℓ
dt .
Although Gn,n (called the Gram matrix of the Bernstein basis) is symmetric positive definite [1], G1n,n and G2n,n are both n d ℓ n dℓ n dℓ n symmetric positive semidefinite since dt ℓ B0 (t ), . . . , dt ℓ Bn (t ) are linearly dependent by i=0 dt ℓ Bi (t ) ≡ 0, ∀t ∈ [0, 1]. 3.3. Curve representation Let ( u0 (t ), . . . , un (t )) defined on the interval [a, b] be a basis of a function space, and represent the curve in this space as n C(t ) = i=0 ui (t )pi , t ∈ [a, b]. Many kinds of bases and curve representations have been proposed for various applications in CAGD. As is well known, the Bernstein–Bézier representation has gained widespread use in CAGD and related fields, mainly because it exhibits elegant geometric properties [3,17]. Although we formulate multiwise merging in terms of Bézier curves, multiwise merging of curves in other representations can also be solved in the same manner. For the extended Chebyshev space
Tα2n = span{cos(kt ), sin(kt ) : t ∈ [0, α]}nk=0 ,
α ∈ (0, π )
of trigonometric polynomials of degree 2n (order n), a trigonometric curve of degree 2n is defined as the convex combination (see [18]) Cαn (t ) =
2n
α T2n ,i (t )pi ,
t ∈ [0, α],
(22)
i=0
where α
α
T2n,i (t ) = τ2n,i sin
2n−i
α−t 2
i
sin
i
t
2
,
α
τ2n,i
⌊2⌋ α i−2r n i−r := 2 cos 2 r sin2n α2 r =0 i − r 1
are the trigonometric normalized B-basis functions. When it comes to implementing multiwise merging of trigonometric curves, matrices Sk are specified by a de Casteljau-type algorithm [19,20], whereas all the entries α α the subdivision 2n,2m α gij = 0 T2n ,i (t )T2m,j (t )dt of G2n,2m = (gij )i=0,j=0 in this case must be evaluated via numerical integration and can be stored in lookup tables in advance. Finally, we would like to mention that multiwise merging can be performed on hyperbolic curves [20,21], algebraictrigonometric curves [22,23], algebraic-hyperbolic curves [23], etc. Note that when the shape parameter α in (22) approaches 0, the trigonometric curve becomes a Bézier curve of degree 2n, see [24, Proposition 2.1]; a similar asymptotic behavior occurs also in the hyperbolic curve case, see [21, Proposition 5]. 4. Examples and applications In this section, we present some examples for the proposed Gk multiwise merging and compare with the previous C k multiwise merging [1,2]. All examples are computed on a PC with Intel Core2 Quad CPU Q9400 2.66 GHz processor and 4 GB RAM. The merging results depend on the choice of the parameters {tk }Kk=+11 . Let t1 = 0, and assign λ,µ,ω
tk+1 = tk +
Vk
V λ,µ,ω
,
k = 1, . . . , K ,
(23)
where λ,µ,ω
Vk
1
:= 0
d k λ P (t ) dt + ω dt
1
|κk (t )|µ dt and V λ,µ,ω :=
0
K
λ,µ,ω
Vk
k=1
with
d k 3 d k d2 k κk (t ) = P (t ) × 2 P (t ) P (t ) dt dt dt denoting the curvature of the kth segment. If ω = 0, it is the commonly used arc length parametrization when λ = 1, the uniform parametrization when λ = 0, and the centripetal parametrization when λ = 1/2. However, if ω > 0, curvature
358
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
.
.
.
.
.
.
Fig. 1. Merging results for the treble clef curve using C k and Gk constraints.
Table 1 Comparative data between C k and Gk multiwise merging for the treble clef curve. (Note: the computation time is measured in milliseconds.) n
C2
C3
G2
G3
E
Time
E
Time
E
Iters
Time
E
Iters
Time
12 14 16 18 20
1.225e−2 3.903e−3 2.145e−3 1.063e−3 5.480e−4
18 21 24 28 32
4.966e−2 2.100e−2 8.801e−3 5.199e−3 2.788e−3
16 18 21 25 29
2.077e−3 6.545e−4 2.434e−4 1.667e−4 8.289e−5
5 6 5 5 6
30 32 35 40 43
2.379e−3 6.708e−4 2.573e−4 2.039e−4 8.940e−5
20 18 21 19 21
73 76 82 83 85
variations of the given Bézier segments are considered as well; this choice may provide aesthetically more pleasing results in case of high curvature segments. The first example is the treble clef curve composed of four Bézier curves of degree 5; see [1] for detailed control points. Fig. 1 illustrates the merged curves of degrees n, in which {t1 , t2 , . . . , t5 } = {0, 0.2387, 0.4872, 0.7626, 1} is taken by arc length parametrization. When applying G2 and G3 multiwise merging, the feasible region is set to D = [0.2, 5]2 . Clearly, owing to the used Gk constraint, better results are generated by Gk multiwise merging. As shown in Table 1, Gk multiwise merging requires more computation time than C k multiwise merging, caused by the three simplification steps and iterative process for solving the constrained minimization problem; and, compared to G2 multiwise merging, G3 multiwise merging increases the number of iterations and thus the computation time, due to the iterative process for solving a more complex problem. The second example is the penguin curve, whose left and right parts are composed of four and three cubic Bézier curves respectively; see [2] for detailed control points. Fig. 2 illustrates the merged curves of degrees n, and D = [0.2, 5]2 is used for the feasible region. Apparently, Gk multiwise merging greatly improves the results. In Fig. 2(a)–(d), {t1 , t2 , . . . , t5 } = {0, 0.08, 0.55, 0.78, 1} for the left part and {t1 , t2 , t3 , t4 } = {0, 0.42, 0.78, 1} for the right part are taken by arc length parametrization. In Fig. 2(e)–(h), {t1 , t2 , . . . , t5 } = {0, 0.18, 0.53, 0.78, 1} for the left part and {t1 , t2 , t3 , t4 } = {0, 0.40, 0.75, 1} for the right part are taken by setting (λ, µ, ω) = (1, 0.5, 0.05) in (23): the results by C k multiwise merging are noticeably improved, referring to the segment at the left bottom; whereas the improvement of Gk multiwise merging is not obvious, possibly owing to the used Gk constraint. We show an application to improve the parametrization of merged curves by using the feasible region. A curve is thought of as having optimal parametrization if it is parametrized by arc length [25], but which is impossible for polynomial curves. While for a merged curve R(t ), the magnitudes of all its tangent vectors, i.e. ∥R′ (t )∥, should be as close to equal as possible in order to achieve an optimal parametrization. Here, we do not intend to find a merged curve with the closest optimal parametrization; instead, we improve the parametrization by interactively adjusting the four bounds in the feasible region. Recalling that these four bounds indirectly specify the ranges of ∥R′ (0)∥ and ∥R′ (1)∥, such interaction is meaningful and flexible. An example is shown in Fig. 3. By making use of the two weights in (21), it increases the flexibility of shape design. As the first term of (21) measures the deviation in positions, small weights may be used for generating satisfactory results from the viewpoint of curve profiles; on the other hand, large weights should be taken if higher-order differential geometric quantities are regarded as important. However, a qualitative criterion for determining suitable values of weights still seems to be unattainable; alternatively, we need to manually assign and adjust their values. As shown in Fig. 4, by setting ω1 = ω2 = 1, the curvature distributions of
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
(a) n = 14.
(b) n = 14.
(c) n = 16.
(d) n = 16.
(e) n = 14.
(f) n = 14.
(g) n = 16.
(h) n = 16.
359
Fig. 2. Merging results for the penguin curve using C k and Gk constraints.
merged curves are improved. We noticed that the improvement is not obvious when the given curves are not C 2 connected, such as the treble clef curve and penguin curve. Finally, Fig. 5 demonstrates an example where multiwise merging is performed on three Bézier curves and three trigonometric curves with α = π /2. Since different curve representations are used, there exist slight differences between the merging results. 5. Conclusion We have developed an iterative algorithm for Gk multiwise merging of Bézier curves. The key idea of the proposed method is to make use of the additional parameters afforded by the Gk constraint, with the ability to produce better merging results. Moreover, the added flexibility can be exploited in many ways, and constitutes an important advance in geometric modeling applications, e.g., improving the parametrization of curves. We also propose to use the feasible region for the purpose of shape preservation and interaction. Acknowledgments This work is supported by the National Natural Science Foundation of China (Nos. 61272307, 61502128). The authors thank the anonymous referees for their constructive comments and for pointing out two extension directions (Sections 3.2 and 3.3).
360
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
(a) n = 9, D = [4, 5] × [2, 2.5].
(b) n = 10, D = [3.5, 4.5] × [3.5, 4.5]. 2
Fig. 3. Applying G multiwise merging to improve the parametrization of merged curves by using the feasible region D . The middle column displays the profiles of ∥R′ (t )∥ with respect to t, while the means and standard deviations of {∥R′ (t¯i )∥}Ni=1 evenly sampled in t are listed in the right column.
(a) n = 10, D = [0.2, 5]2 .
(b) n = 12, D = [0.2, 5]2 .
Fig. 4. Influence of the weights ω1 and ω2 : Red, ω1 = ω2 = 0; Blue, ω1 = ω2 = 1. In (a) and (b), the merged curves with G2 continuity are displayed on the left, while the profiles of curvatures with respect to arc length on the right. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(a) Bézier curves.
(b) Trigonometric curves.
Fig. 5. Merging results for three Bézier curves and three trigonometric curves with α = π/2, using D = [0.2, 5]2 . Note that the control points of the given curves in (a) and (b) are the same.
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
361
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25]
L. Lu, Explicit algorithms for multiwise merging of Bézier curves, J. Comput. Appl. Math. 278 (2015) 138–148. P. Woźny, P. Gospodarczyk, S. Lewanowicz, Efficient merging of multiple segments of Bézier curves, Appl. Math. Comput. 268 (2015) 354–363. G. Farin, Curves and Surfaces for CAGD: A Practical Guide, fifth ed., Morgan Kaufmann, San Francisco, 2002. J. Peters, Geometric continuity, in: G. Farin, J. Hoschek, M.-S. Kim (Eds.), Handbook of Computer Aided Geometric Design, Elsevier, Amsterdam, 2002, pp. 193–227. J. Hoschek, Approximate conversion of spline curves, Comput. Aided Geom. Design 4 (1987) 59–66. S. Hu, R. Tong, T. Ju, J. Sun, Approximate merging of a pair of Bézier curves, Comput. Aided Des. 33 (2001) 125–136. C. Tai, S. Hu, Q. Huang, Approximate merging of B-spline curves via knot adjustment and constrained optimization, Comput. Aided Des. 35 (2003) 893–899. L. Lu, An explicit method for G3 merging of two Bézier curves, J. Comput. Appl. Math. 260 (2014) 421–433. P. Zhu, G. Wang, Optimal approximate merging of a pair of Bézier curves with G2 -continuity, J. Zhejiang Univ. Sci. A 10 (2009) 554–561. M. Cheng, G. Wang, Approximate merging of multiple Bézier segments, Progr. Nat. Sci. 18 (2008) 757–762. H. Pungotra, G.K. Knopf, R. Canas, Merging multiple B-spline surface patches in a virtual reality environment, Comput. Aided Des. 42 (2010) 847–859. S. Lewanowicz, P. Woźny, Bézier representation of the constrained dual Bernstein polynomials, Appl. Math. Comput. 218 (2011) 4580–4586. P. Woźny, S. Lewanowicz, Multi-degree reduction of Bézier curves with constraints, using dual Bernstein basis polynomials, Comput. Aided Geom. Design 26 (2009) 566–579. L. Lu, Some improvements on optimal multi-degree reduction of Bézier curves with geometric constraints, Comput. Aided Des. 59 (2015) 39–42. P. Gospodarczyk, P. Woźny, Merging of Bézier curves with box constraints, J. Comput. Appl. Math. (2016) http://dx.doi.org/10.1016/j.cam.2015.10.005. to appear. M. Schmidt, D. Kim, S. Sra, Projected Newton-type methods in machine learning, in: S. Sra, S. Nowozin, S.J. Wright (Eds.), Optimization for Machine Learning, MIT Press, Cambridge, MA, 2011, pp. 305–330. J.M. Peña, Shape Preserving Representations in Computer-Aided Geometric Design, Nova Science Publishers, Commack, New York, 1999. J. Sánchez-Reyes, Harmonic rational Bézier curves, p-Bézier curves and trigonometric polynomials, Comput. Aided Geom. Design 15 (1998) 909–923. E. Mainar, J.M. Peña, Corner cutting algorithms associated with optimal shape preserving representations, Comput. Aided Geom. Design 16 (1999) 883–906. Á Róth, Control point based exact description of trigonometric/hyperbolic curves, surfaces and volumes, J. Comput. Appl. Math. 290 (2015) 74–91. W. Shen, G. Wang, A class of quasi Bézier curves based on hyperbolic polynomials, J. Zhejiang Univ. Sci. 6A (Suppl. I) (2005) 116–123. J.M. Carnicer, E. Mainar, J.M. Peña, Critical length for design purposes and extended Chebyshev spaces, Constr. Approx. 20 (2004) 55–71. E. Mainar, J.M. Peña, Optimal bases for a class of mixed spaces and their associated spline spaces, Comput. Math. Appl. 59 (2010) 1509–1523. I. Juhász, Á Róth, A scheme for interpolation with trigonometric spline curves, J. Comput. Appl. Math. 263 (2014) 246–261. P. Costantini, R.T. Farouki, C. Manni, A. Sestini, Computation of optimal composite re-parameterizations, Comput. Aided Geom. Design 18 (2001) 875–897.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
An iterative algorithm for G2 multiwise merging of Bézier curves Lizheng Lu ∗ , Chengkai Jiang School of Statistics and Mathematics, Zhejiang Gongshang University, Hangzhou 310018, China
article
abstract
info
This paper presents an iterative algorithm for G2 multiwise merging of Bézier curves. By using the G2 constraint, the L2 distance is represented after simplification as a quartic polynomial in two parameters relating to the magnitudes of end tangents of the merged curve. These two parameters are restricted in a feasible region, in order for the merged curve to preserve the specified directions of end tangents. Then G2 multiwise merging is formulated as a constrained minimization problem, and the classic projected Newton method is applied to find the minimizer. Some extensions of multiwise merging using G3 constraints, other energy functionals and curve representations are also outlined. Several comparative examples are provided to demonstrate the effectiveness of the proposed method. © 2015 Elsevier B.V. All rights reserved.
Article history: Received 4 May 2015 Received in revised form 4 October 2015 Keywords: Bézier curve Merging Geometric continuity L2 distance Iterative method
1. Introduction Assume that the given are K (≥ 2) segments of adjacent Bézier curves of degrees nk , k = 1, . . . , K , each of which is defined in terms of the control points Pknk := (pk0 , . . . , pknk )T as Pk (t ) =
nk
n
Bi k (t )pki = Bnk Pknk ,
t ∈ [0, 1].
(1)
i=0
Throughout, the notation Bn := Bn0 (t ), . . . , Bnn (t ) is used to denote the row vector of the Bernstein basis of degree n, with
Bni (t ) = i t i (1 − t )n−i . The problem of multiwise merging [1,2] aims to find a single degree-n Bézier curve (called merged curve) with control points Rn := (r0 , . . . , rn )T ,
n
R(t ) =
n
Bni (t )ri = Bn Rn ,
t ∈ [0, 1],
(2)
i=0
such that the L2 distance E=
K k=1
∗
tk+1 tk
2 K R(t ) − Pk t − tk dt = ∆ t k ∆t k
Corresponding author. E-mail address: [email protected] (L. Lu).
http://dx.doi.org/10.1016/j.cam.2015.10.007 0377-0427/© 2015 Elsevier B.V. All rights reserved.
k=1
1 0
∥Bn Sk Rn − Bnk Pknk ∥2 dt
(3)
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
353
is minimized. Here, 0 = t1 < t2 < · · · < tK +1 = 1 is a partition of the interval [0, 1], ∆tk = tk+1 − tk , and Sk Rn represents the control points of the sub-curve R(t ), t ∈ [tk , tk+1 ] after it is reparametrized to the standard interval [0, 1]. See [1] for the expression of the (n + 1) × (n + 1) matrix Sk . In geometric modeling, a complex shape is generally constructed by using multiple segments, with the advantage of producing more desirable shapes through local adjustment. It is often required to merge multiple adjacent segments into a single one, while preserving as faithful as possible to the original shape. For this goal, the merging problem is studied as a minimization problem using some suitable distances such as (3). To meet certain smoothness requirement in various applications, it is important to impose continuity conditions on the two endpoints of the merged curve R(t ), which correspond to the initial endpoint of the first segment P1 (t ) and the terminal endpoint of the last segment PK (t ). As multiple segments are merged at a time, multiwise merging does not need to consider continuity conditions at the interior endpoints. It is well known [3,4] that in contrast to the classical parametric continuity (C k ), geometric continuity (Gk ) shares additional parameters that are available for further shape manipulation by direct assignment or optimization techniques. Noting the fact that geometric continuity is a relaxation of parametrization (but not a relaxation of smoothness), it is more suitable to use geometric continuity for the merging problem since the parametrization of unknown merged curve should not be simply taken as that of the first (or last) segment. For instance, G2 continuity (mostly used in applications) is able to retain the positions, tangent directions and curvatures. 1.1. Related work Most of the existing merging approaches can be categorized into matrix-based approaches. The merging problem was firstly introduced in 1987 by Hoschek [5], where many segments of curves are approximately converted into a spline curve through degree reduction and splitting. Although degree reduction was widely investigated by many researchers, merging has not received adequate attention until Hu et al. [6] proposed a pioneering work of pairwise merging where a pair of Bézier curves is merged into one segment instead of applying degree reduction twice. Then, merging has become a special research topic. In terms of matrix representations and operations, many approaches have arisen, including pairwise merging of B-spline curves [7] and Bézier curves [8,9], and multiwise merging of Bézier curves [10,1] and B-spline surfaces [11]. In particular, C r ,s multiwise merging in [1] can preserve C r continuity and C s continuity at the two endpoints respectively. A novel dual approach was proposed by Woźny et al. [2] for C r ,s multiwise merging, which makes use of the advantages of constrained dual Bernstein basis [12,13]. Remarkably, this approach has lower complexity due to the fast evaluation schemes of certain connections between Bernstein and dual Bernstein polynomials. 1.2. Contributions In this paper, we revisit the multiwise merging problem with the constraint of G2 or G3 continuity, and propose an iterative algorithm based on the projected Newton method. According to the G2 constraint, there are four parameters (i.e., α1 , α2 , β1 , β2 in (5)–(8)), and two more parameters by the G3 constraint. Similar to the degree reduction work of Lu [14], we also find that for G2 multiwise merging, α2 and β2 are both expressed as quadratic polynomials of α1 and β1 when the minimum of the L2 distance (3) is reached; therefore, the number of independent parameters is reduced to 2, and the L2 distance is a quartic polynomial in two unknowns. It is important to note that many geometric modeling applications often require the two parameters α1 and β1 to be positive, so that the end tangents of the merged curve coincide with the specified directions, referring to (5) and (7). On the other hand, the values of α1 and β1 indicate the magnitudes of the end tangents. For these reasons, we impose a feasible region on (α1 , β1 ):
D = (α1 , β1 ) ∈ R2 : 0 < l0 ≤ α1 ≤ u0 , 0 < l1 ≤ β1 ≤ u1 .
(4)
By restricting (α1 , β1 ) in the feasible region, the merging problem is changed to a constrained minimization problem which may be efficiently solved by the projected Newton method, as done in [14] for degree reduction. Our contributions are summarized as follows:
• An iterative algorithm is proposed for G2 multiwise merging with the feasible region. • An intuitive interaction tool is provided for shape-preserving merging, in terms of the feasible region. Compared to the previous C r ,s multiwise merging [1,2], our method produces better results due to the used G2 and G3 constraints. Our method can be regarded as a generalization of the previous pairwise merging [8] which uses C 1 G2 and C 1 G3 constraints. Very recently, Gospodarczyk and Woźny [15] proposed a new method for C r ,s multiwise merging by restricting all the interior control points of the merged curve in a specified region; whereas our method restricts (α1 , β1 ) in a specified region. 2. Iterative algorithm for G 2 multiwise merging Given K segments of adjacent Bézier curves Pk (t ) with control points Pknk , k = 1, . . . , K , G2 multiwise merging amounts
to identifying the control points Rn of the merged curve R(t ), such that the L2 distance is minimized and the G2 constraint at the endpoints is satisfied.
354
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
By the G2 constraint, the first three and last three control points in Rn are determined by [3,4]: r0 = p10 ,
n1
r1 = p10 +
r2 = p10 +
2n1 n
rn = pKnK ,
α1 ∆p10 +
α1 ∆p10 , n n1 (n1 − 1) n(n − 1) nK
rn−1 = pKnK − 2nK
rn−2 = pKnK −
n
(5)
α12 ∆2 p10 +
n1 n(n − 1)
α2 ∆p10 ,
(6)
β1 ∆pKnK −1 ,
n
β1 ∆pKnK −1 +
(7)
nK (nK − 1) n( n − 1 )
nK
β12 ∆2 pKnK −2 +
n(n − 1)
β2 ∆pKnK −1 ,
(8)
where α1 , α2 , β1 , β2 ∈ R, and the forward difference operator ∆ is defined by ∆pki = pki+1 − pki and by ∆r pki = ∆r −1 pki+1 − ∆r −1 pki for r > 1. The merged curve is rewritten as R(t ) = Bn0 (t )r0 + Bn1 (t )r1 + Bn2 (t )r2 + Bnn−2 (t )rn−2 + Bnn−1 (t )rn−1 + Bnn (t )rn +
n −3
Bni (t )ri .
i =3
Let Rcn := (r0 , r1 , r2 , rn−2 , rn−1 , rn )T and the L2 distance in (3) is expanded by E=
K
∆tk
:= (r3 , . . . , rn−3 )T be the constrained and free control points respectively. Then
1
0
k=1
f Rn
∥Bn Skc Rcn + Bn Skf Rfn − Bnk Pknk ∥2 dt
=(Rcn )T Hcc Rcn + 2(Rfn )T Hfc Rcn + (Rfn )T Hff Rfn − 2(Rcn )T
K
∆tk (Skc )T Gn,nk Pknk
k=1
− 2(Rfn )T
K
f
∆tk (Sk )T Gn,nk Pknk +
k=1
K
∆tk (Pknk )T Gnk ,nk Pknk
(9)
k=1
n ,m
where Gn,m = (gij )i=0,j=0 ∈ R(n+1)×(m+1) with entries 1
Bni (t )Bm j (t )dt =
gij = 0
1 n+m+1
n
m
n+m
i
j
i+j
and Hcc =
K
∆tk (Skc )T Gn,n Skc ,
k=1
Hfc =
K
f
∆tk (Sk )T Gn,n Skc ,
Hff =
k=1
K
f
f
∆tk (Sk )T Gn,n Sk ,
k=1 f
with Skc := Sk (· · · ; 0, 1, 2, n − 2, n − 1, n) and Sk := Sk (· · · ; 3, . . . , n − 3) denoting the submatrices of Sk formed by selecting the specified columns. f
At first glance, it is natural to determine α1 , β1 , α2 , β2 and Rn using (9) via min
f α1 ,β1 ,α2 ,β2 ,Rn
{E := E (α1 , β1 , α2 , β2 , Rfn )}.
(10)
However, as will be shown in the following subsection, a further analysis reveals that we are able to decrease the complexity of the problem.
2.1. Simplification The expression (9) can be simplified through three steps: 1. Express ri = ri (α1 , β1 , α2 , β2 ), i = 3, . . . , m − 3 as quadratic polynomials, which are quadratic in α1 , β1 and linear in α2 , β2 . 2. Express α2 = α2 (α1 , β1 ) and β2 = β2 (α1 , β1 ) as quadratic polynomials of α1 , β1 . 3. Express E = E (α1 , β1 ) as a quartic polynomial of α1 , β1 .
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
355
f
Firstly, assuming that Rcn and Rn are variables, it follows from (9) that
K ∂E ∂E ∂E ∂E ∂E ∂E ∂E T c T f := , , , , , = 2H R + 2H R − 2 ∆tk (Skc )T Gn,nk Pknk , cc n fc n ∂ Rcn ∂ r0 ∂ r1 ∂ r2 ∂ rn−2 ∂ rn−1 ∂ rn k =1 K ∂E ∂E ∂E T f f c := , . . . , = 2H R + 2H R − 2 ∆tk (Sk )T Gn,nk Pknk . ff fc n n ∂ r3 ∂ rn−3 ∂ Rfn k=1
(11)
(12)
By (12), Rfn = −Hff−1 Hfc Rcn + Hff−1
K
f
∆tk (Sk )T Gn,nk Pknk
(13)
k=1
which implies that ri , i = 3, . . . , m − 3 are quadratic polynomials of α1 , β1 , α2 , β2 . Secondly, note the following facts: r1 is linear in α1 ; r2 is quadratic in α1 and linear in α2 ; rn−1 is linear in β1 ; rn−2 is quadratic in β1 and linear in β2 . After substituting (13) into (9), E is a quartic polynomial of α1 , β1 , α2 , β2 and takes the form E (α1 , β1 , α2 , β2 ) =
1 2
α2
β2
a00
a01 a11
a10
α2 + α2 b0 (α1 , β1 ) + β2 b1 (α1 , β1 ) + c (α1 , β1 ) β2
(14)
where the constants aij and quadratic polynomials bi (α1 , β1 ) will be derived below, and c (α1 , β1 ) = E (α1 , β1 , 0, 0) by (9) and (13). From (5)–(8) and (13), it is easy to show that
∂ Rcn = 0 ∂α2 ∂ Rcn = 0 ∂β2
T
,
T
,
0
v0
0
0
0
0
0
v1
0
0
n
∂ Rfn = −Hff−1 Hfc 0 ∂α2 ∂ Rfn = −Hff−1 Hfc 0 ∂β2
T
,
T
,
0
v0
0
0
0
0
0
v1
0
0
n
K 1 ∆p10 and v1 = n(n− ∆pKnK −1 . Differentiating by the chain rule gives where v0 = n(n− 1) 1)
∂E = 2A · 0 ∂α2
0
v0
0
0
T
0
,
∂E = 2A · 0 ∂β2
0
0
v1
0
T
0
,
(15)
where A = (Hcc − HfcT Hff−1 Hfc )Rcn −
K
∆tk (Skc )T Gn,nk Pknk + HfcT Hff−1
K
f
∆tk (Sk )T Gn,nk Pknk .
k=1
k=1
Using (14) and (15), we obtain b0 (α1 , β1 ) =
∂ E (α1 , β1 , 0, 0) , ∂α2
b1 (α1 , β1 ) =
∂ E (α1 , β1 , 0, 0) ∂β2
(16)
and
∂ 2E ∂α 2 2 = ∂ 2E ∂β2 ∂α2
a00 a10
a01 a11
∂ 2E 2w22 v0 · v0 ∂α2 ∂β2 = 2w32 v1 · v0 ∂ 2E
2w23 v0 · v1 2w33 v1 · v1
(17)
∂β22
5,5 ∂E ∂E where wij are the entries of the 6 × 6 matrix W = wij i=0,j=0 := Hcc − HfcT Hff−1 Hfc . By applying ∂α = ∂β = 0 to (14), α2 2 2 and β2 are solved as α2 =
−a11 b0 (α1 , β1 ) + a01 b1 (α1 , β1 ) , a00 a11 − a201
β2 =
a10 b0 (α1 , β1 ) − a00 b1 (α1 , β1 ) a00 a11 − a201
,
(18)
both of which are linear combinations of quadratic polynomials b0 , b1 and are thus quadratic polynomials of α1 , β1 . Thirdly, after substituting (18) into (14), E is a quartic polynomial of α1 , β1 . Based on the above three steps, we finally arrive at the constrained minimization problem: min {Eˆ (α1 , β1 ) := E (α1 , β1 , α2 (α1 , β1 ), β2 (α1 , β1 ))}.
(α1 ,β1 )∈D
(19)
356
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
2.2. Numerical solution To solve the constrained minimization problem (19), we resort to some numerical methods. Since the objective function is quartic in two unknowns, it is trivial to compute its gradient and Hessian matrix. The projected Newton method [16] is employed, which is an iterative method and achieves a quadratic convergence rate. For numerical optimization, it is advisable to start from an initial point that is not far away from a local minimum. If (1, 1) ∈ D , (α1 , β1 ) = (1, 1) is initially assigned, which corresponds to the C 1 continuity conditions; otherwise, the initial value of (α1 , β1 ) can be taken as the central point of D . 3. Extensions G2 multiwise merging described detailedly in Section 2 can be extended in several ways, with similar derivations. Three directions are outlined in the following. 3.1. G3 constraint Other constraints of geometry continuity can be used, e.g., G1,2 constraint, G2,3 constraint, G3 constraint, etc. We do not suggest to use higher-order geometric continuity of order greater than 3, because G2 continuity and G3 continuity have sufficient smoothness properties for practical applications and more degrees of freedom should be left to the interior control points of merged curves for generating more satisfactory results. For example, G3,2 multiwise merging is briefly described below. In this case, R(t ) and P1 (t ) are G3 continuous at t = 0, while R(t ) and PK (t ) are G2 continuous at t = 1; hence, besides (5)–(8), the control point r3 is expressed using a new parameter α3 ∈ R as 3n1
n 3n1
α1 ∆p10 +
3n1 (n1 − 1)
α12 ∆2 p10 +
n1 (n1 − 1)(n1 − 2)
α13 ∆3 p10 n(n − 1)(n − 2) n1 α2 ∆p10 + α1 α2 ∆2 p10 + α3 ∆p10 . + n( n − 1 ) n(n − 1)(n − 2) n(n − 1)(n − 2)
r3 = p10 +
n(n − 1) 3n1 (n1 − 1)
(20)
In a similar way, the L2 distance is simplified through three steps: 1. Express ri = ri (α1 , α2 , β1 , α3 , β2 ), i = 4, . . . , m − 3 as cubic polynomials, which are cubic in α1 , quadratic in β1 and linear in α2 , α3 , β2 . 2. Express α3 = α3 (α1 , α2 , β1 ) and β2 = β2 (α1 , α2 , β1 ) as cubic polynomials of α1 , α2 , β1 . 3. Express E = E (α1 , α2 , β1 ) as a sextic polynomial of α1 , α2 , β1 . Then, a similar constrained minimization problem like (19) is built, but now the objective function is sextic in three unknowns. 3.2. Energy functional Obviously, the l2 distance in [10,6,8,1] can be used as the energy functional. If geometric information of higher-order derivatives of the given curves is considered, then the L2 distance in (3) is further modified as E˜ =E +
2
ωℓ
ℓ=1
K k=1
tk+1 tk
ℓ 2 d dℓ k t − tk dt R ( t ) − P dt ℓ dt ℓ ∆t k
2 1 K d d 1 k Bn Sk Rn − Bnk Pnk dt = ∆tk ∥Bn Sk Rn − ∥ + ω1 ∆tk 0 dt dt 0 k=1 k=1 K 2 1 2 2 1 d Bn Sk Rn − d Bn Pk dt + ω2 n k k (∆tk )3 0 dt 2 dt 2 k=1 K
1
Bnk Pknk 2 dt
(21)
where the weights ω1 , ω2 ≥ 0 control relative weighting of the three involved terms. Note that E˜ is identical to E when ω1 = ω2 = 0 and that the second derivatives are not included in E˜ when ω1 > 0 and ω2 = 0. We do not suggest to include ˜ since more weights would be introduced that unavoidably increase the difficulty of third (or higher-order) derivatives in E, their adjustments.
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
357
The whole process of Section 2 can be kept to minimize the new energy functional (21), but with the modification that n ,m n ,m for every matrix Gn,m = (gij )i=0,j=0 ∈ R(n+1)×(m+1) , two corresponding matrices Gℓn,m = (gijℓ )i=0,j=0 ∈ R(n+1)×(m+1) , ℓ = 1, 2 are additionally required, where gijℓ =
ℓ m dℓ Bni (t ) d Bj (t )
1
dt ℓ
0
dt ℓ
dt .
Although Gn,n (called the Gram matrix of the Bernstein basis) is symmetric positive definite [1], G1n,n and G2n,n are both n d ℓ n dℓ n dℓ n symmetric positive semidefinite since dt ℓ B0 (t ), . . . , dt ℓ Bn (t ) are linearly dependent by i=0 dt ℓ Bi (t ) ≡ 0, ∀t ∈ [0, 1]. 3.3. Curve representation Let ( u0 (t ), . . . , un (t )) defined on the interval [a, b] be a basis of a function space, and represent the curve in this space as n C(t ) = i=0 ui (t )pi , t ∈ [a, b]. Many kinds of bases and curve representations have been proposed for various applications in CAGD. As is well known, the Bernstein–Bézier representation has gained widespread use in CAGD and related fields, mainly because it exhibits elegant geometric properties [3,17]. Although we formulate multiwise merging in terms of Bézier curves, multiwise merging of curves in other representations can also be solved in the same manner. For the extended Chebyshev space
Tα2n = span{cos(kt ), sin(kt ) : t ∈ [0, α]}nk=0 ,
α ∈ (0, π )
of trigonometric polynomials of degree 2n (order n), a trigonometric curve of degree 2n is defined as the convex combination (see [18]) Cαn (t ) =
2n
α T2n ,i (t )pi ,
t ∈ [0, α],
(22)
i=0
where α
α
T2n,i (t ) = τ2n,i sin
2n−i
α−t 2
i
sin
i
t
2
,
α
τ2n,i
⌊2⌋ α i−2r n i−r := 2 cos 2 r sin2n α2 r =0 i − r 1
are the trigonometric normalized B-basis functions. When it comes to implementing multiwise merging of trigonometric curves, matrices Sk are specified by a de Casteljau-type algorithm [19,20], whereas all the entries α α the subdivision 2n,2m α gij = 0 T2n ,i (t )T2m,j (t )dt of G2n,2m = (gij )i=0,j=0 in this case must be evaluated via numerical integration and can be stored in lookup tables in advance. Finally, we would like to mention that multiwise merging can be performed on hyperbolic curves [20,21], algebraictrigonometric curves [22,23], algebraic-hyperbolic curves [23], etc. Note that when the shape parameter α in (22) approaches 0, the trigonometric curve becomes a Bézier curve of degree 2n, see [24, Proposition 2.1]; a similar asymptotic behavior occurs also in the hyperbolic curve case, see [21, Proposition 5]. 4. Examples and applications In this section, we present some examples for the proposed Gk multiwise merging and compare with the previous C k multiwise merging [1,2]. All examples are computed on a PC with Intel Core2 Quad CPU Q9400 2.66 GHz processor and 4 GB RAM. The merging results depend on the choice of the parameters {tk }Kk=+11 . Let t1 = 0, and assign λ,µ,ω
tk+1 = tk +
Vk
V λ,µ,ω
,
k = 1, . . . , K ,
(23)
where λ,µ,ω
Vk
1
:= 0
d k λ P (t ) dt + ω dt
1
|κk (t )|µ dt and V λ,µ,ω :=
0
K
λ,µ,ω
Vk
k=1
with
d k 3 d k d2 k κk (t ) = P (t ) × 2 P (t ) P (t ) dt dt dt denoting the curvature of the kth segment. If ω = 0, it is the commonly used arc length parametrization when λ = 1, the uniform parametrization when λ = 0, and the centripetal parametrization when λ = 1/2. However, if ω > 0, curvature
358
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
.
.
.
.
.
.
Fig. 1. Merging results for the treble clef curve using C k and Gk constraints.
Table 1 Comparative data between C k and Gk multiwise merging for the treble clef curve. (Note: the computation time is measured in milliseconds.) n
C2
C3
G2
G3
E
Time
E
Time
E
Iters
Time
E
Iters
Time
12 14 16 18 20
1.225e−2 3.903e−3 2.145e−3 1.063e−3 5.480e−4
18 21 24 28 32
4.966e−2 2.100e−2 8.801e−3 5.199e−3 2.788e−3
16 18 21 25 29
2.077e−3 6.545e−4 2.434e−4 1.667e−4 8.289e−5
5 6 5 5 6
30 32 35 40 43
2.379e−3 6.708e−4 2.573e−4 2.039e−4 8.940e−5
20 18 21 19 21
73 76 82 83 85
variations of the given Bézier segments are considered as well; this choice may provide aesthetically more pleasing results in case of high curvature segments. The first example is the treble clef curve composed of four Bézier curves of degree 5; see [1] for detailed control points. Fig. 1 illustrates the merged curves of degrees n, in which {t1 , t2 , . . . , t5 } = {0, 0.2387, 0.4872, 0.7626, 1} is taken by arc length parametrization. When applying G2 and G3 multiwise merging, the feasible region is set to D = [0.2, 5]2 . Clearly, owing to the used Gk constraint, better results are generated by Gk multiwise merging. As shown in Table 1, Gk multiwise merging requires more computation time than C k multiwise merging, caused by the three simplification steps and iterative process for solving the constrained minimization problem; and, compared to G2 multiwise merging, G3 multiwise merging increases the number of iterations and thus the computation time, due to the iterative process for solving a more complex problem. The second example is the penguin curve, whose left and right parts are composed of four and three cubic Bézier curves respectively; see [2] for detailed control points. Fig. 2 illustrates the merged curves of degrees n, and D = [0.2, 5]2 is used for the feasible region. Apparently, Gk multiwise merging greatly improves the results. In Fig. 2(a)–(d), {t1 , t2 , . . . , t5 } = {0, 0.08, 0.55, 0.78, 1} for the left part and {t1 , t2 , t3 , t4 } = {0, 0.42, 0.78, 1} for the right part are taken by arc length parametrization. In Fig. 2(e)–(h), {t1 , t2 , . . . , t5 } = {0, 0.18, 0.53, 0.78, 1} for the left part and {t1 , t2 , t3 , t4 } = {0, 0.40, 0.75, 1} for the right part are taken by setting (λ, µ, ω) = (1, 0.5, 0.05) in (23): the results by C k multiwise merging are noticeably improved, referring to the segment at the left bottom; whereas the improvement of Gk multiwise merging is not obvious, possibly owing to the used Gk constraint. We show an application to improve the parametrization of merged curves by using the feasible region. A curve is thought of as having optimal parametrization if it is parametrized by arc length [25], but which is impossible for polynomial curves. While for a merged curve R(t ), the magnitudes of all its tangent vectors, i.e. ∥R′ (t )∥, should be as close to equal as possible in order to achieve an optimal parametrization. Here, we do not intend to find a merged curve with the closest optimal parametrization; instead, we improve the parametrization by interactively adjusting the four bounds in the feasible region. Recalling that these four bounds indirectly specify the ranges of ∥R′ (0)∥ and ∥R′ (1)∥, such interaction is meaningful and flexible. An example is shown in Fig. 3. By making use of the two weights in (21), it increases the flexibility of shape design. As the first term of (21) measures the deviation in positions, small weights may be used for generating satisfactory results from the viewpoint of curve profiles; on the other hand, large weights should be taken if higher-order differential geometric quantities are regarded as important. However, a qualitative criterion for determining suitable values of weights still seems to be unattainable; alternatively, we need to manually assign and adjust their values. As shown in Fig. 4, by setting ω1 = ω2 = 1, the curvature distributions of
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
(a) n = 14.
(b) n = 14.
(c) n = 16.
(d) n = 16.
(e) n = 14.
(f) n = 14.
(g) n = 16.
(h) n = 16.
359
Fig. 2. Merging results for the penguin curve using C k and Gk constraints.
merged curves are improved. We noticed that the improvement is not obvious when the given curves are not C 2 connected, such as the treble clef curve and penguin curve. Finally, Fig. 5 demonstrates an example where multiwise merging is performed on three Bézier curves and three trigonometric curves with α = π /2. Since different curve representations are used, there exist slight differences between the merging results. 5. Conclusion We have developed an iterative algorithm for Gk multiwise merging of Bézier curves. The key idea of the proposed method is to make use of the additional parameters afforded by the Gk constraint, with the ability to produce better merging results. Moreover, the added flexibility can be exploited in many ways, and constitutes an important advance in geometric modeling applications, e.g., improving the parametrization of curves. We also propose to use the feasible region for the purpose of shape preservation and interaction. Acknowledgments This work is supported by the National Natural Science Foundation of China (Nos. 61272307, 61502128). The authors thank the anonymous referees for their constructive comments and for pointing out two extension directions (Sections 3.2 and 3.3).
360
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
(a) n = 9, D = [4, 5] × [2, 2.5].
(b) n = 10, D = [3.5, 4.5] × [3.5, 4.5]. 2
Fig. 3. Applying G multiwise merging to improve the parametrization of merged curves by using the feasible region D . The middle column displays the profiles of ∥R′ (t )∥ with respect to t, while the means and standard deviations of {∥R′ (t¯i )∥}Ni=1 evenly sampled in t are listed in the right column.
(a) n = 10, D = [0.2, 5]2 .
(b) n = 12, D = [0.2, 5]2 .
Fig. 4. Influence of the weights ω1 and ω2 : Red, ω1 = ω2 = 0; Blue, ω1 = ω2 = 1. In (a) and (b), the merged curves with G2 continuity are displayed on the left, while the profiles of curvatures with respect to arc length on the right. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(a) Bézier curves.
(b) Trigonometric curves.
Fig. 5. Merging results for three Bézier curves and three trigonometric curves with α = π/2, using D = [0.2, 5]2 . Note that the control points of the given curves in (a) and (b) are the same.
L. Lu, C. Jiang / Journal of Computational and Applied Mathematics 296 (2016) 352–361
361
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25]
L. Lu, Explicit algorithms for multiwise merging of Bézier curves, J. Comput. Appl. Math. 278 (2015) 138–148. P. Woźny, P. Gospodarczyk, S. Lewanowicz, Efficient merging of multiple segments of Bézier curves, Appl. Math. Comput. 268 (2015) 354–363. G. Farin, Curves and Surfaces for CAGD: A Practical Guide, fifth ed., Morgan Kaufmann, San Francisco, 2002. J. Peters, Geometric continuity, in: G. Farin, J. Hoschek, M.-S. Kim (Eds.), Handbook of Computer Aided Geometric Design, Elsevier, Amsterdam, 2002, pp. 193–227. J. Hoschek, Approximate conversion of spline curves, Comput. Aided Geom. Design 4 (1987) 59–66. S. Hu, R. Tong, T. Ju, J. Sun, Approximate merging of a pair of Bézier curves, Comput. Aided Des. 33 (2001) 125–136. C. Tai, S. Hu, Q. Huang, Approximate merging of B-spline curves via knot adjustment and constrained optimization, Comput. Aided Des. 35 (2003) 893–899. L. Lu, An explicit method for G3 merging of two Bézier curves, J. Comput. Appl. Math. 260 (2014) 421–433. P. Zhu, G. Wang, Optimal approximate merging of a pair of Bézier curves with G2 -continuity, J. Zhejiang Univ. Sci. A 10 (2009) 554–561. M. Cheng, G. Wang, Approximate merging of multiple Bézier segments, Progr. Nat. Sci. 18 (2008) 757–762. H. Pungotra, G.K. Knopf, R. Canas, Merging multiple B-spline surface patches in a virtual reality environment, Comput. Aided Des. 42 (2010) 847–859. S. Lewanowicz, P. Woźny, Bézier representation of the constrained dual Bernstein polynomials, Appl. Math. Comput. 218 (2011) 4580–4586. P. Woźny, S. Lewanowicz, Multi-degree reduction of Bézier curves with constraints, using dual Bernstein basis polynomials, Comput. Aided Geom. Design 26 (2009) 566–579. L. Lu, Some improvements on optimal multi-degree reduction of Bézier curves with geometric constraints, Comput. Aided Des. 59 (2015) 39–42. P. Gospodarczyk, P. Woźny, Merging of Bézier curves with box constraints, J. Comput. Appl. Math. (2016) http://dx.doi.org/10.1016/j.cam.2015.10.005. to appear. M. Schmidt, D. Kim, S. Sra, Projected Newton-type methods in machine learning, in: S. Sra, S. Nowozin, S.J. Wright (Eds.), Optimization for Machine Learning, MIT Press, Cambridge, MA, 2011, pp. 305–330. J.M. Peña, Shape Preserving Representations in Computer-Aided Geometric Design, Nova Science Publishers, Commack, New York, 1999. J. Sánchez-Reyes, Harmonic rational Bézier curves, p-Bézier curves and trigonometric polynomials, Comput. Aided Geom. Design 15 (1998) 909–923. E. Mainar, J.M. Peña, Corner cutting algorithms associated with optimal shape preserving representations, Comput. Aided Geom. Design 16 (1999) 883–906. Á Róth, Control point based exact description of trigonometric/hyperbolic curves, surfaces and volumes, J. Comput. Appl. Math. 290 (2015) 74–91. W. Shen, G. Wang, A class of quasi Bézier curves based on hyperbolic polynomials, J. Zhejiang Univ. Sci. 6A (Suppl. I) (2005) 116–123. J.M. Carnicer, E. Mainar, J.M. Peña, Critical length for design purposes and extended Chebyshev spaces, Constr. Approx. 20 (2004) 55–71. E. Mainar, J.M. Peña, Optimal bases for a class of mixed spaces and their associated spline spaces, Comput. Math. Appl. 59 (2010) 1509–1523. I. Juhász, Á Róth, A scheme for interpolation with trigonometric spline curves, J. Comput. Appl. Math. 263 (2014) 246–261. P. Costantini, R.T. Farouki, C. Manni, A. Sestini, Computation of optimal composite re-parameterizations, Comput. Aided Geom. Design 18 (2001) 875–897.