C2 Hermite interpolation by Pythagorean-hodograph quintic triarcs

Computer Aided Geometric Design 31 (2014) 412–426

Contents lists available at ScienceDirect

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

C 2 Hermite interpolation by Pythagorean-hodograph quintic triarcs Bohumír Bastl a,b , Michal Bizzarri b,a , Karla Ferjanˇciˇc d , Boštjan Kovaˇc e , Marjeta Krajnc c , Miroslav Láviˇcka a,b , Kristýna Michálková b,a,∗ , Zbynˇek Šír a , Emil Žagar c a

University of West Bohemia, Faculty of Applied Sciences, Department of Mathematics, Univerzitní 8, Plzen, ˇ Czech Republic University of West Bohemia, NTIS, Faculty of Applied Sciences, Univerzitní 8, Plzen, ˇ Czech Republic University of Ljubljana, IMFM and FMF, Jadranska 19, Ljubljana, Slovenia d University of Primorska, IAM, Muzejski trg 2, Koper, Slovenia e Abelium d.o.o., Kajuhova 90, Ljubljana, Slovenia b c

a r t i c l e

i n f o

Article history: Available online 6 September 2014 Keywords: Pythagorean-hodograph curves Hermite interpolation Triarc PH quintic

a b s t r a c t In this paper, the problem of C 2 Hermite interpolation by triarcs composed of Pythagoreanhodograph (PH) quintics is considered. The main idea is to join three arcs of PH quintics at two unknown points – the first curve interpolates given C 2 Hermite data at one side, the third one interpolates the same type of given data at the other side and the middle arc is joined together with C 2 continuity to the first and the third arc. For any set of C 2 planar boundary data (two points with associated first and second derivatives) we construct four possible interpolants. The best possible approximation order is 4. Analogously, for a set of C 2 spatial boundary data we find a six-dimensional family of interpolating quintic PH triarcs. The results are confirmed by several examples. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Shapes (curves, surfaces, volumes) in Computer Aided Geometric Design (CAGD), and in a vast variety of subsequent applications, are often described by piecewise polynomial/rational representations. However, not every shape can be represented using polynomial/rational functions, see Shafarevich (1974) for more details. Another problem in CAGD is that many natural geometric operations, such as offsetting, do not preserve rationality of derived objects. Nonetheless, offsets to certain special classes of shapes admit exact rational representations. In the case of planar curves, the class of Pythagorean-hodograph (PH) curves as polynomial curves possessing rational offset curves and (piecewise) polynomial arc-length functions was introduced in Farouki and Sakkalis (1990). The concept of planar PH curves was later generalized to polynomial spatial PH curves in Farouki and Sakkalis (1994). These curves have the following attractive properties: the arc length of any segment can be determined exactly without numerical approximation and any canal surface based on spatial PH curves as its spine and a rational radius function can be easily parameterized

*

ˇ Czech Republic. Corresponding author at: University of West Bohemia, Faculty of Applied Sciences, Department of Mathematics, Univerzitní 8, Plzen, E-mail addresses: [email protected] (B. Bastl), [email protected] (M. Bizzarri), [email protected] (K. Ferjanˇciˇc), [email protected] (B. Kovaˇc), [email protected] (M. Krajnc), [email protected] (M. Láviˇcka), [email protected] (K. Michálková), [email protected] (Z. Šír), [email protected] (E. Žagar). http://dx.doi.org/10.1016/j.cagd.2014.08.002 0167-8396/© 2014 Elsevier B.V. All rights reserved.

B. Bastl et al. / Computer Aided Geometric Design 31 (2014) 412–426

413

rationally. The case of PH curves in Rn for n > 3 is still not solved satisfactorily. Exploiting recent results from number theory, the structure of PH curves in dimensions n = 5 and n = 9 was characterized in Sakkalis and Farouki (2012). A detailed survey of the literature on shapes with a Pythagorean property and their applications in technical practice can be found in Farouki (2008). Many PH interpolation techniques (yielding C 1 /G 1 or C 2 /G 2 continuity) which form the cornerstone of subsequent approximation algorithms devoted to particular problems originating in technical practice have been formulated in recent years. These methods by planar or spatial polynomial PH curves are usually based on low degree polynomials – mainly 3 and 5 (in some cases also 7 and 9), as only odd degree PH curves are regular. Firstly, Farouki and Neff focused on C 1 Hermite interpolation by PH quintics (cf. Farouki and Neff, 1995). The thorough analysis of other interpolation schemes followed. For planar cubic PH curves one of the first interpolation methods was given in Meek and Walton (1997) where G 1 interpolation of Hermite data was analyzed. These results were generalized in Jakliˇc et al. (2010) to G 2 interpolation by ˇ the same objects. Later the problem was revisited in Byrtus and Bastl (2010), Cernohorská and Šír (2010). For quintic planar PH curves, several results on first and second order continuous spline interpolation are given in Pelosi et al. (2007), Moon et al. (2001), Šír et al. (2006), Albrecht and Farouki (1996), Choi et al. (2008), Farouki (2008). For spatial curves, G 1 Hermite interpolation by PH cubics was thoroughly investigated in Jüttler and Mäurer (1999). Those results were later generalized to some level in Pelosi et al. (2005). The most general results on this type of interpolation can be found in Kwon (2010) and Jakliˇc et al. (2012). The problem of C 1 and C 2 Hermite interpolation by spatial PH curves of degree ≥ 5 has been studied in Šír and Jüttler (2005, 2007), Farouki et al. (2003, 2008), Sestini et al. (2013). G 2 [C 1 ] Hermite interpolation with PH curves of degree 7 has been analyzed in Jüttler (2001). C 1 Hermite interpolation by cubic polynomial spline curves is always uniquely solvable, so it is obvious that the same problem cannot be solved by PH cubic splines. However, cubic (i.e., low degree) polynomial interpolating splines are preferred in many applications. An interesting approach was used in Farouki and Peters (1996), where C 1 Hermite interpolation via double-Tschirnhausen cubics (TC-biarcs) has been considered. The main idea was to join two arcs of planar cubic PH curves at some unknown point – the first curve interpolates C 1 Hermite data at one side, the other one interpolates the same type of data at the other side and the arcs are joined together with C 1 continuity. This approach was recently improved in Bastl et al. (2013), where planar uniform and non-uniform cubic PH biarcs were studied and applied. Using a close analogy between the complex representation of planar PH curves and the quaternion representation of spatial PH curves, this idea was later transformed to the spatial case in Bastl et al. (2014) where the Hermite C 1 interpolation scheme by spatial cubic PH biarcs was presented and a general algorithm for computing interpolants was designed, studied and applied. In this paper, we will extend the ideas from Bastl et al. (2013, 2014) to the C 2 Hermite interpolation by planar/spatial PH curves. Instead of using PH curves of degree 9 as in Šír and Jüttler (2007) we will consider PH triarcs, which enables us to drop the degree of the used PH curve to 5. In other words, the main idea is to join three arcs of PH quintics at two unknown points – the first curve interpolates given C 2 Hermite data at one side, the third one interpolates the same type of given data at the other side and the middle arc is joined together with C 2 continuity to the first and the third arc. The algorithm will be formulated firstly for planar PH quintic triarcs using complex representation. Then, a straightforward generalization to the quaternion representation yields an analogous algorithm for spatial PH quintic triarcs. The functionality of the designed technique will be demonstrated by several planar and spatial examples. The paper is organized as follows. In the next section some preliminaries of the used theory are given. The C 2 Hermite planar interpolation problem is presented and an algorithm for computing such interpolants is given in the third section. In the fourth section the presented approach is generalized to the spatial case and the numerical examples which confirm the theoretical results are shown. The main results of the paper are summarized in the concluding section. 2. Preliminaries In this section, we review fundamentals of the theory of planar and spatial Pythagorean-hodograph curves and their relation to complex numbers and quaternions, respectively. 2.1. Complex numbers and planar PH curves Planar Pythagorean-hodograph curves were introduced in Farouki and Sakkalis (1990). A polynomial parametric curve p = (x, y ) of degree n with the hodograph h = p  = (x , y  ) is a Pythagorean-hodograph curve (shortly PH curve) if the components of its hodograph h fulfill the condition

