四元数

  四元数可以用于表示旋转和方向，它在很多地方都比欧拉角和矩阵表示更加优秀。任何三维方向都可以表示为一个绕特定轴的简单旋转，给定一个旋转轴和旋转角度，可以直接将其转换为一个四元数，或者是从一个四元数中提取出旋转轴和旋转角度；但是对任意方向上的欧拉角进行转换是很困难的。四元数可以用于稳定且恒定速度的方向插值，这是欧拉角很难实现的。

  复数由一个实部和一个虚部组成，每个复数都可以使用两个实数进行表示，其中第二个实数要乘以$\sqrt{-1}$。类似地，四元数由四个部分组成，前三个值与旋转的轴有关，而旋转角度会对四个值都产生影响。由于四元数和向量都有四个分量，因此这里给四元数的符号加上一个小帽子来以示区别，即$\hat{\mathbf{q}}$。在本小节中，将从一些四元数的数学背景开始讲起，然后用它来构建各种有用的变换。

数学背景

  定义：四元数$\hat{\mathbf{q}}$可以使用下面的方法来进行定义，这些定义方法都是等价的。

$$\begin{array}{rl}\hat{\mathbf{q}}& =\left({\mathbf{q}}_v,q_w\right)=i{q}_x+j{q}_y+k{q}_z+q_w={\mathbf{q}}_v+q_w,\\ {\mathbf{q}}_v& =i{q}_x+j{q}_y+k{q}_z=\left(q_x,q_y,q_z\right),\\ i^2& =j^2=k^2=-1,jk=-kj=i,ki=-ik=j,ij=-ji=k.\end{array}\text{\hspace{1em}(4.31)}$$

  其中变量$q_w$是四元数$\hat{\mathbf{q}}$中的实数部分（实部），变量${\mathbf{q}}_{\mathbf{v}}$是四元数$\hat{\mathbf{q}}$中的虚数部分（虚部），$i,j,k$叫做虚数单位。四元数的结构和复数类似，但是复数只有一个虚部，而四元数则包含三个虚部。

  对于虚部${\mathbf{q}}_{\mathbf{v}}$，可以将其看作为一个三维向量，可以对其应用诸如加法、缩放、点乘、叉乘以及其他的向量操作。利用四元数的定义，可以得到两个四元数$\hat{\mathbf{q}},\hat{\mathbf{r}}$之间的乘法运算，如方程4.32所示。这里需要注意的是，虚数单位之间的乘法是不具备交换律的。

$$\begin{array}{rl}\mathbf{乘法}:\hat{\mathbf{q}}\hat{\mathbf{r}}=& \left(i{q}_x+j{q}_y+k{q}_z+q_w\right)\left(i{r}_x+j{r}_y+k{r}_z+r_w\right)\\ =& i\left(q_y{r}_z-q_z{r}_y+r_w{q}_x+q_w{r}_x\right)\\ & +j\left(q_z{r}_x-q_x{r}_z+r_w{q}_y+q_w{r}_y\right)\\ & +k\left(q_x{r}_y-q_y{r}_x+r_w{q}_z+q_w{r}_z\right)\\ & +q_w{r}_w-q_x{r}_x-q_y{r}_y-q_z{r}_z\\ =& \left({\mathbf{q}}_v\times {\mathbf{r}}_v+r_w{\mathbf{q}}_v+q_w{\mathbf{r}}_v,q_w{r}_w-{\mathbf{q}}_v\cdot {\mathbf{r}}_v\right).\end{array}\text{\hspace{1em}(4.32)}$$

  从上述计算公式中可以看出，在计算两个四元数的乘法时，同时使用了点乘和叉乘。

  除了四元数的定义和乘法计算公式之外，下面是四元数的加法、共轭、范数以及其他定义：

$$\begin{array}{rlr}& \mathbf{加法}:& \hat{\mathbf{q}}+\hat{\mathbf{r}}=({\mathbf{q}}_v,q_w)+({\mathbf{r}}_v,r_w)=({\mathbf{q}}_v+{\mathbf{r}}_v,q_w+r_w)\\ & \mathbf{共轭}:& {\hat{\mathbf{q}}}^{\ast }={\left({\mathbf{q}}_v,q_w\right)}^{\ast }=\left(-{\mathbf{q}}_v,q_w\right)\\ & \mathbf{模长}:& \begin{array}{rl}n(\hat{\mathbf{q}})& =\sqrt{\hat{\mathbf{q}}{\hat{\mathbf{q}}}^{\ast }}=\sqrt{{\hat{\mathbf{q}}}^{\ast }\hat{\mathbf{q}}}=\sqrt{{\mathbf{q}}_v\cdot {\mathbf{q}}_v+q_w^2}\\ & =\sqrt{q_x^2+q_y^2+q_z^2+q_w^2}.\end{array}\\ & \mathbf{虚数单位}:& \hat{\mathbf{i}}=(0,1)\end{array}\text{\hspace{1em}(4.33)}$$

  当$n(\hat{\mathbf{q}})=\sqrt{\hat{\mathbf{q}}{\hat{\mathbf{q}}}^{\ast }}$化简之后，可以看到最终结果的虚部消失了，只剩下实部（即一个实数），这个实数叫做虚数$\hat{\mathbf{q}}$的模长，有时候也会使用$\Vert \hat{\mathbf{q}}\Vert =n(\hat{\mathbf{q}})$来表示一个虚数的模长。使用符号来${\hat{\mathbf{q}}}^{-1}$表示一个四元数的逆，四元数的逆有这样一个性质，即${\hat{\mathbf{q}}}^{-1}\hat{\mathbf{q}}=\hat{\mathbf{q}}{\hat{\mathbf{q}}}^{-1}=1$。根据虚数模长的定义，可以推导出:

$$n(\hat{\mathbf{q}}{)}^2=\hat{\mathbf{q}}{\hat{\mathbf{q}}}^{\ast }⟺\frac{\hat{\mathbf{q}}{\hat{\mathbf{q}}}^{\ast }}{n(\hat{\mathbf{q}}{)}^2}=1\text{\hspace{1em}(4.34)}$$

  从方程4.34可以推导出四元数逆的计算公式：

