Skip to content

nicholasferguson/Display.Quaternion.Rotation

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

71 Commits
 
 
 
 
 
 
 
 

Repository files navigation

This is a discussion for using Quaternion math for rotations.
It uses some computer science terminology: 'overloads' and 'interface"

Python version 3.8.6 on Windows 10.

This is documentation for two python files: q3a.py and plotPrims.py

Some computer science terminology will be used to explain quaternions: 'overload' and 'interface'

In the fig below, two quaternions have been used. One quaternion represents a 3D point, displayed as a highlighted vertex in a light blue prism. A 2nd quaternion, used as a rotator, rotates this 3D point. Result of this rotation is displayed as a highlighted vertex in a light purple prism. In fig below, rotation is about 'Z' axis.

from IPython.display import Image 
# A highlighted point displayed as a vertex in light blue prism is rotated to a 
# highlighted point, represented as a vertex in a light purple prism
# Rotation is about a 'Z' axis.

"Result"

Some computer science terminology will be used to explain quaternions: 'overload' and 'interface'

In this write up, quaternions use a format of [Scalar, X, Y, Z] and not [X, Y, Z ,Scalar]. First format is called Hamilton's and 2nd format is called JPL. This article will use Hamilton's format.

We present quaternions as having 'overloads'.

First overload is a quaternion that is a rotator. In our example, rotator is for Z axis only.

[Scalar, X, Y, Z] as a 'rotator' becomes

[Scalar = cos of rotating angle in radians,
X = rotating angle about x axis = 0 radians,
Y = rotating angle about y axis = 0 radians,
Z = rotating angle about z axis = sin of Scalar radians]

Note: This is notebook, so reader can change python code to rotate another axis:X Y Z, combination or all. Also change angle of rotation.

2nd overload is a quaternion as a 'pure quaternion'. In our example it represents coordinates for a 3D point.

[scalar, X, Y, Z] becomes [0, X, Y, Z].

In this example, 'X Y Z' represents a point's coordinates in 3D space.

Note: plotPrism.py is hard coded to properly display a fixed range of X Y Z as coordinate values. In next version it will handle a wider range. Reader can also change values of X Y Z.

Overall Algorithm: A quaternion rotator multiplied with a pure quaternion, as a 3D point, and conjugate of 'quaternion rotator'. This algorithm results in a 3D point rotated to a new position. Result is a 'pure quaternion'of a 3D point.

Wire frames of prisms display these 3D points as highlighted vertices.
Initial (X Y Z) is a vertex of a prism.
Rotated (X Y Z) is a vertex of a second prism.
Both wire frames have a 3D origin, as a vertex (0,0,0).

To prove that wire frames are correct representations, run python code, outside of jupyter, and rotate image along its axes.

# Python version 3.8.6 Windows 10
%matplotlib notebook   
#to get plt.show() from plotPrism.py into notebook
#calling it a second time may prevent some graphics errors
%matplotlib notebook  
import math
import numpy as np
import matplotlib.pyplot as plt
import os
import sys
plotprism_path = os.path.abspath(os.path.join(''))
print("plotPrism.py must be in following directory")
print("Path string will be from your own computer.  In README.MD value is from my computer")
print(plotprism_path)
if plotprism_path not in sys.path:
    sys.path.append(plotprism_path)
    
plotPrism.py must be in following directory
Path string will be from your own computer.  In README.MD value is from my computer
C:\temp\Copula\_MathReview\01.GA\___WriteUp.Quaternion.Rotate\Display.Quaternion.Rotation-main\Quaternion.python
from plotPrism import plot_prism #plotprism.py


def normalize(v, tolerance=0.00001):
    mag2 = sum(n * n for n in v)
    if abs(mag2 - 1.0) > tolerance:
        mag = math.sqrt(mag2)
        v = tuple(n / mag for n in v)
    return np.array(v)

Key math algorithms.