x 2 + y  2 = σ 2 , where σ is a polynomial. The main advantages of PH curves are that they possess rational offset curves and have (piecewise) polynomial arc-length function. With the help of Kubota’s theorem (see Farouki and Sakkalis, 1990; Kubota, 1972) hodographs of all planar polynomial PH curves can be expressed in the form

414

B. Bastl et al. / Computer Aided Geometric Design 31 (2014) 412–426





x = w u 2 − v 2 , y  = w (2uv ) for some real non-zero polynomials u, v, w, where w is usually taken as a constant. In many cases (and also in the follow-up sections), it is advantageous to use complex representation of points and curves in the plane to represent planar PH curves (see Farouki, 1994). More precisely, points in the plane are identified with complex numbers in the complex plane and planar polynomial parametric curves are identified with complex-valued polynomials. Then, p = x + i y is a PH curve iff there exists a complex polynomial w = u + i v, called the preimage, such that the hodograph h of p can be expressed in the form

h = w 2 = u 2 − v 2 + i 2uv . 2.2. Quaternions and spatial PH curves Similarly to the relation of planar PH curves and complex number/polynomials, spatial PH curves can be determined with the help of quaternion algebra. For a definition of quaternions and the description of standard operations in a quaternion space H, the reader is kindly referred, e.g., to Farouki (2008). Here, only some less known facts about quaternions, which will be used in Section 4, are reviewed. A standard multiplication on quaternions is not commutative, but it is possible to define a commutative multiplication. For a pair of quaternions A, B ∈ H, we define

A  B :=

1 2

(A i B¯ + B i A¯ ).

(1)

The result is always a pure quaternion which can be identified with a vector in R3 . A notation A2 := A  A will be used further on. The following lemmas review how linear and quadratic equations with respect to -operation can be solved (see Šír and Jüttler, 2007 and Farouki et al., 2002, respectively). Lemma 1. Let A be a given pure quaternion and B a given non-zero quaternion. Then all solutions of a linear -equation X  B = A form a one-parameter family

X (α ; B, A) :=

(α + A)Bi , BB¯

α ∈ R.

Lemma 2. Let A be a given non-zero pure quaternion. All solutions of a quadratic –equation X 2 = A form a one-parameter family

X := X (ϕ ; A) := X p (A)Qϕ ,

Qϕ := cos ϕ + i sin ϕ ,

ϕ ∈ [−π , π ),

(2)

where X p (A) is a particular solution given by

X p (A) :=

⎧ ⎪ ⎨ √A

A

 A  +i

 A A +i

⎪ ⎩ √A k,

,

A A

= −i,

A A

= −i.

Further, it can be proved that for arbitrary quaternions A and B it holds

AQϕ  BQψ = AQϕ −ψ  B = A  BQψ−ϕ .

(3)

Let us focus on spatial Pythagorean-hodograph curves now. Similarly to the planar case, a polynomial parametric curve p = (x, y , z) of degree n with the hodograph h = p  = (x , y  , z ) is a spatial Pythagorean-hodograph curve if the components of h fulfill the condition

x 2 + y  2 + z 2 = σ 2 , where σ is a polynomial. Analogously to the planar PH curves and their complex representation, it is advantageous to use quaternion algebra to represent spatial PH curves, where points in the space are identified with pure quaternions and spatial parametric polynomials are identified with pure quaternion polynomials. Then, p = x i + y j + z k is a spatial PH curve iff there exists a quaternion polynomial A = u + v i + p j + q k, called the preimage, such that the hodograph h can be expressed in the form

h = A i A¯ = A2 . Note that quaternions A and A Qϕ generate the same hodograph.

B. Bastl et al. / Computer Aided Geometric Design 31 (2014) 412–426

415

3. Planar interpolation problem and an algorithm This section is devoted to the study of C 2 Hermite interpolation with planar quintic PH triarcs. An algorithm for finding such a quintic PH triarc for given Hermite data is provided and the approximation order of the method is analyzed. 3.1. Algorithm In this section the interpolation problem by planar quintic PH triarcs is presented and an algorithm for computing the interpolants is given. Since the previous section shows that a planar PH curve is easily characterized by its complex preimage, vectors in R2 will be identified with complex numbers and vice versa. Suppose that A , B ∈ C are two given points and t A , t B ∈ C are given associated tangent vectors. Additionally, let us prescribe two second derivative vectors c A , c B ∈ C at A , B respectively. The goal is to find a C 2 continuous planar quintic PH triarc interpolant p : [τ0 , τ3 ] → C, composed of three planar quintic PH curves

p i : [τi −1 , τi ] → C, with

τ0 < τ1 < τ2 < τ3 , ⎧ ⎨ p 1 (t ), p (t ) = p 2 (t ), ⎩ p 3 (t ),

i = 1, 2, 3,

i.e.,

t ∈ [τ0 , τ1 ], t ∈ [τ1 , τ2 ],

(4)

t ∈ [τ2 , τ3 ],

satisfying

p  (τ0 ) = t A ,

p (τ0 ) = A ,

p  (τ0 ) = c A ,



p (τ3 ) = B ,

(5)



p (τ3 ) = t B ,

p (τ3 ) = c B .

(6)

Since p should be in C 2 ([τ0 , τ3 ]), we additionally require (k)

(k)

p i (τi ) = p i +1 (τi ),

i = 1 , 2 , k = 0, 1 , 2 .

