By Ilian Bonev, Ph.D., Eng.

The concept of quaternions is fascinating with various applications, but in robotics we care about quaternions mainly as a very efficient parameterization of 3D orientation, primarily that of the robot end-effector. There are dozens of in-depth presentations on quaternions, but none is meant for industrial robot programmers. To fill the void, this tutorial presents the strict minimum about quaternions that users of industrial robots should be familiar with.  

Quaternions as an extension of complex numbers

Quaternions were invented by the Irish mathematician W. R. Hamilton in 1843. They are essentially an extension of complex numbers. As you have learned in high school, a complex number is represented as

a + bi,

where a and b are real numbers, and i is a so-called unit imaginary number such that i² = −1.

“Quaternions are essentially an extension of complex numbers.”

In similar fashion, a quaternion is represented as

q = q₁ + q₂ i + q₃ j + q₄ k,

where i² = j² = k² = −1, ij = k = −ji, jk = i = −kj, and ki = j = −ki. Often, we simply denote quaternions as the quadruple {q₁, q₂, q₃, q₄}. Using the identities just given, if p = {p₁, p₂, p₃, p₄} is another quaternion, it can be easily shown that the product of q times p, called the Hamilton product, is defined by

qp = (q₁p₁ – q₂p₂ – q₃p₃ – q₄p₄) + (q₁p₂ + q₂p₁ + q₃p₄ – q₄p₃)i + (q₁p₃ + q₃p₁ – q₂p₄ + q₄p₂)j + (q₁p₄ + q₄p₁ + q₂p₃ – q₃p₂)k.

It can also be easily demonstrated that quaternion multiplication is not commutative, meaning that in general qp ≠ pq. In contrast, the addition of quaternions is commutative and is performed component-wise:

q + p = p + q = (q₁ + p₁) + (q₂ + p₂)i + (q₃ + p₃)j + (q₄ + p₄)k.

Euler’s rotation theorem

In 1775, the Swiss scientist Leonhard Euler proved that any reorientation of a rigid body is equivalent to a single rotation about an axis. Thus, any orientation in 3D space can be represented by a unit vector u along this axis and the angle of the rotation, θ, where the direction of rotation is determined by the right-hand rule. We will therefore use the term directed axis defined by a unit vector, from now on, and won’t repeat that the direction of positive rotations about this directed axis is measured using the right-hand rule. Thus, if u and θ represent an orientation in space, then −u and −θ represent the same orientation, and so do −u and 360°−θ.

“Any reorientation of a rigid body is equivalent to a single rotation about an axis.”

Imagine that you jog your robot’s end-effector in reorientation mode about its TCP (tool center point). Even after a long sequence of rotations about the x, y, and z axes of the tool reference frame, you can still go back to your original orientation using a single rotation of less than 180° about some axis passing through the TCP. In fact, if you execute a linear motion command (MoveLin, in the case of Mecademic robots) with the same unchanged position of the TCP but the original orientation of the tool reference frame as arguments, the end-effector will perform exactly that single rotation of less than 180° about an axis passing through the TCP. We will get back to this so-called SLERP interpolation later.

Quaternions as rotations in 3D space

A unit quaternion, q, such that |q| = q₁² + q₂² + q₃² + q₄² = 1, can also be represented as

q = {cos(θ/2), u_x sin(θ/2), u_y sin(θ/2), u_z sin(θ/2)},

where u_x, u_y, and u_z are the components of a unit vector u, and θ is some angle. It can be shown, although that’s not trivial, that if we assume the above unit quaternion to represent a rotation of θ about a directed axis defined by vector u, we will end up with some elegant properties. For example, since the same orientation can be represented as a rotation of 360°−θ  about a directed axis defined by vector −u, it can be easily shown that

q = {cos(180°−θ/2), −u_x sin(180°−θ/2), −u_y sin(180°−θ/2), −u_z sin(180°−θ/2)} = {−cos(θ/2), −u_x sin(θ/2), −u_y sin(θ/2), −u_z sin(θ/2)} = −q.

In other words, if we negate all the terms of a quaternion, the new and the old one represent the same orientation.

Unit quaternions as transformations

Now, let R_q be the 3×3 rotation matrix, that you are well familiar with, representing the same orientation as the quaternion q, and let R_p be the rotation matrix representing the same orientation as another unit quaternion p. Then, the orientation represented by the matrix corresponding to the product R_qR_p is the same as the orientation representing the quaternion product qp (which was defined at the very beginning). Similarly, the orientation represented by the matrix corresponding to the product R_pR_q is the same as the orientation represented by the quaternion product pq. Thus, while computing the transformation corresponding to two successive rotations requires 27 multiplications and 18 additions/subtractions when using rotation matrices, the same operation requires only 16 multiplications and 12 additions/subtractions when using quaternions. As you probably remember, the inverse of a rotation matrix is the same as its transpose:

