Motion design with Euler–Rodrigues frames of quintic Pythagorean-hodograph curves

Available online at www.sciencedirect.com

Mathematics and Computers in Simulation 82 (2012) 1696–1711

Original Articles

Motion design with Euler–Rodrigues frames of quintic Pythagorean-hodograph curves Marjeta Krajnc a,b,∗ , Vito Vitrih c,d a c

FMF, University of Ljubljana, Jadranska 19, Ljubljana, Slovenia b IMFM, Jadranska 19, Ljubljana, Slovenia FAMNIT, University of Primorska, Glagoljaˇska 8, Koper, Slovenia d IAM, University of Primorska, Muzejski trg 2, Koper, Slovenia

Received 2 February 2011; received in revised form 29 March 2012; accepted 5 April 2012 Available online 27 April 2012

Abstract The paper presents an interpolation scheme for G1 Hermite motion data, i.e., interpolation of data points and rotations at the points, with spatial quintic Pythagorean-hodograph curves so that the Euler–Rodrigues frame of the curve coincides with the rotations at the points. The interpolant is expressed in a closed form with three free parameters, which are computed based on minimizing the rotations of the normal plane vectors around the tangent and on controlling the length of the curve. The proposed choice of parameters is supported with the asymptotic analysis. The approximation error is of order four and the Euler–Rodrigues frame differs from the ideal rotation minimizing frame with the order three. The scheme is used for rigid body motions and swept surface construction. © 2012 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Pythagorean-hodograph; Euler–Rodrigues frame; Rotation minimizing frame; Motion design; Quaternion; Hermite interpolation

1. Introduction To compute an orthonormal frame of a spatial curve r is an important task in computer animation, motion planning, swept surface construction, etc. Frames determine an orientation of a rigid body as it traverses the curve. Typically an adapted frame (f 1 , f 2 , f 3 ) is searched for, which has the property that f 1 = r˙ /˙r = t is a tangent vector, and the remaining two vectors span the normal plane. A well known adapted frame is the Frenet frame (t, n, b) (see [12]), but it is often unsuitable for practical applications since it is not defined at inflection points and it incurs an unnecessary rotation of the normal plane vectors n and b around t. The most attractive frame in motion design applications and swept surface construction is a rotation minimizing frame (RMF frame), which is characterized through a solution of first-order differential equations (see [11], e.g.). More precisely, there should be no instantaneous rotation of f 2 and f 3 around f 1 = t. The variation of any adapted frame (f 1 , f 2 , f 3 ) along the curve r is determined by the angular velocity vector ω as f˙1 = ω × f 1 , ∗

f˙2 = ω × f 2 ,

f˙3 = ω × f 3 .

Corresponding author at: FMF, University of Ljubljana, Jadranska 19, 1000 Ljubljana, Slovenia. Tel.: +386 1 4766 667; fax: +386 1 2517 281. E-mail address: [email protected] (M. Krajnc).

0378-4754/$36.00 © 2012 IMACS. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.matcom.2012.04.003

M. Krajnc, V. Vitrih / Mathematics and Computers in Simulation 82 (2012) 1696–1711

1697

We can write ω = ω1 f 1 + ω2 f 2 + ω3 f 3 , where ω1 = f˙2 · f 3 = −f 2 · f˙3 , ω2 = f˙3 · f 1 = −f 3 · f˙1 , ω3 = f˙1 · f 2 = −f 1 · f˙2 . The property of the RMF frame is that ω · f 1 = 0, or equivalently ω1 = 0. Operations × and · denote the standard cross and scalar products. An important property for computer aided design applications is that a frame is rationally dependent. Therefore it is necessary that a curve r is a Pythagorean-hodograph curve (PH curve). PH curves are characterized by the property that the Euclidean norm of their hodograph is a piecewise polynomial (see Ref. [3]). As a consequence, they have a rational unit tangent, rational offset, polynomial arc length, etc. But clearly, the PH property does not ensure a rational RMF frame (RRMF frame) of a curve, and a construction of such curves is a difficult task since the nonlinear constraints are involved. For more results on a construction, applications and rational approximation of RMF curves see [4,6,7,9,13], and the citations quoted therein. There exists a special adapted frame defined particularly on spatial PH curves, so called Euler–Rodrigues frame (ERF frame) that has been introduced in Ref. [1]. This frame has the property of being rational by construction and it is always nonsingular at inflection points. Unfortunately, the ERF frame is not in general the RMF frame. In Ref. [1], it is shown that for cubic PH curves ERF and Frenet frames are the same. It is also shown, that a true spatial PH curve with the ERF frame being an RMF frame has to be of degree at least 7. In this paper we present a scheme for interpolation of G1 Hermite motion data, i.e., interpolation of data points and rotations at the points, with quintic PH curves so that the ERF frame of the curve coincides with the rotations at the points. The scheme can be applied for a construction of a rigid body motion with some given positions interpolated, where each position is represented by a position of a center of a rigid body (point in R3 ) and the orientation with respect to some fixed coordinate system (rotation in R3 ). The obtained quintic PH curve is expressed with three free parameters. Two of them are used to control the length of the curve or to provide a C1 interpolant, while the remaining one is computed in such a way that the ERF frame is as close as possible to the corresponding RMF frame. Note that there can exist a quintic PH curve interpolating the G1 Hermite motion data with the exact RRMF frame, but with no degrees of freedom left. This problem is studied in Ref. [6], but it is of a highly nonlinear nature and no a priori guarantee for the existence of the interpolant is provided, except for the asymptotic data. Furthermore, it is easy to construct examples where no interpolant exists or where the shape of the curve is not pleasant to the eye. The advantages of the quintic PH interpolant derived in this paper are that the nonlinear part of the problem has a closed form solution, the frame is rational by construction and has a minimal rotation of normal plane vectors around the tangent among all the possible PH quintics with the same end tangents and with an associated ERF interpolating the assigned end rotations. Furthermore, there are no restrictions on the data for the solution to exist and the main numerical work required is to solve a minimization problem for one variable with a very good starting value provided. The results are supported by the asymptotic analysis. It is shown that the approximation error is of order four and that the ERF frame differs from the RMF frame with the order three. Numerical examples confirm that the interpolant is competitive with interpolants obtained in Ref. [6]. The paper is organized as follows. In Section 2, a quaternion representation of spatial PH curves and a definition of the Euler–Rodrigues frame are briefly recalled. In the next section, the interpolation problem with a spatial PH quintic, which interpolates two end points and whose ERF frame coincides with given rotations at the points, is presented, and a three-parameter family of solutions is derived. Section 4 deals with a C1 interpolation with prescribed lengths of tangents at the interpolation points. An asymptotic analysis is provided, which offers the results also for the nonasymptotic applications. Further, a more general G1 interpolation is considered and several criteria how to choose free parameters are presented. Numerical examples from Section 5 illustrate the theoretical results. The paper is concluded with Section 6 that summarizes the main results of the paper and identifies possible future investigations. 2. Quaternion representation of PH curves Many representations of PH curves were proposed in the last few years. The planar PH curves can be defined through some relations between complex numbers, while the spatial PH curves by relations between quaternions or through a Hopf map representation (see Ref. [3]). Recently, particular algebraic equations, that identify a PH curve

