第4章-变换-4.3-四元数

4.3 四元数

尽管四元数早在1843年就由William Rowan Hamilton爵士发明，作为复数的扩展，但直到1985年Shoemake[1633]才将它们引入计算机图形领域\(^1\)。四元数用于表示旋转和方位。它们在几个方面都优于欧拉角和矩阵。任何三维方向都可以表示为围绕特定轴的单次旋转。给定轴和角度表示，与四元数转换相互转换很简单，然后任一方向的欧拉角转换则具有挑战性。四元数可用于稳定和恒定的方向插值，这是欧拉角无法很好完成的。

复数有实部和虚部。每个复数由两个实数表示，第二个实数乘以\(\sqrt{-1}\)。同样，四元数有四个部分。前三个值与旋转轴密切相关，旋转角度影响所有四个部分（更多可参考第4.3.2节）。每个四元数由四个实数表示，每个实数与不同的部分相关联。由于四元数有四个分量，我们选择将它们表示为向量，但为了区分它们，我们给它们戴上帽子：\(\hat{\rm\pmb{q}}\)。我们从四元数的一些数学背景开始，然后用它来构造各种有用的转换。

4.3.1 数学背景

我们从四元数的定义开始。

定义：四元数\(\hat{\rm\pmb{q}}\)可以通过以下方式定义，都是等价的。

\[\begin{split} \hat{\rm\pmb{q}} = ({\rm\pmb{q}}_v, q_w) = iq_x + jq_y + kq_z + q_w = {\rm\pmb{q}}_v + q_w ,\\ {\rm\pmb{q}}_v = iq_x + jq_y + kq_z = (q_x, q_y, q_z),\\ i^2 = j^2 = k^2 = −1, jk = −kj = i, ki = −ik = j, ij = −ji = k. \end{split} \tag{4.31} \]

变量\(q_w\)称为四元数\(\hat{\rm\pmb{q}}\)的实部。虚部是\({\rm\pmb{q}}_v\)，而\(i\)、\(j\)和\(k\)称为虚数单位。

对于虚部\({\rm\pmb{q}}_v\)，我们可以使用所有法向量运算，例如加法、缩放、点积、叉积等。使用四元数的定义，推导出两个四元数\(\hat{\rm\pmb{q}}\)和\(\hat{\rm\pmb{r}}\)之间的乘法运算，如下所示。请注意，虚数单位的乘法是不可交换的。

\[\begin{equation} \begin{aligned} \pmb{乘法:} \hat{\rm\pmb{q}}\hat{\rm\pmb{r}} &= (iq_x + jq_y + kq_z + q_w )(ir_x + jr_y + kr_z + r_w)\\ &= i(q_yr_z−q_zr_y + r_wq_x + q_wr_x)\\ &\quad + j(q_zr_x − q_xr_z + r_wq_y + q_wr_y )\\ &\quad + k(q_xr_y − q_yr_x + r_wq_z + q_wr_z )\\ &\quad + q_wr_w − q_xr_x − q_yr_y − q_zr_z\\ &= ({\rm\pmb{q}}_v \times {\rm\pmb{r}}_v + r_w{\rm\pmb{q}}_v + q_w{\rm\pmb{r}}_v, q_wr_w − {\rm\pmb{q}}_v \cdot {\rm\pmb{r}}_v ) \end{aligned} \end{equation} \tag{4.32} \]

从这个方程可以看出，我们同时使用叉积和点积来计算两个四元数的乘法。

除了四元数的定义，还需要加法、共轭、范数和单位元素的定义：

