# C2-Continuous Orientation Planning for Robot End-Effector with B-Spline Curve Based on Logarithmic Quaternion

Authors: Yasong Pu; Yaoyao Shi; Xiaojun Lin; Yuan Hu; Zhishan Li
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Smooth orientation planning is beneficial for the working performance and service life of industrial robots, keeping robots from violent impacts and shocks caused by discontinuous orientation planning. Nevertheless, the popular used quaternion interpolations can hardly guarantee C2 continuity for multiorientation interpolation. Aiming at the problem, an efficient quaternion interpolation methodology based on logarithmic quaternion was proposed. Quaternions of more than two key orientations were expressed in the exponential forms of quaternion. These four-dimensional quaternions in space S3, when logarithms were taken for them, could be converted to three-dimensional points in space R3 so that B-spline interpolation could be applied freely to interpolate. The core formulas that B-spline interpolated points were mapped to quaternion were founded since B-spline interpolated point vectors were decomposed to the product of unitized forms and exponents were taken for them. The proposed methodology made B-spline curve applicable to quaternion interpolation through dimension reduction and the high-order continuity of the B-spline curve remained when B-spline interpolated points were mapped to quaternions. The function for reversely finding control points of B-spline curve with zero curvature at endpoints was derived, which helped interpolation curve become smoother and sleeker. The validity and rationality of the principle were verified by the study case. For comparison, the study case was also analyzed by the popular quaternion interpolations, Spherical Linear Interpolation (SLERP) and Spherical and Quadrangle (SQUAD). The comparison results demonstrated the proposed methodology had higher smoothness than SLERP and SQUAD and thus would provide better protection for robot end-effector from violent impacts led by unreasonable multiorientation interpolation. It should be noted that the proposed methodology can be extended to multiorientation quaternion interpolation with higher continuity than the second order.


---

## Body

## 1. Introduction

Trajectory planning of industrial robot end-effector in workspace is a fundamental and important task. Smooth and nonimpact trajectory planning is necessary for the working performance and the service life of industrial robots, otherwise unsmooth trajectory planning will bring added violent impacts and shocks to the robot degrading robot’s work. The trajectory planning of industrial robots is divided into position planning and orientation planning, while the valid and reasonable orientation planning is one of the foundational cores to ensure perfect working completion for the robot mechanical arm. Orientation planning is related to rotation arrangement and is a broad and complicated project [1]. Lots of relevant researchers put great efforts into orientation planning and lots of methodologies were suggested [2–6].The main expressions of orientation of a rigid body are Euler angle, rotation matrix, and quaternion [7]. Quaternion used to describe orientation has the following advantages: (1) avoid the phenomenon of locking universal joints when rotating; (2) higher calculation efficiency; (3) easy to provide smooth interpolation. Because of the remarkable characteristics, quaternion has been widely used in rotating applications. However, quaternion also has its disadvantages in describing orientation. Quaternion belongs to rotation space and is a four-dimensional vector. As a result, the algorithms of the three-dimensional vector in Cartesian space are not suitable for quaternion directly [8], and interpolating theories in the three-dimensional vector space cannot be directly applied to quaternion orientation interpolation. That is, high-order interpolation methods in Cartesian space cannot be directly used for multiorientation interpolation of quaternion. Though many scholars have been proposing kinds of quaternion orientation interpolations, the existing research for quaternion interpolation still needs to be developed further especially for smooth higher-order continuous interpolation between multiple orientations.The most widely used orientation interpolation of quaternion isSpherical Linear Interpolation (SLERP) raised by Shoemake. The reasoning process of SLERP is concise and uncomplicated, and the application of SLERP is simple and convenient. The formula of SLERP is as the following formula [9]:(1)SLERPQ1,Q2,t=sin1−tθsinθQ1+sintθsinθQ2t∈0,1,where Q1 is the original orientation, Q2 is the end orientation, and θ is the angle between the 4-dimesion vectors of Q1 and Q2. The angle value of θ can be gotten by the dot product cos(θ) = Q1.Q2. A series of smooth quaternion orientations between Q1 and Q2 can be obtained when the variable t changes uniformly between 0 and 1.SLERP is very suitable for quaternion interpolation between only two orientations. However, when interpolations between more than two orientations proceed, SLERP can only guarantee C0 continuity, and there would be sharp points at connecting points which connect several segments of interpolations. That means trajectory of orientation interpolation of quaternion is unsmooth at connecting points. As a result, the angular velocity and the angular acceleration would have discontinuous breakpoints which will bring violent impacts and shocks to industrial robots. Therefore, SLERP is not suitable for multiorientation interpolation for industrial robots.In order to acquire smooth quaternion interpolation between more than two orientations, Shoemake presentedSpherical and Quadrangle (SQUAD) developed on the base of SLERP. Although the derivation and application of SQUAD seem more complicated, it is still a matured algorithm and is often used for smooth quaternion interpolation between more than two orientations.The formula of SQUAD is as the following formula [10]:(2)SQUADQi,Qi+1,Si,Si+1,h=SLERPSLERPQi,Qi+1,h,SLERPSi,Si+1,h,2h1−h,i=0,1,…,n−1h∈0,1,Si=Qi∗exp−lnQi−1∗Qi−1+lnQi−1∗Qi+14,i=0,1,…,n−1.Formula (2) has two functions. The first function is the interpolation function, and the second is called auxiliary function. Q0, Q1,…, Qn are series of key orientations given in the form of quaternion, h is the only variable, Qi−1 is the inverse quaternion of Qi, and Si and Si+1 are called auxiliary quaternions which can be gotten by the second function in formula 2. If we want to obtain the values of S0 and Sn in auxiliary function, the values of Q-1 and Qn+1 are needed. However, we only have the values of Q0, Q1,…, Qn. So, we set the value of S0 equal to that of Q0, and the value of Sn equal to that of Qn. This setting has little influence on the resulting interpolation curve. After the SLERP formula is used three times, SQUAD is finished.Dam et al. proved that the formula SQUAD has C1 continuity of quaternion interpolation [10] and showed by the figure that both the angular and the angular velocity in the interpolation examples by SQUAD were continuous. However, they did not mention whether SQUAD has the C2-continuity characteristic. The subsequent case example in Section 6 of this paper shows that there are broken points at connecting nodes in the angular acceleration figure obtained by SQUAD. That means SQUAD is not C2 continuous at connecting points for multiorientation quaternion interpolation. Discontinuous angular acceleration will also bring violent impacts on industrial robots and may degrade its working performance and service life as well.In order to improve high-order smoothness of multiorientation interpolation of quaternion, some scholars introduce quaternion spline curves into quaternion interpolation. Myoung-jun Kim et al. presented a kind of quaternion spline curve with higher-order derivability [11]. They proposed the concept of cumulative basis function of quaternion, and the cumulative basis function of quaternion is taken as the exponent to perform quaternion interpolation. The operation of cumulative basis function becomes more complicated since its calculation is involved quaternion. Moreover, when quaternion B-spline curve is used for orientation interpolation, the control vertexes need be found reversely by the given key orientation quaternions. However, since the quaternion curve is generated in non-Euclidean space, the control vertexes obtained by reverse solution needs to solve nonlinear equations, so the numerical solution of the control vertexes of the quaternion B-spline curve can only be obtained approximately.Yan Xing et al. presented a class of the generalized quaternion curve of B-spline [12]. They constructed series of nonlinear core functions based on the research findings of Juhász and Róth. More paths of rotation can be attained by means of choosing different core functions. They concluded the classical unit quaternion curves of B-spline in these rotating paths and considered the classical unit quaternion curves of B-spline as special cases. And the generalized B-spline quaternion curves can achieve higher order of continuity. Applications and calculations of the research findings seem a little more complex.Although there are also other several ways to construct the high-order quaternion curve [13–18], it is still difficult to tell which one method is better than another, and there is still no perfect solution so far. Relevant researchers have been striving to study and explore in this direction. So, the purpose of this paper is to find a compact and effective method of smooth multiorientation interpolation of quaternion.


