将四元数旋转应用于向量时间序列答案

【问题标题】：Applying quaternion rotation to a vector time series将四元数旋转应用于向量时间序列
【发布时间】：2021-03-07 08:44:26
【问题描述】：

我在 Python numpy 数组中有一个 3D 向量的时间序列，类似于以下内容：

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

我想将每个向量围绕指定轴旋转指定角度theta。我一直在使用四元数来为一个向量实现这一点，正如在 henneray 的回答中找到的 here。

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

我的问题是，在 v1 中对每个向量应用相同旋转的最快方法是什么？

如果v1 包含二维向量数组，则无法从v1 创建四元数，因此我可以使用循环依次旋转每个数组元素；但是，在上面链接中的henneray的回答中，提到四元数可以应用于“适当矢量化的numpy数组”。有人对如何实施有任何建议吗？

（附带问题：如果我的theta 和axis 变量是长度与v1 相等的数组，是否也可以使用相同的方法通过相应的旋转来旋转v1 中的每个向量？）

【问题讨论】：

标签： python numpy vector vectorization quaternions

【解决方案1】：

需要先将[x,y,z]笛卡尔向量转换成第一个分量为零[0,x,y,z]的4-向量。然后您可以将其转换为四元数数组以进行矢量化计算。

下面的这个函数接受一个笛卡尔向量数组并围绕一个旋转轴旋转它们。您需要确保该轴的范数等于您的旋转角度 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:]

作为奖励，此函数采用旋转轴数组来将每个向量旋转相应的轴。

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

请注意，由于轴角和四元数表示之间的转换如此之多，与旋转矩阵代数相比，这几乎不会给您带来什么性能改进。只有当您通过多次连续旋转来旋转向量时，四元数才会真正受益，这样您就可以堆叠四元数乘法。

【讨论】：

【解决方案2】：

进行旋转计算本身的一种“快速”方法是将四元数转换为 3x3 方向余弦矩阵，将向量放在单个 3xN 连续矩阵中，然后调用 BLAS 库例程（例如 dgemm）做一个标准矩阵乘法。具有大 N 的良好 BLAS 库可以多线程执行此计算。

【讨论】：