随笔- 2  文章- 0  评论- 1  阅读- 210 

四元数旋转:从一维到四维

Main

1、前言

本文讨论四元数与三维空间中的旋转关系,也就是四元数从四维空间对三维空间里的向量起了什么作用。众所周知,四元数存在于四维中,这也让它蒙上了一层神秘的色彩;很多文章把它描述得更加神秘,而有些文章为了使其“去神秘化”又写得太过简短,导致很多四元数的细节都被略过了。网上关于四元数的优秀的文献有很多,复制粘贴的也非常多,更有甚者连偏手性都没声明;本人在使用四元数的过程中也产生过一些疑惑,在得到解决后也想写一篇总结性的文章。

以上种种,成为本文的写作动机。三维空间中的旋转其实比想象中要复杂得多,比如二维旋转的次序无关紧要而三维的则不然(数学上的表达就是不满足交换律);再比如我们生活中习以为常的“前后左右”等概念是相对的和模糊的(对此,数学上的严格表达就是偏手性和多坐标系);还有像使用欧拉角导致的万向节锁定的问题存在(四元数则没有这个问题),等等。正因为三维旋转具有这些特性,本文的论述从一维空间开始,逐步升维,在当前维度中阐述对应的数学概念,这有利于更好地理解旋转的本质。由于四元数就是升维了的复数,通过对升维的阐述,把前后相关的知识联系起来。

本文用通俗易懂的方式论述,仅需最基本的线性代数知识即可阅读。但由于是应用于游戏、图形学方向,本文更强调的是四元数(以及其他维度中的数学概念)带来的几何意义;再加上笔者本人才疏学浅的事实,所以如果您是想从更高层次的抽象代数来阅读的话,那么本文并不适合您。本文的结构为:

  • 前面的部分,用一些基础的线性代数和几何学知识,来推导和理解四元数的乘积公式;
  • 后面的部分,用可视化的方式,看看四元数究竟是怎么旋转三维中的向量的。由于四元数存在于四维空间中,对三维的我们来说已经不可见,但是我们可以借用“球极投影”的方法,看看它在三维中的直观投影。有了前置的知识,更有利于这一部分的理解;
  • 最后的部分,我们让四元数在代码中实现。
  • 另外,为了加速计算过程,我将省略一些公式的推导或者证明过程。大多数时候即使不清楚它们的推导过程也不会有什么问题,但如果读者感兴趣,可以移步到扩展阅读部分阅读。

1.1 重要约定

1.1.1 偏手性

偏手性和左右手定则大多数时候只是数学上的惯用约定,实际上可以自由组合。由于unity大多数的坐标空间是左手坐标系,因此本文将采用:

  • 左手坐标系;
  • 左手定则。这意味着,从对着旋转轴的方向看,向量顺时针旋转;影响叉乘方向。

1.1.2 符号约定

  • 向量:使用小写粗体字母,如:v
  • 四元数:使用小写细体字母,如:p
  • 纯四元数:使用小写细体字母,如:v
  • 矩阵:使用大写粗体字母,如:M
  • 集合:使用黑板报粗体字母,如:H

2、一维

一维中只有点和线,可以用一根实数轴和依附于上面的单位刻度来表示数,如下图所示:

加法可以看作是把 0 向左或右移动到相加的结果的位置,比如 0+1=1​ 就是把 0 移动到右边原本1的位置,其他的位置点也同样向右移动相同的距离:

乘法可以看作是把 0 固定住,并以它为原点将 1 拉伸(缩放)到相乘的结果的位置,比如 1×2=2 就是把 1 拉伸到原本 2 的位置。可以看到空间已经被拉伸,或者说缩放了:

比起代数运算结果,我们更应该关注的是运算带来的移动、缩放的动作本身,它就是抽象结果在具象化的过程。在其他的维度空间也有着类似的这种动作。

3、二维

由于二维空间是个平面,到了这里终于可以有旋转的存在了。我们可以像上一个章节那样构造一个由两个数轴围起来的平面,在这里我们用复平面来表示,也就是一个实轴和一个虚轴。复平面我们放到复数的后面来论述,在此之前先看下二维旋转的分解。

3.1 旋转的分解

假设向量 v 从原本 x 轴方向绕原点旋转 θ 度到当前所处的位置,它可以被表示为 v=vx+vy . 也就是说,要求出旋转后的向量,可以将它分解为两个向量:向 x 轴投影的向量和向 y 轴投影的向量,然后将两者相加即可。根据向量投影公式[1] uv|u|2u , 可得:

(3.1.1)v=xv|x|2x+yv|y|2y

如下图所示:

3.2 旋转矩阵

也可以利用三角函数构造旋转矩阵。设旋转前的向量为 v ,它的坐标是 (x,y) ,模长记为 r ,并且与 x 轴的夹角是 α 度。旋转后的向量为 v ,旋转后的坐标是 (x,y) ,并且与 y 轴的夹角是 θ 度。可以得出:

x=r cos(α)y=r sin(α)

而对于旋转后的向量 v ,它的坐标为:

x=r cos(α+θ)y=r sin(α+θ)

根据三角公式可以变形为:

x=r (cos(α)cos(θ)sin(α)sin(θ))y=r (sin(α)cos(θ)+cos(α)sin(θ))

用原向量的 xy 坐标来代替 r cos(α)r sin(α) ,上式可以变形为:

x=x cos(θ)y sin(θ)y=x sin(θ)+y cos(θ)

用矩阵乘法来表达就是:

(3.2.1)[xy]=[cos(θ)sin(θ)sin(θ)cos(θ)][xy]

因此,这个旋转矩阵就是 M=[cos(θ)sin(θ)sin(θ)cos(θ)] .它表示向量绕着原点逆时针旋转。有时候会看到一些文献把旋转矩阵表达为 [cos(θ)sin(θ)sin(θ)cos(θ)]​ ,这表示它是顺时针旋转的。

如下图所示:

在二维旋转中,施加旋转的次序是无关紧要的。比如说原向量 v 以逆时针方向旋转 90 4次,最后再顺时针旋转 90 1次,所得结果是在 y 方向上。无论怎么组合这个旋转序列,所得结果都是一样的。也就是说,二维旋转满足交换律。

旋转矩阵表示的旋转,与式子 3.1.1 的向量投影表示法是等价的。很快我们还能看到复数与它的联系。

3.3 复数

数学上用 z=a+bi 表示复数,zC ,a,bR .其中的 a 是它的实部,b 是它的虚部,并且规定 i2=1 .复数的运算关系如下:

模长:

|z|=a2+b2

乘法运算:

(3.3.1)(1)z1z2=(a+bi)(c+di)(2)=ac+adi+bci+bdi2(3)=ac+adi+bcibd(4)=acbd+(5)(bc+ad)i

请注意,本文将不会对复数进行详细地阐述,只会对一些跟旋转有关以及与前后文有联系的部分进行论述。如果想了解更多复数的细节,可以到网上参考更多资料。

类似于在一维中构造数轴的方式,可以用复数来构造一个二维平面,如下图所示,横轴代表实轴,而竖轴则代表虚轴:

可以看到,复数 z 就是以 {1,i} 为基底进行的操作,这与矩阵非常类似。因此,式子 3.3.1 中的 z1 可以看作是一个矩阵,而 z2 可以看作一个向量,两者相乘相当于对 z2 执行了一次 z1 变换。式子 3.3.1 可以变形为:

z1z2=[abba][cd]

再来看下 z1 所代表的这个矩阵。可以在矩阵左侧乘以 a2+b2 ,矩阵里的元素都除以 a2+b2

(3.3.2)(6)z1=[abba](7)=a2+b2[aa2+b2ba2+b2ba2+b2aa2+b2](8)=|z|[aa2+b2ba2+b2ba2+b2aa2+b2]

观察上图,可以看出 aa2+b2 其实就是 cos(θ) ,而 ba2+b2 就是 sin(θ) .所以 3.3.2 的式子等价于:

z1=|z|[cos(θ)sin(θ)sin(θ)cos(θ)]

这个式子解释了复数乘法与二维旋转的关系。观察式子的右边,它跟我们在上个章节的式子 3.2.1 中推导出的旋转矩阵是一致的。也就是说, z1z2 这个运算是线性的,相当于把 z2 这个向量先逆时针旋转了 θ 度,再缩放了 |z| 倍。而如果 z1 的模长是1,那么 z1 代表的就是纯旋转

如果计算一下,会发现 z1z2z2z1 的结果是一样的[2],也就是说,两个复数相乘是满足交换律的。这也验证了 3.2 章节中提到的“二维旋转中施加旋转的次序是无关紧要的”说法。

通过观察可以得出这个旋转矩阵的复数形式是 z=cos(θ)+i sin(θ) ,可以把它概括成更一般的形式:

(3.3.3)(9)v=(cos(θ)+i sin(θ))v(10)=zv

3.4 复平面

复数在几何上有非常形象的解释,如下图所示:

由于规定了 i2=1 ,1 乘以 i 的平方(也就是连续执行乘以 i 这个动作两次)之后来到了 1 的位置,相当于旋转了 180 .那么,执行一次乘以 i 的动作就相当于旋转 90 ,此时 i=1 .由于 1​ 在实数内无解,只能将它放置在另一个位置——虚轴上。这种几何图称为复平面,也叫阿尔甘图,是以瑞士数学家 Robert Argand 的名字命名的,据说他是最早对复数进行几何描述的人之一(存疑)。

我们还可以参照一维那个章节里展示的思想,构造一张网格,代表复平面中的复数。加法可以看作是平移整张网格,并把 0 移到相加的结果;而乘法则是把 0 固定住,可以想象有一枚钉子钉在了 0 处,并以它为中心点旋转和拉伸整张网格,这与矩阵对空间的作用非常类似。如下图所示:

4、三维

4.1 万向节锁定的本质

与二维旋转不同,三维的旋转顺序是不能变换的,数学上的表达就是运算不满足交换律。比如一个六面不同颜色的魔方,按照 unity 场景视图中的坐标轴为旋转轴,按照 xyz 的顺序旋转,与按照 yxz​ 的顺序旋转,两者所得的结果是不同的。并且,二维坐标系没有偏手性,那是因为无论怎样构造的坐标系我们都可以通过旋转、镜像等变换让它们的指向性变成一致,但是三维坐标系却不能做到这一点,因此三维中的描述要声明偏手性。

使用欧拉角有万向节锁定的风险。在 unity 的场景视图中就可以“人为制造”一个万向锁。具体的做法是,在场景中放置一个立方体,然后在 Inspector 窗口的 Transform 中把 Rotation 的 x 设为 90 ,接着用鼠标在 y 或者 z 上拖动,会发现两者的旋转变成一样了。