## 2. Quaternion Orientation Interpolation Principle with Spline Curve Based on Logarithmic Quaternion

### 2.1. Preliminaries of Quaternion

When a rigid body rotates purely, both the orientation of the rigid body changes and the position of the particle on the rigid body change. A rotation means a new orientation of the rigid body. Conversely, knowing the orientation of a rigid body means knowing what rotation the rigid body has made. That is to say, a rotation equals an orientation.The main expressions of orientation (i.e., rotation) are as follows: rotation matrix, Euler angle, and quaternion. The rotation matrix is composed of nine elements, which results in large data and inefficiency in calculation compared with other two methods. What is more, the rotation matrix is not convenient in smooth interpolation. The big trouble with the Euler angle is the phenomenon ofGimbal Lock. That is, there will be singular points under some conditions when the second Euler angle is the multiple of 90° [19]. When Gimbal Lock happens, movement will not follow the expected path. By comparison, quaternion has less computational complexity, has not the phenomenon of Gimbal Lock, and is convenient for smooth interpolation so that it is widely used in the representation of rotation (i.e., orientation). Although quaternion has many advantages, it also has its weaknesses. Quaternion is a four-dimensional vector and cannot be visually represented in Cartesian space. Moreover, the algorithms of three-dimensional vector operation and interpolation laws cannot be directly applied for quaternion. So, these problems bring some difficulties to quaternion orientation interpolation, especially for high-order continuous orientation interpolation. Thus, lots of research are involved in it.Quaternion can be generally defined asQ=q0+q1i+q2j+q3k, where, q0, q1, q2, and q3 are all real numbers and i, j, and k are the imaginary axes of the complex number. Quaternion can also be represented as the vector form Q=q0,q1,q2,q3=w,v, where w is called the real part of quaternion and v is called the imaginary part of quaternion. The modulus of quaternion can be expressed as Q=q02+q12+q22+q32.When the real part of a quaternionQ is zero, such quaternion is called pure quaternion. Pure quaternion can be expressed as Q=0,q1,q2,q3=0,v; when the imaginary part of the quaternion Q becomes the opposite values, such quaternion is named conjugate quaternion and is expressed in Q∗, and Q∗=q0,-q1,q2,q3=w,-v. The modulus of quaternion can also be expressed by Q=QQ∗; when the multiplication of two quaternions Q1 and Q2 is equal to [1,0,0,0], Q2 is defined as the inverse quaternion of Q1. Inverse quaternion of Q is expressed in Q−1, and Q−1=Q∗/Q2. When the modulus of a quaternion is equal to 1 (Q=1), such quaternion is called unit quaternion. Nonunit quaternion can be converted to unit quaternion by Q/Q. For a unit quaternion, its conjugate quaternion Q∗ is equal to its inverse quaternion Q−1; that is, Q∗=Q−1. Because the modulus of any unit quaternion is equal to 1 (Q=q02+q12+q22+q32=1), unit quaternion belongs to S3 field. Quaternion used to describe rotation is usually unit quaternion.Compared with the mathematical operations of real number, complex number, and vector, the mathematical operations of quaternion is more complicated. And the commutative law is not appropriate for quaternion operation. If there are two arbitrary quaternionsH (H=h0,u) and Q (Q=q0,v), the multiplication operation is as the following equation:(3)HQ=h0,u∗q0,v=h0q0−u⋅v+h0v+q0u+u×v.Notice that HQ is usually not equal to QH.In general, quaternion used for rotation is mostly unit quaternion. Furthermore, unit quaternion can be written as formula (3), which is more popular than original definition expression:(4)Q=cosθ2,sinθ2∗a,b,c=cosθ2,sinθ2∗n,where θ is the rotation angle and n=(a,b,c) is the vector of the rotation axis. In addition, n is the unit vector (n=1). Function (4) shows the rotation characteristic of quaternion; that is, a quaternion means a rotation or an orientation.


### 2.2. Deduction of Quaternion Interpolation Principle Based on Logarithmic Quaternion

As we know, the common complex numberF can be defined as F = A + Bi. The angle of a complex number can be expressed by φ=arctanB/A. The complex, whose modulus equals 1 (F=1), is called the unit complex. The unit complex can be written in the exponential form as F=eϕi. Quaternion is also a kind of the complex number, which is more complicated than the plane complex number. In the same way, unit quaternion can be expressed in exponential form similar to the common complex:(5)Q=cosθ2,sinθ2∗n=eθn=eθ∗a,b,c,where the meanings of θ and n are the same as that in function (4).When we take logarithms on both sides of equation (5), we get(6)lnQ=θn=θ∗a,b,c=x,y,z=P.As can be seen from equation (6), if we take a logarithm for a unit quaternion, the four-dimensional quaternion in S3 space can be converted to three-dimensional points in R3 space. At this time, high-order interpolation theories in Cartesian space become available, such as B-spline interpolation especially. Suppose there are start orientation and end orientation of robot end-effector, which are described as quaternion orientations Qs and Qe, respectively. If the logarithm is taken for Qs and Qe according to equation (5), we can accordingly get three-dimensional start point Ps and end point Pe in Cartesian space. Then, we can use various interpolation methods in R3 space to interpolate and get series of interpolated points marked as Pi=(xi, yi, zi) and i = 0, 1, 2, …, n. According to relevant interpolation equation, the interpolated points can be written as follows:(7)Pt=xt,yt,zt,t∈0,1.The interpolated points are three-dimensional vectors inR3 space and need to be converted to unit quaternions in S3 space. To ensure the rotation axes of the quaternions are unit vectors, the rotation axis vectors should be united. Thus, the interpolation equation should be decomposed to the product of unitized forms as follows:(8)Pt=x2t+y2t+z2t⋅xt,yt,ztx2t+y2t+z2t,t∈0,1.SubstituteP(t) in equation (8) into equation (6), and take exponents on both sides of the equation. We can get(9)Qt=ept=ex2t+y2t+z2t⋅xt,yt,ztx2t+y2t+z2t.By comparing equation (9) with equation (5), the rotation angle and rotation axis vector of interpolated unit quaternion can be obtained as follows:(10)θt=x2t+y2t+z2t,at=xtx2t+y2t+z2t,bt=ytx2t+y2t+z2t,ct=ztx2t+y2t+z2t,t∈0,1,where θ(t) and [a(t), b(t), c(t)], respectively, are the rotation angle and the rotation axis vector of interpolated unit quaternion in S3. [x(t), y(t), z(t)] is the vector of the interpolated point in R3.Finally, the real part and complex part of interpolated unit quaternion can be deduced as follows:(11)Qt=q0t=cosx2t+y2t+z2t2,q1t=sinx2t+y2t+z2t2⋅xtx2t+y2t+z2t,q2t=sinx2t+y2t+z2t2⋅ytx2t+y2t+z2t,q3t=sinx2t+y2t+z2t2⋅ztx2t+y2t+z2t,t∈0,1,where [q0(t), q1(t), q2(t), q3(t)] is interpolated unit quaternion in S3. [x(t), y(t), z(t)] is a vector of the interpolated point in R3. As shown in formula (11), the interpolated points in R3 space can be mapped to the corresponding interpolated unit quaternions in S3 space.Make a summery here. The four-dimensional unit quaternions in S3 space, when logarithms were taken for them, could be converted to three-dimensional points in R3 space. Then, three-dimensional spline curve interpolation, such high-order smooth curve interpolation as the B-spline curve, can be applied freely. When interpolation is completed in R3 space, the three-dimensional interpolated points in R3 space can be mapped to the interpolated unit quaternions in S3 space according to formula (11). Thus, quaternion interpolations are realized.It should be pointed out that the rotation angleθ in formula (10) is positive. When the rotation angle θ of the key orientation given in advance is negative, it can be dealt as follows. The rotation angle θ is taken as absolute value, and the rotation axis vector is taken as the opposite number before interpolation. That is, θ∗(a,b,c) are treated as θ∗−a,−b,−c.From equations (4) and (5), we can obtain the following equation:(12)Q=cosθ2,sinθ2∗n=eθn=ex,y,z.Equation (12) can be written as follows:(13)Qt=ext,yt,zt.In equation (13), [x(t), y(t), z(t)] are interpolated points of the spline curve in Cartesian space and can be expressed as r(t). When we take exponent for r(t) and get a new equation, the continuity of the new function is consistent with the one of r(t); that is, the continuity of the new function will not change. For example, if we choose the cubic B-spline curve to be the interpolation method, the interpolated points are C2 continuous at knot points. When exponent is taken for the function of the cubic B-spline curve and new function is obtained, continuity of new function will not change and remain second-order C2 continuity. Taking exponent is the process of mapping the spline curve in Cartesian space to quaternion space, so the corresponding quaternions also keep the corresponding continuity, and the order of continuity remains unchanged.


