# Rotations with quaternions

A quaternion is a 4 dimensional complex-like number, it has four components, three of which are the "imaginary" part.

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 but I like to put 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):

 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 ) 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 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 (). If we want to apply a rotation then a rotation , the resulting rotation is:

Quaternion multiplication looks like this:

 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:

• 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:

 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); 
q1 = (, , )
q2 = (, , )

## 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:

The model matrix itself looks like this:

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

We convert a quaternion to a rotation matrix like this:

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 is denoted :

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

## Quaternion Inverse

The inverse of a quaternion , denoted , is the conjugate divided by the magnitude squared:

 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:

## Quaternion difference

The difference of two quaternions and is another quaternion that rotates from to :

 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 and its vector part , the exponential of that quaternion is also a quaternion, and it's given by this formula:

 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:

 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. represents twice the rotation of , and represents half of that rotation.

 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 , a quaternion , and an interpolation parameter that goes from to . It gives us an intermediate rotation depending on the value of .

 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 to :

t =

Well that doesn't look quite right, the cube is doing a rotation clockwise on the Y axis, when a 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.

 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 and instead of :

 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 to , this time taking the shortest path for the rotation:

t =