Every robot that carries an IMU has to answer one question many times a second: which way am I facing? There are three common ways to write the answer down: a rotation matrix, three Euler angles, or a unit quaternion. This post builds all three from scratch. It starts in the plane, where complex numbers already know how to rotate, and then follows the same idea into three dimensions until quaternions fall out on their own.

Along the way we will meet the single most useful habit for this topic: reading the same matrix in two ways, as a rotation that moves a vector (active) and as a change of coordinates between two frames (passive). The outline follows Ben-Ari's tutorial [1].

1. Rotations in the plane

1.1 Rotating a vector

Take a vector \(\mathbf{v} = (x, y)\) of length \(r\) that makes an angle \(\varphi\) with the \(x\)-axis, so that \(x = r\cos\varphi\) and \(y = r\sin\varphi\). Rotating it counterclockwise by \(\theta\) keeps the length and adds \(\theta\) to the angle.

xy vr v′r φ θ
Figure 1. Rotating \(\mathbf{v}\) by \(\theta\) adds \(\theta\) to its angle and keeps its length.

With the angle-sum formulas,

\[ \begin{aligned} x' &= r\cos(\varphi + \theta) = r\cos\varphi\cos\theta - r\sin\varphi\sin\theta = x\cos\theta - y\sin\theta, \\ y' &= r\sin(\varphi + \theta) = r\sin\varphi\cos\theta + r\cos\varphi\sin\theta = x\sin\theta + y\cos\theta, \end{aligned} \]

or, as a matrix,

