四元数插值 Interpolating Quaternions

四元数插值 Interpolating Quaternions

对四元数进行 插值，本质上是在多个 旋转姿态 做 平滑过渡 或 权重混合，四元数插值有两种常见做法

这里所说的 “扭矩最小”（torque-minimal 或 minimum-torque） 是从 “刚体动力学” 角度来评价一条旋转插值曲线的物理性质：在理想均匀刚体、忽略外力的前提下，这条插值轨迹所需的外加扭矩 $\tau$（或其二范数的时间积分）是所有可行路径里最小的

Nlerp 归一化线性插值

对任何 标量 或 向量 做 lerp 线性插值 都是

四元数也不例外，只是 $\mathbf{a}, \ \mathbf{b}$ 必须保持单位长度，而进行 lerp 之后，长度一般不再等于 $1$，所以我们可以再最后进行一次 归一化 Normalize，这就是 Nlerp 归一化线性插值

可能不是最短路径

因为 单位四元数 和它的 相反数 $-q$ 在 三维空间 中表示 同一姿态，但却位于 四维超球面 的 对立半球，当我们直接执行 lerp 时，如果起点和终点分属两侧，会沿半个大圆以外的 “长弧” 进行插值，角度将多转 $180^\circ$，为了避免这种情况，我们可以用 点乘 判断是否同属一个半球，如果

那么说明它们分属 四维超球面 的两侧，我们需要将终点进行 取反

典型代码

Quaternion Nlerp(Quaternion start, Quaternion end, float t)
{
    if(Dot(start,end) < 0.0f)         // ① 取最短弧
        end = Negate(end);
    Quaternion res = start + (end - start) * t;   // ② 线性插值
    return Normalize(res);            // ③ 单位化
}

Nlerp 非常快速，只有加减乘＋一次 normalize

带权混合 Mixing

动画系统常用另一种 等价形式 来表示 lerp 和 Nlerp，即用 权重混合 来把两个（或多个）Pose/骨骼姿态叠加，对于四元数，也就是按 权重 对姿态求平均

• 去掉 $\mathrm{Normalize}$ 就是普通的 lerp 带权混合形式

当然，我们仍然需要用 点乘 判断是否同属一个半球用以取反来保证插值按 最短弧 进行

Quaternion Nlerp(Quaternion start, Quaternion end, float t)
{
    if(Dot(start,end) < 0.0f)         // ① 取最短弧
        end = Negate(end);
    Quaternion res =  a*(1-t) + b*t;   // ② 线性插值
    return Normalize(res);            // ③ 单位化
}

这在 动画层叠 / additive blending 里特别常用：你可以给若干动作片段各分配 0.3、0.5、0.2…… 的权重，把它们 “加权平均” 后再归一化，得到最终骨骼姿态

Slerp 球面线性插值

我们先来看看 三维向量 的 Slerp，其完全基于 大圆（测地线） 上做等速运动 的 几何约束

首先，若有 夹角 为 $\theta$ 的 $\vec{q}_0, \ \vec{q}_1$ 均是 三维单位球体 上的点，那么它们生成的 二维子空间 $P = \mathrm{span}\{\vec{q}_0, \ \vec{q}_1 \}$ 与 原点 相交于一个 单位圆，整条 最短弧 都落在这个平面里，也就是说我们将问题降维到一个 二维单位圆

我们首先需要构造 正交基，因为 $\vec{q}_0, \ \vec{q}_1$ 都是 三维单位球体 上的点，所以它们都是 单位向量，那么

而要构造 $\mathbf{e}_2$，我们先使用 Gram-Schmidt 正交化 $\vec{q}_1 - (\vec{q}_1 \cdot \vec{q}_0)\vec{q}_0$，然后再对其执行归一化，也就是除以它的长度 $\| \vec{q}_1 - (\vec{q}_1 \cdot \vec{q}_0)\vec{q}_0 \|$

我们知道 点乘 可以写为