$$\mathbf{四元数的逆}:{\hat{\mathbf{q}}}^{-1}=\frac{1}{n(\hat{\mathbf{q}}{)}^2}{\hat{\mathbf{q}}}^{\ast }\text{\hspace{1em}(4.35)}$$

  方程4.35使用了四元数的标量乘法，其计算方法可以从方程4.32中推导出来，即

$$\begin{gathered} s\hat{\mathbf{q}}=(0,s)\left({\mathbf{q}}_v,q_w\right)=\left(s{\mathbf{q}}_v,s{q}_w\right) \\ \hat{\mathbf{q}}s=\left({\mathbf{q}}_v,q_w\right)(0,s)=\left(s{\mathbf{q}}_v,s{q}_w\right) \end{gathered}$$

  这也意味着复数的标量乘法具备交换律，即：$s\hat{\mathbf{q}}=\hat{\mathbf{q}}s=\left(s{\mathbf{q}}_v,s{q}_w\right)$。

  利用上面提到的基本定义，可以推导出以下的四元数运算规则：

$$\begin{array}{rl}{\left({\hat{\mathbf{q}}}^{\ast }\right)}^{\ast }& =\hat{\mathbf{q}},\\ \mathbf{共轭法则}:(\hat{\mathbf{q}}+\hat{\mathbf{r}}{)}^{\ast }& ={\hat{\mathbf{q}}}^{\ast }+{\hat{\mathbf{r}}}^{\ast },\\ (\hat{\mathbf{q}}\hat{\mathbf{r}}{)}^{\ast }& ={\hat{\mathbf{r}}}^{\ast }{\hat{\mathbf{q}}}^{\ast }.\end{array}\text{\hspace{1em}(4.36)}$$

$$\begin{array}{rl}\mathbf{模长法则}:n\left({\hat{\mathbf{q}}}^{\ast }\right)& =n(\hat{\mathbf{q}})\\ n(\hat{\mathbf{q}}\hat{\mathbf{r}})& =n(\hat{\mathbf{q}})n(\hat{\mathbf{r}})\end{array}\text{\hspace{1em}(4.37)}$$

$$\begin{array}{rl}乘法分配律：& \begin{array}{l}\hat{\mathbf{p}}(s\hat{\mathbf{q}}+t\hat{\mathbf{r}})=s\hat{\mathbf{p}}\hat{\mathbf{q}}+t\hat{\mathbf{p}}\hat{\mathbf{r}}\\ (s\hat{\mathbf{p}}+t\hat{\mathbf{q}})\hat{\mathbf{r}}=s\hat{\mathbf{p}}\hat{\mathbf{r}}+t\hat{\mathbf{q}}\hat{\mathbf{r}}\end{array}\\ 乘法结合律：& \hat{\mathbf{p}}(\hat{\mathbf{q}}\hat{\mathbf{r}})=(\hat{\mathbf{p}}\hat{\mathbf{q}})\hat{\mathbf{r}}\end{array}\text{\hspace{1em}(4.38)}$$

  单位四元数$\hat{\mathbf{q}}=\left({\mathbf{q}}_v,q_w\right)$的长度（模长）为$n\left(\hat{\mathbf{q}}\right)=1$。因此对于某些长度为1的三维向量${\mathbf{u}}_q$，可以将这个四元数改写成如下形式：

$$\hat{\mathbf{q}}=\left(\sin \varphi {\mathbf{u}}_q,\cos \varphi \right)=\sin \varphi {\mathbf{u}}_q+\cos \varphi =\cos \varphi +\sin \varphi ({\mathbf{u}}_{qx}+{\mathbf{u}}_{qy}+{\mathbf{u}}_{qz})\text{\hspace{1em}(4.39)}$$

  当且仅当${\mathbf{u}}_q\cdot {\mathbf{u}}_q=1={\Vert {\mathbf{u}}_q\Vert }^2$的时候，才能够这样改写，这是因为：

$$\begin{array}{rl}n(\hat{\mathbf{q}})& =n\left(\sin \varphi {\mathbf{u}}_q,\cos \varphi \right)=\sqrt{{\sin }^2\varphi \left({\mathbf{u}}_q\cdot {\mathbf{u}}_q\right)+{\cos }^2\varphi }\\ & =\sqrt{{\sin }^2\varphi +{\cos }^2\varphi }=1\end{array}\text{\hspace{1em}(4.40)}$$

  对于复数而言，一个二维单位向量可以被写成$\cos \varphi +i\sin \varphi =e^{i\varphi }$的形式（复变函数中的欧拉公式）；而对于一个四元数来说，其等价形式如下：

$$\hat{\mathbf{q}}=\sin \varphi {\mathbf{u}}_q+\cos \varphi =e^{\varphi {\mathbf{u}}_q}\text{\hspace{1em}(4.41)}$$

  对于一个单位四元数而言，其对数运算和幂运算遵循以下规则，这些规则可以从方程4.41中推导出来：

$$\begin{array}{rl}& \mathbf{Logarithm}:\log (\hat{\mathbf{q}})=\log \left(e^{\varphi {\mathbf{u}}_q}\right)=\varphi {\mathbf{u}}_q\\ & \mathbf{Power}:{\hat{\mathbf{q}}}^t={\left(\sin \varphi {\mathbf{u}}_q+\cos \varphi \right)}^t=e^{\varphi t{\mathbf{u}}_q}=\sin (\varphi t){\mathbf{u}}_q+\cos (\varphi t).\end{array}\text{\hspace{1em}(4.42)}$$

四元数变换

  现在将研究四元数中的一个子集，即单位四元数（unit quaternion），它们的模长为1。单位四元数可以用于表示任何的三维旋转，而且这种表示方式非常紧凑和简单。

  这里将介绍为什么单元四元数对于旋转和方向如此有用。假设现在有一个点（或者向量）$\mathbf{p}={\left(\begin{array}{llll}{p}_x& p_y& p_z& p_w\end{array}\right)}^T$，首先将这个点坐标的四个分量对应放入一个四元数$\hat{\mathbf{p}}$的各个分量中；假设现在还有一个单位四元数$\hat{\mathbf{q}}=(\sin \varphi {\mathbf{u}}_q+\cos \varphi )$，我们对这两个四元数进行如下操作：

