Rotation Representations

A 3D rotation has three degrees of freedom, yet the common ways of writing one down use three, four, six or nine numbers, each with its own failure mode. The extrinsic rotation $R$ of the camera parameters can be stored in any of them; optimizers, interpolators and neural networks each prefer a different one, and mixed-up conventions cause silent errors.

Definition

A rotation is a linear map of 3D space that preserves lengths, angles and handedness. A rotation representation is a set of parameters together with a rule that turns them into that map. The main ones are:

• the rotation matrix, a $3 \times 3$ orthogonal matrix with determinant $+1$;
• Euler angles, three successive rotations about coordinate axes;
• the axis–angle pair and the rotation vector (axis times angle);
• the unit quaternion, four numbers on the unit sphere in $\mathbb{R}^4$;
• the 6D representation, two columns of the rotation matrix, used for learning;
• the Lie algebra $\mathfrak{so}(3)$, the space of small rotations used for optimization.

This article uses right-handed frames and writes quaternions as $(w, x, y, z)$, scalar first, except where a library’s own convention is described.

Intuition

A camera’s orientation can be described by the directions its image axes point in (the rows of $R$), by “turn left 30°, tilt down 20°, roll 5°” (Euler angles), or by “turn 40° about this axis” (axis–angle); every 3D rotation is a turn about some axis.

The difficulty is topological. The rotations form a closed, curved three-dimensional manifold, and Zhou et al. (2019) show that no representation in four or fewer dimensions is continuous in the sense they define. Every compact representation has a weak spot: gimbal lock for Euler angles, the 180° boundary for rotation vectors, the sign ambiguity for quaternions.

Formal Definition

Rotation matrices and SO(3)

The rotation matrices form the special orthogonal group

$$SO(3) = \{ R \in \mathbb{R}^{3 \times 3} : R^\top R = I,\ \det R = +1 \}.$$

Orthogonality preserves lengths and angles; $\det R = +1$ excludes reflections. The nine entries obey six independent constraints, leaving three degrees of freedom (Solà, 2017, arXiv v1, Table 1). Rotations compose by matrix multiplication, $R_2 R_1$ applies $R_1$ first, and the inverse is the transpose. Composition is not commutative.

Euler angles

Euler angles write $R$ as a product of three rotations about coordinate axes. For the yaw–pitch–roll sequence that REP 103 describes,

$$R(\psi, \theta, \phi) = R_z(\psi)\, R_y(\theta)\, R_x(\phi).$$

Read as intrinsic rotations, this turns about $z$, then about the new $y$, then about the newest $x$; read as extrinsic rotations about fixed axes, the same product applies $x$, then $y$, then $z$. Twelve axis sequences are possible (no axis may follow itself), each read intrinsically or extrinsically; REP 103 counts 24 conventions.

Gimbal lock. At $\theta = \pm 90°$ the first and third rotation axes coincide, so yaw and roll turn the camera about the same axis. For this sequence at $\theta = +90°$ the matrix depends only on $\psi - \phi$: one degree of freedom is lost, and the angles cannot be recovered uniquely. It follows that near the lock, small changes of orientation need large changes of angles, which makes Euler angles a poor choice of parameters for optimization.

Axis–angle and the Rodrigues formula

A rotation by angle $\theta$ about a unit axis $\mathbf{u}$ is stored compactly as the rotation vector $\mathbf{r} = \theta \mathbf{u}$. The Rodrigues formula converts it to a matrix; OpenCV’s documentation of Rodrigues (4.14.0 header source) states it as

$$R = \cos\theta\, I + (1 - \cos\theta)\, \mathbf{u}\mathbf{u}^\top + \sin\theta\, [\mathbf{u}]_\times,$$

where $[\mathbf{u}]_\times$ is the cross-product matrix, $[\mathbf{u}]_\times \mathbf{v} = \mathbf{u} \times \mathbf{v}$. The inverse uses $\sin\theta\, [\mathbf{u}]_\times = (R - R^\top)/2$ and $\cos\theta = (\operatorname{tr} R - 1)/2$. OpenCV’s cv::Rodrigues converts in both directions and optionally returns the Jacobian. In the 4.14.0 source, the matrix-to-vector direction first replaces the input by $UV^\top$ from its singular value decomposition, and it handles angles near 180° by a separate branch, because $R - R^\top$ vanishes there and no longer determines the axis.

