Oregon Curriculum Network

Discovering Math with Python

Chapter 11: QUATERNIONS

Quaternions were invented by Sir William Rowan Hamilton around 1843 and were considered a breakthrough.

In subsequent chapters Willard Gibbs and Oliver Heaviside came up with a vector concept that proved easier to use for many of the physics applications for which quaternions had originally been proposed.

However quaternions and vectors together have become chief tools in accomplishing rotation, in computer graphics and games, robotics, rocketry. We may have entered the realm of that proverbially most difficult of disciplines: rocket science.

Quaternions have some advantages over rotation matrices. It's easier to slice up a rotation in a process called SLERP.

You might think of quaternions as "vectors on steroids" if that helps. They have some properties in common with vectors, and are conceived to have a "vector part" somewhat like complex numbers have a "real part". However, like complex numbers, they're considered numbers in their own right, actually a superset of the complex, which in turn contain the reals and so on down to N, the counting numbers (natural numbers).

Unit quaternions (w, x, y, z) with w**2 + x**2 + y**2 + z**2 == 1, form a group under multiplication. Any two unit quaternions, when multiplied, produce a unit quaternion (Closure), and every unit quaternion has an inverse (w, -x, -y, -z), such that q * q**-1 gives the unit quaternion (1, 0, 0, 0).

Indeed the elements i, j, k from which quaternions are made, abet the i of complex number fame, with two more 2nd roots of -1. All three, and their three inverses, engage in a kind of dance.

Group elements are: {i, j, k, 1, -1, -i, -j, -k}. Every product of two of these elements, is in this set (closure); associativity holds; every element has an inverse such that the two give a product of 1, and 1 serves as the neutral (identity) element.

Just the identities we've been given: i**2 = j**2 = k**2 = i * j * k = -1 are suffient to derive the above Cayley Table. This table, in turn enables use to flesh out the __mul__ method, by taking sixteen products (w, x, y, z times each w', x', y', z') and substituting for products such as -i * k and k * j. Some of the derivations for are table are shown above.

What we'll plan to do, when rotating point vectors of polyhedron P, by n degrees around unit rotation vector q, is initialize the corresponding unit quaternion rQ using function rotation(n, q) and then multiply every point vector of P in a "sandwich" between rQ and its inverse:

new_Pv = rQ * Pv * ~rQ (how we apply rotations to each position vector Pv).

Remember our Polyhedron objects, as data structures, start with a set of faces pegged to vertexes, expressed as a dict of position Vectors. We used both xyz and quadray notation to map the same set of points.

Below is the source code for doing that. Understanding multiplication involves developing the multiplication table for i, j, k which play the role of basis vectors in this language game. i**2 == j**2 == k**2 == -1, and i * j * k == -1 as well.


In [1]:
from math import cos, sin, radians, pi
from qrays import Vector
import unittest

class Quaternion:
    
    def __init__(self, w, x, y, z):
        """
        w is the scalar part;
        x,y,z comprise a vector part as 
        the coefficients of i,j,k respectively
        where i,j,k are the three basis vectors
        described by Hamilton such that 
        i**2 == j**2 == k**2 == -1, and 
        i * j * k == -1 as well.
        """
        self.w = w
        self.x = x
        self.y = y
        self.z = z
            
    def __mul__(self, other):
        """
        Derived by inter-multiplying all four terms in 
        each Quaternion to get 16 products and then 
        simplifying according to such rules as ij = k,
        ji = -k.
        
        See:  https://youtu.be/jlskQDR8-bY Mathoma
        'Quaternions explained briefly'
        """
        a, b, c, d = self.w, self.x, self.y, self.z
        e, f, g, h = other.w, other.x, other.y, other.z
        w = (a*e - b*f - c*g - d*h)
        x = (a*f + b*e + c*h - d*g)
        y = (a*g - b*h + c*e + d*f)
        z = (a*h + b*g - c*f + d*e)
        return Quaternion(w, x, y, z)
    
    def __invert__(self):
        return Quaternion( self.w,
                          -self.x,
                          -self.y,
                          -self.z)
        
    def vector(self):
        return Vector((self.x, self.y, self.z))
        
    def __eq__(self, other):
        tolerance = 1e-8
        return (abs(self.w - other.w) < tolerance and
                abs(self.x - other.x) < tolerance and 
                abs(self.y - other.y) < tolerance and
                abs(self.z - other.z) < tolerance)
         
    def __repr__(self):
        return "Quaternion({},{},{},{})".format(self.w, 
                          self.x,
                          self.y,
                          self.z)

def rotator(uV, α):
    w = cos(α/2)
    x = uV.x * sin(α/2)
    y = uV.y * sin(α/2)
    z = uV.z * sin(α/2)
    return Quaternion(w, x, y, z)
    
class TestQuaternion(unittest.TestCase):
    
    def test_inverse(self):
        Q = Quaternion(0, 0, -1, 0)
        inverse = ~Q
        self.assertEqual(Q * inverse, Quaternion(1,0,0,0))
            
    def test_x(self):
        rV = Quaternion(0,1,0,0)
        rQ = rotator(Vector((0,1,0)), pi)
        newQ = rQ * rV * ~rQ
        self.assertTrue(newQ.vector() == Vector((-1, 0, 0)))
        
    def test_360(self):
        rV = Quaternion(0,1,0,0)
        one_deg = radians(1)
        rQ = rotator(Vector((0,1,0)), one_deg)
        for _ in range(360):
            rV = rQ * rV * ~rQ
        self.assertTrue(rV.vector() == Vector((1, 0, 0)))
        
a = TestQuaternion()
# unittest.TextTestRunner().run(a)

suite = unittest.TestLoader().loadTestsFromModule(a)
unittest.TextTestRunner().run(suite)


..
----------------------------------------------------------------------
Ran 2 tests in 0.008s

OK
Out[1]:
<unittest.runner.TextTestResult run=2 errors=0 failures=0>

Back to Chapter 10: Complex Numbers
Continue to Chapter 12: The Mandelbrot Set
Introduction / Table of Contents