Rectifying control polygon for planar Pythagorean hodograph curves

Computer Aided Geometric Design 54 (2017) 1–14

Contents lists available at ScienceDirect

Computer Aided Geometric Design www.elsevier.com/locate/cagd

Rectifying control polygon for planar Pythagorean hodograph curves ✩ Soo Hyun Kim, Hwan Pyo Moon ∗ Department of Mathematics, Dongguk University-Seoul, Seoul, 04620, Republic of Korea

a r t i c l e

i n f o

Article history: Received 14 July 2016 Received in revised form 24 March 2017 Accepted 24 March 2017 Available online 31 March 2017 Keywords: Pythagorean-hodograph curve Bézier control polygon Rectifying control polygon Gauss–Legendre quadrature Bernstein–Vandermonde matrix

a b s t r a c t A Bézier control polygon is not appropriate to control a Pythagorean hodograph curve since it has redundant degrees of freedom. So we propose an alternative, which is the rectifying control polygon. A rectifying control polygon of a PH curve has the same degrees of freedom as the PH curve. It interpolates the end points of the PH curve, but not the end tangents. Most of all, it has the same arc length as the PH curve. In this paper, we present the method to compute the rectifying control polygon from the Bézier control polygon of the PH curve. We also present the procedure to compute the PH curves from a given rectifying control polygon. For the development of these algorithms, we employ the Gauss– Legendre quadrature method and the Bernstein–Vandermonde linear system. © 2017 Elsevier B.V. All rights reserved.

1. Introduction Pythagorean hodograph (PH) curves, introduced by Farouki and Sakkalis (1990), are a special type of polynomial curves, which have the unique property that their parametric speed functions are also polynomials of the curve parameter. The PH property enables us to compute the arc length of the curve exactly without numerical integration. The curve construction problem with prescribed arc length has been studied by Roulier (1993), Roulier and Piper (1996) using Bézier curves or even general parametric curves. Their methods rely on iterative numerical computations. Similar problems using PH curves were solved by Farouki (2016), Huard et al. (2014) in much simpler and elegant ways because PH curves have explicit expression of the arc length. Another important advantage of PH curves is that their offset curves are rational curves. So we do not need to rely on approximation algorithms for offset computation. One may consult Farouki (2007) for more details on PH curves from algebraic frameworks to practical applications. Since PH curves are polynomial curves, they can be expressed as Bézier curves. However, not all Bézier curves are PH curves. Recently, Farouki et al. (2015) suggest the method to determine whether a given Bézier curve is in fact a PH curve, and to compute the parameters of the PH curve. They also presented the method of local modification of quintic PH spline curve while maintaining the PH property (Farouki et al., 2016). Since the Bézier control points of a PH curve should satisfy certain algebraic constraints, we cannot move a Bézier control point arbitrarily while preserving the PH condition. Actually a PH curve has fewer degrees of freedom than a Bézier curve of the same degree has. This situation suggests that a Bézier control polygon is not appropriate to “control” a PH curve. So we propose an alternative method to control PH curves, which is the rectifying control polygon.

✩

*

This paper has been recommended for acceptance by Rida Farouki. Corresponding author. E-mail addresses: [email protected] (S.H. Kim), [email protected] (H.P. Moon).

http://dx.doi.org/10.1016/j.cagd.2017.03.016 0167-8396/© 2017 Elsevier B.V. All rights reserved.

2

S.H. Kim, H.P. Moon / Computer Aided Geometric Design 54 (2017) 1–14

A rectifying control polygon of a PH curve will be defined as a polygon which satisfies the properties: (1) it has the same degrees of freedom as the PH curve has, (2) it interpolates both end points of the PH curve, and (3) the sum of the lengths of polygon segments is the same as the length of the PH curve. The third condition above is the rectifying property of this polygon. To achieve the rectifying property, we will utilize the Gauss–Legendre quadrature for the computation of the arc length of a PH curve. The Gauss–Legendre quadrature with n nodes produces the exact integral of a polynomial of degree 2n − 1. We can construct a rectifying control polygon of a PH curve by attaching the tangent vectors of the PH curve, evaluated at the Gauss–Legendre nodes and scaled by the Gauss–Legendre weights. One may find the details of the Gauss–Legendre quadrature in basic texts of numerical analysis such as Atkinson (1989), Atkinson and Han (2004). The rest of the paper is structured as follows. In Section 2, we summarize the fundamental properties of PH curves and the Gauss–Legendre quadrature, then fix the notations. Section 3 is devoted to the construction of a rectifying polygon from a given PH curve. We here prove that the vertices of the rectifying polygon can be obtained by the convex combinations of the Bézier control points of the PH curve. In Section 4, we solve the inverse problems of Section 3. We construct PH curves from a given rectifying control polygon. In fact, there are multiple instances of PH curves for a given rectifying control polygon. We describe the procedure to compute these PH curves. This procedure can be formulated as a linear system whose coefficient is the Bernstein–Vandermonde matrix evaluated at the Gauss–Legendre nodes. Finally we conclude the paper in Section 5. 2. Preliminaries When we deal with planar geometry, it is convenient to employ the complex number field. A point z = (x, y ) in R2 is identified with the complex number z = x + i y in C. Similarly, a planar parametric curve p(t ) = (x(t ), y (t )) can be identified with a complex valued function p(t ) = x(t ) + i y (t ). A Bézier curve p of degree n with complex Bézier control points bk = xk + i yk for k = 0, · · · , n, can then be expressed as

