Aplicar la rotación de cuaterniones a una serie de tiempo vectorial

Question

Tengo una serie temporal de vectores 3D en una matriz numérica de Python similar a la siguiente:

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],
       .......

Quiero rotar cada vector alrededor de un eje específico a través de un ángulo específico theta. He estado usando cuaterniones para lograr esto para un vector como se encuentra aquí en la respuesta 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

Mi pregunta es, ¿cuál es la forma más rápida de aplicar la misma rotación a cada vector en v1?

No puede crear un cuaternión a partir de v1si v1contiene una matriz 2D de vectores, por lo que podría usar un bucle para rotar cada elemento de la matriz a su vez; sin embargo, en la respuesta de Henneray en el enlace anterior, se menciona que los cuaterniones podrían aplicarse a 'matrices numpy adecuadamente vectorizadas'. ¿Alguien tiene alguna sugerencia sobre cómo se podría implementar?

(Una pregunta paralela: si my thetay axisvariables fueran matrices de igual longitud que v1, ¿podría usarse también el mismo método para rotar cada vector en v1 a través de una rotación correspondiente?)

Answer 1

Primero es necesario convertir los vectores cartesianos [x, y, z] en 4 vectores con el primer componente igual a cero [0, x, y, z]. Luego, puede convertir esto en una matriz de cuaterniones para hacer cálculos vectorizados.

Esta función a continuación toma una matriz de vectores cartesianos y los gira alrededor de un solo eje de rotación. Deberá asegurarse de que la norma de este eje sea igual a su ángulo de rotación theta.

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:]

Como beneficio adicional, esta función toma una matriz de ejes de rotación para rotar cada vector por los ejes correspondientes.

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:]

Tenga en cuenta que con tantas conversiones entre el ángulo del eje y la representación del cuaternión, esto le dará una pequeña mejora en el rendimiento con respecto al álgebra matricial de rotación. Los cuaterniones realmente solo se benefician cuando está rotando un vector a través de muchas rotaciones secuenciales, por lo que puede apilar la multiplicación de cuaterniones.

Answer 2

Una forma "rápida" de hacer el cálculo de rotación en sí sería convertir su cuaternión en una matriz de coseno de dirección 3x3, tener sus vectores en una única matriz contigua 3xN y luego llamar a una rutina de biblioteca BLAS (por ejemplo, dgemm) para hacer un estándar multiplicar matriz. Una buena biblioteca BLAS con N grande haría este cálculo multiproceso.