\[ \begin{equation} \begin{bmatrix} x' \\ y' \end{bmatrix} = \underbrace{\begin{bmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{bmatrix}}_{R(\theta)} \begin{bmatrix} x \\ y \end{bmatrix}. \label{eq:R2} \end{equation} \]

1.2 The columns are the rotated axes

There is a faster way to find \(R(\theta)\), and it is the one we will use in 3D. A matrix times the basis vector \(\mathbf{e}_1 = (1, 0)\) gives its first column, and times \(\mathbf{e}_2 = (0, 1)\) gives its second column. So we only need to know where the two axes go:

xy e1e2 Re1 Re2 cos θsin θ −sin θcos θ θ θ
Figure 2. Rotating the basis vectors: \(R\mathbf{e}_1 = (\cos\theta, \sin\theta)\) and \(R\mathbf{e}_2 = (-\sin\theta, \cos\theta)\) are the two columns of \(R(\theta)\).
\[ R\mathbf{e}_1 = \begin{bmatrix} \cos\theta \\ \sin\theta \end{bmatrix}, \qquad R\mathbf{e}_2 = \begin{bmatrix} -\sin\theta \\ \cos\theta \end{bmatrix} \quad\Longrightarrow\quad R(\theta) = \begin{bmatrix} R\mathbf{e}_1 & R\mathbf{e}_2 \end{bmatrix}. \]

A few properties follow directly and carry over to 3D:

  • The columns are orthonormal, so \(R^\top R = I\) and the inverse is just the transpose: \(R(\theta)^{-1} = R(-\theta) = R(\theta)^\top\).
  • \(\det R = \cos^2\theta + \sin^2\theta = 1\): a rotation never mirrors.
  • Rotations add: \(R(\alpha)R(\beta) = R(\alpha + \beta)\). In particular, 2D rotations commute. This will fail in 3D.

1.3 Complex numbers and Euler's formula

Now encode the vector \((x, y)\) as the complex number \(z = x + iy\). (Engineers usually write \(j\) for the imaginary unit; I use \(i\) here because \(j\) will get a different job in Section 4.) The key fact is Euler's formula:

\[ \begin{equation} e^{i\theta} = \cos\theta + i\sin\theta. \label{eq:euler} \end{equation} \]

It follows from the Taylor series of \(e^x\). Because \(i^2 = -1\), the even powers of \(i\theta\) are real and alternate in sign, and the odd powers are imaginary:

\[ e^{i\theta} = \sum_{n=0}^{\infty} \frac{(i\theta)^n}{n!} = \underbrace{\left(1 - \frac{\theta^2}{2!} + \frac{\theta^4}{4!} - \cdots\right)}_{\cos\theta} + i\underbrace{\left(\theta - \frac{\theta^3}{3!} + \frac{\theta^5}{5!} - \cdots\right)}_{\sin\theta}. \]
ReIm1i eiθ cos θi sin θ θ
Figure 3. \(e^{i\theta}\) is the point on the unit circle at angle \(\theta\).

Writing \(z\) in polar form, \(z = r e^{i\varphi}\), a rotation becomes a single multiplication:

\[ e^{i\theta} z = r\,e^{i(\varphi + \theta)}. \]

The length stays \(r\) and the angle grows by \(\theta\), which is exactly Figure 1. Expanding the product gives back the matrix \(\eqref{eq:R2}\):

\[ (\cos\theta + i\sin\theta)(x + iy) = (x\cos\theta - y\sin\theta) + i\,(x\sin\theta + y\cos\theta). \]

Three ideas from this section drive everything that follows:

  1. A vector is a number built from an imaginary unit: \((x, y) \leftrightarrow x + iy\).
  2. A rotation is multiplication by a number of length one, \(e^{i\theta}\). Composing rotations is multiplying numbers.
  3. The imaginary unit generates the rotation. Multiplying by \(i\) itself rotates by \(90^\circ\), and \(e^{i\theta}\) is “\(i\) applied a fraction \(\theta / (\pi/2)\) of the way”.

2. Rotations and frames

2.1 Axis, angle and the right-hand rule

In 3D, every rotation turns space about some axis, a unit vector \(\mathbf{n}\), by some angle \(\theta\). (This is Euler's rotation theorem.) To say which way is positive we use the right-hand rule: point the thumb of your right hand along \(\mathbf{n}\); your fingers curl in the positive direction. Equivalently, looking from the tip of \(\mathbf{n}\) back toward the origin, a positive rotation is counterclockwise.

(a) right-handed frame: e₁ × e₂ = e₃ x y z ω = θ̇ n (b) positive rotation about n: thumb along n
Figure 4. (a) In a right-handed frame, curling from \(x\) to \(y\) points the thumb along \(z\). (b) The angular velocity \(\boldsymbol{\omega}\) points along the rotation axis, and its length is the rotation rate.

The same rule fixes the frames themselves: a frame is right-handed when \(\mathbf{e}_1 \times \mathbf{e}_2 = \mathbf{e}_3\). It also defines the angular velocity \(\boldsymbol{\omega} = \dot\theta\,\mathbf{n}\): its direction is the axis, and its length is how fast the body turns. A gyroscope measures exactly this vector, in the coordinates of the sensor's own frame.

A 3D rotation matrix has orthonormal columns and determinant \(+1\): \(R^\top R = I\), \(\det R = 1\). As in 2D, its inverse is its transpose.

2.2 Two ways to read a rotation matrix

The equation \(\mathbf{y} = R\mathbf{x}\) can mean two different things, and mixing them up is the source of most sign and order bugs.

xy x x′ = Rx θ (a) active: one frame, the vector moves xWyW xByB p θ (b) passive: one point, two frames
Figure 5. (a) Active: \(R\) moves a vector within one frame. (b) Passive: the point \(p\) stays put, and \(R_{WB}\) converts its coordinates from the body frame \(B\) to the world frame \(W\).

Active reading (a rotation). There is one frame. \(\mathbf{x}\) and \(\mathbf{x}' = R\mathbf{x}\) are coordinates in that same frame, and \(R\) moves the vector to a new place. This is what Section 1 did.

Passive reading (a frame transformation). There is one point \(p\) and two frames: a world frame \(W\), with axes \(\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3\), and a body frame \(B\) attached to the robot, with axes \(\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3\). The body coordinates \(\mathbf{x}_B\) of the point say how to reach it along the body axes:

\[ p = x_{B,1}\,\mathbf{b}_1 + x_{B,2}\,\mathbf{b}_2 + x_{B,3}\,\mathbf{b}_3. \]

To get world coordinates, write each body axis in world coordinates, \([\mathbf{b}_k]_W\), and collect them as columns:

\[ \begin{equation} \mathbf{x}_W = x_{B,1}[\mathbf{b}_1]_W + x_{B,2}[\mathbf{b}_2]_W + x_{B,3}[\mathbf{b}_3]_W = \underbrace{\begin{bmatrix} [\mathbf{b}_1]_W & [\mathbf{b}_2]_W & [\mathbf{b}_3]_W \end{bmatrix}}_{R_{WB}} \mathbf{x}_B . \label{eq:passive} \end{equation} \]

So the columns of \(R_{WB}\) are the body axes written in world coordinates. Its rows are the world axes written in body coordinates, because \(R_{BW} = R_{WB}^{-1} = R_{WB}^\top\). Nothing moves here: \(\mathbf{x}_B\) and \(\mathbf{x}_W\) describe the same point, and \(R_{WB}\) only rewrites its coordinates.

How the two readings connect. Apply \(R_{WB}\) actively to the world axes: \(R_{WB}\mathbf{e}_k\) is the \(k\)-th column, which is \([\mathbf{b}_k]_W\). So \(R_{WB}\) is also the rotation that carries the world frame onto the body frame. “The body is rotated by \(R_{WB}\)” and “\(\mathbf{x}_W = R_{WB}\mathbf{x}_B\)” are the same matrix, read two ways.

The subscripts make the passive reading almost mechanical. \(R_{AB}\) takes \(B\)-coordinates to \(A\)-coordinates, and adjacent subscripts cancel:

\[ \mathbf{x}_A = R_{AB}\,\mathbf{x}_B, \qquad R_{AC} = R_{AB}\,R_{BC}, \qquad R_{BA} = R_{AB}^{-1} = R_{AB}^\top . \]

Example. A robot has turned \(90^\circ\) to the left, so \(R_{WB} = R(90^\circ) = \left[\begin{smallmatrix} 0 & -1 \\ 1 & 0 \end{smallmatrix}\right]\) (in the plane). A point one meter straight ahead of the robot has \(\mathbf{x}_B = (1, 0)\), so \(\mathbf{x}_W = R_{WB}\mathbf{x}_B = (0, 1)\): one meter along the world \(y\)-axis, which is where the robot is now facing.

Why the passive reading is the useful one. On a robot, sensors measure in the body frame and goals live in the world frame. A gyroscope gives \(\boldsymbol{\omega}_B\). An accelerometer at rest measures the reaction to gravity in body coordinates, \(\mathbf{a}_B = R_{WB}^\top (0, 0, g)^\top\), and comparing that with the known world vector is how a filter such as Mahony's corrects its attitude. Every one of these steps is a change of frame, not a motion.

3. Euler angles

3.1 Rotations about one axis

We build the three elementary rotations the way we built \(R(\theta)\) in Section 1.2: find where each axis goes, and use the results as columns. In each case, look down the rotation axis from its tip, so that the rotation is counterclockwise and the other two axes behave like the 2D plane.

yz x Rxe2 Rxe3 φ Rx(φ) zx y Rye3 Rye1 θ Ry(θ) xy z Rze1 Rze2 ψ Rz(ψ)
Figure 6. Each elementary rotation, seen from the tip of its axis (pointing at you). The plane is \((y, z)\) for \(R_x\), \((z, x)\) for \(R_y\), and \((x, y)\) for \(R_z\), always in the cyclic order \(x \to y \to z \to x\).

About \(x\) by the roll angle \(\varphi\). The \(x\)-axis does not move. In the \((y, z)\) plane, \(y\) plays the role of the 2D \(x\)-axis and \(z\) the role of the 2D \(y\)-axis:

\[ R_x\mathbf{e}_1 = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix},\quad R_x\mathbf{e}_2 = \begin{bmatrix} 0 \\ \cos\varphi \\ \sin\varphi \end{bmatrix},\quad R_x\mathbf{e}_3 = \begin{bmatrix} 0 \\ -\sin\varphi \\ \cos\varphi \end{bmatrix}, \]

so

\[ R_x(\varphi) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos\varphi & -\sin\varphi \\ 0 & \sin\varphi & \cos\varphi \end{bmatrix}. \]

About \(z\) by the yaw angle \(\psi\). This is the 2D rotation with \(z\) left alone:

\[ R_z\mathbf{e}_1 = \begin{bmatrix} \cos\psi \\ \sin\psi \\ 0 \end{bmatrix},\quad R_z\mathbf{e}_2 = \begin{bmatrix} -\sin\psi \\ \cos\psi \\ 0 \end{bmatrix},\quad R_z\mathbf{e}_3 = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}, \]

so

\[ R_z(\psi) = \begin{bmatrix} \cos\psi & -\sin\psi & 0 \\ \sin\psi & \cos\psi & 0 \\ 0 & 0 & 1 \end{bmatrix}. \]

About \(y\) by the pitch angle \(\theta\). Here the cyclic order matters. After \(y\) comes \(z\), then \(x\), so in the \((z, x)\) plane it is \(z\) that plays the role of the 2D \(x\)-axis. A positive rotation turns \(z\) toward \(x\), and \(x\) toward \(-z\):

\[ R_y\mathbf{e}_3 = \begin{bmatrix} \sin\theta \\ 0 \\ \cos\theta \end{bmatrix},\quad R_y\mathbf{e}_1 = \begin{bmatrix} \cos\theta \\ 0 \\ -\sin\theta \end{bmatrix},\quad R_y\mathbf{e}_2 = \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix}, \]