对于万向节锁定的解释,有两种比较流行的说法。一种是从旋转矩阵出发,用代数计算的方式可以提供一种优雅而严谨的证明,这种解释需要构造三维旋转矩阵,有兴趣的读者可以在扩展阅读部分中阅读[3]。另一种是使用平衡环架的旋转特点来解释,也就是旋转轴会互相嵌套。这种解释有它独特的价值,因为它够直观也够明了。但是如果要深究的话可能仍然会有些疑惑,比如将 x 轴旋转 90 ,是因为 x 轴与 z 轴重合了,从而丢失了一个自由度——这是平衡环架理论的解释,但是这个解释的前提是旋转轴之间是互相嵌套的关系,也就是 z 轴嵌套于 y 轴而 x 轴则不然。而事实是,物体在旋转时,它的本地坐标轴是跟着转动的,不存在某个轴不转动的情况。真正导致问题的地方其实在于欧拉角的特性,即:欧拉角跟矩阵的意义一样,它描述的应该是变换,而不是动态的运动变化过程。它的变换方式,总是从初始姿态按照旋转顺序和给定参数直接变换成结果姿态。用表达式来描述的话,就是:假如用欧拉角系统执行 A,BCD 三次旋转,其中 A 是旋转前的姿态;B 是执行 B 变换后的姿态,其他的同理。它先是从 A 姿态执行 B 变换从而变成 B ;接着又从 A 姿态出发执行 C 变换从而变换成 C ,并不是预期的从 B 姿态出发再进行变换。最后,它仍然从 A 姿态出发执行 D 变换从而变换成 D .假设 B 变换是绕 x 轴旋转 90 ,那么将会:

  • 当我们尝试绕 y 轴旋转时,由于 x 轴旋转已经改变了物体的朝向,原本应该垂直于 x 轴的 y 轴旋转,现在看起来像是在绕 z 轴旋转。

  • 同样地,当我们尝试绕 z 轴旋转时,它实际上会表现出类似绕 y 轴的效果,因为 x 轴旋转已经使 y 轴和 z 轴的旋转效果变得难以区分。

虽然我们在直觉上认为我们是在绕不同的轴旋转,但由于每次旋转都是相对于初始姿态 A 进行的,所以在某些特定角度下,这些旋转的效果会重叠,从而导致我们失去了对其中一个维度的独立控制——这就是万向锁的本质。换句话说,即使我们意图分别绕三个不同的轴旋转,但由于每次都基于初始姿态,所以一旦某个轴的旋转达到了特定角度(如 90 ),后续的旋转就会受到之前旋转的影响,造成两个轴的旋转效果混淆,进而丧失了一个自由度。所以我们说各轴之间并不是像平衡环架一样互相嵌套,只是它的结果看起来跟平衡环架互相嵌套的结构展示出的过程是一致的。

4.2 三维旋转的分解

我们可以借鉴二维里分解旋转向量的思想,把三维向量也同样分解成对应的分量,从而用平面的方式解决问题。如下图所示:

图中是同一个轴角系统的两个不同视角,右边是俯视图。设三维向量 v (红色),绕单位向量 u 旋转 θ 度,旋转后的 vv (蓝色),它在二维平面上的投影是一个圆。v 分别向两个方向投影,得到 vv 两个分量。对旋转后的 v 也是采用同样的方法分解。可以很容易地知道, v 就是它的两个分量 v (图中没有画出)和 v 相加的结果,且 v=v.现在只需要求出 v 即可。

4.3 罗德里格旋转公式

为了求出 v ,让我们聚焦于上图的右侧,在圆所在平面上继续将 v 分解为 vωvv .显然 v=vω+vv ,且 ω=u×v .于是可得下式:

(4.3.1)(11)v=vv+vω(12)=cos(θ)v+sin(θ)ω(13)=cos(θ)v+sin(θ)(u×v)

根据上一节的向量分解,可得:

(4.3.2)(14)v=v+v(15)=v+cos(θ)v+sin(θ)(u×v)

根据 3.1 节提到的向量投影公式,可得 v

(16)v=uv|u|2u(17)=(uv)u

v 则是:

(18)v=vv(19)=v(uv)u

再看 ω

(20)ω=u×v(21)=u×(vv)(22)=u×vu×v(23)=u×v

将以上变量代入式子 4.3.2​ 可得:

(24)v=(uv)u+cos(θ)(v(uv)u)+sin(θ)(u×v)

将这个式子整理一下可得:

(4.3.3)v=cos(θ)v+(1cos(θ)(uv)u)+sin(θ)(u×v)

这就是罗德里格旋转公式,以法国数学家 Benjamin Olinde Rodrigues 的名字命名的。我们在后面论述四元数的章节里将会证明四元数乘法与此处的式子是等价的。

4.4 构造三维复数?

3.2 节中阐述了复平面,它用一实一虚两个轴构造了一个几何空间。那么,在三维里我们能不能也使用类似的方式,把复数升到三维,即用两个虚轴和一个实轴(三者互相垂直)构造起“三维复数”来描述空间旋转变换呢?实际上,爱尔兰数学家哈密顿(William Rowan Hamilton)也曾花费许多时间在寻找这种数字,而他始终找不到。类似于牛顿被苹果砸中脑袋这种浪漫的情节,他的故事也很有趣,据说他的儿子每天都例行公事般问他找到了没有,他总是回答:“还没。”