\[\begin{equation} \begin{aligned} \pmb{加法:} \hat{\rm\pmb{q}} + \hat{\rm\pmb{r}} = ({\rm\pmb{q}}_v, q_w) + ({\rm\pmb{r}}_v ,r_w ) = ({\rm\pmb{q}}_v + {\rm\pmb{r}}_v, q_w + r_w ) \\ \pmb{共轭:} \hat{\rm\pmb{q}}^∗ = ({\rm\pmb{q}}_v, q_w)^∗ = (−{\rm\pmb{q}}_v, q_w) \\ \pmb{范数:} n(\hat{\rm\pmb{q}}) = \sqrt{\hat{\rm\pmb{q}}\hat{\rm\pmb{q}}^*} = \sqrt{\hat{\rm\pmb{q}}^*\hat{\rm\pmb{q}}} = \sqrt{{\rm\pmb{q}}_v \cdot {\rm\pmb{q}}_v + q_w^2}\\ = \sqrt{q^2_x + q^2_y + q^2_z + q^2_w} \\ \pmb{单位元素:} \hat{\rm\pmb{i}} = (\rm\pmb{0},1) \\ \end{aligned} \end{equation} \tag{4.33} \]

当\(n(\hat{\rm\pmb{q}}) = \sqrt{\hat{\rm\pmb{q}}\hat{\rm\pmb{q}}^*}\)被简化时（结果如上所示），虚部抵消，只剩下实部。范数有时表示为\(||\hat{\rm\pmb{q}}| = n(\hat{\rm\pmb{q}})\)[1105]。上面的结果是可以导出一个乘法逆，用\({\hat{\rm\pmb{q}}}^{-1}\)表示。方程\({\hat{\rm\pmb{q}}}^{-1}{\hat{\rm\pmb{q}}}={\hat{\rm\pmb{q}}}^{-1}{\hat{\rm\pmb{q}}}=1\)对逆必须成立（这对于乘法逆是常见的）。我们从范数的定义推导出一个公式：

\[{n(\hat{\rm\pmb{q}})}^2 = \hat{\rm\pmb{q}}\hat{\rm\pmb{q}}^* \Longleftrightarrow \frac{\hat{\rm\pmb{q}}\hat{\rm\pmb{q}}^*}{{n(\hat{\rm\pmb{q}})}^2} = 1 \tag{4.34} \]

这给出了乘法逆，如下所示：

\[\pmb{逆:} {\hat{\rm\pmb{q}}}^{-1} = \frac{1}{{n(\hat{\rm\pmb{q}})}^2}\hat{\rm\pmb{q}}^* \]

逆的公式使用了标量乘法，可以从方程4.3.1中看到的乘法推导出这个运算：\(s{\hat{\rm\pmb{q}}} = (\pmb{0}, s)({\rm\pmb{q}}_v, q_w) = (s{\rm\pmb{q}}_v, sq_w)\), 并且\({\hat{\rm\pmb{q}}}s = ({\rm\pmb{q}}_v, q_w)(\pmb{0}, s) = (s{\rm\pmb{q}}_v, sq_w)\)。这意味着标量乘法是可交换的：\(s{\hat{\rm\pmb{q}}}={\hat{\rm\pmb{q}}}s = (s{\rm\pmb{q}}_v, sq_w)\)。

以下规则集合很容易从定义中推导出来：

\[\begin{equation} \begin{aligned} \pmb{共轭规则:} (\hat{\rm\pmb{q}}^*)^* = \hat{\rm\pmb{q}} \\ (\hat{\rm\pmb{q}} + \hat{\rm\pmb{r}})^* = \hat{\rm\pmb{q}}^* + \hat{\rm\pmb{r}}^* \\ (\hat{\rm\pmb{q}}\hat{\rm\pmb{r}})^* = \hat{\rm\pmb{r}}^*\hat{\rm\pmb{q}}^* \end{aligned} \end{equation} \tag{4.36} \]

\[\begin{equation} \begin{aligned} \pmb{范数规则:} n(\hat{\rm\pmb{q}}^*) = n(\hat{\rm\pmb{q}}) \\ n(\hat{\rm\pmb{q}}\hat{\rm\pmb{r}}) = n(\hat{\rm\pmb{q}})n(\hat{\rm\pmb{r}}) \\ \end{aligned} \end{equation} \tag{4.37} \]