$$\hat{\mathbf{q}}\hat{\mathbf{p}}{\hat{\mathbf{q}}}^{-1}\text{\hspace{1em}(4.43)}$$

  这个操作意味着将四元数$\hat{\mathbf{p}}$（即点$\mathbf{p}$）以${\mathbf{u}}_q$为旋转轴，旋转了$2\varphi$角度。请注意，由于$\hat{\mathbf{q}}$是一个单位四元数，因此根据方程4.35，有${\hat{\mathbf{q}}}^{-1}={\hat{\mathbf{q}}}^{\ast }$。这个绕任意轴旋转的过程如图4.9所示。

图4.9 通过使用一个单位四元数$\hat{\mathbf{q}}=(\sin \varphi {\mathbf{u}}_q+\cos \varphi )$完成的旋转操作，绕轴${\mathbf{u}}_q$旋转了$2\varphi$ 。

  任意非零倍数的$\hat{\mathbf{q}}$都能表示这个相同的旋转变换，也就是说，$\hat{\mathbf{q}}$和$-\hat{\mathbf{q}}$所表示的旋转变换也都是一样的（不考虑缩放），因此可以这样理解，将旋转轴${\mathbf{u}}_q$和实部$q_w$取反，可以创建一个与原始四元数完全相同的旋转变换。由于单位四元数$\hat{\mathbf{q}}$首先将目标绕轴旋转$\varphi$，然后$-\hat{\mathbf{q}}$再将目标绕轴旋转$\varphi$，并且两次旋转的方向是相同的，因此最终旋转了$2\varphi$。同时这也意味着从变换矩阵中提取出的四元数，可能是$\hat{\mathbf{q}}$或者$-\hat{\mathbf{q}}$中的任意一个。

  给定两个单位四元数$\hat{\mathbf{q}}$和$\hat{\mathbf{r}}$，首先对$\hat{\mathbf{p}}$（可以被认为一个点或者是一个向量）应用四元数$\hat{\mathbf{q}}$所代表的旋转变换，然后再应用四元数$\hat{\mathbf{r}}$所代表的变换，整个过程可以用方程4.44进行描述：

$$\hat{\mathbf{r}}\left(\hat{\mathbf{q}}\hat{\mathbf{p}}{\hat{\mathbf{q}}}^{\ast }\right){\hat{\mathbf{r}}}^{\ast }=(\hat{\mathbf{r}}\hat{\mathbf{q}})\hat{\mathbf{p}}(\hat{\mathbf{r}}\hat{\mathbf{q}}{)}^{\ast }=\hat{\mathbf{c}}\hat{\mathbf{p}}{\hat{\mathbf{c}}}^{\ast }\text{\hspace{1em}(4.44)}$$

  上述公式中的$\hat{\mathbf{c}}=\hat{\mathbf{r}}\hat{\mathbf{q}}$也是一个单位四元数，它代表两个单位四元数$\hat{\mathbf{q}}$和$\hat{\mathbf{r}}$的组合变换。

矩阵转换

  很多时候需要将好几种不同的变换组合在一起，这些变换大多数都以矩阵的形式给出，因此需要一种方法来将方程4.43中的四元数旋转，转换为矩阵形式。一个四元数$\hat{\mathbf{q}}$，可以转换为一个变换矩阵${\mathbf{M}}^{\mathbf{q}}$，其形式如下：

$${\mathbf{M}}^q=\left(\begin{array}{cccc}1-s\left(q_y^2+q_z^2\right)& s\left(q_x{q}_y-q_w{q}_z\right)& s\left(q_x{q}_z+q_w{q}_y\right)& 0\\ s\left(q_x{q}_y+q_w{q}_z\right)& 1-s\left(q_x^2+q_z^2\right)& s\left(q_y{q}_z-q_w{q}_x\right)& 0\\ s\left(q_x{q}_z-q_w{q}_y\right)& s\left(q_y{q}_z+q_w{q}_x\right)& 1-s\left(q_x^2+q_y^2\right)& 0\\ 0& 0& 0& 1\end{array}\right)\text{\hspace{1em}(4.45)}$$

  上述方程中的$s=2/(n(\hat{\mathbf{q}}){)}^2$，对于一个单位四元数而言，可以将其化简为：

$${\mathbf{M}}^q=\left(\begin{array}{cccc}1-2\left(q_y^2+q_z^2\right)& 2\left(q_x{q}_y-q_w{q}_z\right)& 2\left(q_x{q}_z+q_w{q}_y\right)& 0\\ 2\left(q_x{q}_y+q_w{q}_z\right)& 1-2\left(q_x^2+q_z^2\right)& 2\left(q_y{q}_z-q_w{q}_x\right)& 0\\ 2\left(q_x{q}_z-q_w{q}_y\right)& 2\left(q_y{q}_z+q_w{q}_x\right)& 1-2\left(q_x^2+q_y^2\right)& 0\\ 0& 0& 0& 1\end{array}\right)\text{\hspace{1em}(4.46)}$$

   当完成这个四元数的构建，在将其转换为矩阵形式的时候，上述转换方程中不包含任何三角函数运算，因此这个转换在实际应用中是很高效的。

  将一个正交矩阵${\mathbf{M}}^{\mathbf{q}}$反向转换为一个单位四元数$\hat{\mathbf{q}}$要更加复杂一点，这个过程的关键在于，可以通过消元法，将方程4.46中正交矩阵${\mathbf{M}}^{\mathbf{q}}$的某几项进行做差，即：

$$\begin{array}{l}{m}_{21}^q-m_{12}^q=4q_w{q}_x\\ m_{02}^q-m_{20}^q=4q_w{q}_y\\ m_{10}^q-m_{01}^q=4q_w{q}_z.\end{array}\text{\hspace{1em}(4.47)}$$

  方程4.47中的每一项都含有$q_w$，也就是说，如果我们知道了$q_w$的值是多少的话，我们就可以计算出向量${\mathbf{v}}_q$的各个分量，那么这个单位四元数$\hat{\mathbf{q}}$也就计算出来了。$q_w$可以通过计算矩阵${\mathbf{M}}^{\mathbf{q}}$的迹（trace：矩阵主对角线上的元素之和，等于矩阵的特征值之和）获得，计算方程如下：

