Quaternions for rotation

Open a trained 3D Gaussian splatting scene (Kerbl et al. 2023) and look at one Gaussian: a position, three scales, an opacity, spherical-harmonics colour — and four numbers rot_0..rot_3. That's a quaternion (w, x, y, z). The optimiser updates those four numbers freely by gradient descent; the renderer normalises them and turns them into a rotation matrix R, and the Gaussian's covariance is Σ = R S Sᵀ Rᵀ. Why four numbers instead of three angles or nine matrix entries? Because quaternions are the parametrisation of 3D rotation that has no singularities (unlike Euler angles), no constraints to maintain beyond "length 1" (unlike matrices, which must stay orthonormal), and cheap, smooth interpolation. This lesson builds them from the product up and ends with the exact formula 3DGS uses.

1. What a quaternion is

A quaternion is q = w + xi + yj + zk, with Hamilton's rules i² = j² = k² = ijk = −1. Those rules give ij = k, jk = i, ki = j and, reversing the order, ji = −k, kj = −i, ik = −j — the same cyclic pattern, and the same sign flip, as the cross product of the basis vectors. Write it as a scalar plus a vector, q = (w, v), and the product is

(w₁, v₁)(w₂, v₂) = (w₁w₂ − v₁·v₂, w₁v₂ + w₂v₁ + v₁ × v₂)

The v₁ × v₂ term is why the product doesn't commute. Other operations:

• conjugate q* = (w, −v); q q* = |q|², so q⁻¹ = q* / |q|² — just q* for a unit quaternion;
• norm |q| = √(w² + x² + y² + z²); rotations are the unit quaternions, the 3-sphere in 4D.

#include <print>

struct quat { float w = 1, x = 0, y = 0, z = 0; };

constexpr quat operator*(quat a, quat b) {
    return {a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z,
            a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y,
            a.w * b.y - a.x * b.z + a.y * b.w + a.z * b.x,
            a.w * b.z + a.x * b.y - a.y * b.x + a.z * b.w};
}

constexpr quat i{0, 1, 0, 0}, j{0, 0, 1, 0}, k{0, 0, 0, 1};
static_assert((i * j).z == 1 && (j * i).z == -1);   // ij = k, ji = -k
static_assert((i * i).w == -1 && (i * j * k).w == -1);

int main() {
    quat a{1, 2, 3, 4}, b{5, 6, 7, 8};
    quat ab = a * b, ba = b * a;
    std::println("ab = ({}, {}, {}, {})", ab.w, ab.x, ab.y, ab.z);
    std::println("ba = ({}, {}, {}, {})", ba.w, ba.x, ba.y, ba.z);
}

2. Rotating a vector: q v q* and the half angle

The unit quaternion for a rotation by θ about the unit axis k is

q = (cos(θ/2), sin(θ/2)·k)

and a vector v is rotated by the sandwich v' = q (0, v) q*. The result is again a pure quaternion (w = 0), and its length is unchanged. The half angle is there because q appears twice: each side turns the component of v perpendicular to k by θ/2. Forget the ½ and every rotation goes twice as far — the single most common quaternion bug.

Expanding the sandwich gives a cheaper formula with no quaternion products (u is the vector part):

t = 2 (u × v), v' = v + w·t + u × t

#include <cmath>
#include <numbers>
#include <print>

struct vec3 { float x, y, z; };
struct quat { float w, x, y, z; };

quat operator*(quat a, quat b) {
    return {a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z,
            a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y,
            a.w * b.y - a.x * b.z + a.y * b.w + a.z * b.x,
            a.w * b.z + a.x * b.y - a.y * b.x + a.z * b.w};
}
vec3 cross(vec3 a, vec3 b) { return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x}; }

quat from_axis_angle(vec3 k, float theta) {   // k unit
    float s = std::sin(theta / 2);
    return {std::cos(theta / 2), k.x * s, k.y * s, k.z * s};
}

vec3 rotate_sandwich(quat q, vec3 v) {
    quat r = q * quat{0, v.x, v.y, v.z} * quat{q.w, -q.x, -q.y, -q.z};
    return {r.x, r.y, r.z};
}

vec3 rotate_fast(quat q, vec3 v) {
    vec3 u{q.x, q.y, q.z}, t = cross(u, v);
    t = {2 * t.x, 2 * t.y, 2 * t.z};
    vec3 c = cross(u, t);
    return {v.x + q.w * t.x + c.x, v.y + q.w * t.y + c.y, v.z + q.w * t.z + c.z};
}

int main() {
    quat q = from_axis_angle({0, 0, 1}, std::numbers::pi_v<float> / 2);   // 90° about z
    vec3 a = rotate_sandwich(q, {1, 0, 0}), b = rotate_fast(q, {1, 0, 0});
    std::println("q = ({:.4f}, 0, 0, {:.4f})", q.w, q.z);
    std::println("sandwich: ({:.3f}, {:.3f}, {:.3f})  fast: ({:.3f}, {:.3f}, {:.3f})",
                 a.x, a.y, a.z, b.x, b.y, b.z);
}