p(t ) =

n 

B nk (t )bk ,

t ∈ [0, 1]

(1)

k =0

where B nk (t ) is the Bernstein polynomial defined by

 

B nk (t ) =

n k

(1 − t )n−k t k .

Let bk denote the forward difference of the k-th control point, i.e., bk = bk+1 − bk . Then the derivative of the curve p is given by

p (t ) = n

n −1 

B nk −1 (t )bk .

k =0

A planar polynomial curve p(t ) = x(t ) + i y (t ) is called a Pythagorean hodograph (PH) curve (Farouki and Sakkalis, 1990) if and only if its derivative p (t ) = x (t ) + i y  (t ) satisfies the Pythagorean condition

x (t )2 + y  (t )2 = σ (t )2 for some polynomial σ (t ). The algebraic structures of Pythagorean hodograph curves have been studied in many different contexts. We here refer to a few important results. Theorem 1. A planar polynomial curve p(t ) = x(t ) + i y (t ) of odd degree is a Pythagorean hodograph curve if and only if there exist polynomials u (t ) and v (t ) which satisfies

x (t ) = u (t )2 − v (t )2 , y  (t ) = 2u (t ) v (t ). Theorem 2. A planar polynomial curve p(t ) = x(t ) + i y (t ) of odd degree is a Pythagorean hodograph curve if and only if there exists a complex valued polynomial z(t ) such that

p (t ) = z(t )2 . One may consult Farouki (1994) for the proof of the above theorems. For a regular PH curve p(t ) with p (t ) = 0, the polynomial σ (t ) can be chosen as a positive polynomial which is the speed function of p(t ), i.e.,

S.H. Kim, H.P. Moon / Computer Aided Geometric Design 54 (2017) 1–14

3



σ (t ) = p (t ) = x (t )2 + y  (t )2 . The arc-length L of the PH curve over [0, 1] is then

1 L=

σ (t )dt , 0

which can be computed exactly without numerical approximation, since σ (t ) is a polynomial function. If p(t ) is a PH curve of degree n, then σ (t ) is a polynomial of degree n − 1. Thus it can be expressed as a Bernstein polynomial

σ (t ) =

n −1 

σk B nk −1 (t )

k =0

with the coefficients σk ’s. Since the integral of the Bernstein polynomial is arc-length L of the PH curve p(t ) is

1 0

B nk −1 (t )dt = 1/n for k = 0, · · · , n − 1, the

n −1

L=

1 n

σk ,

k =0

which is the average of the coefficients σk ’s. The arc length of a PH curve can be obtained by the integral of a polynomial function, and we can compute the exact value of the polynomial integral by using quadrature rules such as Gauss–Legendre quadrature. Let f (t ) be an integrable function defined on the interval [−1, 1]. A quadrature method is the way to approximate the integral of f (t ) as a weighted sum of values of f (t ) at some specific nodes. Thus the quadrature Q ( f ) with nodes τk and weights ωk for k = 0, · · · , n is

Q(f)=

n 

ωk f (τk ).

k =0

The Gauss–Legendre quadrature is a special type of quadrature which evaluates the exact integral of the polynomials of the maximal degree with the given number of nodes. For the Gauss–Legendre quadrature with n + 1 node points, the nodes τk are chosen as the roots of the Legendre polynomials P n+1 (t ), which are defined as follows:

P 0 (t ) = 1

n  (−1)n dn  2 1 − t , n = 1, 2, · · · , n!2n dt n and the corresponding weights ω0 , · · · , ωn are given by P n (t ) =

ωk =

−2 , n+1 (τk ) P n+2 (τk )

(n + 2) P 

k = 0, · · · , n .

We adopt the formulation in Atkinson (1989), with slightly different indexing scheme for the clarity of later development. Our indexing scheme for the nodes and the weights starts from 0 instead of 1. Due to the orthogonality of Legendre polynomials, the Gauss–Legendre quadrature is the optimal quadrature method in the sense that it requires the minimal number of nodes to get the exact integral of polynomials of the maximal degree. According to the error analysis (Atkinson and Han, 2004; Kythe and Schäferkotter, 2005), the Gauss–Legendre quadrature with n + 1 pairs of nodes τk and weights ωk for k = 0, · · · , n produces the exact integral of polynomial of degree up to 2n + 1. The nodes {τk } and weights {ωk } of the Gauss–Legendre quadrature up to order 5 can be computed algebraically and listed (Farouki et al., 2015; Kythe and Schäferkotter, 2005) in Table 1. For higher order quadrature, those values should be computed numerically. Throughout this paper, τk stands for the Gauss–Legendre node in increasing order and ωk stands for the corresponding Gauss–Legendre weight. 3. Rectifying polygon of PH curves The arc-length formula of a differentiable curve is motivated by the piecewise linear approximation. Let p be a regular curve defined on [0, 1]. If we choose a partition 0 = t 0 < t 1 < · · · < tn+1 = 1 of [0, 1], then the polygon connecting the sequence of points p(t 0 ), p(t 1 ), · · · , p(tn+1 ) approximates the given curve p. The length of this polygon converges to the arc length of p as the partition size tends to 0. However, this process does not provide any specific polygon whose length is exactly the same as the arc length of p. Our goal in this section is to construct a polygon which realizes the exact arc length of a given PH curve. We here introduce the bracket notation [p0 p1 · · · pn+1 ] for the polygon which connects the sequence of points p0 , · · · , pn+1 .