so

\[ R_y(\theta) = \begin{bmatrix} \cos\theta & 0 & \sin\theta \\ 0 & 1 & 0 \\ -\sin\theta & 0 & \cos\theta \end{bmatrix}. \]

The minus sign sits below the diagonal in \(R_y\) but above it in \(R_x\) and \(R_z\). That is not a typo: listing the rows in the order \(x, y, z\) writes the \((z, x)\) plane backwards. (A side effect: with \(z\) pointing up, a positive pitch tips the nose, the \(x\)-axis, down.)

3.2 The order matters

Unlike in the plane, rotations about different axes do not commute. Follow \(\mathbf{e}_1\) through two quarter turns:

\[ R_z(90^\circ)\,R_x(90^\circ)\,\mathbf{e}_1 = R_z(90^\circ)\,\mathbf{e}_1 = \mathbf{e}_2, \qquad R_x(90^\circ)\,R_z(90^\circ)\,\mathbf{e}_1 = R_x(90^\circ)\,\mathbf{e}_2 = \mathbf{e}_3 . \]

Same two rotations, different results. So “three Euler angles” means nothing until we also fix the order of the axes. There are 12 possible orders (any sequence of three axes with no axis twice in a row).

3.3 Yaw, pitch, roll: \(R_{WB} = R_z(\psi)\,R_y(\theta)\,R_x(\varphi)\)

