四元数旋转

四元数旋转

四元数在3D几何中有非常重要的应用——插值旋转。但它在图形学里可能算是比较难啃的一根骨头：首先我们上学的时候基本不会学到四元数的原理，其次我们活在一个3D的世界里，很难想象4D空间在发生着什么。理解四元数旋转还是挺困难的一件事。

——先几句题外话。本文中笔者按照数学的习惯采用右手系，X轴向右、Y轴向上时，Z轴朝向屏幕外，旋转角朝逆时针方向增长；矩阵采用列主序，左乘矩阵进行变换。左手系主要是为了让Z轴指向屏幕内，令顺时针方向为叉积的正方向；行主序主要是为了通过右乘矩阵来进行变换，让变换按从左到右的顺序进行，某种意义上更符合“直觉”。GL系的图形API（如OpenGL和Vulkan）都默认遵守这一规则，只有老版本的DX系API使用左手系和行主序。但是令人比较头痛的是，直到现在新版本DX乃至一票游戏引擎的默认手系仍是左手系，这对不熟悉左手系的人非常不友好，在将数学工具翻译到代码时常常会遇到问题。

要讨论旋转的计算，就要先理解旋转究竟是什么。这里就不扯太抽象的东西了，简单定义：旋转是一种变换，使得变换前后空间内任意两点距离不变。具体到其中一种情况，旋转前后任意一点到旋转中心的距离不变。

我们非常熟悉2D空间中的旋转——这几乎是每一本线性代数教材必讲的内容。所有的2D旋转都可以转写为旋转矩阵：

\[M_\theta = \begin{bmatrix} \cos(\theta)&-\sin(\theta)\\ \sin(\theta)&\cos(\theta)\\ \end{bmatrix} \]

其中 \(\theta\) 为逆时针旋转的角度。第一列即 \(\begin{pmatrix}1&0\end{pmatrix}^T\) 旋转后的位置，第二列即 \(\begin{pmatrix}0&1\end{pmatrix}^T\) 旋转后的位置。但不论是什么计算，在被用矩阵写出来的那一刻起，我们就很难看清它的本质。又因为二维旋转过于特殊，我们退一步，不着急思考旋转操作的几何意义，只在数值集合的意义上做讨论。

首先我们来讨论一些比较简单的情况。对一个自然数 \(x\in\mathbb{N}\) 绕0进行旋转操作，整个集合中 \(x\) 唯一的落点就是其自身；而对于整数 \(y\in\mathbb{Z}\) ，却有 \(y\) 和 \(-y\) 两个不同的落点。

观察这两种情况。显然自然数集不存在比0更小的数，从0到一个非0值只有一个“方向”；而整数集含有一半小于0的负数，对于每一 \(y\) 总有一 \(-y\) 与其有相同的幅度，从0到非0值有两个“方向”可以行进。非常合理的，在向量长度固定的情况下，有多少个“方向”，就有多少个不同的取值。也就是说，对于N维向量而言，其所有可能旋转操作的落点，位于该向量所在的一个N-1维球面上（因为半径确定而-1）。

那么，如果我们给定旋转的起点和终点，想要求这两点中最短的旋转轨迹，基本就是在寻找这个球面上的一条弧。但是我们平时操作的向量都存在于线性空间里，在线性空间中我们是很难求出一条非线性轨迹的。所以我们要主动进入非线性空间，找到这条弧，然后再把这条弧上的点逐个映射回线性空间中来使用。

嗯？我在说什么？我们现在来举个例子：对于二维空间中的向量 \(\vec r\) 绕原点进行旋转，所有落点均在以原点为中心的圆上。如果我们直接对旋转的输入输出进行线性插值，我们算出来的“旋转轨迹”是一条直线，但是很显然它应该是沿着那个圆的切向行进的。

图1 沿蓝色圆上的一条弧旋转，对该过程进行插值。红色为线性插值得到的轨迹；绿色为合理的旋转轨迹。