4

S.H. Kim, H.P. Moon / Computer Aided Geometric Design 54 (2017) 1–14

Table 1 Nodes and weights of Gauss–Legendre quadrature up to order 5. No. of nodes

Nodes

1

0

±

2

Weights

ωk

2



1 3

0, ±

3

τk



3 7

4

±

5

0, ± 13

1 3 5

−

2 7



8 5 , 9 9

6 ,± 5



5−2

3 7

+ 27

10 , ± 13 7



√

√

18+ 30 18− 30 , 36 36

6 5



5+2

10 7

√

√

128 322+13 70 322−13 70 , , 225 900 900

Although one may work with PH curves of even degrees, see e.g., Wang and Fang (2009), we here present the result for the PH curve of odd degrees since they are generic types of PH curves. A PH curve p(t ) of degree 2n + 1 can be specified by a complex valued polynomial z(t ) such that p (t ) = z(t )2 and the starting point p(0). Since the polynomial z(t ) is of degree n, its complex degrees of freedom are n + 1. So the complex degrees of freedom of p(t ) are n + 2. Thus the Bézier control points b0 , · · · , b2n+1 have redundant degrees of freedom, and we propose the rectifying polygon, which consists of n + 2 points. Theorem 3. Let p be a PH curve of degree 2n + 1 with Bézier control points b0 , · · · , b2n+1 . The nodes and the weights of the Gauss– Legendre quadrature of order n + 1 are denoted by {τ0 , · · · , τn } and {ω0 , · · · , ωn } respectively. If we define

p0 = b0 , pk+1 = pk +

ωk 2

p



1 + τk



2

(2)

for k = 0, · · · , n,

,

then [p0 p1 · · · pn+1 ] is a rectifying polygon of p. Proof. The given PH curve p and its derivative p are expressed as

p(t ) =

2n +1 

+1 B 2n (t )b j j

and

p (t ) = (2n + 1)

j =0

2n 

B 2n j (t )b j ,

j =0

respectively. Since the Gauss–Legendre quadrature with n + 1 node points provides the exact integral of polynomial up to degree 2n + 1, we have n  ωk k =0

2

 B 2n j

1 + τk



1

2

1

B 2n j (t ) dt =

=

2n + 1

0

,

for each j = 0, · · · , 2n. We modified the nodes and the weights to change the domain of integration from [−1, 1] to [0, 1]. We now prove the end point agreement. From the recursive definition of pk and the Gauss–Legendre quadrature of Bernstein polynomials, we get

pn+1 = b0 +

n  ωk k =0

= b0 +

n  ωk k =0

= b0 +

2

2n 

2

p

2n  j =0

= b2n+1 .

1 + τk

⎣(2n + 1)

 (2n + 1)

2n 

 B 2n j

j =0 n  ωk k =0

b j



2

⎡

j =0

= b0 +



2

 B 2n j

1 + τk 2

1 + τk 2



⎤ b j ⎦

 b j

S.H. Kim, H.P. Moon / Computer Aided Geometric Design 54 (2017) 1–14

5

Next, we show that the length of the polygon realizes the arc length of the PH curve. Since the speed p (t ) of the PH curve is a polynomial of order 2n, we can apply the Gauss–Legendre quadrature again to get n 

   1 n  ωk   1 + τk   p  = p (t ) dt . 2  2

pk  =

k =0

k =0

(3)

0

2

This completes the proof.

Because of the symmetry of the problem, we can construct the rectifying polygon of the PH curve p backward from the end point in such a way that

pn+1 = b2n+1 ,

ωk

pk = pk+1 −

2

p



1 + τk 2



(4)

for k = n, n − 1, · · · , 0.

,

One can easily check that this backward scheme produces the same rectifying polygon. Both Equation (2) and Equation (4) construct the vertices of a rectifying polygon recursively. If we perform this recursive construction, each vertex pk can be explicitly written in terms of the Bézier control points of the PH curve. We here claim that the vertices of the rectifying polygon are obtained by convex combinations of the Bézier control points b j . Theorem 4. For a PH curve of degree 2n + 1 with the Bézier control points b0 , · · · , b2n+1 , all the vertices p0 , · · · , pn+1 of the rectifying polygon (2) are convex combinations of b0 , · · · , b2n+1 . Proof. By expressing p in Equation (2) in terms of the forward differences of the Bézier control points, we get

pk = (2n + 1)

2n ωk 

2

 B 2n j

1 + τk

 b j .

2

j =0

For k = 0, · · · , n, each

pk+1 = b0 +

k 

pl

l =0

is a barycentric combination of the points b j since the sum of coefficients is one. So it suffices to show that all the coefficients are nonnegative. If we set

μl, j = (2n + 1)

ωl 2



B 2n j

1 + τl 2

 ,

for l = 0, · · · , n and j = 0, · · · , 2n, then all n 

μl, j ’s are positive and

1 B 2n j (t ) dt = 1

μl, j = (2n + 1)

l =0

(5)

0

for j = 0, · · · , 2n. The vertex pk+1 of the rectifying polygon can then be written as

pk+1 = b0 +

2n k  

 = 1−

l =0 j =0 k 



μl, j b j+1 − b j

μl,0





  k  k  2n    b0 + (μl, j −1 − μl, j ) b j + μl,2n b2n+1 .

l =0

j =1

l =0

From Equation (5), it is clear that the first coefficient 1 − denote the intermediate coefficients as