(7)

The interval values τ0 < τ1 < τ2 < τ3 are fixed but can be chosen arbitrarily. Parametric curves p i , i = 1, 2, 3, are assumed to be PH curves, so they are characterized by three associated preimage curves

w i (t ) =

2 

w i , j B 2j

t − τ i −1

τ i −1

j =0

,

i = 1, 2, 3,

where τi −1 := τi − τi −1 , i = 1, 2, 3. Let p i be given in the Bernstein–Bézier basis as

p i (t ) =

5 

P i , j B 5j

j =0

t − τ i −1

τ i −1

i = 1, 2, 3.

,

By Farouki (2008), the control points P i , j are expressed by the control points of the preimages as

P i ,1 = P i ,0 + P i ,2 = P i ,1 + P i ,3 = P i ,2 + P i ,4 = P i ,3 + P i ,5 = P i ,4 +

τ i − 1 5

τ i − 1 5

w 2i ,0 , w i ,0 w i ,1 ,



τ i − 1 2 5

τ i − 1 5

τ i − 1 5

3

w 2i ,1 +

1 3

w i ,0 w i ,2 ,

w i ,1 w i ,2 , w 2i ,2 .

(8)

From (5), (6), (8) and some basic properties of Bézier curves we obtain

P 1,0 = A , and

w 21,0 = t A ,

w 1,0 w 1,1 − w 21,0 =

τ0 4

c A,

(9)

416

B. Bastl et al. / Computer Aided Geometric Design 31 (2014) 412–426

w 23,2 = t B ,

P 3,5 = B ,

w 23,2 − w 3,1 w 3,2 =

τ2 4

(10)

cB.

τ1 and τ2 clearly imply P 2,0 = P 1,5 and P 3,0 = P 2,5 . This, together with (8), gives

τ0 2 1 P 2,0 = A + w 21,0 + w 1,0 w 1,1 + w 21,1 + w 1,0 w 1,2 + w 1,1 w 1,2 + w 21,2 , 5 3 3

τ2 2 1 P 2,5 = B − w 23,0 + w 3,0 w 3,1 + w 23,1 + w 3,0 w 3,2 + w 3,1 w 3,2 + w 23,2 .

The continuity conditions at



5

3

Using (8) it follows that

P 2,5 − P 2,0 =

τ1

w 22,0 + w 2,0 w 2,1 +

5

3

2 3

w 22,1 +

1 3

(11) (12)

w 2,0 w 2,2 + w 2,1 w 2,2 + w 22,2 .

τ1 and τ2 further imply   1  2 1  w i ,2 − w i ,1 w i ,2 = w i +1,0 w i +1,1 − w 2i +1,0 , τ i −1 τ i

(13)

The first and the second order continuity conditions at

w 2i ,2 = w 2i +1,0 ,

i = 1, 2.

(14)

Eqs. (9)–(14) form a system of 11 complex equations for 11 complex unknowns w i , j , i = 1, 2, 3, j = 0, 1, 2, P 2,0 and P 2,5 . From Eqs. (9), (10) and (14) we derive

w 1,0 = χ A where each of

w 1,1 =

√

w 3,2 = χ B

t A,

√

w 1,2 = χ1 w 2,0 ,

tB,

w 3,0 = χ2 w 2,2 ,

(15)

χ A , χ B , χ1 and χ2 is either 1 or −1, and consequently tA +

τ0

√4

cA

χA t A

,

w 3,1 =

t B − 4τ2 c B

√

χB t B

(16)

.

From (14) it follows also

τ0 , τ1 1 w 2,1 + τ0 + τ1 τ1 τ1 τ2 1 w 2,1 + w 2,2 = τ1 + τ2 τ1 w 2,0 =



χ1 t A + 4τ0 c A , √ τ0 χ A tA

(17)

χ2 t B − 4τ2 c B . √ τ2 χ B tB

(18)

1



1

Here, we have assumed that w 2,0 = 0 and w 2,2 = 0, which should be fulfilled since we consider only regular quintic PH curves. Considering (15)–(18) in (11)–(13), we end up with one complex quadratic equation for w 2,1 , namely 3 

τ i −1

w 2i ,0

+ w i ,0 w i ,1 +

i =1

2 3

w 2i ,1

+

1 3

w i ,0 w i ,2 + w i ,1 w i ,2 +

w 2i ,2

− 5( B − A ) = 0.

(19)

It is easy to see that the expressions (8) for the control points of the triarc segments involve only the products χ1 χ A or χ2 χ B . Thus we can assume that, e.g., χ1 = χ2 = 1. We still have two possibilities for each of χ A or χ B , which together with two solutions of the quadratic equation for w 2,1 gives eight solutions in general. But multiplying the control points of the preimage by −1 does not change the curve, thus we can, e.g., choose χ A = 1. This leaves four possible choices for the ±,± solution. Let them be denoted by w 2,1 , where the first symbol ± denotes the sign of the χ B and the second one the sign in the solution of the quadratic equation (19). Consequently, we obtain four possible interpolants, which will be denoted accordingly as p ±,± . Example 1. Let us consider C 2 Hermite data

P 0 = (0, 0) ,

P2 =

3 2

 ,0

,

t 0 = (2, 2) ,

t 2 = 2,

2 3

 ,

c 0 = (−1, 1) ,

c 2 = (1, −1) ,

and let τi = 3i , i = 0, . . . , 3. With the help of the method described in the previous paragraphs we can find all four distinct quintic PH triarcs matching given data (see Fig. 1). For practical applications the scheme should automatically return only one interpolant, determined by some suitable criterion. Among four interpolants one can choose the triarc that provides the best approximation order and, as it is proved in the following subsection, only the p +,+ has the approximation order four, while for the other three triarcs, the approximation order falls. This justifies p +,+ as a reasonable choice. Now, two free parameters, namely τ1 and τ2 , remain and can be used to optimize the shape measures of the curve, see Section 3.3.

B. Bastl et al. / Computer Aided Geometric Design 31 (2014) 412–426

417

Fig. 1. Four distinct quintic PH triarcs for input data in Example 1 labelled by p +,+ (top left), p +,− (top right), p −,+ (bottom left) and p −,− (bottom right).

3.2. Asymptotic approximation order In this section the approximation properties of all four interpolating triarcs will be analyzed. Let the data be sampled from an analytic parametric curve f : [0, h] → C, s → f (s), parameterized by the arc-length and let ϕ : [τ0 , τ3 ] → [0, h], t → h (t − τ0 )/(τ3 − τ0 ), be a particular linear reparameterization. We can assume (without loss of generality) that τ0 = 0, τ3 = 1 and f (0) = 0, f  (0) = 1 to simplify the analysis. Using the well-known Frenet–Serret formulae, the curve can be expressed as

f (s) = s −

1 6

+

1 

1