Unit quaternions

A unit quaternion $q = (w, x, y, z)$ with $w^2 + x^2 + y^2 + z^2 = 1$ represents the rotation by $\theta$ about $\mathbf{u}$ as

$$q = \left( \cos\tfrac{\theta}{2},\ \mathbf{u} \sin\tfrac{\theta}{2} \right).$$

Writing $q = (w, \mathbf{v})$, the corresponding matrix is

$$R(q) = (w^2 - \mathbf{v}^\top\mathbf{v})\, I + 2\,\mathbf{v}\mathbf{v}^\top + 2w\,[\mathbf{v}]_\times .$$

Composition is the Hamilton product (Solà, 2017, arXiv v1, §1.2.2 and Table 1),

$$q_1 \otimes q_2 = \left( w_1 w_2 - \mathbf{v}_1^\top \mathbf{v}_2,\ \ w_1 \mathbf{v}_2 + w_2 \mathbf{v}_1 + \mathbf{v}_1 \times \mathbf{v}_2 \right), \qquad R(q_1 \otimes q_2) = R(q_1)\, R(q_2).$$

Because of the half angle, $q$ and $-q$ give the same matrix: the unit quaternions form a double cover of $SO(3)$ (Solà, 2017, arXiv v1, §2.4.6). Solà (2017, arXiv v1, §3) explains that quaternion conventions also differ in the algebra’s handedness and in the direction of the operator, not only in the component order. The formulas here are the Hamilton convention, which he reports Eigen, ROS and Ceres follow; the JPL convention, common in aerospace, is left-handed.

Slerp. Spherical linear interpolation moves along the great arc between $q_0$ and $q_1$ at constant angular speed:

$$q(t) = q_0 \frac{\sin\big((1 - t)\Omega\big)}{\sin\Omega} + q_1 \frac{\sin(t\Omega)}{\sin\Omega}, \qquad \cos\Omega = q_0^\top q_1 .$$

Solà (2017, arXiv v1, §2.7) gives this form, which he credits to Glenn Davis citing Shoemake (1985), and notes the double-cover condition: if $q_0^\top q_1 < 0$, replace $q_1$ by $-q_1$ first, or the interpolation takes the long way round.

6D representation for learning

Zhou et al. (2019) proposed a continuous representation for networks that regress rotations: the first two columns of $R$, six numbers. A network outputs any six numbers $(\mathbf{a}_1, \mathbf{a}_2)$, and a Gram–Schmidt step maps them to a rotation:

$$\mathbf{b}_1 = \frac{\mathbf{a}_1}{\lVert \mathbf{a}_1 \rVert}, \quad \mathbf{b}_2 = \frac{\mathbf{a}_2 - (\mathbf{b}_1^\top \mathbf{a}_2)\,\mathbf{b}_1}{\lVert \mathbf{a}_2 - (\mathbf{b}_1^\top \mathbf{a}_2)\,\mathbf{b}_1 \rVert}, \quad \mathbf{b}_3 = \mathbf{b}_1 \times \mathbf{b}_2, \quad R = [\,\mathbf{b}_1\ \mathbf{b}_2\ \mathbf{b}_3\,].$$

Their argument: $SO(3)$ is homeomorphic to the real projective space $\mathbb{RP}^3$, which does not embed in $\mathbb{R}^d$ for $d < 5$, so Euler angles, rotation vectors and quaternions all have discontinuities. In their autoencoder sanity test (Figure 5, table c, at 500k iterations), the 6D representation had a mean geodesic error of 0.49° and a maximum of 1.98°, against means of 3.32° for quaternions and 6.98° for Euler angles, each with maxima near 180°. In their point-cloud pose test (Figure 5, table f, at 2,600k iterations), the 6D mean was 2.85°, below quaternions (9.03°) and a directly regressed $3 \times 3$ matrix orthogonalized after prediction (4.21°), though its maximum error was still near 180° there.

