四元数与旋转--计算、推导与实现
发布日期:2021-05-09 07:42:04 浏览次数:15 分类:博客文章

本文共 4773 字,大约阅读时间需要 15 分钟。

转:

上一篇说到四元数的旋转表示,即从轴角法表示变化成的单位四元数的形式,即:

[p=(cosfrac{theta}{2}, boldsymbol{n}sinfrac{theta}{2})]

但是并没有说它是怎么来的。区区一个结果怎么能行呢,所以这里我就来细细说一说四元数的运算法则。

为了方便表示,这里规定两个四元数作为我们接下来演示用的小白鼠:

[q_1 = (s, boldsymbol{u}) = a_1+b_1i+c_1j+d_1k]

[q_2 = (t, boldsymbol{v}) = a_2+b_2i+c_2j+d_2k]

四元数基本运算

加减

加减法很简单,没什么好说的,只要把对应位置的元素加减就行:

[q_1pm q_2=(a_1 pm a_2)+(b_1 pm b_2)i+(c_1 pm c_2)j+(d_1 pm d_2)k]

写成矢量、标量部分分开的表示形式也可以:

[q_1 pm q_2=(s pm t, boldsymbol{u} pm boldsymbol{v})]

乘法

之前给出过有关 (i, j, k) 的运算律:

  1. (i^2=j^2=k^2=-1);
  2. (ij=k, jk=i, ki=j);
  3. (ij=-ji, jk=-kj, ki=-ik).

根据这些,我们把两个四元数当做多项式来计算,就可以得出结果。由于结果太长了,我们换个写法:假定 (q=q_1q_2=[w, x, y, z]), 那么

[w = a_1a_2-b_1b_2-c_1c_2-d_1d_2]

[x = a_1b_2+a_2b_1+c_1d_2-c_2d_1]

[y = a_1c_2+a_2c_1+d_1b_2-d_2b_1]

[z = a_1d_2+a_2d_1+b_1c_2-b_2c_1]

好家伙,这是真的长。有没有更简洁的表达方式呢?答案是有的。把四元数写成矢量、标量部分之后,我们会发现,计算结果的标量部分,也就是上面的 (w),可以写成:

[w=st-boldsymbol{u}cdotboldsymbol{v}]

而矢量部分的计算,也就是上面的 (x, y, z),每一个式子都有四项,这四项可以按照是否含有标量部分因子被分为两部分:含有标量因子的组合起来,我们看到它是 (sboldsymbol{v}+tboldsymbol{u}). 而不含有标量因子的部分,看上去像是行列式的形式,这里直接说结论吧,这部分可以写成 (boldsymbol{utimes v}). 于是,计算结果可以表示成

[q_1q_2=(st-boldsymbol{u}cdotboldsymbol{v}, sboldsymbol{v}+tboldsymbol{u}+boldsymbol{u}timesboldsymbol{v})]

这时,我们请出老朋友二元复数,我们可以看到很多相似点,四元数的乘积比普通的复数多了一项:(boldsymbol{u}timesboldsymbol{v}). 这一项决定了四元数的乘法没有交换律,甚至不具有反交换律。有趣的是,如果矢量部分是同向或反向的,这一项就是 0,而对于一般的复数,虚部只有一个单位 (i),它们只能共线。

按照直觉,四元数 (q) 的模 (||q||),应该就是把四个成员取平方加起来再开方,也就是常说的二范数。实际上也是如此。这部分就这样,公式就不写了。是不是很无聊?这可能是本文最无聊的一个小节了 2333。

共轭

四元数有三个虚部,怎么求共轭呢?答案很简单,就是把三个虚部都反过来。举个例子,(q=(s, boldsymbol{u})),那么 (q) 的共轭为

[q^star=(s, -boldsymbol{u})]

这个定义符合在二元复数里我们对共轭的一般认知,也是就是,满足共轭的运算律。比如说:

[qq^star=q^star q=||q||^2]

[q+q^star=2s, q-q^star=2boldsymbol{u}]

具体证明按照上面的乘法运算律算一算就知道了,不算复杂。

旋转的四元数表示究竟是怎么来的

旋转的分解

为了搞清楚这个问题,我们首先要将旋转操作相关的量用四元数表示出来。其实只有三个相关的向量:旋转前向量 (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}timesboldsymbol{u}_perp) 这两个垂直方向上的分量,我们是知道的,我们有:

[boldsymbol{v}_perp =boldsymbol{u}_perp costheta + boldsymbol{n}timesboldsymbol{u}_perpsintheta]

对比一下我们之前说过的四元数乘积形式,不难发现,如果我们用 (q=(costheta, boldsymbol{n}sintheta)) 来左乘 (u_{perp}=(0, boldsymbol{u}_perp)),会得到:

[qu_{perp}=(0-boldsymbol{n}sinthetacdotboldsymbol{u}_perp, boldsymbol{u}_perpcostheta+boldsymbol{n}timesboldsymbol{u}_perpsintheta)=(0, boldsymbol{v}_perp)]

