Quaternion

Quaternions were initially defined by William Hamilton in 1843 to describe a Cayley-Dickson construction in four dimensions.

Since then, many interpretations have appeared for different applications. The most common definition of a quaternion \(\mathbf{q}\) is as an ordered expression of the form:

\[\mathbf{q} = w + xi + yj + zk\]

where \(w\), \(x\), \(y\) and \(z\) are real numbers, and \(i\), \(j\) and \(k\) are three imaginary unit numbers defined so that [Sola] [Kuipers] [WikiQuaternion] :

\[i^2 = j^2 = k^2 = ijk = -1\]

It is helpful to notice that quaternions are arrays with the same structure: a real number \(w\) followed by three pairs of complex numbers. Thus, we can split it to find a second definition:

\[\mathbf{q} = w + v\]

where \(v=(xi + yj + zk)\in\mathbb{Z}\) and \(w\in\mathbb{R}\). Recognizing that \(i\), \(j\) and \(k\) are always the same, we can omit them to redefine the quaternion as a pair of a scalar and a vector:

\[\mathbf{q} = q_w + \mathbf{q}_v\]

where \(q_w\in\mathbb{R}\) is called the real or scalar part, while \(\mathbf{q}_v = (q_x, q_y, q_z)\in\mathbb{R}^3\) is the imaginary or vector part.

For practical reasons, most literature will use the vector definition of a quaternion [1]:

\[\begin{split}\mathbf{q}\triangleq \begin{bmatrix}q_w \\ \mathbf{q}_v\end{bmatrix} = \begin{bmatrix}q_w \\ q_x \\ q_y \\ q_z\end{bmatrix}\end{split}\]

Sadly, many authors use different notations for the same type of quaternions. Some even invert their order, with the vector part first followed by the scalar part, increasing the confusion among readers. Here, the definition above will be used throughout the package.

Let’s say, for example, we want to use the quaternion \(\mathbf{q}=\begin{pmatrix}0.7071 & 0 & 0.7071 & 0\end{pmatrix}\) with this class:

>>> from ahrs import Quaternion
>>> q = Quaternion([0.7071, 0.0, 0.7071, 0.0])
>>> q
Quaternion([0.70710678, 0.        , 0.70710678, 0.        ])

This will extend the values of the quaternion, because it is handled as a rotation operator.

Tip

You can have a pretty formatting of the quaternion if typecasted as a string, showing it with Hamilton’s notation:

>>> str(q)
'(0.7071 +0.0000i +0.7071j +0.0000k)'

As Rotation Operators