\[\begin{equation} \begin{aligned} \pmb{乘法法则:}\\ \pmb{线性度：} \hat{\rm\pmb{p}}(s\hat{\rm\pmb{q}} + t\hat{\rm\pmb{r}}) &= s\hat{\rm\pmb{p}}\hat{\rm\pmb{q}} + t\hat{\rm\pmb{p}}\hat{\rm\pmb{r}} \\ (s\hat{\rm\pmb{p}} + t\hat{\rm\pmb{q}})\hat{\rm\pmb{r}} &= s\hat{\rm\pmb{p}}\hat{\rm\pmb{r}} + t\hat{\rm\pmb{q}}\hat{\rm\pmb{r}} \\ \pmb{结合性：} \hat{\rm\pmb{p}}(\hat{\rm\pmb{q}}\hat{\rm\pmb{r}}) &= (\hat{\rm\pmb{p}}\hat{\rm\pmb{q}})\hat{\rm\pmb{r}} \end{aligned} \end{equation} \tag{4.38} \]

单位四元数\(\hat{\rm\pmb{q}} = ({\rm\pmb{q}}_v, q_w)\)使得\(n(\hat{\rm\pmb{q}}) = 1\)。由此可知\(\hat{\rm\pmb{q}}\)可写作:

\[\hat{\rm\pmb{q}} = ({\rm{sin}}{\phi}{\rm\pmb{u}}_q, {\rm{cos}}{\phi}) = {\rm{sin}}{\phi}{\rm\pmb{u}}_q + {\rm{cos}}{\phi} \tag{4.39} \]

对于某个三维向量\({\rm\pmb{u}}_q\)，使得\(||{\rm\pmb{u}}_q|| = 1\)，因为

\[n(\hat{\rm\pmb{q}}) = n({\rm{sin}}{\phi}{\rm\pmb{u}}_q, {\rm{cos}}{\phi}) = \sqrt{{\rm{sin}}^2{\phi}({\rm\pmb{u}}_q \cdot {\rm\pmb{u}}_q) + {\rm{cos}}^2{\phi}} \\ = \sqrt{{\rm{sin}}^2{\phi} + {\rm{cos}}^2{\phi}} \tag{4.40} = 1 \]

当且仅当\({\rm\pmb{u}}_q \cdot {\rm\pmb{u}}_q = 1 = ||{\rm\pmb{u}}_q||^2\). 正如将在下一节中看到的，单位四元数非常适合以最有效的方式创建旋转和方向。但在此之前，会为单位四元数引入一些额外的操作。

对于复数，二维单位向量可以写成\({\rm{cos}}{\phi} + i{\rm{sin}}{\phi} = e^{i\phi}\)。 四元数的等价物是

\[\hat{\rm\pmb{q}} = {\rm{sin}}{\phi}{\rm\pmb{u}}_q + {\rm{cos}}{\phi} = e^{\phi{\rm\pmb{u}}_q} \tag{4.41} \]

根据公式4.41，单位四元数的对数函数和幂函数为：

\[\begin{equation} \begin{aligned} \pmb{对数:} {\rm{log}}(\hat{\rm\pmb{q}}) = {\rm{log}}(e^{\phi{\rm\pmb{u}}_q}) = \phi{\rm\pmb{u}}_q \\ \pmb{指数:} \hat{\rm\pmb{q}}^t = ({\rm{sin}}{\phi}{\rm\pmb{u}}_q + {\rm{cos}}{\phi}) ^ t = e^{{\phi}t{\rm\pmb{u}}_q} = {\rm{sin}}({\phi}t){\rm\pmb{u}}_q + {\rm{cos}}({\phi}t)\\ \end{aligned} \end{equation} \tag{4.42} \]

4.3.2 四元数变换

我们现在将研究四元数集的一个子类，即单位长度的子类，称为单位四元数。关于单位四元数的最重要的事实是它们可以表示任何三维旋转，而且这种表示非常紧凑和简单。