The usual choice in robotics and aerospace, and the one in my IMU code, is the Z-Y-X order: yaw \(\psi\) about \(z\), pitch \(\theta\) about \(y\), roll \(\varphi\) about \(x\),

\[ \begin{equation} R_{WB} = R_z(\psi)\,R_y(\theta)\,R_x(\varphi). \label{eq:zyx} \end{equation} \]

Note the order of the factors: \(z\) on the left, \(x\) on the right. Writing \(R_x R_y R_z\) gives a different (and much less common) convention. The two readings from Section 2.2 explain the order.

Active reading: rotations about the fixed world axes. Applied to a vector, \(R_z(\psi)\big(R_y(\theta)(R_x(\varphi)\, \mathbf{v})\big)\) acts right to left: first roll about the world \(x\)-axis, then pitch about the world \(y\)-axis, then yaw about the world \(z\)-axis. The axes stay fixed, and the order is \(x\), then \(y\), then \(z\).

Passive reading: a chain of frames. Insert two intermediate frames, \(1\) and \(2\), and read \(\eqref{eq:zyx}\) as the chain rule of Section 2.2:

\[ \mathbf{x}_W = \underbrace{R_z(\psi)}_{R_{W1}}\ \underbrace{R_y(\theta)}_{R_{12}}\ \underbrace{R_x(\varphi)}_{R_{2B}}\ \mathbf{x}_B . \]

Each factor is an elementary rotation about an axis of the frame before it. Frame 1 is the world frame yawed by \(\psi\) about \(z_W\). Frame 2 is frame 1 pitched by \(\theta\) about its own \(y_1\)-axis. The body frame is frame 2 rolled by \(\varphi\) about its own \(x_2\)-axis. The axes move with the body, and the order is \(z\), then \(y\), then \(x\). This is how we naturally describe a robot or a plane: heading first, then tilt, then bank.

xW yW zW start: world frame W x1 y1 z1 yaw ψ about zW x2 y2 z2 pitch θ about y1 xB yB zB roll φ about x2
Figure 7. The passive (moving-axes) reading of \(R_{WB} = R_z(\psi)R_y(\theta)R_x(\varphi)\). The orange axis is the one being rotated about. Dashed lines show the world frame.

Both readings describe the same matrix: rotating about the moving axes in the order \(z, y, x\) equals rotating about the fixed axes in the order \(x, y, z\). Multiplying out, with \(c_\varphi = \cos\varphi\), \(s_\varphi = \sin\varphi\) and so on:

\[ \begin{equation} R_{WB} = \begin{bmatrix} c_\psi c_\theta & c_\psi s_\theta s_\varphi - s_\psi c_\varphi & c_\psi s_\theta c_\varphi + s_\psi s_\varphi \\ s_\psi c_\theta & s_\psi s_\theta s_\varphi + c_\psi c_\varphi & s_\psi s_\theta c_\varphi - c_\psi s_\varphi \\ -s_\theta & c_\theta s_\varphi & c_\theta c_\varphi \end{bmatrix}. \label{eq:Rzyx} \end{equation} \]

The angles can be read back from the entries \(r_{ij}\). The bottom-left entry gives the pitch, and the ratios in the last row and the first column give roll and yaw:

\[ \theta = -\arcsin r_{31}, \qquad \varphi = \operatorname{atan2}(r_{32},\, r_{33}), \qquad \psi = \operatorname{atan2}(r_{21},\, r_{11}). \]

3.4 Gimbal lock

The read-back fails when \(\cos\theta = 0\), i.e. pitch \(\theta = \pm 90^\circ\). Then \(r_{11} = r_{21} = r_{32} = r_{33} = 0\) and both \(\operatorname{atan2}\) calls become \(\operatorname{atan2}(0, 0)\). Substituting \(\theta = 90^\circ\) into \(\eqref{eq:Rzyx}\) shows why:

\[ R_{WB}\big|_{\theta = 90^\circ} = \begin{bmatrix} 0 & \sin(\varphi - \psi) & \cos(\varphi - \psi) \\ 0 & \cos(\varphi - \psi) & -\sin(\varphi - \psi) \\ -1 & 0 & 0 \end{bmatrix}. \]

The matrix depends only on \(\varphi - \psi\). Every pair \((\varphi, \psi)\) with the same difference gives the same orientation, so one degree of freedom has disappeared. (At \(\theta = -90^\circ\) it is \(\varphi + \psi\) instead.)

xW yW zW (a) pitch θ = 0: three independent axes xB yB zB yaw axis zW roll axis xB (b) pitch θ = 90°: xB points along −zW
Figure 8. At \(\theta = 90^\circ\) the body \(x\)-axis (the roll axis) points along \(-z_W\), the yaw axis. Yawing and rolling now turn the body about the same line.

