四元数

定义与运算

$\hat{\mathbf{q}}=(\mathbf{q}_v,q_w)=iq_x+jq_y+kq_z+q_w$，$i^2=j^2=k^2=-1$，$jk=-kj=i$，$ki=-ik=j$，$ij=-ji=k$。

其他规则：$(\hat{\mathbf{q}}\hat{\mathbf{r}})^*=\hat{\mathbf{r}}^*\hat{\mathbf{q}}^*$，$n(\hat{\mathbf{q}}\hat{\mathbf{r}})=n(\hat{\mathbf{q}})n(\hat{\mathbf{r}})$，乘法满足结合律与线性。

单位四元数可写为 $\hat{\mathbf{q}}=(\sin\phi\,\mathbf{u}_q,\cos\phi)=e^{\phi\mathbf{u}_q}$，于是 $\log(\hat{\mathbf{q}})=\phi\mathbf{u}_q$，$\hat{\mathbf{q}}^t=\sin(\phi t)\mathbf{u}_q+\cos(\phi t)$。

表示旋转

• 把点或向量放进四元数 $\hat{\mathbf{p}}$，$\hat{\mathbf{q}}\hat{\mathbf{p}}\hat{\mathbf{q}}^{-1}$ 绕轴 $\mathbf{u}_q$ 旋转 $2\phi$；单位四元数的逆就是共轭。
• 任意非零实数倍表示同一旋转，所以 $\hat{\mathbf{q}}$ 与 $-\hat{\mathbf{q}}$ 等价；从矩阵提取四元数可能得到其中任一个。
• 先 $\hat{\mathbf{q}}$ 再 $\hat{\mathbf{r}}$：$\hat{\mathbf{c}}=\hat{\mathbf{r}}\hat{\mathbf{q}}$。

单位四元数通过夹心积把向量绕轴旋转，步骤如下（角度为 $2\phi$）。

flowchart LR
  V["向量 v"] --> P["写入纯四元数 p = 0,v"]
  P --> Mul["计算 q p q⁻¹"]
  Q["单位四元数 q"] --> Mul
  Mul --> Out["取出虚部得 v'"]

与矩阵互转

• 单位四元数到矩阵（Eq. 4.46）只需乘加，无三角函数，例如左上元素为 $1-2(q_y^2+q_z^2)$。
• 矩阵到四元数：$q_w=\frac12\sqrt{\operatorname{tr}(\mathbf{M}^q)}$，$q_x=\frac{m_{21}-m_{12}}{4q_w}$ 等。为数值稳定、避免除以小数，先比较 $m_{00},m_{11},m_{22}$ 与迹，确定 $q_x,q_y,q_z,q_w$ 中最大者，先算它再算其余分量。也有无分支但用四次开方的变体。

何时把四元数转成矩阵：需要批量变换大量向量或与仿射矩阵连乘时；插值与串接旋转则留在四元数空间更划算。

flowchart TD
  Need["当前任务"] --> Q1{"要插值或串接旋转？"}
  Q1 -->|"是"| Keep["留在四元数"]
  Q1 -->|"否"| Q2{"要批量乘大量向量 / 与 T·S 合并？"}
  Q2 -->|"是"| Mat["转为矩阵再乘"]
  Q2 -->|"否"| Keep

球面线性插值 (slerp)

$\text{slerp}(\hat{\mathbf{q}},\hat{\mathbf{r}},t)=\frac{\sin(\phi(1-t))}{\sin\phi}\hat{\mathbf{q}}+\frac{\sin(\phi t)}{\sin\phi}\hat{\mathbf{r}},\quad \cos\phi=q_xr_x+q_yr_y+q_zr_z+q_wr_w$

• 代数形式为 $(\hat{\mathbf{r}}\hat{\mathbf{q}}^{-1})^t\hat{\mathbf{q}}$。
• 结果是四维单位球面上的最短大圆弧（测地线），对应绕固定轴匀速旋转。 下图对比 slerp 沿单位球大圆最短弧前进，与先线性插值再归一化（nlerp/lerp）走弦的差别：前者角速度恒定，后者在两端更快、中间更慢。

flowchart TD
  Start["起点朝向 q"] --> Choice{"插值方式"}
  Choice -->|"slerp"| Arc["沿大圆最短弧匀速"]
  Choice -->|"lerp 再归一化"| Chord["沿弦采样再投回球面"]
  Arc --> End1["终点朝向 r：角速度恒定"]
  Chord --> End2["终点朝向 r：两端快、中间慢"]