$$\begin{array}{rl}\mathrm{tr}\left({\mathbf{M}}^q\right)& =4-2s\left(q_x^2+q_y^2+q_z^2\right)=4\left(1-\frac{q_x^2+q_y^2+q_z^2}{q_x^2+q_y^2+q_z^2+q_w^2}\right)\\ & =\frac{4q_w^2}{q_x^2+q_y^2+q_z^2+q_w^2}=\frac{4q_w^2}{(n(\hat{\mathbf{q}}){)}^2}.\end{array}\text{\hspace{1em}(4.48)}$$

  通过方程4.48求出矩阵${\mathbf{M}}^{\mathbf{q}}$的迹，然后就可以计算出单位四元数的实部$q_w$，最后计算出四元数虚部的其他几项即可，具体方程如下：

$$\begin{array}{ll}{q}_w=\frac{1}{2}\sqrt{\mathrm{tr}\left({\mathbf{M}}^q\right)},& q_x=\frac{m_{21}^q-m_{12}^q}{4q_w},\\ q_y=\frac{m_{02}^q-m_{20}^q}{4q_w},& q_z=\frac{m_{10}^q-m_{01}^q}{4q_w}.\end{array}\text{\hspace{1em}(4.49)}$$

  为了在计算过程中保持数值稳定，应当尽可能地避免除以一个过小的数值。为了实现这个目的，首先设$t=q_w^2-q_x^2-q_y^2-q_z^2$，然后又因为$q_x^2+q_y^2+q_z^2+q_w^2=1$（单位四元数的模长为1），可得：

$$\begin{array}{rl}{m}_{00}& =t+2q_x^2,\\ m_{11}& =t+2q_y^2,\\ m_{22}& =t+2q_z^2,\\ u& =m_{00}+m_{11}+m_{22}=t+2q_w^2,\end{array}\text{\hspace{1em}(4.50)}$$

  上述方程反过来意味着，$m_{00},m_{11},m_{22},u$中最大的那一项，对应决定了$q_w,q_x,q_y,q_z$中哪一项是最大的。如果计算出来实部$q_w$是最大的那一项，便可以使用方程4.49计算出这个四元数的其他分量，这样就不会因为除以一个过小的数值而带来精度上的误差。如果计算出来实部$q_w$并不是最大的那一项，那么可以使用下列方程进行计算，从而避免精度误差：

$$\begin{array}{rl}4q_x^2& =+m_{00}-m_{11}-m_{22}+m_{33},\\ 4q_y^2& =-m_{00}+m_{11}-m_{22}+m_{33},\\ 4q_z^2& =-m_{00}-m_{11}+m_{22}+m_{33},\\ 4q_w^2& =\mathrm{tr}\left({\mathbf{M}}^q\right).\end{array}\text{\hspace{1em}(4.51)}$$

  总而言之，方程4.50可以用来计算四元数虚部分量$q_x,q_y,q_z,q_w$中的最大值，根据最大值的不同适当选取方程4.51或者方程4.47来计算$\hat{\mathbf{q}}$的其余分量。 Schuler 提出了一种无分支的变体方法，但是该方法计算了四次平方根。

球面线性插值

  球面线性插值（spherical linear interpolation）是指，给定两个单位四元数$\hat{\mathbf{q}}$，$\hat{\mathbf{r}}$以及一个参数$t\in [0,1]$，然后计算出一个插值的四元数。球面线性插值对于物体动画非常有用，但是对于相机的方向插值来说，却并不常用，因为相机的up向量在插值的过程中可能会发生倾斜，这通常是一个干扰的效果。

  球面线性插值的数学表达式如下：

$$\hat{\mathbf{s}}(\hat{\mathbf{q}},\hat{\mathbf{r}},t)={\left(\hat{\mathbf{r}}{\hat{\mathbf{q}}}^{-1}\right)}^t\hat{\mathbf{q}}\text{\hspace{1em}(4.52)}$$

  但是，对于软件实现而言，通常会使用更为合适的slerp函数来代表球面线性插值，其形式如下：

$$\hat{\mathbf{s}}(\hat{\mathbf{q}},\hat{\mathbf{r}},t)=\mathrm{slerp}(\hat{\mathbf{q}},\hat{\mathbf{r}},t)=\frac{\sin (\varphi (1-t))}{\sin \varphi }\hat{\mathbf{q}}+\frac{\sin (\varphi t)}{\sin \varphi }\hat{\mathbf{r}}\text{\hspace{1em}(4.53)}$$

  其中$\cos \varphi =q_x{r}_x+q_y{r}_y+q_z{r}_z+q_w{r}_w$。对于给定的$t\in [0,1]$，slerp函数所计算出的四元数是唯一的（当且仅当$\hat{\mathbf{q}}$不与$\hat{\mathbf{r}}$反向。）。这些插值生成的四元数组合在一起，构成了一个四维单位球面上，从$\hat{\mathbf{q}}(t=0)$到$\hat{\mathbf{r}}(t=1)$的最短弧。如图4.10所示，四元数$\hat{\mathbf{q}}$，$\hat{\mathbf{r}}$与原点一起构成了一个平面，这个平面与四维单位超球相交，投影到三维空间中便形成了一个圆，而这些插值出来的四元数所形成的圆弧便位于这个圆上。随着参数$t$的均匀变化，这个插值出来的四元数会绕着一个固定的轴，以一个恒定的角速度进行旋转；像这样一个具有恒定变化速度（即加速度为0）的曲线被叫做测地线（geodesic curve）。球面上的大圆（great circle）是由经过原点的平面与球体相交形成的，这个圆上的一部分被叫做大弧（great arc）。