κ02 s3 − κ0 κ1 s4 + 8

κ0 s 2 2

+

κ1 s 3

+

6

1  24

120



κ04 − 4κ2 κ0 − 3κ12 s5 

κ2 − κ03 s4 +

1  120





 

κ3 − 6κ02 κ1 s5 i + O s6 ,

(20)

κ where κ (s) = κ0 + κ1 s + κ22! s2 + 33! s3 + O (s4 ) is the Maclaurin series of the curvature of f . The task is to construct the triarc (4) that satisfies (5) and (6) for the data

A = ( f ◦ ϕ )(0) = 0, B = ( f ◦ ϕ )(1) = h −

t A = ( f ◦ ϕ ) (0) = h, 1 6

h3 κ02 −

1 8



h 4 κ0 κ1 +

c A = ( f ◦ ϕ ) (0) = h2 κ0 i, h 2 κ0 2

+

h 3 κ1

+

6

1 24

h4



κ2 − κ03





 

i + O h5 ,

   1  + h4 κ2 − κ03 i + O h5 , 2 2 2 6

    3 1 c B = ( f ◦ ϕ ) (1) = −h3 κ02 − h4 κ0 κ1 + h2 κ0 + h3 κ1 + h4 κ2 − κ03 i + O h5

t B = ( f ◦ ϕ ) (1) = h −

1

h3 κ02 −

1



h 4 κ0 κ1 + h 2 κ0 +

h 3 κ1

2

2

taken from the reparameterized curve f . From (15)–(18) it follows that the control points of the preimage expressed with the unknown w 2,1 expand as

w 1,0 =

√

h,

w 1,1 =

√

h+



1 4

h3/2 κ0 τ1 i,



τ1 w 2,1 √ τ1 h3/2 κ0 τ1 (τ1 − τ2 ) − w 1,2 = + h 1− i, τ2 τ2 4τ2 w 2,0 = w 1,2 ,

(τ2 − 1) w 2,1 w 2,2 = + τ1 − 1

√

h(τ1 − τ2 )χ B

+

h3/2 κ0 (τ1 − τ2 )(τ2 + 1)χ B

τ1 − 1 4(τ1 − 1) 2 5/2   h (κ0 − 2 κ1 i)(τ1 − τ2 )τ2 χ B − + O h7/2 , 8(τ1 − 1)

w 3,0 = w 2,2 ,

i

418

B. Bastl et al. / Computer Aided Geometric Design 31 (2014) 412–426

w 3,1 = w 3,2 =

√

1

√

4 1

hχ B + hχ B +

2

h3/2 κ0 (τ2 + 1)χ B i − h3/2 κ0 χ B i −

1 8

h5/2



1 8

h5/2









κ02 − 2 κ1 i τ2 χ B + O h7/2 , 





κ02 − 2κ1 i χ B + O h7/2 .

Eq. (19) for w 2,1 thus becomes C 2 (h) w 22,1 + C 1 (h) w 2,1 + C 0 (h) = 0, where

C 2 (h) =

+,+ +,−

w 2,1 =

C 1 (h) = O ( h ),

C 0 (h) = O (h).

χ B = 1 are

The solutions for

w 2,1 =

√

τ12 + τ22 + (τ1 − 3)τ2 , 15(τ1 − 1)τ2 √

h+

1 4

κ0 (τ1 + τ2 )h3/2 i −

1 8



(τ12 − 9τ2 τ1 + τ2 (τ2 + 7)) √

τ + τ2 τ1 + (τ2 − 3)τ2 2 1





κ02 − 2 κ1 i τ1 τ2 h5/2 + O h7/2 ,

h + κ0

τ13 − 3τ2 τ12 − 3τ2 (τ2 + 1)τ1 + τ2 (τ22 + 2τ2 + 5) 3/2 h i 4(τ12 + τ2 τ1 + (τ2 − 3)τ2 )

 τ2 (τ22 + (1 − 3τ1 )τ2 + 1) 5/2    − κ02 − 2 κ1 i h + O h7/2 , 2 8(τ1 + τ2 τ1 + (τ2 − 3)τ2 ) and for

χ B = −1 are −,+

√



−,−

√



w 2,1 = c 1+ (τ1 , τ2 ) h + κ0 c 2+ (τ1 , τ2 )h3/2 i + w 2,1 = c 1− (τ1 , τ2 ) h + κ0 c 2− (τ1 , τ2 )h3/2 i +













κ02 c3+ (τ1 , τ2 ) + κ1 c4+ (τ1 , τ2 ) i h5/2 + O h7/2 , κ02 c3− (τ1 , τ2 ) + κ1 c4− (τ1 , τ2 ) i h5/2 + O h7/2 ,

where c i± (τ1 , τ2 ), i = 1, 2, 3, 4, denote some rather complicated expressions involving only +,+

τ1 and τ2 . Suppose first that

w 2,1 = w 2,1 . Then the triarc (4) expands as +,+

pi

1

1

2

6

(t ) = ht + κ0 t 2 h2 i −



 

κ02 − κ1 t 3 h3 + π4,i (t ; τ1 , τ2 , κ0 , κ1 , κ2 )h4 + O h5 , t ∈ [τi−1 , τi ],

for i = 1, 2, 3, where π4,i is a complex polynomial in t of degree 4 with coefficients involving and (20) then follows that

(21)

τ1 , τ2 , κ0 , κ1 , κ2 . From (21)

 +,+    p − f ◦ ϕ  = O h4 , i.e., the approximation order is equal to four. For the other three solutions the triarc (4) expands as

 

+,−

p +,− (t ) = π5,i (t ; τ1 , τ2 )h + O h2 ,

  −,± p −,± (t ) = π5,i (t ; τ1 , τ2 )h + O h2 , +,−

−,+

t ∈ [τi −1 , τi ], i = 1, 2, 3,

(22)

t ∈ [τi −1 , τi ], i = 1, 2, 3,

(23)

−,−

where π5,i , π5,i and π5,i are complex polynomials in t of degree 5 with coefficients involving only τ1 and checked that none of the polynomials can be equal to t h for any 0 < τ1 < τ2 < 1. Therefrom it follows that

 +,−  p − f ◦ ϕ  = O(h),

τ2 . It can be

  −,± p − f ◦ ϕ  = O(h),

which implies the approximation order one. Let us summarize the results in the following theorem. Theorem 1. Let f : [0, h] → C, s → f (s) be an analytic complex curve and ϕ : [0, 1] → [0, h], t → th, a linear reparameterization. Moreover, by p ±,± : [0, 1] → C let us denote four interpolating triarcs defined by (4), that satisfy

p ±,± (i ) = ( f ◦ ϕ )(i ),





p ±,± (i ) = ( f ◦ ϕ ) (i ),



p ±,±



(i ) = ( f ◦ ϕ ) (i ),

Then for any τ1 , τ2 , 0 < τ1 < τ2 < 1, the triarc p +,+ has the approximation order four, i.e.,

 +,+    p − f ◦ ϕ  = O h4 . For the other three triarcs, the approximation order is equal to one, i.e.,

 +,−  p − f ◦ ϕ  = O(h),

  −,± p − f ◦ ϕ  = O(h).

i = 0, 1 .

B. Bastl et al. / Computer Aided Geometric Design 31 (2014) 412–426

419

Fig. 2. The curve defined in (24) and its approximation by C 2 PH quintic triarc (based on p +,+ solution) for n segments: n = 1 (top left), n = 2 (top right), n = 4 (bottom left) and n = 8 (bottom right).

To confirm the theoretical asymptotic results by a numerical example, let us choose a smooth curve f and a reparameterization ϕ as (see Fig. 2)

f : [0, h] → C,



s → log(s + 1) cos (2s),



s2 + 1 log (s + 1) sin (2s)



,

ϕ : [0, 1] → [0, h], t → ht .

(24)

The data becomes

A = ( f ◦ ϕ )(0),

t A = ( f ◦ ϕ ) (0),

c A = ( f ◦ ϕ ) (0),

B = ( f ◦ ϕ )(1),

t B = ( f ◦ ϕ ) (1),

c B = ( f ◦ ϕ ) (1),

(25)

for some h small enough. Let e ±,± (h) :=  p ±,± − f ◦ ϕ ∞,[0,1] denote the distance between the triarc p ±,± and the given ±,±

curve. To determine the approximation order numerically, we assume that e ±,± (h) ∼ const · hγ and we estimate the decay exponent γ ±,± by comparing the errors of interpolating triarcs for different values of h. In Table 1 the errors e ±,± (h) are shown for all four triarcs and h = 2−k , k = 0, 1, . . . , 7. Moreover, the decay exponents computed from two consecutive measurements are shown too. Since γ ±,± tends to the order of approximation as h approaches zero, the results confirm that the approximation order for p +,+ is four and for the other three solutions it is just one. 3.3. Numerical examination of free parameters τ1 and τ2 In this section the numerical investigation of the influence of the break points τ1 and τ2 on the shape of the interpolant is considered. Without loosing generality, we can assume that τ0 = 0 and τ3 = 1. In Farouki (2008), the use of the elastic bending energy functional

1 U= 0





κ 2 (t ) p  (t )dt

(26)

420

B. Bastl et al. / Computer Aided Geometric Design 31 (2014) 412–426

Table 1 Distances e ±,± (h) of the curve (24) from its triarc interpolants p ±,± and corresponding parameters As it is clearly seen from the table, only the solution p+,+ gives fourth order approximation.

γ ±,± showing numerical values of the decay exponent.

h

e +,+ (h)

γ +,+

e +,− (h)

γ +,−

e −,+ (h)

γ −,+

e −,− (h)

1

0.005653

/

0.3403

/

0.3506

/

0.2399

/

1 2 1 4 1 8 1 16 1 32 1 64 1 128

0.0002240

4.657

0.1140

1.578

0.1242

1.497

0.1075

1.158

0.00001375

4.026

0.04577

1.317

0.05314

1.225

0.05426

0.986

1.174 · 10−6

3.550

0.02292

0.998

0.02684

0.986

0.02775

0.967

8.334 · 10−8

3.816

0.01154

0.990

0.01379

0.961

0.01409

0.978

5.422 · 10−9

3.942

0.005798

0.993

0.007019

0.974

0.007103

0.988

3.436 · 10−10

3.980

0.002907

0.996

0.003545

0.985

0.003567

0.994

2.160 · 10−11

3.992

0.001456

0.998

0.001782

0.992

0.001788

0.997

Fig. 3. Values of a functional U with respect to

γ −,−

τ1 and τ2 for triarcs p +,+ (top left), p +,− (top right), p −,+ (bottom left) and p −,− (bottom right).

as the quantity whose minimum identifies the “good” Hermite interpolant is proposed and explicit formulae are derived for the case of PH quintics. In the following example the behaviour of the functional U with respect to τ1 and τ2 is examinated in order to construct triarcs with minimal elastic bending energy. Example 2. Suppose that C 2 Hermite data from Example 1 are given. The values of the functional U in the case 2 3

for the interpolants p +,+ , p +,− , p −,+ and p −,− are equal to U

τ1 = 13 , τ2 =

≈ 11.33, U ≈ 2261.58, U ≈ 162.45 and U ≈ 104.05, respectively. In Fig. 3 graphs of U for 0 = τ0 < τ1 < τ2 < 1 = τ3 for all four triarcs are shown. They indicate that the minimum of U is reached as the two parameters τi tend to the same value. In particular, for the triarc p +,+ the minimum value of U tends to 2.953 as τ1 → 0, τ2 → 1. Choosing τ1 = 0.001 and τ2 = 0.999 we compute that U ≈ 2.96555 and the corresponding triarc is shown in Fig. 4 (top left). For the other three solutions the minimum is reached as τ1 → τ2 , see Fig. 3. We choose the values τ1 , τ2 that approximate the minimum of U as p +,− :

τ1 = 0.379,

τ2 = 0.38,

U ≈ 1164.91,

p

−,+

: τ1 = 0.605,

τ2 = 0.606,

U ≈ 134.898,

p

−,−

: τ1 = 0.855,

τ2 = 0.856,

U ≈ 68.547,

and we present the triarcs in Fig. 4.

B. Bastl et al. / Computer Aided Geometric Design 31 (2014) 412–426

421

Fig. 4. Four distinct quintic PH triarcs obtained by minimizing elastic bending energy: p +,+ (top left), p +,− (top right), p −,+ (bottom left) and p −,− (bottom right). For a comparison, PH triarcs with equidistant parameters are shown in light gray.

Elastic bending energy integral is not the only choice for the functional to be minimized. One can try to choose other functionals like the length, the curvature or the absolute rotation index of a triarc too, but numerical results lead to similar conclusions, i.e., the minimum is reached as the two parts of a triarc approach. Also, the value of the functional for the equidistant parameters does not differ much from the minimal value, which indicates that such a choice could be good for practical applications. 4. Spatial interpolation problem and algorithm In this section, we will show that it is straightforward to generalize the previous algorithm to the spatial case. The approach presented is similar to the generalization of planar C 1 Hermite interpolation by uniform and non-uniform TC-biarcs (see Bastl et al., 2013) to the spatial case (see Bastl et al., 2014). Recall from Section 2 that spatial PH curves can be generated with the help of quaternion polynomials. Therefore, points and vectors in R3 will be identified with pure quaternions and vice versa in the rest of this section. Let A, B, t A , t B , c A , c B ∈ R3 be pairs of given points, associated tangent vectors and second derivatives vectors, respectively. The task is to compute a C 2 continuous spatial quintic PH triarc

