四元数与旋转--计算、推导与实现
上一篇说到四元数的旋转表示,即从轴角法表示变化成的单位四元数的形式,即:
\]
但是并没有说它是怎么来的。区区一个结果怎么能行呢,所以这里我就来细细说一说四元数的运算法则。
为了方便表示,这里规定两个四元数作为我们接下来演示用的小白鼠:
\]
\]
四元数基本运算
加减
加减法很简单,没什么好说的,只要把对应位置的元素加减就行:
\]
写成矢量、标量部分分开的表示形式也可以:
\]
乘法
之前给出过有关 \(i, j, k\) 的运算律:
- \(i^2=j^2=k^2=-1\);
- \(ij=k, jk=i, ki=j\);
- \(ij=-ji, jk=-kj, ki=-ik\).
根据这些,我们把两个四元数当做多项式来计算,就可以得出结果。由于结果太长了,我们换个写法:假定 \(q=q_1q_2=[w, x, y, z]\), 那么
\]
\]
\]
\]
好家伙,这是真的长。有没有更简洁的表达方式呢?答案是有的。把四元数写成矢量、标量部分之后,我们会发现,计算结果的标量部分,也就是上面的 \(w\),可以写成:
\]
而矢量部分的计算,也就是上面的 \(x, y, z\),每一个式子都有四项,这四项可以按照是否含有标量部分因子被分为两部分:含有标量因子的组合起来,我们看到它是 \(s\boldsymbol{v}+t\boldsymbol{u}\). 而不含有标量因子的部分,看上去像是行列式的形式,这里直接说结论吧,这部分可以写成 \(\boldsymbol{u\times v}\). 于是,计算结果可以表示成
\]
这时,我们请出老朋友二元复数,我们可以看到很多相似点,四元数的乘积比普通的复数多了一项:\(\boldsymbol{u}\times\boldsymbol{v}\). 这一项决定了四元数的乘法没有交换律,甚至不具有反交换律。有趣的是,如果矢量部分是同向或反向的,这一项就是 0,而对于一般的复数,虚部只有一个单位 \(i\),它们只能共线。
模
按照直觉,四元数 \(q\) 的模 \(||q||\),应该就是把四个成员取平方加起来再开方,也就是常说的二范数。实际上也是如此。这部分就这样,公式就不写了。是不是很无聊?这可能是本文最无聊的一个小节了 2333。
共轭
四元数有三个虚部,怎么求共轭呢?答案很简单,就是把三个虚部都反过来。举个例子,\(q=(s, \boldsymbol{u})\),那么 \(q\) 的共轭为
\]
这个定义符合在二元复数里我们对共轭的一般认知,也是就是,满足共轭的运算律。比如说:
\]
\]
具体证明按照上面的乘法运算律算一算就知道了,不算复杂。
旋转的四元数表示究竟是怎么来的
旋转的分解
为了搞清楚这个问题,我们首先要将旋转操作相关的量用四元数表示出来。其实只有三个相关的向量:旋转前向量 \(\boldsymbol{u}\), 旋转后向量 \(\boldsymbol{v}\), 和旋转轴 \(\boldsymbol{n}\) 和一个角度,即旋转角 \(\theta\). 旋转前后的向量可以用纯四元数表示,即转前:\(u=(0, \boldsymbol{u})\),转后:\(v=(0, \boldsymbol{v})\). 注意,这里的转轴的表示向量 \(\boldsymbol{n}\) 是一个单位向量。
要解决这种这种绕固定轴旋转的问题,一个自然的思路是将被旋转的向量分解成平行于轴向的 \(\boldsymbol{u}_{\parallel}\) 和垂直于轴向的 \(\boldsymbol{u}_\perp\).
显然,旋转对于平行分量是无作用的,只作用于垂直分量。暂且不考虑与转轴平行的分量,我们将目光放在垂直分量上。把 \(\boldsymbol{u}_\perp\) 旋转 \(\theta\) 角度,得到的结果就是 \(\boldsymbol{v}_\perp\),所以后者相对前者的分量,也就是在 \(\boldsymbol{u}_\perp\) 和 \(\boldsymbol{n}\times\boldsymbol{u}_\perp\) 这两个垂直方向上的分量,我们是知道的,我们有:
\]
对比一下我们之前说过的四元数乘积形式,不难发现,如果我们用 \(q=(\cos\theta, \boldsymbol{n}\sin\theta)\) 来左乘 \(u_{\perp}=(0, \boldsymbol{u}_\perp)\),会得到:
\]
代数小 trick
我们已经很接近最终结果了,这个结论和最终结论的差异就在于,被乘的是整个 \(q\) 而不是它的分量。要做的就是找到一种运算,可以不影响平行分量,只作用于垂直分量。看看在本文开始提出的最终结果,如果把它取平方的话,会发现
\]
这就是我们说的作用于垂直分量的四元数 \(q\).
再进一步,我们把旋转写成一个很好看的形式:
\]
因为 \(p\) 是单位四元数,所以 \(p^{-1}=p^\star\) (不理解的话可以参考二元复数的性质)注意 \(p\) 的矢量部分,是和转轴平行的,那么把下面的引理一说,大家就明白我要做什么了:
引理:
假设 \(u=(0, \boldsymbol{u})\) 是一个纯四元数,而 \(q=(\alpha, \beta\boldsymbol{n})\),其中 \(\boldsymbol{n}\) 是一个单位向量,\(\alpha, \beta \in R\), 那么,根据 \(\boldsymbol{n}\) 和 \(\boldsymbol{u}\) 之间的关系,有以下结论:
- 若 \(\boldsymbol{n} \parallel \boldsymbol{u}\), 则 \(qu = uq\);
- 若 \(\boldsymbol{n} \perp \boldsymbol{u}\),则 \(qu=uq^\star\).
证明同样省略,直接计算就得到了。
有没有发现,旋转后的两项刚好分别符合引理中的两个结论。于是:
&=&pu_\parallel p^\star+pu_\perp p^\star\\
&=&p(u_\parallel+u_\perp)p^\star\;\;\\
&=&pup^\star\qquad\qquad\;}
\]
这就是我们的结论啦。
四元数运算的 Python 实现
四元数这个东西说年轻也不算年轻,那么是不是有相关的库呢?我在 GitHub 上面翻了翻,还真有,而且已经有 300+ stars 了。
具体来说,这个库是为 numpy
提供了一个四元数的 dtype
,而对于简单的单个四元数的运算也是完全可以胜任的。注意,如果你对运行速度有要求(或者单纯有强迫症,不想每次运行都看到 warning),最好安装一下 numba
,这是一个提高 Python 运行速度的工具。
这个库的安装很简单,因为它已经加入 PyPI 了。注意这个包是基于 numpy
的所以要先安装 numpy
.
pip install numpy-quaternion
有了这个工具,我们就可以轻松地应对四元数的运算了:
import numpy as np
import quaternion
q1 = np.quaternion(1, 2, 3, 4)
q2 = np.quaternion(5, 6, 7, 8)
print("q1 = ", q1)
print("\nq2 = ", q2)
print("\nq1 + q2 = ", q1 + q2) # quaternion(6, 8, 10, 12)
print("\nq2 - q1 = ", q2 - q1) # quaternion(4, 4, 4, 4)
print("\nq2 * q1 = ", q1 * q2) # quaternion(-60, 12, 30, 24)
print("\nThe conjugate of q1 is ", q1.conjugate()) # quaternion(1, -2, -3, -4)
print("\nThe square of q1's norm is ", q1.norm()) # 30.0 这里输出的实际上并不是模,而是模的平方,算是 bug 吧:-/
arr = np.array([[q1, q2]])
print("\nThere is an example of array of quaternions:\n", arr)
print("\nIts transpose times itself is:\n", arr.T * arr)
print("\nIts elementwise exponential in base e is:\n", np.exp(arr))
注意几个点:
- 现行版本的内置函数
norm
是有问题的,会返回模的平方,如果需要求模记得取平方根; -
numpy
的一些逐元素执行的运算,如上面例子中的numpy.exp
在这个库中也是实现了的。 - 这里我没有讲四元数的指数运算,不过其实很简单 w 感兴趣的可以自行了解或者在评论区留言。
- 四元数乘法没有交换律,这一点在矩阵运算的时候也有体现。
结尾
例行总结好无聊啊,不想写,那就写一下本文的参考吧。
本文的求解思路参考了 Krasjet 的四元数与三维旋转一文的思路文章讲解十分详实,推荐阅读。
就到这里啦,有什么问题或发现什么错误可以留言或者私信呀。
祝福每一个在努力学习的人 : )