• $\vec{q}_0, \ \vec{q}_1$ 都是 单位向量

所以 分子 可以写为

而 分母 也可以进行修改

所以我们得到最终的 正交基

$\mathbf{e}_1,\ \mathbf{e}_2$ 在平面 $P$ 内 正交归一，因此任何 插值点 都可写成

从 $\vec{p}_0 \to \vec{p}_1$ 的 插值点 必须满足以下 约束条件

1. 端点条件

   $$\begin{align*} \vec{q}(0) &= \vec{q}_0 = \mathbf{e}_1 &\implies x(0) = 1, \ y(0) = 0 \quad \quad \quad \\ \vec{q}(1) &= \vec{q}_1 = \cos{\theta}\mathbf{e}_1 + \sin{\theta}\mathbf{e}_2 &\implies x(1) = \cos{\theta}, \ y(1) = \sin{\theta} \end{align*}$$
   • $\vec{q}_1 = \cos{\theta}\mathbf{e}_1 + \sin{\theta}\mathbf{e}_2$ 是因为 $$\begin{align*} \vec{q}_1 \cdot \mathbf{e}_1 &= \vec{q}_1 \cdot \vec{q}_0 = \cos{\theta} \\ \vec{q}_1 \cdot \mathbf{e}_2 &= \vec{q}_1 \cdot \frac{\vec{q}_1 - \cos{\theta} \ \vec{q}_0}{\sin{\theta}} \\ &= \frac{\vec{q}_1\cdot\vec{q}_1 - (\vec{q}_1\cdot \vec{q}_0)(\vec{q}_1\cdot \vec{q}_0)}{\sin{\theta}} \\ &= \frac{1 - \cos^2{\theta}}{\sin{\theta}} \\ &= \frac{\sin^2{\theta}}{\sin{\theta}} \\ &= \sin{\theta} \end{align*}$$

2. 保持单位长度
   插值向量 的长度需要保持为 $1$

   $$\|\vec{q}(t)\|^2 = x(t)^2 + y(t)^2 = 1$$

3. 角度线性（匀速）
   我们知道 单位向量 $u$ 和 $v$ 的 点乘 与 夹角 关系是

   $$\vec{u} \cdot \vec{v} = \cos{[\angle(\vec{u}, \ \vec{v})]}$$

   我们还知道 $\vec{q}_0$ 和 $\vec{q}_1$ 的 完整夹角 是 $\theta$，而 插值向量 $\vec{q}(t)$ 与 初始向量 $\vec{q}_0$ 的 夹角 是 $\angle(\vec{q}_0,\ \vec{q}(t))$，我们希望它随着插值的进行，一路 线性 的增长到 $\theta$，也就是

   $$\angle(\vec{q}_0,\ \vec{q}(t)) = t\theta,\quad t\in [0,1]$$

   有了这个结论，再结合 单位向量 的 点乘 与 夹角 关系，可得

   $$\vec{q}_0 \cdot \vec{q}(t) = \cos{t\theta}$$

   我们之前还建立了 正交基 $\mathbf{e}_1,\mathbf{e}_2$，并用它表示了 插值向量 $\vec{q}(t)$，我们知道

   $$\begin{align*} \mathbf{e}_1 &= \vec{q}_0 \\ \vec{q}(t) &= x(t)\mathbf{e}_1 + y(t)\mathbf{e}_2 \end{align*}$$

   所以有

   $$\begin{align*} \vec{q}_0 \cdot \vec{q}(t) &= \mathbf{e}_1 \cdot [x(t)\mathbf{e}_1 + y(t)\mathbf{e}_2] \\ &= x(t)\mathbf{e}_1\cdot\mathbf{e}_1 + y(t)\mathbf{e}_2\cdot\mathbf{e}_1 \\ &= x(t)\cdot 1 + y(t)\cdot 0 \\ &= x(t) \end{align*}$$

   回代，最终得到

   $$x(t) = \cos{t\theta}$$