Geometrically, pitching by \(90^\circ\) swings the roll axis onto the yaw axis, so the two rotations do the same thing. Near this pitch, small changes in orientation cause huge jumps in \(\varphi\) and \(\psi\). Nothing is wrong with the orientation itself; the problem is the parametrization. In fact, every description of 3D rotations by three numbers has singular points like this one. Using four numbers avoids it, which brings us to quaternions.

4. Quaternions

4.1 From the plane to space: what must \(i, j, k\) be?

Let us redo Section 1.3 in three dimensions without assuming anything about quaternions. In the plane, one imaginary unit was enough. In space, give each axis its own unit and write a vector as

\[ \mathbf{v} = (x, y, z) \ \longleftrightarrow\ v = x\,i + y\,j + z\,k . \]

We want to multiply such numbers, and we want the rules for \(i, j, k\) to be forced by a few reasonable requirements.

(R1) Each unit is an imaginary unit: \(i^2 = j^2 = k^2 = -1\), just like \(i\) in the plane.

(R2) Squaring a vector gives its length. For complex numbers, \(z\bar z = |z|^2\). Here we ask for \(v^2 = -|\mathbf{v}|^2\), so every unit vector is a square root of \(-1\): no axis is special. Expanding, while keeping the order of factors because we have not assumed commutativity,

\[ v^2 = -(x^2 + y^2 + z^2) + xy\,(ij + ji) + yz\,(jk + kj) + zx\,(ki + ik). \]

For this to equal \(-(x^2 + y^2 + z^2)\) for every vector, the units must anticommute:

\[ ji = -ij, \qquad kj = -jk, \qquad ik = -ki . \]

(R3) What is \(ij\)? It cannot be a real number, because \((ij)^2 = ijij = -iijj = -(-1)(-1) = -1\), and no real number squares to \(-1\). Can it be a combination \(a + b\,i + c\,j\) of what we already have? Multiply \(ij = a + bi + cj\) on the left by \(i\). The left side becomes \(i^2 j = -j\), and the right side becomes

\[ ai - b + c\,(ij) = ai - b + c\,(a + bi + cj) = (ca - b) + (a + cb)\,i + c^2\,j . \]

Matching the \(j\) parts needs \(c^2 = -1\), which is impossible for a real \(c\). So \(ij\) is genuinely new. We call it \(k\): \(k = ij\). This is why the 3D version of a complex number needs four components, \(1, i, j, k\), not three.

The remaining products now follow from \(k = ij\) and anticommutation:

\[ \begin{aligned} jk &= j\,(ij) = (ji)\,j = -(ij)\,j = -i\,(jj) = i, \\ ki &= (ij)\,i = i\,(ji) = -i\,(ij) = -(ii)\,j = j, \\ k^2 &= (ij)(ij) = i\,(ji)\,j = -i\,(ij)\,j = -(ii)(jj) = -1 . \end{aligned} \]

Altogether this is Hamilton's multiplication table:

\[ i^2 = j^2 = k^2 = ijk = -1, \qquad ij = k,\ \ jk = i,\ \ ki = j, \qquad ji = -k,\ \ kj = -i,\ \ ik = -j . \]
i j k ij = k ji = −k
Figure 9. Multiplying two units along the arrows gives the third with a plus sign; against the arrows, with a minus sign.

A quaternion is a number \(q = w + x\,i + y\,j + z\,k\). We write it as a scalar part plus a vector part, \(q = (s, \mathbf{u})\) with \(s = w\) and \(\mathbf{u} = (x, y, z)\). A quaternion with \(s = 0\) is a pure quaternion, which is how vectors live inside the quaternions.

4.2 Multiplication with dot and cross products

Compare the table with the basis vectors \(\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3\). The square \(i^2 = -1\) looks like \(-\mathbf{e}_1\cdot\mathbf{e}_1\), and \(ij = k = -ji\) looks like \(\mathbf{e}_1 \times \mathbf{e}_2 = \mathbf{e}_3 = -\mathbf{e}_2 \times \mathbf{e}_1\). In general, the product of two units is minus their dot product plus their cross product. By linearity, the same holds for any two pure quaternions:

\[ \begin{equation} \mathbf{u}\,\mathbf{v} = -\,\mathbf{u}\cdot\mathbf{v} + \mathbf{u}\times\mathbf{v} . \label{eq:pure} \end{equation} \]

Adding the scalar parts, the product of two general quaternions is

\[ \begin{equation} (s, \mathbf{u})\,(t, \mathbf{v}) = \big(st - \mathbf{u}\cdot\mathbf{v},\ \ s\,\mathbf{v} + t\,\mathbf{u} + \mathbf{u}\times\mathbf{v}\big). \label{eq:product} \end{equation} \]

This is much easier to remember than the sixteen-term expansion, and it shows at a glance why quaternion multiplication does not commute: the cross product changes sign when the factors swap. (Vince [2] develops all of quaternion algebra in this vector form.) In components it matches the usual formula:

\[ \begin{aligned} (q\,p)_w &= q_w p_w - q_x p_x - q_y p_y - q_z p_z, \\ (q\,p)_x &= q_w p_x + q_x p_w + q_y p_z - q_z p_y, \\ (q\,p)_y &= q_w p_y + q_y p_w + q_z p_x - q_x p_z, \\ (q\,p)_z &= q_w p_z + q_z p_w + q_x p_y - q_y p_x . \end{aligned} \]

Some further facts we will need, each a one-line check with \(\eqref{eq:product}\):

  • The conjugate is \(q^* = (s, -\mathbf{u})\), and \((qp)^* = p^*q^*\).
  • The norm satisfies \(qq^* = q^*q = s^2 + |\mathbf{u}|^2 = |q|^2\) and \(|qp| = |q|\,|p|\).
  • The inverse is \(q^{-1} = q^*/|q|^2\). For a unit quaternion (\(|q| = 1\)), \(q^{-1} = q^*\).
  • For a unit vector \(\mathbf{n}\), \(\eqref{eq:pure}\) gives \(\mathbf{n}^2 = -1\). The Taylor-series argument of Section 1.3 then repeats word for word: \(e^{\mathbf{n}\alpha} = \cos\alpha + \mathbf{n}\sin\alpha\). Every unit vector is an “\(i\)” with its own Euler formula.

4.3 Rotating a vector, and the half angle

First attempt. Copy the plane: multiply \(\mathbf{v}\) by \(q = e^{\mathbf{n}\alpha} = \cos\alpha + \mathbf{n}\sin\alpha\). By \(\eqref{eq:product}\),

\[ q\,\mathbf{v} = \big(-\sin\alpha\ \mathbf{n}\cdot\mathbf{v},\ \ \cos\alpha\,\mathbf{v} + \sin\alpha\,\mathbf{n}\times\mathbf{v}\big). \]

If \(\mathbf{v} \perp \mathbf{n}\), the scalar part vanishes and the vector part, \(\cos\alpha\,\mathbf{v} + \sin\alpha\,\mathbf{n}\times \mathbf{v}\), is exactly \(\mathbf{v}\) rotated by \(\alpha\) about \(\mathbf{n}\), because \(\mathbf{n}\times\mathbf{v}\) is \(\mathbf{v}\) turned \(90^\circ\) about \(\mathbf{n}\). That is the plane all over again. But if \(\mathbf{v}\) has a component along \(\mathbf{n}\), the result has a scalar part and is no longer a vector.

The fix: the sandwich \(q\,\mathbf{v}\,q^*\). Split \(\mathbf{v} = \mathbf{v}_\parallel + \mathbf{v}_\perp\) into the parts along and across \(\mathbf{n}\), and use \(\eqref{eq:pure}\) to see how each part commutes with \(\mathbf{n}\):

\[ \mathbf{n}\,\mathbf{v}_\parallel = -\mathbf{n}\cdot\mathbf{v}_\parallel = \mathbf{v}_\parallel\,\mathbf{n}, \qquad \mathbf{n}\,\mathbf{v}_\perp = \mathbf{n}\times\mathbf{v}_\perp = -\,\mathbf{v}_\perp\times\mathbf{n} = -\,\mathbf{v}_\perp\,\mathbf{n} . \]

So \(\mathbf{v}_\parallel\) commutes with \(q = \cos\alpha + \mathbf{n}\sin\alpha\), while \(\mathbf{v}_\perp\) flips the sign of the \(\mathbf{n}\) term as it passes through: \(\mathbf{v}_\perp\,q^* = q\,\mathbf{v}_\perp\). Therefore

\[ \begin{aligned} q\,\mathbf{v}_\parallel\,q^* &= \mathbf{v}_\parallel\,q\,q^* = \mathbf{v}_\parallel, \\ q\,\mathbf{v}_\perp\,q^* &= q\,q\,\mathbf{v}_\perp = e^{2\alpha\mathbf{n}}\,\mathbf{v}_\perp = \cos 2\alpha\ \mathbf{v}_\perp + \sin 2\alpha\ \mathbf{n}\times\mathbf{v}_\perp . \end{aligned} \]

The part along the axis is untouched, and the part across it turns by \(2\alpha\). That is a rotation by \(2\alpha\) about \(\mathbf{n}\), and the result is a pure quaternion again. To rotate by \(\theta\), take \(\alpha = \theta/2\):