Quaternions can be defined in the geometric space as an alternative form of a rotation operator, so that we can find the image \(\mathbf{a}'\) of some vector \(\mathbf{a}\) using a product similar to other euclidean transformations in 3D:

\[\mathbf{a}' = \mathbf{qa}\]

Back in the XVIII century Leonhard Euler showed that any rotation around the origin can be described by a three-dimensional axis and a rotation magnitude. Euler’s Rotation Theorem is the basis for most three-dimensional rotations out there.

Later, in 1840, Olinde Rodrigues came up with a formula using Euler’s principle in vector form:

\[\mathbf{a}' = \mathbf{a}\cos\theta + (\mathbf{v}\times\mathbf{a})\sin\theta + \mathbf{v}(\mathbf{v}\cdot\mathbf{a})(1-\cos\theta)\]

where \(\mathbf{v}=\begin{bmatrix}v_x & v_y & v_z\end{bmatrix}\) is a unit vector [2] describing the axis of rotation about which \(\mathbf{a}\) rotates by an angle \(\theta\).

The Euler-Rodrigues Formula can be further compacted to:

\[\mathbf{a}' = \mathbf{a} + 2\alpha (\vec{\mathbf{v}}\times\mathbf{a}) + 2\big(\vec{\mathbf{v}}\times(\vec{\mathbf{v}}\times\mathbf{a})\big)\]

where \(\alpha = \cos\frac{\theta}{2}\) and \(\vec{\mathbf{v}}\) represents the vector of rotation as half angles:

\[\begin{split}\vec{\mathbf{v}} = \mathbf{v}\sin\frac{\theta}{2} = \begin{bmatrix} v_x \sin\frac{\theta}{2} \\ v_y \sin\frac{\theta}{2} \\ v_z \sin\frac{\theta}{2} \end{bmatrix}\end{split}\]

This rotation is defined by a scalar and vector pair. Now we can see that the structure of a quaternion can be used to describe a rotation with an axis-angle representation such that:

\[\begin{split}\mathbf{q} = \begin{bmatrix}q_w \\ q_x \\ q_y \\ q_z\end{bmatrix} = \begin{bmatrix} \cos\frac{\theta}{2} \\ v_x \sin\frac{\theta}{2} \\ v_y \sin\frac{\theta}{2} \\ v_z \sin\frac{\theta}{2} \end{bmatrix} = \begin{bmatrix} \cos\frac{\theta}{2} \\ \mathbf{v} \sin\frac{\theta}{2} \end{bmatrix}\end{split}\]

This is only valid if the quaternion is normalized, in which case it is also known as a versor. To enforce the versor, we normalize the quaternion:

\[\begin{split}\mathbf{q} = \frac{1}{\sqrt{q_w^2+q_x^2+q_y^2+q_z^2}}\begin{bmatrix}q_w \\ q_x \\ q_y \\ q_z\end{bmatrix}\end{split}\]

In this module the quaternions are considered rotation operators, so they will always be normalized. From the example above:

>>> q = Quaternion([0.7071, 0.0, 0.7071, 0.0])
>>> q
Quaternion([0.70710678, 0.        , 0.70710678, 0.        ])
>>> import numpy as np
>>> np.linalg.norm(q)
1.0

Very convenient conversion methods are reachable in this class. One of them is the representation of direction cosine matrices from the created quaternion.

>>> q.to_DCM()
array([[-2.22044605e-16,  0.00000000e+00,  1.00000000e+00],
       [ 0.00000000e+00,  1.00000000e+00,  0.00000000e+00],
       [-1.00000000e+00,  0.00000000e+00, -2.22044605e-16]])

Some important observations here help us to clarify this further:

  • When \(\theta=0\) (no rotation) the parametrization becomes \(\mathbf{q} = \begin{pmatrix}1 & 0 & 0 & 0\end{pmatrix}\). This is identified as an identity quaternion.
  • When \(\theta=180\) (half-circle), then \(q_w=0\), making \(\mathbf{q}\) a pure quaternion.
  • The negative values of a versor, \(\begin{pmatrix}-q_w & -q_x & -q_y & -q_z\end{pmatrix}\) represent the same rotation [Grosskatthoefer].
>>> q = Quaternion([1.0, 0.0, 0.0, 0.0])
>>> q.is_identity()
True
>>> q = Quaternion([0.0, 1.0, 2.0, 3.0])
>>> q
Quaternion([0.        , 0.26726124, 0.53452248, 0.80178373])
>>> q.is_pure()
True
>>> q1 = Quaternion([1.0, 2.0, 3.0, 4.0])
>>> q2 = Quaternion([-1.0, -2.0, -3.0, -4.0])
>>> np.all(q1.to_DCM()==q2.to_DCM())
True

To summarize, unit quaternions can also be defined using Euler’s rotation theorem [3] in the form [Kuipers]:

\[\begin{split}\mathbf{q} = \begin{bmatrix}\cos\theta \\ \mathbf{v}\sin\theta\end{bmatrix}\end{split}\]

And they have similar algebraic characteristics as normal vectors [Sola] [Eberly]. For example, given two quaternions \(\mathbf{p}\) and \(\mathbf{q}\):

\[\begin{split}\mathbf{p} \pm \mathbf{q} = \begin{bmatrix}p_w\pm q_w \\ \mathbf{p}_v \pm \mathbf{q}_v \end{bmatrix}\end{split}\]
>>> p = Quaternion([1., 2., 3., 4.])
>>> p
Quaternion([0.18257419, 0.36514837, 0.54772256, 0.73029674])
>>> q = Quaternion([-5., 4., -3., 2.])
>>> q
Quaternion([-0.68041382,  0.54433105, -0.40824829,  0.27216553])
>>> p+q
Quaternion([-0.34359264,  0.62769298,  0.09626058,  0.6918667 ])
>>> p-q
Quaternion([ 0.62597531, -0.1299716 ,  0.69342116,  0.33230917])

The quaternion product uses the Hamilton product to perform their multiplication [4], which can be represented in vector form, as a known scalar-vector form, or even as a matrix multiplication:

\[\begin{split}\begin{array}{rl} \mathbf{pq} &= \begin{bmatrix} p_wq_w-\mathbf{p}_v^T\mathbf{q}_v \\ p_w\mathbf{q}_v + q_w\mathbf{p}_v + \mathbf{p}_v\times\mathbf{q}_v \end{bmatrix} \\ &= \begin{bmatrix} p_w q_w - p_x q_x - p_y q_y - p_z q_z \\ p_w q_x + p_x q_w + p_y q_z - p_z q_y \\ p_w q_y - p_x q_z + p_y q_w + p_z q_x \\ p_w q_z + p_x q_y - p_y q_x + p_z q_w \end{bmatrix} \\ &= \begin{bmatrix} p_w & -p_x & -p_y & -p_z \\ p_x & p_w & -p_z & p_y \\ p_y & p_z & p_w & -p_x \\ p_z & -p_y & p_x & p_w \end{bmatrix} \begin{bmatrix} q_w \\ q_x \\ q_y \\ q_z \end{bmatrix} \\ &= \begin{bmatrix} p_w & -\mathbf{p}_v^T \\ \mathbf{p}_v & p_w \mathbf{I}_3 + \lfloor \mathbf{p}_v \rfloor_\times \end{bmatrix} \begin{bmatrix} q_w \\ \mathbf{q}_v \end{bmatrix} \end{array}\end{split}\]

which is not commutative

\[\mathbf{pq} \neq \mathbf{qp}\]
>>> p*q
array([-0.2981424 ,  0.2981424 , -0.1490712 , -0.89442719])
>>> q*p
array([-2.98142397e-01, -5.96284794e-01, -7.45355992e-01,  4.16333634e-17])

The conjugate of the quaternion, defined as \(\mathbf{q}^*=\begin{pmatrix}q_w & -\mathbf{q}_v\end{pmatrix}\), has the interesting property of:

\[\begin{split}\mathbf{qq}^* = \mathbf{q}^*\mathbf{q} = \begin{bmatrix}q_w^2+q_x^2+q_y^2+q_z^2 \\ \mathbf{0}_v\end{bmatrix}\end{split}\]

But for the case, where the quaternion is a versor, like in this module, the conjugate is actually the same as the inverse \(\mathbf{q}^{-1}\), where:

\[\mathbf{qq}^* = \mathbf{qq}^{-1} = \mathbf{q}^{-1}\mathbf{q} = \begin{pmatrix}1 & 0 & 0 & 0\end{pmatrix}\]
>>> q*q.conj
array([1., 0., 0., 0.])
>>> q*q.inv
array([1., 0., 0., 0.])

Rotation groups are normally represented with direction cosine matrices. It is undeniable that they are the best way to express any rotation operation. However, quaternions are also a good representation of it, and even a better ally for numerical operations.

Footnotes

[1]Some authors use different subscripts, but here we mainly deal with geometric transformations and the notation (w, x, y, z) is preferred.
[2]Any unit vector \(\mathbf{v}\in\mathbb{R}^3\) has, per definition, a magnitude equal to 1, which means \(\|\mathbf{v}\|=1\).
[3]This is the reason why some early authors call the quaternions Euler Parameters.
[4]Many authors decide to use the symbol \(\otimes\) to indicate a quaternion product, but here is not used in order to avoid any confusions with the outer product.

References

[Bar-Itzhack]Y. Bar-Itzhack. New method for Extracting the Quaternion from a Rotation Matrix. Journal of Guidance, Control, and Dynamics, 23(6):1085–1087, 2000. (https://arc.aiaa.org/doi/abs/10.2514/2.4654)
[Chiaverini]S. Chiaverini & B. Siciliano. The Unit Quaternion: A Useful Tool for Inverse Kinematics of Robot Manipulators. Systems Analysis Modelling Simulation. May 1999. (https://www.researchgate.net/publication/262391661)
[Dantam]Dantam, N. (2014) Quaternion Computation. Institute for Robotics and Intelligent Machines. Georgia Tech. (http://www.neil.dantam.name/note/dantam-quaternion.pdf)
[Eberly]Eberly, D. (2010) Quaternion Algebra and Calculus. Geometric Tools. https://www.geometrictools.com/Documentation/Quaternions.pdf
[Grosskatthoefer]K. Grosskatthoefer. Introduction into quaternions from spacecraft attitude representation. TU Berlin. 2012. (http://www.tu-berlin.de/fileadmin/fg169/miscellaneous/Quaternions.pdf)
[Hughes]
  1. Hughes. Spacecraft Attitude Dynamics. 1986. p. 18
[Kuipers](1, 2) Kuipers, Jack. Quaternions and Rotation Sequences. Princenton University Press. 1999.
[Markley2007]F. Landis Markley. Averaging Quaternions. Journal of Guidance, Control, and Dynamics. Vol 30, Num 4. 2007 (https://arc.aiaa.org/doi/abs/10.2514/1.28949)
[Sarabandi]Sarabandi, S. et al. (2018) Accurate Computation of Quaternions from Rotation Matrices. (http://www.iri.upc.edu/files/scidoc/2068-Accurate-Computation-of-Quaternions-from-Rotation-Matrices.pdf)
[Sarkka]Särkkä, S. (2007) Notes on Quaternions (https://users.aalto.fi/~ssarkka/pub/quat.pdf)
[Shepperd]S.W. Shepperd. “Quaternion from rotation matrix.” Journal of Guidance and Control, Vol. 1, No. 3, pp. 223-224, 1978. (https://arc.aiaa.org/doi/10.2514/3.55767b)
[Shoemake]K. Shoemake. Uniform random rotations. Graphics Gems III, pages 124-132. Academic, New York, 1992.
[Sola](1, 2) Solà, Joan. Quaternion kinematics for the error-state Kalman Filter. October 12, 2017. (http://www.iri.upc.edu/people/jsola/JoanSola/objectes/notes/kinematics.pdf)
[WikiConversions]https://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles
[WikiQuaternion]https://en.wikipedia.org/wiki/Quaternion
[Wiki_SLERP]https://en.wikipedia.org/wiki/Slerp