Lie group and Lie algebra

$SO(3)$ is a Lie group. Its tangent space at the identity, the Lie algebra $\mathfrak{so}(3)$, consists of the skew-symmetric matrices $[\boldsymbol{\omega}]_\times$, identified with vectors $\boldsymbol{\omega} \in \mathbb{R}^3$. The exponential map sends a tangent vector to a rotation, and for $SO(3)$ its closed form is the Rodrigues formula; the logarithm inverts it (Solà et al., 2018, arXiv v9, Example 4 and Appendix B):

$$\operatorname{Exp}(\boldsymbol{\omega}) = \exp([\boldsymbol{\omega}]_\times), \qquad \operatorname{Log}(R) = \frac{\theta\, (R - R^\top)^\vee}{2 \sin\theta}, \quad \theta = \arccos\frac{\operatorname{tr} R - 1}{2},$$

where $^\vee$ extracts the vector from a skew-symmetric matrix. The rotation vector is therefore the Lie-algebra coordinate of $R$.

On-manifold optimization. Bundle adjustment, SLAM back ends and perspective-n-point refinement keep each rotation in a global form and solve for a small increment $\boldsymbol{\delta} \in \mathbb{R}^3$, applied on the right as $R\, \operatorname{Exp}(\boldsymbol{\delta})$ (the “right-plus” of Solà et al., 2018, arXiv v9, eq. 25) or on the left as $\operatorname{Exp}(\boldsymbol{\delta})\, R$. Each nonlinear least-squares step then has three unconstrained parameters per rotation and, near the current estimate, no singularity.

Properties

Distances between rotations

The natural distance between $R_1$ and $R_2$ is the geodesic angle, the angle of the relative rotation $R_1^\top R_2$:

$$d(R_1, R_2) = \lVert \operatorname{Log}(R_1^\top R_2) \rVert = \arccos\frac{\operatorname{tr}(R_1^\top R_2) - 1}{2} \in [0, \pi].$$

Zhou et al. (2019) report errors in this metric. For quaternions the same angle is $2 \arccos \lvert q_1^\top q_2 \rvert$, with the absolute value absorbing the sign. The chordal distance $\lVert R_1 - R_2 \rVert_F$ is cheaper; expanding the norm gives $2\sqrt{2} \sin(d/2)$, so it ranks pairs the same way as the geodesic angle.

Conversions and numerical pitfalls

• Drift. Products of many rotations accumulate rounding error and leave $SO(3)$. The fix is to project back to the nearest rotation in the Frobenius norm, $U \operatorname{diag}(1, 1, \det(UV^\top))\, V^\top$ from the SVD $M = U \Sigma V^\top$; the determinant term prevents a reflection. Quaternions are renormalized instead.
• Averaging. An element-wise mean of rotation matrices is not a rotation. Projecting it by SVD gives the chordal $L_2$ mean, which SciPy’s Rotation.mean computes.
• Quaternion sign. Align signs before averaging, interpolating or comparing quaternions (for example so that $q_i^\top q_0 \ge 0$).
• Angles from the trace. $\arccos$ is ill-conditioned near 0 and 180°: for tiny rotations it returns exactly zero. Using $\operatorname{atan2}(\lVert (R - R^\top)^\vee \rVert / 2,\ (\operatorname{tr} R - 1)/2)$ avoids this.

Library Conventions

As of October 2026, from each project’s official documentation or source:

Library	Quaternion order	Other conventions
SciPy Rotation 1.18.0	$(x, y, z, w)$; $(w, x, y, z)$ with scalar_first=True, added in 1.14.0	Euler sequences in upper case are intrinsic, lower case extrinsic; as_euler warns at gimbal lock and sets the third angle to zero; from_matrix orthogonalizes its input; mean is the chordal $L_2$ mean
Eigen 3.4.1 (5.0.1 development docs agree)	constructor Quaternion(w, x, y, z); stored and returned by coeffs() as $(x, y, z, w)$	rotation operations are undefined for non-unit quaternions; angularDistance, slerp
ROS 2	geometry_msgs/Quaternion fields x, y, z, w, identity $(0, 0, 0, 1)$	REP 103: right-handed frames, quaternions preferred, Euler angles discouraged; _optical camera frames use $z$ forward, $x$ right, $y$ down
OpenCV 4.14.0	cv::Quat is $[w, x, y, z]$	calibration and pose functions exchange world-to-camera rotations as rotation vectors (rvec), converted by Rodrigues

