第七章 四元数(有点技术性)
在本章中,我们将会讨论四元数,作为一种可选的旋转表达方式来替代旋转矩阵
R=[r001],
我们主要用四元数来帮助我们在方位角之间进行较自然的插值。给空间中运动的物体做动画时会非常有用。如果在某个特定的程序中我们不打算对旋转进行插值,那么这种表示方式也许就不必要了。
7.1 插值
假设我们已有一个物体坐标系 “time=0,” o0t=wtR0,和另一个物体坐标系“time=1",o1t=wtR1,其中R0和R1是 4x4 旋转矩阵。假设我们希望找到一个坐标系序列oαt,其中α∈[0...1],得到一个从o0t到o1t的自然旋转。
一个想法是定义Rα:=(1−α)R0+(α)R1并且设置otα=wtRα。这样做的问题在于,在对矩阵线性插值时,每个向量只是简单的沿着直线移动(如图 7.1)。此时,过渡的Rα并非旋转矩阵,这明显是不可接受的。 此外(尤其是在 3D 中),想要设法找到一个与这一插值矩阵近似,并且能够抵消掉这种挤压效果的真正的旋转矩阵,也并不容易。
图7.1:对于 Ri 矩阵进行线性插值,只能起到让向量沿着直线方向移动的效果,没有旋转运动
另一个不值得深究想法,是用某种方法将R0和R1提取为欧拉角的形式(如图 2.5)。可以用α对得到的三个标量进行线性插值。插值后的欧拉角可以用来对旋转进行过渡。这样做的结果是,得到的插值旋转序列无论在物理上还是几何上都不自然。例如,当基向量变化时它不具有不变性(invariant)(不变性在随后被定义)。
我们真正想做的是,首先创建一个过渡矩阵R1R0−1。这个矩阵能够和任何旋转矩阵一样,被看做是沿着某个轴向 [kx,ky,kz]t旋转了某个角度θ。我们现在想象一下我们有一个操作 (R1R0−1)α可以给我们一个旋转。然后想要定义
以及设置
oαt=wαtRα
就会很自然了。这样做我们可以得到一个坐标系序列,沿着某个轴向持续旋转,我们不断增加α,直到完成后续姿态。很明显,开始的时候有wt(R1R0−1)0R0=wtR0=o0t, 并且最终有 wt(R1R0−1)1R0=wtR1=o1t, 正如我们所期望的一样。
这里最关键的部分是以轴/角对的方式提取出过渡的旋转矩阵,例如R1R0−1。四元数的理念是,始终跟踪记录轴和角的数据,因而就无需再进行这样的提取了。重要的是,以四元数表示旋转时,使用旋转矩阵时需要用到的必须操作,我们仍然能够继续使用。
7.1.1 周期循环(Cycles)
在这里我们需要澄清一个细节。矩阵R1R0−1实际上,可以被看做θ+n2π度的旋转,n可以为任何整数。将这个旋转的影响看做对一个向量进行线性变换时,这个额外的2π是无关的。但当定义一个指数操作符给我们“一个绕着单一轴逐渐旋转的坐标系序列,直到完成后续姿态”时,我们需要决定怎样通过给定的矩阵R1R0−1选择n的值。对插值来说,自然的选择是选择当 ∣θ+n2π∣ 值最小时的n。即意味着,对于这个n来说,θ+n2π∈[−π...π]。另外,这个选择是没有歧义的(unambiguous)(除了 θ=π+n2π时,此时我们需要在−π和π之间随意取值)。实际上,这个关于n的选择在随后的 7.4 章中会提到。
7.1.2 不变性(Invariance)
对于等式(7.1)所表示的单轴、恒定角速度有许多本质属性。首先,一个没有被施加外力的物体飞越空间,它的质心沿着直线前进,并且它的自转绕着固定的轴向。另外,这种类型的旋转插值同时满足左不变性和右不变性,我们现在就来解释一下。
假设我们有一个变化的世界坐标系w′t例如wt=w′tRl,其中Rl是某个固定的“左”旋转矩阵。使用这个坐标系,我们可以重新表示我们原先的插值问题为在 o0t:=w′tRlR0 和 o1t:=w′tRlR1 之间插值。左不变性(left Invariance)就是说,如果这个表示方式的更换没有改变我们插值的结果:如果原始的插值是 otα=wtRα, 而新方法表示的插值是 w′tRlRα=wtRα,即我们得到了刚好同样的 otα。 换句话说,如果一个插值只取决于几何体上的物体坐标系 o0t和o1t, 而与世界坐标系和R0以及R1的结果无关,那么插值过程具有左不变性。左不变性是一个十分自然的属性;只有很少的情况下插值与世界坐标系的选择是有关的。
图7.2:当我们以 w′t替代wt时,如果一个插值不因此发生变化,则该插值就具有左不变性。(本图中的操作只涉及旋转,其中图像的位移是为了更方便我们观察它。)
我们可以看出等式(7.1)满足左不变性如下。仿射变换的“过渡”, 将坐标系o0t映射到另外一个坐标系o1t上,过程总是唯一的。 在本例中,因为我们的坐标系是右手系,正交,并且共用一个原点,这一变换肯定是旋转变换。同时,我们刚刚定义的指数操作符(也就是,保持轴向不变的同时变换角度)是一个隐性的几何体操作并且可以在不参照任何坐标系的情况下被描述。因此,这个插值过程不依赖于世界坐标系的选择并且具有左不变性。
相比之下,右不变性,表示即使我们改变所使用的物体坐标系,该物体的插值也会发生改变。 比如说,假设我们有一个“右”旋转矩阵Rr,并且用它来定义新的物体坐标系,在 time=0 和 time=1 时分别为 o0tRr和o1tRr。 因为我们不希望物体本身的姿态发生改变,我们使用c′=Rr−1c 为物体上的所有点重新指定合适的物体坐标系。如果对物体坐标系的这一改变,对于旋转物体自身的插值没有影响,我们就说这个插值满足右不变性。换句话说,如果原始的插值给我们Rα 然后新的插值给我们 oα=wtRα。 可以看出,物体的插值(使用新的物体坐标系c′)没有变化。右不变性是一个十分自然的属性。但是我们将会看到当我们将位移包括进来时,我们也许希望插值与物体坐标系的原点有关。
图7.3:在此我们将坐标系换为蓝色的那个,但是使用新的物体坐标系向量继续在 time=0 和 time=1 绘制同一个正方形。如果一个插值满足右不变性,我们会得到对于该正方形相同的动画。
我们可以直接看出等式(7.1)表示的插值满足右不变性如下 [(R1Rr−1)(RrR0)]αR0Rr=(R1R0−1)αR0Rr=RαRr
7.2 表示方式
四元数是四个实数的组合,我们很快会给它定义合适的操作。
我们用以下方式表示四元数
[wc^],
其中w是一个标量,而c^是一个三维的坐标系向量。我们加一个 “∧” 符号来标记它,以便与四维坐标系向量区分。
为了表示一个绕着单位长度的轴向k^旋转θ度,我们使用四元数
[cos(2θ)sin(2θ)k^],
角度需要除以2看起来有点奇怪,但它能够让四元数旋转正常工作,我们随后会描述到。注意一个绕着 −k^轴向旋转−θ角度给我们同样的四元数。 一个绕着轴向k^旋转θ+4π度也会给我们同样的四元数。 目前为止一切都很好。奇怪的是,绕着轴向k^旋转θ+2π度, 它实际上是同样的旋转,却会给我们一个负的四元数
[−cos(2θ)−sin(2θ)k^],
这个问题会在我们定义指数运算时让事情变得更复杂。
四元数
[10^],[−10^]
表示单位旋转矩阵。
四元数
[1k^],[−1k^]
表示沿着k^轴旋转180∘。
任何这样形式的四元数
[cos(2θ)sin(2θ)k^]
都有一个归一化的形式(四个元素的平方和的平方根为1)。相反地,任何单位归一化的四元数(和它的负)都可以被解释为一个唯一的旋转矩阵。
7.3 运算法则
四元数(不必须是归一化的)与标量相乘定义为
α[wc^]=[αwαc^]
在两个四元数(不必须是归一化的)之间相乘,由以下看起来有点奇怪的操作定义
其中⋅和×是 3D 向量的点乘和叉乘。这个奇怪的相乘拥有以下有用的属性:如果[wi,c^i]t表示旋转矩阵Ri, 而乘法[w1,c^1]t[w2,c^2]t表示旋转矩阵R1R2。 这个属性可以通过一系列不是很直观的计算来验证。
一个单位归一化四元数的乘法逆是
[cos(2θ)sin(2θ)k^]−1=[cos(2θ)−sin(2θ)k^]
这个四元数简单的绕着同样的轴旋转−θ∘度。(逆也可以由非单位四元数定义,但我们不需要用到这个。)
重要的是,我们可以用四元数相乘来对一个坐标系向量实施一个旋转。假设我们有一个四维的坐标系向量 c=[c^,1]t,并且我们用一个 4x4 的旋转矩阵左乘它来得到
c′=Rc,
在此返回给我们的四维坐标系向量是 c′=[c^′,1]t。想用四元数来实现它, 需要先让R以单位归一化的四元数形式来表示为
[cos(2θ)sin(2θ)k^]
我们用三维坐标系向量c^来创造一个非单位四元数
[0c^]
随后,我们将以下的三个四元数相乘
这个可以再次通过一系列不是很直观的计算来验证,这三个四元数相乘得到的结果实际上是一个以下形式的四元数
[0c^′],
其中c^′是我们想要得到的三维坐标系向量。
因此,四元数一方面显示的记录了旋转的轴向和角度,另一方面允许我们很容易的像矩阵一样来操作它们。
7.4 指数运算(Power)
给定一个单位四元数来表示旋转,我们可以用以下方式进行指数操作。 我们首先通过归一化四元数最后的三个元素来提取单位轴向$\mathbf{\hat{k}}$。 然后,我们使用atan2函数来提取角$\theta$。 这个函数返回给我们一个唯一的值$\theta/2 \in [-\pi \cdots \pi]$, 由此可以得到一个唯一的$\theta \in [-2\pi \cdots 2\pi]$。然后我们定义
[cos(2θ)sin(2θ)k^]α=[cos(2αθ)sin(2αθ)k^]
当$\alpha$从0变化到1时,我们得到在0和$\theta$角之间一系列的旋转。 如果$\cos(\frac{\theta}{2})>0$,我们有$\theta/2 \in [-\pi/2 \cdots \pi/2]$, 且因此$\theta \in [-\pi \cdots \pi]$。 此时,我们使用$\alpha \in [0 \cdots 1]$来在两个旋转间插值, 我们将会在旋转较短的一侧进行插值。 相反地,如果$\cos(\frac{\theta}{2})<0$, 那么$|\theta| \in [-\pi \cdots \pi]$,而我们将会在“较长的一侧”(大于$180^{\circ}$)进行插值。 通常情况下,在较短的一侧进行旋转插值更加自然,并且当一个给定的四元数第一个元素为负值时,我们总是在进行指数运算前将四元数变负。
7.4.1 球面插值与线性插值(Slerp and Lerp)
将所有这些放在一起,如果我们想要在两个坐标系间进行插值,这两个坐标系通过旋转矩阵$R_0$和$R_1$与世界坐标系相关,并且如果这两个矩阵可以转换为以下两个四元数形式
[cos(2θ0)sin(2θ0)k^][cos(2θ1)sin(2θ1)k^]
然后我们只需要简单的计算四元数:
这个插值操作经常被称为球面线性插值,原因如下。单位四元数的四个实数的平方和为1;因此,我们可以将它看做单位四维球体R4上的一个几何体点。如果从两个单位四元数开始并且对它们进行插值,就可以从公式(7.4)看出,在球体R4上得到的轨迹,实际上和将两个单位球上的点连接起来的大圆弧完全相同。而且,在该轨迹上的插值出的弧长与α成比例。
在任何维度n下,都可以用一个三角参数来表示在任何两个单位向量, v0和v1间的球面线性插值,计算如下
其中Ω是R4内向量间的角度。因此,我们可以将等式(7.4)替换为
其中Ω是R4四元数中起始和最终的角度。稍后,我们会讲到验证等式(7.4)和(7.6)所需要的步骤。注意在等式(7.4)中,为了选择“短侧”小于180∘的插值,如果两个四元数之间四维点乘的结果为负,我们必须将两个四元数之一变负。
这样看来,我们可以用R4中更简单的线性插值去近似等式(7.4),以及等式(7.6),由此,我们可以简单的计算得到
(1−α)[cos(2θ0)sin(2θ0)k^0]+(α)[cos(2θ1)sin(2θ1)k^1]
因为这个插值已经不再是一个单位四元数,我们必须将结果归一化,不过这个很容易。重要的是,这个插值过程,称作线性插值(Lerping),可以得到和更复杂的球面插值(Slerp)同样的轨迹(沿着一个固定的单轴旋转),尽管它在α之间旋转的角度不再是均等的(如图 7.4)。
图7.4:
线性插值运算同时具有左不变性和右不变性。例如,左不变性为
(1−α)[wlc^l][w0c^0]+(α)[wlc^l][w1c^1]
=[wlc^l]((1−α)[w0c^0]+(α)[w1c^1])
标量相乘可以在四元数乘法之间互换,并且四元数相乘可以提取到和之上。同样,我们可以直接从等式(7.6)看出它也同时具有左不变性和右不变性。计算中唯一有技巧的部分展示了角度Ω,同时具有左不变性和右不变性。
线性插值可以比球面插值更有效率的执行。更重要的是,它可以让在n个不同旋转之间的混合概括的更容易。我们可以轻易的混合R4中的四元数,并随后将它归一化。这种 n-way 间的混合在构建一系列的旋转曲线(见 section 9.3)以及当为动画蒙皮时(见 section 23.1.2)会很有用。还有一些其他的方法能够在球体内部实现 n-way 混合,尽管这种方法不是一个固定时间的计算方法[10]。
基于指数和基于球体的球面插值的等效性(选修)
这里我们描述了建立基于指数和基于球体的旋转插值的等效性时所需要的步骤。
应该说(正如我们之前做的一样)等式(7.6)具有左不变性。注意我们之前已经为等式(7.4)建立了左不变性。因为两种插值方式都具有左不变性,我们可以在不损失一般性的情况下,只考虑R0是单位矩阵时的情况。
假设R0是单位矩阵,基于指数的插值等式(7.4)给我们(R1)α,即
因为R0是单位矩阵,初始的四元数是[1,0^]t。将其代入到等式(7.6),我们可以验证它也符合等式(7.7)。
一个三角学的参数可以用来表明等式(7.6)在几何上与沿着球体表面的插值一致。
7.5 代码实现
很容易就可以写出一个四元数的类。
我们定义Quat为一个四组的实数。随后我们根据等式(7.2)来定义乘法(q1 * q2)。 给定一个单位四元数,Quat q,我们定义它的逆 inv(q)。 给定一个单位四元数,Quat q,和一个 Cvec4 c,我们定义(q * c), 根据等式(7.3),将旋转应用于坐标系向量c,并返回坐标系向量c′。
我们还需要写代码来实现 MakeRotation 函数,和我们在矩阵中所做的相同。
给定一个单位四元数q和一个实数alpha,我们定义指数运算符:pow(q,alpha)。 给定两个四元数 q0 和q1,我们可以定义四元数的插值运算:slerp(q0,q1,alpha)。 记得在实现 slerp 时,如果它的第一个元素为负, 我们需要将四元数求负(q1 * inv(q0))来得到在“短侧”的插值。
7.6 重新将位移考虑进来
目前为止,我们已经讨论了在表示旋转时四元数的作用,但我们忽略了位移。 现在我们将会讨论怎样将四元数和位移向量一起使用,来表示刚体变换。
一个刚体变换,即RBT,可以用位移和旋转的组合来描述
由此,我们可以这样来表示一个物体:
记得注意,因为 t 表示一个位移向量,它的第四个元素为 0。
7.6.1 插值(Interpolation)
给定两个RBTs$O0=(O_0)_T(O_0)_R$和$O_1=(O_1)_T(O_1)_R$,我们可以在它们之间插值,首先在表示位移的两个位移向量间插值来得到位移$T\alpha$,随后在表示旋转的两个四元数间进行球面插值得到旋转$R\alpha$,并最终设置插值得到的RBT$O\alpha$为$T\alpha R\alpha$。如果 $\vec{\mathbf{o}}0^t=\vec{\mathbf{w}}^tO_0$且 $\vec{\mathbf{o}}_1^t=\vec{\mathbf{w}}^tO_1$。我们随后可以设置 $\vec{\mathbf{o}}\alpha^t=\vec{\mathbf{w}}^tO_\alpha$。 在这种插值下,$\vec{\mathbf{o}}^t$的原点以固定速度沿着直线运动, 而$\vec{\mathbf{o}}^t$的向量基以固定的角速度沿着固定轴向旋转。正如我们已经提到过的,这样的运动十分自然,就像一个没有施加外力的物体穿过空间,并且它的质心沿着直线运动,并且它绕着固定的轴旋转。另外,这种在$\vec{\mathbf{o}}_0^t$和$\vec{\mathbf{o}}_1^t$之间唯一的几何插值可以在任意坐标系下表达。因此,它不依赖于世界坐标系的选择而一定具有左不变性。
需要注意到RBT插值并不具有右不变性。如果你改变物体坐标系并且用这个方法在它们之间插值,得到的新原点会沿着直线运动,但原先的原点会走过一个曲线。因此,这个方法在被插值的物体有明确含义的“中心”时最有意义。在这样的中心不存在时,最自然的答案就不是很明显了(例如,参见[35])。
7.6.2 运算
回到我们 6.2 节的绘图代码,我们现在可以用 RigTform 数据类型替代 Matrix4 来表示 eyeRBT 和 objRBT。
想要创建一个有用的刚体变换,我们需要以下操作:
我们需要实现 RigTForm A 和 Cvec4 c 的相乘,它返回 A.r * c + A.t 。
随后,我们需要两个 RigTForm 相乘。为了理解如何实现这一点,我们来看看两个这样的刚体变换相乘。
[i0t11][r1001][i0t21][r2001]
=[i0t11][r10r1t21][r2001]
=[i0t11][i0r1t21][r1001][r2001]
=[i0t1+r1t21][r1r2001]
从中我们可以看到,得到的结果是一个新的刚体变换,其中位移为 t1+r1t2,而旋转为 r1r2。
随后,我们需要实现这个数据类型的逆。如果尝试求一个刚体变换的逆,我们可以看到
([i0t1][r001])−1
=[r001]−1[i0t1]−1
=[r−1001][i0t−11]
=[r−10−(r−1t)1]
=[i0−(r−1t)1][i001]
由此,我们可以看到,得到的结果是一个新的刚体变换,其中位移为 −(r−1t),而旋转为 r−1。
给定了这一结构,我们现在可以使用我们新的数据类型来重写函数 doMtoOwrtA(RigTForm M, RigTForm O, RigTForm A) 。
练习
7.1 在本书的网站上,我们有 Quat 类的部分实现。使用这部分代码来创建一个刚体变换类。使用这个类替代 Matrix4 类来表示在练习 6.6 中你的 OpenGL 代码的刚体变换部分。
7.2 为 Quat 类实现指数运算,并通过一个立方体在不同方向上旋转的线性插值来测试它。
Last updated