为什么用两个虚轴和一个实轴无法描述三维的旋转呢?这个问题有很多种说法,这里用一种几何上的解释。如下图所示:

图中的球体是单位球,由于在二维中复数的旋转路径是一个圆,所以在三维里“三维复数”的则是一个球体。现在假设它是存在的,让球面(二维)上的点乘以一个单位“三维复数” c ,这样的话每个点都会绕着球心旋转一定的角度,视觉上看着就像是沿着球面向某个方向移动了一定的距离。这里用橙色的箭头代表这个“移动”。请注意,这里是纯旋转,没有平移,这里所谓的“移动”只是一种感性的称呼。

请想象一下,乘法运算之后球面上所有的箭头(这里只画了几个)都会指向自身的唯一方向,除了两处:那就是上下两个极点。这两处的方向是不确定的,这就好比地球的经纬度描述在两个极点处会失效一样,用数学语言来描述就是:使用三个参数无法完整地覆盖三维旋转的参数空间。拓扑学上的“毛球定理”[4] 也可以说明这一点。具体地说明需要涉及到群论的知识,这里就不展开了。因此,这种“三维复数”不存在,或者说存在,但是它的性质没那么“稳定”。所以,想要完整地描述三维旋转还要再继续升维。

5、四维

既然“三维复数”不行,哈密顿转而寻找四维的数。1843年10月16日,在前往爱尔兰皇家科学院参会的路上,他灵光一闪,觉得应该使用三个虚数(i,j,k)加一个实数的方式,而不是两个虚数加一个实数,来描述三维旋转。他在布鲁姆桥上刻下了著名的式子:i2=j2=k2=ijk=1 .四元数诞生了。

5.1 四元数基本运算

对于四元数 qH ,可以表达为 q=ω+xi+yj+zk 的形式,其中 ω 是实部,xi,yj,zk 是虚部,ω,x,y,zR ,并且规定 i2=j2=k2=ijk=1 . 由于四元数乘法不满足交换律(但是满足分配律,结合律),所以左乘和右乘的结果是不同的。通过简单计算可以得出 i,j,k 之间的运算关系如下表格 5.1.1 所示[5]

i j k
i 1 k j
j k 1 i
k j i 1

注:上述表格乘法顺序按照“列_行”的顺序进行。

除了 q=ω+xi+yj+zk 这种记法,还有以下更简洁的形式:

q=[ω,v](v=[xyz], ω,x,y,zR)

其中的 v ,哈密顿将其称之为:向量。像 q=[0,v] 这种实部为0的四元数,称之为纯四元数。因此,任何向量都可以表示为纯四元数的形式:v=[0,v] .

由于三维坐标系中三个轴都被替代为虚数轴,而第四个轴,即实数轴,它存在于四维空间,以我们无法理解的方式同时正交于这三个本身互相正交的虚数轴。所以,四元数是以 {1,i,j,k}​ 作为基底的线性组合。

请注意,本文不会介绍四元数的矩阵形式、对数型、指数型、点乘和叉乘等概念,只会简单地介绍一些接下来要用到的基本运算。以下是四元数的一些基本计算方式:

  • 标准乘法:(由于 ω1,x1,y1,z1ω2,x2,y2,z2 这种变量名会让计算结果显得非常凌乱,此处用 a,b,c,de,f,g,h​​ 来代替)

(25)q1q2=(a+bi+cj+dk)(e+fi+gj+hk)(26)=ae+afi+agj+ahk+bei+bfi2+bgij+bhik+cej+cfji+cgj2+chjk+dek+dfki+dgkj+dhk2

使用表格 5.1.1​ 的计算结果进一步化简上式:

(27)=ae+afi+agj+ahk+beibf+bgkbhj+cejcfkcg+chi+dek+dfjdgidh(28)=(aebfcgdh)+(be+afdg+ch)i+(ce+df+agbh)j+(decf+bg+ah)k

v=[ecd] ,u=[fgh]​ ,利用向量的点乘和叉乘还可以把上式表达为更简洁的形式:

q1q2=[aevu, au+ev+v×u]

  • 模长:

|q|=ω2+x2+y2+z2

  • 共轭(记法为 q):

q=[ω,v]

(29)qq=[ω,v][ω,v](30)=[ω2v(v), ω(v)+ωv+v×(v)](31)=[ω2+vv, 0](32)=ω2+x2+y2+z2(33)=|q|2

(34)(q)=[ω,(v)](35)=[ω, v](36)=q

(37)qq=(q)(q)(38)=|q|2(39)=ω2+x2+y2+z2(40)=|q|2(41)=qq

  • 逆(记为 q1):

qq1=1

q1=q|q|2[6]

q1=q(如果q是单位四元数)

5.2 v的计算

这里到了最关键的环节。我们同样用 4.2​ 节的分解向量的方法来进行计算,不同的是之前使用的向量在这里被替换为纯四元数:

(42)v=v(43)=[0,v],(44)v=v(45)=[0,v],(46)v=v(47)=[0,v],(48)v=v(49)=[0,v],(50)v=v(51)=[0,v],(52)v=v(53)=[0,v],(54)u=u(55)=[0,u].