R_q⁻¹ = R_q^T

The inverse of a quaternion is not only easier to calculate, it is also more intuitive. The inverse of a rotation of θ about a directed axis defined by unit vector u is simply a rotation of −θ about the same directed axis. In other words, the inverse of the unit quaternion q is:

q⁻¹ = {cos(−θ/2), u_x sin(−θ/2), u_y sin(−θ/2), u_z sin(−θ/2)} = {cos(θ/2), −u_x sin(θ/2), −u_y sin(θ/2), −u_z sin(θ/2)}.

Thus, we have

q⁻¹q = qq⁻¹ = {1, 0, 0, 0},

which corresponds to no orientation (similar to the identity matrix). Finally, you are surely familiar with the way coordinates expressed in one reference frame are expressed in another, using rotation matrices. This operation is equivalent to the following triple quaternion product:

qvq⁻¹,

where v = {0, v_x, v_y, v_z} is a so-called pure quaternion, and v_x, v_y, and v_z are the coordinates to be transformed. Multiplying a rotation matrix times a position vector requires 9 multiplications and 6 additions/subtractions. The equivalent operation using quaternions requires more than double that: 18 multiplications and 13 additions/subtractions. As we just saw, the unit quaternion can represent a rotation, while a pure quaternion can represent a translation. But how do we combine these to perform transformations in space, the way we do with 4×4 homogeneous matrices? This operation is performed using so-called dual-quaternions. In this tutorial, we won’t present dual-quaternions, because these are hardly of any use to the typical robot programmer. However, it is worth mentioning their advantages over homogenous matrices. In terms of computational efficiency, both methods are similar, but dual quaternions are more compact and more numerically robust.

Converting between rotation matrices and unit quaternions

While quaternions are very useful — and sometimes even indispensable — when dealing with 3D orientations, to define the quaternion corresponding to the orientation of one reference frame with respect to another, you generally need to first calculate the corresponding rotation matrix. We therefore need the following formula for computing the components of a quaternion, q, from the elements of a rotation matrix, R:

q₁ = 1/2 √1 + r_1,1 + r_2,2 + r_3,3,

q₂ = 1/2sgn(r_3,2 − r_2,3)√1 + r_1,1 − r_2,2 − r_3,3,

q₃ = 1/2sgn(r_1,3 − r_3,1)√1 − r_1,1 + r_2,2 − r_3,3,

q₄ = 1/2sgn(r_2,1 − r_1,2)√1 − r_1,1 − r_2,2 + r_3,3,

where

R =	r1,1 r1,2 r1,3 r2,1 r2,2 r2,3 r3,1 r3,2 r3,3 .

Of course, recall that if q is the unit quaternion corresponding to the rotation matrix R, then −q is the other quaternion corresponding to the same orientation. Similarly, if you have a unit quaternion q and you wish to convert it to a rotation matrix R, the formula is:

R =	2q12 + 2q22 − 1 2q2q3−2q1q4 2q2q4 + 2q1q3 2q2q3 + 2q1q4 2q12 + 2q32 − 1 2q3q4 − 2q1q2 2q2q4 − 2q1q3 2q3q4 − 2q1q2 2q12 + 2q42 − 1 .

Using unit quaternions to measure rotational displacements

As we saw in our tutorial on Euler angles, these angles are almost always the only way orientation is defined by users and represented to them in the case of industrial robots. Thus, it is a very common mistake to draw a parallel between the nature of position coordinates and Euler angles when it comes to measuring a relative displacement, or a “distance.” In fact, even the norm ISO 9283, which dictates how robot performance should be measured and represented, states that orientational accuracy is essentially measured with exactly the same formula as positional accuracy, but with Euler angles replacing the x, y and z position coordinates. Yet, similar to representational singularities of the Euler angles (commonly known as the Gimbal lock effect), Euler angles can vary a lot, even if the orientation changes only slightly.

The one and only measure of rotational displacement, i.e., the “distance” between an initial orientation and the final one, can be easily calculated using unit quaternions. Indeed, let q be the unit quaternion representing the orientation of the tool reference frame at pose 1, and p be the unit quaternion representing the orientation of the tool reference frame at pose 2. The orientation of the tool frame at pose 2 with respect to the tool frame at pose 1 is

d = q⁻¹p = {cos(φ/2), w_x sin(φ/2), w_y sin(φ/2), w_z sin(φ/2)},

where w_x, w_y, and w_z are the components of a unit vector w. In other words, to get from the orientation of the tool frame at pose 1 to the orientation at pose 2, the tool frame has to rotate φ about a directed axis defined by vector w. Now, since we are interested only in the shortest possible rotation, we want 0 < φ ≤ 180°. Therefore, if we see that the first component of the resulting quaternion d is negative, we should simply negate the quaternion. (Recall that d = −d.) Then, to calculate the equivalent shortest rotation between the two orientations, we should calculate the inverse cosine of the first component of d and multiply it by 2. Similarly, many wrongly assume that the angular velocity of a body in space is the vector of the derivatives of the Euler angles. The angular velocity vector, ω, is in fact the unit vector corresponding to the directed instantaneous axis of rotation multiplied by the rate of rotation about this axis.