1698

M. Krajnc, V. Vitrih / Mathematics and Computers in Simulation 82 (2012) 1696–1711

independently of the dimension of a space, have been derived in Ref. [10]. In this paper a quaternion representation is used, thus let us briefly recall some basic properties and notation needed further on. Quaternions form a 4-dimensional vector space H with a standard basis {1, i, j, k}, 1 = (1, 0, 0, 0),

i = (0, 1, 0, 0),

j = (0, 0, 1, 0),

k = (0, 0, 0, 1).

The first component of a quaternion is called a scalar part, while the remaining three components form a vector part of a quaternion. A quaternion with a zero scalar part is called a pure quaternion. Vectors in R3 can be identified with pure quaternions and vice versa. If we write A = (a, a) and B = (b, b), a, b ∈ R, a, b ∈ R3 , then A + B = (a + b, a + b),

AB = (ab − a · b, ab + ba + a × b).

With this associative, but noncommutative multiplication the space of quaternions becomes an algebra. Moreover, A := (a, −a) denotes a conjugate of A and the norm of a quaternion A is defined as A2 = AA = AA = a2 + a2 , √ where a = a · a is the Euclidean norm of the vector a. Furthermore, let us define a commutative multiplication on the space of quaternions as 1 A  B := (A i B + B i A). 2 Note that A  B is a pure quaternion and will be identified with a vector in R3 several times later on. A shorter notation A2 := A  A will be used, and the solution of a ‘quadratic’ -equation which is given in the next lemma will be needed. The proof can be found in Ref. [8]. Lemma 1. Let A be a given pure quaternion which is not a negative multiple of i. Then all the solutions of a quadratic -equation X2 = A form a one-parametric family  (A/A) + i X = X(φ) = A Q(φ), Q(φ) := cos φ + i sin φ, φ ∈ [−π, π). (A/A) + i √ For A/A = −i, the solution is X(φ) = A k Q(φ). Spatial PH curves and quaternions are connected in the following way (see [3]). A spatial PH curve p can be generated from a quaternion polynomial A(t) = a(t) + i b(t) + j c(t) + k d(t), where a, b, c and d are relatively prime polynomials of degree ≤n, as ˙ := h(t) := A(t)2 . p(t)

(1)

Here p˙ denotes the derivative of p with respect to t. Integration of the hodograph h gives a polynomial PH curve p of degree ≤2n + 1. The quaternion A(t) is called the preimage. Note that quaternions A(t) and A(t)Q(φ) generate the same hodograph. A special adapted frame for PH curves called the Euler–Rodrigues frame is defined as follows. Definition 1. The Euler–Rodrigues frame FER = (f 1 , f 2 , f 3 ) of a spatial PH curve generated from a quaternion polynomial A(t) is defined as f 1 (t) =

A(t) i A(t) , A(t)2

f 2 (t) =

A(t) j A(t) , A(t)2

f 3 (t) =

A(t) k A(t) . A(t)2

The ERF frame has a property of being defined at every point of a PH curve. Moreover, if p is of degree 2n + 1, then the ERF frame is rational of degree 2n. Since the quaternion representation of a PH curve is not unique, the ERF frame is not uniquely defined either. However, each pair of two ERF frames has a constant angular difference along the curve.

M. Krajnc, V. Vitrih / Mathematics and Computers in Simulation 82 (2012) 1696–1711

1699

At every point, any adapted frame is defined by three orthonormal vectors which form a rotation in R3 . Note that each rotation can be written in a form ⎛ 2 ⎞ q0 + q12 − q22 − q32 2(q1 q2 − q0 q3 ) 2(q1 q3 + q0 q2 ) ⎜ ⎟

⎜ 2(q1 q2 + q0 q3 ) q02 − q12 + q22 − q32 2(q2 q3 − q0 q1 ) ⎟ (2) Q i Q, Q j Q, Q k Q , ⎜ ⎟= ⎝ ⎠ 2 2 2 2 2(q1 q3 − q0 q2 ) 2(q2 q3 + q0 q1 ) q 0 − q1 − q2 + q3 where Q = (q0 , q1 , q2 , q3 ) ∈ H, Q = 1, and q0 , q1 , q2 , q3 are called normalized Euler parameters (see [2]). Antipodal quaternions (i.e., quaternions that differ only by the sign of their components) represent the same rotation. Any adapted frame can thus be at each point represented by two antipodal quaternions. The bijective map between the space of rotations and the space of unit quaternions with identified antipodal points is called a kinematic map. 3. Interpolation problem In this section a cornerstone problem for a construction of a quintic PH spline that interpolates the positions of a rigid body is presented. Let P 0 , P 1 ∈ R3 be two given points and let F0 = (f 0,1 , f 0,2 , f 0,3 ) ∈ R3×3 ,

F1 = (f 1,1 , f 1,2 , f 1,3 ) ∈ R3×3

be two given frames, such that f ,i  = 1,

f ,i · f ,j = 0,

Furthermore, let 3 Q := q,i i=0 ∈ H,

 = 0, 1,

i, j ∈ {1, 2, 3}, i = / j.

 = 0, 1,

denote two unit quaternions that correspond to the frames F ,  = 0, 1, and satisfy Q0 · Q1 > 0. Here · denotes the standard scalar product on 4-dimensional vectors. The task is to find a quintic PH curve p : [0, 1] → R3 that interpolates the given points P 0 , P 1 , p(0) = P 0 ,

