- Introduction
- Quaternion magnitude
- Quaternion normalization
- Identity quaternion
- Scaling a quaternion
- Quaternion multiplication
- Euler angle to Quaternion
- Axis angle to Quaternion
- Quaternion to rotation matrix
- Quaternion conjugate
- Rotating a vector by a quaternion
- Quaternion inverse
- Quaternion difference
- Quaternion Exp and Log
- Quaternion exponentiation
- Quaternion slerping
- Quaternion dot product
- Fixed quaternion slerping
- Source code
- Further reading

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:

- 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 = euler_to_quat(, , )

q2 = euler_to_quat(, , )

```
```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_.x};
}

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 :

t = 0.00

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:

t = 0.00

- Let's remove Quaternions from every 3D Engine
- How to Fix Gimbal Lock in N-Dimensions
- Understanding Slerp, Then Not Using It
- Representing Attitude: Euler Angles, Unit Quaternions, and RotationVectors
- Visualizing quaternions, an explorable video series