Using quaternions to interpolate between orientations

As we have just seen, the fastest way to get from one orientation to another is to rotate about the equivalent axis of rotation at an angle of no more than 180°. This rotational path is often referred to as the torque-minimal path. Thus, to get an “equidistant” series of orientations that get you from an orientation represented by the quaternion q, to an orientation represented by the quaternion p, we must use the following simple formula:

Slerp(q,p,t) = q(q⁻¹p)^t,

where t is a parameter varying from 0 (when at orientation q) to 1 (when at orientation p). The above procedure is called spherical linear interpolation, or SLERP. Raising the quaternion d = q⁻¹p = {cos(φ/2), w_x sin(φ/2), w_y sin(−φ/2), w_z sin(−φ/2)} to the power t is equivalent to rotating t times at an angle φ, so

d^t = {cos(tφ/2), w_x sin(tφ/2), w_y sin(−tφ/2), w_z sin(−tφ/2)}.

Thus, if t = 0, 1/n, 2/n, …, 1, SLERP is essentially “adding” a rotation of φ/n at every discretization step, where n − 1 is the number of in-between orientations. Of course, φ must be no more than 180°, otherwise you must negate one of the two quaternions and repeat the calculations. It can be shown that the above SLERP formula is equivalent to the following equation, which is a bit faster to compute:

Slerp(q,p,t) = sin((1 − t)φ/2)/sin(φ/2) q + sin(tφ/2)/sin(φ/2) p,

where φ is as just defined, i.e., the angle of rotation between q and p. Note that both SLERP functions are equivalent and the result is a unit quaternion for any t. A faster way to interpolate between orientations along the same rotation path about the equivalent axis of rotation is using the following simple formula:

Lerp(q,p,t) = (1 − t)q + tp,

where t is again a parameter varying between 0 and 1. However, this so-called LERP interpolation works relatively well only for orientations that are close to each other (e.g., φ < 90°). This is because the reorientation is not performed at constant angular speed. Besides, with LERP, you need to normalize the resulting quaternion at each step t (then the method is actually called NLERP). Most industrial robots use a SLERP interpolation when reorienting the end-effector along a linear path. Obviously, if the start and end orientations are 180° apart, then there would be two identical minimum-torque paths. The robot controller can arbitrarily choose one of them, but that means that the user cannot anticipate the motion of the end-effector. Thus, industrial robots always have some kind of limits on the displacements along a linear path. For example, Mecademic robots simply refuse a linear motion towards a pose with an orientation that is exactly 180° away.

Quaternions vs Euler angles

Whenever I posted about Euler angles, there would always be people criticizing them and praising quaternions. Indeed, Euler angles — and any three-parameter representation of orientation in space — have representational singularities whereby an orientation can be represented with infinitely many triplets. (This is similar to the problem of representing the coordinates at the North or South Poles on Earth — the latitude is not uniquely defined at both poles.)

Euler angles are surely inadequate for dealing with orientations that vary. However, in robotics, Euler angles are typically used to define one fixed reference frame with respect to another, such as the tool reference frame with respect to the flange reference frame. As shown in the figure above, it is generally impossible to calculate by hand the quaternion corresponding to the orientation of one reference frame with respect to another, whereas it is often easy to guess the Euler angles by trial and error.

“Euler angles are surely inadequate for dealing with orientations that vary.”

There do exist situations, though, where it is impossible to guess Euler angles also, and you need to calculate these using trigonometric formulas involving the elements of a rotation matrix. For example, in the figure above, for the orientation of frame 1 with respect to frame 0, the Euler angles according to the mobile XYZ convention used by Mecademic Robotics are {0°, 45°, 40°}, but the Euler angles according to the mobile ZYX convention used by KUKA are {49.879°, 32.798°, 32.732°}, and the quaternions as used by ABB are {0.868163, 0.130885, 0.359605, 0.315986}. Nevertheless, it is not surprising that virtually all industrial robots require users to define orientations using Euler angles, as these are easier to understand, though not trivial. ABB is a rare exception, where quaternions are the main way of entering orientations, though Euler angles are an option too. That said, all industrial robots do use quaternions inside their controllers, be it for orientation interpolation, calculation of transformations, or for other purposes.

Summary

To recap, as an industrial robot programmer you should know about quaternions for at least the following reasons:

• be able to anticipate the reorientation of the robot end-effector along a linear path (the SLERP interpolation)
• know how to correctly measure orientational accuracy and repeatability
• better understand the drawbacks of Euler angles
• know how to use quaternions when these are available

© Mecademic. Reproducing this tutorial, in whole or in part, is strictly prohibited.