⎧ ⎨ p 1 (t ), t ∈ [τ0 , τ1 ], p : [τ0 , τ3 ] → R3 : p (t ) = p 2 (t ), t ∈ [τ1 , τ2 ], ⎩ p 3 (t ), t ∈ [τ2 , τ3 ],

τ0 < τ1 < τ2 < τ3 ,

interpolating given Hermite data. The modification of the planar to the spatial case is based on one principle: all complex numbers are replaced by quaternions and the standard multiplication of complex numbers is replaced by -operation defined on quaternions, see (1). The preimage Ai and the associated PH curves p i are represented in the Bernstein–Bézier form as

p i (t ) =

5 

P i , j B 5j

t − τ i −1

j =0

τ i −1

Ai (t ) =

,

2 

Ai , j B 2j

j =0

t − τ i −1

τ i −1

,

i = 1, 2, 3,

and it holds that

τ i − 1 A2i ,0 , P i ,2 = P i ,1 + Ai ,0  Ai ,1 , 5 5

τ i − 1 2 2  1 P i ,3 = P i ,2 + Ai ,1 + Ai ,0  Ai ,2 , P i ,1 = P i ,0 +

τ i − 1

5

P i ,4 = P i ,3 +

τ i − 1 5

3

3

Ai ,1  Ai ,2 ,

P i ,5 = P i ,4 +

τ i − 1 5

A2i ,2 .

The interpolating conditions that determine the spatial triarc are of the same form as the conditions (9)–(14) for the planar case, only analogously modified. Clearly, P 1,0 = A and P 3,5 = B. The first order continuity conditions in (9), (10) and (14) can be solved with the help of Lemma 2. Namely,

422

B. Bastl et al. / Computer Aided Geometric Design 31 (2014) 412–426

A1,0 = X (ϕ A ; t A ),

A3,2 = X (ϕ B ; t B ),

A1,2 = A2,0 Qϕ1 ,

A3,0 = A2,2 Qϕ2 ,

(27)

where ϕ A , ϕ B , ϕ1 and ϕ2 are free angular parameters. The equations for the interpolation of second derivatives at the boundary, analogous to (9) and (10), can be solved using Lemma 1. The unknowns

A1,1 = X α1 ; A1,0 , t A +

τ0 4



A3,1 = X α2 ; A3,2 , t B −

cA ,

are determined with two new free parameters presented in (14), simplify to

A2,0 

1

τ0

A2,2 

1

τ1

(A2,0 − A1,1 Q−ϕ1 ) = A2,0  (A2,2 − A2,1 ) = A2,2 

1

τ2

τ2 4

cB

(28)

α1 and α2 involved. The second order continuity conditions at τ1 and τ2 , 1

τ1

(A2,1 − A2,0 ),

(A3,1 Q−ϕ2 − A2,2 ),

which by following the planar case should determine the unknowns A2,0 and A2,2 . To obtain a regular curve we must assume A2,0 = 0 and A2,2 = 0. By using Lemma 2 they can be expressed with two additional free angle parameters. To simplify further computations, we can take a particular solution which is to equate the quaternions on the right side of the star operation instead. This gives

A2,0 =

1

τ0 + τ1

( τ0 A2,1 + τ1 A1,1 Q−ϕ1 ),

A2,2 =

1

τ1 + τ2

( τ1 A3,1 Q−ϕ2 + τ2 A2,1 ).

(29)

Finally, we end up with one quadratic -equation, analogous to (19), for the remaining control point A2,1 of the preimage: 3 



2 1 τi −1 A2i ,0 + Ai ,0  Ai ,1 + A2i ,1 + Ai ,0  Ai ,2 + Ai ,1  Ai ,2 + A2i ,2 − 5( B − A ) = 0. 3

i =1

3

(30)

Solving this equation with the help of Lemma 2 introduces another free angular parameter ψ . To sum up, the control points of the preimage quaternion curve are completely determined by (27)–(29) and by a solution of (30). They are expressed with seven free parameters ϕ A , ϕ B , ϕ1 , ϕ2 , α1 , α2 , ψ . However, analogously to Bastl et al. (2014), it is possible to prove with the help of (3) that one angular parameter is superfluous. We can choose e.g., ψ = 0. Note that the number of free parameters coincides with the difference between the number of equations and unknowns that are involved in C 2 quintic triarc construction. Namely, each quintic PH curve is determined by 14 parameters, the quintic PH triarc thus has p = 42 parameters of freedom. Interpolation conditions give e = 4 · 9 = 36 equations, and the difference is p − e = 6. To sum up, for any non-degenerate spatial C 2 Hermite data there exists a six-parametric family of interpolating quintic PH triarcs p τ (t ; ϕ A , ϕ B , ϕ1 , ϕ2 , α1 , α2 ). Example 3. Let us consider C 2 Hermite data

P 0 = (2, 1, 1) ,



c0 =

1 1

, ,0

2 2

P 2 = (1, 2, 1) ,





,

c2 =

1 2

1

,− ,0 2

t 0 = (1, 1, 1) ,



t 2 = (1, 1, −2)

.

Fig. 5 (left) shows quintic PH triarcs obtained for pairs ϕ A , ϕ B ∈ {−π , − π2 , 0, π2 } and ϕ1 , ϕ2 , α1 , α2 equal to zero. Fig. 5 (right) shows quintic PH triarcs obtained for pairs α1 , α2 ∈ {0, 2, 4, 6}, and parameters ϕ A , ϕ B , ϕ1 , ϕ2 set to zero. In both cases, uniform quintic PH triarcs (i.e., τ1 = 13 , τ2 = 23 ) are shown. Analysis of this spatial algorithm is very similar to C 1 Hermite interpolation by spatial PH cubic biarcs (see Bastl et al., 2014) and C 2 Hermite interpolation by PH curves of degree 9 (see Šír and Jüttler, 2007). With exactly the same reasoning as in Šír and Jüttler (2007), Bastl et al. (2014) it is possible to prove the following Proposition 1. For any planar C 2 Hermite data P i , t i , c i , i = 0, 2, the four quintic PH triarcs p τ (t ; 0, 0, 0, 0, 0, 0), p τ (t ; 0, −π , 0, 0, 0, 0), p τ (t ; −π , 0, 0, 0, 0, 0), p τ (t ; −π , −π , 0, 0, 0, 0) are planar and correspond to p ±,± solutions of the planar interpolation problem with quintic PH triarcs, presented in Section 3, respectively. Proof. If we set α1 = 0, α2 = 0, ϕ1 = 0 and ϕ2 = 0, then the proof is completely analogous to the proofs of Proposition 1 in Bastl et al. (2014) and Theorem 3.11 in Šír and Jüttler (2007). 2