Composition works like matrices: rotate(a * b, v) = rotate(a, rotate(b, v)) — the right factor acts first, because (ab) v (ab)* = a (b v b*) a*. One product composes two rotations in 16 multiplies (a 3x3 product needs 27).

3. The double cover

q and −q give the same rotation: negate q and both copies in q v q* change sign, which cancels. A full 360° turn has half angle 180°, so its quaternion is (−1, 0, 0, 0) — not the identity quaternion, but the identity rotation. Every rotation is hit exactly twice, by two antipodal points of the 4D sphere.

This matters whenever you compare or interpolate quaternions: two networks (or two keyframes) can store "the same" orientation with opposite signs, and the straight path between q and −q passes through zero.

4. From quaternion to matrix — the 3DGS formula

Expanding q (0, v) q* for the three basis vectors gives the columns of the rotation matrix. For a unit q = (w, x, y, z):

This is exactly build_rotation() in the 3DGS reference code (utils/general_utils.py, with r = w), which first divides the raw parameters by their norm. The CUDA rasterizer's computeCov3D() writes the same nine expressions into a glm::mat3 constructor — which takes them column by column, so glm's R is the transpose of this one; the code then forms M = S·R and Σ = Mᵀ M = Rᵀ S² R (in glm's transposed convention) — the same covariance. Reading research code means tracking these conventions line by line.

#include <cmath>
#include <print>

struct quat { float w, x, y, z; };
using mat3 = float[3][3];

void to_mat3(quat q, mat3& R) {
    float n = std::sqrt(q.w * q.w + q.x * q.x + q.y * q.y + q.z * q.z);
    float r = q.w / n, x = q.x / n, y = q.y / n, z = q.z / n;
    float m[3][3] = {{1 - 2 * (y * y + z * z), 2 * (x * y - r * z), 2 * (x * z + r * y)},
                     {2 * (x * y + r * z), 1 - 2 * (x * x + z * z), 2 * (y * z - r * x)},
                     {2 * (x * z - r * y), 2 * (y * z + r * x), 1 - 2 * (x * x + y * y)}};
    for (int i = 0; i < 3; ++i)
        for (int j = 0; j < 3; ++j) R[i][j] = m[i][j];
}

int main() {
    quat raw{1, 0, 0, 1};             // unnormalised: 90° about z, length √2
    float s[3] = {2.0f, 0.5f, 0.1f};  // a flat, elongated Gaussian
    mat3 R;
    to_mat3(raw, R);
    float cov[3][3] = {};             // Σ = R S Sᵀ Rᵀ
    for (int i = 0; i < 3; ++i)
        for (int j = 0; j < 3; ++j)
            for (int k = 0; k < 3; ++k) cov[i][j] += R[i][k] * s[k] * s[k] * R[j][k];
    for (auto& row : cov) std::println("{:7.3f} {:7.3f} {:7.3f}", row[0], row[1], row[2]);
    // The long axis (scale 2 along local x) now lies along world y: Σ_yy = 4.
}

Two properties make R S Sᵀ Rᵀ a good parametrisation: it's a valid covariance (symmetric positive semi-definite) for any parameter values, and its eigenvectors are the columns of R with eigenvalues sᵢ² — the Gaussian's principal axes. Optimising Σ's six entries directly would need a constraint to stay positive definite.

5. Interpolation: slerp and nlerp

Unit quaternions live on a sphere, so "interpolate two rotations" means "walk along the great circle between them". slerp (Shoemake 1985) does it at constant angular velocity: with cos θ = a·b,

slerp(a, b, t) = (sin((1 − t)θ)·a + sin(tθ)·b) / sin θ

Two guards make it robust:

1. Shortest arc. If a·b < 0, negate b (and the dot). Otherwise slerp takes the long way round — more than 180° of rotation.
2. Nearly parallel. As a → b, θ → 0 and sin θ → 0: the weights become 0/0 = NaN, and a dot product that rounds to 1.0000001 makes acos NaN. For a·b > 0.9995, use nlerp — normalize((1 − t)·a + t·b) — which agrees with slerp to float precision over a tiny arc. Clamp the dot to ≤ 1 before acos regardless.

nlerp follows the same arc but not at constant speed: it's slower near the ends and faster in the middle. For small arcs (animation blending, per-frame smoothing) that's invisible and nlerp is cheaper, commutative, and easy to weight across many inputs.

#include <algorithm>
#include <cmath>
#include <numbers>
#include <print>

struct quat { float w, x, y, z; };
quat operator*(quat a, float s) { return {a.w * s, a.x * s, a.y * s, a.z * s}; }
quat operator+(quat a, quat b) { return {a.w + b.w, a.x + b.x, a.y + b.y, a.z + b.z}; }
float dot(quat a, quat b) { return a.w * b.w + a.x * b.x + a.y * b.y + a.z * b.z; }
quat normalized(quat q) { return q * (1 / std::sqrt(dot(q, q))); }

quat slerp(quat a, quat b, float t) {
    float d = dot(a, b);
    if (d < 0) { b = b * -1; d = -d; }                       // shortest arc
    if (d > 0.9995f) return normalized(a * (1 - t) + b * t); // nearly parallel
    float th = std::acos(std::min(d, 1.0f)), s = std::sin(th);
    return a * (std::sin((1 - t) * th) / s) + b * (std::sin(t * th) / s);
}

float degrees_between(quat a, quat b) {
    return 2 * std::acos(std::min(std::fabs(dot(a, b)), 1.0f)) * 180 / std::numbers::pi_v<float>;
}

int main() {
    float h = 160.0f / 2 * std::numbers::pi_v<float> / 180;
    quat a{1, 0, 0, 0}, b{std::cos(h), 0, 0, std::sin(h)};   // identity → 160° about z
    for (float t : {0.25f, 0.5f, 0.75f})
        std::println("t = {}: slerp {:.1f}°  nlerp {:.1f}°", t,
                     degrees_between(a, slerp(a, b, t)), degrees_between(a, normalized(a * (1 - t) + b * t)));
    quat nb = b * -1;   // same rotation, opposite sign
    std::println("slerp towards -b at 0.5: {:.1f}°", degrees_between(a, slerp(a, nb, 0.5f)));
    std::println("slerp(a, a, 0.5) is finite: {}", std::isfinite(slerp(a, a, 0.5f).w));
}

6. Why not Euler angles?

Yaw, pitch, roll are three rotations about fixed axes in sequence, R = R_z(yaw)·R_y(pitch)·R_x(roll). They're intuitive for a UI, but when the middle angle reaches ±90° the first and last axes line up and only their difference matters: a degree of freedom disappears. That's gimbal lock — a property of the parametrisation, not of rounding — and near it, small orientation changes need huge angle changes, which is poison for an optimiser or an animation curve.

#include <array>
#include <cmath>
#include <numbers>
#include <print>

using M = std::array<std::array<double, 3>, 3>;
M mul(const M& a, const M& b) {
    M r{};
    for (int i = 0; i < 3; ++i)
        for (int j = 0; j < 3; ++j)
            for (int k = 0; k < 3; ++k) r[i][j] += a[i][k] * b[k][j];
    return r;
}
M Rx(double a) { return {{{1, 0, 0}, {0, std::cos(a), -std::sin(a)}, {0, std::sin(a), std::cos(a)}}}; }
M Ry(double a) { return {{{std::cos(a), 0, std::sin(a)}, {0, 1, 0}, {-std::sin(a), 0, std::cos(a)}}}; }
M Rz(double a) { return {{{std::cos(a), -std::sin(a), 0}, {std::sin(a), std::cos(a), 0}, {0, 0, 1}}}; }

double snap(double v) { return std::fabs(v) < 5e-4 ? 0.0 : v; }   // print -0.000 as 0

int main() {
    double d = std::numbers::pi / 180;
    for (auto [yaw, roll] : {std::pair{30.0, 0.0}, {0.0, -30.0}, {20.0, -10.0}}) {
        M r = mul(Rz(yaw * d), mul(Ry(90 * d), Rx(roll * d)));
        std::println("yaw {:>3} roll {:>3}: row0 = ({:.3f}, {:.3f}, {:.3f})", yaw, roll,
                     snap(r[0][0]), snap(r[0][1]), snap(r[0][2]));
    }
    // All three are the same rotation: at pitch 90° only yaw − roll matters.
}

Matrices avoid gimbal lock but have nine numbers for three degrees of freedom: an optimiser or an integrator drifts off the rotation group, and pulling it back (Gram-Schmidt, SVD) is expensive. Quaternions have four numbers and one cheap constraint — |q| = 1, restored with a single normalisation. Repeated q = dq * q updates do drift in length through rounding; renormalise every frame, or (like 3DGS) normalise inside every use.

Storage order varies: 3DGS's PLY stores (w, x, y, z); glm's constructor takes (w, x, y, z) but stores x, y, z, w; Eigen likewise (coeffs() is x-first); ROS, Unity and many file formats use (x, y, z, w). Swapping w and z yields a valid-looking but wrong rotation, so check one known value (90° about z must give w = z ≈ 0.707) whenever data crosses a boundary.

In an interview

Expect "why do quaternions use the half angle?" (q appears twice in the sandwich), "how do you interpolate rotations?" (slerp with the sign flip and the nlerp fallback — say both guards), and "what's gimbal lock?". For a 3DGS-flavoured role: "how is each Gaussian's covariance parametrised, and why not optimise Σ directly?" — R(q) S Sᵀ R(q)ᵀ stays positive semi-definite for any parameters, with q normalised in the forward pass.

In Python (NumPy / PyTorch)

The 3DGS code builds all rotation matrices at once: q = r / r.norm(dim=1, keepdim=True) then nine vectorised expressions into an (N, 3, 3) tensor, and L = R @ S; Σ = L @ L.transpose(1, 2). Autograd differentiates through the normalisation, so the gradient on the raw quaternion has no component along q itself — scaling q changes nothing. SciPy's Rotation.from_quat expects (x, y, z, w) by default; mind the order.