The traps are in the mixtures: ROS messages and SciPy’s default put $w$ last, and the ROS 2 quaternion tutorial warns that Eigen puts it first, yet Eigen stores it last.

Examples

The examples orient a camera toward a synthetic cube and probe each pitfall. They are deterministic and were run with NumPy 2.3.5 and SciPy 1.16.3, whose Rotation calls used here match the 1.18.0 documentation.

Round trips

A look-at construction gives the world-to-camera rotation for OpenCV’s axes ($x$ right, $y$ down, $z$ forward); its rows are the camera axes in world coordinates. Hand-written Rodrigues and quaternion formulas are checked against SciPy.

import numpy as np
from scipy.spatial.transform import Rotation

np.set_printoptions(precision=4, suppress=True)

def look_at(C, target, up=(0.0, 0.0, 1.0)):
    """World-to-camera rotation for an OpenCV-style camera (x right, y down, z forward)."""
    z = (target - C) / np.linalg.norm(target - C)
    x = np.cross(z, up); x /= np.linalg.norm(x)
    y = np.cross(z, x)
    return np.stack([x, y, z])            # rows = camera axes in world coordinates

# A synthetic object: the 8 corners of a cube centred at the world origin.
cube = np.array([[i, j, k] for i in (-1, 1) for j in (-1, 1) for k in (-1, 1)], float) * 0.5
C = np.array([4.0, -3.0, 2.0])            # camera centre
R = look_at(C, cube.mean(0))
t = -R @ C                                # X_c = R X_w + t
K = np.array([[800, 0, 320], [0, 800, 240], [0, 0, 1.0]])

print("R^T R = I:", np.allclose(R.T @ R, np.eye(3)), " det R = %.6f" % np.linalg.det(R))
x = K @ (R @ cube.mean(0) + t)
print("object centre projects to pixel", x[:2] / x[2])

def skew(v):
    return np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])

def rodrigues(r):
    """Rotation vector -> matrix: R = cos(th) I + (1 - cos(th)) u u^T + sin(th) [u]x."""
    th = np.linalg.norm(r)
    if th < 1e-12:
        return np.eye(3) + skew(r)
    u = r / th
    return np.cos(th) * np.eye(3) + (1 - np.cos(th)) * np.outer(u, u) + np.sin(th) * skew(u)

def quat_to_matrix(q):
    """Unit quaternion (w, x, y, z), Hamilton convention -> rotation matrix."""
    w, x, y, z = q / np.linalg.norm(q)
    return np.array([[1 - 2*(y*y + z*z), 2*(x*y - w*z),     2*(x*z + w*y)],
                     [2*(x*y + w*z),     1 - 2*(x*x + z*z), 2*(y*z - w*x)],
                     [2*(x*z - w*y),     2*(y*z + w*x),     1 - 2*(x*x + y*y)]])

rot = Rotation.from_matrix(R)
q = rot.as_quat(scalar_first=True)        # (w, x, y, z); SciPy's default is (x, y, z, w)
rv = rot.as_rotvec()
eul = rot.as_euler("ZYX", degrees=True)   # intrinsic yaw, pitch, roll
print("quaternion (w, x, y, z):", q)
print("rotation vector:", rv, " angle %.3f deg" % np.degrees(np.linalg.norm(rv)))
print("Euler ZYX (deg):", eul)

err = lambda A: np.abs(A - R).max()
print("round trips, max |R' - R|:")
print("  quaternion (own formula)  %.1e" % err(quat_to_matrix(q)))
print("  rotation vector (own)     %.1e" % err(rodrigues(rv)))
print("  Euler ZYX (SciPy)         %.1e" % err(Rotation.from_euler("ZYX", eul, degrees=True).as_matrix()))
print("  scalar-last read as (w,x,y,z) by mistake: %.2f" % err(quat_to_matrix(rot.as_quat())))