## 2.1. Preliminaries of Quaternion

When a rigid body rotates purely, both the orientation of the rigid body changes and the position of the particle on the rigid body change. A rotation means a new orientation of the rigid body. Conversely, knowing the orientation of a rigid body means knowing what rotation the rigid body has made. That is to say, a rotation equals an orientation.The main expressions of orientation (i.e., rotation) are as follows: rotation matrix, Euler angle, and quaternion. The rotation matrix is composed of nine elements, which results in large data and inefficiency in calculation compared with other two methods. What is more, the rotation matrix is not convenient in smooth interpolation. The big trouble with the Euler angle is the phenomenon ofGimbal Lock. That is, there will be singular points under some conditions when the second Euler angle is the multiple of 90° [19]. When Gimbal Lock happens, movement will not follow the expected path. By comparison, quaternion has less computational complexity, has not the phenomenon of Gimbal Lock, and is convenient for smooth interpolation so that it is widely used in the representation of rotation (i.e., orientation). Although quaternion has many advantages, it also has its weaknesses. Quaternion is a four-dimensional vector and cannot be visually represented in Cartesian space. Moreover, the algorithms of three-dimensional vector operation and interpolation laws cannot be directly applied for quaternion. So, these problems bring some difficulties to quaternion orientation interpolation, especially for high-order continuous orientation interpolation. Thus, lots of research are involved in it.Quaternion can be generally defined asQ=q0+q1i+q2j+q3k, where, q0, q1, q2, and q3 are all real numbers and i, j, and k are the imaginary axes of the complex number. Quaternion can also be represented as the vector form Q=q0,q1,q2,q3=w,v, where w is called the real part of quaternion and v is called the imaginary part of quaternion. The modulus of quaternion can be expressed as Q=q02+q12+q22+q32.When the real part of a quaternionQ is zero, such quaternion is called pure quaternion. Pure quaternion can be expressed as Q=0,q1,q2,q3=0,v; when the imaginary part of the quaternion Q becomes the opposite values, such quaternion is named conjugate quaternion and is expressed in Q∗, and Q∗=q0,-q1,q2,q3=w,-v. The modulus of quaternion can also be expressed by Q=QQ∗; when the multiplication of two quaternions Q1 and Q2 is equal to [1,0,0,0], Q2 is defined as the inverse quaternion of Q1. Inverse quaternion of Q is expressed in Q−1, and Q−1=Q∗/Q2. When the modulus of a quaternion is equal to 1 (Q=1), such quaternion is called unit quaternion. Nonunit quaternion can be converted to unit quaternion by Q/Q. For a unit quaternion, its conjugate quaternion Q∗ is equal to its inverse quaternion Q−1; that is, Q∗=Q−1. Because the modulus of any unit quaternion is equal to 1 (Q=q02+q12+q22+q32=1), unit quaternion belongs to S3 field. Quaternion used to describe rotation is usually unit quaternion.Compared with the mathematical operations of real number, complex number, and vector, the mathematical operations of quaternion is more complicated. And the commutative law is not appropriate for quaternion operation. If there are two arbitrary quaternionsH (H=h0,u) and Q (Q=q0,v), the multiplication operation is as the following equation:(3)HQ=h0,u∗q0,v=h0q0−u⋅v+h0v+q0u+u×v.Notice that HQ is usually not equal to QH.In general, quaternion used for rotation is mostly unit quaternion. Furthermore, unit quaternion can be written as formula (3), which is more popular than original definition expression:(4)Q=cosθ2,sinθ2∗a,b,c=cosθ2,sinθ2∗n,where θ is the rotation angle and n=(a,b,c) is the vector of the rotation axis. In addition, n is the unit vector (n=1). Function (4) shows the rotation characteristic of quaternion; that is, a quaternion means a rotation or an orientation.


## 2.2. Deduction of Quaternion Interpolation Principle Based on Logarithmic Quaternion

