Rotations with quaternions

A quaternion is a 4 dimensional complex-like number, it has four components, three of which are the "imaginary" part.
$$ q = w+x\textrm{i}+y\textrm{j}+z\textrm{k} $$ $$ q = (x,y,z,w) $$ $$ \textrm{i}^{2}=\textrm{j}^{2}=\textrm{k}^{2}=\textrm{i}\textrm{j}\textrm{k}=-1 $$

We represent a quaternion with this data structure:

typedef union{
    float q[4];
    struct{
        float x;
        float y;
        float z;
        float w;
    };
} Quaternion;

The four components are usually ordered \(w,x,y,z\) but I like to put \(w\) at the end.

Initializing a quaternion:

Quaternion q = (Quaternion){1, 2, 3, 4};

Quaternion magnitude

A quaternion is basically a 4 dimensional vector, so it has a magnitude (or norm, or length):

$$||q|| = \sqrt{x^{2}+y^{2}+z^{2}+w^{2}}$$
float quat_magnitude(Quaternion q){
    return sqrt(q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w);
}

Quaternion normalization

A quaternion can be normalized by dividing each component by the magnitude:

Quaternion quat_normalize(Quaternion q){
    float m = quat_magnitude(q);
    if(m == 0) return (Quaternion){0, 0, 0, 0};
    return (Quaternion){
        q.x/m,
        q.y/m,
        q.z/m,
        q.w/m
    };
}

A special property of quaternions is that a unit quaternion (a quaternion with magnitude \(1\)) represents a rotation in 3D space.

Identity quaternion

There is a special quaternion called the identity quaternion which corresponds to no rotation:

Quaternion quat_id(){
    return (Quaternion){0, 0, 0, 1};
}

Geometrically, we can also consider \((0, 0, 0, -1)\) to be an identity quaternion since it corresponds to no rotation.

Scaling a quaternion

Scaling a quaternion is multiplying each of its components by a real number (the scalar):

Quaternion quat_scale(Quaternion q, float s){
    return (Quaternion){q.x*s, q.y*s, q.z*s, q.w*s};
}

Quaternion multiplication

Multiplying two unit quaternions represents a composition of two rotations.

Quaternion multiplication isn't commutative (\(q_{1}.q_{2} \ne q_{2}.q_{1}\)). If we want to apply a rotation \(q_{1}\) then a rotation \(q_{2}\), the resulting rotation \(q_{3}\) is:

$$q_{3}=q_{2}.q_{1}$$

Quaternion multiplication looks like this:

$$q_{1} = a+b\textrm{i}+c\textrm{j}+d\textrm{k}$$ $$q_{2} = e+f\textrm{i}+g\textrm{j}+h\textrm{k}$$ $$\begin{align*} q_{1}.q_{2} = (ae-bf-cg-dh)+(af+be+ch-dg)\textrm{i}+\\ (ag-bh+ce+df)\textrm{j}+(ah+bg-cf+de)\textrm{k}\end{align*}$$
Quaternion quat_mul(Quaternion a, Quaternion b){
    return (Quaternion){
        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,
        a.w*b.w - a.x*b.x - a.y*b.y - a.z*b.z
    };
}

Quaternion vs Euler angles

We use quaternions instead of Euler angles to represent rotations for a couple of reasons:

We represent the orientation of an object using only a quaternion, then we multiply that orientation by another quaternion to rotate it.

However writing a rotation directly in quaternion form isn't really intuitive, what we do instead is convert an Euler angle to a quaternion then use it for rotating.

If we have an Euler angle rotation in the order ZYX (Yaw -> Pitch -> Roll, we can chose any order but must stay consistent), we can convert it to a quaternion like this:

$$ q = \begin{bmatrix} \sin(x/2)\cos(y/2)\cos(z/2)-\cos(x/2)\sin(y/2)\sin(z/2) \\ \cos(x/2)\sin(y/2)\cos(z/2)+\sin(x/2)\cos(y/2)\sin(z/2) \\ \cos(x/2)\cos(y/2)\sin(z/2)-\sin(x/2)\sin(y/2)\cos(z/2) \\ \cos(x/2)\cos(y/2)\cos(z/2)+\sin(x/2)\sin(y/2)\sin(z/2) \end{bmatrix} $$
typedef union{
    float v[3];
    struct{
        float x;
        float y;
        float z;
    };
} Vector3;

