Here is an elegant method using quaternions that are blazingly fast; I can calculate 10 million rotations per second with appropriately vectorised numpy arrays. It relies on the quaternion extension to numpy found here.
Quaternion Theory:
A quaternion is a number with one real and 3 imaginary dimensions usually written as q = w + xi + yj + zk
where 'i', 'j', 'k' are imaginary dimensions. Just as a unit complex number 'c' can represent all 2d rotations by c=exp(i * theta)
, a unit quaternion 'q' can represent all 3d rotations by q=exp(p)
, where 'p' is a pure imaginary quaternion set by your axis and angle.
We start by converting your axis and angle to a quaternion whose imaginary dimensions are given by your axis of rotation, and whose magnitude is given by half the angle of rotation in radians. The 4 element vectors (w, x, y, z)
are constructed as follows:
import numpy as np
import quaternion as quat
v = [3,5,0]
axis = [4,4,1]
theta = 1.2 #radian
vector = np.array([0.] + v)
rot_axis = np.array([0.] + axis)
axis_angle = (theta*0.5) * rot_axis/np.linalg.norm(rot_axis)
First, a numpy array of 4 elements is constructed with the real component w=0 for both the vector to be rotated vector
and the rotation axis rot_axis
. The axis angle representation is then constructed by normalizing then multiplying by half the desired angle theta
. See here for why half the angle is required.
Now create the quaternions v
and qlog
using the library, and get the unit rotation quaternion q
by taking the exponential.
vec = quat.quaternion(*v)
qlog = quat.quaternion(*axis_angle)
q = np.exp(qlog)
Finally, the rotation of the vector is calculated by the following operation.
v_prime = q * vec * np.conjugate(q)
print(v_prime) # quaternion(0.0, 2.7491163, 4.7718093, 1.9162971)
Now just discard the real element and you have your rotated vector!
v_prime_vec = v_prime.imag # [2.74911638 4.77180932 1.91629719] as a numpy array
Note that this method is particularly efficient if you have to rotate a vector through many sequential rotations, as the quaternion product can just be calculated as q = q1 * q2 * q3 * q4 * ... * qn and then the vector is only rotated by 'q' at the very end using v' = q * v * conj(q).
This method gives you a seamless transformation between axis angle <---> 3d rotation operator simply by exp
and log
functions (yes log(q)
just returns the axis-angle representation!). For further clarification of how quaternion multiplication etc. work, see here