图4.10 单位四元数可以表示为单位球面（实际上是一个四维空间中的超球，投影到了三维空间中）上的点。slerp函数用于对两个四元数进行插值，插值形成的路径是球面上的一段圆弧。需要注意的是，从${\hat{\mathbf{q}}}_1$插值到${\hat{\mathbf{q}}}_2$，与从${\hat{\mathbf{q}}}_1$插值到${\hat{\mathbf{q}}}_3$，再插值到${\hat{\mathbf{q}}}_2$并不是一回事。虽然它们插值的起点和终点是相同的，但是中间具体的路径和方向是不同的。

  slerp函数非常适合用于在两个方向之间进行插值，这个插值过程会围绕一个固定的旋转轴，且旋转速度是恒定的；而在使用欧拉角进行方向插值的时候则并不是这样的（旋转轴不固定，旋转速度也不恒定）。在实际实现中，由于slerp函数中涉及到了很多三角函数的计算，因此效率是很低的。Malyshau 讨论了如何将四元数集成到渲染管线中，他指出：当直接在像素着色器中对四元数进行归一化，而不是使用slerp函数来调整三角形的朝向时，在90度角内的最大误差为4度，这个误差在光栅化三角形的时候是可以接受的。Li 提供了一种效率更高的增量方法来计算slerp，并且不会对精度产生任何影响。Eberly 提出了一种仅使用加法和乘法来快速计算slerp的方法。

  当现在有两个以上的方向时，例如${\hat{\mathbf{q}}}_0,{\hat{\mathbf{q}}}_1,\dots ,{\hat{\mathbf{q}}}_{n-1}$，希望可以从${\hat{\mathbf{q}}}_1$插值到${\hat{\mathbf{q}}}_2$，再插值到${\hat{\mathbf{q}}}_3$，直到最终插值到${\hat{\mathbf{q}}}_{n-1}$，使用slerp函数是最直接的方式。例如现在我们到达了${\hat{\mathbf{q}}}_{i-1}$，正要向${\hat{\mathbf{q}}}_i$方向前进，我们将使用${\hat{\mathbf{q}}}_i$，${\hat{\mathbf{q}}}_{i-1}$来作为slerp的参数进行插值；当到达了${\hat{\mathbf{q}}}_i$，要向${\hat{\mathbf{q}}}_{i+1}$方向前进时，我们使用${\hat{\mathbf{q}}}_i$，${\hat{\mathbf{q}}}_{i+1}$作为slerp的参数进行插值。但是这样的操作会导致在插值过程中的方向发生突变，如图4.10所示；这和点的线性插值情况相类似。

  一个更好的方法是使用某种样条线来进行插值。在四元数${\hat{\mathbf{q}}}_i$与${\hat{\mathbf{q}}}_{i+1}$之间引入了2个新的中间四元数${\hat{\mathbf{a}}}_i$，${\hat{\mathbf{a}}}_{i+1}$。我们可以使用${\hat{\mathbf{q}}}_i,{\hat{\mathbf{a}}}_i,{\hat{\mathbf{a}}}_{i+1},{\hat{\mathbf{q}}}_{i+1}$这四个四元数，来定义一个球面立方插值（Spherical cubic interpolation）。这两个中间四元数可以通过下面的方程进行计算（Shoemake 给出了另一种推导方式）：

$${\hat{\mathbf{a}}}_i={\hat{\mathbf{q}}}_i\exp \left[-\frac{\log \left({\hat{\mathbf{q}}}_i^{-1}{\hat{\mathbf{q}}}_{i-1}\right)+\log \left({\hat{\mathbf{q}}}_i^{-1}{\hat{\mathbf{q}}}_{i+1}\right)}{4}\right].\text{\hspace{1em}(4.54)}$$

  使用${\hat{\mathbf{q}}}_i,{\hat{\mathbf{a}}}_i$来构建一个光滑的三次样条线，从而对四元数进行光滑的球面插值：

$$\begin{gathered} \mathbf{squad}\left({\hat{\mathbf{q}}}_i,{\hat{\mathbf{q}}}_{i+1},{\hat{\mathbf{a}}}_i,{\hat{\mathbf{a}}}_{i+1},t\right)= \\ \mathbf{slerp}\left(\mathbf{slerp}\left({\hat{\mathbf{q}}}_i,{\hat{\mathbf{q}}}_{i+1},t\right),\mathbf{slerp}\left({\hat{\mathbf{a}}}_i,{\hat{\mathbf{a}}}_{i+1},t\right),2t(1-t)\right)\text{\hspace{1em}(4.55)} \end{gathered}$$

  从上述方程中可以看到，$\mathbf{squad}$函数是通过重复使用slerp进行球面线性插值构建而成的。这个插值形成的路径，会通过每一个${\hat{\mathbf{q}}}_i$，但是并不会通过这些用于构建样条线的${\hat{\mathbf{a}}}_i$——它们用于指示位于初始位置时的切线方向。

将一个向量旋转到另一个向量

  一个常见的操作是将方向$\mathbf{s}$以最短路径变换到另一个方向$\mathbf{t}$，使用四元数可以大大简化这个变换过程，这也显示了四元数与其表示方法的密切关系。首先要将$\mathbf{s}$和$\mathbf{t}$进行归一化；其次通过$\mathbf{u}=(\mathbf{s}\times \mathbf{t})/\Vert \mathbf{s}\times \mathbf{t}\Vert$来计算归一化的旋转轴（同时垂直于向量$\mathbf{s}$和$\mathbf{t}$）；然后根据$e=\mathbf{s}\cdot \mathbf{t}=\cos (2\varphi )，\Vert \mathbf{s}\times \mathbf{t}\Vert =\sin (2\varphi )$，其中$2\varphi$是向量$\mathbf{s}$和$\mathbf{t}$之间的夹角。这时候就可以表示从$\mathbf{s}$变换到$\mathbf{t}$的四元数了，即$\hat{\mathbf{q}}=(\sin \varphi \mathbf{u},\cos \varphi )$；我们将其展开可得$\hat{\mathbf{q}}=\left(\frac{\sin \varphi }{\sin 2\varphi }(\mathbf{s}\times \mathbf{t}),\cos \varphi \right)$，又因为$\cos \varphi =\sqrt{\frac{e+1}{2}}$，进一步化简可得：

$$\hat{\mathbf{q}}=\left({\mathbf{q}}_v,q_w\right)=\left(\frac{1}{\sqrt{2(1+e)}}(\mathbf{s}\times \mathbf{t}),\frac{\sqrt{2(1+e)}}{2}\right)\text{\hspace{1em}(4.56)}$$

  使用这种方式构建四元数（相比于对叉积$\mathbf{s}\times \mathbf{t}$的结果进行标准化），避免了当$\mathbf{s}$和$\mathbf{t}$几乎指向相同方向时造成的数值不稳定。但是当$\mathbf{s}$和$\mathbf{t}$指向相反方向的时候，此时$e=-1$，两种方法都会出现数值不稳定的情况，因为分母为0；当检测到这种特殊情况出现的时候，可以使用任何一个垂直于$\mathbf{s}$的向量作为旋转轴，然后再将$\mathbf{s}$旋转到$\mathbf{t}$。

  有时可能会需要一个表示$\mathbf{s}$旋转到$\mathbf{t}$的矩阵，上文中的方程4.46给出了单位四元数所表示的旋转矩阵，可以从这个矩阵出发，将方程4.56带入其中并进行化简，可得：