Quaternion euler_to_quat(Vector3 e){
    float cx = cos(e.x/2);
    float sx = sin(e.x/2);
    float cy = cos(e.y/2);
    float sy = sin(e.y/2);
    float cz = cos(e.z/2);
    float sz = sin(e.z/2);
    return (Quaternion){
        sx*cy*cz - cx*sy*sz,
        cx*sy*cz + sx*cy*sz,
        cx*cy*sz - sx*sy*cz,
        cx*cy*cz + sx*sy*sz
    };
}
typedef struct Transform{
    Vector3 position;
    Quaternion rotation;
    Vector3 scale;
} Transform;
Transform obj;
obj.position = (Vector3){0, 0, 0};
obj.scale = (Vector3){1, 1, 1};
obj.rotation = quat_id(); // Initially our object isn't rotated

// We rotate the object by PI/4 around the Y axis
obj.rotation = quat_mul(euler_to_quat((Vector3){0, PI/4, 0}), obj.rotation);

// We rotate again by PI/4 making it a PI/2 rotation around Y
obj.rotation = quat_mul(euler_to_quat((Vector3){0, PI/4, 0}), obj.rotation);

Quaternion to rotation matrix

When doing 3D rendering, we pass an MVP (Model View Projection) matrix to a shader to properly display our objects in the scene:

$$\textit{MVP} = M_{\textit{projection}}.M_{\textit{view}}.M_{\textit{model}}$$

The model matrix itself looks like this:

$$M_{\textit{model}} = M_{\textit{scale}}.M_{\textit{rotate}}.M_{\textit{translate}}$$

Each of those matrices is a 4x4 matrix in homogeneous coordinates.

We convert a quaternion to a rotation matrix like this:

$$q = (x, y, z, w)$$
$$ M_{\textit{rotate}} = \begin{bmatrix} 1-2yy-2zz && 2xy-2zw && 2xz+2yw && 0 \\ 2xy+2zw && 1-2xx-2zz && 2yz-2xw && 0 \\ 2xz-2yw && 2yz+2xw && 1-2xx-2yy && 0 \\ 0 && 0 && 0 && 1 \end{bmatrix} $$

Graphics APIs (like OpenGL) usually represent matrices in memory in a column-major notation, so we have to transpose the matrices in our code:

typedef union{
    float m[16];
    struct{
        float m00; float m10; float m20; float m30;
        float m01; float m11; float m21; float m31;
        float m02; float m12; float m22; float m32;
        float m03; float m13; float m23; float m33;
    };
} Mat4;

Mat4 rotate_3d_matrix(Quaternion q){
    float xx = q.x*q.x;
    float yy = q.y*q.y;
    float zz = q.z*q.z;
    return (Mat4){
        1-2*yy-2*zz, 2*q.x*q.y+2*q.z*q.w, 2*q.x*q.z-2*q.y*q.w, 0,
        2*q.x*q.y-2*q.z*q.w, 1-2*xx-2*zz, 2*q.y*q.z+2*q.x*q.w, 0,
        2*q.x*q.z+2*q.y*q.w, 2*q.y*q.z-2*q.x*q.w, 1-2*xx-2*yy, 0,
        0, 0, 0, 1
    };
}

Quaternion Conjugate

The conjugate of a quaternion \(q\) is denoted \(q^{*}\):

$$q^{*} = w-x\textrm{i}-y\textrm{j}-z\textrm{k}$$
Quaternion quat_conjugate(Quaternion q){
    return (Quaternion){-q.x, -q.y, -q.z, q.w};
}

Quaternion Inverse

The inverse of a quaternion \(q\), denoted \(q^{-1}\), is the conjugate divided by the magnitude squared:

$$q^{-1} = \frac{q^{*}}{||q||^{2}}$$
Quaternion quat_inverse(Quaternion q){
    float m = quat_magnitude(q);
    if(m == 0) return (Quaternion){0, 0, 0, 0};
    m *= m;
    return (Quaternion){-q.x/m, -q.y/m, -q.z/m, q.w/m};
}

For unit quaternions, the conjugate is equal to the inverse.
Multiplying a quaternion by its inverse results in the identity quaternion:

$$q.q^{-1} = (0, 0, 0, 1)$$

Quaternion difference

The difference of two quaternions \(q_{1}\) and \(q_{2}\) is another quaternion \(q_{3}\) that rotates from \(q_{1}\) to \(q_{2}\):

$$q_{3} = q_{1}^{-1}.q_{2}$$
Quaternion quat_difference(Quaternion a, Quaternion b){
    return quat_mul(quat_inverse(a), b);
}

Quaternion Exp and Log

The exponential and the logarithm of a quaternion won't be very useful by themselves, but we will use them to compute other functions later.

Given a quaternion \(q = (x,y,z,w)\) and its vector part \(v = (x,y,z)\), the exponential of that quaternion is also a quaternion, and it's given by this formula:

