Parallelize rotation of points around axis

[adriengoeller]

Sadly, scipy rotation cannot be used in that way: do in parallel several different rotations applied to several points groups vectors.

We propose here to implement a parallel rotation based on quaternions and numpy only to be pyodide-compatible.

Here is a code implementing the rotation and comparing values to scipy and another standard quaternion library: pyquaternion.

import numpy as np
from pyquaternion import Quaternion
from scipy.spatial.transform import Rotation


def ham(q1, q2):
    a1, b1, c1, d1 = q1
    a2, b2, c2, d2 = q2

    return np.array(
        [
            a1 * a2 - b1 * b2 - c1 * c2 - d1 * d2,
            a1 * b2 + b1 * a2 + c1 * d2 - d1 * c2,
            a1 * c2 - b1 * d2 + c1 * a2 + d1 * b2,
            a1 * d2 + b1 * c2 - c1 * b2 + d1 * a2,
        ]
    )


vector = np.array([-9.86411084, 0.10916063, -0.68953008])
purequat = np.array([0, -9.86411084, 0.10916063, -0.68953008])
# order: w, x, y, z
quat = np.array([-0.54312134, 0.42388916, -0.45617676, 0.5632019])
# order x, y,  z, w
quat_for_scipy = np.array([0.42388916, -0.45617676, 0.5632019, -0.54312134])
conj = np.array([1, -1, -1, -1])
quatconj = quat * conj  # hand conjugate
Q = Quaternion(quat)
R = Rotation.from_quat(quat_for_scipy)

print("manual:", ham(quat, ham(purequat, quatconj))[1:])
print("pyQuaternion:", Q.rotate(vector))
print("scipyRotation:", R.apply(vector))
#-----------------
manual: [-0.14691211  9.88691296 -0.08305227]
pyQuaternion: [-0.14691852  9.88734378 -0.08305589]
scipyRotation: [-0.14691852  9.88734378 -0.08305589]

This function can be vectorized in this way: (example for the trocky part: hamilton operator)

def quat_multiply(quaternion0, quaternion1):
    x0, y0, z0, w0 = np.split(quaternion0, 4, axis=-1)
    x1, y1, z1, w1 = np.split(quaternion1, 4, axis=-1)
    return np.concatenate(
        (x1*w0 + y1*z0 - z1*y0 + w1*x0,
         -x1*z0 + y1*w0 + z1*x0 + w1*y0,
         x1*y0 - y1*x0 + z1*w0 + w1*z0,
         -x1*x0 - y1*y0 - z1*z0 + w1*w0),
        axis=-1)

# define 
q1= np.random.rand(10,4)
q2 = np.random.rand(10,4)

display(q1,q2)

quat_multiply(q1,q2)

Goal: provide in geometry/tools.py of vectorized function allowing to do:

q_rotation_vectors = np.random.rand(10,4)
q_points_to_rotate = np.random.rand(10,4)

q_points_rotated = rotate(q_rotation_vectors , q_points_to_rotate)