Application de la rotation de quaternions à une série chronologique vectorielle
J'ai une série chronologique de vecteurs 3D dans un tableau numpy Python similaire à ce qui suit:
array([[-0.062, -0.024, 1. ],
[-0.071, -0.03 , 0.98 ],
[-0.08 , -0.035, 0.991],
[-0.083, -0.035, 0.98 ],
[-0.083, -0.035, 0.977],
[-0.082, -0.035, 0.993],
[-0.08 , -0.034, 1.006],
[-0.081, -0.032, 1.008],
.......
Je souhaite faire pivoter chaque vecteur autour d'un axe spécifié selon un angle spécifié theta
. J'ai utilisé des quaternions pour atteindre cet objectif pour un vecteur comme trouvé ici dans la réponse de Henneray.
v1 = np.array ([1, -2, 0])
axis = np.array([-4, -2, 3])
theta = 1.5
rot_axis = np.insert(axis, 0, 0, axis=0)
axis_angle = (theta*0.5) * rot_axis/np.linalg.norm(rot_axis)
vec = quat.quaternion(*v1)
qlog = quat.quaternion(*axis_angle)
q = np.exp(qlog)
v_prime = q * vec * np.conjugate(q)
v_prime_vec = v_prime.imag
Ma question est la suivante: quel est le moyen le plus rapide d'appliquer la même rotation à chaque vecteur dans la v1?
Vous ne pouvez pas créer un quaternion à partir de v1
si v1
contient un tableau 2D de vecteurs, donc je pourrais utiliser une boucle pour faire pivoter chaque élément du tableau à son tour; cependant, dans la réponse de Henneray dans le lien ci-dessus, il est mentionné que les quaternions pourraient être appliqués à des «tableaux numpy correctement vectorisés». Quelqu'un a-t-il des suggestions sur la façon dont cela pourrait être mis en œuvre?
(Une question secondaire: si mes variables theta
et axis
étaient des tableaux de longueur égale à v1, la même méthode pourrait-elle également être utilisée pour faire pivoter chaque vecteur dans v1 par une rotation correspondante?)
Réponses
Il faut d'abord convertir les vecteurs cartésiens [x, y, z] en 4 vecteurs avec la première composante égale à zéro [0, x, y, z]. Ensuite, vous pouvez convertir ceci en un tableau quaternion pour effectuer des calculs vectorisés.
Cette fonction ci-dessous prend un tableau de vecteurs cartésiens et les fait pivoter autour d'un seul axe de rotation. Vous devrez vous assurer que la norme de cet axe est égale à votre angle de rotation thêta.
def rotate_vectors(vecs, axis):
"""
Rotate a list of 3D [x,y,z] vectors about corresponding 3D axis
[x,y,z] with norm equal to the rotation angle in radians
Parameters
----------
vectors : numpy.ndarray with shape [n,3]
list of [x,y,z] cartesian vector coordinates
axis : numpy.ndarray with shape [3]
[x,y,z] axis to rotate corresponding vectors about
"""
# Make an 4 x n array of zeros
vecs4 = np.zeros([vecs.shape[0],vecs.shape[1]+1])
# Fill the imaginary i, j, k components with x, y, z values, leaving the real part w=0
vecs4[:,1:] = vecs
# Convert to quaternion array
vecsq = quat.as_quat_array(vecs4)
# Make a rotation quaternion
qrot = quat.from_rotation_vector(axis)
# Rotate vectors
vecsq_rotated = qrot * vecsq * qrot.conjugate()
# Cast quaternion array to float and return only imaginary components (ignore real part)
return quat.as_float_array(vecsq_rotated)[:,1:]
En prime, cette fonction utilise un tableau d'axes de rotation pour faire pivoter chaque vecteur par les axes correspondants.
def rotate_vectors_each(vecs, axes):
"""
Rotate a list of 3D [x,y,z] vectors about corresponding 3D axes
[x,y,z] with norm equal to the rotation angle in radians
Parameters
----------
vectors : numpy.ndarray with shape [n,3]
list of [x,y,z] cartesian vector coordinates
axes : numpy.ndarray with shape [n,3]
axes to rotate corresponding vectors about
n = pulse shape time domain
3 = [x,y,z]
"""
# Make an 4 x n array of zeros
vecs4 = np.zeros([vecs.shape[0],vecs.shape[1]+1])
# Fill the imaginary i, j, k components with x, y, z values, leaving the real part w=0
vecs4[:,1:] = vecs
# Convert to quaternion array
vecsq = quat.as_quat_array(vecs4)
# Make an 4 x n array of zeros
rots4 = np.zeros([rots.shape[0],rots.shape[1]+1])
# Fill the imaginary i, j, k components with x, y, z values, leaving the real part w=0
rots4[:,1:] = rots
# Convert to quaternion array and take exponential
qrots = np.exp(quat.as_quat_array(0.5 * rots4))
# Rotate vectors
vecsq_rotated = qrots * vecsq * qrots.conjugate()
return quat.as_float_array(vecsq_rotated)[:,1:]
Notez qu'avec autant de conversions entre l'angle d'axe et la représentation quaternion, cela vous donnera peu d'amélioration des performances par rapport à l'algèbre matricielle de rotation. Les quaternions ne bénéficient vraiment que lorsque vous faites pivoter un vecteur sur de nombreuses rotations séquentielles, grâce auxquelles vous pouvez empiler la multiplication de quaternions.
Un moyen "rapide" de faire le calcul de rotation lui-même serait de transformer votre quaternion en une matrice cosinus de direction 3x3, d'avoir vos vecteurs dans une seule matrice contiguë 3xN, puis d'appeler une routine de bibliothèque BLAS (par exemple, dgemm) pour faire un standard la matrice se multiplie. Une bonne bibliothèque BLAS avec un grand N ferait ce calcul multi-thread.