From quaternion to matrix
Write the unit quaternion as q = (w,\ x,\ y,\ z) with
w^2 + x^2 + y^2 + z^2 = 1. We want the matrix
R with R\mathbf{v} = q\,(0,\mathbf{v})\,q^{-1} for
every \mathbf{v}.
Step 1 — the columns are the images of the basis vectors. Any matrix is determined by
what it does to \mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3: column
j of R is
R\mathbf{e}_j = q\,(0,\mathbf{e}_j)\,q^{-1}. So we just sandwich each basis
vector and read off the result.
Step 2 — sandwich one basis vector. Take
\mathbf{e}_1 = (1,0,0), i.e. the pure quaternion p = i,
and compute q\,i\,q^{*} using q^{-1} = q^{*} = (w,-x,-y,-z).
Multiplying out with the rules i^2 = j^2 = k^2 = ijk = -1 and collecting the
i, j, k coefficients gives the first column of R:
q\,i\,q^{*} = \big(1 - 2(y^2 + z^2)\big)\,i + 2(xy + wz)\,j + 2(xz - wy)\,k.
Reading off the i, j, k coefficients gives the column
\big(1 - 2(y^2+z^2),\ 2(xy+wz),\ 2(xz-wy)\big)^{\!\top}. Notice the diagonal
entry: 1 - 2(y^2 + z^2) says "how much of \mathbf{e}_1
survives" — a rotation about an axis tilts \mathbf{e}_1 away by an amount set by
the axis components y and z perpendicular to it.
Step 3 — repeat for the other two basis vectors (sandwich
j and k the same way) and stack the three columns.
The whole matrix is
R = \begin{pmatrix} 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{pmatrix}.
Every entry is a quadratic in (w,x,y,z) — no trig, just multiplies and adds,
which is exactly why it is GPU-friendly. The diagonal carries the
1 - 2(\cdot) "what survives" terms; the off-diagonals come in
2(\,\cdot \pm w\,\cdot) pairs whose antisymmetric part
(\pm 2w\cdot) is the genuine turn and whose symmetric part
(2\,\cdot) is the axis tilt. (Sanity check: set
q = (1,0,0,0) — the no-rotation quaternion — and every off-diagonal vanishes,
every diagonal is 1, so R = I, the identity.)
From matrix back to quaternion
Step 1 — read w off the trace. Add the three diagonal
entries of R above:
\operatorname{tr} R = 3 - 4(x^2 + y^2 + z^2) = 3 - 4(1 - w^2) = 4w^2 - 1,
using x^2 + y^2 + z^2 = 1 - w^2. Solving,
w = \tfrac{1}{2}\sqrt{1 + \operatorname{tr} R}\,.
Step 2 — recover x, y, z from the off-diagonal pairs.
Subtract the symmetric partners to kill the tilt and leave only the w term:
R_{32} - R_{23} = 4wx, \qquad R_{13} - R_{31} = 4wy, \qquad R_{21} - R_{12} = 4wz,
so, once w is known,
x = \frac{R_{32} - R_{23}}{4w}, \qquad y = \frac{R_{13} - R_{31}}{4w}, \qquad z = \frac{R_{21} - R_{12}}{4w}.
Step 3 — mind the degenerate case. If \operatorname{tr} R \approx -1
then w \approx 0 and dividing by 4w is unstable;
real code instead reads the largest diagonal entry first to recover the biggest of
x, y, z, then divides by that. The idea is identical — solve from the
trace and the off-diagonals — only the pivot changes for numerical safety.
Let q = (w, x, y, z) be a unit quaternion.
-
Its rotation matrix is
R = \begin{pmatrix} 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{pmatrix},
with R\mathbf{v} = q\,(0,\mathbf{v})\,q^{-1} for every
\mathbf{v}.
-
The reverse uses the trace:
w = \tfrac12\sqrt{1 + \operatorname{tr} R}, then
x = (R_{32}-R_{23})/4w,
y = (R_{13}-R_{31})/4w,
z = (R_{21}-R_{12})/4w (pivot on the largest diagonal entry when
w \approx 0).
-
Store quaternions, render with matrices: engines keep orientation as a compact,
drift-resistant quaternion and convert to R when handing vertices to the GPU.
Both representations wander under repeated floating-point arithmetic, but they wander differently — and
that asymmetry is the whole reason for the store-as-quaternion convention.
A rotation matrix must stay orthogonal: its columns must remain unit-length and mutually
perpendicular (R^{\top}R = I). Multiply many matrices and rounding error nudges
the columns off this constraint — the matrix quietly starts to shear and scale, distorting your
model. Re-orthogonalising (Gram–Schmidt, or a polar decomposition) is comparatively expensive and easy to
get subtly wrong.
A quaternion has just one constraint to honour: |q| = 1. After
any drift, snapping it back to the unit sphere is a single normalisation
q \mapsto q/|q| — one square root and four divides — and it can never represent
a shear or a scale, only a rotation. So the cheap, safe place to accumulate rotation is the
quaternion; the matrix is generated fresh from it each frame, never left to compound its own error.