v=v+v, v 在旋转前后不变,比较关键的是求出旋转后的 v ,也就是 v .

由式子 4.3.1 可得:

v=cos(θ)v+sin(θ)(u×v)

由于

(56)uv=[0,u][0,v](57)=[uv, u×v](58)=[0, u×v](59)=u×v(60)=u×v

uv 代入上式,可得:

v=cos(θ)v+sin(θ)(uv)

使用分配律变形上式,可得:

(5.2.1)(61)v=cos(θ)v+sin(θ)(uv)(62)=(cos(θ)+sin(θ)u) v

p=cos(θ)+sin(θ)u ,可得:

(63)p=cos(θ)+sin(θ)u(64)=[cos(θ),0]+[0,sin(θ)u](65)=[cos(θ), sin(θ)u]

上式 5.2.1 可变形为:

(5.2.2)v=pv

将式子 5.2.2 代入以下式子可得:

(5.2.3)(66)v=v+v(67)=v+pv

p=qq ,则 q=[cos(θ2),sin(θ2)u] [7] ,上式可变形为:

(5.2.4)(68)v=v+v(69)=v+pv(70)=1v+pv(71)=qq1v+qqv

到这一步,四元数对向量(或者说纯四元数)的旋转作用已经非常直观了:对于 v ,它执行 q 之后又再执行 q1 ,相当于没有操作;而对于 v ,它执行两次 θ2 旋转,相当于旋转了 θ 度。

但是数学公式追求一种最简洁的美,上式还可以继续变形。

因为 q 是单位四元数[8] ,所以 q1=q .将其代入式子 5.2.4 可得:

(5.2.5)v=qqv+qqv

又因为 qv=vq [9], qv=vq [10],式子 5.25 可变形为:

(5.2.6)(72)v=qqv+qqv(73)=qvq+qvq(74)=q(v+v)q(分配律)(75)=qvq      (v=v+v)(76)=qvq1

至此,得出了上式形式的四元数一般性的乘法公式,非常简洁。

总结一下,设 v=[0,v], q=[cos(θ2),sin(θ2)u] ,u 为单位向量,则旋转前的 vu 轴旋转 θ后的 v 可用四元数乘法表示为: v=qvq=qvq1​ .

5.3 与罗德里格旋转公式的比较

以上的四元数乘法,其结果等价于之前论述的罗德里格旋转公式(式子 4.3.3)。下面我们就来论证这一点。

证明:

q=[cos(θ2),sin(θ2)u],v=[0,v],q=[cos(θ2),sin(θ2)u] .则

(77)qv=[cos(θ2),sin(θ2)u] [0,v](78)=[sin(θ2)uv, cos(θ2)v+sin(θ2)u×v].

(1)qvq=[sin(θ2)uv, cos(θ2)v+sin(θ2)u×v] [cos(θ2),sin(θ2)u].

由于罗德里格旋转公式是旋转后的向量表示,对应于四元数里的虚部。因此,只需要计算式子 (1) 的虚部即可。

(2)(79)Im=(sin(θ2)uv)(sin(θ2)u)+(80)cos(θ2)(cos(θ2)v+sin(θ2)u×v)+(81)(cos(θ2)v+sin(θ2)u×v)×(sin(θ2)u)(82)=sin2(θ2)(uv)u+(83)cos2(θ2)v+(84)cos(θ2)sin(θ2)u×v+(85)cos(θ2)v×(sin(θ2)u)+(86)(sin(θ2)u×v)×(sin(θ2)u)

其中,

(3)(87)cos(θ2)v×(sin(θ2)u)(88)=cos(θ2)sin(θ2)(v×u)

根据向量三重积公式 (a×b)×c=(ac)b(bc)a 可得:

(4)(89)(sin(θ2)u×v)×(sin(θ2)u)(90)=sin2(θ2)v+sin2(θ2)(vu)u

(3),(4) 代入 (2) ,可得:

(5)(91)Im=sin2(θ2)(uv)u+(92)cos2(θ2)v+(93)cos(θ2)sin(θ2)(u×v)(94)cos(θ2)sin(θ2)(v×u)(95)sin2(θ2)v+(96)sin2(θ2)(vu)u

因为v×u=u×v(反交换律), (5) 中的 cos(θ2)sin(θ2)(v×u) 等于:

(6)cos(θ2)sin(θ2)(u×v)

(6) 代入 (5)​ 并整理可得:

(7)(97)Im=cos2(θ2)vsin2(θ2)v(98)+2sin2(θ2)(uv)u(99)+2cos(θ2)sin(θ2)(u×v)

利用三角公式化简 (7) ,可得:

Im=cos(θ)v+(1cos(θ))(uv)u+sin(θ)(u×v)

则所得结果与罗德里格旋转公式一致。

证毕。

6、可视化

四元数藏匿于四维空间,想实现真正的可视化是不现实的,我们只能通过其在三维空间中的投影来建立起对它的几何直觉。关于四元数可视化的这部分,建议观看 3Blue 1Brown 的视频,更直观:四元数的可视化_哔哩哔哩_bilibili 。但是应该提醒自己,关注点不是作者的酷炫动画展示,而是他在视频中提到的左乘和右乘的不同影响。直观去感受四元数旋转在三维中的投影变化方式对理解它没有什么大的帮助,因此本文将补充一点视频中没提到的部分。以下部分图片来源于他视频中的截图。

