我认为我只需要一个欧氏四元数就可以完成这个工作,即:
>>> q * Vector3(x, y, z)
Vector3(x', y', z')
但是为了做到这一点,我相信我需要一个旋转轴向量和一个旋转角度。但我不知道如何从i'、j'和k'计算出它们。这似乎是一个可以从头开始编写的简单过程,但我怀疑这需要线性代数才能自己解决。非常感谢指导。
我认为我只需要一个欧氏四元数就可以完成这个工作,即:
>>> q * Vector3(x, y, z)
Vector3(x', y', z')
但是为了做到这一点,我相信我需要一个旋转轴向量和一个旋转角度。但我不知道如何从i'、j'和k'计算出它们。这似乎是一个可以从头开始编写的简单过程,但我怀疑这需要线性代数才能自己解决。非常感谢指导。
使用四元数来表示旋转在代数角度上并不困难。就我个人而言,我发现很难从可视化的角度来理解四元数,但是使用它们进行旋转所涉及的公式并不过于复杂。在此,我提供了一组基本的参考函数1。(请参见hosolmaz的精彩答案,在其中,他将这些打包在一起创建了一个方便的四元数类。)
你可以将四元数(对我们而言)看作一个标量加上一个三维向量——抽象地说,w + xi + yj + zk
,这里用一个普通的元组(w, x, y, z)
来表示。 3D旋转空间在四元数中由一部分子空间完全表示,即单位四元数空间,因此确保四元数已经被标准化是很重要的。您可以按照任何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度。(想象一下右手法则;关于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. 四元数乘法公式起初看起来像是一个可怕的老鼠窝,但推导很容易,尽管繁琐。用铅笔和纸,你可以将两个四元数表示为: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
的公式也是同样的原理。其余的项都是普通标量乘法的例子。
qv_mult
有用,它应该返回一个向量,而不是四元数,因此我删除了第一个分量(它本来就是零!)。 - Hooked这个问题和@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的绝妙答案一样。我希望这能帮助那些想要理解和尝试使用四元数进行新事物的人。
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
也很快。在我的普通电脑上,1百万个点的转换速度非常快:
In []: P= matrix(rand(2, 1e6))
In []: %timeit T* P
10 loops, best of 3: 37.7 ms per loop
(0, 0, 1)
、(0, 1, 0)
和(1, 0, 0)
应该翻译成什么? - Anon.