Output:

R^T R = I: True  det R = 1.000000
object centre projects to pixel [320. 240.]
quaternion (w, x, y, z): [ 0.5014  0.7406  0.3703 -0.2507]
rotation vector: [ 1.7901  0.895  -0.606 ]  angle 119.809 deg
Euler ZYX (deg): [ 26.344   47.9689 123.6901]
round trips, max |R' - R|:
  quaternion (own formula)  3.3e-16
  rotation vector (own)     2.2e-16
  Euler ZYX (SciPy)         2.2e-16
  scalar-last read as (w,x,y,z) by mistake: 1.73

Every round trip returns the matrix to rounding error. The last line is the classic bug: reading SciPy’s default scalar-last quaternion as $(w, x, y, z)$ yields a valid but wrong rotation, with entries off by up to 1.73 and no error raised.

Gimbal lock and the double cover

import warnings
import numpy as np
from scipy.spatial.transform import Rotation

def euler_zyx(yaw, pitch, roll):
    return Rotation.from_euler("ZYX", [yaw, pitch, roll], degrees=True).as_matrix()

def geodesic_deg(A, B):
    """Angle of the rotation A^T B, in degrees (atan2 form, see the last example)."""
    M = A.T @ B
    s = np.linalg.norm([M[2, 1] - M[1, 2], M[0, 2] - M[2, 0], M[1, 0] - M[0, 1]]) / 2
    return np.degrees(np.arctan2(s, (np.trace(M) - 1) / 2))

# Gimbal lock: at pitch = 90 deg, yaw and roll turn about the same axis.
A = euler_zyx(30, 90, 10)
B = euler_zyx(50, 90, 30)                    # yaw and roll both +20 deg
print("pitch 90: (30, 90, 10) vs (50, 90, 30), angle between them: %.1e deg" % geodesic_deg(A, B))

def jacobian_svals(angles, h=1e-6):
    """Singular values of d vec(R) / d(yaw, pitch, roll) per radian, by central differences."""
    J = np.column_stack([(euler_zyx(*(angles + h * e)) - euler_zyx(*(angles - h * e))).ravel() / (2 * h)
                         for e in np.eye(3)]) * 180 / np.pi
    return np.linalg.svd(J, compute_uv=False)

for pitch in (20.0, 89.0, 90.0):
    print(f"pitch {pitch:4.0f}: Jacobian singular values", np.array2string(jacobian_svals(np.array([30.0, pitch, 10.0])), precision=3))

with warnings.catch_warnings(record=True) as w:
    warnings.simplefilter("always")
    back = Rotation.from_matrix(A).as_euler("ZYX", degrees=True)
print("SciPy as_euler at the lock:", np.round(back, 4), "| warning:", w[0].message if w else None)

# Double cover: q and -q are the same rotation.
q = Rotation.from_euler("ZYX", [30, 40, 10], degrees=True).as_quat(scalar_first=True)
Rq = Rotation.from_quat(q, scalar_first=True).as_matrix()
Rmq = Rotation.from_quat(-q, scalar_first=True).as_matrix()
print("q =", np.round(q, 4), " |R(q) - R(-q)| max: %.1e" % np.abs(Rq - Rmq).max())
print("Euclidean distance |q - (-q)| = %.1f, rotation angle between them = %.1e deg"
      % (np.linalg.norm(2 * q), geodesic_deg(Rq, Rmq)))

Output:

pitch 90: (30, 90, 10) vs (50, 90, 30), angle between them: 2.3e-14 deg
pitch   20: Jacobian singular values [1.638 1.414 1.147]
pitch   89: Jacobian singular values [2.    1.414 0.017]
pitch   90: Jacobian singular values [2.000e+00 1.414e+00 1.034e-08]
SciPy as_euler at the lock: [20. 90.  0.] | warning: Gimbal lock detected. Setting third angle to zero since it is not possible to uniquely determine all angles.
q = [ 0.9119 -0.0091  0.3503  0.2135]  |R(q) - R(-q)| max: 0.0e+00
Euclidean distance |q - (-q)| = 2.0, rotation angle between them = 0.0e+00 deg