$$\exp(q) = \exp(w)\begin{pmatrix} \frac{v_{x}}{||v||}\sin(||v||)\\ \frac{v_{y}}{||v||}\sin(||v||)\\ \frac{v_{z}}{||v||}\sin(||v||)\\ \cos(||v||) \end{pmatrix}$$
Quaternion quat_exp(Quaternion q){
    Vector3 v = (Vector3){q.x, q.y, q.z};
    float v_m = Vector3_magnitude(v);
    Vector3 v_n = Vector3_normalize(v);
    float sin_v = sin(v_m);
    float exp_w = exp(q.w);
    return (Quaternion){
        v_n.x*sin_v*exp_w,
        v_n.y*sin_v*exp_w,
        v_n.z*sin_v*exp_w,
        cos(v_m)*exp_w
    };
}

The logarithm of a quaternion is also a quaternion and is given by this formula:

$$\log(q) = \begin{pmatrix} \frac{v_{x}}{||v||}\arccos(\frac{w}{||q||})\\ \frac{v_{y}}{||v||}\arccos(\frac{w}{||q||})\\ \frac{v_{z}}{||v||}\arccos(\frac{w}{||q||})\\ \log(||q||) \end{pmatrix}$$
Quaternion quat_log(Quaternion q){
    Vector3 v = (Vector3){q.x, q.y, q.z};
    Vector3 v_n = Vector3_normalize(v);
    float m = quat_magnitude(q);
    float a = acos(q.w/m);
    return (Quaternion){
        v_n.x*a,
        v_n.y*a,
        v_n.z*a,
        log(m)
    };
}

Quaternion exponentiation

Raising a quaternion to a power results in either a fraction or a multiple of that quaternion. \(q^{2}\) represents twice the rotation of \(q\), and \(q^{\frac{1}{2}}\) represents half of that rotation.

$$q^{n} = \exp(n\log(q))$$
Quaternion quat_pow(Quaternion q, float n){
    return quat_exp(quat_scale(quat_log(q), n));
}

Quaternion slerping

Arguably one of the most important advantages of quaternions, "Slerp" stands for spherical linear interpolation. It's a function that takes three parameters: a quaternion \(q_{1}\), a quaternion \(q_{2}\), and an interpolation parameter \(t\) that goes from \(0\) to \(1\). It gives us an intermediate rotation depending on the value of \(t\).

$$\textrm{slerp}(q_{1}, q_{2}, t) = q_{1}(q_{1}^{-1}q_{2})^{t}$$
Quaternion quat_slerp(Quaternion q1, Quaternion q2, float t){
    t = t < 0 ? 0 : t;
    t = t > 1 ? 1 : t;
    return quat_mul(q1, quat_pow(quat_mul(quat_inverse(q1), q2), t));
}

Here is an interactive demo showing a cube slerping from Euler angle \((0, -\frac{3\pi}{2}, 0)\) to \((0, 0, 0)\):

t =


Well that doesn't look quite right, the cube is doing a \(\frac{3\pi}{2}\) rotation clockwise on the Y axis, when a \(\frac{\pi}{2}\) counterclockwise rotation would have suffised. To fix that we will use the quaternion dot product.

Quaternion dot product

The quaternion dot product (not to be confused with quaternion multiplication) is analogous to the vector dot product, in the sense that it is a measure of how similar two quaternions are.

$$q_{1} = a+b\textrm{i}+c\textrm{j}+d\textrm{k}$$ $$q_{2} = e+f\textrm{i}+g\textrm{j}+h\textrm{k}$$ $$ q_{1} \bullet q_{2} = ae+bf+cg+dh $$
float quat_dot(Quaternion q1, Quaternion q2){
    return q1.x*q2.x + q1.y*q2.y + q1.z*q2.z + q1.w*q2.w;
}

Fixed quaternion slerping

With the dot product we check if two quaternions are on opposite hemispheres, if they are we interpolate between \(q_{1}\) and \(-q_{2}\) instead of \(q_{2}\):

$$ \textrm{slerp}(q_{1}, q_{2}, t) = \begin{cases} q_{1}(q_{1}^{-1}q_{2})^{t},& \text{if } q_{1} \bullet q_{2} > 0\\ q_{1}(q_{1}^{-1}(-q_{2}))^{t},& \text{otherwise} \end{cases} $$
Quaternion quat_slerp(Quaternion q1, Quaternion q2, float t){
    t = t < 0 ? 0 : t;
    t = t > 1 ? 1 : t;
    if(quat_dot(q1, q2) < 0) q2 = quat_scale(q2, -1);
    return quat_mul(q1, quat_pow(quat_mul(quat_inverse(q1), q2), t));
}

Here is another demo showing the cube slerping from Euler angle \((0, -\frac{3\pi}{2}, 0)\) to \((0, 0, 0)\), this time taking the shortest path for the rotation:

t =


Source code

Further reading