$$\mathbf{R}(\mathbf{s},\mathbf{t})=\left(\begin{array}{cccc}e+h{v}_x^2& h{v}_x{v}_y-v_z& h{v}_x{v}_z+v_y& 0\\ h{v}_x{v}_y+v_z& e+h{v}_y^2& h{v}_y{v}_z-v_x& 0\\ h{v}_x{v}_z-v_y& h{v}_y{v}_z+v_x& e+h{v}_z^2& 0\\ 0& 0& 0& 1\end{array}\right)\text{\hspace{1em}(4.57)}$$

  在这个计算过程中，会使用到以下的中间计算过程：

$$\begin{array}{rl}\mathbf{v}& =\mathbf{s}\times \mathbf{t}\\ e& =\cos (2\varphi )=\mathbf{s}\cdot \mathbf{t},\\ h& =\frac{1-\cos (2\varphi )}{{\sin }^2(2\varphi )}=\frac{1-e}{\mathbf{v}\cdot \mathbf{v}}=\frac{1}{1+e}.\end{array}\text{\hspace{1em}(4.58)}$$

  可以看到，通过对方程进行简化，避免了平方根运算和三角函数运算，因此这是一种效率很高的计算方式。请注意方程4.57（使用四元数将一个向量旋转到另一个向量）和方程4.30（使用欧拉变换将一个向量旋转到另一个向量），二者的结构是类似的，但是使用四元数变换的方程4.57中，并没有涉及对三角函数的计算，因此其效率更高。

  当向量$\mathbf{s}$和$\mathbf{t}$平行或者接近平行的时候，此时$\Vert \mathbf{s}\times \mathbf{t}\Vert \approx 0$，需要格外小心。当$\varphi \approx 0$的时候，可以返回单位矩阵；当时当$2\varphi \approx \pi$的时候，可以绕任意轴旋转$\pi$弧度，这个轴可以是向量$\mathbf{s}$与任何其他不平行于$\mathbf{s}$的向量之间的叉乘。Moller和Hughes使用了Householder矩阵（初等反射矩阵，根据一个向量和平面，生成一个反射的向量），来以一种不同的方法来处理这个特殊情况。

顶点混合

  想象现在有一个数字人的手臂，他由两部分组成，分别是前臂和上臂，这个模型可以使用刚体变换来进行动画处理，如图4.11左半边所示。但是这两部分之间的关节并不像一个真正的手臂肘关节，这是因为使用了两个独立的子模型，即这个肘关节由两个独立模型的重叠部分所组成。这个由静态模型组合而成的物体，并没有解决如何使得关节灵活柔韧（flexible）的问题，很明显，使用一个单独的物体进行动画处理，可以获得更好的效果。

  顶点混合（vertex blending）是一种解决这个问题的常见方法，它还有几个其他的名称，例如：线性蒙皮混合（linear-blend skinning）、包络（enveloping）或者骨架子空间变形（skeleton-subspace deformation）。虽然该算法的确切来源尚不清楚，但是定义骨骼，并让模型皮肤随着骨骼的变化产生相应的反应则是计算机动画中的基本概念。在最简单的形式中，数字人的前臂和上臂可以像以前一样分别进行动画处理，但是在关节处，这两个独立的部分会通过一个具有弹性的“皮肤（skin）”连接起来。这部分皮肤由两组不同的顶点组成，其中一部分顶点由前臂的变换矩阵所驱动，另一部分顶点则被上臂的变换矩阵所驱动。这意味着位于关节处的“皮肤”，其三角形顶点不再使用统一的矩阵进行变换，而是每个三角形都使用不同的矩阵来进行变换，如图4.11所示。

图4.11 左侧的手臂由前臂和上臂两部分组成，并使用了刚体变换进行动画处理，其关节处表现的并不是很真实自然。在右侧，对整个手臂模型使用了顶点混合方法。其中，位于中间的手臂展示了，当一个连接两部分的简单皮肤直接覆盖在肘关节时会发生什么。最右侧的手臂展示了使用顶点混合进行处理的结果，其中一些顶点使用了不同的权重进行混合：$(2/3,1/3)$意味着应用在该顶点上的变换有$2/3$来自上臂，有$1/3$来自前臂。在最右侧的图中，还展示了使用顶点混合的一个缺点，即可以看见肘关节内部的折叠。还可以使用更多的骨骼节点，以及更加合理的权重选择来获得更好的动画效果。

  更进一步说，可以允许单个顶点被若干个不同的矩阵所变换，并将这些变换所得到的结果加权混合在一起。这是通过为动画物体设置一个骨架来完成的，通过用户自定义的权重，每个骨骼的变换都可能会对其他的顶点产生影响。由于整个手臂都可能是“有弹性的”，也就是说，所有的顶点都可能会受到多个矩阵变换的影响，因此也可以将整个网格模型称作为包裹骨骼的蒙皮（skin），如图4.12所示。许多商业化的建模软件都具有这样的骨架-骨骼建模功能，尽管名字中含有骨骼，但是这个骨骼也不一定是刚性的。例如：Mohr和Gleicher 提出了一种想法，可以通过增加额外关节的方法，来实现诸如肌肉膨胀的效果。James和Twigg 讨论了如何使用可拉伸与可挤压骨骼来进行动画蒙皮。

图4.12 顶点混合的一个真实案例。左上角的图展示了这个手臂的两个骨骼，此时手臂处于伸展状态。右上角的图展示了手臂的顶点和网格模型，其中使用颜色标注了该顶点属于哪个骨骼。在底部的图中，被着色处理的手臂处于一个稍微不同的位置，呈弯曲状态。

  顶点混合的数学表达式如方程4.59所示，其中$\mathbf{p}$是原始的顶点坐标，$\mathbf{u}(t)$是变化后的顶点位置，$t$代表了动画时间。方程的具体形式如下：

