Application de la rotation de quaternions à une série chronologique vectorielle

Nov 24 2020

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 v1si v1contient 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 thetaet 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

1 henneray Nov 26 2020 at 11:52

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.

1 JamesTursa Nov 25 2020 at 01:21

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.