As we know, the common complex numberF can be defined as F = A + Bi. The angle of a complex number can be expressed by φ=arctanB/A. The complex, whose modulus equals 1 (F=1), is called the unit complex. The unit complex can be written in the exponential form as F=eϕi. Quaternion is also a kind of the complex number, which is more complicated than the plane complex number. In the same way, unit quaternion can be expressed in exponential form similar to the common complex:(5)Q=cosθ2,sinθ2∗n=eθn=eθ∗a,b,c,where the meanings of θ and n are the same as that in function (4).When we take logarithms on both sides of equation (5), we get(6)lnQ=θn=θ∗a,b,c=x,y,z=P.As can be seen from equation (6), if we take a logarithm for a unit quaternion, the four-dimensional quaternion in S3 space can be converted to three-dimensional points in R3 space. At this time, high-order interpolation theories in Cartesian space become available, such as B-spline interpolation especially. Suppose there are start orientation and end orientation of robot end-effector, which are described as quaternion orientations Qs and Qe, respectively. If the logarithm is taken for Qs and Qe according to equation (5), we can accordingly get three-dimensional start point Ps and end point Pe in Cartesian space. Then, we can use various interpolation methods in R3 space to interpolate and get series of interpolated points marked as Pi=(xi, yi, zi) and i = 0, 1, 2, …, n. According to relevant interpolation equation, the interpolated points can be written as follows:(7)Pt=xt,yt,zt,t∈0,1.The interpolated points are three-dimensional vectors inR3 space and need to be converted to unit quaternions in S3 space. To ensure the rotation axes of the quaternions are unit vectors, the rotation axis vectors should be united. Thus, the interpolation equation should be decomposed to the product of unitized forms as follows:(8)Pt=x2t+y2t+z2t⋅xt,yt,ztx2t+y2t+z2t,t∈0,1.SubstituteP(t) in equation (8) into equation (6), and take exponents on both sides of the equation. We can get(9)Qt=ept=ex2t+y2t+z2t⋅xt,yt,ztx2t+y2t+z2t.By comparing equation (9) with equation (5), the rotation angle and rotation axis vector of interpolated unit quaternion can be obtained as follows:(10)θt=x2t+y2t+z2t,at=xtx2t+y2t+z2t,bt=ytx2t+y2t+z2t,ct=ztx2t+y2t+z2t,t∈0,1,where θ(t) and [a(t), b(t), c(t)], respectively, are the rotation angle and the rotation axis vector of interpolated unit quaternion in S3. [x(t), y(t), z(t)] is the vector of the interpolated point in R3.Finally, the real part and complex part of interpolated unit quaternion can be deduced as follows:(11)Qt=q0t=cosx2t+y2t+z2t2,q1t=sinx2t+y2t+z2t2⋅xtx2t+y2t+z2t,q2t=sinx2t+y2t+z2t2⋅ytx2t+y2t+z2t,q3t=sinx2t+y2t+z2t2⋅ztx2t+y2t+z2t,t∈0,1,where [q0(t), q1(t), q2(t), q3(t)] is interpolated unit quaternion in S3. [x(t), y(t), z(t)] is a vector of the interpolated point in R3. As shown in formula (11), the interpolated points in R3 space can be mapped to the corresponding interpolated unit quaternions in S3 space.Make a summery here. The four-dimensional unit quaternions in S3 space, when logarithms were taken for them, could be converted to three-dimensional points in R3 space. Then, three-dimensional spline curve interpolation, such high-order smooth curve interpolation as the B-spline curve, can be applied freely. When interpolation is completed in R3 space, the three-dimensional interpolated points in R3 space can be mapped to the interpolated unit quaternions in S3 space according to formula (11). Thus, quaternion interpolations are realized.It should be pointed out that the rotation angleθ in formula (10) is positive. When the rotation angle θ of the key orientation given in advance is negative, it can be dealt as follows. The rotation angle θ is taken as absolute value, and the rotation axis vector is taken as the opposite number before interpolation. That is, θ∗(a,b,c) are treated as θ∗−a,−b,−c.From equations (4) and (5), we can obtain the following equation:(12)Q=cosθ2,sinθ2∗n=eθn=ex,y,z.Equation (12) can be written as follows:(13)Qt=ext,yt,zt.In equation (13), [x(t), y(t), z(t)] are interpolated points of the spline curve in Cartesian space and can be expressed as r(t). When we take exponent for r(t) and get a new equation, the continuity of the new function is consistent with the one of r(t); that is, the continuity of the new function will not change. For example, if we choose the cubic B-spline curve to be the interpolation method, the interpolated points are C2 continuous at knot points. When exponent is taken for the function of the cubic B-spline curve and new function is obtained, continuity of new function will not change and remain second-order C2 continuity. Taking exponent is the process of mapping the spline curve in Cartesian space to quaternion space, so the corresponding quaternions also keep the corresponding continuity, and the order of continuity remains unchanged.


## 3. B-Spline Curve Interpolation with Zero Curvature at Endpoints

### 3.1. Preliminaries of B-Spline Curve

The B-spline curve plays an important role in the freeform curve and surface modeling design and is a significant branch for smooth interpolation. The B-spline curve is developed on the basis of the Bezier curve. It retains the advantages of the Bezier curve and overcomes the defect that the Bezier curve has not the ability of local modification. Furthermore, it has advantages in dealing with complicated shapes. Because of its smooth interpolation and convenience, the B-spline curve has been one of the most popular methods of mathematical description of freeform curve and surface.The mathematical definition of the B-spline interpolation curve is as follows [20]:(14)rt=∑i=0nNi,kudi,where di are the control points and Ni,k (u) are the basis functions. The label k is the power of basis function. The label i is the sequence label. Basis function of de Boor Cox is commonly used and is defined as follows:(15)Ni,0u=1ui≤u≤ui+10other cases,Ni,ku=u−uiui+k−uiNi,k−1u+ui+k+1−uui+k+1−ui+1Ni+1,k−1u,Regulation:00=0,where ui are the nodes of the basis function.The Cubic uniform B-spline curve is most widely used in engineering among kinds of B-spline curves because its expression is not very complicated and it has the properties of high smoothness and C2 continuity, which can meet the needs of most engineering. At the same time, the cubic uniform B-spline also has the advantages of low power, simple calculation, fast speed, and the ability of constructing various shapes. As a result, numerous engineers and designers are favored of the cubic uniform B-spline curve when modeling freeform curve and surface.The basis function of the cubic uniform B-spline curve can be described further in its effective regions of basis function nodes as follows [21]:(16)Ni,3u=N0,3u=16−u3+3u2−3u+1,N1,3u=163u3−6u2+4,N2,3u=16−3u3+3u2+3u+1,N3,3u=16u3,0≤u≤1.By synthesizing formula (14) and (16), practical formula of the cubic uniform B-spline curve can be obtained as the following formula [21]. Unless the other is specified, the following B-spline curves all refer to cubic uniform B-spline curves:(17)riu=∑j=03Nj,3udi+j=161uu2u31410-30303-630-13-31didi+1di+2di+3,u∈0,1,where di+j are the control points, Nj,3(u) are the basis functions, and u is the variable.Formula (17) is really practical for usage of the cubic uniform B-spline curve. The formula is easy to understand and convenient to use. As can be seen from formula (17), every four control points hold a B-spline curve. For example, if there are five control points, the two B-spline curves can be obtained, and so on. B-spline curves are C2 continuous at knots between every two curves [22]. It should be pointed out that no matter how many control points there are, the basis functions are the same four functions as far as the cubic uniform B-spline curve is concerned.As described above, the B-spline curve can be obtained if control points are given. However, in many engineering problems, we are often not given the control points but the knots which the spline curve would pass through. Knots are the key points belonging to the B-spline curve while control points are not a part of the curve so that we need to find control points for the B-spline curve. The process of finding control points by the use of knots is called inverse solution.For the cubic uniform B-spline curve, the formula of inverse solution of control points is expressed as follows [21]:(18)16di-1+4di+di+1=Pi,i=0,1,…,n,where di-1, di, and di+1 are the control points and Pi are the key points called as knots through which the B-spline curve must pass.In function (18), the number of unknown control points is more than the number of equations, so there are numerous solutions for the equations. In order to obtain unique solutions of this function, we need to add two constraint equations. There are several common methods to add constraint equations. We need to search a better method suitable for the characteristic of quaternion interpolation.


### 3.2. The Unique Solution of Control Points with Zero Curvature at Endpoints

For the purpose that the interpolation curves are prevented from changing violently, we could consider that the curvatures at the beginning knot and the terminal knot of the interpolation curve are all set zero. We need add two constrain equations when the curvatures at the beginning and the terminal knots of the B-spline curve are all set to zero and can obtain the unique set of solutions of control points.When the first and second control points are coincident, the curvature at the beginning knot of the B-spline curve will be zero, and the equation as below should hold [21]:(19)diPi⟶=16×didi+1⟶,where di−1, di, and di+1 are the first, the second, and the third control points, respectively. Control points di−1 and di are coincident. Pi is the beginning knot of the B-spline curve. These relationships are clearly shown in Figure 1.Figure 1
Cubic B-spline curve with free endpoint conditions.From equation (19), we can deduce one of the constraint equations as follows:(20)56di+16di+1=Pi.Similarly, when the last and the penultimate control points are coincident, the curvature at the terminal knot of the B-spline curve will be zero, and the equation as below should hold:(21)dnPn⟶=16×dndn−1⟶,where dn-1, dn, and dn+1 are the last control point, the penultimate control point, and the third control from the bottom, respectively. Control points dn+1 and dn are coincident. Pn is the terminal knot of the B-spline curve.From equation (21), we can also deduce the other constrain equation as follows;(22)16dn−1+56dn=Pn.By synthesizing equations (18), (20), and (22), we can finally obtain the following formula, which can be used to solve control points inversely by key knots of the B-spline curve under the condition that the curvatures at the beginning and the terminal knots of the B-spline curve are zero:(23)1651000⋯014100⋯001410⋯0001410⋯⋮00⋯0141000⋯015d1d2d3d4⋮dn−1dn=P1P2P3P4⋮Pn−1Pn.In formula (23), d1, d2,…, dn are a series of control points. The control point d1 has a coincidence point d0, which means the first control point has two identical control points d1 and d0. Similarly, the control point dn has a coincidence point dn+1, and the last control point has two identical control points dn and dn+1. The key knots P1, P2,…,Pn are given beforehand and must be passed through by the B-spline interpolation curve. The control points are the unknown and need to be solved. We can see from function (23) that the quantity of the unknown control points is equal to the equation quantity of the matrix so that function (23) has a unique set of solutions of control points. Then, we can interpolate between the key knots with cubic uniform B-spline curve after the values of control points are found.