$$\mathbf{u}(t)=\sum _{i=0}^{n-1}{w}_i{\mathbf{B}}_i(t){\mathbf{M}}_i^{-1}\mathbf{p},\text{ where}\sum _{i=0}^{n-1}{w}_i=1,w_i\ge 0\text{\hspace{1em}(4.59)}$$

  上述方程的含义是：有$n$个骨骼会对点$\mathbf{p}$的世界空间坐标产生影响；$w_i$的值代表了骨骼$i$对于顶点$\mathbf{p}$坐标的影响权重；矩阵${\mathbf{M}}_i^{-1}$代表骨骼$i$的模型变换矩阵，它负责将骨骼变换到世界空间中。通常来说，骨骼的控制关节位于自身坐标系的原点位置，例如对于前臂骨骼而言，其控制关节为肘关节，则肘关节应当位于这个前臂骨骼模型空间的原点处；这样可以直接对模型的旋转矩阵进行动画处理，从而可以对整个前臂进行移动，并且移动结束之后，肘关节仍然位于模型空间的原点处。${\mathbf{B}}_i(t)$是一个世界变换矩阵，代表了第$i$个骨骼的世界变换，它会随着时间进行变化，从而让物体动起来；它通常是由若干个矩阵连接而成的，例如上一层级（父级）的变换以及局部的动画矩阵等。

  Woodland对一种维护和更新动画矩阵${\mathbf{B}}_i(t)$的方法进行了深入讨论。首先找到所有会对顶点产生影响的骨骼，然后每个骨骼都会将这个顶点坐标转换到自身的坐标系中，从而生成一个变换后的新顶点，而该顶点的最终位置，是通过在这些新顶点之间进行插值获得的。在一些关于蒙皮的讨论中，模型变换矩阵${\mathbf{M}}_i$被认为是世界变换矩阵${\mathbf{B}}_i(t)$的一部分，因此并不会将${\mathbf{M}}_i$显式展示出来。在这里将其展示了出来，因为这是一个十分有用的矩阵，它总会成为矩阵连接过程的一部分。

  实际上，矩阵${\mathbf{B}}_i(t)$和矩阵${\mathbf{M}}_i^{-1}$在每一帧中都会进行连接组合，最终生成的变换矩阵会被用于对顶点进行变换。点$\mathbf{p}$会被不同骨骼的变换矩阵进行变换，并通过使用权重$w_i$，来将这些变换进行加权混合——这也是顶点混合的名称由来。这些权重的值通常都是非负的，且它们的和为1，即点$\mathbf{p}$会被变换到若干个位置上，并在这些位置之间进行插值。换句话说，对于一个固定的时刻$t$，最终插值出来的点$\mathbf{u}$位于点集${\mathbf{B}}_i(t){\mathbf{M}}_i^{-1}\mathbf{p},i=0\cdots n-1$的凸包（convex hull）内部。通常我们也可以使用方程4.59来对顶点的法线进行变换，当然这也取决于所使用的变换类型：如果骨骼被大量拉伸或者挤压（非均匀缩放）的话，就需要使用矩阵${\mathbf{B}}_i(t){\mathbf{M}}_i^{-1}$的逆来对法线进行变换。

  顶点混合非常适合在GPU上运行，可以将网格中的顶点集合放入GPU的静态缓存中（只需要发送一次），之后重复使用这个缓存即可。在每一帧中，只有骨骼节点发生了变化（顶点数量和连接关系不会发生变化，即GPU上的静态缓存也不会发生变化），使用一个顶点着色器，来计算骨骼节点对缓存中存储网格的影响。通过这种方式，可以最大程度的减少在CPU上进行的处理，以及与GPU之间的数据传输（会带来大量延迟），这使得GPU可以高效地渲染这个网格。最简单的情况是可以同时使用模型的整套骨骼矩阵，否则的话就需要拆分模型，并复制一些骨骼。还可以将骨骼的变换存储在顶点可以访问的纹理中，从而避免触及寄存器的限制。通过使用四元数来表示旋转变换，那么每个变换都可以被存储在两张纹理中（一张存储顶点的位置，一张存储用于旋转的单位四元数）。如果可以使用无序访问视图（UAV）的话，还可以对蒙皮的动画结果进行重复使用。

  还可以指定超出$[0,1]$范围之外的权重值，或者是让权重之和不为1。当然这只在使用一些其他的混合算法时才有意义。

  刚才介绍的是最基础的顶点混合算法，它当然会存在一些缺点，例如不希望出现的折叠、扭曲或者是自相交等现象，如图4.13所示。一个更好的方案是使用双四元数混合（dual quaternion），这个技术可以让蒙皮保持原始形状的刚性，从而避免四肢出现像“糖纸”一样的扭曲；这个技术的效果良好，并且计算成本小于线性蒙皮混合的1.5倍，因此被广泛采用。但是双四元数混合技术会导致蒙皮发生膨胀，Le和Hodings 提出了中心旋转蒙皮（center-of-rotation）作为双四元数的替代，他们认为：局部变换应当是一种刚体变换，并且具有相似权重$w_i$的顶点，应当也具有相似的变换。这种技术预先计算了每个顶点的旋转中心，并对其施加了正交（刚体）约束，从而防止肘部折叠以及糖纸弯曲现象。在运行时，该技术类似于线性蒙皮混合，它运行在GPU上，在顶点的旋转中心进行线性蒙皮混合，紧接着再进行四元数混合的操作。

图4.13 左图展示了在使用线性蒙皮混合算法时关节会出现的问题；右图则使用了双四元数混合算法，它明显改善了关节处的表现。

变形

  在动画中，从一个三维模型变形到另一个三维模型是非常有用的。假设在$t_0$时刻有一个正在展示的模型，希望它能够在$t_1$时刻变形为另一个模型。在$t_0$和$t_1 $时刻之间，我们会使用某种插值方法，得到一个不断变化的“混合”模型。图4.14展示了变形的一个例子。