现在我们将描述是什么使单位四元数对旋转和方向如此有用。首先，将点或向量的四个坐标\(\rm\pmb{p} = (p_x\ p_y\ p_z\ p_w)^T\)代入四元数\(\hat{\rm\pmb{p}}\)的分量，并假设我们有一个单位四元数\(\hat{\rm\pmb{q}} = ({\rm{sin}}{\phi}{\rm\pmb{u}}_q, {\rm{cos}}{\phi})\)。可以证明

\[\hat{\rm\pmb{q}}\hat{\rm\pmb{p}}\hat{\rm\pmb{q}}^{-1} \tag{4.43} \]

将\(\hat{\rm\pmb{p}}\)（以及点\(\rm\pmb{p}\)）绕轴\({\rm\pmb{u}}_q\)旋转角度\(2\phi\)。注意，由于\(\hat{\rm\pmb{q}}\)是一个单位四元数，\(\hat{\rm\pmb{q}}^{-1} = \hat{\rm\pmb{q}}^{*}\)。参见图4.9。

图4.9. 由单位四元数表示的旋转变换的图示，\(\hat{\rm\pmb{q}} = ({\rm{sin}}{\phi}{\rm\pmb{u}}_q, {\rm{cos}}{\phi})\)。变换围绕轴\({\rm\pmb{u}}_q\)旋转 \(2\phi\)弧度。

\(\hat{\rm\pmb{q}}\)的任何非零实数倍数也表示相同的变换，这意味着\(\hat{\rm\pmb{q}}\)和\(-\hat{\rm\pmb{q}}\)表示相同的旋转。也就是说，取反轴\(\rm\pmb{u}_q\)和实部\(q_w\)会创建一个与原始四元数完全一样旋转的四元数。这也意味着从矩阵中提取四元数可以返回\(\hat{\rm\pmb{q}}\)或\(-\hat{\rm\pmb{q}}\)。

给定两个单位四元数\(\hat{\rm\pmb{q}}\)和\(\hat{\rm\pmb{r}}\)，首先应用\(\hat{\rm\pmb{q}}\)再应用\(\hat{\rm\pmb{r}}\)到四元数\(\hat{\rm\pmb{p}}\)（可以解释为点\({\rm\pmb{p}}\)）的级联由公式 4.44 给出：

\[\hat{\rm\pmb{r}}(\hat{\rm\pmb{q}}\hat{\rm\pmb{p}}\hat{\rm\pmb{q}}^{∗})\hat{\rm\pmb{r}}^{∗} = (\hat{\rm\pmb{r}}\hat{\rm\pmb{q}})\hat{\rm\pmb{p}}(\hat{\rm\pmb{r}}\hat{\rm\pmb{q}})^{∗} = \hat{\rm\pmb{c}}\hat{\rm\pmb{p}}\hat{\rm\pmb{c}}^{∗} \tag{4.44} \]

这里，\(\hat{\rm\pmb{c}} = \hat{\rm\pmb{r}}\hat{\rm\pmb{q}}\)是单位四元数，表示单位四元数\(\hat{\rm\pmb{q}}\)和\(\hat{\rm\pmb{r}}\)的级联。

矩阵转换

由于经常需要将几种不同的变换组合起来，而且大部分都是矩阵形式，因此需要一种方法将式4.43转化为矩阵。四元数 \(\hat{\rm\pmb{q}}\)可以转换为矩阵\(\rm\pmb{M}^q\)，如公式4.45[1633,1634]所示：

\[{\rm\pmb{M}}^q = \left( \begin{matrix} 1 − s(q^2_y + q^2_z) & s(q_xq_y − q_wq_z) & s(q_xq_z + q_wq_y) & 0\\ s(q_xq_y + q_wq_z ) & 1 − s(q^2_x + q^2_z ) & s(q_yq_z − q_wq_x ) & 0\\ s(q_xq_z − q_wq_y) & s(q_yq_z + q_wq_x ) & 1 − s(q^2_x + q^2_y) & 0\\ 0 & 0 & 0 & 1 \end{matrix} \right) \tag{4.45} \]