代数小 trick

我们已经很接近最终结果了,这个结论和最终结论的差异就在于,被乘的是整个 (q) 而不是它的分量。要做的就是找到一种运算,可以不影响平行分量,只作用于垂直分量。看看在本文开始提出的最终结果,如果把它取平方的话,会发现

[p^2=(costheta, boldsymbol{n}sintheta)]

这就是我们说的作用于垂直分量的四元数 (q).

再进一步,我们把旋转写成一个很好看的形式:

[v=u_{parallel}+qu_{perp}=pp^{-1}u_{parallel}+ppu_{perp}]

因为 (p) 是单位四元数,所以 (p^{-1}=p^star) (不理解的话可以参考二元复数的性质)注意 (p)矢量部分,是和转轴平行的,那么把下面的引理一说,大家就明白我要做什么了:

引理:

假设 (u=(0, boldsymbol{u})) 是一个纯四元数,而 (q=(alpha, betaboldsymbol{n})),其中 (boldsymbol{n}) 是一个单位向量,(alpha, beta in R), 那么,根据 (boldsymbol{n})(boldsymbol{u}) 之间的关系,有以下结论:

  1. (boldsymbol{n} parallel boldsymbol{u}), 则 (qu = uq);
  2. (boldsymbol{n} perp boldsymbol{u}),则 (qu=uq^star).

证明同样省略,直接计算就得到了。

有没有发现,旋转后的两项刚好分别符合引理中的两个结论。于是:

[array{v&=&pp^star u_parallel+ppu_perp;;\&=&pu_parallel p^star+pu_perp p^star\&=&p(u_parallel+u_perp)p^star;;\&=&pup^starqquadqquad;}]

这就是我们的结论啦。

四元数运算的 Python 实现

四元数这个东西说年轻也不算年轻,那么是不是有相关的库呢?我在 GitHub 上面翻了翻,还真有,而且已经有 300+ stars 了。

> 链接在这 <

具体来说,这个库是为 numpy 提供了一个四元数的 dtype,而对于简单的单个四元数的运算也是完全可以胜任的。注意,如果你对运行速度有要求(或者单纯有强迫症,不想每次运行都看到 warning),最好安装一下 numba,这是一个提高 Python 运行速度的工具。

这个库的安装很简单,因为它已经加入 PyPI 了。注意这个包是基于 numpy 的所以要先安装 numpy.

pip install numpy-quaternion

有了这个工具,我们就可以轻松地应对四元数的运算了:

import numpy as npimport quaternionq1 = 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 的四元数与三维旋转一文的思路文章讲解十分详实,推荐阅读。

就到这里啦,有什么问题或发现什么错误可以留言或者私信呀。

祝福每一个在努力学习的人 : )

转:

上一篇:PXE配合Kickstart无人值守——批量装机简单如喝水(详细)
下一篇:Docker学习笔记

发表评论

最新留言

表示我来过!
[***.240.166.169]2025年04月21日 04时38分50秒

关于作者

    喝酒易醉,品茶养心,人生如梦,品茶悟道,何以解忧?唯有杜康!
-- 愿君每日到此一游!

推荐文章

bat 命令返回结果_【批处理】带你入门命令行 2023-01-24
c++ string取子串_Integer与String的设计哲学 2023-01-24
c++ 数组批量赋值_数组之间不能赋值?穿个马甲吧! 2023-01-24
cad模糊查询符号_mysql 正则模式和like模糊查询 2023-01-24
continue可以用if判断里面吗_谁能说说if()else()里的continue是干嘛的? 2023-01-24
ctrl c 和 ctrl v 不能用了_神奇操作,原来CTRL键还能这么用 2023-01-24
cytoscape安装java_Cytoscape史上最全攻略 2023-01-24
c语言程序设计年历显示,C语言程序设计报告《万年历》.doc 2023-01-24
C语言程序设计梁海英答案,1.5 习题 2023-01-24
c语言编写单片机中断,C语言AVR单片机中断程序写法 2023-01-24
#pragma region、{} 2023-01-24
ddr2的上电顺序_S5PV210 DDR2初始化 28个步骤总结 2023-01-24
deque stack java_「集合系列」- 初探 java 集合框架图 2023-01-24
easyexcel 导出 代码翻译converter_【starter推荐】简单高效Excel 导出工具 2023-01-24
echarts 如何在一条柱形显示两个数字_干货 | 如何快速制作数据地图?让你的可视化逼格再高一级!... 2023-01-24
eclipse设置utf8编码_记住没:永远不要在 MySQL 中使用 UTF8 2023-01-24
eclipse里source的快捷方法_Eclipse快捷键/快捷操作汇总 2023-01-24
elasticsearch 查询_Elasticsearch地理信息存储及查询之Geo_Point 2023-01-24
embedding层_【预估排序】Embedding+MLP: 深度学习预估排序通用框架(一) 2023-01-24
excel中最常用的30个函数_Excel玩转数据分析常用的43个函数! 2023-01-24