B. Bastl et al. / Computer Aided Geometric Design 31 (2014) 412–426

423

Fig. 5. PH triarcs obtained for spatial input Hermite data in Example 3 (tangent vectors are scaled down by factor 3).

Fig. 6. Uniform quintic PH triarcs obtained for planar input Hermite data in Example 4 (tangent vectors are scaled down by factor 3).

Example 4. Let us consider C 2 Hermite data 

P 0 = (0, 0, 0) ,

t 2 = 2,

2 3

 ,0

P2 =

3 2

 , 0, 0

c 0 = (−1, 1, 0) ,

t 0 = (2, 2, 0) ,

,

c 2 = (1, −1, 0) .

Fig. 6 shows four quintic PH triarcs obtained for

τ2 =

2 . 3

ϕ A , ϕ B ∈ {−π , 0}, parameters ϕ1 , ϕ2 , α1 , α2 equal to zero and τ1 = 13 ,

These solutions exactly correspond to the solutions of the planar problem.

Similarly as in the planar case one can analyze the asymptotic behaviour of the interpolating triarc in order to determine the set of parameters providing the best approximation order. By direct computation, analogously to the analysis of an approximation order in Bastl et al. (2014), the following theorem can be proved. Nevertheless, we will skip the details of the proof here, because expressions become much more complicated. Theorem 2. Let f : [0, h] → H, s → f (s) be an analytic quaternion curve and ϕ : [0, 1] → [0, h], t → th, a linear reparameterization. Moreover, by p τ (t ; ϕ A , ϕ B , ϕ1 , ϕ2 , α1 , α2 ) : [0, 1] → H let us denote interpolating quintic PH triarcs that satisfy

p τ (t ; ϕ A , ϕ B , ϕ1 , ϕ2 , α1 , α2 )|t =i = ( f ◦ ϕ )(i ),

 



 p τ (t ; ϕ A , ϕ B , ϕ1 , ϕ2 , α1 , α2 ) t =i = ( f ◦ ϕ ) (i ),

 

 p τ (t ; ϕ A , ϕ B , ϕ1 , ϕ2 , α1 , α2 ) t =i = ( f ◦ ϕ ) (i ),

i = 0, 1 .

424

B. Bastl et al. / Computer Aided Geometric Design 31 (2014) 412–426

Table 2 The approximation errors e and the corresponding decay exponents order approximation for p τ (t ; π2 , − π2 , 0, 0, 0, 0).

γ indicating the fourth order approximation for pτ (t ; 0, 0, 0, 0, 0, 0) and just the first p τ (t ; π2 , − π2 , 0, 0, 0, 0)

p τ (t ; 0, 0, 0, 0, 0, 0)

h

e (h)

γ

e (h)

γ

1

0.002474

/

0.1519

/

1 2 1 4 1 8 1 16 1 32 1 64 1 128

0.0002292

3.432

0.03733

2.024

0.00001522

3.913

0.01301

1.521

8.388 · 10−7

4.181

0.006332

1.038

4.724 · 10−8

4.150

0.003142

1.011

2.786 · 10−9

4.084

0.001567

1.004

1.691 · 10−10

4.042

0.0007826

1.001

1.042 · 10−11

4.021

0.0003912

1.001

Then for any τ1 , τ2 , 0 < τ1 < τ2 < 1, the triarc p τ (t ; 0, 0, 0, 0, 0, 0) has the approximation order four, i.e.,

     p τ (t ; 0, 0, 0, 0, 0, 0) − f ◦ ϕ  = O h4 . In the following example, we also give a numerical evidence of the above mentioned fact and we show that for other, non-zero, choices of the parameters ϕ A , ϕ B , the approximation order is generally only one. Example 5. We are given an analytic curve f and a linear reparameterization

f : [0, h] → R3 ,



s → log(s + 1) cos(s), log(s + 1) sin(s), 1 + s2



ϕ, ,

ϕ : [0, 1] → [0, h], t → th. The data are sampled as (25) for some h small enough. Let τ1 = 13 , τ2 = 23 and let p τ (t ; 0, 0, 0, 0, 0, 0) denote the interpolating quintic PH triarc with all six parameters, ϕ A , ϕ B , ϕ1 , ϕ2 , α1 , α2 set to zero. Furthermore, let e (h) :=  p τ − f ◦ ϕ ∞,[0,1] denote the distance between the triarc p τ and the given curve. In Table 2 the errors e and the corresponding decay exponents γ are shown. The third column numerically indicates that for p τ (t ; 0, 0, 0, 0, 0, 0) the approximation order of the curve approximation is four. In case when all free parameters are not set to zero, the approximation order falls. For example, let ϕ A = π2 , ϕ B = − π2 , ϕ1 = ϕ2 = α1 = α2 = 0. The associated errors and the decay exponents are presented in Table 2 which confirms that the approximation order for p τ (t ; π2 , − π2 , 0, 0, 0, 0) is just one. As in the planar case, the free parameters τ1 and τ2 can be used to optimize the shape of the triarc. Minimization of the elastic bending energy integral is considered in the next example. Example 6. Suppose that the C 2 Hermite data from Example 3 are given. Let us consider the PH triarcs p τ (t ; 0; 0; 0; 0; 0; 0) that ensure the best approximation order only. The value of the functional U , defined in (26), is for the choice τ0 = 0, τ1 = 13 , τ2 = 23 , τ3 = 1 equal to U ≈ 16.505. Similarly as in the planar case, the minimum of U is reached as τ1 → 0 and τ2 → 1. Choosing τ1 = 0.001 and τ2 = 0.999 we compute that U ≈ 13.6895. The corresponding triarc p τ (t ; 0; 0; 0; 0; 0; 0) is shown in Fig. 7 (right) together with the triarc obtained by the equidistant choice of parameters (left). 5. Conclusion In this paper, the problem of C 2 Hermite interpolation of planar and spatial data (two points with associated first and second derivatives) by PH quintic triarcs was considered. More precisely, a PH quintic triarc curve was constructed by joining three PH quintics such that the first curve interpolates given C 2 Hermite data at one side, the third one interpolates the same type of given data at the other side and the middle arc is joined together with C 2 continuity to the first and the third arc. The method was formulated firstly for planar PH quintic triarcs using complex representation and consequently using the quaternion representation also for spatial PH quintic triarcs. For any set of C 2 planar boundary data we constructed four possible interpolants and, similarly, for any set of C 2 spatial boundary data we found a six-dimensional family of interpolating quintic PH triarcs. It was shown that the best possible approximation order is 4. Several numerical examples were presented which confirmed that the constructed triarcs are of good shape and can be used for the approximation of planar/spatial parametric curves by PH quintic spline curves. One of the advantages of the proposed interpolation scheme is that it enables to drop the degree of the used PH curve from 9 to 5.