通常采用“球极投影”(Stereographic Projection)的方法来实现,如下图所示:

从球面的极点出发引一条射线,穿过球面上的某一个点最后投射到平面上,则球面上的点被投影到平面。如上图黄色的线是投影线,C 点被投影为平面上的 E​​​ 点。按照这个方式,我们可以把球面上的每个点都投影到平面上,除了作为射线出发点的极点——它被投影到无穷远的地方。出发极点会把自身所在的半球投影到平面上赤道的投影之外的地方,球体的赤道则会保持完整的形态不变。

视频中作者从 1 点,也就是南极点向上投射投影线。

如上图所示:通过 {k,i,k,i} 的圆会被投影成通过 {i,i} 的直线,因此当这个圆发生旋转时,对应地,这条直线上出现了在一维的章节中提到过的轴上的点的移动。通过 {k,k}{j,j} 的也是如此。所以当看到这三条轴上的点在朝一个方向“前进”时,我们应该提醒自己它实际上对应的是一个圆的旋转。

三维球体投影到二维平面上保持不变的部分是一个圆(赤道),而四维超球体投影到我们三维空间中,保持不变的部分则是一个同时穿过 {i,i,j,j,k,k} 的三维球体。这个球体上的球面,代表的是实部为 0 的四元数,即纯四元数,要发生旋转的向量需要沿着这个球面进行才是纯旋转,如果脱离这个球面就是混合着缩放或平移等的其他变换。脱离这个球面的部分代表的是实部不为 0​ 的四元数,它们在三维中的投影以一种我们看不见的或者我们不好理解的扭曲方式存在,正如投影线出发的极点所对立方向的半球被投影到赤道之外而发生扭曲一样。

让我们先考虑三个虚轴上的情景。当一个单位四元数左乘一个四元数时,如 ip , i 可以看作是一个会产生变换的函数p 则可以看作一个被发生变换的。此时 {j,k,j,k} 圆发生旋转,{i,i} 的点在移动(实际上是一个圆的旋转),也就是把 1 移动到所乘的结果的位置,因此它所代表的是一个互相垂直且同时发生的双旋转。请注意,这是四维空间中发生的旋转,当前的画面只是它投影在三维中的情形。之所以会有双旋转,是因为 {i,i} 作为旋转轴可以确定两个平面:一个是以它为法向的过原点的平面;另一个是它与实轴张成的复平面。这两个平面内的旋转是同时发生的。它遵循右手法则,拇指指向作为函数的运动方向;四指弯曲的方向则是 {j,k,j,k} 圆旋转的方向,它作为被函数变换了的的结果,也就是表格 5.1.1 中的计算结果。请注意,视频中作者使用的是右手坐标系。如下图所示:

但是如果是右乘,则是相反的操作,此时遵循左手定则。

其他的两个虚轴也是类似的。

除了三个虚轴,还可以扩展到任意轴:

现在回到四元数乘法公式上:qvq 中,可以把 v 看作是上述中的,左乘一个作为函数作用的 q 和右乘一个作为函数作用的 q .设要旋转的角度是 θ 度,旋转轴是 1 和球面上 b 的连线,当左乘 q 时,即 qv ,会遵循右手定则绕着 {1,b} 旋转轴(就是 1q 的连线,它代表的是一个圆,就像上面双旋转中的情景一样)把 v 拽到脱离单位球面的位置(实部小于 0 的在球面之外,实部在 [0,1) 之间的在球面之内),注意此时的操作是混合着缩放和旋转的复合变换,因为它的实部已经不为 0 ;由于它是一个点,所以投影在单位球面上的是一个旋转了 θ2 度的效果。但如果是一个向量,此时投影在单位球面上的将会是被缩放过的向量,而不是纯旋转。如下图所示(以下图片来自 B站 Up主“开心expz”的视频截图):

当右乘 q 时,即 vq ,会遵循左手定则把原本被拽到脱离单位球面的 v 重新拉回单位球面上,这意味着它的实部重新变为 0 ,此时不带缩放效果并且再次旋转 θ2 度,两次的叠加正好是旋转了 θ 度。如下图所示:

以上的左乘和右乘对点(向量)的影响,才是四元数可视化中最重要的地方。它在三维中炫酷的展开方式只是一种几何直觉,对理解四元数没有多大的帮助。并且,从视频里也能看出,当进行四元数乘法计算时,每个点都在移动(旋转),这也说明其不会遇到 4.4 章节中提到的构造点的问题。

如果以上图中展示的方法还是感觉太麻烦,还有一种十分简便的方法来理解:用您的双手,保持每只手的大拇指都和其余四指垂直,弯曲四指,然后让双手的小拇指弯曲的形状重合。此时,两个大拇指的方向是相反的,因此,大拇指方向的旋转互相抵消,而四指弯曲方向的旋转叠加了。

7、代码实现

得益于四元数乘法公式的简洁性,它在代码中实现非常地简单,比较需要注意的地方是在归一化的操作中要注意避免除零的计算风险。在本节代码中只计算纯旋转的变换,也没有提供四元数的求负、加减等操作,而且没有涉及多个四元数的复合变换,也没有提供双倍覆盖的处理方式。所以这并不是一个很健壮的代码。本文没有对四元数插值进行论述(也许以后有机会会写一写这方面的文章),因此本节代码也将不加入插值的函数。

using System;