Hamilton solved several issues to make this work. He devised quaternion structures, their overloads and a complex number interface for quaternion multiplication. We have covered two quaternion overloads, that are used in rotations.

This complex number interface reduces a multiplication of two quaternions to a quaternion result. In our example it's used to multiply a 'rotator quaternion' and 'pure quaternion' and conjugate of a 'rotator quaternion'. Result is a 'pure quaternion'

'Complex Number Interface' to generate multiplication algorithm for quaternions

When a quaternion multiplies a 2nd quaternion, Hamilton used 'Hamilton product'.
Quaternion 1: [a0,a1,a2,a3]
Quaternion 2: [b0,b1,b2,b3]

Result of multiplying Quaternion 2 by Quaternion 1:
a0b0 + a0b1 + a0b2+ a0b3 +
a1b0 + a1b1 + a1b2+ a1b3 +
a2b0 + a2b1 + a2b2+ a2b3 +
a3b0 + a3b1 + a3b2+ a3b3

In this example, with a 'complex number interface' we could then reduce these multiplication results to a [X Y Z] or 'pure quaternion'.

Generate a 'Complex Number Interface' to get a quaternion multiplication function

A complex number interface also has three axes: i, j, k. They interface to X,Y,Z axis. In our example, this interface will allow us to simplify a quaternion multiplication to a pure quaternion.

First we need some rules for complex number multiplication. With these rules we will substitute pairs of complex numbers with a single complex number. This yields a new quaternion [Scalar, X, Y, Z]. In our example, eventually it yields 3D coordinates for a rotated point [0,X,Y,Z]

Hamilton discovered that
ijk = -1 and of course ii, jj, kk = -1.
It's not commutative. kj neq jk. But kj eq -jk (neq is not equal)

First step in this interface is to create a translation where pairs of complex number equal a single complex number.

ijk = -1 # Hamilton's discovery
iijk = -i # we multiply each side by 'i'
-jk = -i # ii = -1
jk = i AND kj = -i

ijk = -1
-ikj = -1 # jk = -kj
kij = -i # -ik = ki
kkij = -k # we multiply both sides by k
-ij = -k ij=k AND ji = -k

ijk = -1
-jik = -1
-jjik = -j
ik = -j AND ki = j

Now we redo multiplication of two quaternions (regardless of overload)
Note: we need additional layers of indices to work with our complex number interface.
Quaternion 1: [a0,a1i,a2j,a3k] This is [Scalar, i, j, k]
Quaternion 2: [b0,b1i,b2j,b3k]

Result of multiplying Quaternion 2 by Quaternion 1:

a0b0 + a0b1i + a0b2j+ a0b3k + # we multiply Quaternion 2 by a0
a1ib0 + a1ib1i + a1ib2j+ a1ib3k +
a2jb0 + a2jb1i + a2jb2j+ a2jb3k +
a3kb0 + a3kb1i + a3kb2j+ a3kb3k

We restate this result by collecting complex variables together.

a0b0 + (a0b1)i + (a0b2)j+ (a0b3)k +
(a1b0)i + (a1b1)ii + (a1b2)ij+ (a1b3)ik +
(a2b0)j + (a2b1)ji + (a2b2)jj+ (a2b3)jk +
(a3b0)k + (a3b1)ki + (a3b2)kj+ (a3b3)kk

Then we substitute pairs of complex variables with a single complex variable, from our translations.

a0b0 + (a0b1)i + (a0b2)j+ (a0b3)k +
(a1b0)i + (a1b1)-1 + (a1b2)k+ (a1b3)-j +
(a2b0)j + (a2b1)-k + (a2b2)-1+ (a2b3)i +
(a3b0)k + (a3b1)j + (a3b2)-i+ (a3b3)-1

Finally, we group by scalar or i, j, k

