Technology

Rotations with quaternions




Rotations with quaternions


A quaternion is a 4 dimensional complex-like number, it has four components, three of which are the “imaginary” part.
$$ q = a+btextrm{i}+ctextrm{j}+dtextrm{k} $$
$$ q = (b,c,d, a) $$
$$ 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};

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);
}

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

Quaternion quat_normalize(Quaternion q){
    float m = quat_magnitude(q);
    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.

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 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};
}

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+btextrm{i}+ctextrm{j}+dtextrm{k}$$
$$q_{2} = e+ftextrm{i}+gtextrm{j}+htextrm{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
    };
}

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

  • Euler angles suffer from gimbal lock
  • Interpolating between two Euler angles lead to weird results

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);

When doing 3D rendering, we usually pass an MVP (Model View Projection) matrix to a shader to proprely 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 4×4 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
    };
}

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

$$q^{*} = a-btextrm{i}-ctextrm{j}-dtextrm{k}$$

Quaternion quat_conjugate(Quaternion q){
    return (Quaternion){-q.x, -q.y, -q.z, q.w};
}

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}; // avoid division by 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)$$

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);
}

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)
    };
}

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(nlog(q))$$

Quaternion quat_pow(Quaternion q, float n){
    return quat_exp(quat_scale(quat_log(q), n));
}

Arguably one of the most important advantages of quaternions, “Slerp” stands for spherical linear interpolation. It’s a function thats 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  1 ? 1 : t;
    return quat_mul(q1, quat_pow(quat_mul(quat_inverse(q1), q2), t));
}

Here is an animation showing a cube slerping from Euler angle ((0, frac{5pi}{4}, frac{pi}{4})) to ((0, 0, 0)):

Transform obj;

Quaternion start_rotation = euler_to_quat((Vector3){0, 5*PI/4, PI/4});
Quaternion target_rotation = euler_to_quat((Vector3){0, 0, 0})
obj.rotation = start_rotation;

float t = 0;

// Main loop
while(1){
    ...
    obj.rotation = quat_slerp(start_rotation, target_rotation, t);
    t += delta_time;
    ...
}

Related Articles

Back to top button