Quaternion Quick Facts ====================== This is not meant to be a treatise on quaternions, but a summary relevant for the code implementation. It also assumes familiarity with linear algebra. See also the `Matrix and Quaternions FAQ `_ for a practical compendium of equations for implementing matrix and quaternion operations. Quaternions ----------- The space of quaternions is a vector space generated by the basis :math:`{\boldsymbol{1}, \boldsymbol{i}, \boldsymbol{j}, \boldsymbol{k}}`. In this space, quaternions are the sum of a *scalar* part generated by :math:`{\boldsymbol{1}}` and a vector component generated by :math:`{\boldsymbol{i}, \boldsymbol{j}, \boldsymbol{k}}`. A multiplication operation is defined over the space such that: * it is distributive with respect to the usual sum * the identity is the element :math:`\boldsymbol{1}` * the products of the basis elements is given by :math:`{\boldsymbol{i}, \boldsymbol{j}, \boldsymbol{k}}`: .. math:: :label: quaternion_ops \begin{eqnarray} i^2 &=& j^2 = k^2 = -1 \\ i j &=& - j i = k \\ j k &=& - k j = i \\ k i &=& - i k = j \end{eqnarray} **Example**: The product of two pure vectors :math:`\boldsymbol{v} = \vec{v} = v_0 \boldsymbol{i} + v_1 \boldsymbol{j} + v_2 \boldsymbol{k}` and :math:`\vec{u} = u_0 \boldsymbol{i} + u_1 \boldsymbol{j} + u_2 \boldsymbol{k}` is: .. math:: \begin{eqnarray} \boldsymbol{v} \boldsymbol{u} &=& - (v_0 u_0 + v_1 u_1 + v_2 u_2)\boldsymbol{1} + (v_1 u_2 - v_2 u_1)\boldsymbol{i} + (v_2 u_0 - v_0 u_2)\boldsymbol{j} + (v_0 u_1 - v_1 u_0)\boldsymbol{k} \\ &=& - \vec{v} \cdot \vec{u} + \vec{v} \times \vec{u} \end{eqnarray} Rotations ---------- Most often, 3-d rotations are represented as :math:`3\times3` matrices .. math:: :label: matrix_def \begin{eqnarray} M &=& \left[\begin{matrix} M_{00} & M_{01} & M_{02} \\ M_{10} & M_{11} & M_{12} \\ M_{20} & M_{21} & M_{22} \end{matrix}\right] \end{eqnarray} Some properties of rotations: * they are orthogonal: :math:`M^{-1} = M^T` or :math:`M^T M = M M^T = \mathbb{1}` * det(M) = 1 In what follows we use the following elemental (right-handed) rotations of :math:`\alpha` around the axes :math:`R_{\alpha}^z`: .. math:: :label: elemental_rotations \begin{eqnarray} R_{\alpha}^x &=& \left[\begin{matrix} 1 & 0 & 0 \\ 0 & \cos\alpha & -\sin\alpha \\ 0 & \sin\alpha & \cos\alpha \end{matrix}\right] \\ R_{\alpha}^y &=& \left[\begin{matrix} \cos\alpha & 0 & \sin\alpha \\ 0 & 1 & 0 \\ -\sin\alpha & 0 & \cos\alpha \end{matrix}\right] \\ R_{\alpha}^z &=& \left[\begin{matrix} \cos\alpha & -\sin\alpha & 0 \\ \sin\alpha & \cos\alpha & 0 \\ 0 & 0 & 1 \end{matrix}\right] \end{eqnarray} Angle-axis Representation ^^^^^^^^^^^^^^^^^^^^^^^^^ A 3-d rotation can be uniquely defined by a unit vector and an angle :math:`(\vec{u}, \phi)`. * there is one and only one real-valued eigenvalue: :math:`\lambda = 1`. * the rotation axis is given by a unit eigenvector of M with eigenvalue :math:`\lambda = 1`. * the rotation angle is given by the trace of M: :math:`Tr(M) = 1 + 2 \cos \phi`. Equatorial Representation ^^^^^^^^^^^^^^^^^^^^^^^^^ Any rotation can be decomposed as a sequence of three rotations around different cartesian axes. Using the equatorial coordinates related to a *fixed* intertial reference frame, this can be expressed as :math:`R_{\alpha}^z R_{-\delta}^y R_{\rho}^x`, where :math:`{\rho}`, :math:`{\delta}` and :math:`{\alpha}` are roll, declination and right ascension respectively. Quaternion Representation ^^^^^^^^^^^^^^^^^^^^^^^^^ A rotation can be represented by a quaternion operator: .. math:: :label: quaternion_operator L_q(\boldsymbol{v}) = \boldsymbol{q}^* \boldsymbol{v} \boldsymbol{q}, where :math:`\boldsymbol{q}` is a unit quaternion (i.e.: :math:`\boldsymbol{v}^2 = \boldsymbol{1}`). The quaternion is related to the rotation axis and angle: .. math:: :label: quaternion_axis_angle \boldsymbol{q} = \cos \frac{\phi}{2} + \sin \frac{\phi}{2}\vec{u} Switching Representations ------------------------- .. _equatorialmatrix: Equatorial -> Matrix ^^^^^^^^^^^^^^^^^^^^ Let's write :math:`M = R_{\alpha}^z R_{-\delta}^y R_{\rho}^x` in matrix form: .. math:: :label: eq_matrix_repr \begin{eqnarray} M &=& \left[\begin{matrix}\cos{\alpha} & - \sin{\alpha} & 0\\\sin{\alpha} & \cos{\alpha} & 0\\0 & 0 & 1\end{matrix}\right] \left[\begin{matrix}\cos{\delta} & 0 & - \sin{\delta}\\0 & 1 & 0\\ \sin{\delta} & 0 & \cos{\delta}\end{matrix}\right] \left[\begin{matrix}1 & 0 & 0\\0 & \cos{\rho} & - \sin{\rho}\\0 & \sin{\rho} & \cos{\rho}\end{matrix}\right]\\ &=& \left[\begin{matrix}\cos{\alpha} \cos{\delta} & - \sin{\alpha} \cos{\rho} - \sin{\delta} \sin{\rho} \cos{\alpha} & \sin{\alpha} \sin{\rho} - \sin{\delta} \cos{\alpha} \cos{\rho}\\\sin{\alpha} \cos{\delta} & - \sin{\alpha} \sin{\delta} \sin{\rho} + \cos{\alpha} \cos{\rho} & - \sin{\alpha} \sin{\delta} \cos{\rho} - \sin{\rho} \cos{\alpha}\\ \sin{\delta} & \sin{\rho} \cos{\delta} & \cos{\delta} \cos{\rho}\end{matrix}\right] \end{eqnarray} Matrix -> Equatorial ^^^^^^^^^^^^^^^^^^^^ From Eq. :eq:`eq_matrix_repr` we can invert and get: .. math:: \begin{eqnarray} \tan \alpha &=& \frac{M_{10}}{M_{00}} \\ \tan \rho &=& \frac{M_{21}}{M_{22}} \\ \sin \delta &=& - M_{20} \end{eqnarray} Quaternion -> Matrix ^^^^^^^^^^^^^^^^^^^^ Expanding Eq. :eq:`quaternion_operator` using the multiplication rules in Eq. :eq:`quaternion_ops`, and then grouping factors of :math:`v_*`, we can write this in matrix form: .. math:: :label: quat_matrix_repr \begin{eqnarray} M &=& \left[\begin{matrix} 1 - 2(q_1^2 + q_2^2) & 2(q_0 q_1 - q_2 q_3) & 2(q_0 q_2 + q_1 q_3) \\2(q_0 q_1 + q_2 q_3) &1 - 2(q_0^2 + q_2^2) &2(q_1 q_2 - q_0 q_3) \\2(q_0 q_2 - q_1 q_3) &2(q_1 q_2 + q_0 q_3) &1 - 2(q_0^2 + q_1^2)\end{matrix}\right] \end{eqnarray} For practical reasons, some times we want to express the transform in terms of a quaternion with norm :math:`\left\| \boldsymbol{q} \right\| \neq 1`, only requiring that :math:`\left\| \boldsymbol{q} \right\| \neq 0`: .. math:: :label: quat_matrix_repr_2 \begin{eqnarray} M &=& \frac{1}{\left\| q \right\|^2} \left[\begin{matrix} \left\| q \right\|^2 - 2(q_1^2 + q_2^2) & 2(q_0 q_1 - q_2 q_3) & 2(q_0 q_2 + q_1 q_3) \\ 2(q_0 q_1 + q_2 q_3) & \left\| q \right\|^2 - 2(q_0^2 + q_2^2) &2(q_1 q_2 - q_0 q_3) \\ 2(q_0 q_2 - q_1 q_3) & 2(q_1 q_2 + q_0 q_3) & \left\| q \right\|^2 - 2(q_0^2 + q_1^2) \end{matrix}\right] \end{eqnarray} Matrix -> Quaternion ^^^^^^^^^^^^^^^^^^^^ This is a bit more involved. The idea is to invert from Eq. :eq:`quat_matrix_repr_2`. From the diagonal elements we can already get the squares of the quaternion components: .. math:: :label: quat_matrix_den S = \frac{\left\| q \right\|^2}{4} \left[\begin{matrix} 1 + M_{00} - M_{11} - M_{22} \\ 1 - M_{00} + M_{11} - M_{22} \\ 1 - M_{00} - M_{11} + M_{22} \\ 1 + M_{00} + M_{11} + M_{22} \end{matrix}\right] = \left[\begin{matrix} q_0^2 \\ q_1^2 \\ q_2^2 \\ q_3^2 \end{matrix}\right], The next step depends on which entry in Eq. :eq:`quat_matrix_den` is the largest: .. math:: :label: transform2quat \begin{eqnarray} 0 &\rightarrow& \left[ \sqrt{S_0}, \frac{(M_{01} + M_{10})}{4 \sqrt{S_0}}, \frac{(M_{02} + M_{20})}{4 \sqrt{S_0}}, \frac{(M_{21} - M_{12})}{4 \sqrt{S_0}} \right] \\ 1 &\rightarrow& \left[ \frac{(M_{01} + M_{10})}{4 \sqrt{S_1}}, \sqrt{S_1}, \frac{(M_{12} + M_{21})}{4 \sqrt{S_1}}, \frac{(M_{02} - M_{20})}{4 \sqrt{S_1}} \right] \\ 2 &\rightarrow& \left[ \frac{(M_{20} + M_{02})}{4 \sqrt{S_2}}, \frac{(M_{12} + M_{21})}{4 \sqrt{S_2}}, \sqrt{S_2}, \frac{(M_{10} - M_{01})}{4 \sqrt{S_2}} \right] \\ 3 &\rightarrow& \left[ \frac{(M_{21} - M_{12})}{4 \sqrt{S_3}}, \frac{(M_{02} - M_{20})}{4 \sqrt{S_3}}, \frac{(M_{10} - M_{01})}{4 \sqrt{S_3}}, \sqrt{S_3} \right] \end{eqnarray} Note that the denominator is always :math:`4 \sqrt{S_i}`, so we always choose the case with the largest denominator in order to minimize round-off errors. **Derivation**. As an example, we derive one of the cases in Eq. :eq:`transform2quat`. The rest are similar. The idea is to scale the quaternion so its largest component is equal to 1. In other words, if the i-th entry in Eq. :eq:`quat_matrix_den` is the largest, then multiply q by a factor :math:`1/S_i` to make :math:`q_i = 1`. The magnitude of the resulting quaternion is :math:`\left\| \boldsymbol{q} \right\| = 1/q_i`. If :math:`S_1 = 1 - M_{00} + M_{11} - M_{22}` is the largest, we scale q by :math:`1/\sqrt{S_1}`. The matrix in Eq. :eq:`quat_matrix_repr_2` then takes a simpler form: .. math:: \begin{eqnarray} M &=& S_1 \left[\begin{matrix} 1/S_1 - 2(1 + q_2^2) & 2(q_0 - q_2 q_3) & 2(q_0 q_2 + q_3) \\ 2(q_0 + q_2 q_3) & 1/S_1 - 2(q_0^2 + q_2^2) &2(q_2 - q_0 q_3) \\ 2(q_0 q_2 - q_3) & 2(q_2 + q_0 q_3) & 1/S_1 - 2(q_0^2 + 1)\end{matrix}\right] \end{eqnarray} and we can set: .. math:: \begin{eqnarray} q_0 &=& \frac{(M_{01} + M_{10})}{4 S_1} \\ q_1 &=& 1 \\ q_2 &=& \frac{(M_{12} + M_{21})}{4 S_1} \\ q_3 &=& \frac{(M_{02} - M_{20})}{4 S_1} \\ \left\| q \right\| &=& \frac{1}{\sqrt{S_1}} \end{eqnarray} and after normalizing: .. math:: q = \left[ \frac{(M_{01} + M_{10})}{4 \sqrt{S_1}}, \sqrt{S_1}, \frac{(M_{12} + M_{21})}{4 \sqrt{S_1}}, \frac{(M_{02} - M_{20})}{4 \sqrt{S_1}} \right] Equatorial -> Quaternion ^^^^^^^^^^^^^^^^^^^^^^^^ This code first transforms from equatorial to matrix, and then from matrix to quaternion. Quaternion -> Equatorial ^^^^^^^^^^^^^^^^^^^^^^^^ Transform from quaternion to matrix, and then from matrix to equatorial.