图形学中的旋转(Rotation)非常重要,用欧拉角(Euler Angle)旋转非常简单,但是却存在 万向锁(Gimbal Lock) 的问题;同时,用欧拉角插值也不尽方便。基于四元数(Quaternion)的旋转既解决了万向锁的问题,又能非常便利地插值。本文将从概念出发,对四元数的定义、推导、性质、应用,以及它与欧拉角之间的联系进行介绍。本文主要参考了Krasjet的《四元数与三维旋转》1,其他参考内容包括2345678。
复数与二维旋转
首先我们将简要介绍复数及它与二维旋转的关系,从而可以更加自然地过度到四元数及其与三维旋转的关系。
定义
复空间
从线性代数的角度来看,复数
复数的加法与乘法
复数的加法
复数的加法很简单,分别把实部和虚部相加即可。减法类似,不多赘述。
复数的乘法
对复数
利用上面讲过的复数的向量形式,可以把上式写成矩阵乘向量的形式:
右侧的向量
所以,我们可以把
所以,我们现在可以用至少三种形式去表达一个复数了(具体用哪种形式取决于我们的需要):
- 代数形式:
; - 向量形式:
; - 矩阵形式:
.
有两个特殊的复数
对复数
复数的模长
复数
这可以用它的共轭(Conjugate)表示:
复数与二维旋转
从复数矩阵形式导出旋转关系
下面将重点引出复数与二维旋转的关系。
注意观察复数
矩阵中的每一项都被模长
这样一来,我们就可以把上式写成:
这样一来,任意复数
所以,对复数
我们可以构造一个复数
复数的极坐标形式
当我们看到
再根据欧拉公式(Euler's Formula):
可以立得复数的极坐标形式:
现在,复数
如果仅需表示旋转,则令
Proof of Euler's Formula:
令
复数的几种表示形式总结
现在我们已经得到了如下几种复数
- 代数形式一:
; - 向量形式:
; - 矩阵形式一:
; - 矩阵形式二:
; - 代数形式二:
; - 极坐标形式:
.
这几种形式是完全等价的。
旋转复合
假设我们现在有两个仅表示旋转的复数
所以,用两个复数进行旋转时,所得到的结果仍然是一个旋转,且旋转的角度是两次旋转角度之和。
基于Axis Angle的旋转
在这篇博客中,我们已经介绍了三维空间中一个向量绕任意轴旋转的公式(被称为Rodrigues' formula)。因为下文会涉及大量此公式,因此把结论在此重新记述,推导过程请参考原文。
记
或者矩阵形式:
因为这种旋转方式是通过一个旋转轴和一个旋转角度定义的,所以我们称之为Axis Angle旋转。 请记住此公式,我们将看到它与四元数之间紧密的联系。
欧拉角与万向锁(Gimbal Lock)问题
Unity中的Gimbal Lock现象
在正式进入四元数之前,我们先介绍欧拉角及其引发的万向锁问题。网上有很多关于Gimbal Lock的文章,但是找了一通发现大都表述得非常抽象。下面我把这个现象放到Unity中可视化展示,读者应该就能立刻知道问题所在。
注意左上角红框圈出的部分,X/Y/Z=1表示这个Cube正在绕世界坐标系的x/y/z轴旋转,在我们一开始就把Cube绕x轴旋转90度时,我们惊异地发现,旋转y轴或者z轴的效果都是一样的!显然这不是我们想要的结果。
1 | using System.Collections; |
那么为什么会出现Gimbal Lock现象呢?核心原因在于欧拉角旋转顺序是被固定的。要理解这句话,我们就要首先知道什么是欧拉角,什么叫固定顺序旋转。
欧拉角、坐标系与旋转顺序
如果你没有任何基础知识取看Wikipedia的Euler Angles页面,你一定会感到一脸懵逼。从严格的定义来看,Wikipedia的介绍比较准确,但是在这里我们只需要简单地把欧拉角理解为点绕着世界坐标轴旋转的角度即可。
欧拉证明了三维空间中的任意旋转都可以拆分为沿着自身的三个坐标轴的旋转,欧拉角就是这三次旋转的角度。但是在这里我们使用相对于世界坐标系的旋转角度,这是因为相对自身坐标系的旋转与相对世界坐标系的旋转是可以互相转换的,只需要把旋转矩阵按照相反顺序相乘即可。
我们把绕着自身坐标系旋转称为intrinsic rotation,把绕着世界坐标系旋转称为extrinsic rotation。
假定点
证明:
我们知道在世界坐标系中绕
用
现在考虑在点
这里要稍微做一些解释。粉色的
可以看到在绕
可以看到上图进行了两次旋转操作,第一次是绕着自身(其实也是世界坐标系)的
那么,绕自身
这个过程可以用下图表示:
首先我们绕
如果你还不理解,我们再来看一下最后公式的粉色部分
有了上述结论,我们在这里就不再关心是局部坐标系还是全局坐标系了,我们默认使用全局坐标系。
此外,顺序对基于欧拉角的旋转来说至关重要,这是因为矩阵往往不具备交换律,即
Unity中的欧拉角
Unity中的欧拉角是采用
这是因为当我们变动欧拉角的时候,Unity会重新按照
Unity中的Gimbal Lock问题
有了上面的铺垫,我们就知道为什么会出现开头那样的Gimbal
Lock现象了。当我们固定
我们发现,这个旋转矩阵竟然仅仅是旋转了
所以,由此可见,Gimbal
Lock产生的根本原因在于采用固定的旋转顺序,因为旋转顺序固定,所以就会在某些情况(如此处
除了Gimbal Lock之外,Euler Angles还存在插值困难和歧义性(与旋转矩阵之间的对应关系复杂)等问题。
四元数
定义
四元数可以定义为如下形式:
这个定义和复数非常相似,只不过四元数有三个“虚部” (
上述的
如果把上式称为四元数的代数形式,那么同样地,它也可以有对应的向量形式一:
如果把四元数
所以,仅仅从定义来看,四元数就有如下三种等价的形式:
- 代数形式:
; - 向量形式一:
; - 向量形式二:
.
加法、乘法、模长
加法
给定两个四元数
用四元数的向量形式
乘法
四元数的标量乘法和实数、复数一样,可以定义为:
然而四元数之间的乘积可以通过如下步骤定义:
我们可以从性质
于是能够得到下表:
容易观察到交换律并不成立,即
现在,得到的这个新四元数
或者可以写成下述的矩阵形式:
这样一来,我们就把两个四元数的乘积写成了矩阵与向量的乘积,并且由于四元数的乘法不可交换,
Graßmann积
显然矩阵形式又臭又长,我们希望用一种更简单的形式——Graßmann积——去表示四元数乘法。
重新整理之前的乘法结果:
如果令
这样一来,我们就可以用点乘和叉乘的形式表示四元数的乘积:
这就是所谓的Graßmann积(Graßmann Product)。所以,从Graßmann积也可能看到两个四元数的乘积一般是不可交换的,因为叉乘一般不可交换。
模长
四元数的模长和向量的模长定义一样:
纯四元数
定义“实部”为零的四元数为纯四元数,即:
对于两个纯四元数
共轭
四元数的共轭(Conjugation)与复数的共轭类似,可以定义为:
四元数与它自己的共轭相乘可以得到一个标量:
正好是四元数模长的平方。同时:
这表明四元数与它共轭的乘法是可交换的。这是一个非常好的性质。
逆
有了共轭,我们最后就可以来定义四元数的逆(Inverse)了。记四元数
利用共轭,我们可以进行如下推导:
所以,四元数的逆就是它的共轭除以它模长的平方。如果
四元数与三维旋转
从三维旋转公式到四元数表示
现在我们正式考察四元数与三维旋转之间的关系。首先,让我们回忆一下普通3D旋转的推导过程:
A sketch derivation of Rodrigues' formula:
记
现在我们想要找到四元数与上式的关系,我们自然地想到能不能像推导四元数乘法的矩阵形式那样,把上面所有的向量都用四元数去表示。于是我们可以定义:
从而我们有:
前面两项我们已经顺利地用四元数表示出来了,但是第三项的叉乘怎么办?我们先尝试一下用对应的四元数相乘看看:
这里就用到了纯四元数的概念。不管怎样,我们很“幸运”地用四元数的乘法把最后的叉乘表示了出来,代入
这里的
四元数旋转公式的简化
现在我们也可以像之前那样,分别把
引理一
如果
证明从略。
现在我们引入一个新的四元数
所以,之前的结果现在就可以重写为:
这里用
但是这种形式还是不能简化右式,我们需要下面的引理二。
引理二
如果
该引理证明如下:
根据引理二,我们就可以把上式的第一项最后两个因子交换,得到:
我们再使用引理三进行最后的化简。
引理三
如果
证明与引理二过程类似,不多赘述。根据这个引理,我们就能对原式做最后的化简:
最后我们得到了一个非常优美的结果,在这个公式中,我们不再把待旋转向量拆分为两个分量,而是直接通过四元数之间的乘法实现了旋转。该结论可以总结为一个定理:
三维旋转定理(四元数型)
3D Rotation Theorem (Quaternion): 任意三维向量
现在我们进行验证:
正好就是Rodrigues' formula!
并且,从上面的推导我们还知道:所有的单位四元数都对应了一个三维旋转(但这二者之间并非双射)。那么给定单位四元数
因此,我们得到了:
四元数乘积的性质
下面我们介绍几个关于四元数乘积的性质,这些性质与我们的主题没有直接关联,但对我们深入理解四元数大有帮助,供感兴趣的读者参考。
Property 1:
令
证明:
此外,还有
Property 1告诉我们:对四元数的旋转不改变标量和向量的模长。这与四元数表示旋转的本质是符合的。
Property 2:
令
证明:
其中第四个等号和第七个等号可以通过Property 1推知。
Property 2告诉我们:对一个单位四元数的
接下来的四个性质与单位四元数的求导有关。
Property 3 (Derivative of exponential):
令
证明:
Property 4 (Product rule): 设
证明:
Property 5 (Chain rule): 设
证明:
Property 6: 设
证明:
最后一个Property联系了四元数的乘积与点乘。
Property 7: 设
证明:
使用恒等式
最后,我们使用下述等式:
现在我们有:
三维旋转定理(矩阵型)
通过之前的讨论(尽管不完全相同,但可以类似推导)我们知道左乘一个四元数
而右乘
现在我们就可以利用这个结论把
实际上我们不需要存这个完整的
3D Rotation Theorem (Matrix):
任意三维向量
其中
尽管看起来比较复杂,但是在实际中如果旋转轴
三维旋转定理(指数型)
类似复数,我们也可以把四元数写成指数的形式。如果
也就是说,
3D Rotation Theorem (Exponential):
任意三维向量
其中
有了指数形式,我们就可以很自然地运用关于指数的一些运算了。比如对单位四元数
下面是幂运算;
即一个单位四元数的
旋转复合
现在来考虑基于四元数的旋转复合,假设分别在四元数
上式隐含了
需要注意的是,复合后旋转
双倍覆盖
单位四元数与三维旋转并不是双射(一对一)的,一个单位四元数对应一个三维旋转,但是一个三维旋转可以找到两个单位四元数去对应它,这就是所谓的2对1满射同态,或者说单位四元数双倍覆盖了三维旋转。
具体来说,
即
现在我们已经知道四元数相比欧拉角而言的优点了:
- 插值非常简单(将在下一节介绍);
- 几何意义明显,双倍覆盖三维旋转;
- 与旋转的顺序无关,一次四元数乘法即可完成旋转;
- 约束少,仅要求为单位四元数;
- 不存在Gimbal lock;
- 复合旋转简单。
四元数插值
基于幂运算的插值
假设有两个旋转变换
然而直接求这个中间旋转似乎比较困难,我们不妨从结果出发,期望原始向量
根据旋转定理的指数型,我们只要对
这个插值方法称为Slerp,但它涉及到多个四元数的乘法,而且还包含幂运算,实际应用中效率很低。我们一般的插值都是使用线性组合,而不是幂运算。
四元数在四维空间内的夹角
公式
我们有;
如果我们把
与此同时,我们又知道
Lerp, Nlerp, Slerp
前文已经提到,我们希望对四元数进行线性插值,即表示为
Lerp
第一种是我们熟知的线性插值Lerp:
但是,这种插值方法得到的四元数并不是单位四元数,也就是并不能真正表示三维旋转。现在我们来验证一下:
显然,当
Nlerp
既然Lerp的结果不是一个单位四元数,那么除以它的模长不就可以变成单位四元数了吗?这就是所谓的Nlerp(Normalized Lerp):
但是Nlerp的问题在于,随着
如下图所示,
Slerp
所以,我们想要直接对角度插值。如果
因为对角度线性插值相当于是让向量在球面上的一个弧上旋转,所以这种插值被称为球面线性插值Slerp(Spherical Lerp)。我们最开始说到的基于幂运算的插值就是我们接下来要介绍的Slerp的一种等价形式。
对于向量
为了求解
同样地,对两边同乘
现在,把
把
于是,我们能就得到最终Slerp的公式:
现在的Slerp已经比最开始介绍的基于幂运算的Slerp快很多,但是由于它涉及了三个三角函数、一个反三角函数和一次点乘,所以效率仍然比Nlerp低一些。此外,如果夹角
需要注意的是,Slerp并不是一个“万能答案”,它的运行效率显著低于Nlerp。此外,Nlerp的角速度变化不一致并不一定是个坏处,在不同的情景下可能有不同的用处。所以,在使用插值方法的时候要根据实际情况选择,大多数时候Nlerp已经足够了。
双倍覆盖问题
之前我们介绍过,四元数
如下图所示,虽然可以从
这告诉我们,在对四元数插值之前,要检查
Squad
Slerp已经是我们理想中的插值方式了:它直接对角度插值,插值角度匀速变化,运算效率尚可。但是它还有一个小问题:角度变化的速率等于夹角,即
为此,我们希望插值的曲线能够在高阶导处也连续,下面介绍的Squad(Spherical and quadrangle)就是一种解决方法。
下图是Slerp和Squad的插值结果比较,曲线是每个四元数插值结果的角速度。可以看到:Slerp在四元数切换时角速度会发生突变,而Squad的角速度则是光滑变化的。
Bézier曲线
假设我们有一个向量序列
显然这个曲线是
关于Bézier曲线的相关知识可以在Wikipedia上找到介绍。
三次Bézier曲线
为了解决这个问题,我们想让生成的曲线(1)经过每个点和(2)一阶导数连续。因此,我们可以分段对每两个向量
上述做法可以保证曲线经过每个点,但是不能保证曲线在每个点处都是
蓝色的线就是曲线在点
下面我们来具体了解如何构造一个三次Bézier曲线。
de Casteljau算法
de Casteljau算法可以非常直观简单地构造Bézier曲线,在这里我们只关注三次Bézier曲线的情况。
假设我们有四个向量
之后,我们对
最后,我们对
整个过程可以总结为下图(
当我们把
三次Bézier曲线的过程可以表述为:
这里的
Squad
三次Bézier曲线实际上是嵌套了三层一次(one-order)插值,而Squad则使用的是一层二次插值嵌套了一层一次插值。
我们首先是分别对
之后,我们使用
上述过程可以通过下图阐明(
当然,我们也可以把它写成递归形式:
我们成功地将原先的七次Lerp减少到了三次,且生成的曲线与原先的曲线差别也不大,可以看下图的两个对比(左边为Bézier曲线,右边为Squad曲线:
按照Lerp的定义对上式展开,就能得到:
它仍然是一个三次曲线,只不过系数不同。
把它应用到四元数,就能得到我们想要的四元数插值版本:
我们又知道
这个形式对我们后面解出控制点非常有帮助。
Squad应用
现在我们可以使用上式对多个单位四元数进行插值了。假设我们有一个四元数序列
接下来就是要找出控制点
找出
注意到上面的
通过上面的等式,我们可以求出
下面是上述结果的推导过程,篇幅较长,不感兴趣的读者粗看即可。所以,和原始版本的三次贝塞尔曲线(包含
A derivation of control point
为方便起见,我们记
因为
所以我们就能求出
从而,我们把
我们再把这个结果代回
同样地,我们可以求出
我们令两式相等,即得:
最后我们再根据
你可以验证最后一个等式的正确性。
与两个四元数插值一样,Squad在多个四元数插值时仍然会有双重覆盖的影响,在计算控制点和插值之前,需要先选中一个四元数,检测它与其他三个四元数之间的夹角,如果是钝角就翻转。
Angle、Euler、Matrix、Quaternion之间的相互转换
我们首先总结一下表示旋转的四种方法:
- Angle Axis(简称为Angle):任何一个旋转都可以表示为绕一个轴
旋转 度,这里的旋转轴 并不一定是坐标轴。 - Euler Angles(简称为Euler):任何一个旋转都可以表示为绕标准坐标系的旋转,旋转顺序可以自定。
- Matrix Multiplication(简称为Matrix):任何一个旋转都可以用矩阵乘法表示。
- Quaternion Rotation(简称为Quaternion):任何一个旋转都可以用四元数的乘法表示。
这几种表示方法实际上可以相互转化,下面将进行介绍。
Angle转化为其他形式
Angle转化为Euler
我们可以将四元数作为中间媒介进行转换,即先使用Angle转化为Quaternion将Axis Angle转化为四元数:
然后再通过Quaternion转化为Angle将四元数转化为Euler Angles:
其中
Angle转化为Matrix
在向量绕任意轴旋转的简单推导这篇文章中,我们已经推导出了AngleAxis对应的矩阵形式:
Angle转化为Quaternion
Angle与Quaternion在本质上是一致的,我们在三维旋转定理(四元数型)一节中已经证明过了。所以,对于轴为
Euler转化为其他形式
Euler转化为Angle
首先通过Euler转化为Matrix把欧拉角转化为对应的矩阵
然后再通过Matrix转化为Angle求解对应的(unnormalized)旋转轴
Euler转化为Matrix
在欧拉角、坐标系与旋转顺序一节中我们已经证明了:对于以
Euler转化为Quaternion
我们知道,Euler角的旋转是有顺序的,且每一次旋转都是绕标准坐标轴的旋转,所以,我们可以相应地构造出三个四元数去表示这些旋转(旋转顺序为
将这三个四元数相乘,即得到结果:
Matrix转化为其他形式
Matrix转化为Angle
给定旋转矩阵
从特征向量的角度来看,
从而我们有:
所以
既然我们已经知道了旋转轴
我们也可以从特征值的角度出发。在三维空间内,旋转矩阵的三个特征值分别为
Matrix转化为Euler
对于旋转矩阵
上式仅在
Matrix转化为Quaternion
给定矩阵
所以我们就能解出
从而得到四元数
Quaternion转化为其他形式
Quaternion转化为Angle
Quaternion转化为Axis Angle的结果我们已经在三维旋转定理(四元数型)一节中给出。给定单位四元数
得到它的等价形式
Quaternion转化为Euler
使用Quaternion转化为Matrix和Matrix转化为Euler的结果,我们能求出以
Quaternion转化为Matrix
结果已经由三维旋转定理(矩阵型)一节给出,结果是:
其中
https://krasjet.github.io/quaternion/quaternion.pdf↩︎
https://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles↩︎
http://www.euclideanspace.com/maths/geometry/rotations/conversions/index.htm↩︎
http://number-none.com/product/Hacking%20Quaternions/index.html↩︎
http://number-none.com/product/Understanding%20Slerp,%20Then%20Not%20Using%20It/↩︎
https://web.mit.edu/2.998/www/QuaternionReport1.pdf↩︎
https://www.geometrictools.com/Documentation/RotationRepresentations.pdf↩︎
https://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions#Rotation_matrix_%E2%86%94_Euler_axis/angle↩︎