At 90° pitch, two angle triples with the same $\psi - \phi$ give the same orientation, and the smallest singular value of the Jacobian of $R$ with respect to the angles falls to finite-difference noise: a lost degree of freedom. Negating a quaternion moves it a Euclidean distance of 2 but leaves the rotation unchanged.

Interpolation and continuity

Naive Euler-angle interpolation and SciPy’s Slerp are compared at 101 points with the geodesic $R_0 \operatorname{Exp}(t \operatorname{Log}(R_0^\top R_1))$, computed with the matrix exponential and logarithm. Then two rotations 0.2° apart, either side of 180°, are compared in each representation.

import numpy as np
from scipy.linalg import expm, logm
from scipy.spatial.transform import Rotation, Slerp

def geodesic_deg(A, B):
    """Angle of A^T B in degrees, via atan2 (accurate near 0 and 180 degrees)."""
    M = A.T @ B
    s = np.linalg.norm([M[2, 1] - M[1, 2], M[0, 2] - M[2, 0], M[1, 0] - M[0, 1]]) / 2
    return np.degrees(np.arctan2(s, (np.trace(M) - 1) / 2))

def path_length(Rs):
    return sum(geodesic_deg(a, b) for a, b in zip(Rs[:-1], Rs[1:]))

ts = np.linspace(0, 1, 101)
# Pairs of camera orientations as intrinsic ZYX (yaw, pitch, roll) in degrees.
for e0, e1 in [((10, 20, 0), (80, 70, 40)),          # no angle wraps around
               ((170, 50, 0), (-170, 60, 30))]:     # yaw crosses +-180
    e0, e1 = np.array(e0, float), np.array(e1, float)
    R0 = Rotation.from_euler("ZYX", e0, degrees=True).as_matrix()
    R1 = Rotation.from_euler("ZYX", e1, degrees=True).as_matrix()
    W = np.real(logm(R0.T @ R1))                      # reference geodesic R0 exp(t log(R0^T R1))
    geo = [R0 @ expm(t * W) for t in ts]
    slerp = Slerp([0, 1], Rotation.from_matrix([R0, R1]))(ts).as_matrix()
    euler = Rotation.from_euler("ZYX", (1 - ts)[:, None] * e0 + ts[:, None] * e1,
                                degrees=True).as_matrix()
    err_e = max(geodesic_deg(a, g) for a, g in zip(euler, geo))
    err_s = max(geodesic_deg(a, g) for a, g in zip(slerp, geo))
    print(f"{e0} -> {e1}: endpoints {geodesic_deg(R0, R1):.2f} deg apart")
    print(f"  max error from geodesic: Euler-lerp {err_e:.2f} deg, slerp {err_s:.1e} deg")
    print(f"  path length: Euler-lerp {path_length(euler):.2f} deg, slerp {path_length(slerp):.2f} deg")

# Continuity: two rotations 0.2 deg apart, either side of 180 deg about the same axis.
axis = np.array([1.0, 2.0, 2.0]) / 3
Ra = Rotation.from_rotvec(np.radians(179.9) * axis)
Rb = Rotation.from_rotvec(np.radians(180.1) * axis)
print("angle between Ra and Rb: %.2f deg" % geodesic_deg(Ra.as_matrix(), Rb.as_matrix()))
d = lambda u, v: np.linalg.norm(u - v)
print("distance between representations:")
print("  rotation vector       %.3f" % d(Ra.as_rotvec(), Rb.as_rotvec()))
print("  quaternion (w >= 0)   %.3f" % d(Ra.as_quat(canonical=True), Rb.as_quat(canonical=True)))
six = lambda r: r.as_matrix()[:, :2].T.ravel()      # first two columns
print("  6D (two columns)      %.3f" % d(six(Ra), six(Rb)))

