python - 通过四元数旋转坐标系

Question

我们有无数个空间坐标(x、y 和 z)表示 3d 空间中的原子，我正在构建一个函数，将这些点转换为新的坐标系。将坐标移动到任意原点很简单，但我无法理解下一步:3d 点旋转计算。换句话说，我试图将点从 (x, y, z) 转换为 (x', y', z')，其中 x', y' 和 z' 是根据 i', j'和 k'，我在 euclid python module 的帮助下制作的新轴向量.

我认为我需要的只是一个欧几里得四元数来做到这一点，即

>>> q * Vector3(x, y, z)
Vector3(x', y', z')

但要做到这一点，我相信我需要一个旋转轴矢量和一个旋转角度。但我不知道如何从 i'、j' 和 k' 计算这些。这似乎是一个从头开始编写代码的简单过程，但我怀疑这样的事情需要线性代数来自己解决。非常感谢您对正确方向的插入。

Answer 1

从代数的角度来看，使用四元数来表示旋转并不困难。就我个人而言，我发现很难从视觉上推理四元数，但使用它们进行旋转所涉及的公式并不过分复杂。我将在此处提供一组基本的引用函数。1(另请参阅 hosolmaz 的这个可爱的答案，他将它们打包在一起以创建一个方便的 Quaternion 类。)
您可以将四元数(出于我们的目的)视为标量加上 3 维向量——抽象地说，w + xi + yj + zk ，这里用一个普通的元组 (w, x, y, z) 表示. 3-d 旋转空间完全由四元数的子空间表示，单位四元数的空间，所以你要确保你的四元数是归一化的。您可以按照对任何 4 向量进行归一化的方式来执行此操作(即幅度应接近 1；如果不是，则按幅度缩小值):

def normalize(v, tolerance=0.00001):
    mag2 = sum(n * n for n in v)
    if abs(mag2 - 1.0) > tolerance:
        mag = sqrt(mag2)
        v = tuple(n / mag for n in v)
    return v

请注意，为简单起见，以下函数假设四元数值已经标准化。在实践中，您需要不时地重新规范化它们，但处理这种情况的最佳方法将取决于问题域。这些功能仅提供非常基础的功能，仅供引用。
每个旋转都由一个单位四元数表示，旋转的串联对应于单位四元数的乘法。其公式2如下:

def q_mult(q1, q2):
    w1, x1, y1, z1 = q1
    w2, x2, y2, z2 = q2
    w = w1 * w2 - x1 * x2 - y1 * y2 - z1 * z2
    x = w1 * x2 + x1 * w2 + y1 * z2 - z1 * y2
    y = w1 * y2 + y1 * w2 + z1 * x2 - x1 * z2
    z = w1 * z2 + z1 * w2 + x1 * y2 - y1 * x2
    return w, x, y, z

要将向量旋转四元数，您还需要四元数的共轭:

def q_conjugate(q):
    w, x, y, z = q
    return (w, -x, -y, -z)

四元数向量乘法就是将向量转换为四元数(通过设置 w = 0 并保留 x 、 y 和 z 相同)然后乘以 q * v * q_conjugate(q) :

def qv_mult(q1, v1):
    q2 = (0.0,) + v1
    return q_mult(q_mult(q1, q2), q_conjugate(q1))[1:]
    

最后，您需要知道如何从轴角旋转转换为四元数并返回。这也出人意料地直截了当。通过调用 normalize 来“清理”输入和输出是有意义的。 .

def axisangle_to_q(v, theta):
    v = normalize(v)
    x, y, z = v
    theta /= 2
    w = cos(theta)
    x = x * sin(theta)
    y = y * sin(theta)
    z = z * sin(theta)
    return w, x, y, z

然后回来:

def q_to_axisangle(q):
    w, v = q[0], q[1:]
    theta = acos(w) * 2.0
    return normalize(v), theta

这是一个快速使用示例。围绕 x、y 和 z 轴旋转 90 度的序列会将 y 轴上的向量返回到其原始位置。此代码执行这些轮换:

x_axis_unit = (1, 0, 0)
y_axis_unit = (0, 1, 0)
z_axis_unit = (0, 0, 1)
r1 = axisangle_to_q(x_axis_unit, numpy.pi / 2)
r2 = axisangle_to_q(y_axis_unit, numpy.pi / 2)
r3 = axisangle_to_q(z_axis_unit, numpy.pi / 2)

v = qv_mult(r1, y_axis_unit)
v = qv_mult(r2, v)
v = qv_mult(r3, v)

print v
# output: (0.0, 1.0, 2.220446049250313e-16)

请记住，此旋转序列不会将所有向量返回到相同位置；例如，对于 x 轴上的向量，它将对应于绕 y 轴旋转 90 度。 (想想这里的 right-hand-rule；绕 y 轴的正旋转会将 x 轴上的向量插入负 z 区域。)

v = qv_mult(r1, x_axis_unit)
v = qv_mult(r2, v)
v = qv_mult(r3, v)

print v
# output: (4.930380657631324e-32, 2.220446049250313e-16, -1.0)

与往常一样，如果您在这里发现任何问题，请告诉我。

1. 这些改编自 OpenGL 教程 archived here .
2.四元数乘法公式乍一看像个可怕的老鼠窝，但推导很容易，虽然很乏味。使用铅笔和纸，您可以像这样表示两个四元数:w + xi + yj + zk .然后注意ii = jj = kk = -1 ;那个ij = k , jk = i , ki = j ;还有那个ji = -k , kj = -i , ik = -j .最后，将两个四元数相乘，分配项并根据 16 次乘法中的每一个的结果重新排列它们。这有助于说明为什么可以使用四元数来表示旋转；最后六个恒等式遵循右手定则，在来自 i 的旋转之间创建双射至 j和旋转 k ， 等等。
如果你这样做，你会看到身份 ii = jj = kk = -1解释 w 中的最后三个术语公式;那个jk = i解释 x 中的第三项公式;还有那个kj = -i解释了 x 中的第四项公式。 y和 z公式的工作方式相同。其余项都是普通标量乘法的例子。