public struct Vector3 
{
    public float X, Y, Z;

    //Vector3构造函数
    public Vector3(float x, float y, float z)
    {
        X = x;
        Y = y;
        Z = z;
    }
}

public class Quaternion 
{
    //四元数实部、虚部属性
    public double W { get; private set; }
    public double X { get; private set; }
    public double Y { get; private set; }
    public double Z { get; private set; }

    //Quaternion构造函数
    public Quaternion(double w, double x, double y, double z)
    {
        W = w;
        X = x;
        Y = y;
        Z = z;
        Normalize();
    }

    //轴角对的计算函数
    public static Quaternion FromAxisAngle(Vector3 axis, double angleInRadians) 
    {
        double halfAngle = angleInRadians / 2;	//注意公式qvq*中一个q的角度是θ/2
        double sinHalfAngle = Math.Sin(halfAngle);
        return new Quaternion(
            Math.Cos(halfAngle),
            axis.X * sinHalfAngle,
            axis.Y * sinHalfAngle,
            axis.Z * sinHalfAngle
        );
    }

    //归一化四元数
    private void Normalize() 
    {
        double lengthSquared = W * W + X * X + Y * Y + Z * Z;
    
        // 检查是否为零向量
        if (Math.Abs(lengthSquared) < 1e-6)     //判断模长平方的绝对值是否小于0.000001
        {
            // 如果是零向量,可以选择不同的处理方式。
            // 这里简单粗暴地将其设为单位四元数 (1, 0, 0, 0)
            W = 1.0;
            X = Y = Z = 0.0;
            return;
        }

        double length = Math.Sqrt(lengthSquared);

        // 检查长度是否接近1
        if (Math.Abs(length - 1.0) > 1e-6)  //大于1必须初始化
        {
            double invLength = 1.0 / length;
            W *= invLength;
            X *= invLength;
            Y *= invLength;
            Z *= invLength;
        }
    }

    //共轭,也就是q*
    public Quaternion Conjugate() 
    {
        return new Quaternion(W, -X, -Y, -Z);
    }

    //执行四元数乘法公式:qvq*
    public Vector3 RotateVector(Vector3 v)
    {
        // 将向量转换成纯四元数
        Quaternion qv = new Quaternion(0, v.X, v.Y, v.Z);
        Quaternion qConj = this.Conjugate();
        
        // 应用旋转: q * v * q*
        Quaternion result = Multiply(Multiply(this, qv), qConj);
        
        // 提取旋转后的向量部分
        return new Vector3((float)result.X, (float)result.Y, (float)result.Z);
    }

    //四元数标准乘法
    private static Quaternion Multiply(Quaternion a, Quaternion b) 
    {
        return new Quaternion(
            a.W * b.W - a.X * b.X - a.Y * b.Y - a.Z * b.Z,
            a.W * b.X + a.X * b.W + a.Y * b.Z - a.Z * b.Y,
            a.W * b.Y - a.X * b.Z + a.Y * b.W + a.Z * b.X,
            a.W * b.Z + a.X * b.Y - a.Y * b.X + a.Z * b.W
        );
    }
}

使用示例:

 		// 定义一个要旋转的向量
        Vector3 vectorToRotate = new Vector3(1, 0, 0);

        // 创建一个绕Y轴旋转90度的四元数
        Quaternion rotation = Quaternion.FromAxisAngle(new Vector3(0, 1, 0), Math.PI / 2);

        // 应用旋转
        Vector3 rotatedVector = rotation.RotateVector(vectorToRotate);

        Console.WriteLine($"Rotated Vector: ({rotatedVector.X}, {rotatedVector.Y}, 			{rotatedVector.Z})");

8、扩展阅读

8.1 推导部分

[1] 向量投影公式的推导

​ 以其中 v 为例,因为 uv 同向,所以 v 可以看作是一个标量乘以 u ,即:

v=ku(1)

​ 现在的问题就是找到这个标量k .因为 uv 同向(cosθ=1),所以

(2)(100)vu=|v||u|cos0(101)=|v||u|

​ 将(1) 代入(2) ,可得 knu=|v||u| .因为 kuu=k(uu)=k|u|2 ,所以

(102)k|u|2=|v||u|(103)k=|v||u||u|2(104)=|v||u|

​ 把这个 k 代入(1) ,可得:

(3)v=u|v||u|

​ 根据三角函数的定义可得 |v|=cosθ|v| ,将其代入 (3) ,可得:

v=ucosθ|v||u|

​ 将其上下同时乘以 |u| ,可得

(105)v=ucosθ|v||u||u|2(106)=uvuuu

[2] 复数满足交换律

(107)z1z2=(a+bi)(c+di)(108)=acbd+(ad+bc)i

(109)z2z1=(c+di)(a+bi)(110)=acbd+(ad+bc)i

​ 可见两者的结果是一致的。

[3] 首先构造三维旋转矩阵。在 3.2 节中我们已经构造了一个二维旋转矩阵,在三维中绕某个轴的旋转一样可以看作 是在二维平面内(向量与旋转轴围起来的平面)的旋转,如此一来二维旋转矩阵仍然适用。我们已经知道绕着某 个轴旋转,这个轴是不变的,因此可以分别为 x,y,z 轴构造绕着它们的旋转矩阵:

​ 首先是 x 轴,因为旋转前后它是不变的,因此是