• 适合动画物体；不太适合相机朝向插值，因为相机的"上"向量可能在插值中倾斜。
• 直接计算 slerp 要调用三角函数，较贵。书中提到：在像素着色器里只做归一化代替 slerp，90° 时误差最大 4°，对光栅化三角形可以接受；也有不损精度的增量算法和只用加乘的快速算法。

多个关键帧：squad

逐段 slerp 在经过关键帧时会突然变向。用样条更好：在 $\hat{\mathbf{q}}_i$ 与 $\hat{\mathbf{q}}_{i+1}$ 之间引入 $\hat{\mathbf{a}}_i=\hat{\mathbf{q}}_i\exp\left[-\frac{\log(\hat{\mathbf{q}}_i^{-1}\hat{\mathbf{q}}_{i-1})+\log(\hat{\mathbf{q}}_i^{-1}\hat{\mathbf{q}}_{i+1})}{4}\right]$ $\text{squad}(\hat{\mathbf{q}}_i,\hat{\mathbf{q}}_{i+1},\hat{\mathbf{a}}_i,\hat{\mathbf{a}}_{i+1},t)=\text{slerp}(\text{slerp}(\hat{\mathbf{q}}_i,\hat{\mathbf{q}}_{i+1},t),\text{slerp}(\hat{\mathbf{a}}_i,\hat{\mathbf{a}}_{i+1},t),2t(1-t))$ 插值经过各 $\hat{\mathbf{q}}_i$，但不经过 $\hat{\mathbf{a}}_i$（它们表示切向朝向）。

从一个方向转到另一个方向

单位向量 $\mathbf{s}$ 转到 $\mathbf{t}$，$e=\mathbf{s}\cdot\mathbf{t}=\cos 2\phi$： $\hat{\mathbf{q}}=\left(\frac{1}{\sqrt{2(1+e)}}(\mathbf{s}\times\mathbf{t}),\ \frac{\sqrt{2(1+e)}}{2}\right)$

• 这样直接构造（而不是先归一化 $\mathbf{s}\times\mathbf{t}$）在两者几乎同向时数值稳定。
• 两者反向时会除以零，此时取任意垂直于 $\mathbf{s}$ 的轴旋转。
• 也可直接写出对应矩阵（Eq. 4.57，$h=1/(1+e)$），不含开方和三角函数；接近平行时，$\phi\approx0$ 返回单位矩阵，$2\phi\approx\pi$ 绕任意垂直轴转 $\pi$。

对比

GEA 的视角：游戏引擎中的四元数（卷 I 5.4）

GEA 第 5 章从游戏引擎实践出发讲四元数，结论与上文 RTR4 一致，但在记法、乘法顺序的表述、矩阵转换代码和 lerp/slerp 取舍上有补充。GEA 使用行向量约定。

为什么不用矩阵表示旋转（卷 I 5.4）

1. 9 个浮点数表示只有 3 个自由度的旋转，太浪费；
2. 旋转向量要做向量-矩阵乘法：三个点积，共 9 次乘法和 6 次加法；
3. 很难在两个矩阵表示的朝向之间求中间旋转（如相机在几秒内从朝向 A 平滑转到 B）。

原书问题：GEA 说四元数「克服了这三个问题」，但第 2 条（旋转单个向量的开销）并不成立：用四元数直接旋转一个向量，即使用优化公式 $\mathbf{v}'=\mathbf{v}+2\mathbf{r}\times(\mathbf{r}\times\mathbf{v}+w\mathbf{v})$，也要约 15 次乘法 + 15 次加法，而 3×3 矩阵只要 9 次乘法 + 6 次加法；要批量旋转大量向量，通常先把四元数转成矩阵。四元数真正便宜的是串接旋转（16 乘 12 加，矩阵为 27 乘 18 加）、存储和插值。来源：Wikipedia「Quaternions and spatial rotation」§Performance comparisons（引 D. Eberly, Rotation Representations and Performance Issues）。 四元数由 Hamilton 于 1843 年作为复数的推广提出，可视为一个实轴加三个虚轴 $i,j,k$ 的「四维复数」，复数形式 $\mathbf{q}=iq_x+jq_y+kq_z+q_w$；技术上它服从实数上的四维赋范可除代数。游戏中只需知道：单位长度四元数表示三维旋转。

单位四元数与轴角（卷 I 5.4.1）