mk, j =

l =0

k

l=0

μl,0 and the last coefficient

k  (μl, j −1 − μl, j ) l =0

for k = 0, · · · , n and j = 1, · · · , 2n. Then m0, j = 0, and mn, j = 0 by Equation (5).

k

l=0

μl,2n are positive. Let us

6

S.H. Kim, H.P. Moon / Computer Aided Geometric Design 54 (2017) 1–14

We now fix the index j and consider the sequence {mk, j }nk=0 . This sequence is the cumulative summation of (μl, j −1 −

μl, j ), so its increment is

mk−1, j = mk, j − mk−1, j = μk, j −1 − μk, j      ωk 2n 1 + τk 1 + τk − B 2n = (2n + 1) B j −1 j 2

2

2



(2n)! 1 + τk = (2n + 1) 1− 2 (2n − j + 1)! j ! 2 2

ωk

2n− j 

1 + τk

 j −1 

j 2n + 1

2

−

1 + τk 2

 .

This increment mk−1, j is positive while (1 + τk )/2 < j /(2n + 1), and becomes negative after that. So the sequence {mk, j }nk=0 starting from m0, j = 0 increases while (1 + τk )/2 < j /(2n + 1), then it decreases to mn, j = 0. Therefore all the intermediate coefficients mn, j ’s are nonnegative. 2 We now provide a few illustrative examples for the rectifying polygon of low degree PH curves. The above theorem is valid for n ≥ 0. When n = 0, the above theorem states that a PH curve of degree 1, which is a straight line, has a rectifying polygon [p0 p1 ] = [b0 b1 ]. The first nontrivial case occurs when n = 1. Let us apply the above theorem to the case n = 1. The cubic PH curve

p(t ) =

3 

B 3j (t )b j

j =0

has a rectifying polygon [p0 p1 p2 ] which is defined by p0 = b0 , p2 = b3 and

p1 = b0 +

= b0 +

ω0 2





p

3ω0

1 + τ0



2 2 

2

 B 2j

1 + τ0

 b j .

2

j =0

Since the Gauss–Legendre quadrature of order 2 has the nodes middle point p1 can be explicitly written as

p1 =

2−

4

√

√

√

3

√

τ0 = −1/ 3, τ1 = 1/ 3 and the weights ω0 = ω1 = 1, the

(b0 + b3 ) +

3

4

(b1 + b2 ).

(6)

The following example illustrates a cubic PH curve and its rectifying polygon. In order to show the exact arc-length realization of rectifying polygon, we express all real numbers as algebraic numbers and perform all computations using algebraic arithmetics without any floating point approximation. Example 5. A cubic Bézier curve p(t ) = points satisfy the condition b0 b2 = condition as follows:

b0 = 0,

b1 =

2 3

+

1 2

i,

b2 =

3

3 j =0 B j (t )b j with complex control points b j becomes a PH curve 2 (b1 ) . We choose the control points of a cubic Bézier curve which

√

4+5 3 12

+

1 4

i,

if its control abide by this

b3 = 1.

This control polygon is shown as a blue polygon in Fig. 1. The corresponding cubic PH curve is illustrated as a blue curve. We now compute the rectifying polygon [p0 p1 p2 ] of this PH curve. Both end points of the rectifying polygon should agree with the end points of the PH curve, i.e., p0 = b0 = 0, p2 = b3 = 1. The middle point p1 computed by using Equation (6) is

p1 =

13 16

√

+

3 3 16

i.

The corresponding rectifying polygon is shown as the yellow polygon in Fig. 1. Let us check that the rectifying polygon realizes the arc length of the PH curve p. The arc-length L of p can be obtained by the Gauss–Legendre quadrature as

1 L= 0

p (t )dt

S.H. Kim, H.P. Moon / Computer Aided Geometric Design 54 (2017) 1–14

7