## 3.1. Preliminaries of B-Spline Curve

The B-spline curve plays an important role in the freeform curve and surface modeling design and is a significant branch for smooth interpolation. The B-spline curve is developed on the basis of the Bezier curve. It retains the advantages of the Bezier curve and overcomes the defect that the Bezier curve has not the ability of local modification. Furthermore, it has advantages in dealing with complicated shapes. Because of its smooth interpolation and convenience, the B-spline curve has been one of the most popular methods of mathematical description of freeform curve and surface.The mathematical definition of the B-spline interpolation curve is as follows [20]:(14)rt=∑i=0nNi,kudi,where di are the control points and Ni,k (u) are the basis functions. The label k is the power of basis function. The label i is the sequence label. Basis function of de Boor Cox is commonly used and is defined as follows:(15)Ni,0u=1ui≤u≤ui+10other cases,Ni,ku=u−uiui+k−uiNi,k−1u+ui+k+1−uui+k+1−ui+1Ni+1,k−1u,Regulation:00=0,where ui are the nodes of the basis function.The Cubic uniform B-spline curve is most widely used in engineering among kinds of B-spline curves because its expression is not very complicated and it has the properties of high smoothness and C2 continuity, which can meet the needs of most engineering. At the same time, the cubic uniform B-spline also has the advantages of low power, simple calculation, fast speed, and the ability of constructing various shapes. As a result, numerous engineers and designers are favored of the cubic uniform B-spline curve when modeling freeform curve and surface.The basis function of the cubic uniform B-spline curve can be described further in its effective regions of basis function nodes as follows [21]:(16)Ni,3u=N0,3u=16−u3+3u2−3u+1,N1,3u=163u3−6u2+4,N2,3u=16−3u3+3u2+3u+1,N3,3u=16u3,0≤u≤1.By synthesizing formula (14) and (16), practical formula of the cubic uniform B-spline curve can be obtained as the following formula [21]. Unless the other is specified, the following B-spline curves all refer to cubic uniform B-spline curves:(17)riu=∑j=03Nj,3udi+j=161uu2u31410-30303-630-13-31didi+1di+2di+3,u∈0,1,where di+j are the control points, Nj,3(u) are the basis functions, and u is the variable.Formula (17) is really practical for usage of the cubic uniform B-spline curve. The formula is easy to understand and convenient to use. As can be seen from formula (17), every four control points hold a B-spline curve. For example, if there are five control points, the two B-spline curves can be obtained, and so on. B-spline curves are C2 continuous at knots between every two curves [22]. It should be pointed out that no matter how many control points there are, the basis functions are the same four functions as far as the cubic uniform B-spline curve is concerned.As described above, the B-spline curve can be obtained if control points are given. However, in many engineering problems, we are often not given the control points but the knots which the spline curve would pass through. Knots are the key points belonging to the B-spline curve while control points are not a part of the curve so that we need to find control points for the B-spline curve. The process of finding control points by the use of knots is called inverse solution.For the cubic uniform B-spline curve, the formula of inverse solution of control points is expressed as follows [21]:(18)16di-1+4di+di+1=Pi,i=0,1,…,n,where di-1, di, and di+1 are the control points and Pi are the key points called as knots through which the B-spline curve must pass.In function (18), the number of unknown control points is more than the number of equations, so there are numerous solutions for the equations. In order to obtain unique solutions of this function, we need to add two constraint equations. There are several common methods to add constraint equations. We need to search a better method suitable for the characteristic of quaternion interpolation.


## 3.2. The Unique Solution of Control Points with Zero Curvature at Endpoints

For the purpose that the interpolation curves are prevented from changing violently, we could consider that the curvatures at the beginning knot and the terminal knot of the interpolation curve are all set zero. We need add two constrain equations when the curvatures at the beginning and the terminal knots of the B-spline curve are all set to zero and can obtain the unique set of solutions of control points.When the first and second control points are coincident, the curvature at the beginning knot of the B-spline curve will be zero, and the equation as below should hold [21]:(19)diPi⟶=16×didi+1⟶,where di−1, di, and di+1 are the first, the second, and the third control points, respectively. Control points di−1 and di are coincident. Pi is the beginning knot of the B-spline curve. These relationships are clearly shown in Figure 1.Figure 1
Cubic B-spline curve with free endpoint conditions.From equation (19), we can deduce one of the constraint equations as follows:(20)56di+16di+1=Pi.Similarly, when the last and the penultimate control points are coincident, the curvature at the terminal knot of the B-spline curve will be zero, and the equation as below should hold:(21)dnPn⟶=16×dndn−1⟶,where dn-1, dn, and dn+1 are the last control point, the penultimate control point, and the third control from the bottom, respectively. Control points dn+1 and dn are coincident. Pn is the terminal knot of the B-spline curve.From equation (21), we can also deduce the other constrain equation as follows;(22)16dn−1+56dn=Pn.By synthesizing equations (18), (20), and (22), we can finally obtain the following formula, which can be used to solve control points inversely by key knots of the B-spline curve under the condition that the curvatures at the beginning and the terminal knots of the B-spline curve are zero:(23)1651000⋯014100⋯001410⋯0001410⋯⋮00⋯0141000⋯015d1d2d3d4⋮dn−1dn=P1P2P3P4⋮Pn−1Pn.In formula (23), d1, d2,…, dn are a series of control points. The control point d1 has a coincidence point d0, which means the first control point has two identical control points d1 and d0. Similarly, the control point dn has a coincidence point dn+1, and the last control point has two identical control points dn and dn+1. The key knots P1, P2,…,Pn are given beforehand and must be passed through by the B-spline interpolation curve. The control points are the unknown and need to be solved. We can see from function (23) that the quantity of the unknown control points is equal to the equation quantity of the matrix so that function (23) has a unique set of solutions of control points. Then, we can interpolate between the key knots with cubic uniform B-spline curve after the values of control points are found.


## 4. Steps of Orientation Interpolation with B-Spline Curve Based on Logarithmic Quaternion

We sum up the steps of orientation interpolation with the B-spline curve based on logarithmic quaternion as follows:Step 1: depending on the tasks of industrial robots, several key orientations of end-effector are obtained. The key orientations are expressed by unit quaternionsQi, and rotation angles θi and rotation axes ni can be gotten as well.Step 2: according to functions (5) and (6), the four-dimensional unit quaternions Qi are expressed in the exponential forms of the complex number and are transformed to three-dimensional vectors Pi in R3 space when logarithms are taken for these quaternions. Then, interpolation theory of the B-spline curve becomes available.Step 3: the obtained pointsPi are considered as the key knots of B-spline interpolation curves, and we can find the control points di of the B-spline curve with zero curvature at end knot points by resolving function (23).Step 4: according to the obtained control pointsdi, we can use formula (17) to apply the B-spline curve to interpolate between every two key knots Pi, and get series of interpolated points P(t) in R3.Step 5: all the interpolated pointsP(t) of the B-spline curve in R3 space are mapped to interpolated quaternions in S3 space via formula (11). Quaternion interpolation ends.The interpolation process can be demonstrated clearly by the following flowchart in Figure2.Figure 2
Flowchart of quaternion interpolation with B-spline based on logarithmic quaternion.