$\mathbf{q}=[\mathbf{q}_v\ \ q_s]=[\mathbf{a}\sin\tfrac{\theta}{2}\ \ \cos\tfrac{\theta}{2}]$ $\mathbf{a}$ 为单位旋转轴，$\theta$ 为转角，方向遵循右手定则。它很像轴角表示 $[\mathbf{a}\ \theta]$，但数学上更方便。这与上文 RTR4「$\hat{\mathbf{q}}\hat{\mathbf{p}}\hat{\mathbf{q}}^{-1}$ 旋转 $2\phi$」是同一件事：RTR4 的 $\phi$ 即 GEA 的半角 $\theta/2$。

运算（卷 I 5.4.2）

• 两个单位四元数之和不是单位四元数，不表示旋转；引擎里看不到未经缩放的四元数和。
• Grassman 积（原书拼写；人名应为 Grassmann）：$\mathbf{pq}=[(p_s\mathbf{q}_v+q_s\mathbf{p}_v+\mathbf{p}_v\times\mathbf{q}_v)\ \ (p_sq_s-\mathbf{p}_v\cdot\mathbf{q}_v)]$，向量部分放在 x/y/z、标量部分在 w。四元数乘法有多种，GEA 只讨论与三维旋转配合的这一种。$\mathbf{pq}$ 表示先 Q 后 P 的复合旋转。
• 共轭 $\mathbf{q}^*=[-\mathbf{q}_v\ q_s]$，逆 $\mathbf{q}^{-1}=\mathbf{q}^*/|\mathbf{q}|^2$；单位四元数时逆等于共轭，免去除法——求四元数的逆通常比求 3×3 矩阵的逆快得多，可用于优化。$[0\ 0\ 0\ 1]$ 表示零旋转。
• $(\mathbf{pq})^*=\mathbf{q}^*\mathbf{p}^*$，$(\mathbf{pq})^{-1}=\mathbf{q}^{-1}\mathbf{p}^{-1}$，与矩阵乘积的转置、求逆反序类似。

旋转向量与串接（卷 I 5.4.3）

• 向量写成标量部分为 0 的四元数 $[\mathbf{v}\ 0]$，$\mathbf{v}'=\mathbf{q}\mathbf{v}\mathbf{q}^{-1}=\mathbf{q}\mathbf{v}\mathbf{q}^*$。
• 例：游戏约定模型空间前向 $\mathbf{f}_M=[0\ 0\ 1]$，飞机朝向四元数 $\mathbf{q}$ 给出世界前向 $\mathbf{f}_W=\mathbf{q}[0\ 0\ 1\ 0]\mathbf{q}^{-1}$。
• 串接顺序：依次施加旋转 1、2、3，矩阵（行向量）为 $\mathbf{R}_{net}=\mathbf{R}_1\mathbf{R}_2\mathbf{R}_3$，四元数却是 $\mathbf{q}_{net}=\mathbf{q}_3\mathbf{q}_2\mathbf{q}_1$——因为四元数总在向量两侧相乘，未取逆的在左、从右往左读，取逆的在右、从左往右读。上文 RTR4 的「先 $\hat{\mathbf{q}}$ 再 $\hat{\mathbf{r}}$ 得 $\hat{\mathbf{r}}\hat{\mathbf{q}}$」与此一致；差别在于 RTR4 的矩阵是列向量，所以矩阵与四元数顺序相同，而 GEA 的行向量矩阵与四元数顺序相反。

四元数与矩阵互转（卷 I 5.4.4）

设 $\mathbf{q}=[x\ y\ z\ w]$，行向量约定下 $\mathbf{R}=\begin{bmatrix}1-2y^2-2z^2&2xy+2zw&2xz-2yw\\2xy-2zw&1-2x^2-2z^2&2yz+2xw\\2xz+2yw&2yz-2xw&1-2x^2-2y^2\end{bmatrix}$ 它是 RTR4 列向量形式（Eq. 4.46）的转置。矩阵到四元数，GEA 给出改编自 Nick Bobic 1998 年 Gamasutra 文章的 C 代码（q[0..2] 为向量部分，q[3] 为标量），思路与 RTR4 相同：

• 迹 $>0$：$s=\sqrt{\text{tr}+1}$，$q_w=s/2$，$t=0.5/s$，$q_x=(R_{21}-R_{12})t$，$q_y=(R_{02}-R_{20})t$，$q_z=(R_{10}-R_{01})t$；
• 否则找出对角线最大元素的下标 $i$，按 NEXT = {1, 2, 0} 得到 $j,k$，$s=\sqrt{R_{ii}-(R_{jj}+R_{kk})+1}$，$q_i=s/2$，$t=0.5/s$（$s=0$ 时 $t=s$），$q_w=(R_{kj}-R_{jk})t$，$q_j=(R_{ji}+R_{ij})t$，$q_k=(R_{ki}+R_{ik})t$。 GEA 另指出存在利用矩阵性质假设的更快转换方法。