由 角度线性 约束，我们已经得到了 $x(t)$，接下来我们只需要结合 保持单位长度 约束就可以得到 $y(t)$

沿 最短弧 时 $y(t) \geq 0$ (与 $\mathbf{e}_2$ 同向)，故取正号，得到

所以我们得出

然后我们将 $\mathbf{e}_1, \ \mathbf{e}_2$ 展开写回为 $\vec{q}_0, \ \vec{q}_1$

我们可以利用 三角恒等式

所以我们有

回代，得到

这正是 正弦权重 版的 Slerp 公式，整条 插值曲线 因此可以视为在以 $\vec{q}_0, \ \vec{q}_1$ “按弦长正弦比” 进行的匀速运动，当 $\theta\to0$ 时，$\sin{\theta} \approx \theta$，此时公式滑退化为普通 lerp 线性插值 $\vec{q}_0 + (\vec{q}_1 - \vec{q}_0)t$

四元数 Slerp

要推导 四元数 $q_0$ 到 $q_1$ 的 Slerp，核心思想是在 四维单位球 $S^3$ 上寻找经过 $q_0 \to q_1$ 的最短大圆弧，并要求匀速（角速度恒定）地走这条弧线

单位四元数 对应 三维旋转，其集合构成 三维流形 $S^3$ （半径 $1$ 的 四维超球面），两个 四元数 点 $q_0,\ q_1 \in S^3$ 间的测地线就是穿过两点的 大圆（Great Circle），设

是两点在 $S^3$ 上的 夹角 (用四元数点积即可求得)，在大圆坐标系里，想要 匀速 行走意味着 插值点 必须满足

即 插值角度 $\theta(t)$ 与参数 $t$ 成 线性 关系

根据我们前面的推导，在二维圆上我们有以下 slerp 插值公式

我们将 $\vec{p}_0$ 和 $\vec{p}_1$ 换成 四元数 $q_0,\ q_1$ 就得到

四元数 Slerp - 指数形式

我们知道，四元数的相乘表示旋转的叠加，并且我们已知

• 起点姿态 $q_0$
• 终点姿态 $q_1$

我们需要找到一个 “中间旋转” $\Delta$ 使得

可得

可以认为我们先用 $q_0^{-1}$ 把当前姿态 “复原” 到基准坐标，再用 $q_1$ 把它转到目标姿态，整个组合就等价于一次旋转 $\Delta$

• 注意，$\Delta$ 是 “姿态到姿态” 的 差旋转 输入/输出都是完整的 四元数，包含 朝向 + 自旋，而我们前面章节介绍的，输入 from 向量 和 to 向量 来构造 from-to 四元数 是 “方向到方向” 的最短旋转，输入仅仅是两个 3D 向量，输出一个把第一个向量转到第二支向量的 四元数，不含轴向自旋
  • 旋转差 3 DoF 自由度：$\Delta = q_0^{-1}q_1$ 完整 SO(3) 旋转；既包含 “把坐标轴指向哪儿”，也包含 “绕该方向自旋了多少”
  • from-to 2 DoF 自由度：只能把一根方向矢量对齐到另一根；无法决定 绕该方向自旋（roll）

由于 $\Delta$ 也是表示旋转的 单位四元数，所以我们可以将它写成 轴-角 表示以及等价的 指数表示

如此一来，我们执行插值将会变的容易，我们可以取 对数 以后对角度按 $t$ 进行 线性插值

如果我们再对左右取 指数，可得

也就是说，要对 $\Delta$ 进行 $t \in [0,1]$ 且保持旋转 轴 不变，仅对 角度 进行 线性插值，我们只需要对原四元数 $\Delta$ 作 实数幂 操作 $\Delta^t$，所以 四元数幂 的实现是