def from_6d(a):
    """Zhou et al.'s Gram-Schmidt map from 6 numbers to a rotation matrix."""
    b1 = a[:3] / np.linalg.norm(a[:3])
    b2 = a[3:] - (b1 @ a[3:]) * b1
    b2 /= np.linalg.norm(b2)
    return np.column_stack([b1, b2, np.cross(b1, b2)])

noisy = six(Ra) + np.array([0.05, -0.02, 0.03, 0.01, 0.04, -0.03])   # an unconstrained "network output"
Rn = from_6d(noisy)
print("6D round trip error %.1e; noisy 6D -> R: orthogonality %.1e, det %.6f"
      % (np.abs(from_6d(six(Ra)) - Ra.as_matrix()).max(), np.abs(Rn.T @ Rn - np.eye(3)).max(), np.linalg.det(Rn)))

Output:

[10. 20.  0.] -> [80. 70. 40.]: endpoints 66.49 deg apart
  max error from geodesic: Euler-lerp 11.95 deg, slerp 2.5e-14 deg
  path length: Euler-lerp 71.54 deg, slerp 66.49 deg
[170.  50.   0.] -> [-170.   60.   30.]: endpoints 20.24 deg apart
  max error from geodesic: Euler-lerp 179.63 deg, slerp 2.8e-14 deg
  path length: Euler-lerp 365.09 deg, slerp 20.24 deg
angle between Ra and Rb: 0.20 deg
distance between representations:
  rotation vector       6.280
  quaternion (w >= 0)   2.000
  6D (two columns)      0.004
6D round trip error 2.2e-16; noisy 6D -> R: orthogonality 3.0e-17, det 1.000000

Slerp follows the geodesic to rounding error. Euler-angle interpolation strays by 12° even without wrap-around, and when yaw crosses ±180° it turns the camera through 365° to cover a 20° change. Near 180°, the rotation vectors are nearly $2\pi$ apart and the canonical quaternions 2 apart, while the 6D vectors differ by 0.004: the discontinuity Zhou et al. (2019) identify. The Gram–Schmidt map turns a perturbed 6D vector into an exact rotation.

Averaging and re-orthonormalization

import numpy as np
from scipy.spatial.transform import Rotation

def geodesic_deg(A, B):
    M = A.T @ B
    s = np.linalg.norm([M[2, 1] - M[1, 2], M[0, 2] - M[2, 0], M[1, 0] - M[0, 1]]) / 2
    return np.degrees(np.arctan2(s, (np.trace(M) - 1) / 2))

def project_to_so3(M):
    """Nearest rotation in the Frobenius norm: U diag(1, 1, det(U V^T)) V^T."""
    U, _, Vt = np.linalg.svd(M)
    D = np.diag([1.0, 1.0, np.sign(np.linalg.det(U @ Vt))])
    return U @ D @ Vt

# Five estimates of one camera's orientation, each off by 10 to 30 deg about some axis.
R_true = Rotation.from_euler("ZYX", [40, -30, 15], degrees=True)
axes = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 1, 0], [0, -1, 1]], float)
angles = np.radians([30, -20, 25, 10, -15])
offsets = Rotation.from_rotvec(angles[:, None] * axes / np.linalg.norm(axes, axis=1)[:, None])
Rs = (R_true * offsets).as_matrix()
print("estimates' errors (deg):", np.round([geodesic_deg(R, R_true.as_matrix()) for R in Rs], 1))

M = Rs.mean(axis=0)                                   # element-wise average
print("element-wise mean: max |M^T M - I| = %.3f, det = %.3f, singular values %s"
      % (np.abs(M.T @ M - np.eye(3)).max(), np.linalg.det(M), np.round(np.linalg.svd(M, compute_uv=False), 3)))
P = project_to_so3(M)
print("after SVD projection: max |P^T P - I| = %.1e, det = %.6f, error %.2f deg"
      % (np.abs(P.T @ P - np.eye(3)).max(), np.linalg.det(P), geodesic_deg(P, R_true.as_matrix())))
print("SciPy Rotation.mean agrees: %.1e deg" % geodesic_deg(P, Rotation.from_matrix(Rs).mean().as_matrix()))

