第七章 四元数(有点技术性)

在本章中,我们将会讨论四元数,作为一种可选的旋转表达方式来替代旋转矩阵

我们主要用四元数来帮助我们在方位角之间进行较自然的插值。给空间中运动的物体做动画时会非常有用。如果在某个特定的程序中我们不打算对旋转进行插值,那么这种表示方式也许就不必要了。

7.1 插值

R_\alpha:=(R_1R_0^{-1})^\alpha R_0 \tag{7.1}

以及设置

7.1.1 周期循环(Cycles)

7.1.2 不变性(Invariance)

对于等式(7.1)所表示的单轴、恒定角速度有许多本质属性。首先,一个没有被施加外力的物体飞越空间,它的质心沿着直线前进,并且它的自转绕着固定的轴向。另外,这种类型的旋转插值同时满足左不变性和右不变性,我们现在就来解释一下。

图7.3:在此我们将坐标系换为蓝色的那个,但是使用新的物体坐标系向量继续在 time=0 和 time=1 绘制同一个正方形。如果一个插值满足右不变性,我们会得到对于该正方形相同的动画。

7.2 表示方式

四元数是四个实数的组合,我们很快会给它定义合适的操作。

我们用以下方式表示四元数

这个问题会在我们定义指数运算时让事情变得更复杂。

四元数

表示单位旋转矩阵。

四元数

任何这样形式的四元数

都有一个归一化的形式(四个元素的平方和的平方根为1)。相反地,任何单位归一化的四元数(和它的负)都可以被解释为一个唯一的旋转矩阵。

7.3 运算法则

四元数(不必须是归一化的)与标量相乘定义为

在两个四元数(不必须是归一化的)之间相乘,由以下看起来有点奇怪的操作定义

\left[ \begin{array}{ccc} w_1 \\ \mathbf{\hat{c}}_1 \end{array}\right] \left[ \begin{array}{ccc} w_2 \\ \mathbf{\hat{c}}_2 \end{array}\right] = \left[ \begin{array}{ccc} (w_1w_2-\mathbf{\hat{c}}_1 \cdot \mathbf{\hat{c}}_2) \\ (w_1\mathbf{\hat{c}}_2+w_2\mathbf{\hat{c}}_1+ \mathbf{\hat{c}}_1 \times \mathbf{\hat{c}}_2) \end{array}\right] \tag{7.2}

一个单位归一化四元数的乘法逆是

随后,我们将以下的三个四元数相乘

\left[ \begin{array}{ccc} \cos(\frac{\theta}{2}) \\ \sin(\frac{\theta}{2})\mathbf{\hat{k}} \end{array}\right] \left[ \begin{array}{ccc} 0 \\ \mathbf{\hat{c}} \end{array}\right] \left[ \begin{array}{ccc} \cos(\frac{\theta}{2}) \\ \sin(\frac{\theta}{2})\mathbf{\hat{k}} \end{array}\right]^{-1} \tag{7.3}

这个可以再次通过一系列不是很直观的计算来验证,这三个四元数相乘得到的结果实际上是一个以下形式的四元数

因此,四元数一方面显示的记录了旋转的轴向和角度,另一方面允许我们很容易的像矩阵一样来操作它们。

7.4 指数运算(Power)

给定一个单位四元数来表示旋转,我们可以用以下方式进行指数操作。 我们首先通过归一化四元数最后的三个元素来提取单位轴向$\mathbf{\hat{k}}$。 然后,我们使用atan2函数来提取角$\theta$。 这个函数返回给我们一个唯一的值$\theta/2 \in [-\pi \cdots \pi]$, 由此可以得到一个唯一的$\theta \in [-2\pi \cdots 2\pi]$。然后我们定义

当$\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$与世界坐标系相关,并且如果这两个矩阵可以转换为以下两个四元数形式

然后我们只需要简单的计算四元数: \Bigg(\left[ \begin{array}{ccc} \cos(\frac{\theta_1}{2}) \\ \sin(\frac{\theta_1}{2})\mathbf{\hat{k}} \end{array}\right] \left[ \begin{array}{ccc} \cos(\frac{\theta_0}{2}) \\ \sin(\frac{\theta_0}{2})\mathbf{\hat{k}} \end{array}\right]\Bigg)^\alpha \left[ \begin{array}{ccc} \cos(\frac{\theta_0}{2}) \\ \sin(\frac{\theta_0}{2})\mathbf{\hat{k}} \end{array}\right] \tag{7.4}