图4.14 顶点变形（vertex morphing）。每个顶点上都定义了两组位置和法线，在每一帧中，顶点着色器会对这两组位置和法线进行线性插值，从而获得中间的位置和法线。

  变形主要会涉及到两个问题：顶点对应（vertex correspondence）问题和插值（interpolation）问题。给定两个任意的模型，它们具有不同的拓扑结构，不同数量的顶点以及不同的网格连通性。通常要从设置两个模型顶点之间的对应关系开始，这是一个十分困难的问题。

  如果两个模型之间已经有了一一对应的顶点了，即对于第一个模型中的每个顶点，都可以在第二个模型中找到唯一相对应的顶点。这种情况下的变形操作是很简单的，只需要在对应的顶点属性之间进行插值就可以获得中间数据了，例如直接使用线性插值。对于$t\in [t_0,t_1]$，为了计算这个变形的顶点，首先我们需要计算$s=\left(t-t_0\right)/\left(t_1-t_0\right)$（即当前时刻$t$占总时间的比例），然后将两个对应顶点进行线性混合：

$$\mathbf{m}=(1-s){\mathbf{p}}_0+s{\mathbf{p}}_1,\text{\hspace{1em}(4.60)}$$

  其中${\mathbf{p}}_0$与${\mathbf{p}}_1$对应了同一个顶点不同时刻（$t_0,t_1$）的状态。

  变形目标（morph target）或者形状混合（blend shape）是变形技术的其中一种，用户可以对其进行更加直观的控制，图4.15展示了其基本思想。

图4.15 给定两个嘴巴的姿态，可以计算一组差分向量，用于控制插值。在变形目标中，我们可以通过调整差分向量的权重，从而在中性面部（标准面部姿态，类似于T Pose）上“添加”运动。在这个例子中，当权重为正的时候，我们可以获得一个嘴角上扬的笑脸；当权重为负的时候，我们可以获得一个相反的结果，即嘴角下垂的哭脸。

  从一个标准的中性模型（neutral model）开始，在这个例子中，这个中性模型是一个没有表情的脸，将其标记为$\mathcal{N}$。此外还需要有一组不同的面部姿态，而在这个例子中只有一个笑脸作为额外的面部姿态。一般来说，可以允许$k\ge 1$个不同的姿态，我们使用${\mathcal{P}}_i,i\in [1,\dots ,k]$来标记它们。在预处理阶段，通过使用每个不同的姿态来减去中性姿态，从而计算出姿态之间的差分向量，即${\mathcal{D}}_i={\mathcal{P}}_i-\mathcal{N}$。

  此时有了一个中性模型以及一组不同姿态${\mathcal{D}}_i$，可以使用下列方程计算得到一个变形的模型$\mathcal{M}$：

$$\mathcal{M}=\mathcal{N}+\sum _{i=1}^k{w}_i{\mathcal{D}}_i\text{\hspace{1em}(4.61)}$$

  在中性模型的基础上，通过使用不同的权重$w_i$，根据需要来将不同姿态的特征进行加权混合。在图4.15中，如果设定此时的权重为$w_1=1$，那么就会获得一个笑脸，如中间的面部表情；如果设定此时的权重为$w_1=0.5$，那么就可以获得一个半微笑的脸，以此类推。同时还可以使用一个负数权重，或者是使用和不为1的权重。

  对于这个简单的面部模型而言，还可以再添加一个带有“悲伤”眉毛的面部姿态，通过使用一个负数权重，从而获得一个“开心”的眉毛。由于这些顶点的位移都是可叠加且相互独立的，因此这个眉毛姿态可以和和嘴巴姿态同时使用。

  变形目标是一个强大的技术，它为动画师提供了很多控制能力，因为模型中的不同特征，可以进行独立操控。Lewis等人介绍了姿态空间变形技术（pose-space deformation），该技术将顶点混合和变形目标结合在一起。Senior 使用了预计算的顶点纹理来存储和检索目标姿态之间的位移。对于支持流式输出和顶点ID的硬件，允许在单个模型上使用多个变形目标，并在GPU上专门计算这个效果。此外，还可以使用一个低分辨率的网格，然后通过曲面细分和位移映射（displacement mapping）来生成一个高分辨率的网格，这样避免了在一个高分辨率模型上，对每个顶点进行蒙皮和处理的开销。

  图4.16展示了一个同时使用蒙皮和变形的例子。Weronko和Andreason 在《教团：1886》中使用了蒙皮和变形。

图4.16 《声名狼藉：次子》的游戏画面，角色Delsin的面部使用了形状混合来制作表情动画。所有的这些面部表情，都是基于一个没有表情的面部（中性），然后通过修改不同的权重，从而使得面部表情看起来不同。

几何缓存回放

  在一些过场动画中，可能会希望使用一些极高质量的动画，这些高质量动画使用上述所提到的任何方法和技术都无法实时实现。一种简单的方法是，将所有帧的顶点数据预先计算并存储起来， 然后在游戏运行的时候，从磁盘上读取这些数据并对顶点进行更新。但是对于一个包含30000个顶点的模型而言，一段简单的动画可能就会占据50 MB/s的带宽。Gneiting 提供了若干种方法，可以将内存消耗降低到原来的10%，其具体做法如下：

  首先会对数据进行量化处理，例如：使用16 位整数来存储每个顶点的位置和纹理坐标，这一步操作会损失一些精度和信息，因为这个压缩过程是不可逆的，在压缩之后无法恢复原始的数据。为了进一步对数据进行压缩，还使用了基于空域（spatial）和时域（temporal）的预测，并对这些差异信息（同一顶点的不同帧，同一帧中的相邻顶点）进行了编码。对于空域压缩而言，可以使用平行四边形预测（parallelogram prediction）技术。对于三角形条带上的一个顶点，我们可以使用该顶点所在的三角形，以当前边为基准，将其反射到当前边的另一侧；这两个三角形构成了一个平行四边形，然后使用这个平行四边形来预测下一个顶点的位置，然后对原顶点和新顶点之间的差异进行编码。在良好预测的情况下，大多数差异数值都会接近于0，这对于许多压缩方案来说都是很理想的。与MPEG压缩类似，还可以在时域上进行预测，即每隔$n$帧进行一次空域压缩。还可以在时域上对两帧之间进行预测，例如：如果一个顶点被一个增量向量从第$n-1$帧移动到了第$n$帧，那么它很有可能会以类似的增量，再次移动到第$n+1$帧。这些技术大大减少了所需要的存储空间，从而使得动画系统可以实时使用流式数据。