为什么插值出来的路径会是一条直线呢？这和线性插值的原理有关：假设我们有 \(\vec r_0,\vec r_1\) 两个已知数据，求两端点之间 \(\alpha\) 比例位置我们所不知道的点 \(\vec r_\alpha\) ，线性插值通过计算 \((1 - \alpha)\vec r_0 + (\alpha)\vec r_1\) 得到这个点坐标的一个近似。如果我们重整这条式子，会得到一条直线的表达式

\[\vec r_\alpha = \vec r_0 + \alpha(\vec r_1 - \vec r_0) \]

既然线性插值不能达到我们所期望的效果，如果我们不使用线性插值，改用别的插值方法，是否可以呢？当然可以。我们有一个极其优雅的公式可以帮助我们从线性空间迁移到非线性空间中，得到这条合理的轨迹，即著名的欧拉公式：

\[e^{ix} = \cos x + i\sin x \]

把我们的二维平面看作复平面，我们的X轴变成了实轴，Y轴变成了虚轴。空间中任意一个向量 \(\vec r\) 都可以用 \(c\) 这个复数来转述：

\[c=Ae^{i\phi} \]

其中

\[\begin{aligned} A&=\|\vec r\|\\ \phi&=\arctan(\frac{r_y}{r_x})\\ \end{aligned}\]

随后若要将这个复数在复平面中逆时针旋转 \(\theta\) 角，可以通过复数乘法完成：

\[\begin{aligned} \mathrm{rotate}(c, \theta)&=Ae^{i\phi}e^{i\theta}\\ &=Ae^{i(\phi + \theta)}\\ &=A\cos(\phi + \theta)+Ai\sin(\phi + \theta)\\ \end{aligned}\]

可能有同学会发现了，这就是极坐标嘛。显然，向量到原点的距离 \(A\) 不变，我们对 \(\phi\) 进行插值能得到一条完美的弧线，而且叠加旋转不涉及复杂的三角学，这在数字计算中无疑是非常有优势的。

现在我们来讨论三维空间中的旋转。传统上，3D旋转通过欧拉角完成。欧拉角涉及三次旋转，分别绕X、Y、Z三轴进行，矩阵形式如下：

\[M=M_xM_yM_z \]

其中

\[\begin{aligned} M_x &= \begin{bmatrix}1&0&0\\0&\cos\psi&-\sin\psi\\0&\sin\psi&\cos\psi\end{bmatrix}\\ M_y &= \begin{bmatrix}\cos\theta&0&\sin\theta\\0&1&0\\-\sin\theta&0&\cos\theta\end{bmatrix}\\ M_z &= \begin{bmatrix}\cos\phi&-\sin\phi&0\\\sin\phi&\cos\phi&0\\0&0&1\end{bmatrix} \end{aligned}\]

非常有趣的是，欧拉角并没有规定旋转轴必须分别是X、Y、Z轴，也没有规定矩阵相乘的顺序，所以你可能会在不同的文献里看到各种各样的欧拉角。

欧拉角在图形学场景下不是很好用。首先，我们很难将一个旋转操作分解成三个独立的操作，欧拉角更适合应用在飞机操控那类本来就有三个控制变量的场景；其次，欧拉角存在一个叫做万向锁的问题，在第二次旋转角恰为90°时，旋转整体会丢失一个自由度；再者，计算三个旋转矩阵，再逐个相乘，要是有叠加的旋转操作还要做更多的乘法计算，这对实时图形来说太耗时了…那么有什么好的方法可以替代欧拉角呢？要是三维空间有什么东西和复数一样好用就好了。

于是学者们就将二元复数扩展成了四元数，虚部从一个量 \(i\) 扩增为三元的向量 \(\vec u = (i,j,k)^{\mathrm T}\) 。和复数相似的，欧拉公式对四元数也成立：

\[e^{\vec ux} = \cos x + \vec u\sin x \]

也就有

\[q = Ae^{\vec u\phi} \]

通常标记四元数时，四元数的实部叫做 \(w\) ，而其虚部的三个成分分别为 \(x,y,z\) ，虚部整体作为一个向量用 \(\vec u\) 标记 。