p(1) = P 1 ,

(3)

and has the ERF frames at t = 0 and t = 1 equal to F0 and F1 , respectively. The curve, its hodograph and the preimage can be written in a Bernstein form as p(t) :=

5 

pi Bi5 (t),

h(t) :=

i=0

4 

hi Bi4 (t),

i=0

A(t) :=

2 

Ai Bi2 (t),

(4)

i=0

where pi , hi (pure quaternions) and Ai (quaternions) are the control points, and n j n Bj (t) = t (1 − t)n−j , j = 0, 1, . . . , n, j are the Bernstein basis polynomials. From the interpolation of the frames and from interpolation conditions (3) we obtain p0 = P 0 ,

p5 = P 1 ,

A0 = λ0 Q0 ,

A2 = λ1 Q1 ,

(5)

where λ0 , λ1 ∈ R are the unknown parameters. The control points of the curve and its hodograph are connected as hi = 5 pi ,

i = 0, 1, . . . , 4,

where pi := pi+1 − pi is a forward difference. Conditions (3) imply the equation 4  i=0

hi = 5 P 0 .

(6)

1700

M. Krajnc, V. Vitrih / Mathematics and Computers in Simulation 82 (2012) 1696–1711

Furthermore, the relation between the hodograph and the preimage h(t) = A(t)2 can be expressed as 1 h3 = λ1 A1  Q1 , h4 = λ21 Q2 (2A2 1 + λ0 λ1 Q0  Q1 ), 1 . (7) 3 Let A1 := (a1 , b1 , c1 , d1 ). Substituting (7) into (6) gives the nonlinear system of three equations for six unknowns a1 , b1 , c1 , d1 , λ0 and λ1 . It can be written as h0 = λ20 Q2 0 ,

h1 = λ0 Q0  A1 ,

h2 =

A2 1 + A1  B + c = 0, where

3 3 1 2 2 2 2 λ0 Q0 + λ1 Q1 + λ0 λ1 Q0  Q1 − 5 P 0 . B := (λ0 Q0 + λ1 Q1 ), c := 2 2 3 By transforming Eq. (8) into

2 1 1 A1 + B = B2 − c, 4 2 and by using Lemma 1, we obtain a one-parametric family of solutions for the unknown quaternion A1 :    1 2  (1/4B2 − c)/(1/4B2 − c) + i 1  − c B A1 = − B + Y Q(φ), Y :=    (1/4B2 − c)/(1/4B2 − c) + i , φ ∈ [−π, π). 2 4

(8)

(9)

The coefficients a1 , b1 , c1 and d1 are thus expressed with three free (shape) parameters λ0 , λ1 and φ and their choice clearly has a great influence on the shape of the resulting curve and on the properties of the ERF frame. Since the ERF frame of a quintic PH curve can be an RMF frame only if the curve is planar, our first aim is that the resulting ERF frame would be as close as possible to the ideal RMF frame. Another criteria that can be used to choose the parameters is to control the length of a curve p. The following theorem reveals that the length of p is independent of the parameter φ. Theorem 1. The length of the interpolating PH curve p defined by (4), (5) and (9) depends only on the parameters λ0 and λ1 and not on the parameter φ. Proof. From a quaternion representation of a PH curve the norm of the hodograph can be expressed as 4  σi Bi4 (t), σ(t) := h(t) = A(t)A(t) =: i=0

where

1 A0 A1 + A1 A0 , 2

1 σ2 = A0 A2 + 4A1 A1 + A2 A0 , 6

1 σ3 = A 1 A 2 + A 2 A 1 , σ4 = A 2 A 2 , 2 and the length of p simplifies to  1 4 1 h(t) dt = σi . (p) := 5 0 i=0     Recall (9) and note that Y2 =  41 B2 − c. It is straightforward to check that (10) is equal to σ0 = A0 A0 ,

(p) =

σ1 =

 1  15λ20 Q0 2 − 10λ0 λ1 Q0 · Q1 + 15λ21 Q1 2 + 16Y2 120    1 2  1 2 2 2 2  = 15λ0 Q0  − 10λ0 λ1 Q0 · Q1 + 15λ1 Q1  + 16  B − c  . 120 4

(10)

(11)

M. Krajnc, V. Vitrih / Mathematics and Computers in Simulation 82 (2012) 1696–1711

1701

Expression (11) is independent of φ which completes the proof.  Remark 1. Note that Theorem 1 is closely related to results obtained in [5, Section 7]. The RMF frame has a property that it does not rotate about the tangent of a curve. More precisely, the first component ω1 = f˙ 2 (t) · f 3 (t) of the angular velocity vector ω is identically equal to zero. The component ω1 of the ERF frame FER = f 1 , f 2 , f 3 of a curve determined through the quaternion A can be written as ω1 (t) = −4 where

α(1 − t)3 + (α − γ)(1 − t)2 t + (β − γ)(1 − t)t 2 + β t 3 , A(t)2

⎛ α = λ0 ⎝det



q0,0

a1

q0,1

b1

⎛



q0,2

c1

(12)

⎞⎞

⎠⎠ , q d 0,3 1 ⎛ ⎛ ⎞ ⎞⎞ ⎛ a1 q1,0 c1 q1,2 ⎠ − det ⎝ ⎠⎠ , β = λ1 ⎝det ⎝ b q d q 1 1,1 1 1,3      q1,0 q0,0 q1,2 q0,2 γ = λ0 λ1 det − det . q1,1 q0,1 q1,3 q0,3 − det ⎝

In order to approach the ERF frame of a curve as close as possible to the RMF frame of the same curve, we will minimize the integral  1 ω12 (t) dt. (13) min φ∈[−π,π) 0

Similar minimization has been done in [5,8], e.g. For a numerical computation of the minimum (13), the integral must first be approximated using some quadrature rule (Simpson, Gauss, etc.). A minimum of the approximated integral can then be computed using a gradient method or some other Newton type method. Asymptotic analysis from Section 4 provides very good starting values for this nonlinear minimization problem. Eq. (9) shows that there exists a three-parameter family of solutions of the interpolation problem considered. The following lemma reveals a symmetry between choosing ±λi . Lemma 2. Let p (t ; λ0 , λ1 ) denote the PH quintic interpolant with the parameter φ = φ(λ0 , λ1 ) computed by minimizing the integral (13). Then p (t; −λ0 , −λ1 ) = p (t; λ0 , λ1 ) and

 φ(−λ0 , −λ1 ) =