Quaternion Pow(Quaternion q, float t)
{
    //---------------------------------------
    // 0)  若 q 不是严格单位四元数，先归一化
    //---------------------------------------
    q.normalize();   // 可省略，但最好保守处理
 
    //---------------------------------------
    // 1)  提取半角 θ/2
    //---------------------------------------
    float cosHalf = q.w;              // w = cos(θ/2)
    // 若 |cosHalf| ≈ 1 → 角度极小，需要特殊处理
    if (fabs(cosHalf) > 0.9999f)
    {
        // 与 Lerp 近似：q^t ≈ (1-t) + t*q
        return Quaternion::Lerp(Quaternion::Identity(), q, t).normalized();
    }
    float halfAngle = acos(cosHalf);  // θ/2  ∈ (0,π)
 
    //---------------------------------------
    // 2)  求旋转轴 n̂
    //---------------------------------------
    float sinHalf = sqrtf(1.0f - cosHalf * cosHalf);   // = |v|, > 0 因为上面排除了极小角
    Vector3 axis(q.x / sinHalf, q.y / sinHalf, q.z / sinHalf);  // 已单位化
 
    //---------------------------------------
    // 3)  把角度乘 t，再算新的 sin/cos
    //---------------------------------------
    float newHalf = t * halfAngle;     // = tθ/2
    float newCos = cosf(newHalf);
    float newSin = sinf(newHalf);
 
    //---------------------------------------
    // 4)  组合回四元数
    //---------------------------------------
    return Quaternion(axis * newSin, newCos);   // (x,y,z,w)
}

• 若平台有 sincosf() 或 SIMD sin_cos(), 可以一次性同时取 sin 和 cos，更省性能

最终，我们可以得到四元数 Slerp 的 指数幂形式

它与 三角函数 写法完全等价，只需将 $(q_0^{-1}q_1)^t$ 展开成 $(\cos{\frac{\Omega}{2}},\ \hat{n}\sin{\frac{\Omega}{2}})$ 并用 倍角公式 即可，但是 指数 写法更紧凑，且与实现 Pow() 函数（幂运算）天然耦合

struct Quaternion { float x, y, z, w; /* …ctor & ops… */ };
 
Quaternion Slerp(const Quaternion& a_, const Quaternion& b_, float t)
{
    Quaternion a = a_, b = b_;
 
    // ① 最短弧：同半球
    if (Dot(a, b) < 0.0f)
        b = -b;
 
    // ② 差旋转
    Quaternion delta = b * Inverse(a);
 
    // ④ 小角退化（利用 delta.w ≈ cos(θ/2)）
    if (fabs(delta.w) > 0.9999f)
        return Normalize(Lerp(a, b, t));   // Nlerp 近似 
 
    // ③ Δ^t —— 直接用四元数 Pow() 函数
    Quaternion deltaPow = Pow(delta, t);
 
    // ⑤ 左乘起点
    return deltaPow * a;
}
 

以上是 指数形式 实现的 四元数 Slerp，但是不管是 指数形式 还是 三角函数形式，都无法避免 acosf，所以四元数的 Slerp 存在一定的性能消耗

四元数 Slerp vs 向量 Slerp

四元数保存了 “绕方向自转” 的信息；为了保持 匀速自旋，就得处理 $q_0^{-1}q_1$ 再做幂运算，而向量只关心方向角度，公式更简单，如果只旋转镜头指向，而不想处理 roll，向量 Slerp 就够；要连同 roll 一起平滑过渡，必须用 四元数 Slerp

双覆盖问题

在上面的分析中我们提到，在两个四元数之间进行插值时，插值可以采用两者之间最短或最长的圆弧，因为这两个四元数都代表 同一姿态

我们一般会想让插值沿着 最短路径 进行，如果你曾经见过游戏中角色的弯曲错误，但动作却正确，那么很可能是没有沿着 最短路径 执行插值，我们在上面的代码中已经加入了 半球检测

// ① 最短弧：同半球
if (Dot(a, b) < 0.0f)
    b = -b;