四元数旋转坐标系

Question

我们有无数的空间坐标（x，y和z）表示三维空间中的原子，我正在构造一个函数，将这些点转换成一个新的坐标系。将坐标转换为任意原点很简单，但我无法在下一步：三维点旋转计算中全神贯注。换句话说，我试图将点从（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

这个问题和@senderle给出的答案真的帮助了我的一个项目。答案是最小的，涵盖了一个人可能需要执行的大多数四元数计算的核心。

对于我自己的项目，我发现为所有操作分别使用函数并在每次需要时逐个导入它们是很乏味的，所以我实现了一个面向对象的版本。

四元数.py：

import numpy as np
from math import sin, cos, acos, sqrt

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 np.array(v)

class Quaternion:

    def from_axisangle(theta, v):
        theta = theta
        v = normalize(v)

        new_quaternion = Quaternion()
        new_quaternion._axisangle_to_q(theta, v)
        return new_quaternion

    def from_value(value):
        new_quaternion = Quaternion()
        new_quaternion._val = value
        return new_quaternion

    def _axisangle_to_q(self, theta, v):
        x = v[0]
        y = v[1]
        z = v[2]

        w = cos(theta/2.)
        x = x * sin(theta/2.)
        y = y * sin(theta/2.)
        z = z * sin(theta/2.)

        self._val = np.array([w, x, y, z])

    def __mul__(self, b):

        if isinstance(b, Quaternion):
            return self._multiply_with_quaternion(b)
        elif isinstance(b, (list, tuple, np.ndarray)):
            if len(b) != 3:
                raise Exception(f"Input vector has invalid length {len(b)}")
            return self._multiply_with_vector(b)
        else:
            raise Exception(f"Multiplication with unknown type {type(b)}")

    def _multiply_with_quaternion(self, q2):
        w1, x1, y1, z1 = self._val
        w2, x2, y2, z2 = q2._val
        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

        result = Quaternion.from_value(np.array((w, x, y, z)))
        return result

    def _multiply_with_vector(self, v):
        q2 = Quaternion.from_value(np.append((0.0), v))
        return (self * q2 * self.get_conjugate())._val[1:]

    def get_conjugate(self):
        w, x, y, z = self._val
        result = Quaternion.from_value(np.array((w, -x, -y, -z)))
        return result

    def __repr__(self):
        theta, v = self.get_axisangle()
        return f"((%.6f; %.6f, %.6f, %.6f))"%(theta, v[0], v[1], v[2])

    def get_axisangle(self):
        w, v = self._val[0], self._val[1:]
        theta = acos(w) * 2.0

        return theta, normalize(v)

    def tolist(self):
        return self._val.tolist()

    def vector_norm(self):
        w, v = self.get_axisangle()
        return np.linalg.norm(v)

在这个版本中，我们可以使用四元数向量乘法和四元数向量乘法的重载运算符

from quaternion import Quaternion
import numpy as np

x_axis_unit = (1, 0, 0)
y_axis_unit = (0, 1, 0)
z_axis_unit = (0, 0, 1)

r1 = Quaternion.from_axisangle(np.pi / 2, x_axis_unit)
r2 = Quaternion.from_axisangle(np.pi / 2, y_axis_unit)
r3 = Quaternion.from_axisangle(np.pi / 2, z_axis_unit)

# Quaternion - vector multiplication
v = r1 * y_axis_unit
v = r2 * v
v = r3 * v

print(v)

# Quaternion - quaternion multiplication
r_total = r3 * r2 * r1
v = r_total * y_axis_unit

print(v)

我不打算实现一个完整的四元数模块，所以这也是出于教学目的，就像@senderle的伟大回答一样。我希望这对那些想了解和尝试四元数新事物的人有帮助。

Answer 2

从代数的角度来看，用四元数来表示旋转并不困难。就我个人而言，我发现很难从视觉上解释四元数，但使用它们进行旋转的公式非常简单。我将在这里提供一组基本的引用函数。¹（另请参见hosolmaz给出的这个可爱的答案，他将这些函数打包在一起，创建了一个方便的四元数类。）

您可以将四元数（出于我们的目的）视为标量加上三维向量——抽象地说，w + xi + yj + zk，这里用一个简单的元组(w, x, y, z)表示。三维旋转的空间完全由四元数的子空间表示，即单位四元数的空间，因此要确保四元数是标准化的。您可以按照将任何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

请注意，为了简单起见，下面的函数假设四元数值已经规范化。实际上，您需要不时地对它们进行重新规范化，但处理这一问题的最佳方法将取决于问题域。这些函数只提供最基本的功能，仅供参考。

每个旋转都由一个单位四元数表示，旋转的串联对应于单位四元数的乘法。公式²如下：

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度旋转。（此处请记住右手法则；围绕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。四元数乘法公式看起来像老鼠窝，但推导起来很简单（如果繁琐的话）。先注意ii = jj = kk = -1；然后是ij = k，jk = i，ki = j；最后是ji = -k，kj = -i，ik = -j。然后将这两个四元数相乘，然后根据16次相乘的结果将这些项分布并重新排列。这也有助于说明为什么可以使用四元数来表示旋转；最后六个恒等式遵循右手规则，在从i到j的旋转和围绕k的旋转之间创建双射，以此类推。}

Answer 3

注意，矩阵的求逆一点也不简单！首先，所有n点（其中n是空间的维数）必须处于一般位置（即，没有单个点可以表示为其余点的线性组合[注意：这似乎确实是一个简单的要求，但在数值线性代数领域，它是非平凡的；当这种配置真的存在或不存在时，最终的决定将基于“实际领域”的特定知识）。

此外，新点和旧点的“对应关系”可能不准确（然后，您应该使用“真正对应关系”的最佳近似值，即：）。当lib提供伪逆（而不是试图使用纯逆）时，建议使用所有方法。

伪逆的优点是可以使用更多的点进行变换，因此增加了至少n个点处于一般位置的概率。

这里是一个例子，单位平方旋转90度，在二维逆时针方向（但显然这个测定在任何一个维度都有效），用numpy：

In []: P=  matrix([[0, 0, 1, 1],
                   [0, 1, 1, 0]])
In []: Pn= matrix([[0, -1, -1, 0],
                   [0,  0,  1, 1]])
In []: T= Pn* pinv(P)
In []: (T* P).round()
Out[]:
matrix([[ 0., -1., -1.,  0.],
        [ 0.,  0.,  1.,  1.]])

p.S.numpy也很快。在我那台普通的电脑上转换100万点：

In []: P= matrix(rand(2, 1e6))
In []: %timeit T* P
10 loops, best of 3: 37.7 ms per loop