φ(λ0 , λ1 ) + π, φ(λ0 , λ1 ) ∈ [−π, 0) φ(λ0 , λ1 ) − π, φ(λ0 , λ1 ) ∈ [0, π)

.

Proof. It follows straightforwardly from (1), (5), (9) and (12).  4. Computation of free parameters 4.1. C1 interpolation The first column of the ERF frame at each point is the normalized tangent vector of a curve p. This implies that the spline composed of PH curves defined in the previous section is globally G1 . In some applications it might be more

1702

M. Krajnc, V. Vitrih / Mathematics and Computers in Simulation 82 (2012) 1696–1711

appropriate to have a C1 smoothness. Therefore the interpolation problem must be extended in such a way that a curve p also satisfies ˙ p(0) = d0,

˙ p(1) = d1,

where d 0 and d 1 are prescribed derivatives at interpolation points. By (7), this determines the values of λ0 and λ1 as   d i  λi = = d i , i = 0, 1. (14) Q2∗ i  It remains only to compute the parameter φ by minimizing (13). Let us now analyze this C1 interpolation scheme for the data taken from a smooth parametric curve. Let f : [0, h] → 3 R , s → f (s), be a smooth parametric curve parameterized by the arc-length with nonvanishing curvature. Further, let ψ : [0, 1] → [0, h],

t → ht,

be the reparameterization, and R(s) := (R1 (s), R2 (s), R3 (s)) ∈ R3×3 the RMF frame of f at s. Our goal is to compute a quintic PH polynomial curve p that interpolates the points p(0) = f (ψ(0)) = f (0) =: P 0 ,

p(1) = f (ψ(1)) = f (h) =: P 1 ,

(15)

the frames F0 = R(ψ(0)),

F1 = R(ψ(1)),

(16)

and also the derivatives dp(0) d(f ◦ ψ)(0) df (0) = =h , dt dt ds

df (h) dp(1) d(f ◦ ψ)(1) = =h . dt dt ds

(17)

Interpolation of derivative directions is already included in (16), but conditions (17) additionally imply the values for λ0 and λ1 . Without loos of generality, we may assume ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ 0 1 0 

 ⎜ ⎟ ⎜ ⎟ ⎜ ⎟



  (18) f (0) = ⎝ 0 ⎠ , f (0) = ⎝ 0 ⎠ , f (0) = f (0) ⎝ 1 ⎠ . 0 0 0 Suppose that the curvature κ and the torsion τ of the curve f at s = 0 expand as κ(s) = κ0 + κ1 s +

κ2 2 κ3 3 s + s + O(s4 ), 2 3!

τ(s) = τ0 + τ1 s +

τ2 2 τ3 3 s + s + O(s4 ), 2 3!

where κ0 > 0. Then (see [12]) the Frenet–Serret formulae give an expansion of the curve, simplified by the assumptions (18) to ⎞ ⎛ 1 1 s − κ02 s3 − κ0 κ1 s4 ⎟ ⎜ 6 8 ⎜  ⎟ ⎜1 2 1 3 1  3 2 4⎟ 5 ⎟ f (s) = ⎜ (19) ⎜ 2 κ0 s + 6 κ1 s + 24 −κ0 − τ0 κ0 + κ2 s ⎟ + O(s ). ⎟ ⎜ ⎠ ⎝ 1 1 3 4 (2κ ) κ0 τ0 s + 1 τ0 + κ0 τ1 s 6 24 The first component ω1 of the angular velocity vector ω of the Frenet  s frame is in the case of the arc-length parameterization equal to the torsion. From the angular deviation θ(s) := 0 τ(t)dt we obtain the RMF frame R as R1 = t,

R2 = cos θ n − sin θ b,

R3 = sin θ n + cos θ b,

M. Krajnc, V. Vitrih / Mathematics and Computers in Simulation 82 (2012) 1696–1711

where (t, n, b) is the Frenet frame of f . From R(0) and R(h), ⎛ ⎞ 1 0 0 ⎜ ⎟ R(0) = ⎝ 0 1 0 ⎠ , 0 0 1 ⎛ ⎞ 1 1 1 − κ0 h − κ 1 h 2 − κ 0 τ0 h 2 1 − κ02 h2 ⎜ ⎟ 2 2 2 ⎜ ⎟   ⎜ ⎟ 1 1 2 2 2 ⎜ ⎟ + O h3 , R(h) = ⎜ κ0 h + κ1 h 1 − κ0 h 0 ⎟ 2 2 ⎜ ⎟ ⎝1 ⎠ 2 0 1 κ0 τ0 h 2

1703

(20)

(21)

we can by (2) compute the corresponding unit quaternions Q0 and Q1 which expand as

1 1 1 1 1 − κ02 h2 − κ0 κ1 h3 , κ02 τ0 h3 , − κ0 τ0 h2 8 8 24 4

   1 1 1 1 1  3 −κ0 − 4τ02 κ0 + 4κ2 h3 + O h4 . + − κ1 τ0 − κ 0 τ1 h 3 , κ 0 h + κ 1 h 2 + 6 12 2 4 48

Q0 = (1, 0, 0, 0) ,

Q1 =

From (14) and (17) we obtain the parameters λ0 and λ1 as √ √ λ0 = h, λ1 = h.

(22)

(23)

The unknown quaternion A1 is by (9) given as A1 = A˜ 1 + O(h7/2 ), where ⎛ 1 ⎞ √ 3 − (5 sin(φ) + 3) h + (sin(φ) + 1)κ02 h5/2 ⎜ 2 ⎟ 32 ⎜ ⎟ ⎜5 ⎟ √ 3 ⎜ ⎟ 2 5/2 cos(φ)κ0 h ⎜ cos(φ) h − ⎟ ⎜ ⎟ 2 32 A˜ 1 = ⎜ ⎟, ⎜5 ⎟ 3 5/2 ⎜ cos(φ)κ0 h3/2 + ⎟ (cos(φ)κ1 + (sin(φ) + 1)κ0 τ0 ) h ⎜8 ⎟ 16 ⎜ ⎟ ⎝ ⎠ 1 3 ((sin(φ) + 1)κ1 − cos(φ)κ0 τ0 ) h5/2 − (5 sin(φ) + 3)κ0 h3/2 − 8 16 and ω1 , given in (12), expands as ω1 (t) =

