A quaternion is a four dimensional complex-like number. It has two parts: an imaginary (or vector) part with three components, and a real (or scalar) part with one component.
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};
A quaternion with zero real (or scalar) term is called a pure quaternion.
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);
}
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.
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 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 (). 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
};
}
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:
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 axis_angle_to_quat(Vector3 axis, float angle){
float s = sin(angle/2);
return (Quaternion){
axis.x*s,
axis.y*s,
axis.z*s,
cos(angle/2)
};
}
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 quat_to_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 is denoted :
Quaternion quat_conjugate(Quaternion q){
return (Quaternion){-q.x, -q.y, -q.z, q.w};
}
We can use a quaternion to rotate a vector directly without converting it to a matrix. First we convert the 3D vector to a pure quaternion, then we pre-multiply it by and post-multiply it by the conjugate :
Vector3 rotate_vector(Vector3 v, Quaternion q){
Quaternion v_ = (Quaternion){v.x, v.y, v.z, 0};
v_ = quat_mul(quat_mul(q, v_), quat_conjugate(q));
return (Vector3){v_.x, v_.y, v_.z};
}
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:
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);
}
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)
};
}
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));
}
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 :
Well that doesn't look quite right, the cube is doing a rotation counterclockwise on the Y axis, when a clockwise rotation would have suffised. To fix that we will use the 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;
}
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: