通过四元数旋转坐标系

Question

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

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

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

但是为了使我相信我需要一个旋转轴向量和一个旋转角度.但我不知道如何从i',j'和k'计算这些.这似乎是一个从头开始编码的简单程序,但我怀疑这样的东西需要线性代数来自行解决.非常感谢你在正确的方向上轻推.

Answer 1

从代数的角度来看,使用四元数来表示旋转并不困难.就个人而言,我发现很难在视觉上对四元数进行推理,但使用它们进行旋转所涉及的公式非常简单.我将在这里提供一组基本的参考函数.¹(请参阅hosolmaz的这个可爱的答案,他将这些包装在一起以创建一个方便的Quaternion类.)

您可以将四元数(出于我们的目的)视为标量加上三维向量 - 抽象地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教程.}

^{2.四元数乘法公式看起来像一个疯狂的老鼠的巢,但推导很简单(如果单调乏味).先注意一下ii = jj = kk = -1; 那么ij = k,jk = i,ki = j; 最后是ji = -k,kj = -i,ik = -j.然后乘以两个四元数,分配出术语并根据16次乘法中的每一次的结果重新排列它们.这也有助于说明为什么可以使用四元数来表示旋转; 过去六年身份遵循右手法则,创造旋转之间的双射来自 i于j和旋转左右 k,依此类推.}

Answer 2

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

对于我自己的项目，我发现为所有操作都具有单独的功能并在每次需要时逐一导入它们很麻烦，因此我实现了一个面向对象的版本。

quaternion.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的一个很好的答案一样。我希望这对希望了解并尝试使用四元数的新事物有所帮助。