Quaternion
Quaternions are a number system which can be used to represent rotations in 3D space.
Double quaternion allows representing 4D rotations using two quaternions.
Dual quaternion allows represent both rotations and translations in 3D space.
The algebra of quaternions, double quaternions, and dual quaternions is called Clifford algebras.
Background
Quaternions are a number space that can be represented as a 4D point \(\displaystyle q = (q_0, q_1, q_2, q_3) = (q_0, \mathbf{q})\) with unit norm.
Here, \(\displaystyle \mathbf{q}\) represents the imaginary part: \(\displaystyle q_1 \mathbf{i}+q_2 \mathbf{j}+q_3 \mathbf{k}\).
The conjugate is \(\displaystyle \bar{q} = (q_0, -\mathbf{q})\).
Multiplication
Multiplication is \(\displaystyle q*p = (q_0 p_0 - \mathbf{q} \cdot \mathbf{p}, q_0 \mathbf{p} + p_0 \mathbf{q} + \mathbf{q} \times \mathbf{p})\).
Similar to imaginary numbers, \(\displaystyle \mathbf{i}^2=\mathbf{j}^2=\mathbf{k}^2=\mathbf{i}\mathbf{j}\mathbf{k}=-1\).
Multiplication between bases follow cross product rules: \(\displaystyle \mathbf{i} * \mathbf{j} = \mathbf{k}\) and \(\displaystyle \mathbf{j} * \mathbf{i} = -\mathbf{k}\).
In general, quaternion multiplication does not commute: \(\displaystyle q_1 * q_2 \neq q_2 * q_1\)
Intuition: Left multiplication by a quaternion is a 4D rotation plus a scaling operation.
Rotations
Quadratically conjugating a quaternion as follows is equivalent to applying a rotation on \(\displaystyle \mathbf{x}=(x,y,z)\):
\[\mathbf{q} * (c,x,y,z)*\bar{\mathbf{q}} = (c, R(\mathbf{q}) \cdot \mathbf{x})\]
Here, \(\displaystyle R(q)=R(-q)\) is a two-to-one mapping from quaternions to 3x3 rotation matrices.
Quaternion multiplication corresponds to composing rotations:
\(\displaystyle R(q*p) = R(q) \cdot R(p)\)
The quaternion \(\displaystyle q=(\cos(\frac{\theta}{2}), \hat{n}_1 \sin(\frac{\theta}{2}), \hat{n}_2 \sin(\frac{\theta}{2}), \hat{n}_3 \sin(\frac{\theta}{2}))=\cos(\frac{\theta}{2}) + \sin(\frac{\theta}{2}) \mathbf{\hat{n}}\) is equivalent to the rotation around axis \(\displaystyle \hat{n}\) by angle \(\displaystyle \theta\).
Slurp
Slerp, or spherical linear interpolation is done as: \[ \operatorname{slerp}(q_0, q_1, s) \equiv q(s)[q_0, q_1] = q_0 \frac{\sin((1-s)\phi)}{\sin \phi} + q_1 \frac{\sin(s\phi)}{\sin \phi} \] Here \(\displaystyle \phi\) is the angle between the two unit quaternions: \(\displaystyle \cos(\phi) = q_0 \cdot q_1\).
Distance Metrics
There are multiple ways to calculate the distance between two quaternions. See Stackexchange.
- \(\displaystyle \theta = \cos^{-1}( \langle q_1, q_2 \rangle)\)
- Note that this loss can be misleading for rotations because of double cover:
- \(\displaystyle R(q_1) = R(-q_1)\) but \(\displaystyle \cos^{-1}(x) \neq \cos^{-1}(-x)\)
- Thus some people use the distance below.
- \(\displaystyle \theta = \cos^{-1}(2 \langle q_1, q_2 \rangle^2 - 1)\)
Spatial Alignment Problem
Also known as the RMSD problem.
Let \(\displaystyle \{y_k\}\) be a data array with N columns of D-dimensional points, known as the reference structure.
Let \(\displaystyle \{x_k\}\) be a data array of matched points, known as the test structure.
The goal is to rotation \(\displaystyle \{x_k\}\) by a rotation matrix \(\displaystyle R_{D} \in SO(D)\) to minimize the mean squared distance:
\[
\mathbf{S}_D = \sum_{k=1}^{N} \Vert R_D \cdot x_k - y_x \Vert^2
\]
The paper by Hanson focuses on 3D points.
Cross-term maximization
It can be shown that the least squares minimization problem can be converted to a cross-term maximization:
\(\displaystyle
\begin{aligned}
\mathbf{S}_D &= \sum_{k=1}^{N} \Vert R_D x_k - y_k \Vert^2\\
&= \sum_{k=1}^{N} (R_D x_k - y_x)^T(R_D x_k - y_k)\\
&= \sum_{k=1}^{N} ( x_k ^T R_D^T - y_k^T)(R_D x_k - y_k)\\
&= \sum_{k=1}^{N} ( x_k ^T R_D^T - y_k^T)(R_D x_k - y_k)\\
&= \sum_{k=1}^{N} - 2(R_D x_k) \cdot y_k\\
\end{aligned}
\)
So minimization of \(\displaystyle \mathbf{S}_D\) is equivalent to maximization of \(\displaystyle \sum_{k=1}^{N} (R_D x_k) \cdot y_k = tr R_d \cdot E\) where \(\displaystyle E_{ab} = \sum_{k=1}^{N} [x_k]_a [y_k]_b = [X \cdot Y^T]_{ab}\). This is also equivalent to solving \(\displaystyle \delta(q) = q M(E) q^t\).
Algebraic Solutions
There are expressions for solving the eigenvalues of \(\displaystyle M\).
Orientation-Frame Alignment Problem
An orientation frame are the columns of the \(\displaystyle SO(3)\) orientation matrix.
In this problem, we assume we have \(\displaystyle n\) paired orientation frames, \(\displaystyle \{p_k\}\) and \(\displaystyle \{r_k\}\) which we want to align using quaternion \(\displaystyle q\).
Geodesic arc-length distance
For two unit quaternions \(\displaystyle q_1\) and \(\displaystyle q_2\), the geodesic arc length is \(\displaystyle \alpha\) where \(\displaystyle q_1 \cdot q_2 = \cos \alpha\).
Because of sign ambiguity, \(\displaystyle R(q_2) = R(-q_2)\), we take the minimum value:
\[
d_{geodesic}(q_1, q_2) = \min(\alpha, \pi-\alpha) \text{ where }
0 \leq d_{geodesic}(q_1, q_2) \leq \frac{\pi}{2}
\]
The least-squares optimization is: \[ \begin{align*} \mathbf{S}_{geodesic} &= \sum_{k=1}^{N}(\operatorname{arccos} |(q*p_k) \cdot r_k|)^2\\ &= \sum_{k=1}^{N}(\operatorname{arccos} |q \cdot (r_k * \bar{p}_k)|)^2 \end{align*} \]
Chord distance
The geodesic least squares cannot be solved using linear algebra.
Instead you can use chord distance:
\(\displaystyle \begin{aligned} d_{\text{chord}}(q_1, q_2) &= \min(\Vert q_1 - q_2 \Vert, \Vert q_1 + q_2 \Vert)\\ & 0 \leq d_{\text{chord}}(q_1, q_2) \leq \sqrt{2} \end{aligned} \)
In this case, our least squares problem is: \[ \mathbf{S}_{\text{chord}} = \sum_{k=1}^{N}\left( \min( \Vert (q*p_k) - r_k \Vert, \Vert (q*p_k)+r_k \Vert ) \right)^2 \]
Double Quaternions
Note that double quaternions are different from dual quaternions.
Double quaternions are written as \(\displaystyle q_1 + \epsilon q_2\) with \(\displaystyle \epsilon^2 = 1\) and applied as \(\displaystyle q_1 p q_2\).
Dual quaternions are written as \(\displaystyle q_1 + \varepsilon q_2\) with \(\displaystyle \varepsilon^2 = 0\) and applied as \(\displaystyle (q_1 + \varepsilon q_2) p (\bar{q}_1 + \varepsilon \bar{q}_2)\)
Cayley Factorization
See [Federico Thomas].
Any 4D rotation matrix can be decomposed into a right and a left isoclinic rotation matrix:
\(\displaystyle R = R^L R^R = R^R R^L\)
\(\displaystyle R^R\) and \(\displaystyle R^L\) can be viewed as a matrix representation of single double quaternion \(\displaystyle (q_1, q_2)\).
For a double quaternion, the 4D rotation written is \(\displaystyle x' = q_1 x q_2\).
The product of left and right isoclinic rotation matrices commute.
Furthermore, the product of two left isoclinic rotation matrices is a left isoclinic rotation matrix. Same with right.
Thus, \(\displaystyle R_1 R_2 = (R_1^L R_1^R) (R_2^L R_2^R) = (R_1^L R_2^L) (R_1^R R_2^R) = R\).
This shows that the composition of two double quaternions will be a double quaternion.
Approximating 3D Translations using Double quaternions
Dual Quaternions
Combined Point + Frame Alignment Problem
Given a spatial profile matrix \(\displaystyle M(E)\) and orientation frame profile matrix \(\displaystyle U(S)\), we can normalize the scale to have unit-eigenvalue and then apply linear interpolation with dimensional constant \(\displaystyle \sigma\): \[ \Delta_{xf}(t, \sigma) = q \cdot \left[ (1-t) \frac{M(E)}{\epsilon_{x}} + t \sigma \frac{U(S)}{\epsilon_{f}} \right] \cdot q \]
Or alternatively, we can use slerp: \[ q(t) = \operatorname{slerp}(q_{x:opt}, q_{f:opt}, t) \]
Resources
- The Quaternion-Based Spatial Coordinate and Orientation Frame Alignment Problems by Andrew Hanson Edited Paper
- Approaching Dual Quaternions From Matrix Algebra by Federico Thomas
- Double Quaternions for Motion Interpolation by Q.J. Ge, Amitabh Varshney, Jai P. Menon, Chu-Fei Chang