3 - 使用对偶四元数插值

在看完GAMES104动画系统那一部分后，我也想写一个简单的动画系统，又要开一个坑了！本文将简要介绍如何用对偶四元数给蒙皮动画插值。

对偶四元数

背景

我们在之前用的是直接四元数蒙皮（Direct Quaternion Blending）方法，它将变换拆成两部分，平移和缩放用矩阵表示，旋转用四元数表示，可能会导致一些严重的断裂现象；而使用对偶四元数插值平移和旋转，缩放用矩阵表示则会使插值结果流畅。

数学原理

一个对偶四元数$dq$由两个四元数$p,q$和一个$\epsilon$组成：

$$\nonumber dq=p+\epsilon q\,,\text{其中}\ \epsilon^2=0\ \text{且} \ \epsilon \neq0$$

其中$p$是实部(real)，$q$是对偶部(dual)。

然后对于蒙皮动画，只需知道对偶四元数的加法即可，有对偶四元数$dq_1$和$dq_2$如下：

$$\begin{align} \nonumber dq_1=p_1+\epsilon q_1 \\ \nonumber dq_2=p_2+\epsilon q_2 \end{align}$$

那么它俩的和如下：

$$\nonumber dq_1+dq_2=(p_1p_2)+\epsilon(q_1+q_2)$$

存储数据

对于两个四元数，有如下性质：

• 进行加法后归一化，就能得到两个四元数的均值。很适合顶点的旋转。
• 进行加法后不归一化，相当于两个四维向量相加。很适合顶点的平移。

因此可以将旋转和平移存储为四元数，用一个对偶四元数来管理。首先是作为实部的$p$，用于存储旋转：

$$\nonumber p(\phi)=\cos(\frac{\phi}{2})+\sin(\frac{\phi}{2})i+\sin(\frac{\phi}{2})j+\sin(\frac{\phi}{2})k$$

然后将平移存储到对偶四元数的虚部$q$中：

$$\nonumber q(t)=(\frac{t_x}{2}i+\frac{t_y}{2}j+\frac{t_z}{2}k)*p$$

那么平移信息可由如下计算得到：

$$\nonumber t=2*q*p$$

需要注意的是，在最后别忘了归一化对偶四元数，否则会出现一些模型上的变形。

此外，对偶四元数可以用矩阵表示：

$$\nonumber \epsilon= \begin{bmatrix} 0 & 1\\ 0 & 1 \end{bmatrix} \quad\text{及}\quad a+b\epsilon= \begin{bmatrix} a & 0\\ b & a \end{bmatrix}$$

代码实现

glm数学库

glm在<glm/gtx/dual_quaternion.hpp>中实现了对偶四元数：

glm::dualquat dq;

可用数组索引的方式来获取对偶四元数的实部和虚部：

glm::quat p = dq[0];
glm::quat q = dq[1];

由于GLSL不支持四元数和对偶四元数，需要在CPU侧将四元数转换成2x4矩阵，然后传给GPU。glm通过glm::mat2x4_cast()将一个对偶四元数转换为一个2x4矩阵：

glm::mat2x4 dqMat = glm::mat2x4_cast(dq);

使用对偶四元数

在Animator类中，添加一个新的成员变量，用于存储关节的对偶四元数：

std::vector<glm::mat2x4> m_boneDualQuaternions;

这里使用glm::mat2x4的原因是Shader只接受这种类型的四元数。

然后别忘了在构造函数中初始化它的大小：

m_boneDualQuaternions.resize(100);

由于使用uniform传递变量，数组元素最多100个，对于一些动画关键帧很多的模型无法使用。在下一篇文章中我们将使用 着色器存储缓冲对象SSBO 来解除这一限制，此外也能预先将关键帧信息存入纹理，然后采样。

接下来需要通过glm::decompose()将计算好的蒙皮矩阵拆解，结果会存储到下列临时变量中，其中我们只用到旋转和平移：