a0b0 - (a1b1) - (a2b2) -(a3b3) # Scalars
((a0b1) +a1b0) + (a2b3) - (a3b2))i # i's
((a0b2) - (a1b3) + (a2b0) + (a3b1)j # j's
((a0b3) + (a1b2) - (a2b1) + (a3b0))k # k's

This give us an interface to generate our quaternion multiplication function for our python file:

We substitute pq[0] for a0 and r[0] for b0 etc... This multiplication function will now return a quaternion of [Scalar, X, Y ,Z]
When we multiply that result by a 'conjugate of rotator', we get a pure quaternion of [Scalar==0, X, Y ,Z]

def quaternion_mult(r,pq):
    return [pq[0]*r[0]-pq[1]*r[1]-pq[2]*r[2]-pq[3]*r[3], # Scalar
            pq[0]*r[1]+pq[1]*r[0]+pq[2]*r[3]-pq[3]*r[2], # X
            pq[0]*r[2]-pq[1]*r[3]+pq[2]*r[0]+pq[3]*r[1], # Y
            pq[0]*r[3]+pq[1]*r[2]-pq[2]*r[1]+pq[3]*r[0]] # Z

We complete our rotation algorithm. In below code block, you can comment out r_conj in quaternion_mult to examine results in plot_prism

# pq is pure quaternion
# rq is rotator quaternion
def point_rotation_by_quaternion(pq,rq):
    r = normalize(rq)
    r_conj = [r[0],-1*r[1],-1*r[2],-1*r[3]]
    return quaternion_mult(quaternion_mult(r,pq),r_conj)

Below are our two overloaded quaternions, with their angle of rotation.
These values can be changed by reader.
Note that plotPrism.py, currently can only handle a narrow range of values. This will be extended on a next version of code.
Note: A formula for rotating a point (x, y) in 2D is given by:

x': = x * cos (angle) - y * sin (angle)
y': = y * cos (angle) + x * sin (angle)

degrees = math.pi/180;
rot = 180*degrees;
# for rotating about an axis.
w = math.cos(rot/2.);
ax = math.sin(rot/2.);
# quaternion format is [scalar, x, y ,z]
pq = [0, 1, 2, 3]  # pure quaternion as 3D coordinates.  Scalar is zero.
rq = [w, 0, 0, ax]  # play with ax on different axis: x,y,z  Change values of w and ax.

We now apply rotation algorithms

pq2 = point_rotation_by_quaternion(pq,rq)
print(pq2)
[0.0, -0.9999999999999998, -2.0, 3.0]

We now prep for display purposes

# Fill in a start prism wire frame, using point and a 3D origin as vertices
# Highlight point as vertex

p1 = np.array([0.,0.,0.])
p2 = np.array([pq[1],0.,0.])
p3 = np.array([0.,pq[2],0.])
p4 = np.array([0.,0.,pq[3]])

prism1 = [
    (p1[0],p1[1],p1[2]), (p2[0],p2[1],p2[2]), (p3[0],p3[1],p3[2]), (p4[0],p4[1],p4[2])
]

# Fill in a rotated prism wire frame, using point and a 3D origin as vertices
# Highlight point as vertex

p1A = np.array([0.,0.,0.])
p2A = np.array([pq2[1],0.,0.])
p3A = np.array([0.,pq2[2],0.])
p4A = np.array([0.,0.,pq2[3]])

prism2 = [
    (p1A[0],p1A[1],p1A[2]), (p2A[0],p2A[1],p2A[2]), (p3A[0],p3A[1],p3A[2]), (p4A[0],p4A[1],p4A[2])
]

# to prove that plot is correct, you need to rotate image along its axes.

plot_prism(prism1,prism2)
<IPython.core.display.Javascript object>

This image can be rotated by running files q3a.py and plotPrism.py in a python ide like 'Sublime Text'

Notes: This point being rotated could be a significant point of a geometric object. For a sphere, it could represent a center. Quaternion rotation is meant to preserve shape.