# Quaternion sign: averaging q's without aligning their signs.
q = Rotation.from_matrix(Rs).as_quat(scalar_first=True)
q[::2] *= -1                                          # same rotations, three of five sign-flipped
naive = q.mean(axis=0)
aligned = (q * np.sign(q @ q[0])[:, None]).mean(axis=0)
for name, m in [("naive", naive), ("sign-aligned", aligned)]:
    R_m = Rotation.from_quat(m, scalar_first=True).as_matrix()
    print(f"quaternion mean, {name:12s}: |mean| = {np.linalg.norm(m):.3f}, "
          f"error {geodesic_deg(R_m, R_true.as_matrix()):.2f} deg")

# Drift: compose 100,000 small camera increments in float32.
step = Rotation.from_rotvec([1e-3, 2e-3, -1e-3]).as_matrix().astype(np.float32)
R = np.eye(3, dtype=np.float32)
for k in range(1, 100_001):
    R = R @ step
    if k in (1_000, 10_000, 100_000):
        print(f"after {k:6d} products: max |R^T R - I| = {np.abs(R.T @ R - np.eye(3)).max():.1e}")
R = project_to_so3(R.astype(np.float64))
print("re-orthonormalized:            max |R^T R - I| = %.1e" % np.abs(R.T @ R - np.eye(3)).max())

# Angle from the trace: arccos loses small angles; atan2 does not.
small = Rotation.from_rotvec([1e-8, 0, 0]).as_matrix()
print("1e-8 rad rotation: arccos form %.2e rad, atan2 form %.2e rad"
      % (np.arccos(np.clip((np.trace(small) - 1) / 2, -1, 1)), np.radians(geodesic_deg(np.eye(3), small))))

Output:

estimates' errors (deg): [30. 20. 25. 10. 15.]
element-wise mean: max |M^T M - I| = 0.082, det = 0.888, singular values [0.967 0.959 0.957]
after SVD projection: max |P^T P - I| = 3.9e-16, det = 1.000000, error 7.99 deg
SciPy Rotation.mean agrees: 1.5e-13 deg
quaternion mean, naive       : |mean| = 0.200, error 36.04 deg
quaternion mean, sign-aligned: |mean| = 0.985, error 8.01 deg
after   1000 products: max |R^T R - I| = 1.9e-05
after  10000 products: max |R^T R - I| = 1.6e-04
after 100000 products: max |R^T R - I| = 1.7e-03
re-orthonormalized:            max |R^T R - I| = 8.9e-16
1e-8 rad rotation: arccos form 0.00e+00 rad, atan2 form 1.00e-08 rad

The element-wise mean of five estimates is not orthogonal (determinant 0.888); its SVD projection is a rotation and matches SciPy’s Rotation.mean (its 8° error comes from offsets that do not cancel). Averaging quaternions with mixed signs gives a 36° error; aligning the signs first gives 8.01°, close to the chordal mean. After 100,000 single-precision products, $R^\top R$ is off the identity by $1.7 \times 10^{-3}$ until one SVD projection restores it. For a $10^{-8}$ rad rotation, the $\arccos$ form returns zero.

Common Misconceptions

• “A quaternion is four numbers, so $(q_0, q_1, q_2, q_3)$ is unambiguous.” The scalar may come first or last, and Hamilton and JPL quaternions differ in handedness.
• “Euler angles are fine away from 90° pitch.” Their wrap-around at ±180° still breaks interpolation, distances and losses.
• “$q$ and $-q$ are different orientations.” They are the same rotation; distances and losses must take the absolute value of $q_1^\top q_2$.
• “A rotation vector is a vector, so it can be added.” Adding rotation vectors is not composition except for rotations about the same axis; composition goes through the matrices or quaternions.

Where It Is Used

Every estimate of a camera’s orientation needs a representation: perspective-n-point, the essential matrix decomposition, iterative closest point (whose closed-form solutions use either the SVD or unit quaternions), bundle adjustment, visual odometry, SLAM and the error-state extended Kalman filter, which estimates orientation errors as rotation vectors. Networks that regress object or camera pose can output the 6D form of Zhou et al. (2019). Camera pose combines a rotation with a translation and a transform direction.