// 使用对偶四元数
//    1. 准备临时变量
glm::quat orientation;
glm::vec3 scale;
glm::vec3 translation;
glm::vec3 skew;
glm::vec4 perspective;
//    2. 拆解蒙皮矩阵
if (glm::decompose(m_finalBoneMatrices[index], scale, orientation, translation, skew, perspective))
{
    glm::dualquat dq;
    dq[0] = orientation;
    dq[1] = glm::quat(0.0, translation.x, translation.y, translation.z) * orientation * 0.5f;
    m_boneDualQuaternions[index] = glm::mat2x4_cast(dq);
}
else
{
    LOG_WARN(std::format("[{}] Could not decompose skinning matrix for bone {}, use direct quat instead for animation...", __FUNCTION__, index));
    m_useDualQuaternion = false;
}

修改着色器

CPU端有数据以后，就该修改着色器代码了，首先添加对应的uniform变量：

uniform mat2x4 u_FinalBonesDQs[MAX_BONES];

然后写一个新的函数getBoneTransform()获取加权且插值后的对偶四元数：

// 获取加权且插值后的对偶四元数
mat2x4 getBoneTransform()
{
    // 获取影响该顶点的对偶四元数
    mat2x4 dq0 = u_FinalBonesDQs[boneIds.x];
    mat2x4 dq1 = u_FinalBonesDQs[boneIds.y];
    mat2x4 dq2 = u_FinalBonesDQs[boneIds.z];
    mat2x4 dq3 = u_FinalBonesDQs[boneIds.w];

    // 根据最短旋转路径加权, 防止动作发生突变
    vec4 shortestWeights = weights;
    shortestWeights.y *= sign(dot(dq0[0], dq1[0]));
    shortestWeights.z *= sign(dot(dq0[0], dq2[0]));
    shortestWeights.w *= sign(dot(dq0[0], dq3[0]));

    // 得到加权插值后的对偶四元数, 别忘了标准化
    mat2x4 result = shortestWeights.x * dq0 +
                    shortestWeights.y * dq1 +
                    shortestWeights.z * dq2 +
                    shortestWeights.w * dq3;
    float norm = length(result[0]);
    return result / norm;
}

有了对偶四元数后，还需要将其转换为矩阵，通过新的函数getSkinMatFromDQ()实现：

mat4 getSkinMatFromDQ()
{
    mat2x4 boneDQ = getBoneTransform();

    vec4 r = boneDQ[0];        // 旋转
    vec4 t = boneDQ[1];        // 平移

    return mat4(
        1.0 - (2.0 * r.y * r.y) - (2.0 * r.z * r.z),
              (2.0 * r.x * r.y) + (2.0 * r.w * r.z),
              (2.0 * r.x * r.z) - (2.0 * r.w * r.y),
        0.0,

              (2.0 * r.x * r.y) - (2.0 * r.w * r.z),
        1.0 - (2.0 * r.x * r.x) - (2.0 * r.z * r.z),
              (2.0 * r.y * r.z) + (2.0 * r.w * r.x),
        0.0,

              (2.0 * r.x * r.z) + (2.0 * r.w * r.y),
              (2.0 * r.y * r.z) - (2.0 * r.w * r.x),
        1.0 - (2.0 * r.x * r.x) - (2.0 * r.y * r.y),
        0.0,

        2.0 * (-t.w * r.x + t.x * r.w - t.y * r.z + t.z * r.y),
        2.0 * (-t.w * r.y + t.x * r.z + t.y * r.w - t.z * r.x),
        2.0 * (-t.w * r.z - t.x * r.y + t.y * r.x + t.z * r.w),
        1
    );
}

其中，矩阵转换的格式如下：

$$\nonumber T= \begin{bmatrix} & & & t_x \\ & R & & t_y \\ & & & t_z \\ 0 & 0 & 0 & 1 \end{bmatrix}$$

R是旋转矩阵，由四元数转换而成。

那么就能通过如下形式变换顶点了：

mat4 skinMat = getSkinMatFromDQ();

vs_objData.Normal = u_Normal * normalize(mat3(skinMat) * normal);
vs_objData.FragPos = vec3(u_Model * skinMat * vec4(position, 1.0));
vs_objData.TexCoords = texCoords;

gl_Position = u_MVP * skinMat * vec4(position, 1.0);

这样应该就能用对偶四元数进行蒙皮动画插值了，效果应该会更好。在下一篇文章中我们将打破uniform数组的限制，使用着色器存储缓存对象SSBO来存储更多的关键帧信息！

参考资料

• 蒙皮算法之对偶四元数 (Skinning algorithm for dual Quaternions) - 知乎 (zhihu.com)
• 《C++ Game Animation Programming》