B. Bastl et al. / Computer Aided Geometric Design 31 (2014) 412–426

Fig. 7. Two interpolating triarcs p τ (t ; 0; 0; 0; 0; 0; 0). The left one corresponds to τ1 = 13 , which ensure the minimal bending energy (tangent vectors are scaled down by factor 3).

τ2 =

2 3

and the right one corresponds to

425

τ1 = 0.001, τ2 = 0.999,

Acknowledgements All authors were supported by the Czech/Slovenian project KONTAKT MEB 091105. B. Bastl and K. Michálková were supported by Technology agency of the Czech Republic through the project TA03011157. We thank the referees for their valuable comments which have helped us to improve the manuscript. References Albrecht, G., Farouki, R.T., 1996. Construction of C 2 Pythagorean-hodograph interpolating splines by the homotopy method. Adv. Comput. Math. 5 (4), 417–442. Bastl, B., Slabá, K., Byrtus, M., 2013. Planar C 1 Hermite interpolation with uniform and non-uniform TC-biarcs. Comput. Aided Geom. Des. 30 (1), 58–77. Bastl, B., Bizzarri, M., Krajnc, M., Láviˇcka, M., Slabá, K., Šír, Z., Vitrih, V., Žagar, E., 2014. C 1 Hermite interpolation with spatial Pythagorean-hodograph cubic biarcs. J. Comput. Appl. Math. 257, 65–78. Byrtus, M., Bastl, B., 2010. G 1 Hermite interpolation by PH cubics revisited. Comput. Aided Geom. Des. 27 (8), 622–630. ˇ Cernohorská, E., Šír, Z., 2010. Support function of Pythagorean hodograph cubics and G 1 Hermite interpolation. In: Mourrain, B., Schaefer, S., Xu, G. (Eds.), Advances in Geometric Modeling and Processing. In: Lect. Notes Comput. Sci., vol. 6130. Springer, Berlin, Heidelberg, pp. 29–42. Choi, H.I., Farouki, R.T., Kwon, S.-H., Moon, H.P., 2008. Topological criterion for selection of quintic Pythagorean-hodograph Hermite interpolants. Comput. Aided Geom. Des. 25 (6), 411–433. Farouki, R.T., 1994. The conformal map z → z2 of the hodograph plane. Comput. Aided Geom. Des. 11 (4), 363–390. Farouki, R.T., 2008. Pythagorean-Hodograph Curves: Algebra and Geometry Inseparable. Geom. Comput., vol. 1. Springer. Farouki, R., Neff, C., 1995. Hermite interpolation by Pythagorean-hodograph quintics. Math. Comput. 64 (212), 1589–1609. Farouki, R.T., Peters, J., 1996. Smooth curve design with double-Tschirnhausen cubics. Ann. Numer. Math. 3 (1–4), 63–82. Farouki, R., Sakkalis, T., 1990. Pythagorean hodographs. IBM J. Res. Dev. 34 (5), 736–752. Farouki, R., Sakkalis, T., 1994. Pythagorean-hodograph space curves. Adv. Comput. Math. 2 (1), 41–66. Farouki, R.T., al Kandari, M., Sakkalis, T., 2002. Hermite interpolation by rotation-invariant spatial Pythagorean-hodograph curves. Adv. Comput. Math. 17 (4), 369–383. Farouki, R.T., Manni, C., Sestini, A., 2003. Spatial C 2 PH quintic splines. In: Lyche, T., Mazure, M., Schumaker, L. (Eds.), Curve and Surface Design. St. Malo 2002. Nashboro Press, Brentwood, pp. 147–156. Farouki, R.T., Giannelli, C., Manni, C., Sestini, A., 2008. Identification of spatial PH quintic Hermite interpolants with near-optimal shape measures. Comput. Aided Geom. Des. 25 (4–5), 274–297. Jakliˇc, G., Kozak, J., Krajnc, M., Vitrih, V., Žagar, E., 2010. On interpolation by planar cubic G 2 Pythagorean-hodograph spline curves. Math. Comput. 79 (269), 305–326. Jakliˇc, G., Kozak, J., Krajnc, M., Vitrih, V., Žagar, E., 2012. An approach to geometric interpolation by Pythagorean-hodograph curves. Adv. Comput. Math. 37 (1), 123–150. Jüttler, B., 2001. Hermite interpolation by Pythagorean hodograph curves of degree seven. Math. Comput. 70 (235), 1089–1111. Jüttler, B., Mäurer, C., 1999. Cubic Pythagorean-hodograph spline curves and applications to sweep surface modeling. Comput. Aided Des. 31 (1), 73–83. Kubota, K., 1972. Pythagorean triples in unique factorization domains. Am. Math. Mon. 79 (5), 503–505. Kwon, S.-H., 2010. Solvability of G 1 Hermite interpolation by spatial Pythagorean-hodograph cubics and its selection scheme. Comput. Aided Geom. Des. 27 (2), 138–149. Meek, D.S., Walton, D.J., 1997. Geometric Hermite interpolation with Tschirnhausen cubics. J. Comput. Appl. Math. 81 (2), 299–309. 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. Pelosi, F., Farouki, R.T., Manni, C., Sestini, A., 2005. Geometric Hermite interpolation by spatial Pythagorean-hodograph cubics. Adv. Comput. Math. 22 (4), 325–352. Pelosi, F., Sampoli, M., Farouki, R., Manni, C., 2007. A control polygon scheme for design of planar C 2 PH quintic spline curves. Comput. Aided Geom. Des. 24 (1), 28–52. Sakkalis, T., Farouki, R.T., 2012. Pythagorean-hodograph curves in Euclidean spaces of dimension greater than 3. J. Comput. Appl. Math. 236 (17), 4375–4382. Sestini, A., Landolfi, L., Manni, C., 2013. On the approximation order of a space data-dependent PH quintic Hermite interpolation scheme. Comput. Aided Geom. Des. 30 (1), 148–158. Shafarevich, I.R., 1974. Basic Algebraic Geometry. Springer-Verlag.

426

B. Bastl et al. / Computer Aided Geometric Design 31 (2014) 412–426

Šír, Z., Jüttler, B., 2005. Spatial Pythagorean hodograph quintics and the approximation of pipe surfaces. In: Martin, R., Bez, H., Sabin, M. (Eds.), Mathematics of Surfaces XI. In: Lect. Notes Comput. Sci., vol. 3604. Springer, Berlin, Heidelberg, pp. 364–380. Šír, Z., Jüttler, B., 2007. C 2 Hermite interpolation by Pythagorean hodograph space curves. Math. Comput. 76 (259), 1373–1391. Šír, Z., Feichtinger, R., Jüttler, B., 2006. Approximating curves and their offsets using biarcs and Pythagorean hodograph quintics. Comput. Aided Des. 38 (6), 608–618.