在这里，标量是\(s = 2/(n(\hat{\rm\pmb{q}}))^2\)。 对于单位四元数，上式可简化为：

\[{\rm\pmb{M}}^q = \left( \begin{matrix} 1 − 2(q^2_y + q^2_z) & 2(q_xq_y − q_wq_z) & 2(q_xq_z + q_wq_y) & 0\\ 2(q_xq_y + q_wq_z ) & 1 − 2(q^2_x + q^2_z ) & 2(q_yq_z − q_wq_x ) & 0\\ 2(q_xq_z − q_wq_y) & 2(q_yq_z + q_wq_x ) & 1 − 2(q^2_x + q^2_y) & 0\\ 0 & 0 & 0 & 1 \end{matrix} \right) \tag{4.46} \]

一旦构造了四元数，就不需要计算三角函数，因此转换过程在实践中很有效率。

从正交矩阵\(\rm\pmb{M}^q\)到单位四元数\(\hat{\rm\pmb{q}}\)的逆转换要复杂一些。此过程的关键是与公式4.46中的矩阵的以下差异：

\[\begin{equation} \begin{aligned} m^q_{21} − m^q_{12} = 4q_wq_x \\ m^q_{02} − m^q_{20} = 4q_wq_y \\ m^q_{10} − m^q_{01} = 4q_wq_z \\ \end{aligned} \end{equation} \tag{4.47} \]

这些方程的含义是，如果\(q_w\)已知，则可以计算向量\({\rm\pmb{v}}_q\)的值，从而推导出\(\hat{\rm\pmb{q}}\)。\(\rm\pmb{M}^q\)的迹由下式计算:

\[\begin{equation} \begin{aligned} tr(\rm\pmb{M}^q) &= 4 - 2s(q^2_x + q^2_y + q^2_z) = 4(1-\frac{q^2_x + q^2_y + q^2_z}{q^2_x + q^2_y + q^2_z+q^2_w}) \\ &= \frac{4q^2_w}{q^2_x + q^2_y + q^2_z+q^2_w} = \frac{4q^2_w}{(n(\hat{\rm\pmb{q}}))^2} \end{aligned} \end{equation} \tag{4.48} \]

对于单位四元数来说，可得出以下转换：

\[\begin{equation} \begin{aligned} q_w &= \frac{1}{2}\sqrt{tr({\rm\pmb{M}}^q)} \\ q_x &= \frac{(m^q_{21} − m^q_{12})}{4q_w} \\ q_y &= \frac{(m^q_{02} − m^q_{20})}{4q_w} \\ q_z &= \frac{(m^q_{10} − m^q_{01})}{4q_w} \\ \end{aligned} \end{equation} \tag{4.49} \]

为了有一个数值稳定的程序[1634]，小数除法应该被避免。因此，首先设置\(t = q^2_w − q^2_x − q^2_y − q^2_z\)，由此得出:

\[\begin{equation} \begin{aligned} m_{00} &= t + 2q^2_x \\ m_{11} &= t + 2q^2_y \\ m_{22} &= t + 2q^2_z \\ u &= m_{00} + m_{11} + m_{22} = t + 2q^2_w\\ \end{aligned} \end{equation} \tag{4.50} \]

这反过来意味着\(m_{00}\)、\(m_{11}\)、\(m_{22}\)和\(u\)中的最大值，确定了\(q_x\)、\(q_y\)、\(q_z\)和\(q_w\)中最大的。如果\(q_w\)最大，则使用公式4.49导出四元数。否则，我们注意到以下内容成立：

\[\begin{equation} \begin{aligned} 4q^2_x &= +m_{00} − m_{11} − m_{22} + m_{33} \\ 4q^2_y &= −m_{00} + m_{11} − m_{22} + m_{33} \\ 4q^2_z &= −m_{00} − m_{11} + m_{22} + m_{33} \\ 4q^2_w &= tr({\rm\pmb{M}}^q)\\ \end{aligned} \end{equation} \tag{4.51} \]