25(t

10(2t − 1) cos(φ) + O(h2 ). − 5t + 1 + 5(t − 1)t sin(φ))2

− 1)2 t 2 cos2 (φ) + (5t 2

In order to have ω1 as small as possible for every t ∈ [0, 1] and h small enough, we have to choose φ = ± π/2. For φ = − π/2, ω1 (t) = and

 0

1

1 τ0 κ2 (6t 2 − 10t + 3)h3 + O(h4 ) 24 0

ω12 (t) dt =

23 4 2 6 κ τ h + O(h7 ). 8640 0 0

For φ = π/2, ω1 (t) =

(14t 2 − 10t + 3) 3 1 h + O(h4 ), τ0 κ02 24 (10t 2 − 10t + 1)2

(24)

1704

M. Krajnc, V. Vitrih / Mathematics and Computers in Simulation 82 (2012) 1696–1711

1 but the integral 0 ω12 (t)dt diverges, therefore this choice is not appropriate. For √ √ π λ0 = h, λ1 = h, φ = − , 2 the interpolant p expands as ⎞ ⎛ 1 th − t 3 κ02 h3 ⎟ ⎜ 6 ⎟ ⎜ ⎟ ⎜1 1 4 2 2 3 3 ⎜ p(t) = ⎜ t κ0 h + t κ1 h ⎟ ⎟ + O(h ), 6 ⎟ ⎜2 ⎠ ⎝1 t 3 κ0 τ 0 h 3 6 and the parametric distance reads as  h4 p − f ◦ ψ∞ = 36κ02 κ12 + (κ03 + 4τ02 κ0 − 4κ2 )2 + 16(2κ1 τ0 + κ0 τ1 )2 + O(h5 ). 1536 The results are summarized in the next theorem. Theorem 2. Let f : [0, h] → R3 , s → f (s), be a smooth parametric curve parameterized by the arc-length with the expansion (19), and let ψ : [0, 1] → [0, h], t → ht. The PH curve p that satisfies (15)–(17) and is defined by (4), (5) and (9), with parameters chosen as √ √ π λ0 = h, λ1 = h, φ = − , 2 approximates f with the asymptotic approximation order h4 . Furthermore, the first component ω1 of the angular 1 velocity vector ω tends to zero with the order h3 , and 0 ω12 (t) dt = O(h6 ). The curve f defined by (19) is rotated in such a way that the RMF frame at s = 0 is the identity matrix. Suppose that the curve is rotated by some matrix U, f → Uf , and let QU = (qU,0 , qU,1 , qU,2 , qU,3 ) be the corresponding quaternion. The data then change to P 0 → U P 0 ,

Q0 → QU Q0 ,

Q1 → QU Q1 ,

and the unknown quaternion A1 must therefore be multiplied by QU too. Unfortunately, the solution (9) is not invariant with respect to a quaternion multiplication and as a consequence the optimal angle φ changes. The asymptotic expansions in this case imply

qU,0 φ = c0 + c1 h + c2 h2 + O(h3 ), c0 = − arctan , (25) qU,1 where c1 and c2 are some constants that depend only on qU,i , κ0 , κ1 and τ 0 . Similarly as before, it can √more involved √ be proven that for λ0 = h, λ1 = h and φ chosen by (25) the conclusions of Theorem 2 hold. Remark 2. Asymptotic analysis provides a very good starting value φ for the nonlinear minimization problem (13) for general data. Namely,

q0,0 φstart = − arctan . (26) q0,1 The asymptotic results will now be confirmed with a numerical example. Let the data be sampled from a spatial curve T   f (s) = log(1 + s) cos s, log(1 + s) sin s, 1 + s2 , s ∈ [0, h], (27)

M. Krajnc, V. Vitrih / Mathematics and Computers in Simulation 82 (2012) 1696–1711

1705

Table 1 Approximation error between the curve (27) and the C1 interpolant p together with values of free parameters and the value of the integral (13). h

λ0

λ1

φ

φ − φstart

1

1 0.7071 0.5 0.3536 0.25 0.1768

1.0532 0.6706 0.4651 0.3364 0.2431 0.1742

−1.3258 −1.3369 −1.3387 −1.3389 −1.3390 −1.3390

0.0131 2.0372 × 10−3 3.1926 × 10−4 4.4659 × 10−5 5.8587 × 10−6 7.4821 × 10−7

1 2 1 4 1 8 1 16 1 32

h 1 1 2 1 4 1 8 1 16 1 32

1 0

ω12 (t)dt

5.6055 × 10−4

2.0905 × 10−5 5.7286 × 10−7 1.0861 × 10−8 1.8012 × 10−10 2.8712 × 10−12

Decay exponent

Approx. error

Decay exponent

– 4.7449 5.1895 5.7210 5.9141 5.9712

5.4300 × 10−3

– 3.0712 3.8960 4.1507 4.1447 4.0909

6.4607 × 10−4 4.3397 × 10−5 2.4433 × 10−6 1.3813 × 10−7 8.1064 × 10−9

at s = 0 and s = h. Recall that the frames F0 and F1 are computed as RMF frames of f at s = 0 and s = h. The parameters λi for the interpolant p are determined by (14) and (17), and φ is obtained by minimizing (13) using (26). The first part of Table 1 shows the values of λi , φ and the distance between the optimal φ from the starting value (26) for different values of h. In the second part of the table the values of the integral (13) and the parametric distance between f and p together with the decay exponent are shown. Note that the decay exponent is computed as the binary logarithm of the quotient of two consecutive measurements. It tends to the order of approximation as h approaches zero. The numbers numerically confirm the results of Theorem 2. Fig. 1(left) shows the RMF frame of the curve f and the ERF frame of the interpolant p for h = 1. The curves f and p are almost indistinguishable and so are the frames. The graph of ω1 at the optimal φ is shown in Fig. 1 (right). 4.2. G1 interpolation Let us now return to the problem from Section 3. A PH curve p is determined by three shape parameters λ0 , λ1 and φ ∈ [− π, π), that will be chosen in two steps. First the parameters λ0 and λ1 will be computed and secondly the optimal φ will be determined. The spline composed of such PH curves is globally G1 smooth. As stated in Theorem 1 the parameters λ0 and λ1 influence the length of the resulting curve p. Following the asymptotic results from Section 4.1 we can select λ0 and λ1 accordingly to (23) as  λi = P 1 − P 0 , i = 0, 1. (28)

0.04 0.02

0.2

0.4

0.6

0.8

1.0

0.02

Fig. 1. In the left figure, the RMF frame of the curve (27) for h = 1 (gray) and the ERF frame of the C1 interpolant (black) are shown. The right figure shows the graph of ω1 at the optimal φ.

1706

M. Krajnc, V. Vitrih / Mathematics and Computers in Simulation 82 (2012) 1696–1711

Another way to select λ0 and λ1 is to choose in advance the length (p) of the interpolant p. By prescribing the length of p, Eq. (11) determines the relation between parameters λ0 and λ1 that must be satisfied. The next lemma gives a closed form solution of Eq. (11) under the assumption λ0 = λ1 . Lemma 3. Let us assume that λ0 = λ1 and let L ∈ R, L≥  P 0 . Then (p) = L iff 



η + 2 ξ L − (η + 2 ξ L)2 − 4 ξ 2 − ζ L2 − P 0 2 √ λ0 = λ1 = ± Λ, Λ := , 2(ξ 2 − ζ)

(29)

where 1 1 S0,3 , ζ := (ρ2 + ρ22 + (S0,1 − S2,3 )2 ), 24 576 1 1 η := (S2,3 − S0,1 , ρ1 , ρ2 ) · P 0 , 12 ξ :=

and Sj,k :=

k 

2 2 (2 q0,i + 2 q1,i + (q0,i − q1,i )2 ),

i=j

ρ(x, y, z, w) := det



x −z

y w



 − 3 det

0 ≤ j < k ≤ 3,

x −w

y z

 ,

ρ1 := 2(ρ(q0,0 , q1,0 , q0,3 , q1,3 ) + ρ(q0,1 , q1,1 , q0,2 , q1,2 )), ρ2 := 2(ρ(q0,1 , q1,1 , q0,3 , q1,3 ) − ρ(q0,0 , q1,0 , q0,2 , q1,2 )). Proof. It is easy to see that ξ ≥ 0, ζ ≥ 0, η2 − 4 ζ  P 0  2 ≤ 0 and ξ 2 − ζ ≥ 0. The relation (11) can be simplified to  (p) − ξλ = ζλ2 + η λ + P 0 2 , (30) where λ := λ20 . It is straightforward to see that inserting λ = Λ into (30) implies (p) = L. Let us now assume that (p) = L. By squaring both sides of (30) we obtain a quadratic equation in λ. Using the cylindrical decomposition, one can prove that only the solution λ = Λ satisfies (30). In order to have λ0 = λ1 well defined, it remains to prove that Λ ∈ R and Λ ≥ 0. The term under the square root in Λ, X := (η + 2 ξ L)2 − 4(ξ 2 − ζ)(L2 − P 0 2 ), can be considered as a quadratic function in L having both roots complex or equal, since (ξ 2 − ζ)(η2 − 4ζ  P 0  2 ) ≤ 0. The positive leading coefficient implies that X ≥ 0 and thus Λ ∈ R. Under the assumption L≥  P 0 , it follows X ≤ (η + 2 ξ L)2 , which implies Λ ≥ 0. The proof is completed.  √ Remark 3. Using the result of Lemma 2, one can always choose λ0 = λ1 = + Λ. The influence of the length L on the shape of the interpolant p is shown in Fig. 2. For more examples see Section 5. Once the parameters λ0 and λ1 are determined, it remains to compute φ so that the ERF frame of a curve p will be as close as possible to the RMF frame of p. As suggested in Section 3, we have to minimize the integral (13). This minimization problem is clearly nonlinear. To find the solution numerically, we can approximate the integral with some Gaussian quadrature and then apply some minimization solver. The initial value for an iterative procedure is provided by (26). Let us conclude this section by considering the asymptotic analysis. It turns out that the result for G1 interpolation is similar as in the C1 case. Theorem 3. Suppose that the assumptions of Theorem 2 hold. The PH curve p that satisfies (15)–(16) and is defined by (4), (5) and (9), with parameters chosen as  π (31) λi = f (h) − f (0), i = 0, 1, φ = − , 2

M. Krajnc, V. Vitrih / Mathematics and Computers in Simulation 82 (2012) 1696–1711

1707

Fig. 2. The interpolant p (black curve) with parameters λ0 and λ1 chosen by (28) and three different solution curves (gray curves) with prescribed lengths 2.6, 3, 3.5, and λ0 = λ1 given by (29).

approximates f with the asymptotic approximation order h4 . Furthermore, the first component ω1 of the angular 1 velocity vector ω tends to zero with the order h3 , and 0 ω12 (t) dt = O(h6 ). Proof. Recall the expansions (19)–(22) from Section 4.1. The parameters λ0 and λ1 defined by (31) then expand as λi

=

√ 1 1 h − κ02 h5/2 − κ0 κ1 h7/2 + O(h9/2 ), 48 48

i = 0, 1.

Similarly as for the C1 interpolation in Section 4.1, one can see that ω1 , defined in (12), is of order h3 for φ = − π2 and that (24) holds. In this case the PH curve p is equal to ⎞ ⎛ 1 1 th − t(6t 2 − 3t + 1)κ02 h3 − t(3t 3 + 10t 2 − 9t + 2)κ0 κ1 h4 ⎟ ⎜ 24 48 ⎟ ⎜ 1 3 1 2 ⎜1 2 2 3 2 3 2 4 ⎟ + O(h5 ). ⎜ t κ0 h + t κ1 h + t ((1 − 3t )κ0 + 2(1 − 2t)(τ0 κ0 − κ2 ))h ⎟ ⎟ ⎜2 6 48 ⎠ ⎝1 1 2 3 3 4 t κ0 τ0 h + t (2t − 1)(2κ1 τ0 + κ0 τ1 )h 6 24 The distance between curves p and f is measured as a parametric distance dist(p, f ) = inf p − f ◦ Ψ , Ψ

where Ψ : [0, 1] → [0, h] is a regular reparameterization. For   1  1  Ψ (t) = th − t 2t 2 − 3t + 1 κ02 h3 − t 2t 2 − 3t + 1 κ0 κ1 h4 + O(h5 ), 24 24 curves p and f ◦ Ψ have a second order contact at t = 0 and t = 1, and the distance is  h4 dist(p, f ) ≤ 9κ02 κ12 + (3κ03 + 2τ02 κ0 − 2κ2 )2 + 4(2κ1 τ0 + κ0 τ1 )2 + O(h5 ). 768 This concludes the proof. 

1708

M. Krajnc, V. Vitrih / Mathematics and Computers in Simulation 82 (2012) 1696–1711

Table 2 Approximation error between the curve (27) and the G1 interpolant p together with values of free parameters and the value of the integral (13). h

λ0

λ1

φ

φ − φstart

1

0.8986 0.6498 0.4746 0.3436 0.2463 0.1754

0.8986 0.6498 0.4746 0.3436 0.2463 0.1754

−1.3252 −1.3368 −1.3386 −1.3389 −1.3390 −1.3390

0.0138 2.2068 × 10−3 3.4378 × 10−4 4.6762 × 10−5 6.0067 × 10−6 7.5793 × 10−7

1 2 1 4 1 8 1 16 1 32

h 1 1 2 1 4 1 8 1 16 1 32

1 0

ω12 (t)dt

4.5732 × 10−4

1.6673 × 10−5 5.2661 × 10−7 1.0625 × 10−8 1.7914 × 10−10 2.8673 × 10−12

Decay exponent

Approx. error

Decay exponent

– 4.7776 4.9846 5.6313 5.8902 5.9653

0.0329 4.9251 × 10−3 3.7743 × 10−4 2.2303 × 10−5 1.2858 × 10−6 7.6091 × 10−8

– 2.7382 3.7059 4.0809 4.1166 4.0787

To confirm the asymptotic results numerically, let us use the same example as for the C1 interpolation in Section 4.1. The lengths are computed by (28) and φ is obtained by minimizing (13) using (26). The results are shown in Table 2 and in Fig. 3. Again the numbers numerically confirm the results of Theorem 3. 5. Numerical examples Let us conclude the paper with some numerical examples, which illustrate the theoretical results obtained in previous sections. First let us compare our G1 interpolant p with parameters λ0 , λ1 chosen by (28) and φ computed by (13), with quintic PH curves having the exact RRMF frame on three different examples considered in [6, pp. 18–23]. Let P 0 = (1, 0, 0)T and let first ⎛√ √ ⎞ 2 2 ⎛ ⎞ 0 − ⎟ ⎜ 0.3106 0.5059 0.8047 2 ⎟ ⎜ 2 √ ⎟ ⎜√ ⎜ −0.3106 −0.5059 0.8047 ⎟ ⎜ 2 2 ⎟ ⎟ ⎟ , F1 = ⎜ F0 = ⎜ (32) ⎜ ⎟. 0 ⎜ 2 ⎝ 0.5059 −0.8047 −0.3106 ⎠ 2 ⎟ ⎟ ⎜ ⎟ ⎜ −1 0 ⎠ ⎝ 0

0.02 0.01 0.2

0.4

0.6

0.8

1.0

0.01 0.02 0.03

Fig. 3. In the left figure, the RMF frame of the curve (27) for h = 1 (gray) and the ERF frame of the G1 interpolant (black) are shown. The right figure shows the graph of ω1 at the optimal φ.

M. Krajnc, V. Vitrih / Mathematics and Computers in Simulation 82 (2012) 1696–1711

1709

Fig. 4. G1 quintic interpolants with the RRMF frame (top), interpolant p with the ERF frame (bottom left), and all three solution curves together (bottom right), for data (32).

The solution curves are shown in Fig. 4. On the top, both solutions having RRMF frame are presented, bottom left figure shows the interpolant p with the ERF frame, and on the bottom right, all three curves are presented together. It can be concluded that our G1 quintic is competitive with interpolants obtained in Ref. [6]. The same comparison for √ √ ⎞ ⎛ √ 3 5 5 − − ⎜ 2 ⎟ ⎛ ⎞ 10 5 ⎜ √ ⎟ 0.5833 −0.1869 0.7904 ⎜ ⎟ ⎜ ⎟ 5 ⎜ −0.6238 0.5202 0.5833 ⎟ ⎜ −0.0536 0.8928 ⎟ ⎟ ⎟ , F1 = ⎜ F0 = ⎜ (33) 5 ⎜ ⎟ ⎜ √ ⎟ −0.5202 ⎝ −0.8333 0.1869 ⎠ ⎜ ⎟ 5 ⎜ ⎟ ⎜− −0.9732 0.0536 ⎟ ⎝ 10 ⎠ is shown in Fig. 5. The conclusions are the same as for the previous example. Finally, let √ ⎛ ⎞ ⎛ √ ⎞ 1 2 1 1 3 ⎜ 2 2 2 ⎟ 0 − ⎜ ⎜ √ ⎟ √ ⎟ 2 ⎟ ⎜ 2 ⎜ ⎟ ⎜ ⎜ ⎟ 2 2⎟ ⎜ 0 ⎜− ⎟ ⎟ 1 0 0 ⎟ , F1 = ⎜ 2 F0 = ⎜ 2 ⎟ ⎜ √3 ⎜ ⎟ ⎟. √ 1 ⎟ ⎜ ⎜ ⎟ 0 2 1 ⎟ ⎜ ⎜ 1 ⎟ ⎜ ⎟ 2 ⎠ ⎝ 2 − ⎝ 2 2 2 ⎠

(34)

Fig. 5. G1 quintic interpolants with the RRMF frame (top), interpolant p with the ERF frame (bottom left), and all three solution curves together (bottom right), for data (33).

1710

M. Krajnc, V. Vitrih / Mathematics and Computers in Simulation 82 (2012) 1696–1711

Fig. 6. A G1 quintic PH interpolant together with its ERF frame for data (34).

Fig. 7. G1 Hermite splines with parameters λk0 , λk1 and φk , k = 1, 2, 3, given in Table 3.

For these data no quintic PH interpolant having the RRMF frame exists [6], while our interpolant with the ERF frame is shown in Fig. 6. As another example, let us consider a G1 spline interpolation. The data are taken from a curve g : R → R3 , g(t) = 4(log(1 + t) sin(πt), log(1 + t) cos(πt), t 2 )T , as P i := g(i/3), i = 0, 1, 2, 3. Further let Fi , i = 0, 1, 2, 3, be the Frenet frames of g at points P i . In order to construct a quintic G1 Hermite spline, which interpolates the given points and frames, we have to determine three free parameters λk0 , λk1 and φk for each polynomial segment pk , k = 1, 2, 3. Three different cases are considered. For each of them, the parameters φk are chosen in such a way that the integral (13) on every polynomial segment is minimal, but the choice of λk0 , λk1 , k = 1, 2, 3, is different for each case. Firstly, let us follow the asymptotic analysis and choose  λk0 = λk1 = P k − P k−1 , k = 1, 2, 3. The obtained G1 spline is shown in Fig. 7(left). The parameters λk0 and λk1 can be selected also in such a way that the lengths k :=  (pk ), k = 1, 2, 3, of all three polynomial segments are prescribed in advance. For the second case let us choose ( 1 ,  2 ,  3 ) = (1.6, 2.5, 3.6). The parameters λk0 and λk1 are then computed from (29). This G1 spline is shown in Fig. 7(middle). As the last case, let ( 1 ,  2 ,  3 ) = (1.3, 2.3, 3.4). The solution is presented in Fig. 7(right). The corresponding values for shape parameters are given in Table 3. It is straightforward to use the presented interpolation scheme for rigid body motions and swept surface construction. In Fig. 8(left) a motion of a rectangular parallelepiped which interpolates the given positions (black) is shown. A center of mass traverses the interpolating curve and rotations are specified by the ERF frame of the interpolant. An example of swept surface is shown in Fig. 8(right). Table 3 The values of shape parameters λk0 , λk1 and φk , k = 1, 2, 3, for all three cases.

Case 1 Case 2 Case 3

λ10 = λ11

φ1

λ20 = λ21

φ2

λ30 = λ31

φ3

1.11066 1.62246 0.69293

−0.11059 −0.12313 −0.09727

1.48977 2.28030 1.23509

−0.23447 −0.22636 −0.23066

1.82674 2.79703 1.38596

−0.74604 −0.72479 −0.74154

M. Krajnc, V. Vitrih / Mathematics and Computers in Simulation 82 (2012) 1696–1711

1711

Fig. 8. A motion of a rectangular parallelepiped (left) and an example of a swept surface (right) for the given G1 Hermite motion data.

6. Conclusion A method for computing rational motions of a rigid body has been presented. It is based on the construction of spatial quintic Pythagorean-hodograph spline curves that interpolate 3D points and corresponding frame orientations, and posses rational orthogonal frames called Euler–Rodrigues frames. A solution of a polynomial Hermite interpolation problem is expressed in a closed form with three free parameters. They are used for minimizing the length of the curve and to approach the Euler–Rodrigues frame as much as possible to the rotation minimizing one. The method can be used in several applications connected to motions of a rigid body, such as computer animation, robot manipulation, construction of a smooth camera motion, and spatial path planning in manufacturing. Another important application is a swept surface construction that attaches to the isogeometric analysis for solving partial differential equations. Similar interpolation schemes can be derived using Pythagorean-hodograph spline curves of degree n > 5. An advantage of using higher degrees is that the interpolants can have Euler–Rodrigues frames equal to rotation minimizing frames. Since the equations are highly nonlinear the analysis becomes more complicated and this is a topic for a future research. References [1] H.I. Choi, C.Y. Han, Euler–Rodrigues frames on spatial Pythagorean-hodograph curves, Computer Aided Geometric Design 19 (2002) 603–620. [2] G. Farin, J. Hoschek, M.S. Kim, Handbook of Computer Aided Geometric Design, first ed., Elsevier, Amsterdam, 2002. [3] R.T. Farouki, Pythagorean-Hodograph Curves: Algebra and Geometry Inseparable, Volume 1 of Geometry and Computing, Springer, Berlin, 2008. [4] R.T. Farouki, Quaternion and hopf map characterizations for the existence of rational rotation-minimizing frames on quintic space curves, Advances in Computational Mathematics 33 (2010) 331–348. [5] R.T. Farouki, C. Giannelli, C. Manni, A. Sestini, Identification of spatial PH quintic Hermite interpolants with near-optimal shape measures, Computer Aided Geometric Design 25 (2008) 274–297. [6] R.T. Farouki, C. Giannelli, C. Manni, A. Sestini, Design of rational rotation-minimizing rigid body motions by hermite interpolation, Mathematics of Computation 81 (2012) 879–903. [7] R.T. Farouki, C.Y. Han, Rational approximation schemes for rotation-minimizing frames on Pythagorean-hodograph curves, Computer Aided Geometrics Design 20 (2003) 435–454. [8] R.T. Farouki, M. al Kandari, T. Sakkalis, Hermite interpolation by rotation-invariant spatial Pythagorean-hodograph curves, Advances in Computational Mathematics 17 (2002) 369–383. [9] R.T. Farouki, T. Sakkalis, Rational rotation-minimizing frames on polynomial space curves of arbitrary degree, Journal of Symbolic Computation 45 (2010) 844–856. ˇ [10] G. Jakliˇc, J. Kozak, M. Krajnc, V. Vitrih, E. Zagar, An approach to geometric interpolation by Pythagorean-hodograph curves, Advances in Computational Mathematics, http://dx.doi.org/10.1007/s10444-011-9209-0, in press. [11] F. Klok, Two moving coordinate frames for sweeping along a 3D trajectory, Computer Aided Geometric Design 3 (1986) 217–229. [12] E. Kreyszig, Differential Geometry, Dover Publications Inc., New York, 1991 (Reprint of the 1963 Edition). [13] C. Mäurer, B. Jüttler, Rational approximation of rotation minimizing frames using Pythagorean-hodograph cubics, Journal of Geometry and Graphics 3 (1999) 141–159.