Fig. 1. Rectifying polygon for a cubic PH curve. (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

 

=

ω0 

 p

 



1 + τ0  

+

ω1 

 p



1 + τ1  

 2 2 2  2   2     √ √   1   1 3   2 3− 3    2 3+ Bk Bk bk  + 3 bk  = 3  2  2 6 6 2

k =0

k =0

5

= . 4

On the other hand, the length of the rectifying polygon is

5

p0  + p1  = . 4

So the lengths of the PH curve and its rectifying polygon are exactly the same. We now construct the rectifying polygon of a quintic PH curve using Theorem 3 with n = 2. Let p be a quintic PH curve expressed as a Bézier curve

p(t ) =

5 

B 5j (t )b j .

j =0

The Bézier control points of the quintic PH curve p cannot be chosen arbitrary. One can find the sufficient and necessary condition for the control points to make p a PH curve in Farouki (2007). The rectifying polygon of the quintic PH curve has 3 legs, and can be expressed as [p0 p1 p2 p3 ]. Both end points should be p0 = b0 and p3 = b5 . To get symmetric expression for two middle points p1 and p2 , we apply Equation (2) to p1 , and Equation (4) to p2 , respectively. Then we get

p1 = b0 +

p2 = b5 −

4 5ω0 

2



B 4j

j =0

4 5ω2 

2

j =0

√



b j ,

2

 B 4j

1 + τ0

1 + τ2

 b j

2

√

where τ0 = − 3/5, τ2 = 3/5, ω0 = 5/9, above equations can be written as

√

p1 =

41 − 8 15 72

+ p2 =

−10 + 4 15 72

√

72

+

√

b0 +

15 + 4 15 72

√

31 − 8 15

√

b3 +

b0 +

10 + 4 15 72

ω2 = 5/9. By substituting the values of nodes τ0 , τ2 and weights ω0 , ω2 , the 10 + 4 15

−15 + 4 15 72

√

72

72

√

−15 + 4 15

b3 +

√

b1 +

√

15 + 4 15 72

b4 +

b1 + b4 +

b2

√

31 − 8 15 72

√

−10 + 4 15 72

√

41 − 8 15 72

b5 , (7)

b2

b5 .

As shown in Theorem 4, these two middle points p1 and p2 are convex combinations of the Bézier control points b0 , · · · , b5 .

8

S.H. Kim, H.P. Moon / Computer Aided Geometric Design 54 (2017) 1–14

Fig. 2. Rectifying polygon of a quintic PH curve. (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

Example 6. We prepare a quintic PH curve using Hermite interpolation method (Farouki and Neff, 1995). For a quintic Bézier curve

p(t ) =

5 

B 5j (t )b j ,

j =0

we set the end points as b0 = 0 and b5 = 1. To specify the derivatives at both end points. We choose two more control points b1 and b4 as

b1 = 0.1 + 0.25 i,

b4 = 0.9 + 0.4 i.

Then the intermediate control points b2 and b3 can be determined by the PH condition. This Hermite interpolation problem has 4 different solutions. Several algorithms for the selection of the “best” quintic PH curve have been proposed (Choi et al., 2007; Choi and Kwon, 2008; Moon et al., 2001). For “reasonable” Hermite data, most of the selection schemes suggest the same best solution. For the above Hermite data, the best quintic PH curve has the intermediate control points

b2 = 0.29102072 + 0.46044857 i,

b3 = 0.57948601 + 0.54478958 i,

which are computed numerically up to 8 significant digits. Fig. 2 illustrates this quintic PH curve and its Bézier control polygon as blue polygon. The rectifying polygon [p0 p1 p2 p3 ] for this quintic PH curve can be computed by Equation (7) as

p0 = 0, p1 = 0.19596175 + 0.31318651 i, p2 = 0.74831322 + 0.39911457 i,

(8)

p3 = 1. This rectifying polygon is shown as a yellow polygon in Fig. 2. Remark 7. The rectifying polygon constructed by either Equation (2) or Equation (4) has the end point agreement property such as

p0 = b0 ,

pn = b2n+1 .

In fact, any Bézier curve has this property. The proof of the end point agreement property in Theorem 3 does not rely on the PH condition of p. So if we construct a polygon [p0 p1 p2 p3 ] using Equation (2) or Equation (4), from an arbitrary Bézier curve, then we still get the end point agreement. However, this polygon cannot rectify a non-PH Bézier curve. 4. PH curves with given rectifying polygon The rectifying polygon [p0 · · · pn+1 ] has the same complex degrees of freedom as the PH curve of degree 2n + 1. So we suggest using the rectifying polygon as a kind of control polygon for PH curves. We call this polygon a rectifying control polygon of PH curves. We discovered that the rectifying control polygons have many advantages over Bézier control polygons to guide the shape of PH curves. We here want to solve the inverse problem of Theorem 3: For a given polygon [p0 · · · pn+1 ], find PH curves which have [p0 · · · pn+1 ] as their rectifying polygon. To answer this question, we need to exploit the Bernstein–Vandermonde matrix.

S.H. Kim, H.P. Moon / Computer Aided Geometric Design 54 (2017) 1–14

9

For a set of distinct parameters t 0 , t 1 , · · · , tn , the Vandermonde matrix is defined by

⎛

1 t0 ⎜ 1 t1 ⎜

Vn = ⎜ . .

.. .

⎝.

1 tn

⎞ · · · t n0 n ⎟ · · · t1 ⎟ . ⎟. .. . .. ⎠ · · · tnn

Each row of V n is the values of the power basis {1, t , t 2 , · · · , t n } at each parameter t i . The Bernstein–Vandermonte matrix M n is defined by using the Bernstein basis instead of the power basis such as

⎛

⎞

··· ··· Mn = ⎜ .. .. .. ⎝ . . . B n0 (tn ) B n1 (tn ) · · · B n0 (t 0 ) ⎜ B n0 (t 1 ) ⎜

B n1 (t 0 ) B n1 (t 1 )

B nn (t 0 ) B nn (t 1 ) ⎟ ⎟

.. ⎟ . . ⎠ n B n (tn )

It is well known that V n has the determinant

#

det ( V n ) =

(t j − t i ),

0≤i < j ≤n

which is nonzero for distinct parameters t 0 , · · · , tn . The determinant of the Bernstein–Vandermonde matrix is given by

  

det ( M n ) =

n

n

0

1

 

···

n n

det ( V n ).

Thus the Bernstein–Vandermonde matrix is also nonsingular if the parameters t 0 , · · · , tn are distinct. One may consult Marco and Martìnez (2007) for more information on the Bernstein–Vandermonde matrix. For a given rectifying control polygon, there are multiple instances of PH curves that share the same rectifying polygon. The following theorem presents the method to identify these PH curves. Theorem 8. For a given polygon [p0 p1 · · · pn+1 ] of n + 1 segments, there exist 2n PH curves of degree 2n + 1 whose rectifying polygon is [p0 p1 · · · pn+1 ]. These PH curves can be specified by

p(0) =p0 , 

p (t ) =

 n 

2

(9)

B ln (t )zl

l =0

where z0 , z1 , · · · , zn are the solution of the linear system

⎞ ⎛  ± ω20 p0 z0 ⎟ ⎜  ⎟ ⎜ z1 ⎟ ⎜ ± ω21 p1 ⎟ ⎜ ⎟ ⎜ ⎟ Mn ⎜ . ⎟ = ⎜ ⎟ .. ⎝ .. ⎠ ⎜ ⎟ ⎜ . ⎠ ⎝  zn 2 ± ωn pn ⎛

⎞

(10)

with the Bernstein–Vandermonde matrix

⎛

⎜ ⎜ ⎜ Mn = ⎜ ⎜ ⎜ ⎝





1 +τ 0 2

B n0 1+2τ1

B n0

B n0

.. .

1+τn 2







1 +τ 0 2

B n1 1+2τ1

B n1

B n1

.. .

1+τn 2



··· ··· .. . ···



⎞

.. .

⎟ ⎟ ⎟ ⎟. ⎟ ⎟

⎠

1 +τ 0 2

B nn 1+2τ1

B nn

B nn

1+τn 2

Proof. Let p be a PH curve of degree 2n + 1. Then its hodograph p is a square of a complex-valued polynomial of degree n, which can be expressed as 

p (t ) =

 n 

2 B ln (t )zl

.

l =0

If the rectifying polygon of p is [p0 p1 · · · pn+1 ], then Equation (2) yields

10

S.H. Kim, H.P. Moon / Computer Aided Geometric Design 54 (2017) 1–14

ωk

pk =

2





p





1 + τk

=

2



n ωk 

2

B ln

1 + τk

 2 zl

2

l =0

for k = 0, · · · , n. Each of these equations can be expressed as



n 

B ln

1 + τk

2

zk = ±

2

l =0

$



ωk

pk .

These equations constitute the linear system (10). Since the signs of the right hand side of Equation (10) can be selected arbitrary, we have 2n+1 linear systems. A pair of linear systems with opposite sign combination has the solutions z0 , z1 , · · · , zn with opposite sign combination. This pair of solutions produces the same hodograph p by Equation (9). So the total number of PH curves with the given rectifying polygon is 2n . 2 Remark 9. Because of the nonlinear nature of PH curves, the rectifying control polygon does not inherit some of the nice properties of Bézier control polygons. The rectifying control polygon does not have end tangent interpolation property. However, the computation of the end tangent from the rectifying control polygon is not a complicated process. One can easily see from (9) that the end tangents of the PH curve in Theorem 8 are

p (0) = z20 ,

p (1) = zn2

where z0 , zn are the solution of (10). We here illustrate a few examples of PH curves specified by the given rectifying control polygons. The first example is the inverse problem of Example 5. Example 10. Suppose a rectifying control polygon [p0 p1 p2 ] is given by

p0 = 0,

p1 =

13 16

√

3 3

+

16

p2 = 1,

i,

which is the rectifying polygon of the PH curve in Example 5. The above theorem claims that there are two cubic PH curves whose rectifying polygon is [p0 p1 p2 ]. Equation (10) can be written as

1



6

√

√ 

3 + √3 3 − √3 3− 3 3+ 3

whose solutions are



z0 z1

√ 



=

2

2

z0 z1



=

 √  ±√2p0 , ± 2p1

√ √ √ √  ±(1 + √3) p0 ± (1 − √3) p1 √ √ . ±(1 − 3) p0 ± (1 + 3) p1

So we can compute the Bézier√control points √ b0 , · · · , b3 from Equation (9) after we fix the signs in the above expression. If we choose the same sign for p0 and p1 , then the corresponding Bézier control points are

b0 = 0,

b1 =

2 3

+

1 2

√

i,

b2 =

4+5 3 12

+

1 4

i,

b3 = 1.

These points of the original cubic PH curve in Example 5. On the other hand, if we choose different signs √ are the control √ for p0 and p1 , we get

√

b0 = 0,

b1 =

8+5 3 12

+

1 4

i,

b2 =

1 3

+

1 2

i,

b3 = 1.

Fig. 3 (a) illustrates this cubic PH curve, and Fig. 3 (b) shows this curve together with the original PH curve in Example 5. The next example illustrates quintic PH curves with the given rectifying control polygon. This is the inverse problem of Example 6. Example 11. We want to find quintic PH curves for the rectifying control polygon [p0 p1 p2 p3 ] given by Equation (8) in Example 6. There exist 4 quintic PH curves which share the rectifying control polygon [p0 p1 p2 p3 ]. To compute these curves, we need to solve the linear systems

⎛  ⎞  p ± 18 0 8 + 15 4 8 − 15 z0 ⎜ 8 ⎟ 1 ⎜ ⎟ ⎝ 5√ 10 5√ ⎠ ⎝ z1 ⎠ = ⎜ ± 18 p1 ⎟ , 5 ⎝ ⎠ 20  z2 8 − 15 4 5 + 15 ± 18 p2 8 ⎛

√

√

⎞⎛

⎞

(11)

S.H. Kim, H.P. Moon / Computer Aided Geometric Design 54 (2017) 1–14

Fig. 3. (a) Cubic PH curve defined by the given rectifying control polygon with different sign combination for cubic PH curves with the same rectifying control polygon.

11

√ √ p0 and p1 . (b) Comparison of two

Table 2 Intermediate control points of the quintic PH curves with the same rectifying control polygon for different sign configurations. (0,0)

= 0.1 + 0.25 i

b1

(0,0)

= 0.29102072 + 0.46044857 i

b2

(0,0)

= 0.57948601 + 0.54478958 i

b3

(0,0)

= 0.9 + 0.4 i

b4

(1,0)

= 0.61756386 + 0.81100984 i

b1

(1,0)

= −0.29947217 − 0.25849981 i

b2

(1,0)

= 0.42656732 + 0.77347242 i

b3

b1 b2 b3 b4 b1 b2 b3

(1,0)

b4

which has the solution

= 1.12584773 + 0.32925569 i

(0,1)

= −0.10072651 + 0.18878062 i

(0,1)

= 0.36309317 + 0.65645730 i

(0,1)

= 1.42348283 − 0.06793782 i

(0,1)

= 0.18465724 + 0.87793804 i

(1,1)

= 1.08486144 + 0.75096640 i

(1,1)

= −1.22943586 − 0.06425497 i

(1,1)

= 2.27260028 + 0.16250892 i

(1,1)

= −0.25751913 + 0.80601780 i

b4

 √ ⎞ √ √ 2 p1 ± (5 − 15) p2 5 z0 ⎜ ⎟  √ √ √ ⎟ ⎝ z1 ⎠ = 1 ⎜ ⎜ ⎟. ∓5 p0 ± 32 25 p1 ∓ 5 p2 ⎠ 4⎝  √ √ √ √ z2 2√ ±(5 − 15) p0 ∓ 8 5 p1 ± (5 + 15) p2 ⎛

⎛

⎞

±(5 +

√

√

15) p0 ∓ 8

The above solution has 3 multiple signs, which produces 8 different sign configurations. Since the√opposite sign configuration produces the same PH curve, √ we can√freely fix one of the multiple signs. We fix the sign of p1 in Equation (11) as plus. The choice of signs for p0 and p2 are represented by (−1)σ0 and (−1)σ2 respectively, where σi ∈ {0, 1} for i = 0, 2. Then the above solution can be rewritten as

 √ % √ % √ % (−1)σ0 (25 + 5 15) p0 − 8 10 p1 + (−1)σ2 (25 − 5 15) p2 , 20  √ % % % 1  (σ0 ,σ2 ) −(−1)σ0 25 p0 + 32 10 p1 − (−1)σ2 25 p2 , z1 = 20  √ % √ % √ % 1  (σ1 ,σ2 ) z2 = (−1)σ0 (25 − 5 15) p0 − 8 10 p1 + (−1)σ2 (25 + 5 15) p2 . (σ0 ,σ2 )

z0

=

1 

20

We can compute the control points of the quintic PH curves using Equation (9), or equivalently

b0 = p0 , b1 = b0 + b2 = b1 +

b5 = p3 , 1 5 1 5

z20 ,

b4 = b5 −

z0 z1 ,

b3 = b4 −

1 5 1 5

z22 , z1 z2 .

Table 2 exhibits the results of numerical computation for each sign configuration (σ0 , σ2 ). We here skip both end points since they are fixed as b0 = 0 and b5 = 1. For the first sign configuration (σ0 , σ2 ) = (0, 0), the above process recovers the original quintic PH curve in Example 6. Although this process involves numerical errors of the floating point arithmetic, it recovers the correct values for b2 and b3 up to 8 significant digits. The other sign configurations produce different PH curves. The PH curve with sign configuration

12

S.H. Kim, H.P. Moon / Computer Aided Geometric Design 54 (2017) 1–14

Fig. 4. (a), (b), (c) The quintic PH curves for the given rectifying control polygon with the sign configuration (0, 1), (1, 0), and (1, 1), respectively. (d) All 4 quintic PH curves with the same rectifying control polygon.

Fig. 5. (a) Variation of the cubic PH curve as p1 moves to the right. (b) Variation of the quintic PH curve as p2 moves upward.

(0, 1), (1, 0), and (1, 1) are shown in Fig. 4 (a), (b), and (c), respectively. We put all these 4 quintic PH curves on the same plane in Fig. 4 (d). The rectifying control polygon developed so far is well-suited to control a PH curve, because it has the same degrees of freedom as the corresponding PH curve. Thus an individual point of the rectifying control polygon of a PH curve can be moved arbitrarily for the design purposes. The continuous movement of a rectifying control point results in the continuous deformation of corresponding PH curve. Fig. 5 illustrates this situation. Starting from the cubic PH curve with the same sign √ √ for p0 and p1 in Example 10, we move the middle control point p1 to the right to get the sequence of cubic PH curves shown in Fig. 5 (a). Similarly, Fig. 5 (b) shows the sequence of quintic PH curves produced by upward movement of the control point p2 in Example 11 with the same signs in (11).

S.H. Kim, H.P. Moon / Computer Aided Geometric Design 54 (2017) 1–14

Fig. 6. The quintic PH curves and the Bézier control polygons for given rectifying control polygons (a) p0 = 0, p1 = (b) p0 = 0, p1 =

4 5

+

4 5

i, p2 =

1 5

−

4 5

i, p3 = 1.

13

2 3

+

3 5

i, p2 =

1 3

+

1 10

i, p3 = 1 and

The rectifying control polygon for a PH curve has some drawbacks compared with the Bézier control polygon for a polynomial curve. First, a single rectifying control polygon produces multiple instances of PH curves. However this is a typical situation for a PH curve construction problems, and several selection methods (Choi et al., 2007; Farouki and Neff, 1995) for the best solution have been proposed. Another drawback of the rectifying control polygon is that it does not have the shape preserving properties (Carnicer, 1999; Carnicer and Peña, 1993; Peña, 1999). As one can see in the previous examples, PH curves lie outside of the convex hull of the rectifying control polygon. So it does not have the convex hull property. The rectifying control polygon [p0 p1 p2 p3 ] in Example 11 has monotonicity along the positive x axis. It also rotates along the clockwise orientation. The corresponding PH curve with the positive signs in Equation (11) has the same monotonicity and the same orientation. However, other instances of PH curves have tiny loops that change the monotonicity and the orientations. So the rectifying control polygon has neither the monotonicity preserving property nor the orientation preserving property in general. The only shape preserving property of the rectifying control polygon is the length diminishing property, which means the length of the curve is less than or equal to the length of the control polygon. The length diminishing property of Bézier polygon sometimes makes it difficult to predict the shape of the curve from the shape of the control polygon. It is because the length of the Bézier polygon is too long compared with the length of the curve and the Bézier control polygon tends to exaggerate the behavior of the curve. On the other hand, the rectifying control polygon has the length preserving property such that the length of the control polygon is the same as the length of the curve. Fig. 6 illustrates the comparison between the Bézier control polygons and the rectifying control polygons of quintic PH curves. 5. Conclusion We proposed to use the rectifying control polygon instead of the Bézier polygon to guide PH curves. The rectifying control polygon has exactly the same degrees of freedom as the corresponding PH curve. It also has the same arc length as the PH curve. So we expect that the rectifying control polygon can be used to develop simpler PH curve manipulation algorithms than using the standard Bézier polygon. Acknowledgement This work was supported by the National Research Foundation of Korea (NRF) Grant funded by the Korea Government (MEST) (No. 20120001978). References Atkinson, K.E., 1989. An Introduction to Numerical Analysis, 2nd edition. John Wiley & Sons, Hoboken, NJ.

14

S.H. Kim, H.P. Moon / Computer Aided Geometric Design 54 (2017) 1–14

Atkinson, K.E., Han, W., 2004. Elementary Numerical Analysis, 3rd edition. John Wiley & Sons, Inc. Carnicer, J.M., 1999. Interpolation, shape control and shape properties. In: Shape Preserving Representations in Computer-Aided Geometric Design, pp. 15–43. Carnicer, J.M., Peña, J.M., 1993. Shape preserving representations and optimality of the Bernstein basis. Adv. Comput. Math. 1 (2), 173–196. Choi, H.I., Farouki, R.T., Kwon, S.H., Moon, H.P., 2007. Topological criterion for selection of quintic Pythagorean-hodograph Hermite interpolants. Comput. Aided Geom. Des. 25, 411–433. Choi, H.I., Kwon, S.H., 2008. Absolute hodograph winding number and planar PH quintic splines. Comput. Aided Geom. Des. 25, 230–246. Farouki, R.T., 1994. The conformal map z → z2 of the hodograph plane. Comput. Aided Geom. Des. 11 (4), 363–390. Farouki, R.T., 2007. Pythagorean-Hodograph Curves: Algebra and Geometry Inseparable. Springer. Farouki, R.T., 2016. Construction of G 1 planar Hermite interpolants with prescribed arc lengths. Comput. Aided Geom. Des. 46, 64–75. Farouki, R.T., Giannelli, C., Sestini, A., 2015. Identification and “reverse engineering” of Pythagorean-hodograph curves. Comput. Aided Geom. Des. 34, 21–36. Farouki, R.T., Giannelli, C., Sestini, A., 2016. Local modification of Pythagorean-hodograph quintic spline curves using the B-spline form. Adv. Comput. Math. 42 (1), 199–225. Farouki, R.T., Neff, C.A., 1995. Hermite interpolation by Pythagorean hodograph quintics. Math. Comput. 64 (212), 1589–1609. Farouki, R.T., Sakkalis, T., 1990. Pythagorean hodographs. IBM J. Res. Dev. 34, 736–752. Huard, M., Farouki, R.T., Sprynski, N., Biard, L., 2014. C 2 interpolation of spatial data subject to arc-length constraints using Pythagorean–hodograph quintic splines. Graph. Models 76 (1), 30–42. Kythe, P.K., Schäferkotter, M.R., 2005. Handbook of Computational Methods for Integration. Chapman & Hall/CRC. Marco, A., Martìnez, J.-J., 2007. A fast and accurate algorithm for solving Bernstein–Vandermonde linear systems. Linear Algebra Appl. 422, 616–628. Moon, H.P., Farouki, R.T., Choi, H.I., 2001. Construction and shape analysis of PH quintic Hermite interpolants. Comput. Aided Geom. Des. 18 (2), 93–115. Peña, J.M., 1999. Bases with optimal shape preserving properties. In: Shape Preserving Representations in Computer-Aided Geometric Design. Roulier, J.A., 1993. Specifying the arc length of Bézier curves. Comput. Aided Geom. Des. 10 (1), 25–56. Roulier, J.A., Piper, B., 1996. Prescribing the length of parametric curves. Comput. Aided Geom. Des. 13 (1), 3–22. Wang, G., Fang, L., 2009. On control polygons of quartic Pythagorean-hodograph curves. Comput. Aided Geom. Des. 26, 1006–1015.