然后使用上述适当方程计算\(q_x\)、\(q_y\)和\(q_z\)中的最大值，然后使用方程4.47计算\(\hat{\rm\pmb{q}}\)的其余分量。Schüler[1588]提出了一种无分支但使用四个平方根的变体。

球面线性插值

球面线性插值是一种运算，在给定两个单位四元数\(\hat{\rm\pmb{q}}\)和\(\hat{\rm\pmb{r}}\)以及参数\(t∈[0,1]\)的情况下，计算插值四元数。例如，这对于动画对象很有用。它对于插值相机方向没有用，因为相机的“向上”矢量在插值过程中可能会倾斜，通常是一种干扰效果。

此运算的代数形式由复合四元数\(\hat{\rm\pmb{s}}\)表示，如下所示：

\[\hat{\rm\pmb{s}}(\hat{\rm\pmb{q}},\hat{\rm\pmb{r}},t) = (\hat{\rm\pmb{r}}\hat{\rm\pmb{q}}^{−1})^t\hat{\rm\pmb{q}} \tag{4.52} \]

但是，对于软件实现，以下形式更合适，其中slerp代表球面线性插值：

\[\hat{\rm\pmb{s}}(\hat{\rm\pmb{q}},\hat{\rm\pmb{r}},t) = slerp(\hat{\rm\pmb{q}},\hat{\rm\pmb{r}},t) = \frac{sin(\phi(1-t))}{sin\phi}\hat{\rm\pmb{q}} + \frac{sin(\phi t)}{sin\phi}\hat{\rm\pmb{r}} \tag{4.53} \]

为了计算这个方程所需的φ，可以使用以下式子：\(cos\phi = q_xr_x + q_yr_y + q_zr_z + q_wr_w\)[325]。对于 \(t∈[0,1]\)，slerp函数计算（唯一的\(^2\)）插值四元数，它们共同构成四维单位球面上从\(\hat{\rm\pmb{q}}(t = 0)\)到\(\hat{\rm\pmb{r}}(t = 1)\)的最短弧。圆弧位于由\(\hat{\rm\pmb{q}}\),\(\hat{\rm\pmb{r}}\)和原点给出的平面与四维单位球面的交点形成的圆上。如图 4.10 所示。计算出的旋转四元数以恒定速度绕固定轴旋转。像这样的曲线具有恒定的速度，因此加速度为零，称为测地曲线[229]。过原点的平面与球面的交点在球面上形成一个大圆，这个圆的一部分称为大圆弧。

图4.10. 单位四元数表示为单位球面上的点。函数slerp用于四元数之间的插值，插值的路径是球体上的一个大圆弧。请注意，从\(\hat{\rm\pmb{q}}_1\)到\(\hat{\rm\pmb{q}}_2\)和从\(\hat{\rm\pmb{q}}_1\)到\(\hat{\rm\pmb{q}}_3\)到\(\hat{\rm\pmb{q}}_2\)的插值不是一回事，即使它们到达相同的方向。

slerp函数非常适合在两个方位之间进行插值，并且有良好的效果（固定轴，恒速）。使用多个欧拉角进行插值时，情况并非如此。实际上，直接计算slerp是一项涉及调用三角函数的昂贵操作。Malysau[1114]讨论了将四元数集成到渲染管线中。他指出，当不使用slerp而只是在像素着色器中对四元数进行归一化时，对于90度角，三角形方向的误差最大为 4度。光栅化三角形时，此误差率是可以接受的。Li[1039, 1040]提供了更快的增量方法来计算slerps，并且不会牺牲任何准确性。Eberly[406]提出了一种仅使用加法和乘法计算slerps的快速技术。

