Quaternions
3. Quaternions
Quaternions are seen as a generalization of complex numbers. A quaternion is a 4-tuple \( q = w + xi + yj + zk \) where \( i^2 = j^2 = k^2 = ijk = -1 \). Quaternions are used to represent rotations in 3D space. But before we get to that, let's understand some algebriac properties of quaternions.
3.1 A little bit of algebra with quaternions
The product of two quaternions \( a = a_0 + a_1i + a_2j + a_3k \) and \( b = b_0 + b_1i + b_2j + b_3k \) is:
\( ab = (a_0b_0 - a_1b_1 - a_2b_2 - a_3b_3) + (a_0b_1 + a_1b_0 + a_2b_3 - a_3b_2)i + (a_0b_2 - a_1b_3 + a_2b_0 + a_3b_1)j + (a_0b_3 + a_1b_2 - a_2b_1 + a_3b_0)k \)
Now if \( a_0 = b_0, a_1 = -b_1, a_2 = -b_2, a_3 = -b_3 \), then \( ab = a_0^2 + a_1^2 + a_2^2 + a_3^2 \) which is equal to the square of the norm of the vector \(
\begin{bmatrix}
a_0 \\ a_1 \\ a_2 \\ a_3
\end{bmatrix}
\).
Now we will create two quaternions from a 3D unit vector \( q = \)\(
\begin{bmatrix}
q_1 \\ q_2 \\ q_3
\end{bmatrix}
\) and an angle \( θ \) like so:Note that \( q_θ . q_θ^* = 1 \) because if \(
\begin{bmatrix}
q_1 \\ q_2 \\ q_3
\end{bmatrix}
\) is a unit vector, then \(
\begin{bmatrix}
cos(θ/2) \\ q_1.sin(θ/2) \\ q_2.sin(θ/2) \\ q_3.sin(θ/2)
\end{bmatrix}
\) is also a unit vector for any \( θ \). Thus, \( q_θ^* \) is the inverse of \( q_θ \) and we can write it as \( q_θ^* = q_θ^{-1} \).
\( q_θ = cos(θ/2) + sin(θ/2).(q_1i + q_2j + q_3k) \)
\( q_θ^* = cos(θ/2) - sin(θ/2).(q_1i + q_2j + q_3k) \)
3.2 Rotating a vector using quaternions
Let's say we want to rotate a vector \( x = \)\(
\begin{bmatrix}
x_1 \\ x_2 \\ x_3
\end{bmatrix}
\) around an axis represented by a unit vector \( q = \)\(
\begin{bmatrix}
q_1 \\ q_2 \\ q_3
\end{bmatrix}
\) by an angle \( \theta \). Let's denote the output vector as \( y = \)\(
\begin{bmatrix}
y_1 \\ y_2 \\ y_3
\end{bmatrix}
\). Then:
\(
\begin{bmatrix}
0 \\ y_1 \\ y_2 \\ y_3
\end{bmatrix}
\)\( = q_θ \)\(
\begin{bmatrix}
0 \\ x_1 \\ x_2 \\ x_3
\end{bmatrix}
\)\( q_θ^{-1} \)
Here we took the liberty to represent a quaternion \( a + bi + cj + dk \) as a vector \(
\begin{bmatrix}
a \\ b \\ c \\ d
\end{bmatrix}
\). You can verify the above equation by multiplying the quaternions on the right hand side.
It would be ideal to spend some time understanding the above equation. However we will move on for now. We will get an understanding of this equation when we cover rotations using geometric algebra.
3.3 Matrix representation of a quaternion
A quaternion \( q = w + xi + yj + zk \) can be represented as a matrix:\(
\begin{bmatrix}
1-2s(y^2+z^2) & 2s(xy - zw) & 2s(xz+wy) \\ 2s(xy+zw) & 1-2s(x^2+z^2) & 2s(yz-wx) \\ 2s(xz-wy) & 2s(yz+wx) & 1-2s(x^2+y^2)
\end{bmatrix}
\)where \( s = ||q||^{-2} \) and \( ||q|| = \sqrt{w^2 + x^2 + y^2 + z^2} \).
Here is a basic implementation of a quaternion:
class Quaternion:
def __init__(self, w:float, x:float, y:float, z:float):
self.w = w
self.x = x
self.y = y
self.z = z
def norm(self):
return np.sqrt(self.w**2 + self.x**2 + self.y**2 + self.z**2)
def unit(self):
n = self.norm()
return Quaternion(self.w/n, self.x/n, self.y/n, self.z/n)
def conjugate(self):
return Quaternion(self.w, -self.x, -self.y, -self.z)
def __mul__(self, other):
w = self.w*other.w - self.x*other.x - self.y*other.y - self.z*other.z
x = self.w*other.x + self.x*other.w + self.y*other.z - self.z*other.y
y = self.w*other.y - self.x*other.z + self.y*other.w + self.z*other.x
z = self.w*other.z + self.x*other.y - self.y*other.x + self.z*other.w
return Quaternion(w, x, y, z)
def to_matrix(self):
u = self.unit()
w,x,y,z = u.w, u.x, u.y, u.z
xx, yy, zz = x*x, y*y, z*z
xy, xz, yz = x*y, x*z, y*z
wx, wy, wz = w*x, w*y, w*z
m = [[1 - 2*(yy + zz), 2*(xy - wz), 2*(xz + wy)],
[2*(xy + wz), 1 - 2*(xx + zz), 2*(yz - wx)],
[2*(xz - wy), 2*(yz + wx), 1 - 2*(xx + yy)]]
return np.array(m)
Here is how one would use it to rotate a vector:
def rotate_vector(
axis_x:float, axis_y:float, axis_z:float,
angle:float,
vector_x:float, vector_y:float, vector_z:float
):
cos_angle, sin_angle = np.cos(angle/2), np.sin(angle/2)
# make the axis a unit vector
axis_norm = np.sqrt(axis_x**2 + axis_y**2 + axis_z**2)
axis_x, axis_y, axis_z = axis_x/axis_norm, axis_y/axis_norm, axis_z/axis_norm
q = Quaternion(cos_angle, sin_angle*axis_x, sin_angle*axis_y, sin_angle*axis_z)
v = Quaternion(0, vector_x, vector_y, vector_z)
y = q * v * q.conjugate()
return y.x, y.y, y.z