图形学基础（3）——模型变换动画

模型变换动画

　　当场景中的物体进行运动时，有时候不能每一帧都由人来控制，可以使用相近的动作为两个关键帧，然后进行插值。每帧动画其实就是模型特定姿态的一个“快照”，通过在帧之间插值的方法，从而得到平滑的动画效果。
　　一种比较好的运动动画是矩阵分解（matrix decomposition）——给予任意变换矩阵 M，可以将它分解为缩放矩阵 S，旋转矩阵 R 和平移矩阵 T，即

M=SRT

　　然后每一个矩阵进行独立插值，然后将插值完成的矩阵再次用乘法组合起来使用。
　　平移和缩放矩阵的插值非常简单，可以使用线性插值达到非常精确的效果，对渲染的插值就要困难很多。为了得到平滑的插值，我们需要使用四元数，关于四元数的定义不再赘述，只介绍稍后要用到的四元数插值。
　　
　　四元数其实是四维超球面上的点，直接使用线性插值 LERP 运算实际上是沿超球的弦上进行插值，而不是在超球面上插值，这样会导致旋转动画并非以恒定角速度进行。旋转在两端看似较慢，但在动画中间就会较快。
　　可以使用 LERP 的变体——球面线性插值（spherical linear interpolation，SLERP）解决这个问题。SLERP 使用正弦和余弦在四维超球面的大圆上进行插值，即

Slerp(p,q,β)=ωpp+ωqq

　　其中：

wp=sin((1−β)θ)sinθ

wq=sin(βθ)sinθ

　　两个单位四元数之间的夹角，可以使用四维点积求得 cosθ，再求反余弦：

cosθ=p⋅q=pxqx+pyqy+pzqz+pwqw

θ=cos−1(p⋅q)

　　有了四元数的基础，可以尝试实现 AnimatedTransform 类了，以下代码来自 pbrt 的实现，矩阵变换分解为缩放，旋转和平移的函数是 Decompose()。

　　<AnimatedTransform 类构造函数> ≡    AnimatedTransform::AnimatedTransform(const Transform *startTransform,                   startTime, const Transform *endTransform, Float endTime)             : startTransform(startTransform), endTransform(endTransform),                   startTime(startTime), endTime(endTime),                   actuallyAnimated(*startTransform != *endTransform) {        Decompose(startTransform->m, &T[0], &R[0], &S[0]);        Decompose(endTransform->m, &T[1], &R[1], &S[1]);        <如果需要选择最短路径，需要将 R 取相同的符号>        if (Dot(R[0], R[1]) < 0) R[1] = -R[1];        hasRotation = Dot(R[0], R[1]) < 0.9995f;        <然后计算需要的运动微分函数项>    }    <私有成员>    const Transform *startTransform, *endTransform;    const Float startTime, endTime;    const bool actuallyAnimated;    Vector3f T[2];    Quaternion R[2];    Matrix4x4 S[2];    bool hasRotation;

　　如果得到一个合成好的变换矩阵，其实合成它的单个变换矩阵的细节已经消失了，同样的矩阵可以使用不同数量的分解矩阵来合成，所以需要规范分解的顺序，一种推荐的分解形式是：

M=TRS

　　其中，M是给定的变换，T 是平移，R 是旋转，S 是缩放，S是广义上的缩放，并不是当前坐标的缩放，但无论如何，都可以实现精确的线性插值。

    <Decompose 函数定义>    void AnimatedTransform::Decompose(const Matrix4x4 &m, Vector3f *T,        Quaternion *Rquat, Matrix4x4 *S) {        <首先将平移 T 提取出来>        <然后计算没有平移分量的矩阵 M>        <然后从变换矩阵中提取旋转分量 R>        <最后用旋转分量和最初的变换矩阵计算缩放分量 S >    }

　　提取平移变换 T 只需要将第 4 列取出即可。

    <首先将平移 T 提取出来>    T->x = m.m[0][3];    T->y = m.m[1][3];    T->z = m.m[2][3];

　　然后将第四列用(0,0,0,1) 代替即可，此时左上的 3x3矩阵就是旋转缩放混合矩阵。比较有挑战的是提取旋转分量，这里采用极分解（polar decomposition），重复给 M 和 M 的逆转置取平均来得到旋转分量 R 和缩放分量S，直到收敛 Mi=R。

Mi+1=12(Mi+(MTi)−1)

　　如果 M 是纯旋转矩阵，那么 (MTi)−1 与 M 相同，那么就可以立即退出运算。极分解是已经被证明的定理，详情请见Polar decomposition。迭代运算最终会收敛到一个特别小的范围或者定值，实践证明，这个运算收敛得也很快，当计算迭代 100 次或者两个矩阵之间的插值小于 0.0001 也就可以停止循环。

    <然后从变换矩阵中提取旋转分量 R>    Float norm;    int count = 0;    Matrix4x4 R = M;    do {        <计算下一个矩阵 Rnext>        <计算两个矩阵之间的相差的定值>        R = Rnext;    } while (++count < 100 && norm > .0001);    *Rquat = Quaternion(R);

    <计算下一个矩阵 Rnext>    Matrix4x4 Rnext;    Matrix4x4 Rit = Inverse(Transpose(R));    for (int i = 0; i < 4; ++i)    for (int j = 0; j < 4; ++j)    Rnext.m[i][j] = 0.5f * (R.m[i][j] + Rit.m[i][j]);

    <计算两个矩阵之间的相差的定值>    norm = 0;    for (int i = 0; i < 3; ++i) {        Float n = std::abs(R.m[i][0] - Rnext.m[i][0]) +        std::abs(R.m[i][1] - Rnext.m[i][1]) +        std::abs(R.m[i][2] - Rnext.m[i][2]);        norm = std::max(norm, n);    }

　　提取出旋转矩阵之后，就需要找到满足 M=RS 的缩放分量 S，所以有 S=R−1M。

    <计算缩放矩阵>    *S = Matrix4x4::Mul(Inverse(R), M);

　　对于旋转矩阵的四元数，正负表示同一个旋转，如果点乘两个旋转后得到的值是负的，那么进行球面插值 Slerp 获得的不是最短路径（请理解四维超球面），所以要有一个取负值参与运算。
　　最后是计算运动微分方程，代码比较复杂，请见 pbrt 源码。

　　AnimatedTransform 中另一个必须要有的函数是插值函数 Interpolate()，使用输入的时间得到输出的变换矩阵。

    void AnimatedTransform::Interpolate(Float time, Transform *t) const {        <获得变换的边界条件>        Float dt = (time - startTime) / (endTime - startTime);        <在 dt 点插值平移>        <在 dt 点插值旋转>        <在 dt 点插值缩放>        <合成变换矩阵并返回>    }

　　如果给予的时间值超出范围则立即返回，如果构造 AnimatedTransform 类时起始变换和终止变换相同，则 actuallyAnimated 为 true，也就没有插值的必要了。

    <获得变换的边界条件>    if (!actuallyAnimated || time <= startTime) {        *t = *startTransform;        return;    }    if (time >= endTime) {        *t = *endTransform;        return;    }

　　平移和缩放使用线性插值，旋转使用球面插值。

    <在 dt 点插值平移>    Vector3f trans = (1 - dt) * T[0] + dt * T[1];    <在 dt 点插值旋转>    Quaternion rotate = Slerp(dt, R[0], R[1]);    <在 dt 点插值缩放>    Matrix4x4 scale;    for (int i = 0; i < 3; ++i)        for (int j = 0; j < 3; ++j)            scale.m[i][j] = Lerp(dt, S[0].m[i][j], S[1].m[i][j]);

　　最后合成矩阵并返回。

    *t = Translate(trans) * rotate.ToTransform() * Transform(scale);

　　在物体进行旋转动画时，也需要保证包围盒的变换，始终将物体包含在包围盒内。计算旋转的困难在于运动的时候难以确定极值，许多渲染器的做法是对这段时间内对包围盒进行多次插值，计算每一个包围盒转换后的位置，再对所有的包围盒做 Union 操作。
　　pbrt 采用微分方程计算极值的方法，算法比较复杂，有时间会再进行专题解析，感兴趣的可以翻看 pbrt 源码。

参考：
　　Physically Based Rendering
　　游戏引擎架构