A great video explaining and visualizing quaternions by 3Blue1Brown: Visualizing quaternions (4d numbers) with stereographic projection
It’s nice to see geometric algebra being used for practical applications. It’s a very elegant mathematical framework for doing a wide variety of analysis in multidimensional spaces. I don’t often get to use it for my daily work, but I try to keep a working knowledge of it in the back of my mind just in case I ever come across an application that could possibly be expressed a little more clearly in this framework.
You made me think of this Visual Basic code I had to write in 2002, after finding out the hard way that the usual Trig rotation subroutines did not work for articulating a skeleton. If I recall correctly then it’s like after a couple of rotations there is no easy way of knowing where body parts are located. The main subroutines were:
Sub XYZtoQuat()
'Euler To Quaternions: Roll = XA, Pitch = YA, Yaw = ZA
'Calculate trig identities
CR = Cos(XA / 2)
CP = Cos(YA / 2)
CY = Cos(ZA / 2)
SR = Sin(XA / 2)
SP = Sin(YA / 2)
SY = Sin(ZA / 2)
CPCY = CP * CY
SPSY = SP * SY
QW = CR * CPCY + SR * SPSY
QX = SR * CPCY - CR * SPSY
QY = CR * SP * CY + SR * CP * SY
QZ = CR * CP * SY - SR * SP * CY
End Sub
Sub MulQuat()
'If object is rotated by r1, then by r2, then the composed rotation (r2 o r1) is represented
'by the quaternion (q2 * q1), or by the matrix (m1*m2).
'Note that matrix composition is reversed because matrices and vectors are represented rowwise.
'
'You have to be careful that your quaternion represents an absolute and not a relative rotation.
'You can think of a relative rotation as a rotation from the previous (intermediate) orientation
'and an absolute rotation as the rotation from the initial orientation.
'This becomes clearer if you think of the q2 quaternion orientation in Figure 3 as a relative rotation,
'since it moved with respect to the q1 orientation. To get an absolute rotation of a given quaternion,
'you can just multiply the current relative orientation by a previous absolute one.
'The initial orientation of an object can be represented as a multiplication identity [1, (0, 0, 0)].
'This means that the first orientation is always an absolute one, because q = qidentity q
'From: Rotating Objects Using Quaternions, by Nick Bobick
'
'Efficient quaternion multiplication.
'QuatMul(QUAT *q1, QUAT *q2, QUAT *res){
'
A = (Q1W + Q1X) * (Q2W + Q2X)
B = (Q1Z - Q1Y) * (Q2Y - Q2Z)
C = (Q1W - Q1X) * (Q2Y + Q2Z)
D = (Q1Y + Q1Z) * (Q2W - Q2X)
E = (Q1X + Q1Z) * (Q2X + Q2Y)
F = (Q1X - Q1Z) * (Q2X - Q2Y)
G = (Q1W + Q1Y) * (Q2W - Q2Z)
H = (Q1W - Q1Y) * (Q2W + Q2Z)
'
ResW = B + (-E - F + G + H) / 2
ResX = A - (E + F + G + H) / 2
ResY = C + (E - F + G - H) / 2
ResZ = D + (E - F - G + H) / 2
QW = B + (-E - F + G + H) / 2
QX = A - (E + F + G + H) / 2
QY = C + (E - F + G - H) / 2
QZ = D + (E - F - G + H) / 2
End Sub
Sub QuatToMatrix()
'Picture1.Print QW; QX; QY; QZ
' Action: Converts quaternions representation of a rotation to a matrix representation
' Params: QW, QX, QY, QZ
' Returns: M(4x4 matrix) is filled
' Comments: remember matrix (in OGL) is represented in COLUMN major form
X2 = QX + QX
Y2 = QY + QY
Z2 = QZ + QZ
XX = QX * X2
XY = QX * Y2
XZ = QX * Z2
YY = QY * Y2
YZ = QY * Z2
ZZ = QZ * Z2
WX = QW * X2
WY = QW * Y2
WZ = QW * Z2
M(0, 0) = 1 - (YY + ZZ)
M(0, 1) = XY - WZ
M(0, 2) = XZ + WY
M(0, 3) = 0
M(1, 0) = XY + WZ
M(1, 1) = 1 - (XX + ZZ)
M(1, 2) = YZ - WX
M(1, 3) = 0
M(2, 0) = XZ - WY
M(2, 1) = YZ + WX
M(2, 2) = 1 - (XX + YY)
M(2, 3) = 0
M(3, 0) = 0
M(3, 1) = 0
M(3, 2) = 0
M(3, 3) = 1
End Sub
Sub MatrixToQuat()
'Matrix to quaternion conversion
'MatToQuat(float m[4][4], QUAT * quat)
Nxt(0) = 1
Nxt(1) = 2
Nxt(2) = 0
'{1, 2, 0};
TR = M(0, 0) + M(1, 1) + M(2, 2)
'Check the diagonal for Positive or Negative
If TR > 0 Then GoTo DiagPositive
I = 0
If M(1, 1) > M(0, 0) Then I = 1
If M(2, 2) > M(I, I) Then I = 2
J = Nxt(I)
K = Nxt(J)
S = Sqr((M(I, I) - (M(J, J) + M(K, K))) + 1)
Q(I) = S * 0.5
If S <> 0 Then S = 0.5 / S
Q(3) = (M(J, K) - M(K, J)) * S
Q(J) = (M(I, J) + M(J, I)) * S
Q(K) = (M(I, K) + M(K, I)) * S
QX = Q(0)
QY = Q(1)
QZ = Q(2)
QW = Q(3)
GoTo EndCnv
DiagPositive:
S = Sqr(TR + 1)
QW = S / 2
S = 0.5 / S
QX = (M(1, 2) - M(2, 1)) * S
QY = (M(2, 0) - M(0, 2)) * S
QZ = (M(0, 1) - M(1, 0)) * S
EndCnv:
End Sub
Sub DrawMatrix()
Picture1.Print ""
Picture1.Print "matrix"; N
For Y = 0 To 3
For X = 0 To 3
Picture1.Print Format(M(Y, X), "###0.####"),
Next X
Picture1.Print ""
Next Y
End Sub
Sub WorldToMatrix()
M(0, 0) = X
M(1, 0) = 0
M(2, 0) = 0
M(3, 0) = 0
M(0, 1) = 0
M(1, 1) = Y
M(2, 1) = 0
M(3, 1) = 0
M(0, 2) = 0
M(1, 2) = 0
M(2, 2) = Z
M(3, 2) = 0
M(0, 3) = 0
M(1, 3) = 0
M(2, 3) = 0
M(3, 3) = 1
End Sub
Sub NormalizeQuat()
If QW = 0 Then QW = 1
Magnitude = (QW * QW) + (QX * QX) + (QY * QY) + (QZ * QZ)
QW = QW / Magnitude
QX = QX / Magnitude
QY = QY / Magnitude
QZ = QZ / Magnitude
End Sub
Private Sub SaveSettings_Click()
Open "skeleton.set" For Output As #2
For N = 0 To TotalNodes
Print #2, NodeRoll(N)
Print #2, NodeRollMin(N)
Print #2, NodeRollMax(N)
Print #2, NodePitch(N)
Print #2, NodePitchMin(N)
Print #2, NodePitchMax(N)
Print #2, NodeYaw(N)
Print #2, NodeYawMin(N)
Print #2, NodeYawMax(N)
Print #2, NodeX(N)
Print #2, NodeY(N)
Print #2, NodeZ(N)
Print #2, NodeQW(N)
Print #2, NodeQX(N)
Print #2, NodeQY(N)
Print #2, NodeQZ(N)
Next N
Close #2
End Sub
Private Sub LoadSettings_Click()
On Error GoTo ErrorRet
Open "skeleton.set" For Input As #2
For N = 0 To TotalNodes
Input #2, NodeRoll(N)
Input #2, NodeRollMin(N)
Input #2, NodeRollMax(N)
Input #2, NodePitch(N)
Input #2, NodePitchMin(N)
Input #2, NodePitchMax(N)
Input #2, NodeYaw(N)
Input #2, NodeYawMin(N)
Input #2, NodeYawMax(N)
Input #2, NodeX(N)
Input #2, NodeY(N)
Input #2, NodeZ(N)
Input #2, NodeQW(N)
Input #2, NodeQX(N)
Input #2, NodeQY(N)
Input #2, NodeQZ(N)
Next N
Close #2
Call ChangeNode
Call Reposition
ErrorRet:
End Sub
Sub MatrixToWorld()
XW = (X * M(0, 0)) + (Y * M(1, 0)) + (Z * M(2, 0))
YW = (X * M(0, 1)) + (Y * M(1, 1)) + (Z * M(2, 1))
ZW = (X * M(0, 2)) + (Y * M(1, 2)) + (Z * M(2, 2))
End Sub
Can you recommend any Python code for molecular dynamics simulations? To simulate close packing geometry on a non-flat surface I wrote python code that keeps randomly placed particles/cells bound to a sphere, and I now need to apply Lennard-Jones potential between points:
#from random import random
import sys
import math
import pygame
import random
import numpy as np
Pi=3.14159265358979323846264338327950288412
HalfPi=Pi/2
Radian=Pi*2
SphereRadius=200
Screen_Size = Width, Height = (640, 480)
BLACK = (0, 0, 0)
WHITE = (255, 255, 255)
RED = (255, 50, 50)
GREEN = (50, 255, 50)
CIRCLE_RADIUS = 30
# Initialization
pygame.init()
screen = pygame.display.set_mode(Screen_Size)
xc= int(Screen_Size[0]/2)
yc= int(Screen_Size[1]/2)
pygame.display.set_caption('Spherical Membrane')
fps = pygame.time.Clock()
paused = False
# Initialize Environment xe,ye,ze arrays, with ne points.
ne=100
xyze = [0 for i in range(ne)]
sphe = [0 for i in range(ne)]
rval=0.1
# Fill Environment arrays with randomly placed points.
for i in range(ne):
xyze[i] = random.randrange(-100,100),random.randrange(-100,100),random.randrange(-100,100)
# Cartesian to Spherical coordinates.
def cart2sphere(xyz):
hxy = np.hypot(xyz[0], xyz[1])
r = np.hypot(hxy, xyz[2])
el = np.arctan2(xyz[2], hxy)
az = np.arctan2(xyz[1], xyz[0])
return az, el, r
# Spherical to Cartesian coordinates.
def sphere2cart(sp):
theta = sp[2] * np.cos(sp[1])
x = theta * np.cos(sp[0])
y = theta * np.sin(sp[0])
z = sp[2] * np.sin(sp[1])
return x, y, z
def update():
for i in range(ne):
# Keep on a sphere by converting to Spherical then using its radius back again.
s = cart2sphere(xyze[i]) # Cartesian to Spherical.
sphe[i] = s[0], s[1], SphereRadius # Change r to SphereRadius.
xyze[i] = sphere2cart(sphe[i]) # Spherical to Cartesian again.
# Change x,y,z according to Lennard-Jones attraction/repulsion.
# for j in range(1):
# if i != j:
def render():
screen.fill(BLACK)
for i in range(ne):
xyze[i] = sphere2cart(sphe[i]) # Spherical to Cartesian again.
s=sphe[i]
s=[s[0]+rval,s[1],s[2]]
xyz = sphere2cart(s) # Spherical to Cartesian again.
xe=int(xyz[2]+xc) # X location in Environment
ye=int(xyz[1]+yc) # Y location in Environment
pygame.draw.circle(screen, WHITE, [xe, ye], 2, 0)
pygame.display.update()
fps.tick(60)
while True:
for event in pygame.event.get():
if event.type == pygame.QUIT:
pygame.quit()
sys.exit()
if event.type == pygame.KEYUP:
if event.key == pygame.K_SPACE:
paused = not paused
if not paused:
update()
render()
rval += .01I don’t have an answer for python, but I am using an extremely fast Common Lisp library (intended originally for graphics games) that supports quaternions and spherical coordinates: https://github.com/cbaggers/rtg-math (as well as vectors/matricies etc) We are experiencing amazing performance in commercial code with quaternion transformations with this library.
You could look at the lisp code for quaternions and port to python: https://github.com/cbaggers/rtg-math/tree/master/quaternions
The last time I programmed in LISP it was installed into an IBM-XT clone from a 4 inch paper floppy. I experimented with it for awhile but was soon back to programming in an early Basic.
Now that the era of VB being the worlds most used programming language is long over I’m forced to experiment with the now popular Python, which looks like it will maybe remain that way long enough to make it worth becoming proficient in. But perhaps I may soon need to experiment with LISP again.
I took a look for what was available for Python. I found that there is a “pyquaternion” and it properly installed into PyCharm, so I’m set to go with this:
http://kieranwynn.github.io/pyquaternion/
I’m not a natural born math wiz, and am still unsure how you would most efficiently use it to sum forces in the “j” loop, I commented out in code while getting all the rest of the program working:
def update():
for i in range(ne):
# Keep on a sphere by converting to Spherical then using its radius back again.
s = cart2sphere(xyze[i]) # Cartesian to Spherical.
sphe[i] = s[0], s[1], SphereRadius # Change r to SphereRadius.
xyze[i] = sphere2cart(sphe[i]) # Spherical to Cartesian again.
# Change x,y,z according to Lennard-Jones attraction/repulsion.
# for j in range(1):
# if i != j:
For calculating spherical distances between the points I adapted the Haversine Formula, and can adapt to return a 0 to 1 fractional distance.
# ------------------------------------------------------------------------------+
#
# Nathan A. Rooy
# Haversine Formula
# June, 2016
# https://nathanrooy.github.io/posts/2016-09-07/haversine-with-python/
# ------------------------------------------------------------------------------+
# Radian Haversine to calculate distance between two points on a sphere.
# Distance is given and returned in Radians, max distance away is Pi
# Haversine(-HalfPi, 0, HalfPi, 0) = pole to pole = Pi = half around sphere.
import math
Pi=3.141592653589793238462643383279502884197169399375105820974944592307816406286
HalfPi=Pi/2
Radian=Pi*2
def Haversine(lat1, lon1, lat2, lon2):
phi_1 = lat1
phi_2 = lat2
delta_phi = lat2 - lat1
delta_lambda = lon2 - lon1
a = math.sin(delta_phi / 2.0)**2 + math.cos(phi_1) * math.cos(phi_2) * math.sin(delta_lambda / 2.0)**2
return 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
print(Haversine(-HalfPi, 0, HalfPi, 0))
print(Haversine(-HalfPi, 0, HalfPi/2, 0))
print(Haversine(-HalfPi, 0, 0, 0))
print(Haversine(-HalfPi/2, 0, 0, HalfPi))
The above code might be replaced by a quaternion too. But since I’m not a natural born math wiz I’m not sure what the proper way of adding a 3D displacement vector to the x,y,z points. I have a long (re)learning curve ahead. All help up the hill is appreciated.
In this example a small vector of the same value should draw all the points to one side of the sphere. For bubble-like cell migration behavior this “ljp=” line of code would be used alone for displacement/magnitude, without complicating things by adding accelerations as in gasses in a vacuum.
epsilon = 1 # Energy minimum
sigma = 1 # Distance to zero crossing point
print("Lennard-Jones potential at radius of 0.01 to 5")
print(" i r ljp")
print("------------------------------------")
# Lennard-Jones potential
for i in range(1,501):
r=i/100
ljp = 4*epsilon*((sigma/r)**12-(sigma/r)**6)
print("{:3}{:7} {}".format(i, r, ljp))
From experience I learned that at some point along the way towards modeling a human is the frustrating problem of the usual methods not working as expected, to overcome. An easy quaternion example of cells packing together on a spherical cortex of a small brain or “fold” of large ones might become useful, to others in this forum.