当有两个以上的方向时，比如\(\hat{\rm\pmb{q}}_0,\hat{\rm\pmb{q}}_1,...,\hat{\rm\pmb{q}}_n\)，并且我们想要从\(\hat{\rm\pmb{q}}_0\)到\(\hat{\rm\pmb{q}}_1\)到\(\hat{\rm\pmb{q}}_2\)，依此类推直到\(\hat{\rm\pmb{q}}_{n-1}\)，可以直接使用slerp。现在，当我们接近\(\hat{\rm\pmb{q}}_i\)时，我们将使用\(\hat{\rm\pmb{q}}_{i-1}\)和\(\hat{\rm\pmb{q}}_i\)作为 slerp 的参数。通过\(\hat{\rm\pmb{q}}_i\)后，我们将使用\(\hat{\rm\pmb{q}}_i\)和\(\hat{\rm\pmb{q}}_{i+1}\)作为 slerp 的参数。这将导致方向插值中出现突然的抖动，如图4.10所示。这类似于当点被线性插值时发生的情况；参见第720页图17.3的右上部分。有些读者可能希望在阅读第17章中的样条曲线后重新阅读以下段落。

更好的插值方法是使用某种样条。我们在\(\hat{\rm\pmb{q}}_i\)和\(\hat{\rm\pmb{q}}_{i+1}\)之间引入了四元数\(\hat{\rm\pmb{a}}_i\)和\(\hat{\rm\pmb{a}}_{i+1}\)。球面三次插值可以在四元数\(\hat{\rm\pmb{q}}_i\)，\(\hat{\rm\pmb{a}}_i\)，\(\hat{\rm\pmb{a}}_{i+1}\)和\(\hat{\rm\pmb{q}}_{i+1}\)的集合内定义。令人惊讶的是，这些额外的四元数计算如下[404]\(^3\)：

\[\hat{\rm\pmb{a}}_i = \hat{\rm\pmb{q}}_i \, {\rm{exp}} \, [-\frac{{\rm{log}}(\hat{\rm\pmb{q}}^{-1}_{i}\hat{\rm\pmb{q}}_{i-1})+{\rm{log}}(\hat{\rm\pmb{q}}^{-1}_{i}\hat{\rm\pmb{q}}_{i+1})}{4}]\tag{4.54} \]

通过三次样条平滑算法，\(\hat{\rm\pmb{q}}_i\)和\(\hat{\rm\pmb{a}}_i\)用来对四元数进行球面插值，如公式4.55所示：

\[{\rm{squad}}(\hat{\rm\pmb{q}}_i,\hat{\rm\pmb{q}}_{i+1},\hat{\rm\pmb{a}}_i,\hat{\rm\pmb{a}}_{i+1},t) = \\ {\rm{slerp}}({\rm{slerp}}(\hat{\rm\pmb{q}}_i,\hat{\rm\pmb{q}}_{i+1},t),{\rm{slerp}}(\hat{\rm\pmb{a}}_i,\hat{\rm\pmb{a}}_{i+1},t),2t(1 − t)) \tag{4.55} \]

从上面可以看出，squad函数是使用slerp从同样的球面插值构建的（第17.1.1节有关点的重复线性插值的信息）。插值将通过初始方向\(\hat{\rm\pmb{q}}_i,i ∈ [0,...,n − 1]\)，但不通过\(\hat{\rm\pmb{a}}_i\)——这些用于指示初始方向的切线方向。

从一个向量到另一个向量的旋转

一个常见的操作是通过尽可能最短的路径从一个方向\(\rm\pmb{s}\)转换到另一个方向\(\rm\pmb{t}\)。四元数的数学运算极大地简化了这个过程，并显示了四元数与这种表示的密切关系。首先，将\(\rm\pmb{s}\)和\(\rm\pmb{t}\)归一化。然后计算单位旋转轴，称为\(\rm\pmb{u}\)，计算为\(\rm\pmb{u} = (\pmb{s} × \pmb{t})/||\pmb{s} × \pmb{t}||\)。 接下来，\(e = {\rm\pmb{s}} · {\rm\pmb{t}} = cos(2\phi)\)和\(||{\rm\pmb{s}} × {\rm\pmb{t}}|| = sin(2\phi)\)，其中\(2\phi\)是\({\rm\pmb{s}}\)和\({\rm\pmb{t}}\)之间的角度。表示从\(\rm\pmb{s}\)到\(\rm\pmb{t}\)的旋转的四元数是\(\hat{\rm\pmb{q}} = (sin{\phi}{\rm\pmb{u}}, cos{\phi})\)。事实上，使用半角关系和三角恒等式简化\(\hat{\rm\pmb{q}} = (\frac{sin{\phi}}{sin2{\phi}}({\rm\pmb{s}}×{\rm\pmb{t}}), cos{\phi})\)，可得出[1197]：

