"Give me three angles and I will eventually lose a degree of freedom at the worst possible moment. Give me four numbers on a sphere and I will never look down and find the floor missing."
A Well-Oriented AI Agent
Prerequisites
This section assumes the angular-rate signal \(\boldsymbol\omega\) and the body-frame idea from Chapter 23: a gyroscope reports rotational velocity in the sensor's own axes, and to speak about the world we must relate that moving frame to a fixed one. Vectors, matrices, and the cross product from the math reference in Part I are enough; no prior exposure to quaternions is needed. The integration-and-drift consequences of the kinematics derived here are developed in Section 24.2 and Section 24.4. Python and NumPy suffice for the code.
The Big Picture
Everything in this chapter, complementary filters, attitude and heading reference systems, dead reckoning, learned inertial odometry, rests on one object: the orientation of a body in space, and how it changes when a gyroscope hands you an angular rate. Choose the wrong way to write that orientation down and you inherit singularities, expensive re-orthogonalization, or interpolation that visibly stutters. The unit quaternion is the representation that estimation code overwhelmingly settles on, not for elegance but because it integrates cleanly, composes cheaply, never suffers gimbal lock, and stays valid with a single division. Master the four-number object here and the rest of Part VI becomes bookkeeping over rotations rather than a fight with the coordinate system.
Four ways to write down a rotation
A rigid-body orientation is a single geometric fact, but it has several algebraic clothes, and each trades away something. A rotation matrix \(R \in SO(3)\) is a \(3\times 3\) orthogonal matrix with determinant \(+1\); it rotates a vector by plain multiplication, \(\mathbf v' = R\,\mathbf v\), and composes by matrix product. It is unambiguous and fast to apply, but it spends nine numbers on three degrees of freedom, and after many multiplications floating-point error makes it drift away from orthogonality, forcing a periodic and costly re-orthogonalization. Euler angles (roll, pitch, yaw) use the minimum three numbers and read intuitively, which is exactly why they persist in user interfaces and datasheets. They carry a fatal defect for estimation: at a pitch of \(\pm 90^\circ\) two of the three axes align and one degree of freedom vanishes, a singularity called gimbal lock, where the orientation is still well defined but its rate becomes undefined and small motions demand infinite angular speed in the parameters.
The axis-angle form states the geometric truth directly: any orientation is a rotation by angle \(\theta\) about a unit axis \(\hat{\mathbf n}\) (Euler's rotation theorem). It has no gimbal lock but composes awkwardly. The unit quaternion is the axis-angle written so that composition and integration become algebra. It is a four-vector \[ q = \begin{bmatrix} w \\ \mathbf v \end{bmatrix} = \begin{bmatrix} \cos(\theta/2) \\ \hat{\mathbf n}\,\sin(\theta/2) \end{bmatrix}, \qquad \lVert q \rVert = 1, \] where \(w\) is the scalar part and \(\mathbf v \in \mathbb{R}^3\) the vector part. The half-angle is not decoration: it is what makes the quaternion double-cover the rotation group, so \(q\) and \(-q\) name the identical physical orientation, a subtlety we return to.
Key Insight
Pick the representation by the operation you do most. If you mostly apply a fixed rotation to many vectors, cache a rotation matrix; the multiply is one dense \(3\times 3\). If you mostly estimate and integrate orientation from a gyroscope, hundreds of times a second, use quaternions: they update with a small multiply, renormalize with one division, and cannot hit a singularity. If a human must read or set the value, convert to Euler angles at the very last moment, at the boundary of your system, and never let them into the filter state. Most orientation bugs in inertial code trace to Euler angles leaking inward.
The quaternion algebra you actually use
Three operations carry almost all the weight. Composition of rotations is the Hamilton product. For \(q_1=(w_1,\mathbf v_1)\) and \(q_2=(w_2,\mathbf v_2)\), \[ q_1 \otimes q_2 = \begin{bmatrix} w_1 w_2 - \mathbf v_1 \cdot \mathbf v_2 \\ w_1 \mathbf v_2 + w_2 \mathbf v_1 + \mathbf v_1 \times \mathbf v_2 \end{bmatrix}, \] which applies \(q_2\) first, then \(q_1\), just like stacking matrices. Because it contains a cross product, quaternion multiplication does not commute, correctly reflecting that rotations do not either. Rotating a vector \(\mathbf u\) embeds it as a pure quaternion \((0,\mathbf u)\) and conjugates: \(\mathbf u' = q \otimes (0,\mathbf u) \otimes q^{-1}\), and for a unit quaternion the inverse is just the conjugate \(q^{-1}=(w,-\mathbf v)\). Inversion is therefore free. These three, multiply, conjugate, rotate, are the entire working vocabulary; you rarely need anything else in an estimator.
The operation that makes quaternions the natural home for inertial work is kinematics: how \(q\) evolves under a measured angular rate \(\boldsymbol\omega\) (body frame, from your gyroscope in Chapter 23). The continuous law is beautifully compact, \[ \dot q = \tfrac{1}{2}\, q \otimes (0,\boldsymbol\omega), \] a single quaternion product, with none of the trigonometric case-splitting that Euler-angle rates require near their singularity. Over one sample of duration \(\Delta t\) the exact update is a rotation by \(\lVert\boldsymbol\omega\rVert\,\Delta t\) about \(\boldsymbol\omega/\lVert\boldsymbol\omega\rVert\), which you build as a small quaternion \(\Delta q\) and compose. This exact-exponential step is why a quaternion integrator stays numerically clean where a matrix integrator slowly forgets it is a rotation.
import numpy as np
def quat_mul(a, b):
"""Hamilton product a (x) b, quaternions as [w, x, y, z]."""
w1, x1, y1, z1 = a
w2, x2, y2, z2 = b
return np.array([
w1*w2 - x1*x2 - y1*y2 - z1*z2,
w1*x2 + x1*w2 + y1*z2 - z1*y2,
w1*y2 - x1*z2 + y1*w2 + z1*x2,
w1*z2 + x1*y2 - y1*x2 + z1*w2,
])
def integrate_gyro(q, omega, dt):
"""Advance orientation q by body-frame rate omega (rad/s) over dt.
Uses the exact axis-angle exponential, robust as ||omega|| -> 0."""
theta = np.linalg.norm(omega) * dt # rotation angle this step
if theta < 1e-9: # avoid 0/0 for tiny rates
dq = np.array([1.0, *(0.5 * omega * dt)])
else:
axis = omega / np.linalg.norm(omega)
dq = np.array([np.cos(theta/2), *(axis * np.sin(theta/2))])
q_new = quat_mul(q, dq) # body-frame update: q (x) dq
return q_new / np.linalg.norm(q_new) # one division keeps it valid
# 200 Hz gyro, a steady 30 deg/s yaw for 3 seconds -> expect ~90 deg heading.
q = np.array([1.0, 0.0, 0.0, 0.0]) # identity orientation
omega = np.deg2rad([0.0, 0.0, 30.0])
for _ in range(600):
q = integrate_gyro(q, omega, 1/200)
yaw = np.rad2deg(2 * np.arctan2(q[3], q[0])) # yaw from a pure-z rotation
print(f"integrated yaw = {yaw:.2f} deg") # ~ 90.00
The printout lands on roughly \(90^\circ\), as a steady \(30^\circ/\mathrm{s}\) yaw for three seconds should. Note what the loop does not need: no special case at pitch \(90^\circ\), no re-orthogonalization, no trigonometric unwinding. That economy, one multiply and one divide per sample, is why quaternion integration runs comfortably on a microcontroller sampling an IMU at hundreds of hertz, a point we lean on again in the edge-deployment discussion of Chapter 61.
Interpolation, the double cover, and staying on the sphere
Two quaternion pitfalls bite in practice. First, interpolation. Blending two orientations by averaging their four components and renormalizing (nlerp) is cheap but moves at non-uniform angular speed and can swing the long way around. The correct constant-rate blend is spherical linear interpolation, slerp: \[ \mathrm{slerp}(q_0,q_1,t) = \frac{\sin((1-t)\Omega)}{\sin\Omega}\,q_0 + \frac{\sin(t\Omega)}{\sin\Omega}\,q_1, \qquad \cos\Omega = q_0\cdot q_1, \] which traces the shortest great-circle arc at uniform angular velocity, exactly what you want to resample an orientation trajectory or to smooth a camera path. Second, the double cover: because \(q\) and \(-q\) are the same rotation, \(q_0\cdot q_1\) can be negative even when the orientations are close. Before slerping (or before comparing, or averaging) flip one sign so the dot product is non-negative; skip this and slerp will take the \(340^\circ\) detour instead of the \(20^\circ\) one. The same sign ambiguity is why you never compute orientation error as a raw component difference: use \(q_{\text{err}} = q_{\text{ref}}^{-1}\otimes q_{\text{est}}\) and read the angle from its scalar part, a habit that carries straight into the trajectory metrics of Section 24.7.
In Practice: the VR headset that snapped through the floor
A standalone VR headset fused its gyroscope into a quaternion orientation at 1 kHz and streamed it to the renderer at 90 Hz, interpolating between fused samples to hit exact frame times. On fast head turns some users reported a brief, sickening snap where the world spun almost fully around before settling. The cause was textbook double cover: the interpolator blended consecutive quaternions without the sign check, so whenever integration produced a sample near \(-q\) of its predecessor, the slerp chose the long arc. The fix was two lines, negate the incoming quaternion when its dot product with the previous one was negative, before interpolating. No filter change, no added latency; the artifact vanished. It is the canonical reminder that a quaternion pipeline needs sign hygiene as much as it needs the right math.
Right Tool: scipy.spatial.transform.Rotation
Everything above, Hamilton products, vector rotation, conversions to and from matrices, Euler angles, and axis-angle, plus a numerically careful slerp, ships in SciPy. Rotating a batch of vectors is Rotation.from_quat(q).apply(vectors); converting for a human is .as_euler('xyz', degrees=True); interpolating a trajectory is Slerp(times, rotations)(query_times). That collapses the roughly 40 lines of hand-written product, normalization, sign-fix, and conversion code a full pipeline accumulates into three or four calls, and SciPy uses the scalar-last \([x,y,z,w]\) convention, so confirm the ordering when you hand quaternions across a library boundary. Hand-roll the algebra once to understand it (as above), then let the library carry it in production.
Research Frontier: rotations that neural networks can regress
Quaternions are ideal for integrating a gyroscope, but they are a poor output target for a neural network: the double cover and the unit-norm constraint make the map from network activations to rotations discontinuous, and regression to it trains badly. Zhou et al. (2019), "On the Continuity of Rotation Representations in Neural Networks," showed that quaternions and Euler angles are discontinuous representations for learning, and proposed a continuous 6D representation (two 3-vectors passed through Gram-Schmidt to form a rotation matrix) that networks regress far more accurately. This choice is now standard in learned pose and orientation heads, including the learned inertial odometry models of Section 24.6 and the pose estimators of Chapter 27. The lesson: use quaternions in the filter, but reach for a continuous parameterization at a network's output layer.
When each representation earns its place
Put the choices in one frame. Store and integrate orientation state as a unit quaternion: it is singularity-free, cheap to propagate from angular rate, and trivially renormalized, which is why the filters of Section 24.2 and the quaternion error-state Kalman filters that build on Chapter 9 carry a quaternion in their state. Apply orientation with a rotation matrix when you transform many vectors at once, for example projecting a lidar sweep into a world frame in Chapter 52. Expose orientation to people as Euler angles only at the interface. Regress orientation from a network with a continuous representation. One geometric fact, four algebraic clothes, and the engineering skill is matching the garment to the task rather than forcing one everywhere.
Exercise
Extend the integrator above. (1) Feed it a time-varying rate that ramps yaw up and back down, and confirm the net orientation returns to identity, so that \(q\) reads \([1,0,0,0]\) up to sign. (2) Add a constant gyroscope bias of \(0.5^\circ/\mathrm{s}\) on one axis and integrate for ten minutes; plot the growing orientation error and connect the slope to the drift analysis you will meet in Section 24.4. (3) Implement slerp with the sign fix, then deliberately remove the sign fix and construct two nearby orientations whose interpolation takes the long way around; measure the peak angular speed in both cases.
Self-Check
1. Why does the half-angle appear in the quaternion definition, and what does the resulting \(q\) versus \(-q\) ambiguity mean physically? 2. A device pitches through vertical (\(90^\circ\)) during a maneuver. Explain why an Euler-angle filter state can blow up there while a quaternion state does not. 3. You must interpolate between two orientations whose quaternion dot product is \(-0.98\). What is the shortest-arc angle between them, and what one step must you take before calling slerp?
What's Next
In Section 24.2, we turn the clean quaternion kinematics of this section into a working estimator. Integrating the gyroscope alone drifts (a preview of Section 24.4), so we blend that smooth-but-drifting orientation with the absolute-but-noisy tilt and heading references from the accelerometer and magnetometer. The complementary filter and the Mahony and Madgwick filters are exactly that fusion, expressed as small, elegant corrections applied to the quaternion state you now know how to propagate.