Rx(θ)=[1000cos(θ)sin(θ)0sin(θ)cos(θ)]

​ 同理可以构造 yz :

Ry(θ)=[cos(θ)0sin(θ)010sin(θ)0cos(θ)]

Rz(θ)=[cos(θ)sin(θ)0sin(θ)cos(θ)0001]

​ 假设现在 x 轴旋转 α 度, y 轴旋转 90 度, z 轴旋转 θ 度:

(111)Rz(θ)Ry(π2)Rx(α)(112)=[0cos(θ)sin(α)cos(α)sin(θ)sin(α)sin(θ)+cos(α)cos(θ)0sin(α)sin(θ)+cos(α)cos(θ)cos(α)sin(θ)cos(θ)sin(α)100](113)=[0sin(αθ)cos(αθ)0cos(αθ)sin(αθ)100](114)=Ry(π2)Rx(αθ)

​ 可见表达式中丢失了 Rz .

[4] 毛球定理(Hairy Ball Theorem)是拓扑学中的一个定理,它指出在二维球面上,不存在连续的、非零的切向量 场。换句话说,如果你试图将毛发均匀地分布在球体表面(想象成在球上刷油漆或梳理毛球),总会有一个地方 毛发会竖起来或者聚集在一起,形成所谓的“旋涡”或“牛眼”。

​ 具体来说,如果我们尝试定义一个从实数到 SO(3)(三维空间中所有可能旋转组成的特殊正交群)的连续映射, 使得每个点对应于一个特定的角度和轴线上的旋转,那么根据毛球定理,这样的映射必然会在某些点失效。这是 因为 SO(3) 的空间拓扑结构本质上是与三维球面 S3 相关的,而 S3 又可以通过去掉一对对径点后的 S2 来近似 理解。因此,在尝试构造这样一个连续映射时,总会有某些方向或角度无法被平滑地覆盖,导致出现不连续性或 奇异性。

[5] i2=j2=k2=1 ,ijk=1

​ 两边同时左乘 i​ ,即:

(1)(115)iijk=i(1)(116)jk=i(117)jk=i

​ 两边同时右乘 k ,即:

(2)(118)ijkk=(1)k(119)ij=k(120)ij=k

​ 同理,由(1)可以继续变形:

(121)jjk=ji(122)k=ji

(123)jkk=ik(124)j=ik

​ 同理,由(2)可以继续变形:

(125)iij=ik(126)j=ik

(127)ijj=kj(128)i=kj

[6] 逆定义为 qq1=1 ,表达式两边同时左乘 q​ ,即:

(129)qqq1=q1(130)(qq)q1=q(131)|q|2q1=q因为qq=|q|2(132)q1=q|q|2

[7] 设 p=[cos(θ),sin(θ)u] ,则

(133)pp=[cos(θ),sin(θ)u] [cos(θ),sin(θ)u](134)=[cos2(θ)sin2(θ)(uu), 2cos(θ)sin(θ)u+sin(θ)u×sin(θ)u]

​ 因为 u×u=0 ,uu=|u|2=1​ ,根据倍角公式可得:

pp=[cos(2θ), sin(2θ)u]

​ 可知乘一次 p 就是对 θ​ 追加一次累加。

[8] q​​ 是单位四元数

​ 设 q=[cos2(θ2)+sin2(θ2)u] ,也就是 q=cos(θ2)+sin(θ2)uxi+sin(θ2)uyj+sin(θ2)uzk ,则

(135)|q|=cos2(θ2)+sin(θ2)usin(θ2)u(136)=cos2(θ2)+sin2(θ2)(uu)(137)=cos2(θ2)+sin2(θ2)(138)=1

[9] qv=vq

​ 为简便起见,设 q=[ω,αu] ,v=[0,v] .因为 uv 同向,所以 u×v=0 ,点积满足交换律,则

(139)qv=[ω,αu] [0,v](140)=[αuv, ωv]

(141)vq=[0,v] [ω,αu](142)=[αuv, ωv]

​ 可见两者结果是一致的。

[10] qv=vq

​ 同样设 q=[ω,αu] ,v=[0,v] .因为 uv=0 ,叉乘满足反交换律,则

(143)qv=[ω,αu] [0,v](144)=[0, ωv+αu×v]

(145)vq=[0,v] [ω,αu](146)=[0, ωv+αu×v]

​ 可见两者结果是一致的。

8.2​ 参考资料

[1] Ian Parberry, Fletcher Dunn. 3D数学基础:图形与游戏开发(第2版)[M]. 北京: 清华大学出版社, 2020-05.

[2] Eric Lengyel. 3D游戏与计算机图形学中的数学方法(第3版)[M]. 北京: 清华大学出版社, 2016-05.

[3] HMMNRST. Unityで3Dモデルの法線を反転させる方法[EB/OL].[2019-05-22].https://qiita.com/HMMNRST/items/0a4ab86ed053c770ff6a.

[4] 视频:四元数的可视化_哔哩哔哩_bilibili

[5] 视频:三维旋转需要经过四维空间?关于四元数_哔哩哔哩_bilibili

 posted on   小花猫喵喵  阅读(182)  评论(1编辑  收藏  举报
相关博文:
阅读排行:
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· C# 集成 DeepSeek 模型实现 AI 私有化(本地部署与 API 调用教程)
· spring官宣接入deepseek,真的太香了~
点击右上角即可分享
微信分享提示