n v∥ v v′ v⊥ n × v⊥ v′⊥ θ
Figure 10. Rotating \(\mathbf{v}\) about \(\mathbf{n}\): \(\mathbf{v}_\parallel\) stays, \(\mathbf{v}_\perp\) turns by \(\theta\) toward \(\mathbf{n}\times\mathbf{v}_\perp\).
\[ \begin{equation} \boxed{\ q = e^{\mathbf{n}\theta/2} = \Big(\cos\frac{\theta}{2},\ \ \sin\frac{\theta}{2}\,\mathbf{n}\Big), \qquad \mathbf{v}' = q\,\mathbf{v}\,q^* = \mathbf{v}_\parallel + \cos\theta\ \mathbf{v}_\perp + \sin\theta\ \mathbf{n}\times\mathbf{v}\ } \label{eq:sandwich} \end{equation} \]

(Note \(\mathbf{n}\times\mathbf{v} = \mathbf{n}\times\mathbf{v}_\perp\). The right-hand side is Rodrigues' rotation formula.) The half angle appears because \(q\) acts twice, once from each side. Two consequences:

  • \(q\) and \(-q\) give the same rotation, since the two minus signs cancel in \(q\,\mathbf{v}\,q^*\). Every rotation has exactly two unit quaternions.
  • Every unit quaternion is a rotation: any \(q\) with \(|q| = 1\) can be written as \((\cos\frac\theta2, \sin\frac\theta2\,\mathbf{n})\).

Back to the plane. Take \(\mathbf{n} = k\) and a vector in the \(xy\)-plane, \(\mathbf{v} = x\,i + y\,j \perp k\). Then \(q\,\mathbf{v}\,q^* = q^2\,\mathbf{v} = (\cos\theta + k\sin\theta)(x\,i + y\,j)\), and with \(ki = j\) and \(kj = -i\),

\[ (\cos\theta + k\sin\theta)(x\,i + y\,j) = (x\cos\theta - y\sin\theta)\,i + (x\sin\theta + y\cos\theta)\,j , \]

which is the 2D rotation \(\eqref{eq:R2}\). The quaternion \(k\) plays exactly the role that \(i\) played in the complex plane.

Example. Rotate \(\mathbf{e}_1 = i\) by \(90^\circ\) about \(z\). Here \(q = \cos 45^\circ + k\sin 45^\circ = \tfrac{\sqrt2}{2}(1 + k)\). Since \(i \perp k\), \(q\,i\,q^* = q^2\,i = k\,i = j\), so the result is \(\mathbf{e}_2\), as it should be.

4.4 Composition and frames

Quaternions compose like matrices. Rotating by \(q_1\) and then by \(q_2\):

\[ q_2\,(q_1\,\mathbf{v}\,q_1^*)\,q_2^* = (q_2 q_1)\,\mathbf{v}\,(q_2 q_1)^* , \]

so the combined rotation is \(q_2 q_1\), in the same order as \(R_2 R_1\). Everything from the passive reading carries over. Let \(q_{WB}\) be the quaternion of \(R_{WB}\). Then

\[ \mathbf{x}_W = q_{WB}\,\mathbf{x}_B\,q_{WB}^*, \qquad q_{WC} = q_{WB}\,q_{BC}, \qquad q_{BW} = q_{WB}^* . \]

From Euler angles. Each elementary rotation is a rotation about a basis axis, so by \(\eqref{eq:sandwich}\), \(q_x(\varphi) = \cos\frac\varphi2 + i\sin\frac\varphi2\), \(q_y(\theta) = \cos\frac\theta2 + j\sin\frac\theta2\), and \(q_z(\psi) = \cos\frac\psi2 + k\sin\frac\psi2\). Multiplying them in the same order as \(\eqref{eq:zyx}\), \(q_{WB} = q_z(\psi)\,q_y(\theta)\,q_x(\varphi)\), and writing \(c_\varphi = \cos\frac\varphi2\), \(s_\varphi = \sin\frac\varphi2\) (half angles this time):

\[ \begin{aligned} w &= c_\varphi c_\theta c_\psi + s_\varphi s_\theta s_\psi, & x &= s_\varphi c_\theta c_\psi - c_\varphi s_\theta s_\psi, \\ y &= c_\varphi s_\theta c_\psi + s_\varphi c_\theta s_\psi, & z &= c_\varphi c_\theta s_\psi - s_\varphi s_\theta c_\psi . \end{aligned} \]

Integrating a gyroscope. The passive reading also tells us how to update the orientation from gyroscope readings. Over a short time \(\Delta t\), the body turns by the small angle \(|\boldsymbol\omega|\Delta t\) about \(\boldsymbol\omega_B\), expressed in its own frame. That small rotation is \(q_{BB'} \approx 1 + \tfrac12\boldsymbol\omega_B\Delta t\), and the chain rule puts it on the right: \(q_{WB'} = q_{WB}\,q_{BB'}\). Letting \(\Delta t \to 0\) gives the quaternion kinematics used in attitude filters [3]:

\[ \dot q_{WB} = \tfrac12\, q_{WB} \otimes (0,\ \boldsymbol\omega_B). \]

4.5 From a quaternion to a rotation matrix

Expanding the sandwich with \(\eqref{eq:product}\) gives a formula that needs no trigonometry. For a unit \(q = (s, \mathbf{u})\), first \(q\,\mathbf{v} = (-\mathbf{u}\cdot\mathbf{v},\ s\mathbf{v} + \mathbf{u}\times\mathbf{v})\). Multiplying on the right by \(q^* = (s, -\mathbf{u})\), the scalar part cancels and the vector part is

\[ (\mathbf{u}\cdot\mathbf{v})\,\mathbf{u} + s^2\,\mathbf{v} + 2s\,\mathbf{u}\times\mathbf{v} - (\mathbf{u}\times\mathbf{v})\times\mathbf{u} . \]

With the identity \((\mathbf{u}\times\mathbf{v})\times\mathbf{u} = |\mathbf{u}|^2\mathbf{v} - (\mathbf{u}\cdot\mathbf{v})\,\mathbf{u}\),

\[ \begin{equation} q\,\mathbf{v}\,q^* = (s^2 - |\mathbf{u}|^2)\,\mathbf{v} + 2\,(\mathbf{u}\cdot\mathbf{v})\,\mathbf{u} + 2s\,\mathbf{u}\times\mathbf{v} . \label{eq:qvq} \end{equation} \]

Each term is linear in \(\mathbf{v}\), so this is a matrix: \(R = (s^2 - |\mathbf{u}|^2)\,I + 2\,\mathbf{u}\mathbf{u}^\top + 2s\,[\mathbf{u}]_\times\), where \([\mathbf{u}]_\times\) is the matrix of \(\mathbf{v} \mapsto \mathbf{u}\times\mathbf{v}\). With \(q = (w, x, y, z)\) and \(w^2 + x^2 + y^2 + z^2 = 1\):

\[ \begin{equation} R(q) = \begin{bmatrix} 1 - 2(y^2 + z^2) & 2(xy - wz) & 2(xz + wy) \\ 2(xy + wz) & 1 - 2(x^2 + z^2) & 2(yz - wx) \\ 2(xz - wy) & 2(yz + wx) & 1 - 2(x^2 + y^2) \end{bmatrix}. \label{eq:Rq} \end{equation} \]

Going back is just as direct. The trace is \(\operatorname{tr}R = 3 - 4(x^2 + y^2 + z^2) = 4w^2 - 1\), and the off-diagonal differences isolate the vector part:

\[ w = \tfrac12\sqrt{1 + \operatorname{tr}R}, \qquad x = \frac{r_{32} - r_{23}}{4w}, \qquad y = \frac{r_{13} - r_{31}}{4w}, \qquad z = \frac{r_{21} - r_{12}}{4w} . \]

This works as long as \(w\) is not close to zero. For rotations near \(180^\circ\), robust code starts from the largest diagonal entry instead.

4.6 From a quaternion to Euler angles

Now match \(\eqref{eq:Rq}\) against \(\eqref{eq:Rzyx}\) and reuse the read-back formulas from Section 3.3. With \(r_{31} = 2(xz - wy)\), \(r_{32} = 2(yz + wx)\), and so on:

\[ \begin{aligned} \text{roll: }\ \varphi &= \operatorname{atan2}\big(2(wx + yz),\ 1 - 2(x^2 + y^2)\big), \\ \text{pitch: }\ \theta &= \arcsin\big(2(wy - xz)\big), \\ \text{yaw: }\ \psi &= \operatorname{atan2}\big(2(wz + xy),\ 1 - 2(y^2 + z^2)\big). \end{aligned} \]

In floating point, clamp the argument of \(\arcsin\) to \([-1, 1]\). Gimbal lock still appears here, as \(2(wy - xz) = \pm 1\), but only in the angles. The quaternion itself is perfectly well behaved there, which is why filters store a quaternion and convert to Euler angles only for display.

5. Summary

Representation Numbers Composition Singularities Good for
Rotation matrix \(R\) 9 (with 6 constraints) \(R_2 R_1\) none transforming many vectors
Euler angles \((\psi, \theta, \varphi)\) 3 through matrices gimbal lock at \(\theta = \pm 90^\circ\) reading and setting by humans
Unit quaternion \(q\) 4 (with 1 constraint) \(q_2 q_1\) none (\(q\) and \(-q\) are the same) integration, filtering, storage

The thread through all of it:

  • In the plane, \(e^{i\theta}\) rotates by multiplication. In space, \(e^{\mathbf{n}\theta/2}\) rotates by the sandwich \(q\,\mathbf{v}\,q^*\).
  • The columns of a rotation matrix are the rotated axes. The columns of \(R_{WB}\) are the body axes in world coordinates.
  • Read \(R_{WB}\) and \(q_{WB}\) passively, and let the subscripts cancel: \(R_{AC} = R_{AB}R_{BC}\), \(q_{AC} = q_{AB}\,q_{BC}\).

References

[1] Moti Ben-Ari. “A Tutorial on Euler Angles and Quaternions.” Weizmann Institute of Science, version 2.0.2, 2014–2017.

[2] John Vince. Quaternions for Computer Graphics. Springer, 2011.

[3] Joan Solà. “Quaternion kinematics for the error-state Kalman filter.” arXiv:1711.02508, 2017.