## 5. Case Study

In order to verify the correctness and rationality of the proposed methodology, the following study case is employed. The usage of the proposed methodology is further described in detail during the study case analysis.Case 1: In order to accomplish a certain task for an industrial robot, the end-effector has to experience four necessary orientations. The first orientation is to rotate 60° around theZ axis in the fixed frame. The second orientation is to rotate 45° around the Y axis in the fixed frame. The third orientation is to rotate 90° around the X axis in the fixed frame. The last orientation is to rotate 90° around the axis whose vector is n=(2,3,5) in the fixed frame. Using the interpolation methodology proposed in this paper, complete quaternion orientation interpolation with guarantee of C2 continuity.Case analysis: there are four key orientations given in the case as seen in Figure3. The multiorientation planning has three interpolation segments: the first segment between the 1st key orientation and the 2nd key orientation, the second segment between the 2nd key orientation and the 3rd key orientation, and the third segment between the 3rd key orientation and the 4th key orientation. The orientation interpolation process is as below. It should be pointed out that the B-spline curves used in the case are the cubic uniform B-spline curve.(1)
According to the four necessary orientations of industrial robots, we can get four key quaternionsQi and then get the corresponding rotation angles θi and axes ni, as seen in Table 1.(2)
According to formula (6), we can calculate the key knots Pi of the cubic uniform B-spline curve in R3 Cartesian space. Depending on these key knots Pi obtained, we can find the resolutions of control points di by the use of formula (23), which ensures that the curvatures at the beginning and the terminal points of the B-spline curve are zero to avoid the curve changing violently. The values of key knots Pi and control points di are filled in Table 2.(3)
Using the control pointsdi obtained, interpolations can be carried out between ever two key knots with cubic B-spline curves. According to B-spline interpolation formula (17), the interpolation functions of cubic B-spline curves can be obtained as below.The first B-spine curve function:(24)r1u=x1u=0.791u3−0.396u2−0.396u,y1u=−0.914u3+0.850u2+0.850u,z1u=0.510u3−0.779u2−0.779u+1.04,7,u∈0,1.The second B-spine curve function:(25)r2u=x2u=−1.595u3+1.978u2+1.187u,y2u=1.299u3−1.892u2−0.193u+0.785z2u=0.053u3+0.752u2−0.805u,,u∈0,1.The third B-spine curve function:(26)r3u=x3u=1.343u3−2.805u2+0.360u+1.571,y3u=−0.989u3+2.005u2−0.079u,z3u=−0.599u3+0.911u2+0.858u,u∈0,1.When the variableu changes within the range [0,1], series of interpolated points of cubic B-spline curves can be obtained. Connecting these interpolated points, we can get three cubic B-spline curves which must pass through the key knots P1, P2, P3, and P4, as can be seen in Figure 4.It is obvious that curver1(u), curve r2(u), and curve r3(u) are C3 continuous, respectively, according the above functions. It can also be verified that the whole spline composed of these three B-spline curves is C2 continuous at the key knots P2 and P3. As seen in function (27), the second derivative of r1(u) with u = 1 is equal to the second derivative of r2(u) with u = 0. That is, the whole spline is C2 continuous at the knots P2. Similarly, the second derivative of r2(u) with u = 1 is equal to the second derivative of r3(u) with u = 0 in function (28), and the whole spline is C2 continuous at the knots P3. Thus, the whole B-spline r(u) is C2 continuous:(27)d2r1udu2u=1=d2r2udu2u=0=3.957−3.7841.505,(28)d2r2udu2u=1=d2r3udu2u=0=−5.6104.0111.822.According to equation (5), the interpolated quaternion can be marked as the below equation to explain the continuity of interpolation quaternion clearly:(29)Qu=eru.In function (29), on the left side of the equation is the interpolated quaternions. On the right side of the equation is the exponential form of interpolated B-spine. Exponent is taken for r(u), and the new function obtained would have the same order continuity as r(u). Because r(u) in the case study described above is C2 continuous, especially at the key knots P2 and P3 showed above, exp[r(u)] is surely also C2 continuous at same knots points as checked by following equations:(30)d2er1udu2u=1=d2er2udu2u=0=5.366−8.2172.153,(31)d2er2udu2u=1=d2er3udu2u=0=−26.3644.0172.559.The right side of equation (29), exp[r(u)], is C2 continuous. Then, the left side of equation (29), Q(u), must be C2 continuous because the continuity on the left and right sides of the equation is the same. So, the interpolated quaternions are C2 continuous, which is demonstrated with detail in Section 6.(4)
According to formula (11), these interpolated points of B-spline curves in R3 Cartesian space can be mapped to the unit quaternions in S3 space. Thus, we eventually complete the orientation interpolations of unit quaternion. Taking u = 0, 0.2, 0.5, 0.8, and 1 partly as examples, the relevant data of the quaternion interpolations are listed in Tables 3–5, respectively.Figure 3
The four key orientations in working.Figure 4
Interpolated curve with cubic B-spline.As can been seen in Table3, when the variable u is equal to 0 and 1, respectively, the interpolated quaternion is [0.866, 0, 0, 0.5] and [0.924, 0, 0.383, 0] correspondingly, which are precisely the values of the quaternions corresponding to the first and second key orientations.Table 1
Initial information of the four key orientations.
No.QuaternionQiRotation angleθiRotation axisni1[0.866, 0, 0, 0.500]π/3(0, 0, 1)2[0.924, 0, 0.383, 0]π/4(0, 1, 0)3[0.707, 0.707, 0, 0]π/2(1, 0, 0)4[0.707, 0.211, 0.422, 0.527]π/2(0.298, 0.596, 0.745)Table 2
Control points founded by the key knots.
Key knot (Pi)Control point (di)P1 (0, 0, 1.047)d0 (0.132, −0.283, 0.307)d1 (0.132, −0.283, 0.307)P2 (0, 0.785, 0)d2 (−0.659, 1.416, −0.251)P3 (1.571, 0, 0)d3 (2.506, −0.668, −0.304)P4 (0.468, 0.937, 1.171)d5 (0.061, 1.258, 1.466)d5 (0.061, 1.258, 1.466)Table 3
The first B-spline curve is mapped to unit quaternions (data partly).
SNVariablePosition of B-spline in Cartesian spaceRotation angleVector of rotation axisUnit quaternionReal partImaginary partuxyzθabcq0q1q2q310001.047`1.0470010.866000.520.2−0.0790.1750.8850.905−0.0880.1930.9770.899−0.0380.0850.42730.5−0.1940.5030.5500.770−0.2510.6530.7150.927−0.0940.2450.26840.8−0.1740.7460.2080.794−0.2190.9400.2620.922−0.0850.3630.1015100.78500.7850100.92400.3830Table 4
The second B-spline curve is mapped to unit quaternions (data partly).
SNVariablePosition of B-spline in Cartesian spaceRotation angleVector of rotation axisUnit quaternionReal partImaginary partuxyzθabcq0q1q2q31000.78500.7850100.92400.383020.20.2690.697−0.1200.7570.3550.921−0.1590.9290.1310.340−0.05930.50.8490.401−0.2070.9620.8830.417−0.2160.8870.4090.193−0.10040.81.3730.101−0.1451.3840.9920.073−0.1050.7700.6330.046−0.067511.571001.5711000.7070.70700Table 5
The third B-spline curve is mapped to unit quaternions (data partly).
SNVariablePosition of B-spline in Cartesian spaceRotation angleVector of rotation axisUnit quaternionReal partImaginary partuxyzθabcq0q1q2q3101.571001.5711000.7070.7070020.21.5530.0450.1811.5640.9930.0290.1150.7100.7000.0200.08130.51.2460.3150.5561.4000.8900.2250.3970.7650.5730.1450.25640.80.7820.6890.9401.4030.5570.4910.6700.7640.3600.3170.432510.4680.9371.1711.5710.2980.5960.7450.7070.2110.4220.527In Table4, when the variable u is equal to 0 and 1, respectively, the interpolated quaternion is [0.924, 0, 0.383, 0] and [0.707, 0.707, 0, 0], which are precisely the values of the quaternions corresponding to the second and third key orientations.Similarly, in Table5, when the variable u is equal to 0 and 1, respectively, the interpolated quaternion is [0.707, 0.7071, 0, 0] and [0.707, 0.211, 0.422, 0.527] which are precisely the values of the quaternions corresponding to the third and fourth key orientations. The results are consistent with the expected values of interpolation theory proposed.When we take evenly 50 times variableu within the range of [0,1], we can get 150 unit quaternions among the three segments of interpolation. In consideration of paper space, only partial interpolation data are listed in the above tables. By now, quaternion orientation interpolation ends.


## 6. Comparison and Discussion

In order to compare with the proposed quaternion interpolation methodology, the most commonly used quaternion interpolation methods SLERP and SQUAD described in the previous section are applied for the same case study above as well. The results of the research and analysis below reflect there are more or less deficiencies for these two methods to be used in multiorientation quaternion interpolation. The comparison and analysis mainly focus on these items: quaternion interpolation curve, rotation angle routing, angular velocity routing, angular acceleration routing, and rotating axis trajectory. The detailed comparison and analysis process are as follows.


### 6.1. Quaternion Interpolation Curves by Proposed Methodology, SQUAD and SLERP

Unit quaternion belongs to S3 space and can be turned by Hopf Mapping to the vector in S2 space [23], which is on the surface of unit sphere in Cartesian space. Then, the interpolated unit quaternions of end-effector orientation can be plotted graphically for observation. By this way, we can intuitively draw the interpolation trajectory diagram for the interpolated four-dimensional unit quaternions on the unit sphere in three-dimensional space. Figure 5 shows the quaternion interpolation curves for the above study case obtained by the methods of proposed methodology, SQUAD and SLERP, respectively.Figure 5
Quaternion interpolation curves. (a) Quaternion interpolation curves by proposed methodology. (b) Quaternion interpolation curves by SQUAD. (c) Quaternion interpolation curves by SLERP.
(a)(b)(c)In Figure5(a), we can see the whole quaternion interpolation curve by proposed methodology is continuous and smooth. There is neither sharp point nor broken point even at the segment connecting points k2 or k3, which means the proposed multiorientation quaternion interpolation is reasonable. The same is shown in Figure 5(b) by the SQUAD method. However, we cannot tell whether the proposed methodology or SQUDA is C2 continuous only from the quaternion interpolation curves. Because both of their quaternion interpolation curves seem smooth enough in appearance, further analysis and discussion are needed in below section.As for the method SLERP, the situation becomes a little simpler. We can clearly see that sharp points occur at the segment connecting pointsk2 or k3 in Figure 5(c) by SLERP. This result fully shows that the quaternion interpolation of SLERP for multiple orientations is not smooth despite being continuous and will bring violent impacts and shocks to end-effector. So, we can immediately come to a conclusion that SLERP is not suitable for multiorientation quaternion interpolation.


### 6.2. Rotation Angle Routings by Proposed Methodology, SQUAD and SLERP

We can get the values of rotation angle, angular velocity, and angular acceleration according to these three interpolation formulas above. Figure6 shows rotation angle routing curves of the above study case obtained by the methods of proposed methodology, SQUAD and SLERP, respectively. The angle curve in Figure 6(c) has sharp points θc1 and θc2, which implies the curve is continuous but not smooth enough and rotation angle arrangement is not reasonable. The defect led by unsmooth rotation angle planning of multiorientation interpolation of SLERP is more obvious in later analysis.Figure 6
Angular routings. (a) Angular routing by proposed methodology. (b) Angular routing by SQUAD. (c) Angular routing by SLERP.
(a)(b)(c)The curves in Figures6(a) and 6(b) are very similar in appearance. Both of them look smooth and continuous, which means we cannot judge from rotation angle routing curves which one is the better choice in proposed methodology and SQUAD. Analysis below is needed.


### 6.3. Angular Velocity Routings by Proposed Methodology, SQUAD and SLERP

Angular velocity routing curves of these three methods for the case are shown in Figure7. There are very obvious broken angular velocity points ωc1 and ωc2 in Figure 7(c), which shows clearly that jump changes take place in angular velocity and will bring violent impacts and shocks harmful to industrial robots. Figure 7(c) demonstrates exactly once more that SLERP is not suitable for multiorientation quaternion interpolation.Figure 7
Angular velocity routings. (a) Angular velocity routing by proposed methodology. (b) Angular velocity routing by SQUAD. (c) Angular velocity routing by SLERP.
(a)(b)(c)The angular velocity routing curves in Figure7(a) and 7(b) are still similar. Because both of them are continuous and have no obvious sharp point, it is still difficult to determine which one is better in the proposed methodology and SQUAD by angular velocity.


### 6.4. Angular Acceleration Routings by Proposed Methodology, SQUAD and SLERP

Angular acceleration routing curves of these three methods for the case are shown in Figure8. The angular acceleration routing curve in Figure 8(a) obtained by proposed methodology is continuous while the curve in Figure 8(b) by SQUAD is broken by the points εc1 and εc2. In other words, the quaternion interpolation SQUAD for multiple orientations is not C2 continuous and the angular acceleration is discontinuous, which will also lead to violent shocks and impacts more or less and degrade the working performance and service life of industrial robot as well. Therefore, SQUAD is not perfect quaternion interpolation for multiple orientation. The angular acceleration routing curves in Figure 8(c) obtained by SLERP are discontinuous expectedly. The previous analysis has made it very clear that SLERP is not suitable for multiorientation quaternion interpolation.Figure 8
Angular acceleration routings. (a) Angular acceleration routing by proposed methodology. (b) Angular acceleration routing by SQUAD. (c) Angular acceleration routing by SLERP.
(a)(b)(c)Up to now, we can tell evidentially that the proposed methodology is superior to SLERP and SQUAD in multiorientation quaternion interpolation.


### 6.5. Rotation Axis Trajectories by Proposed Methodology, SQUAD and SLERP

We can also observe which method is more advantageous from rotation axis trajectories as shown in Figure9. The rotation axis trajectory curve is from the starting orientation to the ending orientation of end-effector during the whole orientation planning.Figure 9
Rotation axis trajectories. (a) Rotation axis trajectory by proposed methodology. (b) Rotation axis trajectory by SQUAD. (c) Rotation axis trajectory by SLERP.
(a)(b)(c)The rotation axis trajectory by the proposed method in Figure9(a) is smooth enough. Although the rotation axis trajectory curve by SQUAD in Figure 9(b) is continuous and smooth, we can still observe that the rotation axis trajectory curve by proposed methodology in Figure 9(a) is more fluid. Moreover, the rotation axis trajectory by SLERP in Figure 9(c) has apparent sharp points n2 and n3 at the connecting corners and is obviously not smooth. It can be concluded as well that SLERP is not suitable for multiorientation quaternion interpolation, and proposed methodology is more advantageous than SQUAD and SLERP in multiorientation quaternion interpolation.


## 6.1. Quaternion Interpolation Curves by Proposed Methodology, SQUAD and SLERP

Unit quaternion belongs to S3 space and can be turned by Hopf Mapping to the vector in S2 space [23], which is on the surface of unit sphere in Cartesian space. Then, the interpolated unit quaternions of end-effector orientation can be plotted graphically for observation. By this way, we can intuitively draw the interpolation trajectory diagram for the interpolated four-dimensional unit quaternions on the unit sphere in three-dimensional space. Figure 5 shows the quaternion interpolation curves for the above study case obtained by the methods of proposed methodology, SQUAD and SLERP, respectively.Figure 5
Quaternion interpolation curves. (a) Quaternion interpolation curves by proposed methodology. (b) Quaternion interpolation curves by SQUAD. (c) Quaternion interpolation curves by SLERP.
(a)(b)(c)In Figure5(a), we can see the whole quaternion interpolation curve by proposed methodology is continuous and smooth. There is neither sharp point nor broken point even at the segment connecting points k2 or k3, which means the proposed multiorientation quaternion interpolation is reasonable. The same is shown in Figure 5(b) by the SQUAD method. However, we cannot tell whether the proposed methodology or SQUDA is C2 continuous only from the quaternion interpolation curves. Because both of their quaternion interpolation curves seem smooth enough in appearance, further analysis and discussion are needed in below section.As for the method SLERP, the situation becomes a little simpler. We can clearly see that sharp points occur at the segment connecting pointsk2 or k3 in Figure 5(c) by SLERP. This result fully shows that the quaternion interpolation of SLERP for multiple orientations is not smooth despite being continuous and will bring violent impacts and shocks to end-effector. So, we can immediately come to a conclusion that SLERP is not suitable for multiorientation quaternion interpolation.


## 6.2. Rotation Angle Routings by Proposed Methodology, SQUAD and SLERP

We can get the values of rotation angle, angular velocity, and angular acceleration according to these three interpolation formulas above. Figure6 shows rotation angle routing curves of the above study case obtained by the methods of proposed methodology, SQUAD and SLERP, respectively. The angle curve in Figure 6(c) has sharp points θc1 and θc2, which implies the curve is continuous but not smooth enough and rotation angle arrangement is not reasonable. The defect led by unsmooth rotation angle planning of multiorientation interpolation of SLERP is more obvious in later analysis.Figure 6
Angular routings. (a) Angular routing by proposed methodology. (b) Angular routing by SQUAD. (c) Angular routing by SLERP.
(a)(b)(c)The curves in Figures6(a) and 6(b) are very similar in appearance. Both of them look smooth and continuous, which means we cannot judge from rotation angle routing curves which one is the better choice in proposed methodology and SQUAD. Analysis below is needed.


## 6.3. Angular Velocity Routings by Proposed Methodology, SQUAD and SLERP

Angular velocity routing curves of these three methods for the case are shown in Figure7. There are very obvious broken angular velocity points ωc1 and ωc2 in Figure 7(c), which shows clearly that jump changes take place in angular velocity and will bring violent impacts and shocks harmful to industrial robots. Figure 7(c) demonstrates exactly once more that SLERP is not suitable for multiorientation quaternion interpolation.Figure 7
Angular velocity routings. (a) Angular velocity routing by proposed methodology. (b) Angular velocity routing by SQUAD. (c) Angular velocity routing by SLERP.
(a)(b)(c)The angular velocity routing curves in Figure7(a) and 7(b) are still similar. Because both of them are continuous and have no obvious sharp point, it is still difficult to determine which one is better in the proposed methodology and SQUAD by angular velocity.


## 6.4. Angular Acceleration Routings by Proposed Methodology, SQUAD and SLERP

Angular acceleration routing curves of these three methods for the case are shown in Figure8. The angular acceleration routing curve in Figure 8(a) obtained by proposed methodology is continuous while the curve in Figure 8(b) by SQUAD is broken by the points εc1 and εc2. In other words, the quaternion interpolation SQUAD for multiple orientations is not C2 continuous and the angular acceleration is discontinuous, which will also lead to violent shocks and impacts more or less and degrade the working performance and service life of industrial robot as well. Therefore, SQUAD is not perfect quaternion interpolation for multiple orientation. The angular acceleration routing curves in Figure 8(c) obtained by SLERP are discontinuous expectedly. The previous analysis has made it very clear that SLERP is not suitable for multiorientation quaternion interpolation.Figure 8
Angular acceleration routings. (a) Angular acceleration routing by proposed methodology. (b) Angular acceleration routing by SQUAD. (c) Angular acceleration routing by SLERP.
(a)(b)(c)Up to now, we can tell evidentially that the proposed methodology is superior to SLERP and SQUAD in multiorientation quaternion interpolation.


## 6.5. Rotation Axis Trajectories by Proposed Methodology, SQUAD and SLERP

We can also observe which method is more advantageous from rotation axis trajectories as shown in Figure9. The rotation axis trajectory curve is from the starting orientation to the ending orientation of end-effector during the whole orientation planning.Figure 9
Rotation axis trajectories. (a) Rotation axis trajectory by proposed methodology. (b) Rotation axis trajectory by SQUAD. (c) Rotation axis trajectory by SLERP.
(a)(b)(c)The rotation axis trajectory by the proposed method in Figure9(a) is smooth enough. Although the rotation axis trajectory curve by SQUAD in Figure 9(b) is continuous and smooth, we can still observe that the rotation axis trajectory curve by proposed methodology in Figure 9(a) is more fluid. Moreover, the rotation axis trajectory by SLERP in Figure 9(c) has apparent sharp points n2 and n3 at the connecting corners and is obviously not smooth. It can be concluded as well that SLERP is not suitable for multiorientation quaternion interpolation, and proposed methodology is more advantageous than SQUAD and SLERP in multiorientation quaternion interpolation.


## 7. Conclusions

High-order continuous orientation planning is beneficial to working performance and service life of robot end-effector. However, the most commonly used quaternion interpolation algorithms, SLERP and SQUAD, cannot guarantee C2 continuity of multiorientation interpolation planning. When SLERP is applied for quaternion multiorientation interpolation, it will generate intermittent angular velocity and angular acceleration at the connecting point, which will bring additional impacts and vibrations to the robot and is very harmful to the working performance of robots. When SQUAD is employed for quaternion multiorientation interpolation, it will also produce discontinuous angular acceleration at the connection point, which will also lead to additional shocks and vibrations more or less and are not conducive to the robot’s work either. Therefore, SLERP and SQUAD are not suitable for quaternion multiorientation planning of robot end-effector.The proposed quaternion interpolation method with the cubic B-spline curve based on logarithmic quaternion is able to guarantee C2 continuity. When the proposed method is applied to the multiorientation planning of robot end-effector, the planned angular velocity and angular acceleration are both continuous and smooth avoiding additional impacts and vibrations unfavourable for robot working. By the way, in the procedure of reversely finding control points of the cubic B-spline curve, the constraint function condition is that the curvatures at the beginning and the terminal knots of the B-spline curve are set to zero. It helps the quaternion interpolation become smoother and sleeker. Moreover, the formula expression of the proposed methodology is concise and its application is simple and convenient. The study case verified the C2 continuity of quaternion interpolation and its rationality. Therefore, it is suitable for quaternion multiorientation planning of robot end-effector. Furthermore, the proposed methodology can be extended to multiorientation quaternion interpolation with higher continuity than the second order.



---

*Source: 2543824-2020-07-22.xml*