\frac{\sin[(1-\alpha)\Omega]}{\sin(\Omega)} \vec{v_0} + \frac{sin(\alpha\Omega)}{sin(\Omega)} \vec{v_1} \tag{7.5}

\frac{\sin[(1-\alpha)\Omega]}{\sin(\Omega)} \left[ \begin{array}{ccc} \cos(\frac{\theta_0}{2}) \\ \sin(\frac{\theta_0}{2})\mathbf{\hat{k}}_0 \end{array}\right] + \frac{sin(\alpha\Omega)}{sin(\Omega)} \left[ \begin{array}{ccc} \cos(\frac{\theta_1}{2}) \\ \sin(\frac{\theta_1}{2})\mathbf{\hat{k}}_1 \end{array}\right] \tag{7.6}

图7.4:

线性插值运算同时具有左不变性和右不变性。例如,左不变性为

基于指数和基于球体的球面插值的等效性(选修)

这里我们描述了建立基于指数和基于球体的旋转插值的等效性时所需要的步骤。

\left[ \begin{array}{ccc} \cos(\frac{\alpha\theta_1}{2}) \\ \sin(\frac{\alpha\theta_1}{2})\mathbf{\hat{k}}_1 \end{array}\right] \tag{7.7}

  • 一个三角学的参数可以用来表明等式(7.6)在几何上与沿着球体表面的插值一致。

7.5 代码实现

很容易就可以写出一个四元数的类。

我们还需要写代码来实现 MakeRotation 函数,和我们在矩阵中所做的相同。

给定一个单位四元数q和一个实数alpha,我们定义指数运算符:pow(q,alpha)。 给定两个四元数 q0q1,我们可以定义四元数的插值运算:slerp(q0,q1,alpha)。 记得在实现 slerp 时,如果它的第一个元素为负, 我们需要将四元数求负(q1 * inv(q0))来得到在“短侧”的插值。

7.6 重新将位移考虑进来

目前为止,我们已经讨论了在表示旋转时四元数的作用,但我们忽略了位移。 现在我们将会讨论怎样将四元数和位移向量一起使用,来表示刚体变换。

一个刚体变换,即RBT,可以用位移和旋转的组合来描述

\begin{eqnarray} A &=& TR \\ \left[ \begin{array}{ccc} r & t \\ 0 & 1 \end{array}\right] &=& \left[ \begin{array}{ccc} i & t \\ 0 & 1 \end{array}\right] \left[ \begin{array}{ccc} r & 0 \\ 0 & 1 \end{array}\right] \end{eqnarray}

由此,我们可以这样来表示一个物体:

class RigTform{
		Cvec4 t;
		Quat r;
	};

记得注意,因为 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 来表示 eyeRBTobjRBT

想要创建一个有用的刚体变换,我们需要以下操作:

RigTForm identity();
RigTForm makeXRotation(const double ang);
RigTForm makeYRotation(const double ang);
RigTForm makeYRotation(const double ang);
RigTForm makeTranslation(const double tx, const double ty,const double tz);

我们需要实现 RigTForm ACvec4 c 的相乘,它返回 A.r * c + A.t

随后,我们需要两个 RigTForm 相乘。为了理解如何实现这一点,我们来看看两个这样的刚体变换相乘。

随后,我们需要实现这个数据类型的逆。如果尝试求一个刚体变换的逆,我们可以看到

给定了这一结构,我们现在可以使用我们新的数据类型来重写函数 doMtoOwrtA(RigTForm M, RigTForm O, RigTForm A)

练习

7.1 在本书的网站上,我们有 Quat 类的部分实现。使用这部分代码来创建一个刚体变换类。使用这个类替代 Matrix4 类来表示在练习 6.6 中你的 OpenGL 代码的刚体变换部分。

7.2 为 Quat 类实现指数运算,并通过一个立方体在不同方向上旋转的线性插值来测试它。

Last updated