Aplicação de rotação de quatérnio a uma série temporal de vetor

Nov 24 2020

Eu tenho uma série temporal de vetores 3D em uma matriz numpy Python semelhante à seguinte:

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

Quero girar cada vetor em torno de um eixo especificado por meio de um ângulo especificado theta. Tenho usado quatérnios para conseguir isso para um vetor, conforme encontrado aqui na resposta 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

Minha pergunta é: qual é a maneira mais rápida de aplicar a mesma rotação a cada vetor na v1?

Você não pode criar um quatérnion a partir de v1se v1contém uma matriz 2D de vetores, então eu poderia usar um loop para girar cada elemento da matriz por vez; entretanto, na resposta de Henneray no link acima, é mencionado que os quatérnios poderiam ser aplicados a 'matrizes numpy vetorizadas apropriadamente'. Alguém tem alguma sugestão sobre como isso poderia ser implementado?

(Uma questão lateral: se minhas variáveis thetae axisvariáveis ​​fossem matrizes de comprimento igual a v1, o mesmo método também poderia ser usado para girar cada vetor em v1 por meio de uma rotação correspondente?)

Respostas

1 henneray Nov 26 2020 at 11:52

É necessário primeiro converter os vetores cartesianos [x, y, z] em 4 vetores com o primeiro componente igual a zero [0, x, y, z]. Em seguida, você pode lançar isso para uma matriz de quatérnio para fazer cálculos vetorizados.

Esta função abaixo pega uma matriz de vetores cartesianos e os gira em torno de um único eixo de rotação. Você precisará certificar-se de que a norma desse eixo é igual ao seu ângulo de rotação teta.

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 um bônus, esta função usa uma matriz de eixos de rotação para girar cada vetor pelos eixos correspondentes.

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

Observe que, com tantas conversões entre o ângulo do eixo e a representação do quatérnio, isso proporcionará pouca melhoria no desempenho em relação à álgebra da matriz de rotação. Os quatérnios realmente só se beneficiam quando você gira um vetor por meio de muitas rotações sequenciais, por meio das quais você pode empilhar a multiplicação do quatérnio.

1 JamesTursa Nov 25 2020 at 01:21

Uma maneira "rápida" de fazer o cálculo de rotação em si seria transformar seu quatérnion em uma matriz cosseno de direção 3x3, ter seus vetores em uma única matriz contígua 3xN e, em seguida, chamar uma rotina de biblioteca BLAS (por exemplo, dgemm) para fazer um padrão multiplicação da matriz. Uma boa biblioteca BLAS com N grande faria esse cálculo multi-threaded.