慢着，为什么我们旋转一个三维向量需要用到第四个维度？这其实和球面的性质有关系。和前面提到的一样，我们N维空间中的球面只有N-1维，也就是说，三维空间中的球面，只有两个维度。听上去有些反直觉，但是只要想想我们GPS用的经纬坐标系马上就能理解。因为半径位置已经固定了，我们只需要两个坐标就可以定位球面上的一个点，但如果我们想用向量的三个坐标来定位一个球面上的点，我们就必须进入四维空间，然后“废掉”一个维度，原来的球面就变成了3-球面（3-sphere）。为了计算的简洁，这个被“废掉”的维度就是我们的四元数实部。我们称这种实部为0的四元数为纯四元数（Pure quaternion）。

相比复数乘法来说，四元数乘法有一些特殊的性质。比如说，四元数乘法是反交换的，也就是说 \(qp=-pq\) ，因此我们不可以随意交换乘数的顺序，必须按照从左到右的顺序计算。此外，四元数乘法的计算方式也更奇怪一些。给定 \(q=a+\vec u\) 和 \(p=t+\vec v\) ，有

\[qp = at+a\vec v+t\vec u-\vec u\cdot\vec v+\vec v\times\vec u \]

这被称作格拉斯曼积。如果我们将其和复数 \(z=x+ui\) 和 \(w=y+vi\) 间的乘法比较：

\[zw = xu + xvi + yui - yv \]

我们会发现四元数的乘法多出来了一个项，四元数虚数的虚数乘积部分从一项变成了两项。要理解这一点，我们可以将这虚数乘积部分展开：

\[\begin{aligned} \vec u\vec v=&(ai+bj+ck)(xi+yj+zk)\\ =&ax(i^2)+ay(ij)+az(ik)\\ &+bx(ji)+by(j^2)+bz(jk)\\ &+cx(ki)+cy(kj)+cz(k^2)\\ =&\underbrace{(bz-cy)i+(cx-az)j+(ay-bx)k}\_{\vec u\times\vec v}-\underbrace{(ax+by+cz)}_{\vec u\cdot\vec v} \end{aligned}\]

如果四元数乘数虚部的相同两个分量为0，四元数行为上就会退化成虚数，这里限于篇幅不作演示。对了，说出来你可能不信，向量叉积的定义和 \(R^3\) 中使用 \(i,j,k\) 标记三个轴向单位的习惯都是从四元数拿来的。

既然四元数是复数的扩展，那么四元数旋转应当也能通过乘法完成。假设我们有旋转轴 \((x,y,z)^\mathrm T\) 和旋转角 \(\phi\) ，我们可以通过欧拉公式得到旋转四元数

\[q=e^{(xi+yj+zk)\phi}=\cos\phi+\underbrace{(xi+yj+zk)}_\vec u\sin\phi \]

与一个表示位置的纯四元数 \(p=\underbrace{ai+bj+ck}_\vec v\) 相乘，有

\[\begin{aligned} p'&=qp\\ &=(\cos\phi+\vec u\sin\phi)\vec v\\ &=\vec v\cos\phi\vec+(\vec u\times\vec v)\sin\phi-(\vec u\cdot\vec v)\sin\phi \end{aligned}\]

注意，我们操作后的位置 \(p'\) 应当也是一个纯四元数，要不然就脱离我们的三维子空间了；但上面这个式子只有在 \(\phi\) 为180°的整数倍时实部才为0。为了解决这个问题，前人想到了一个非常巧妙的方法：将旋转分两次进行。

重新定义

\[q=\cos\frac\phi 2+\vec u\sin\frac\phi 2 \]

有

\[\begin{aligned} p'&=qpq^*\\ &=(\cos\frac\phi 2+\vec u\sin\frac\phi 2)\vec v(\cos\frac\phi 2-\vec u\sin\frac\phi 2) \end{aligned}\]

其中 \(q^*\) 是 \(q\) 的共轭，四元数的共轭和复数的相似，都可以通过改变虚部符号得到。先让向量绕轴正方向逆时针旋转一半的角度，再让它绕轴反方向顺时针旋转一半的角度。经过两次四元数乘法，原本被带入实部的成分又被带回了虚部。

图2 利用四元数乘法进行180°旋转的两个步骤。

关于四元数旋转等距的相关证明可以在这里找到，此处不再赘述。