\[\hat{\rm\pmb{q}} = ({\rm\pmb{q}}_v,q_w) = (\frac{1}{\sqrt{2(1+e)}}({\rm\pmb{s}}×{\rm\pmb{t}}), \frac{\sqrt{2(1+e)}}{2}) \tag{4.56} \]

以这种方式直接生成四元数（相对于标准化叉积\({\rm\pmb{s}}×{\rm\pmb{t}}\)）避免了当\({\rm\pmb{s}}\)和\({\rm\pmb{t}}\)指向几乎相同的方向时的数值不稳定[1197]。当\({\rm\pmb{s}}\)和\({\rm\pmb{t}}\)指向相反的方向时，两种方法都会出现稳定性问题，因为发生了除以零的情况。当检测到这种特殊情况时，任何垂直于\({\rm\pmb{s}}\)的旋转轴都可以用来旋转到\({\rm\pmb{t}}\)。

有时我们需要从\({\rm\pmb{s}}\)到\({\rm\pmb{t}}\)旋转的矩阵表示。在对等式4.46进行一些代数和三角简化后，旋转矩阵变为[1233]:

\[{\rm\pmb{R}}({\rm\pmb{s}}, {\rm\pmb{t}}) = \left( \begin{matrix} e + hv_x^2 & hv_xv_y − v_z & hv_xv_z + v_y & 0 \\ hv_xv_y + v_z & e + hv_y^2 & hv_yv_z − v_x & 0 \\ hv_xv_z − v_y & hv_yv_z + v_x & e + hv_z^2 & 0 \\ 0 & 0 & 0 & 1 \\ \end{matrix} \right) \tag{4.57} \]

在这个等式中，我们使用了以下中间计算：

\[\begin{equation} \begin{aligned} {\rm\pmb{v}} &= {\rm\pmb{s}} × {\rm\pmb{t}} \\ e &= cos(2{\phi}) = {\rm\pmb{s}} \cdot {\rm\pmb{t}} \\ h &= \frac{1-cos(2\phi)}{sin^2(2\phi)} = \frac{1-e}{{\rm\pmb{v}} \cdot {\rm\pmb{v}}} = \frac{1}{1+e}\\ \end{aligned} \end{equation} \tag{4.58} \]

可以看出，由于简化，所有平方根和三角函数都消失了，因此这是创建矩阵的有效方法。请注意，公式4.57的结构类似于公式4.30的结构，并注意后一种形式如何不需要三角函数。

请注意，当\({\rm\pmb{s}}\)和\({\rm\pmb{t}}\)平行或接近平行时必须小心，因为这时\(||{\rm\pmb{s}} × {\rm\pmb{t}}|| ≈ 0\)。如果\(\phi ≈ 0\)，那么我们可以返回单位矩阵。然而，如果\(2\phi ≈ \pi\)，那么我们可以绕任意轴旋转\(\pi\)弧度。该轴可以是\({\rm\pmb{s}}\)和任何其他不平行于\({\rm\pmb{s}}\)的向量的叉积（第4.2.4节）。Möller和Hughes使用Householder矩阵以不同的方式处理这种特殊情况[1233]。

---

\(^1\)公平地说，Robinson[1502]在1958年使用四元数进行刚体模拟。
\(^2\)当且仅当\(\hat{\rm\pmb{q}}\)和\(\hat{\rm\pmb{r}}\)不相反。
\(^3\)Shoemake[1633]给出了另一个推导。