原书问题：上面这段代码与 GEA 自己给出的行向量矩阵 $\mathbf{R}$ 符号不一致。代码沿用列向量矩阵的下标（$q_x=(R_{21}-R_{12})t$ 等），但对上式的行向量 $\mathbf{R}$ 有 $R_{21}-R_{12}=(2yz-2xw)-(2yz+2xw)=-4xw$，于是得到 $q_x=-x$；同理 $q_y,q_z$ 以及 else 分支里的 $q_w=(R_{kj}-R_{jk})t$ 也都变号，结果是**共轭（逆旋转）**而不是 $\mathbf{q}$。（本 wiki 用数值验证过：取任意单位 $\mathbf{q}$，按上式构造 $\mathbf{R}$，确认 $\mathbf{v}\mathbf{R}=\mathbf{q}\mathbf{v}\mathbf{q}^*$，再代入迹 $>0$ 分支，得到 $[-x\ -y\ -z\ w]$。）要配合 GEA 的行向量 $\mathbf{R}$，应交换这几处差值的两项（如 $q_x=(R_{12}-R_{21})t$），或者把矩阵转置后再传入。上文 RTR4 的公式用于列向量矩阵，本身没有问题；问题在于 GEA 把它原样用在行向量矩阵上。

记法分歧：[x y z w] 还是 [w x y z]

GEA 写作 $[x\ y\ z\ w]$，与齐次向量 $[x\ y\ z\ 1]$ 的习惯一致；很多学术论文按复数类比写成 $q=w+ix+jy+kz$，即 $[w\ x\ y\ z]$。读论文前先确认约定。（本页上文 RTR4 写作 $(\mathbf{q}_v,q_w)$，与 GEA 顺序相同。）

lerp 与 slerp（卷 I 5.4.5）

• 四元数 lerp：把四元数当四维向量线性插值再重新归一化：$\mathbf{q}_{lerp}=\dfrac{(1-\beta)\mathbf{q}_A+\beta\mathbf{q}_B}{|(1-\beta)\mathbf{q}_A+\beta\mathbf{q}_B|}$。最简单、最便宜。
• lerp 沿四维超球面的弦插值，$\beta$ 匀速变化时角速度不恒定：两端慢、中间快。slerp 沿大圆弧插值，角速度恒定：$\text{slerp}(\mathbf{p},\mathbf{q},\beta)=\frac{\sin(1-\beta)\theta}{\sin\theta}\mathbf{p}+\frac{\sin\beta\theta}{\sin\theta}\mathbf{q}$，$\cos\theta=\mathbf{p}\cdot\mathbf{q}$（四维点积）。
• 补充（GEA 未提及的实现注意事项）：因为 $\mathbf{q}$ 与 $-\mathbf{q}$ 表示同一旋转（见上文 RTR4），若 $\mathbf{p}\cdot\mathbf{q}<0$，应先把其中一个取负再做 lerp/slerp，否则会走「远路」，转过大于 180° 的角度；当 $\theta\to0$ 时 $\sin\theta\to0$，slerp 的权重会变成除以零，实践中此时退回归一化 lerp。
• 用不用 slerp 仍无定论：Jonathan Blow 撰文认为 slerp 太贵而 lerp 质量不差，主张理解 slerp 但不在引擎中使用；Naughty Dog 的同事发现好的 slerp 实现几乎和 lerp 一样快——PS3 SPU 上 Ice 团队的 slerp 每关节 20 个周期，lerp 每关节 16.25 个周期。作者建议先分别做性能分析：slerp 开销可接受就用（动画可能略好看），否则 lerp 对多数用途足够。
• 与上文 RTR4 对照：RTR4 也指出直接计算 slerp 要调三角函数、较贵，并给出在像素着色器里只做归一化代替 slerp 的做法（90° 时误差最大 4°）。两书结论相容：归一化 lerp 是廉价近似，是否值得用 slerp 取决于实测。RTR4 还提醒 slerp 不太适合相机朝向插值，GEA 未讨论这一点。

旋转表示的整体比较（欧拉角、矩阵、轴角、四元数、SRT、对偶四元数）见 旋转表示的比较。