- Email: [email protected]
Construction of low degree rational motions
Journal of Computational and Applied Mathematics 256 (2014) 92–103
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Construction of low degree rational motions Marjeta Krajnc a,b , Karla Počkaj c , Vito Vitrih c,d,∗ a
FMF, University of Ljubljana, Jadranska 19, Ljubljana, Slovenia
b
IMFM, Jadranska 19, Ljubljana, Slovenia
c
IAM, University of Primorska, Muzejski trg 2, Koper, Slovenia
d
FAMNIT, University of Primorska, Glagoljaška 8, Koper, Slovenia
article
info
Article history: Received 10 September 2012 Received in revised form 11 July 2013 Keywords: Motion design Geometric interpolation Rational spline motion Geometric continuity
abstract Construction of rational spline motions is an important issue in robotics, animations and related fields. In this paper a geometric approach to interpolate given sequence of rigid body positions is considered, which, in contrast to standard approaches, is free of choosing parameter values in advance and it enables the lowest possible degree of the motion. A general solution to the problem how to interpolate 2n given positions by rational motion of degree 2n is presented and two particular cases, motions of degree six and eight, are studied in more detail. This interpolation scheme is further generalized to a method for constructing first order geometric continuous rational spline motions of degree six. Numerical examples are given which confirm the presented theoretical results. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Rational spline motions are important in robotics, computer graphics, animations and several related fields. They are defined by the property that the trajectories of points of a moving rigid body are piecewise rational curves. An important task is to construct a motion through given positions of a rigid body. If some standard interpolation techniques from curve design are applied, the resulting motion heavily depends on a particular parameterization, which has to be chosen in advance. Furthermore, standard algorithms imply rational motions of a relatively high degree with a large discrepancy in degrees of freedom involved in the spherical and the translational part of a motion. To construct rational motions of lower degree, geometric interpolation techniques (see e.g. [1] and references therein) that generate spline motions only from given geometric data (points, tangents, etc.) can be used instead. This approach yields at least two important advantages: an automatically chosen parameterization and a higher asymptotic approximation order. The first steps in this direction were done in [2,3]. Geometric interpolation by parabolic splines was used to construct G1 quartic rational spline motions. The generalization of this approach to cubic interpolation, which leads to G2 rational spline motions of degree six, was considered in [4]. More precisely, general conditions that determine Gk motions were derived and the problem how to interpolate two orientations and the corresponding velocity and curvature data was studied in detail. The difficulty with geometric interpolation methods is that nonlinear equations are involved and the analysis of the existence and uniqueness of the interpolant can be quite complicated. In this paper, we apply the geometric interpolation approach to interpolate a given sequence of rigid body positions (center positions and rotations) with a rational motion of the lowest possible degree. We will name such motions Lagrange motions, since in the geometric interpolation, the curves that match only the given points are called Lagrange interpolants. It is well known that Lagrange geometric interpolation problems are usually much more difficult to handle than the Hermite Gk interpolation problems (see e.g., [1,5,6]). However, we conjecture that n positions can be interpolated with a rational motion
∗
Corresponding author at: IAM, University of Primorska, Muzejski trg 2, Koper, Slovenia. Tel.: +386 5 611 75 84. E-mail addresses: [email protected] (M. Krajnc), [email protected] (K. Počkaj), [email protected] (V. Vitrih).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.07.014
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
93
of the degree n for even n, and degree n + 1 for odd n. Such motions are of a smaller degree compared to motions obtained by standard approaches (see e.g., [7–9]). The first nontrivial and practically very important case, i.e., rational motion of degree six interpolating six given positions, is studied in detail. The solution of a nonlinear system involved is given in an explicit form and is proven to be unique. Additionally, partial analysis is done for motion of degree eight where the nonlinear system becomes much more complicated. The Lagrange interpolation scheme is further generalized to a method for constructing first order geometric continuous rational spline motions of degree six, that locally interpolate four given positions and first order boundary data. The paper is organized as follows. In the next section a general construction of a rational motion is presented. In Section 3, the interpolation problem is stated and the equations are derived. Two particular cases, interpolation of six (eight) positions with rational motions of degree six (eight), are considered in the following two sections. In Section 6, a method for constructing G1 continuous rational spline motion of degree six is analyzed. The paper is concluded with some numerical examples and conclusions section, which summarizes the main results of the paper and identifies possible future investigations. 2. Rational motions Motion of a rigid body is described by a coordinate transformation between fixed and the moving coordinate system in three dimensional Euclidean space (see [7]). It consists of a translational and a rotational (spherical) part. The translational part is described by a trajectory c = (c1 , c2 , c3 )⊤ of the origin of the moving coordinate system with axes determined by columns of a rotation matrix R = R(t ) ∈ R3×3 which defines the spherical part of the motion. The trajectory of an arbitrary point p on the rigid body is given as p(t ) = c (t ) + R(t ) p. Rational spline motions are obtained by choosing rational spline functions for the elements of the matrix R and for the components of the trajectory c. The degree of the motion is the maximal degree of the functions involved. Particular motions where c is determined by R and motions satisfy some additional constraints are considered e.g., in [10–12]. A difficult part is to construct a spherical motion since the matrix R must lie in a special orthogonal group SO3 . There is a bijective correspondence between the space of rotations and the unit quaternion sphere S3 with identified antipodal points (see [7]). It is given by the kinematic mapping χ : H \ {0} → SO3 , Q = (q0 , q1 , q2 , q3 ) → χ (Q ) q20 + q21 − q22 − q23 2(q1 q2 + q0 q3 ) := 2 q0 + q21 + q22 + q23 2(q1 q3 − q0 q2 )
1
2(q1 q2 − q0 q3 ) − q21 + q22 − q23 2(q2 q3 + q0 q1 )
q20
2(q1 q3 + q0 q2 ) 2(q2 q3 − q0 q1 ) , q20 − q21 − q22 + q23
(1)
that maps any nonzero quaternion Q to a rotation matrix R = χ (Q ), where H denotes the space of quaternions. Note that all nonzero quaternions
λQ ,
λ ∈ R, λ ̸= 0,
(2)
lying on the same line passing through the origin, represent the same rotation. This equivalence relation defines a 3dimensional projective space, described by homogeneous quaternion coordinates. Spherical motions can thus be identified with two antipodal curves on the unit quaternion sphere S3 . By applying the kinematic mapping to a polynomial (spline) curve q = (q0 , q1 , q2 , q3 ) of degree n, one obtains a rational spherical motion of degree 2n. Rational functions ci should then be chosen as ci =
wi r
,
with r = q20 + q21 + q22 + q23 , i = 1, 2, 3,
(3)
where w := (w1 , w2 , w3 ) is a polynomial (spline) curve of degree ≤ 2n. 3. Geometric interpolation problem Suppose that a sequence of N positions Posi , i = 0, 1, . . . , N − 1, of a rigid body is given. The task is to construct a rational motion of a low degree that interpolates these positions. The positions are described by center positions Ci and by unit quaternions Qi ∈ H that represent the rotations. Note that the neighboring quaternions must satisfy
⟨Qi , Qi+1 ⟩ ≥ 0,
i = 0, 1, . . . , N − 2,
(4)
where ⟨·, ·⟩ is the standard inner product in R4 . In standard interpolation approaches the parameter values at which the positions are interpolated at must be prescribed in advance, which leads to rational motions of degree 2(N − 1). But, in the geometric interpolation the parameterization is used as an additional freedom, and as a result the degree of the motion can be reduced. The rotational part of the motion is determined from a quaternion curve q of degree n, which by (2) satisfies q(ti ) = λi Qi ,
λi ̸= 0, i = 0, 1, . . . , N − 1.
(5)
94
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
To require the interpolation of positions Posi in a particular order, which is usually a very important issue in practice, the parameters ti must be ordered as 0 = t0 < t1 < · · · < tN −2 < tN −1 = 1.
(6)
Similarly as for rational Bézier curves, we assume that the quaternion curves are written in a standard form, i.e., λ0 = λN −1 = 1, which can always be obtained by a bilinear reparameterization (see e.g., [13]). Moreover, by (4), the parameters λi must all be positive for
⟨q(ti ), q(ti+1 )⟩ > 0,
i = 0, 1, . . . , N − 2,
(7)
to hold. System (5) is a nonlinear system of 4N equations for 4(n + 1) unknown coefficients of q, N − 2 unknown parameters (ti )Ni=−12 and N − 2 unknowns (λi )Ni=−12 . The numbers of equations and unknowns are equal iff 4N = 4(n + 1) + 2(N − 2), which yields n = ⌈N /2⌉. This implies the following conjecture. Conjecture 1. Even number of positions, N = 2n, can be geometrically interpolated by a quaternion curve of degree n. Odd number of positions, N = 2n − 1, can be interpolated by a quaternion curve of degree n too, but with two additional degrees of freedom. Note that in case all the parameters (ti )Ni=−12 and (λi )Ni=−12 are chosen in advance, the degree of a quaternion curve is N − 1. If some of the parameters are left free (see e.g. [8]), the degree can be reduced but it is always ≥ n. Another advantage of our geometric method is the rotation invariance. More precisely, the same result is obtained if one rotates the data first and then computes the geometric interpolant or if first the interpolant is computed and then rotated. Note that this may not be true if some of the parameters are chosen in advance (see e.g. [8]). In this paper we assume that N = 2n. The existence and uniqueness of the solution of the system (5) is a difficult problem in general. The first step towards the solution is to separate the linear part of the problem from the nonlinear −2 −2 one. More precisely, once the unknown parameters t := (ti )2n and λ := (λi )2n are determined, the coefficients i=1 i=1 of q can be computed using any standard interpolation scheme (Newton, Lagrange, . . . ) componentwise. To obtain the equations for unknowns t and λ only, we apply the divided differences [tℓ , tℓ+1 , . . . , tℓ+n+1 ] , ℓ = 0, 1, . . . , n − 2, that map any polynomial of degree ≤ n to zero, to the system (5). Note that the divided difference is a linear functional and [tℓ , tℓ+1 , . . . , tℓ+n+1 ] f prescribes the leading coefficient of the interpolating polynomial of degree n + 1, that matches the function f at t = ti , i = ℓ, ℓ + 1, . . . , ℓ + n + 1. This leads to 4(n − 1) nonlinear equations 0 = [tℓ , tℓ+1 , . . . , tℓ+n+1 ] q =
ℓ+ n+1
1
i=ℓ
ω˙ ℓ,n+1 (ti )
λi Qi ,
ℓ = 0, 1, . . . , n − 2,
(8)
for 4(n − 1) unknowns t and λ, where
ωℓ,k (t ) :=
ℓ+k (t − ti ),
ω˙ ℓ,k (t ) :=
i=ℓ
d dt
ωℓ,k (t ).
Further notation will be used later on. By A[j1 ,j2 ,...,jk ] we will denote a matrix A with columns indexed by j1 , j2 , . . . , jk omitted. Moreover let An := Q0 ,
Q1 ,
···,
Q2n−1 ∈ R4,2n .
Relations (6) and (7) imply the following definition. Definition 1. The solution of the nonlinear system (8) is called admissible if it lies in
Dn := {(t , λ) : 0 = t0 < t1 < · · · < t2n−1 = 1, λi > 0, i = 0, 1, . . . , 2n − 1} . The corresponding q is called the admissible interpolant. Note that the restriction of the solution to the open set Dn makes the analysis even more complicated and usually requires some extra conditions to be satisfied. The first case to be considered is the parabolic one. But, every parabolic quaternion curve lies in a 2-dimensional subspace of H, and so should the data. In practice this condition is usually not fulfilled for 4 consecutive data. Thus, the first interesting case is the cubic one, which leads to motions of degree 6. It is analyzed in the next section where a closed form solution of the nonlinear system (8) is derived and the conditions for the solution to be admissible are stated. 4. Lagrange motion of degree six In this section, a motion of degree six that interpolates six positions Posi of a rigid body is derived. Interpolating parameters ti are determined from the rotational part of the motion. For the construction of the trajectory of the center three more degrees of freedom are left.
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
95
Consider first the rotational part of the motion. Let us define constants
αj,i := det A3 −5j,i+1] , [6
i = j, j + 1, . . . , j + 4, j = 0, 1,
that depend only on given quaternions. Note also that α0,0 = α1,5 . By applying linear functionals Qi1 , Qi2 , Qi3 ,
det ·,
1 ≤ i1 < i2 < i3 ≤ 4,
to (8), the system transforms to
(−1)i+j+1
λi λ5j α0,0 + αj,i = 0, ω˙ j,4 (ti ) ω˙ j,4 (t5j )
i = 1, 2, 3, 4, j = 0, 1.
(9)
Recall that λ0 = λ5 = 1. From (9) it is straightforward to eliminate unknowns (λi )4i=1 . We obtain four equations
ω˙ 0,4 (t0 )ω˙ 1,4 (ti ) α0,i =− =: −αi , ω˙ 1,4 (t5 )ω˙ 0,4 (ti ) α1,i
i = 1, 2, 3, 4,
for four unknowns (ti )4i=1 , which further simplify to 4
tj
= αi
1 − tj j =1
ti 1 − ti
,
i = 1, 2, 3, 4.
(10)
It is easy to see that (10) has a unique real solution
√ 3
ui
α1 α2 α3 α4
, ui := , i = 1, 2, 3, 4. 1 + ui αi Once the parameters ti are determined, the solution for λi follows from (9): ti =
λi =
4 (−1)i α0,i (ti − tj ), t1 t2 t3 t4 α0,0 j=0
i = 1, 2, 3, 4.
(11)
(12)
j̸=i
Let us summarize the results. Theorem 2. Let Qi ∈ H, i = 0, 1, . . . , 5, be a sequence of normalized quaternions, such that αj,i ̸= 0, i = j, j + 1, . . . , j + 4, j = 0, 1. Then there exists a unique cubic interpolating quaternion q, that satisfies q(0) = Q0 ,
q(1) = Q5 ,
q(ti ) = λi Qi ,
i = 1, 2, 3, 4,
for pairwise distinct ti and nonzero λi , determined by (11) and (12). Remark 3. Note that the uniqueness of the solution follows also from the known result in projective and algebraic geometry which says that there exists a unique geometric rational cubic curve, that interpolates six unordered spatial points in a general position (see [14]). However, this result does not ensure the existence of the solution if the order of the points is prescribed. Remark 4. Geometrically, conditions αj,i ̸= 0, i = 1, 2, 3, 4, j = 0, 1, are equivalent to Q0 , Q5 ̸∈ Πi1 ,i2 ,i3 , where Πi1 ,i2 ,i3 are hyperplanes through the origin, defined by
Πi1 ,i2 ,i3 := det ·, Qi1 , Qi2 , Qi3 = 0,
1 ≤ i1 < i2 < i3 ≤ 4.
The obtained solution does not necessarily lie in D3 . The conditions that imply the solution to be admissible are stated in the next theorem. Theorem 5. Suppose that the assumptions of Theorem 2 hold. Then the interpolating quaternion q is admissible iff
α1 > α2 > α3 > α4 > 0 and
α0,i > 0, α0,0
i = 1, 2, 3, 4.
Proof. The result follows directly from (11) and (12). The kinematic mapping χ defined by (1) maps a quaternion curve q = (q0 , q1 , q2 , q3 ) of degree ≤ 3 to a rotational matrix
R = χ (q) with rational coefficients of degree ≤ 6. The trajectory c of the center of the rigid body must satisfy c (ti ) = Ci ,
i = 0, 1, . . . , 5.
1 w, r
q20
(13) q21
q22
q23
Recall that c = where r = + + + and the unknown polynomial curve w must be of degree ≤ 6. Since w has 21 unknown coefficients and (13) is a linear system of only 18 equations, three more degrees of freedom are left for the construction of the translational part. If one chooses the degree of w to be ≤ 5 then c is uniquely determined.
96
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
Fig. 1. Rational motion of degree six interpolating six positions (left) and the spherical part of this motion (right). Table 1 The parametric distances between trajectory of the original and the rational motion for different values h. h
Parametric distance
Decay exponent
1
0.909483 0.0273524
/ 5.05531
0.0013242
4.36845
0.0000412
5.00488
9.7394880 × 10−7
5.40413
1.9224901 × 10−8
5.6628
3.4065213 × 10−10
5.81853
1 2 1 4 1 8 1 16 1 32 1 64
In the following example we are given a quaternion curve q
⊤ πt πt πt q(t ) = t , t + cos , sin , cos 4
4
(14)
10
and a trajectory of the center
c (t ) = (3 log(t + 1) cos(t ), 3 log(t + 1) sin(t ), 3(t + 1))⊤ .
(15)
We sample the positions
q(ti ) , Ci = c (ti ), ∥ q(ti )∥ where ti = ih, i = 0, 1, . . . , 5. Fig. 1(left) shows a rational motion of a cuboid interpolating these positions. On the right hand side of Fig. 1 the spherical part of this motion together with the trajectory R p of a cuboid center p, which is a spherical rational curve of degree 6 lying on the sphere centered at the origin with radius ∥ p∥, is shown for h = 1. Six interpolation Qi =
rotations of a cuboid (left) and interpolation positions (right) are denoted by bold cuboids. In Table 1 the parametric distances between trajectory of a point p of the original and the rational motion of degree six for different values h are shown. Moreover, the decay exponents computed as the binary logarithm of the quotient of two consecutive measurements are shown too. Since they tend to the order of approximation as h approaches zero, the last column numerically indicates that the approximation order is optimal, i.e. six. 5. Lagrange motion of degree eight In this section, the Lagrange interpolation problem for n = 4 is shortly considered. A nonlinear system (8) is now a system of 12 equations for unknown parameters t1 , t2 , . . . , t6 and unknown scalar values λ1 , λ2 , . . . , λ6 . It turns out that the analysis of equations and numerical computation of the solutions becomes much more complicated than in the cubic case. The first step towards the simplification is to separate unknown parameters ti from the unknown scalar values λi . Therefore Eqs. (8) are replaced by the equivalent ones given by
[t0 , t3 , t4 , t5 , t6 , t7 ]q = 0,
[t0 , t1 , t2 , t5 , t6 , t7 ]q = 0,
[t0 , t1 , t2 , t3 , t4 , t7 ]q = 0.
These equations can be written as 7 i=0 i̸=2j−1,2j
1
ω˙ j (ti )
λi Qi = 0,
ωj (t ) :=
ω0,7 (t ) , (t − t2j−1 )(t − t2j )
j = 1, 2, 3.
(16)
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
97
Fig. 2. Two rational spherical motions of degree 8 that interpolate 8 positions.
Let us define Ij := {1, 2, . . . , 6} \ {2j − 1, 2j} for j = 1, 2, 3. By applying linear functionals
det ·,
Qi1 , Qi2 , Qi3 ,
1 ≤ i1 < i2 < i3 ≤ 6, i1 , i2 , i3 ∈ Ij ,
(17)
to (16), the system transforms to
[2j−1,2j,i,7]
det A4
[0,2j−1,2j,7]
det A4
+ (−1)i+1 λi ω˙ j (t0 ) ω˙ j (ti ) for i ∈ Ij and j = 1, 2, 3. If we define the constants λ0
[2j−1,2j,i,7]
γj,i = (−1)i+1
det A4
[0,2j−1,2j,7]
det A4
,
δj,i = (−1)i
− λ7
[0,2j−1,2j,i]
det A4
ω˙ j (t7 )
[0,2j−1,2j,i]
[0,2j−1,2j,7]
det A4 det A4
= 0,
(18)
that depend only on the given data, and we use λ0 = λ7 = 1, then Eqs. (18) simplify to
γj,i + λi
ω˙ j (t0 ) ω˙ j (t0 ) + δ j ,i = 0, ω˙ j (ti ) ω˙ j (t7 )
i ∈ Ij , j = 1, 2, 3.
It is now straightforward to eliminate the unknowns (λi )6i=1 . Note that every equation contains only one unknown λi and every λi occurs in precisely two equations. If we eliminate λi from these two equations we end up with six equations t3 t4 (tℓ − t5 )(tℓ − t6 ) σ (2, ℓ) = t5 t6 (tℓ − t3 )(tℓ − t4 ) σ (3, ℓ), t5 t6 (tℓ − t1 )(tℓ − t2 ) σ (3, ℓ) = t1 t2 (tℓ − t5 )(tℓ − t6 ) σ (1, ℓ), t1 t2 (tℓ − t3 )(tℓ − t4 ) σ (1, ℓ) = t3 t4 (tℓ − t1 )(tℓ − t2 ) σ (2, ℓ),
ℓ = 1, 2, ℓ = 3, 4, ℓ = 5, 6,
for unknown parameters (ti )6i=1 , where
σ (j, ℓ) := γj,ℓ −
i∈Ij
ti 1 − ti
δj,ℓ ,
ℓ ∈ Ij , j = 1, 2, 3.
Although we managed to reduce the nonlinear system to only six equations, the problem is highly nonlinear and difficult to handle even numerically. Direct methods like Groebner basis are in this case computationally too expensive. Solutions can be computed by iterative methods only, but the equations are quite sensitive and require very good starting values. Clearly, more than one admissible solution can exist. As an example let us sample the positions Qi from the quaternion curve (14) at values 2i , i = 0, 1, . . . , 7. Two solutions of a nonlinear system are found numerically that yield quite similar motions. Namely, the unknown parameter values ti are equal to
(ti )6i=1 = {0.239946, 0.428939, 0.579601, 0.704642, 0.812773, 0.909874} in the first solution and to
(ti )6i=1 = {0.254982, 0.431498, 0.571357, 0.690916, 0.798808, 0.90052} in the second. Corresponding rational spherical motions of degree 8 that interpolate 8 positions are almost indistinguishable and are shown in Fig. 2. This short analysis of the case n = 4 shows that Lagrange interpolation schemes for higher degrees seem out of reach for today’s computer power at will. Therefore the usual approach for interpolating a large number of data is to use splines composed of low degree polynomial/rational segments. Construction of a rational G1 continuous spline motion of degree six is considered in the next section.
98
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
6. G 1 continuous spline motion of degree six When the number of interpolation data is very large, a general approach is to construct an interpolating spline, composed of polynomial (rational) segments of a low degree. If the presented Lagrange interpolation scheme is used for construction of spline segments, then the spline motion is only C 0 continuous. To obtain a G1 continuous spline motion that can be constructed locally, tangent directions and velocity quaternions should additionally be interpolated at breakpoints. For this purpose, the Lagrange interpolation problem from Section 3 must slightly be modified. Namely, interpolation of two center positions and two quaternions is replaced with interpolation of two tangent directions and two velocity quaternions. Unfortunately, in the geometric interpolation, the solution cannot simply be obtained by computing the limit of the solution of the Lagrange interpolation problem (see [15,16], e.g.). By replacing points with tangent information, the type of the unknowns changes. Additionally, one has to consider the positivity requirements separately. The interpolation problem is stated as follows. For a given sequence of N = 3m + 1 ≫ 6 positions Posi ∼ (Ci , Qi ), i = 0, 1, . . . , N − 1, together with prescribed tangent vectors t3ℓ−3 , t3ℓ and velocity quaternions U3ℓ−3 , U3ℓ , ℓ = 1, 2, . . . , m, the task is to construct a G1 continuous rational spline motion of degree 6, composed of m segments that interpolate these data in the geometric sense. Spherical part of the motion is defined by a cubic spline curve qS : [0, N − 1] → H with integer knots, composed of polynomial segments qℓ t ℓ := qS (u)|[ℓ−1,ℓ] ,
t ℓ := u − ℓ + 1 ∈ [0, 1], ℓ = 1, 2, . . . , m.
Similarly, cℓ t ℓ := cS (u)|[ℓ−1,ℓ]
cS : [0, N − 1] → R3 ,
is a rational spline representing the translational part. A motion is G1 continuous if a trajectory of every point is a G1 continuous curve, i.e. there exists a reparameterization, such that a reparameterized trajectory is C 1 continuous. The conditions that determine G1 continuous motions are derived in [4]. Therefrom the equations for the rotational part of the motion follow as qℓ tiℓ = λℓi Q3ℓ−3+i , i = 0, 1, 2, 3, ℓ = 1, 2, . . . , m, q˙ℓ (0) = µℓ0 Q3ℓ−3 + λℓ0 ϕ0ℓ U3ℓ−3 , q˙ℓ (1) = µℓ1 Q3ℓ + λℓ3 ϕ1ℓ U3ℓ ,
(19)
where we again assume that λℓ0 = λℓ3 = 1 and t0ℓ = 0, t3ℓ = 1. As in (5), the unknowns at every segment ℓ are the coefficients of a quaternion curve qℓ , parameters t1ℓ , t2ℓ and scalar coefficients λℓ1 , λℓ2 . Additionally, we have four new unknowns µℓ0 , µℓ1 , ϕ0ℓ and ϕ1ℓ . As explained in [4], the first two correspond to derivatives of some nonzero scalar function and the last two are the derivatives of the reparameterization function at t ℓ = 0 and t ℓ = 1. The number of unknowns is 24 m and it is equal to the number of equations. Definition 2. The solution of the nonlinear system (19) is admissible if the parameters tiℓ are ordered as 0 < t1ℓ < t2ℓ < 1,
ℓ = 1, 2, . . . , m,
(20)
and the other unknowns satisfy
λℓ1 > 0,
λℓ2 > 0,
ϕ0ℓ > 0,
ϕ1ℓ > 0,
ℓ = 1, 2, . . . , m.
(21)
Let us arrange the given data into the matrices Bℓ := U3ℓ−3 ,
Q3ℓ−3 ,
Q3ℓ−2 ,
Q3ℓ−1 ,
Q3ℓ ,
U3ℓ ∈ R4,6 ,
ℓ = 1, 2, . . . , m,
and let
βjℓ,i := det Bℓ −5j,i+1] , [6
i = j, j + 1, . . . , j + 4, j = 0, 1.
Note that β0ℓ,0 = β1ℓ,5 and assume that this determinants are nonzero. Since the problem is completely local, it is enough to consider only the case ℓ = 1. In order to simplify the notation, we will skip the superscript ℓ = 1 in all unknowns. In the first step we eliminate in (19) (for ℓ = 1) the equations for coefficients of q := q1 from the ones which involve the rest of the unknowns. Therefore we apply particular divided differences, [t0 , t1 , t2 , t3 ] ([t0 , ·]q) =
3 [t0 , ti ]q i=0
[t0 , t1 , t2 , t3 ] ([t3 , ·]q) =
ω( ˙ ti )
3 [t3 , ti ]q i=0
ω( ˙ ti )
= 0, (22)
= 0,
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
99
where ω(t ) := ω0,3 (t ). More precisely, Eqs. (22) are of the form
3 µ0 λ0 λ0 ϕ0 − U0 = 0, Q0 + ω( ˙ 0) t ω( ˙ ti ) ω( ˙ 0) i =1 i =1 i 2 2 λi µ1 λ3 ϕ1 λ3 Qi + − Q3 + U3 = 0. ( t − 1 ) ω( ˙ t ) ω( ˙ 1 ) ( t − 1 ) ω( ˙ t ) ω( ˙ 1) i i i i i =0 i=0 λi Qi + ti ω( ˙ ti )
3
(23)
Recall that λ0 = λ3 = 1. By applying linear functionals Qi1 , Qi2 , Qi3 ,
0 ≤ i1 < i2 < i3 ≤ 3,
det ·,
to (23), the system transforms to
ϕj (−1)i λi−j βj,i−j+1 + 1−j β0,0 = 0, i = 1, 2, 3, ω( ˙ j) ti (ti−j − 1)j ω( ˙ ti−j ) 3−j 1 ϕj µj βj,3j+1 + − β0,0 = 0, ω( ˙ j) ω( ˙ j) ℓ=1−j tℓ1−j (tℓ − 1)j ω( ˙ tℓ )
(24)
(25)
for j = 0, 1, where βj,i := βj1,i . From (24) it is straightforward to express ϕj in terms of t1 and t2 ,
ϕj = (−1)j+1
β0,0 βj,4−3j
t1 t2
1−2j
(1 − t1 )(1 − t2 )
,
j = 0, 1.
(26)
Further, from (25) one obtains
µj = (−1)
j
3−j
ω( ˙ j)
βj,3j+1 + 1 −j j β j,4−3j ˙ tℓ ) ℓ=1−j tℓ (1 − tℓ ) ω(
t1 t2
1−2j
(1 − t1 )(1 − t2 )
,
(27)
for j = 0, 1, and from (24) it follows
λj =
β0,j+1 (t2 − t1 )tj2 , β0,4 (1 − t3−j )
j = 1, 2.
(28)
By slightly modifying and dividing particular equations in (24) and applying (26), we obtain tj (1 − t1 )(1 − t2 ) β1,1 βj+1 = , β0,4 1 − tj t1 t2
j = 1, 2,
(29)
β
where βj := β0,j . The system (29) has precisely one real solution 1,j tj =
uj 1 + uj
,
uj :=
3
√ β0,4 3 β2 β3 , β1,1 βj+1
j = 1, 2.
(30)
Since the construction was completely local, let us summarize the results for the whole spline qS . Theorem 6. Let Q3ℓ−3+i ∈ H, i = 0, 1, 2, 3, be a sequence of normalized quaternions, and let U3ℓ−3 and U3ℓ be given velocity quaternions at every segment ℓ = 1, 2, . . . , m, such that βjℓ,i ̸= 0, i = 0, 1, 2, 3, j = 0, 1. Then there exists a unique cubic
interpolating quaternion spline qS , satisfying (19) with λℓ0 = λℓ3 = 1, where ϕjℓ , µℓj , λℓj+1 and tjℓ+1 , j = 0, 1, are determined by (26)–(28) and (30), where the superscript ℓ in all unknowns is omitted. The obtained solution is not admissible if some of the relations in (20) and (21) are not fulfilled. Theorem 7. Suppose that the assumptions of Theorem 6 hold. Then the interpolating quaternion spline qS is admissible iff
β2ℓ < β3ℓ < 0,
β0ℓ,j+2 β0ℓ,4
> 0,
(−1)j β0ℓ,0 βjℓ,4−3j
< 0, j = 0, 1, ℓ = 1, 2, . . . , m.
Proof. The result follows directly from (26), (28), (30) and the fact that the last two conditions imply β0ℓ,4 /β1ℓ,1 < 0.
100
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103 Table 2 Ten given positions, i.e. unit quaternions and center positions of the cuboid. Posi i i i i i i i i i i
=0 =1 =2 =3 =4 =5 =6 =7 =8 =9
Qi⊤
Ci⊤
(0, 0.7071, 0, 0.7071) (0.3039, 0.5188, 0.5188, 0.6078) (0.6489, 0.3244, 0.4867, 0.4867) (0.8356, 0.2129, 0.3442, 0.3714) (0.9147, 0.1715, 0.2287, 0.2858) (0.9464, 0.1625, 0.1625, 0.2271) (0.9601, 0.1600, 0.1333, 0.1867) (0.9677, 0.1522, 0.1243, 0.1580) (0.9735, 0.1369, 0.1217, 0.1369) (0.9787, 0.1173, 0.1173, 0.1208)
(1, 1, 2) (−1, 1, 3) (−3, 0, 4) (−3, −1, 5) (1, −2, 6) (6, −1, 7) (5, 1, 8) (−1, 2, 9) (−8, 1, 10)
(0, 0, 1)
Table 3 Velocity quaternions and tangent vectors determined by the method presented in [17]. Ui⊤
Posi i i i i
=0 =3 =6 =9
ti⊤
(0, 0, 1, 0) (0.12698, −0.076360, −0.14218, −0.11011) (0.00501, −0.003232, −0.00818, −0.01716) (0.01, −0.01, −0.01, −0.01)
(1, 0, 1)
(−0.81650, −1.11536, 1.11536) (0.55400, 1.00895, 0.60070)
(−5, −2, 1)
Fig. 3. Ten positions of a cuboid interpolated by the spatial G1 continuous spline motion (left) and ten rotations of a cuboid interpolated by the spherical part of the motion (right).
Let us now construct the trajectory of the origin cS . By (3), polynomials (wℓ )i , i = 1, 2, 3, of degree at most 6, have to be determined. As in the Lagrange case, these polynomials are not uniquely determined and we will restrict the degrees to 5. Therefore the polynomials (wℓ )i are uniquely defined by the conditions wℓ (tiℓ ) = rℓ (tiℓ ) C3ℓ−3+i , i = 0, 1, 2, 3, ℓ w˙ℓ (t ℓ3i ) = ϕiℓ rℓ (t3i ) t3(ℓ+i−1) + r˙ℓ (t3iℓ ) C3(ℓ+i−1) ,
i = 0, 1,
ℓ = 1, 2, . . . , m,
ℓ ℓ ℓ 2 where rℓ = i=0 (qℓ )i and ϕ0 , ϕ1 , ti , i = 0, 1, 2, 3, are the parameters which have already been determined in the spherical part of the motion. Polynomials wℓ can thus be computed by the standard Newton interpolation scheme componentwise. In case when only the positions Posi are given, tangent vectors and velocity quaternions must be estimated from the given data. This can be done, e.g., by using the method presented in [17] as shown in the next example. In Table 2, positions Posi , i = 0, 1, . . . , 9, are listed and in Table 3 the tangent vectors t3ℓ−3 , t3ℓ and the velocity quaternions Q3ℓ−3 , Q3ℓ , ℓ = 1, 2, 3, obtained by the method proposed in [17] are presented. Note that at the first and at the last position the derivative data are chosen arbitrarily. Fig. 3 (left) shows the corresponding G1 spline motion of a cuboid, and in Fig. 3 (right) the spherical part of the motion together with the trajectory R p of a cuboid center p is presented. The interpolation positions are denoted by bold cuboids.
3
7. Examples In this section, two more numerical examples are presented. In the first one the Lagrange interpolation scheme, discussed in Section 4, is applied to interpolate larger number of data, i.e., more than 6 positions. The data are sampled from the quaternion curve
q(t ) =
√
t + 1, cos
πt 5
, sin
πt 5
,t
⊤
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
101
Fig. 4. Eleven initial positions of the cuboid center and its positions after one, two and six steps of the subdivision scheme.
and are determined by Qi =
q(ti ) , ∥ q(ti )∥
where ti = ih, i = 0, 1, . . . , k, k ≫ 5. The spatial motion which interpolates these k + 1 rotation data is constructed by using the following 6-point nonlinear interpolatory subdivision scheme Q2ik+1 = Qik , 1 Q2ik+ +1
=
qki
tik,2 + tik,3 2
,
i ∈ Z,
where the initial quaternions are defined by Qi0 := Qi and qki is a cubic Lagrange interpolating quaternion, defined in Section 4, which satisfies qki (tik,j−i+2 ) = λki,j−i+2 Qjk ,
j = i − 2, i − 1, . . . , i + 3.
Fig. 4 shows the interpolation positions of the cuboid center, its positions after one, two and six steps of our subdivision scheme. In Fig. 5, the spherical motion of a cuboid which corresponds to six steps of the subdivision scheme is shown. As a second example, let us construct a G1 continuous spherical spline motion, where the input data are only the unit quaternions Qi =
q(i) , ∥ q(i)∥
i = 0, 1, . . . , 9,
(31)
with q defined in (14). Instead of applying the method presented in [17] as done in one of the previous examples, the velocity quaternions U3i−3 , i = 1, 2, 3, 4, are now estimated using local normalized quartic polynomials through five consecutive quaternions. The corresponding rotational motion is presented in Fig. 6 (left). Furthermore, in Fig. 6 (right) comparison with the standard C 1 motion based on quintic quaternion curves interpolating the same data, constructed using standard Hermite interpolation techniques, is shown. Both motions are quite similar, but of course the degree of the geometrically continuous motion is much smaller.
102
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
Fig. 5. The spherical motion of a cuboid which interpolates eleven rotation data and corresponds to six steps of the subdivision scheme.
Fig. 6. Comparison between spherical parts of a G1 motion interpolating ten given rotations (31) with the velocity quaternions estimated by using local normalized quartic polynomials (left) and a standard C 1 motion interpolating the same data (right).
8. Conclusion In this paper the geometric approach has been used to interpolate a given sequence of rigid body positions (center positions and rotations) with a rational motion of the lowest possible degree. In comparison to the classical approaches, the degree of a motion has been reduced from 4n − 2 to 2n, when 2n positions are interpolated. The first nontrivial case, interpolation of six positions in a prescribed order by a motion of degree six, has been studied in detail. Under some natural restrictions, such a motion exists and it is unique. Moreover it is given in a closed form. On the other hand, interpolation of eight positions by motions of degree eight is a highly nonlinear problem, which is difficult to solve even numerically. One can find examples where more than one admissible solutions exist. Thus the usual approach to interpolate a large number of data using splines composed of low degree segments has been considered and a construction of rational G1 continuous spline motions of degree six has been analyzed. Many interesting examples have been given to illustrate practical importance and advantages of the presented methods. To uniquely construct the translational part of a motion the degree of the numerator has been reduced by one. Instead, we could use the remaining free parameters and try to minimize some particular quantities of a motion. This is perhaps an issue which would be worth considering in the future. References [1] G. Jaklič, J. Kozak, M. Krajnc, E. Žagar, On geometric interpolation by planar parametric polynomial curves, Math. Comp. 76 (260) (2007) 1981–1993. [2] H.-P. Schröcker, B. Jüttler, Motion interpolation with Bennett biarcs, in: A. Kecskemethy, A. Müller (Eds.), Proc. Computational Kinematics, Springer, 2009, pp. 141–148. [3] B. Jüttler, M. Krajnc, E. Žagar, Geometric interpolation by quartic rational spline motions, in: Advances in Robot Kinematics: Motion in Man and Machine, Springer, New York, 2010, pp. 377–384. [4] G. Jaklič, B. Jüttler, M. Krajnc, V. Vitrih, E. Žagar, Hermite interpolation by rational Gk motions of low degree, J. Comput. Appl. Math. 240 (2013) 20–30. [5] J. Kozak, M. Krajnc, Geometric interpolation by planar cubic polynomial curves, Comput. Aided Geom. Design 24 (2) (2007) 67–78. [6] J. Kozak, E. Žagar, On geometric interpolation by polynomial curves, SIAM J. Numer. Anal. 42 (3) (2004) 953–967. [7] O. Bottema, B. Roth, Theoretical Kinematics, Dover Publications Inc., New York, 1990.
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
103
[8] B. Jüttler, Visualization of moving objects using dual quaternion curves, Comput. Graph. 18 (3) (1994) 315–326. [9] A. Gfrerrer, Rational interpolation on a hypersphere, Comput. Aided Geom. Design 16 (1) (1999) 21–37. [10] G. Jaklič, M.L. Sampoli, A. Sestini, E. Žagar, C 1 rational interpolation of spherical motions with rational rotation-minimizing directed frames, Comput. Aided Geom. Design 30 (1) (2013) 159–173. [11] R.T. Farouki, C. Giannelli, C. Manni, A. Sestini, Design of rational rotation-minimizing rigid body motions by Hermite interpolation, Math. Comp. 81 (278) (2012) 879–903. [12] R.T. Farouki, Exact rotation-minimizing frames for spatial pythagorean-hodograph curves, Graph. Models 64 (6) (2002) 382–395. [13] G. Farin, Curves and surfaces for computer-aided geometric design, in: Computer Science and Scientific Computing, fourth ed., Academic Press Inc., San Diego, CA, 1997. [14] J. Harris, Algebraic Geometry: A First Course, Springer, New York, 1995. [15] J. Kozak, M. Krajnc, Geometric interpolation by planar cubic G1 splines, BIT 47 (3) (2007) 547–563. [16] M. Krajnc, Geometric Hermite interpolation by cubic G1 splines, Nonlinear Anal. 70 (7) (2009) 2614–2626. [17] T. Horsch, B. Jüttler, Cartesian spline interpolation for industrial robots, Comput. Aided Des. 30 (1998) 217–224.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Construction of low degree rational motions Marjeta Krajnc a,b , Karla Počkaj c , Vito Vitrih c,d,∗ a
FMF, University of Ljubljana, Jadranska 19, Ljubljana, Slovenia
b
IMFM, Jadranska 19, Ljubljana, Slovenia
c
IAM, University of Primorska, Muzejski trg 2, Koper, Slovenia
d
FAMNIT, University of Primorska, Glagoljaška 8, Koper, Slovenia
article
info
Article history: Received 10 September 2012 Received in revised form 11 July 2013 Keywords: Motion design Geometric interpolation Rational spline motion Geometric continuity
abstract Construction of rational spline motions is an important issue in robotics, animations and related fields. In this paper a geometric approach to interpolate given sequence of rigid body positions is considered, which, in contrast to standard approaches, is free of choosing parameter values in advance and it enables the lowest possible degree of the motion. A general solution to the problem how to interpolate 2n given positions by rational motion of degree 2n is presented and two particular cases, motions of degree six and eight, are studied in more detail. This interpolation scheme is further generalized to a method for constructing first order geometric continuous rational spline motions of degree six. Numerical examples are given which confirm the presented theoretical results. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Rational spline motions are important in robotics, computer graphics, animations and several related fields. They are defined by the property that the trajectories of points of a moving rigid body are piecewise rational curves. An important task is to construct a motion through given positions of a rigid body. If some standard interpolation techniques from curve design are applied, the resulting motion heavily depends on a particular parameterization, which has to be chosen in advance. Furthermore, standard algorithms imply rational motions of a relatively high degree with a large discrepancy in degrees of freedom involved in the spherical and the translational part of a motion. To construct rational motions of lower degree, geometric interpolation techniques (see e.g. [1] and references therein) that generate spline motions only from given geometric data (points, tangents, etc.) can be used instead. This approach yields at least two important advantages: an automatically chosen parameterization and a higher asymptotic approximation order. The first steps in this direction were done in [2,3]. Geometric interpolation by parabolic splines was used to construct G1 quartic rational spline motions. The generalization of this approach to cubic interpolation, which leads to G2 rational spline motions of degree six, was considered in [4]. More precisely, general conditions that determine Gk motions were derived and the problem how to interpolate two orientations and the corresponding velocity and curvature data was studied in detail. The difficulty with geometric interpolation methods is that nonlinear equations are involved and the analysis of the existence and uniqueness of the interpolant can be quite complicated. In this paper, we apply the geometric interpolation approach to interpolate a given sequence of rigid body positions (center positions and rotations) with a rational motion of the lowest possible degree. We will name such motions Lagrange motions, since in the geometric interpolation, the curves that match only the given points are called Lagrange interpolants. It is well known that Lagrange geometric interpolation problems are usually much more difficult to handle than the Hermite Gk interpolation problems (see e.g., [1,5,6]). However, we conjecture that n positions can be interpolated with a rational motion
∗
Corresponding author at: IAM, University of Primorska, Muzejski trg 2, Koper, Slovenia. Tel.: +386 5 611 75 84. E-mail addresses: [email protected] (M. Krajnc), [email protected] (K. Počkaj), [email protected] (V. Vitrih).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.07.014
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
93
of the degree n for even n, and degree n + 1 for odd n. Such motions are of a smaller degree compared to motions obtained by standard approaches (see e.g., [7–9]). The first nontrivial and practically very important case, i.e., rational motion of degree six interpolating six given positions, is studied in detail. The solution of a nonlinear system involved is given in an explicit form and is proven to be unique. Additionally, partial analysis is done for motion of degree eight where the nonlinear system becomes much more complicated. The Lagrange interpolation scheme is further generalized to a method for constructing first order geometric continuous rational spline motions of degree six, that locally interpolate four given positions and first order boundary data. The paper is organized as follows. In the next section a general construction of a rational motion is presented. In Section 3, the interpolation problem is stated and the equations are derived. Two particular cases, interpolation of six (eight) positions with rational motions of degree six (eight), are considered in the following two sections. In Section 6, a method for constructing G1 continuous rational spline motion of degree six is analyzed. The paper is concluded with some numerical examples and conclusions section, which summarizes the main results of the paper and identifies possible future investigations. 2. Rational motions Motion of a rigid body is described by a coordinate transformation between fixed and the moving coordinate system in three dimensional Euclidean space (see [7]). It consists of a translational and a rotational (spherical) part. The translational part is described by a trajectory c = (c1 , c2 , c3 )⊤ of the origin of the moving coordinate system with axes determined by columns of a rotation matrix R = R(t ) ∈ R3×3 which defines the spherical part of the motion. The trajectory of an arbitrary point p on the rigid body is given as p(t ) = c (t ) + R(t ) p. Rational spline motions are obtained by choosing rational spline functions for the elements of the matrix R and for the components of the trajectory c. The degree of the motion is the maximal degree of the functions involved. Particular motions where c is determined by R and motions satisfy some additional constraints are considered e.g., in [10–12]. A difficult part is to construct a spherical motion since the matrix R must lie in a special orthogonal group SO3 . There is a bijective correspondence between the space of rotations and the unit quaternion sphere S3 with identified antipodal points (see [7]). It is given by the kinematic mapping χ : H \ {0} → SO3 , Q = (q0 , q1 , q2 , q3 ) → χ (Q ) q20 + q21 − q22 − q23 2(q1 q2 + q0 q3 ) := 2 q0 + q21 + q22 + q23 2(q1 q3 − q0 q2 )
1
2(q1 q2 − q0 q3 ) − q21 + q22 − q23 2(q2 q3 + q0 q1 )
q20
2(q1 q3 + q0 q2 ) 2(q2 q3 − q0 q1 ) , q20 − q21 − q22 + q23
(1)
that maps any nonzero quaternion Q to a rotation matrix R = χ (Q ), where H denotes the space of quaternions. Note that all nonzero quaternions
λQ ,
λ ∈ R, λ ̸= 0,
(2)
lying on the same line passing through the origin, represent the same rotation. This equivalence relation defines a 3dimensional projective space, described by homogeneous quaternion coordinates. Spherical motions can thus be identified with two antipodal curves on the unit quaternion sphere S3 . By applying the kinematic mapping to a polynomial (spline) curve q = (q0 , q1 , q2 , q3 ) of degree n, one obtains a rational spherical motion of degree 2n. Rational functions ci should then be chosen as ci =
wi r
,
with r = q20 + q21 + q22 + q23 , i = 1, 2, 3,
(3)
where w := (w1 , w2 , w3 ) is a polynomial (spline) curve of degree ≤ 2n. 3. Geometric interpolation problem Suppose that a sequence of N positions Posi , i = 0, 1, . . . , N − 1, of a rigid body is given. The task is to construct a rational motion of a low degree that interpolates these positions. The positions are described by center positions Ci and by unit quaternions Qi ∈ H that represent the rotations. Note that the neighboring quaternions must satisfy
⟨Qi , Qi+1 ⟩ ≥ 0,
i = 0, 1, . . . , N − 2,
(4)
where ⟨·, ·⟩ is the standard inner product in R4 . In standard interpolation approaches the parameter values at which the positions are interpolated at must be prescribed in advance, which leads to rational motions of degree 2(N − 1). But, in the geometric interpolation the parameterization is used as an additional freedom, and as a result the degree of the motion can be reduced. The rotational part of the motion is determined from a quaternion curve q of degree n, which by (2) satisfies q(ti ) = λi Qi ,
λi ̸= 0, i = 0, 1, . . . , N − 1.
(5)
94
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
To require the interpolation of positions Posi in a particular order, which is usually a very important issue in practice, the parameters ti must be ordered as 0 = t0 < t1 < · · · < tN −2 < tN −1 = 1.
(6)
Similarly as for rational Bézier curves, we assume that the quaternion curves are written in a standard form, i.e., λ0 = λN −1 = 1, which can always be obtained by a bilinear reparameterization (see e.g., [13]). Moreover, by (4), the parameters λi must all be positive for
⟨q(ti ), q(ti+1 )⟩ > 0,
i = 0, 1, . . . , N − 2,
(7)
to hold. System (5) is a nonlinear system of 4N equations for 4(n + 1) unknown coefficients of q, N − 2 unknown parameters (ti )Ni=−12 and N − 2 unknowns (λi )Ni=−12 . The numbers of equations and unknowns are equal iff 4N = 4(n + 1) + 2(N − 2), which yields n = ⌈N /2⌉. This implies the following conjecture. Conjecture 1. Even number of positions, N = 2n, can be geometrically interpolated by a quaternion curve of degree n. Odd number of positions, N = 2n − 1, can be interpolated by a quaternion curve of degree n too, but with two additional degrees of freedom. Note that in case all the parameters (ti )Ni=−12 and (λi )Ni=−12 are chosen in advance, the degree of a quaternion curve is N − 1. If some of the parameters are left free (see e.g. [8]), the degree can be reduced but it is always ≥ n. Another advantage of our geometric method is the rotation invariance. More precisely, the same result is obtained if one rotates the data first and then computes the geometric interpolant or if first the interpolant is computed and then rotated. Note that this may not be true if some of the parameters are chosen in advance (see e.g. [8]). In this paper we assume that N = 2n. The existence and uniqueness of the solution of the system (5) is a difficult problem in general. The first step towards the solution is to separate the linear part of the problem from the nonlinear −2 −2 one. More precisely, once the unknown parameters t := (ti )2n and λ := (λi )2n are determined, the coefficients i=1 i=1 of q can be computed using any standard interpolation scheme (Newton, Lagrange, . . . ) componentwise. To obtain the equations for unknowns t and λ only, we apply the divided differences [tℓ , tℓ+1 , . . . , tℓ+n+1 ] , ℓ = 0, 1, . . . , n − 2, that map any polynomial of degree ≤ n to zero, to the system (5). Note that the divided difference is a linear functional and [tℓ , tℓ+1 , . . . , tℓ+n+1 ] f prescribes the leading coefficient of the interpolating polynomial of degree n + 1, that matches the function f at t = ti , i = ℓ, ℓ + 1, . . . , ℓ + n + 1. This leads to 4(n − 1) nonlinear equations 0 = [tℓ , tℓ+1 , . . . , tℓ+n+1 ] q =
ℓ+ n+1
1
i=ℓ
ω˙ ℓ,n+1 (ti )
λi Qi ,
ℓ = 0, 1, . . . , n − 2,
(8)
for 4(n − 1) unknowns t and λ, where
ωℓ,k (t ) :=
ℓ+k (t − ti ),
ω˙ ℓ,k (t ) :=
i=ℓ
d dt
ωℓ,k (t ).
Further notation will be used later on. By A[j1 ,j2 ,...,jk ] we will denote a matrix A with columns indexed by j1 , j2 , . . . , jk omitted. Moreover let An := Q0 ,
Q1 ,
···,
Q2n−1 ∈ R4,2n .
Relations (6) and (7) imply the following definition. Definition 1. The solution of the nonlinear system (8) is called admissible if it lies in
Dn := {(t , λ) : 0 = t0 < t1 < · · · < t2n−1 = 1, λi > 0, i = 0, 1, . . . , 2n − 1} . The corresponding q is called the admissible interpolant. Note that the restriction of the solution to the open set Dn makes the analysis even more complicated and usually requires some extra conditions to be satisfied. The first case to be considered is the parabolic one. But, every parabolic quaternion curve lies in a 2-dimensional subspace of H, and so should the data. In practice this condition is usually not fulfilled for 4 consecutive data. Thus, the first interesting case is the cubic one, which leads to motions of degree 6. It is analyzed in the next section where a closed form solution of the nonlinear system (8) is derived and the conditions for the solution to be admissible are stated. 4. Lagrange motion of degree six In this section, a motion of degree six that interpolates six positions Posi of a rigid body is derived. Interpolating parameters ti are determined from the rotational part of the motion. For the construction of the trajectory of the center three more degrees of freedom are left.
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
95
Consider first the rotational part of the motion. Let us define constants
αj,i := det A3 −5j,i+1] , [6
i = j, j + 1, . . . , j + 4, j = 0, 1,
that depend only on given quaternions. Note also that α0,0 = α1,5 . By applying linear functionals Qi1 , Qi2 , Qi3 ,
det ·,
1 ≤ i1 < i2 < i3 ≤ 4,
to (8), the system transforms to
(−1)i+j+1
λi λ5j α0,0 + αj,i = 0, ω˙ j,4 (ti ) ω˙ j,4 (t5j )
i = 1, 2, 3, 4, j = 0, 1.
(9)
Recall that λ0 = λ5 = 1. From (9) it is straightforward to eliminate unknowns (λi )4i=1 . We obtain four equations
ω˙ 0,4 (t0 )ω˙ 1,4 (ti ) α0,i =− =: −αi , ω˙ 1,4 (t5 )ω˙ 0,4 (ti ) α1,i
i = 1, 2, 3, 4,
for four unknowns (ti )4i=1 , which further simplify to 4
tj
= αi
1 − tj j =1
ti 1 − ti
,
i = 1, 2, 3, 4.
(10)
It is easy to see that (10) has a unique real solution
√ 3
ui
α1 α2 α3 α4
, ui := , i = 1, 2, 3, 4. 1 + ui αi Once the parameters ti are determined, the solution for λi follows from (9): ti =
λi =
4 (−1)i α0,i (ti − tj ), t1 t2 t3 t4 α0,0 j=0
i = 1, 2, 3, 4.
(11)
(12)
j̸=i
Let us summarize the results. Theorem 2. Let Qi ∈ H, i = 0, 1, . . . , 5, be a sequence of normalized quaternions, such that αj,i ̸= 0, i = j, j + 1, . . . , j + 4, j = 0, 1. Then there exists a unique cubic interpolating quaternion q, that satisfies q(0) = Q0 ,
q(1) = Q5 ,
q(ti ) = λi Qi ,
i = 1, 2, 3, 4,
for pairwise distinct ti and nonzero λi , determined by (11) and (12). Remark 3. Note that the uniqueness of the solution follows also from the known result in projective and algebraic geometry which says that there exists a unique geometric rational cubic curve, that interpolates six unordered spatial points in a general position (see [14]). However, this result does not ensure the existence of the solution if the order of the points is prescribed. Remark 4. Geometrically, conditions αj,i ̸= 0, i = 1, 2, 3, 4, j = 0, 1, are equivalent to Q0 , Q5 ̸∈ Πi1 ,i2 ,i3 , where Πi1 ,i2 ,i3 are hyperplanes through the origin, defined by
Πi1 ,i2 ,i3 := det ·, Qi1 , Qi2 , Qi3 = 0,
1 ≤ i1 < i2 < i3 ≤ 4.
The obtained solution does not necessarily lie in D3 . The conditions that imply the solution to be admissible are stated in the next theorem. Theorem 5. Suppose that the assumptions of Theorem 2 hold. Then the interpolating quaternion q is admissible iff
α1 > α2 > α3 > α4 > 0 and
α0,i > 0, α0,0
i = 1, 2, 3, 4.
Proof. The result follows directly from (11) and (12). The kinematic mapping χ defined by (1) maps a quaternion curve q = (q0 , q1 , q2 , q3 ) of degree ≤ 3 to a rotational matrix
R = χ (q) with rational coefficients of degree ≤ 6. The trajectory c of the center of the rigid body must satisfy c (ti ) = Ci ,
i = 0, 1, . . . , 5.
1 w, r
q20
(13) q21
q22
q23
Recall that c = where r = + + + and the unknown polynomial curve w must be of degree ≤ 6. Since w has 21 unknown coefficients and (13) is a linear system of only 18 equations, three more degrees of freedom are left for the construction of the translational part. If one chooses the degree of w to be ≤ 5 then c is uniquely determined.
96
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
Fig. 1. Rational motion of degree six interpolating six positions (left) and the spherical part of this motion (right). Table 1 The parametric distances between trajectory of the original and the rational motion for different values h. h
Parametric distance
Decay exponent
1
0.909483 0.0273524
/ 5.05531
0.0013242
4.36845
0.0000412
5.00488
9.7394880 × 10−7
5.40413
1.9224901 × 10−8
5.6628
3.4065213 × 10−10
5.81853
1 2 1 4 1 8 1 16 1 32 1 64
In the following example we are given a quaternion curve q
⊤ πt πt πt q(t ) = t , t + cos , sin , cos 4
4
(14)
10
and a trajectory of the center
c (t ) = (3 log(t + 1) cos(t ), 3 log(t + 1) sin(t ), 3(t + 1))⊤ .
(15)
We sample the positions
q(ti ) , Ci = c (ti ), ∥ q(ti )∥ where ti = ih, i = 0, 1, . . . , 5. Fig. 1(left) shows a rational motion of a cuboid interpolating these positions. On the right hand side of Fig. 1 the spherical part of this motion together with the trajectory R p of a cuboid center p, which is a spherical rational curve of degree 6 lying on the sphere centered at the origin with radius ∥ p∥, is shown for h = 1. Six interpolation Qi =
rotations of a cuboid (left) and interpolation positions (right) are denoted by bold cuboids. In Table 1 the parametric distances between trajectory of a point p of the original and the rational motion of degree six for different values h are shown. Moreover, the decay exponents computed as the binary logarithm of the quotient of two consecutive measurements are shown too. Since they tend to the order of approximation as h approaches zero, the last column numerically indicates that the approximation order is optimal, i.e. six. 5. Lagrange motion of degree eight In this section, the Lagrange interpolation problem for n = 4 is shortly considered. A nonlinear system (8) is now a system of 12 equations for unknown parameters t1 , t2 , . . . , t6 and unknown scalar values λ1 , λ2 , . . . , λ6 . It turns out that the analysis of equations and numerical computation of the solutions becomes much more complicated than in the cubic case. The first step towards the simplification is to separate unknown parameters ti from the unknown scalar values λi . Therefore Eqs. (8) are replaced by the equivalent ones given by
[t0 , t3 , t4 , t5 , t6 , t7 ]q = 0,
[t0 , t1 , t2 , t5 , t6 , t7 ]q = 0,
[t0 , t1 , t2 , t3 , t4 , t7 ]q = 0.
These equations can be written as 7 i=0 i̸=2j−1,2j
1
ω˙ j (ti )
λi Qi = 0,
ωj (t ) :=
ω0,7 (t ) , (t − t2j−1 )(t − t2j )
j = 1, 2, 3.
(16)
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
97
Fig. 2. Two rational spherical motions of degree 8 that interpolate 8 positions.
Let us define Ij := {1, 2, . . . , 6} \ {2j − 1, 2j} for j = 1, 2, 3. By applying linear functionals
det ·,
Qi1 , Qi2 , Qi3 ,
1 ≤ i1 < i2 < i3 ≤ 6, i1 , i2 , i3 ∈ Ij ,
(17)
to (16), the system transforms to
[2j−1,2j,i,7]
det A4
[0,2j−1,2j,7]
det A4
+ (−1)i+1 λi ω˙ j (t0 ) ω˙ j (ti ) for i ∈ Ij and j = 1, 2, 3. If we define the constants λ0
[2j−1,2j,i,7]
γj,i = (−1)i+1
det A4
[0,2j−1,2j,7]
det A4
,
δj,i = (−1)i
− λ7
[0,2j−1,2j,i]
det A4
ω˙ j (t7 )
[0,2j−1,2j,i]
[0,2j−1,2j,7]
det A4 det A4
= 0,
(18)
that depend only on the given data, and we use λ0 = λ7 = 1, then Eqs. (18) simplify to
γj,i + λi
ω˙ j (t0 ) ω˙ j (t0 ) + δ j ,i = 0, ω˙ j (ti ) ω˙ j (t7 )
i ∈ Ij , j = 1, 2, 3.
It is now straightforward to eliminate the unknowns (λi )6i=1 . Note that every equation contains only one unknown λi and every λi occurs in precisely two equations. If we eliminate λi from these two equations we end up with six equations t3 t4 (tℓ − t5 )(tℓ − t6 ) σ (2, ℓ) = t5 t6 (tℓ − t3 )(tℓ − t4 ) σ (3, ℓ), t5 t6 (tℓ − t1 )(tℓ − t2 ) σ (3, ℓ) = t1 t2 (tℓ − t5 )(tℓ − t6 ) σ (1, ℓ), t1 t2 (tℓ − t3 )(tℓ − t4 ) σ (1, ℓ) = t3 t4 (tℓ − t1 )(tℓ − t2 ) σ (2, ℓ),
ℓ = 1, 2, ℓ = 3, 4, ℓ = 5, 6,
for unknown parameters (ti )6i=1 , where
σ (j, ℓ) := γj,ℓ −
i∈Ij
ti 1 − ti
δj,ℓ ,
ℓ ∈ Ij , j = 1, 2, 3.
Although we managed to reduce the nonlinear system to only six equations, the problem is highly nonlinear and difficult to handle even numerically. Direct methods like Groebner basis are in this case computationally too expensive. Solutions can be computed by iterative methods only, but the equations are quite sensitive and require very good starting values. Clearly, more than one admissible solution can exist. As an example let us sample the positions Qi from the quaternion curve (14) at values 2i , i = 0, 1, . . . , 7. Two solutions of a nonlinear system are found numerically that yield quite similar motions. Namely, the unknown parameter values ti are equal to
(ti )6i=1 = {0.239946, 0.428939, 0.579601, 0.704642, 0.812773, 0.909874} in the first solution and to
(ti )6i=1 = {0.254982, 0.431498, 0.571357, 0.690916, 0.798808, 0.90052} in the second. Corresponding rational spherical motions of degree 8 that interpolate 8 positions are almost indistinguishable and are shown in Fig. 2. This short analysis of the case n = 4 shows that Lagrange interpolation schemes for higher degrees seem out of reach for today’s computer power at will. Therefore the usual approach for interpolating a large number of data is to use splines composed of low degree polynomial/rational segments. Construction of a rational G1 continuous spline motion of degree six is considered in the next section.
98
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
6. G 1 continuous spline motion of degree six When the number of interpolation data is very large, a general approach is to construct an interpolating spline, composed of polynomial (rational) segments of a low degree. If the presented Lagrange interpolation scheme is used for construction of spline segments, then the spline motion is only C 0 continuous. To obtain a G1 continuous spline motion that can be constructed locally, tangent directions and velocity quaternions should additionally be interpolated at breakpoints. For this purpose, the Lagrange interpolation problem from Section 3 must slightly be modified. Namely, interpolation of two center positions and two quaternions is replaced with interpolation of two tangent directions and two velocity quaternions. Unfortunately, in the geometric interpolation, the solution cannot simply be obtained by computing the limit of the solution of the Lagrange interpolation problem (see [15,16], e.g.). By replacing points with tangent information, the type of the unknowns changes. Additionally, one has to consider the positivity requirements separately. The interpolation problem is stated as follows. For a given sequence of N = 3m + 1 ≫ 6 positions Posi ∼ (Ci , Qi ), i = 0, 1, . . . , N − 1, together with prescribed tangent vectors t3ℓ−3 , t3ℓ and velocity quaternions U3ℓ−3 , U3ℓ , ℓ = 1, 2, . . . , m, the task is to construct a G1 continuous rational spline motion of degree 6, composed of m segments that interpolate these data in the geometric sense. Spherical part of the motion is defined by a cubic spline curve qS : [0, N − 1] → H with integer knots, composed of polynomial segments qℓ t ℓ := qS (u)|[ℓ−1,ℓ] ,
t ℓ := u − ℓ + 1 ∈ [0, 1], ℓ = 1, 2, . . . , m.
Similarly, cℓ t ℓ := cS (u)|[ℓ−1,ℓ]
cS : [0, N − 1] → R3 ,
is a rational spline representing the translational part. A motion is G1 continuous if a trajectory of every point is a G1 continuous curve, i.e. there exists a reparameterization, such that a reparameterized trajectory is C 1 continuous. The conditions that determine G1 continuous motions are derived in [4]. Therefrom the equations for the rotational part of the motion follow as qℓ tiℓ = λℓi Q3ℓ−3+i , i = 0, 1, 2, 3, ℓ = 1, 2, . . . , m, q˙ℓ (0) = µℓ0 Q3ℓ−3 + λℓ0 ϕ0ℓ U3ℓ−3 , q˙ℓ (1) = µℓ1 Q3ℓ + λℓ3 ϕ1ℓ U3ℓ ,
(19)
where we again assume that λℓ0 = λℓ3 = 1 and t0ℓ = 0, t3ℓ = 1. As in (5), the unknowns at every segment ℓ are the coefficients of a quaternion curve qℓ , parameters t1ℓ , t2ℓ and scalar coefficients λℓ1 , λℓ2 . Additionally, we have four new unknowns µℓ0 , µℓ1 , ϕ0ℓ and ϕ1ℓ . As explained in [4], the first two correspond to derivatives of some nonzero scalar function and the last two are the derivatives of the reparameterization function at t ℓ = 0 and t ℓ = 1. The number of unknowns is 24 m and it is equal to the number of equations. Definition 2. The solution of the nonlinear system (19) is admissible if the parameters tiℓ are ordered as 0 < t1ℓ < t2ℓ < 1,
ℓ = 1, 2, . . . , m,
(20)
and the other unknowns satisfy
λℓ1 > 0,
λℓ2 > 0,
ϕ0ℓ > 0,
ϕ1ℓ > 0,
ℓ = 1, 2, . . . , m.
(21)
Let us arrange the given data into the matrices Bℓ := U3ℓ−3 ,
Q3ℓ−3 ,
Q3ℓ−2 ,
Q3ℓ−1 ,
Q3ℓ ,
U3ℓ ∈ R4,6 ,
ℓ = 1, 2, . . . , m,
and let
βjℓ,i := det Bℓ −5j,i+1] , [6
i = j, j + 1, . . . , j + 4, j = 0, 1.
Note that β0ℓ,0 = β1ℓ,5 and assume that this determinants are nonzero. Since the problem is completely local, it is enough to consider only the case ℓ = 1. In order to simplify the notation, we will skip the superscript ℓ = 1 in all unknowns. In the first step we eliminate in (19) (for ℓ = 1) the equations for coefficients of q := q1 from the ones which involve the rest of the unknowns. Therefore we apply particular divided differences, [t0 , t1 , t2 , t3 ] ([t0 , ·]q) =
3 [t0 , ti ]q i=0
[t0 , t1 , t2 , t3 ] ([t3 , ·]q) =
ω( ˙ ti )
3 [t3 , ti ]q i=0
ω( ˙ ti )
= 0, (22)
= 0,
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
99
where ω(t ) := ω0,3 (t ). More precisely, Eqs. (22) are of the form
3 µ0 λ0 λ0 ϕ0 − U0 = 0, Q0 + ω( ˙ 0) t ω( ˙ ti ) ω( ˙ 0) i =1 i =1 i 2 2 λi µ1 λ3 ϕ1 λ3 Qi + − Q3 + U3 = 0. ( t − 1 ) ω( ˙ t ) ω( ˙ 1 ) ( t − 1 ) ω( ˙ t ) ω( ˙ 1) i i i i i =0 i=0 λi Qi + ti ω( ˙ ti )
3
(23)
Recall that λ0 = λ3 = 1. By applying linear functionals Qi1 , Qi2 , Qi3 ,
0 ≤ i1 < i2 < i3 ≤ 3,
det ·,
to (23), the system transforms to
ϕj (−1)i λi−j βj,i−j+1 + 1−j β0,0 = 0, i = 1, 2, 3, ω( ˙ j) ti (ti−j − 1)j ω( ˙ ti−j ) 3−j 1 ϕj µj βj,3j+1 + − β0,0 = 0, ω( ˙ j) ω( ˙ j) ℓ=1−j tℓ1−j (tℓ − 1)j ω( ˙ tℓ )
(24)
(25)
for j = 0, 1, where βj,i := βj1,i . From (24) it is straightforward to express ϕj in terms of t1 and t2 ,
ϕj = (−1)j+1
β0,0 βj,4−3j
t1 t2
1−2j
(1 − t1 )(1 − t2 )
,
j = 0, 1.
(26)
Further, from (25) one obtains
µj = (−1)
j
3−j
ω( ˙ j)
βj,3j+1 + 1 −j j β j,4−3j ˙ tℓ ) ℓ=1−j tℓ (1 − tℓ ) ω(
t1 t2
1−2j
(1 − t1 )(1 − t2 )
,
(27)
for j = 0, 1, and from (24) it follows
λj =
β0,j+1 (t2 − t1 )tj2 , β0,4 (1 − t3−j )
j = 1, 2.
(28)
By slightly modifying and dividing particular equations in (24) and applying (26), we obtain tj (1 − t1 )(1 − t2 ) β1,1 βj+1 = , β0,4 1 − tj t1 t2
j = 1, 2,
(29)
β
where βj := β0,j . The system (29) has precisely one real solution 1,j tj =
uj 1 + uj
,
uj :=
3
√ β0,4 3 β2 β3 , β1,1 βj+1
j = 1, 2.
(30)
Since the construction was completely local, let us summarize the results for the whole spline qS . Theorem 6. Let Q3ℓ−3+i ∈ H, i = 0, 1, 2, 3, be a sequence of normalized quaternions, and let U3ℓ−3 and U3ℓ be given velocity quaternions at every segment ℓ = 1, 2, . . . , m, such that βjℓ,i ̸= 0, i = 0, 1, 2, 3, j = 0, 1. Then there exists a unique cubic
interpolating quaternion spline qS , satisfying (19) with λℓ0 = λℓ3 = 1, where ϕjℓ , µℓj , λℓj+1 and tjℓ+1 , j = 0, 1, are determined by (26)–(28) and (30), where the superscript ℓ in all unknowns is omitted. The obtained solution is not admissible if some of the relations in (20) and (21) are not fulfilled. Theorem 7. Suppose that the assumptions of Theorem 6 hold. Then the interpolating quaternion spline qS is admissible iff
β2ℓ < β3ℓ < 0,
β0ℓ,j+2 β0ℓ,4
> 0,
(−1)j β0ℓ,0 βjℓ,4−3j
< 0, j = 0, 1, ℓ = 1, 2, . . . , m.
Proof. The result follows directly from (26), (28), (30) and the fact that the last two conditions imply β0ℓ,4 /β1ℓ,1 < 0.
100
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103 Table 2 Ten given positions, i.e. unit quaternions and center positions of the cuboid. Posi i i i i i i i i i i
=0 =1 =2 =3 =4 =5 =6 =7 =8 =9
Qi⊤
Ci⊤
(0, 0.7071, 0, 0.7071) (0.3039, 0.5188, 0.5188, 0.6078) (0.6489, 0.3244, 0.4867, 0.4867) (0.8356, 0.2129, 0.3442, 0.3714) (0.9147, 0.1715, 0.2287, 0.2858) (0.9464, 0.1625, 0.1625, 0.2271) (0.9601, 0.1600, 0.1333, 0.1867) (0.9677, 0.1522, 0.1243, 0.1580) (0.9735, 0.1369, 0.1217, 0.1369) (0.9787, 0.1173, 0.1173, 0.1208)
(1, 1, 2) (−1, 1, 3) (−3, 0, 4) (−3, −1, 5) (1, −2, 6) (6, −1, 7) (5, 1, 8) (−1, 2, 9) (−8, 1, 10)
(0, 0, 1)
Table 3 Velocity quaternions and tangent vectors determined by the method presented in [17]. Ui⊤
Posi i i i i
=0 =3 =6 =9
ti⊤
(0, 0, 1, 0) (0.12698, −0.076360, −0.14218, −0.11011) (0.00501, −0.003232, −0.00818, −0.01716) (0.01, −0.01, −0.01, −0.01)
(1, 0, 1)
(−0.81650, −1.11536, 1.11536) (0.55400, 1.00895, 0.60070)
(−5, −2, 1)
Fig. 3. Ten positions of a cuboid interpolated by the spatial G1 continuous spline motion (left) and ten rotations of a cuboid interpolated by the spherical part of the motion (right).
Let us now construct the trajectory of the origin cS . By (3), polynomials (wℓ )i , i = 1, 2, 3, of degree at most 6, have to be determined. As in the Lagrange case, these polynomials are not uniquely determined and we will restrict the degrees to 5. Therefore the polynomials (wℓ )i are uniquely defined by the conditions wℓ (tiℓ ) = rℓ (tiℓ ) C3ℓ−3+i , i = 0, 1, 2, 3, ℓ w˙ℓ (t ℓ3i ) = ϕiℓ rℓ (t3i ) t3(ℓ+i−1) + r˙ℓ (t3iℓ ) C3(ℓ+i−1) ,
i = 0, 1,
ℓ = 1, 2, . . . , m,
ℓ ℓ ℓ 2 where rℓ = i=0 (qℓ )i and ϕ0 , ϕ1 , ti , i = 0, 1, 2, 3, are the parameters which have already been determined in the spherical part of the motion. Polynomials wℓ can thus be computed by the standard Newton interpolation scheme componentwise. In case when only the positions Posi are given, tangent vectors and velocity quaternions must be estimated from the given data. This can be done, e.g., by using the method presented in [17] as shown in the next example. In Table 2, positions Posi , i = 0, 1, . . . , 9, are listed and in Table 3 the tangent vectors t3ℓ−3 , t3ℓ and the velocity quaternions Q3ℓ−3 , Q3ℓ , ℓ = 1, 2, 3, obtained by the method proposed in [17] are presented. Note that at the first and at the last position the derivative data are chosen arbitrarily. Fig. 3 (left) shows the corresponding G1 spline motion of a cuboid, and in Fig. 3 (right) the spherical part of the motion together with the trajectory R p of a cuboid center p is presented. The interpolation positions are denoted by bold cuboids.
3
7. Examples In this section, two more numerical examples are presented. In the first one the Lagrange interpolation scheme, discussed in Section 4, is applied to interpolate larger number of data, i.e., more than 6 positions. The data are sampled from the quaternion curve
q(t ) =
√
t + 1, cos
πt 5
, sin
πt 5
,t
⊤
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
101
Fig. 4. Eleven initial positions of the cuboid center and its positions after one, two and six steps of the subdivision scheme.
and are determined by Qi =
q(ti ) , ∥ q(ti )∥
where ti = ih, i = 0, 1, . . . , k, k ≫ 5. The spatial motion which interpolates these k + 1 rotation data is constructed by using the following 6-point nonlinear interpolatory subdivision scheme Q2ik+1 = Qik , 1 Q2ik+ +1
=
qki
tik,2 + tik,3 2
,
i ∈ Z,
where the initial quaternions are defined by Qi0 := Qi and qki is a cubic Lagrange interpolating quaternion, defined in Section 4, which satisfies qki (tik,j−i+2 ) = λki,j−i+2 Qjk ,
j = i − 2, i − 1, . . . , i + 3.
Fig. 4 shows the interpolation positions of the cuboid center, its positions after one, two and six steps of our subdivision scheme. In Fig. 5, the spherical motion of a cuboid which corresponds to six steps of the subdivision scheme is shown. As a second example, let us construct a G1 continuous spherical spline motion, where the input data are only the unit quaternions Qi =
q(i) , ∥ q(i)∥
i = 0, 1, . . . , 9,
(31)
with q defined in (14). Instead of applying the method presented in [17] as done in one of the previous examples, the velocity quaternions U3i−3 , i = 1, 2, 3, 4, are now estimated using local normalized quartic polynomials through five consecutive quaternions. The corresponding rotational motion is presented in Fig. 6 (left). Furthermore, in Fig. 6 (right) comparison with the standard C 1 motion based on quintic quaternion curves interpolating the same data, constructed using standard Hermite interpolation techniques, is shown. Both motions are quite similar, but of course the degree of the geometrically continuous motion is much smaller.
102
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
Fig. 5. The spherical motion of a cuboid which interpolates eleven rotation data and corresponds to six steps of the subdivision scheme.
Fig. 6. Comparison between spherical parts of a G1 motion interpolating ten given rotations (31) with the velocity quaternions estimated by using local normalized quartic polynomials (left) and a standard C 1 motion interpolating the same data (right).
8. Conclusion In this paper the geometric approach has been used to interpolate a given sequence of rigid body positions (center positions and rotations) with a rational motion of the lowest possible degree. In comparison to the classical approaches, the degree of a motion has been reduced from 4n − 2 to 2n, when 2n positions are interpolated. The first nontrivial case, interpolation of six positions in a prescribed order by a motion of degree six, has been studied in detail. Under some natural restrictions, such a motion exists and it is unique. Moreover it is given in a closed form. On the other hand, interpolation of eight positions by motions of degree eight is a highly nonlinear problem, which is difficult to solve even numerically. One can find examples where more than one admissible solutions exist. Thus the usual approach to interpolate a large number of data using splines composed of low degree segments has been considered and a construction of rational G1 continuous spline motions of degree six has been analyzed. Many interesting examples have been given to illustrate practical importance and advantages of the presented methods. To uniquely construct the translational part of a motion the degree of the numerator has been reduced by one. Instead, we could use the remaining free parameters and try to minimize some particular quantities of a motion. This is perhaps an issue which would be worth considering in the future. References [1] G. Jaklič, J. Kozak, M. Krajnc, E. Žagar, On geometric interpolation by planar parametric polynomial curves, Math. Comp. 76 (260) (2007) 1981–1993. [2] H.-P. Schröcker, B. Jüttler, Motion interpolation with Bennett biarcs, in: A. Kecskemethy, A. Müller (Eds.), Proc. Computational Kinematics, Springer, 2009, pp. 141–148. [3] B. Jüttler, M. Krajnc, E. Žagar, Geometric interpolation by quartic rational spline motions, in: Advances in Robot Kinematics: Motion in Man and Machine, Springer, New York, 2010, pp. 377–384. [4] G. Jaklič, B. Jüttler, M. Krajnc, V. Vitrih, E. Žagar, Hermite interpolation by rational Gk motions of low degree, J. Comput. Appl. Math. 240 (2013) 20–30. [5] J. Kozak, M. Krajnc, Geometric interpolation by planar cubic polynomial curves, Comput. Aided Geom. Design 24 (2) (2007) 67–78. [6] J. Kozak, E. Žagar, On geometric interpolation by polynomial curves, SIAM J. Numer. Anal. 42 (3) (2004) 953–967. [7] O. Bottema, B. Roth, Theoretical Kinematics, Dover Publications Inc., New York, 1990.
M. Krajnc et al. / Journal of Computational and Applied Mathematics 256 (2014) 92–103
103
[8] B. Jüttler, Visualization of moving objects using dual quaternion curves, Comput. Graph. 18 (3) (1994) 315–326. [9] A. Gfrerrer, Rational interpolation on a hypersphere, Comput. Aided Geom. Design 16 (1) (1999) 21–37. [10] G. Jaklič, M.L. Sampoli, A. Sestini, E. Žagar, C 1 rational interpolation of spherical motions with rational rotation-minimizing directed frames, Comput. Aided Geom. Design 30 (1) (2013) 159–173. [11] R.T. Farouki, C. Giannelli, C. Manni, A. Sestini, Design of rational rotation-minimizing rigid body motions by Hermite interpolation, Math. Comp. 81 (278) (2012) 879–903. [12] R.T. Farouki, Exact rotation-minimizing frames for spatial pythagorean-hodograph curves, Graph. Models 64 (6) (2002) 382–395. [13] G. Farin, Curves and surfaces for computer-aided geometric design, in: Computer Science and Scientific Computing, fourth ed., Academic Press Inc., San Diego, CA, 1997. [14] J. Harris, Algebraic Geometry: A First Course, Springer, New York, 1995. [15] J. Kozak, M. Krajnc, Geometric interpolation by planar cubic G1 splines, BIT 47 (3) (2007) 547–563. [16] M. Krajnc, Geometric Hermite interpolation by cubic G1 splines, Nonlinear Anal. 70 (7) (2009) 2614–2626. [17] T. Horsch, B. Jüttler, Cartesian spline interpolation for industrial robots, Comput. Aided Des. 30 (1998) 217–224.