Learned Motion Matching

_-ARCHER-_

已于 2022-05-23 15:58:16 修改

阅读量2k

点赞数 6

分类专栏：神经网络文章标签：算法 python 动画

于 2022-05-16 11:10:02 首次发布

本文链接：https://blog.csdn.net/qq_41668705/article/details/124793351

版权

神经网络专栏收录该内容

1 篇文章 0 订阅

订阅专栏

Learned Motion Matching

文章链接: https://dl.acm.org/doi/pdf/10.1145/3386569.3392440
Daniel Holden 大佬主页：https://theorangeduck.com/page/publications
unity项目AI4Animation地址(17-19年论文的实现)：https://github.com/sebastianstarke/AI4Animation

一、Motion Matching

输入： 预期轨迹和其他特征的 feature

输出： 最匹配的 pose

在这里插入图片描述
具体来说一个传统的motion matching分为以下几步：

Projection: 当用户输入的轨迹和查找到的目标 Frame 的轨迹不一致（或是出路很大），则需要重新搜索，搜索的过程主要是通过用户模拟的轨迹和当前 pose 的状态来去 database 中进行 match，然后作为修正的 Frame
Decompression: 要把查找到的目标 Frame 重新映射回反正的目标 Pose
Stepping: index++ 逐步推出下一个的 Frame

*注：此时的 Frame 只是用户输入的信息与之匹配的 x，并没有做任何映射改变。

按这样的思路，每帧去从动画数据库中去查询，得到最优 Pose，就是motion matching的基本思路。

问题：

Projection过程中匹配计算开销大，并且database也非常大，此外根据查找的 Frame 重新解压缩Pose的过程也不一定能完美的进行还原
不会组合数据，如果有步行动画和饮酒动画，Motion Match 无法生成角色边喝酒边走路的场景，Motion Match仅播放提供的数据

基于这些缺点，就有了《Learned Motion Matching》这篇 paper。

二、Learned Motion Matching

输入： 预期轨迹和其他特征的 feature

输出： 最匹配的 pose

在这里插入图片描述

LMM具有和MM一样的输入和输出，值得注意的是大体思路没有变，以上所描述的MM具体步骤同样可以适用于LMM，不同之处在于：用网络对数据库进行压缩，来代替数据库，进而达到降低模型的大小

在 LMM 中将 x 叫做 feature，并对其增广，增加得到一个隐变量 latent z

网络结构

Decompressor D

MM：通过 Frame index 直接映射至pose
LMM：额外添加了隐变量 latent z 来作为新添加的特征，这个 latent 由将 Pose y 作为输入，通过一个auto encoder学习到的，作者将这个encoder过程叫做Compressor C ，这也是LMM最关键的一步（正是因为这一步，其他匹配全变成了网络）

Projection P

MM：Motion Matching Search，去查询 Frame
LMM：通过sample各种可能的用户输入，然后寻找最近的 feature 和 latent ，组成pair的数据，从而可以 train 对应的projector。通过用户输入，Projection 生成匹配的动作（在这个部分可以进行 “模态” 转换），比如跑步到跳这个动作，可以用一个变量来对 feature 进行增广，使其有0 -> 1的剧烈变化

stepper S

MM：根据 x 对应的motion库中的index，index++ 就能得到下一帧的x
LMM：因为把这个motion库给去掉，所以作者train了一个自回归的网络，能够根据当前的 feature 和 latent ,来推测下一时刻的 feature 和 latent

Compressor C

训练 auto encoder 可提供每帧额外特征的矢量 latent，这些矢量会自动选择以提高Decompressor 的准确性。所以必须将其与Maching Feature Database一起存储到它们自己的额外的database中

MM的匹配方式：

LMM的匹配方式：

可以看出将原来的Animation Dataset 换成了Matching Feature Database 和一个Extra Freature Database ，

Matching Feature Database：为了减小时间空间开销
Extra Freature Database：Matching Feature不能很好地携带足够的信息，使我们能够在所有情况下重建姿势，这些Extra Freature会自动选择以提高 Decompressor 的准确性，而在论文中将这一部分叫做 Latent

节省的内存：

2022.5.18 继续更新

三、TOOL

阻尼器

举个栗子：在打游戏的时候，我们先往前快跑，然后瞬间再往后走，人物按照惯性不可能是直接往后跑的，所以在程序设定了一个 阻尼器 来满足用户 Input 至网络 Input 的过渡。如果需要原理链接在这里：

Spring-It-On: The Game Developer’s Spring-Roll-Call

这里以一个 float 类型的阻尼器为例

static inline void simple_spring_damper_implicit(
    float& x, 				// 位置
    float& v, 				// 速度
    const float x_goal, 	// 目标位置
    const float halflife,	// 阻尼
    const float dt)			// 时间间隔 1/fps
{
    float y = halflife_to_damping(halflife) / 2.0f; 
    float j0 = x - x_goal;
    float j1 = v + j0*y;
    float eydt = fast_negexpf(y*dt);

    x = eydt*(j0 + j1*dt) + x_goal;
    v = eydt*(v - j1*y*dt);
}

在实现过程中不仅仅只有一个位置信息，可能这是一个角度信息，比如人物的朝向，都需要用到该功能，所以通过复用，比如说输入位置信息 x 为 vec3 向量，或者输入角度信息 x 为 quat 四元数，这几个函数作用都是一样的，都是加上阻尼，防止输入信息变化过快。

Quaternion

在这里四元数是工具，我只想要知道怎么用就够了，所以理论部分请参考大佬的文章：

四元数与三维旋转

首先，四元数的 标向量形式定义如下
$\ q = \left[ \begin{matrix} w, \bm s \end{matrix} \right],( \bm s = \left[ \begin{matrix} x \\ y \\ z \end{matrix} \right])$
四元数的乘法可以表示为

$\ q_aq_b = \left[ w_aw_b - \bm s_a \bm s_b , w_b\bm s_a +w_a\bm s_b+ \bm s_a \times \bm s_b \right]$

应用于旋转的四元数的 标向量形式可以表示为
$\ q = \left[ \begin{matrix} cos(\frac{\theta}{2}),sin(\frac{\theta}{2}) \bm u \end{matrix} \right],( \bm u = \left[ \begin{matrix} x \\ y \\ z \end{matrix} \right])$

通过一个四元数 $q$ 对一个三维点 $p=\left[ 0,\bm v \right]$ (把 $p$ 当成一个实部为 $0$ 的四元数)施加旋转变换，只需要做如下计算:
$p'=qpq^{-1} \\ qpq^{-1}=\left[ 0, \bm v + 2cos(\frac{\theta}{2})(\bm u\times \bm v) +2\bm u\times (\bm u\times \bm v)\right]$

输入： $p$ 为对应父节点的 Local 坐标
输出： $p^{'}$ 为以父节点为0点位置的的世界坐标

注： $p^{'}$ 也是一个实部为 $0$ 的四元数

四元数相乘的代码为：

static inline quat quat_mul(quat q, quat p){
  return quat(
    p.w*q.w - p.x*q.x - p.y*q.y - p.z*q.z,
    p.w*q.x + p.x*q.w - p.y*q.z + p.z*q.y,
    p.w*q.y + p.x*q.z + p.y*q.w - p.z*q.x,
    p.w*q.z - p.x*q.y + p.y*q.x + p.z*q.w);
}

四元数旋转的代码为：

static inline vec3 quat_mul_vec3(quat q, vec3 v){
    vec3 t = 2.0f * cross(vec3(q.x, q.y, q.z), v);
    return v + q.w * t + cross(vec3(q.x, q.y, q.z), t);
}

static inline vec3 cross(vec3 v, vec3 w){
    return vec3(
        v.y*w.z - v.z*w.y,
        v.z*w.x - v.x*w.z,
        v.x*w.y - v.y*w.x);
}

之后设计Local & Globel 坐标之间的转换全部都是用四元数，使用的时候只需要知道输入是啥输出是啥就好了

Forward Kinematics

定义： Forward Kinematics 是指使用机器人的运动学方程根据关节参数的指定值计算末端执行器的位置

简单来说就是已知每一个关节的位置角度信息，来计算每一个关节的位置。在 LMM 中和普通的 FK 不同之处在于旋转使用的是四元数，此时的 FK 可以递归定义如下：

$\{ q^{t} ,q^{r},\dot q^{t},\dot q^{r} \}$
where
$q_{j}^{t} =q_{p}^{t} + q_{p}^{r} \otimes y_{j}^{t} \\ q_{j}^{r} =q_{p}^{r} \otimes y_{j}^{r} \\ \dot q_{j}^{t} =\dot q_{p}^{t} + q_{p}^{r} \otimes \dot y_{j}^{t} + \dot q_{p}^{r} \times(q_{p}^{r} \otimes y_{j}^{t}) \\ \dot q_{j}^{r} =\dot q_{p}^{r} + q_{p}^{r} \otimes y_{j}^{r} \\$

基于以上公式，设root 节点的 index 为 0，对于所有关节的 KF 可以写为函数（只包含 positions 和 rotations ）：

void forward_kinematics_full(
    slice1d<vec3> global_bone_positions,		// 输出全局坐标
    slice1d<quat> global_bone_rotations,		// 输出全局旋转角度
    const slice1d<vec3> local_bone_positions,	// 输入局部坐标
    const slice1d<quat> local_bone_rotations,	// 输入局部旋转角度
    const slice1d<int> bone_parents)			// 父节点（ ROOT 的父节点为 -1）
{
    for (int i = 0; i < bone_parents.size; i++) {
        // Assumes bones are always sorted from root onwards
        assert(bone_parents(i) < i);
        
        if (bone_parents(i) == -1)  {
            global_bone_positions(i) = local_bone_positions(i);
            global_bone_rotations(i) = local_bone_rotations(i);
        }
        else  {
            vec3 parent_position = global_bone_positions(bone_parents(i));
            quat parent_rotation = global_bone_rotations(bone_parents(i));
            global_bone_positions(i) = quat_mul_vec3(parent_rotation, local_bone_positions(i)) + parent_position;
            global_bone_rotations(i) = quat_mul(parent_rotation, local_bone_rotations(i));
        }
    }
}

如果是确定某一个关节点的位置，那么可以使用递归来完成，这个程序中包括了position velocity rotation 和 angular_velocity 四个变量的传递

// Forward kinematics but also compute the velocities
void forward_kinematics_velocity(
    vec3& bone_position,							// 输出全局坐标
    vec3& bone_velocity,							// 输出全局速度
    quat& bone_rotation,							// 输出全局旋转角度
    vec3& bone_angular_velocity,					// 输出全局旋转角速度
    const slice1d<vec3> bone_positions,				// 输入局部坐标
    const slice1d<vec3> bone_velocities,			// 输入局部速度
    const slice1d<quat> bone_rotations,				// 输入局部旋转角度
    const slice1d<vec3> bone_angular_velocities,	// 输入局部旋转角速度
    const slice1d<int> bone_parents,				// 父节点索引集
    const int bone)									// 需要确定的关节点位置
{
    if (bone_parents(bone) != -1) {					// root 节点
        vec3 parent_position;
        vec3 parent_velocity;
        quat parent_rotation;
        vec3 parent_angular_velocity;
        
        forward_kinematics_velocity(...);			// 递归
        
        bone_position = quat_mul_vec3(parent_rotation, bone_positions(bone)) + parent_position;
        bone_velocity = 
            parent_velocity + 
            quat_mul_vec3(parent_rotation, bone_velocities(bone)) + 
            cross(parent_angular_velocity, quat_mul_vec3(parent_rotation, bone_positions(bone)));
        bone_rotation = quat_mul(parent_rotation, bone_rotations(bone));
        bone_angular_velocity = quat_mul_vec3(parent_rotation, bone_angular_velocities(bone) + parent_angular_velocity);
    }
    else{
        bone_position = bone_positions(bone);
        bone_velocity = bone_velocities(bone);
        bone_rotation = bone_rotations(bone);
        bone_angular_velocity = bone_angular_velocities(bone); 
    }
}

四、LMM Test

工欲善其事，必先利其器，按照以上对 LMM 的描述，不妨先来看看 LMM 是怎样对动作进行预测匹配的，其代码为 controller.cpp 但其中包括着许多游戏视角、游戏控制有关的代码，为了简洁描述我将不必要的代码删光了，只保留了必要的部分。

简单介绍

首先，LMM 的运行需要提供以下.bin 文件

database db;
database_load(db, "./resources/database.bin");   			// 包含了 X Y 的信息

nnet decompressor, stepper, projector;
nnet_load(decompressor, "./resources/decompressor.bin");	// D network 的参数
nnet_load(stepper, "./resources/stepper.bin");				// S network 的参数
nnet_load(projector, "./resources/projector.bin");			// P network 的参数

controller.cpp 中只实现了跑步走路，还有一些时间转动的 function，万事开头难，开头简单万事就简单，controller.cpp 中一些重要的参数能够方便理解整个程序

玩家控制的参数：

float desired_gait = 0.0f;			// 步态，0~1 代表着从走路到跑步的过度

vec3 gamepadstick_left = vec3();	// 控制移动方向，是一个二维（x,0,z）的单位向量

在controller.cpp 中还有控制视角的右摇杆变量，就不做考虑

网络中重要的参数：

float feature_weight_foot_position;			// 足部 位置 权重
float feature_weight_foot_velocity;			// 足部 速度 权重
float feature_weight_hip_velocity;			// 髋部 速度 权重
float feature_weight_trajectory_positions;	// 轨迹 位置 权重
float feature_weight_trajectory_directions;	// 轨迹 速度 权重

// 每个变量组成：关节数 * 3D向量
// local	（相较于父节点的坐标，ROOT节点除外）
array1d<vec3> curr_bone_positions; 			// 每个关节的局部位置 	size: 23 * 3
array1d<vec3> curr_bone_velocities;			// 每个关节的局部速度 	size: 23 * 3	
array1d<quat> curr_bone_rotations;			// 每个关节的局部四元数 	size: 23 * 4
array1d<vec3> curr_bone_angular_velocities;	// 用四元数的轴角表示速度 size: 23 * 3
array1d<bool> curr_bone_contacts;			// 两脚是否接地（锁足部分 size:  2 * 1	

// global	（世界坐标，ROOT节点还是原来的坐标）
array1d<vec3> global_bone_positions;		// 每个关节的全局位置 	size: 23 * 3
array1d<vec3> global_bone_velocities;		// 每个关节的全局速度 	size: 23 * 3
array1d<quat> global_bone_rotations;		// 每个关节的全局四元数 	size: 23 * 4
array1d<vec3> global_bone_angular_velocities;// 用四元数的轴角表示速度 size: 23 * 3

Evaluate Network

Projector

projector_evaluate() 的作用就是通过网络 Projector P 对原姿势进行更新。MM 或是 LMM 其原理都是从数据库匹配，只不过 LMM 是用网络在数据库中匹配。具体的步骤可以表示为：

标准化 feature 数据，生成 Projector 的输入
运行 Projector
计算更新之后的 feature 和更新前的 feature 之间的距离，差距过大则输出更新之后的 feature；差距不大则输出原 feature

void projector_evaluate(
    bool& transition,						// 是否需要转移（Projector 匹配结果和现阶段比较）
    float& best_cost,						// 判断未标准化feature之间差距
    					    				// trns_dist_squared 判断标准化feature之间差距
    slice1d<float> proj_features,			// Projector 生成 feature
    slice1d<float> proj_latent,				// Projector 生成 latent
    nnet_evaluation& evaluation,			// Projector 结构
    const slice1d<float> query,				// 当前 feature 值 （非标准化的值）
    const slice1d<float> features_offset,	// 当前 feature 值的偏移
    const slice1d<float> features_scale,	// 当前 feature 值的放大倍数
    const slice1d<float> curr_features,		// 当前 feature 值 （标准化的值）
    const nnet& nn,							// Projector 参数
    const float transition_cost = 0.0f)
{
    slice1d<float> input_layer = evaluation.layers.front();
    slice1d<float> output_layer = evaluation.layers.back();
    
    // 标准化 生成网络的输入
    for (int i = 0; i < query.size; i++){
        input_layer(i) = (query(i) - features_offset(i)) / features_scale(i);      
    }
    
    // 根据 Input_layer 跑 Projector 网络
    nnet_evaluate(evaluation, nn);
    
    // 拷出 Projector 网络输出的
    for (int i = 0; i < proj_features.size; i++){
        proj_features(i) = output_layer(i);
    }
    for (int i = 0; i < proj_latent.size; i++){
        proj_latent(i) = output_layer(proj_features.size + i);
    }
    
    // 计算更新之后的 feature 和更新前的 feature 之间的距离（标准化对比）
    best_cost = 0.0f;
    for (int i = 0; i < proj_features.size; i++){
        best_cost += squaref(query(i) - proj_features(i));
    }
    best_cost = sqrtf(best_cost);
    
    // 计算更新之后的 feature 和更新前的 feature 之间的距离（非标准化对比）
    float trns_dist_squared = 0.0f;
    for (int i = 0; i < proj_features.size; i++){
        trns_dist_squared += squaref(curr_features(i) - proj_features(i));
    }
    
    // 差距过大，则使用 P 输出的结果（这里用哪个差距都一样，都得转）
    if (trns_dist_squared > squaref(transition_cost)){
        // transition and add the transition cost
        transition = true;
        best_cost += transition_cost;
    }
    else {   
        // Don't transition and use current features as-is
        transition = false;
        
        for (int i = 0; i < proj_features.size; i++){
            proj_features(i) = curr_features(i);
        }
        
        // Re-compute the projection cost
        best_cost = 0.0f;
        for (int i = 0; i < curr_features.size; i++) {
            best_cost += squaref(query(i) - curr_features(i));
        }
        best_cost = sqrtf(best_cost);
    }
}

可以看出这里使用了两个更新前后 feature 之间的距离（非标准化和标准化），但实际上 transition_cost = 0，故无论如何都会置transition = true，并且两个距离变量 best_cost，transition_cost，在后边都没有使用到，所以这块代码是冗余的（我个人认为这是作者的失误）

如果要改我认为：由于 proj_features 和 curr_features 都为标准化的 feature 变量，故 trns_dist_squared 能更好的比较输出与输入之间的距离，故可以在函数输入 transition_cost 变量取一个阈值

Decompressor

decompressor_evaluate() 功能：输入 Stepper 或者Projector 生成的 features 和 latent ；输出此时应该做的 pose。具体的步骤可以表示为：

设 features 和 latent 为 Input
运行 Decompressor
备份输出

void decompressor_evaluate(
    slice1d<vec3> bone_positions,			// 输出骨架 positions 信息
    slice1d<vec3> bone_velocities,			// 输出骨架 velocities 信息
    slice1d<quat> bone_rotations,			// 输出骨架 rotations信息
    slice1d<vec3> bone_angular_velocities,	// 输出骨架 angular_velocities信息
    slice1d<bool> bone_contacts,			// 双脚离地了，病毒就关闭了，聪明的智商又占领高地了
    nnet_evaluation& evaluation,			// Decompressor 结构
    const slice1d<float> features,			// Decompressor 输入的 feature
    const slice1d<float> latent,			// Decompressor 输入的 latent
    const vec3 root_position,				// 上一时刻 root position  	（目标位置）
    const quat root_rotation,				// 上一时刻 root rotation	（目标位置）
    const nnet& nn,							// Decompressor 参数
    const float dt = 1.0f / 60.0f)
{
    slice1d<float> input_layer = evaluation.layers.front();
    slice1d<float> output_layer = evaluation.layers.back();
  
    // 设置网络的输入 features + latent
    for (int i = 0; i < features.size; i++) {
        input_layer(i) = features(i);
    }
    for (int i = 0; i < latent.size; i++) {
        input_layer(features.size + i) = latent(i);
    }
    
    // Decompressor 生成 pose
    nnet_evaluate(evaluation, nn);
    
    // 提取输出骨架 positions 信息
    int offset = 0;
    for (int i = 0; i < bone_positions.size - 1; i++) {
        bone_positions(i + 1) = vec3(
            output_layer(offset+i*3+0),
            output_layer(offset+i*3+1),
            output_layer(offset+i*3+2));
    }
    offset += (bone_positions.size - 1) * 3;
    
    // 提取输出骨架 rotations 信息，并转换为 4元数形式
    for (int i = 0; i < bone_rotations.size - 1; i++) {   
        bone_rotations(i + 1) = quat_from_xform_xy(
            vec3(output_layer(offset+i*6+0),
                 output_layer(offset+i*6+2),
                 output_layer(offset+i*6+4)),
            vec3(output_layer(offset+i*6+1),
                 output_layer(offset+i*6+3),
                 output_layer(offset+i*6+5)));
    }
    offset += (bone_rotations.size - 1) * 6;
    
    // 提取输出骨架 velocities 信息
    for (int i = 0; i < bone_velocities.size - 1; i++) {
        bone_velocities(i + 1) = vec3(
            output_layer(offset+i*3+0),
            output_layer(offset+i*3+1),
            output_layer(offset+i*3+2));
    }
    offset += (bone_velocities.size - 1) * 3;
    
    // 提取输出骨架 angular velocities 信息，并转换为 4元数形式
    for (int i = 0; i < bone_angular_velocities.size - 1; i++) {
        bone_angular_velocities(i + 1) = vec3(
            output_layer(offset+i*3+0),
            output_layer(offset+i*3+1),
            output_layer(offset+i*3+2));
    }
    offset += (bone_angular_velocities.size - 1) * 3;
    
    // 提取输出 root 位置的 velocities 和 angular velocity 信息，并转换为世界坐标 
    // 即两个速度信息代表着目标所在的位置
    vec3 root_velocity = quat_mul_vec3(root_rotation, vec3(
        output_layer(offset+0),
        output_layer(offset+1),
        output_layer(offset+2)));     
    vec3 root_angular_velocity = quat_mul_vec3(root_rotation, vec3(
        output_layer(offset+3),
        output_layer(offset+4),
        output_layer(offset+5)));
    offset += 6;

	// 通过 root 的 velocity 信息更新 root_position
    bone_positions(0) = dt * root_velocity + root_position;
    bone_velocities(0) = root_velocity;
    
    // 通过 root 的 angular velocity 信息更新 bone_rotations
    bone_rotations(0) = quat_mul(quat_from_scaled_angle_axis(root_angular_velocity * dt), root_rotation);
    bone_angular_velocities(0) = root_angular_velocity;    
    
    // 提取双脚状态 	输出 > 1/2 认为接触地面 
    //				输出 < 1/2 认为未接触地面 
    if (bone_contacts.data != nullptr) {
        bone_contacts(0) = output_layer(offset+0) > 0.5f;
        bone_contacts(1) = output_layer(offset+1) > 0.5f;
    }
    offset += 2;
    
    // Check we got everything!
    assert(offset == nn.output_mean.size);
}

Stepper

stepper_evaluate() 功能：更新输入的 features 和 latent

设 features 和 latent 为 Input
运行 Stepper
更新 features 和 latent

void stepper_evaluate(
    slice1d<float> features,			// 输入的 features 并且在函数中更新
    slice1d<float> latent,				// 输入的 latent   并且在函数中更新
    nnet_evaluation& evaluation,		// Stepper 结构
    const nnet& nn,						// Stepper 参数
    const float dt = 1.0f / 60.0f)
{
    slice1d<float> input_layer = evaluation.layers.front();
    slice1d<float> output_layer = evaluation.layers.back();
  
    // 设置 input 为 features 和 latents
    for (int i = 0; i < features.size; i++) {
        input_layer(i) = features(i);
    }    
    for (int i = 0; i < latent.size; i++) {
        input_layer(features.size + i) = latent(i);
    }
    
    // Stepper 生成输出 features 和 latent 的变化量
    nnet_evaluate(evaluation, nn);
    
    // 更新 features 和 latents 
    for (int i = 0; i < features.size; i++){
        features(i) += dt * output_layer(i);
    }    
    for (int i = 0; i < latent.size; i++) {
        latent(i) += dt * output_layer(features.size + i);
    }
}

可以看到，Stepper 的输出不是 features 和 latent 而是其变化量，这个在训练的时候会提到，这里只是使用它，搞工程就是这样哈哈哈哈哈哈

主循环

开始我们的主循环，我设置的是1s 右走 1s 左跑，上代码

首先（加速标志位）判断该时刻是加速 or 减速来对 desired_gait 进行更新（存在阻尼）
这个（加速标志位）在程序中就是阻尼器的 “目标位置” ，可以在函数 simple_spring_damper_implicit() 输入的第三位找到

// Get the desired gait (walk / run)
desired_gait_update(
	desired_gait,
	desired_gait_velocity,
	dt,
	tt);

其中

void desired_gait_update(
	float& desired_gait,
	float& desired_gait_velocity,
	const float dt,
	int tt,
	const float gait_change_halflife = 0.1f)
{
	simple_spring_damper_implicit(
		desired_gait,
		desired_gait_velocity,
		(tt / FPS) % 2,    //  这里来判断是否加速   我这里写的是单数秒的时候加速
		gait_change_halflife,
		dt);
}

这是阻尼器的一个最简单的实现，对于单个 float 变量加上了阻尼，在 “单数秒的时候加速” 这个定义下，desired_gait 随时间变化如下所示：

现在有了 gait 步态信息，直接在0~1之间线性插值就能计算出来 desired 的速度了

// Get the desired simulation speeds based on the gait
// 方向速度 =  (跑步速度 - 走路速度) * desired_gait + 走路速度 
float simulation_fwrd_speed = lerpf(simulation_run_fwrd_speed, simulation_walk_fwrd_speed, desired_gait);
float simulation_side_speed = lerpf(simulation_run_side_speed, simulation_walk_side_speed, desired_gait);
float simulation_back_speed = lerpf(simulation_run_back_speed, simulation_walk_back_speed, desired_gait);

通过 desired 速度计算世界坐标下的 desired 速度和方向（我们简化版本在不考虑视角旋转的情况下其实可以不用干）

// Get the desired velocity
vec3 desired_velocity_curr = desired_velocity_update(...);
// Get the desired rotation/direction
quat desired_rotation_curr = desired_rotation_update(...);

使用生成的 desired 变量： desired_rotation 和 desired_velocity
来生成未来20 40 60 帧的轨迹： trajectory_desired_xxxx 和 trajectory_xxxx

值得一提的是，trajectory_desired_xxxx 的预测不需要阻尼，纯粹的由用户 Input 计算生成，但是 trajectory_xxxx 是需要阻尼的，因为真实人物的运动是有惯性的。

// Predict Future Trajectory
// 期望未来20 40 60 帧的面朝方向
trajectory_desired_rotations_predict(...);
// 实际算出未来20 40 60 帧的面朝方向应该是什么样的（存在阻尼）
trajectory_rotations_predict(...);
// 期望未来20 40 60 帧的速度
trajectory_desired_velocities_predict(...);
// 实际算出未来20 40 60 帧的应该在的位置 （存在阻尼）
// 这个计算根据的是 模拟的位置、速度、加速度、轨迹上期望的速度
trajectory_positions_predict(...);

到这里就完成了准备工作，即动画上显示的运动方向

图中紫色的圆圈和方向就是我们算出来用户输入的 trajectory_desired_xxxx ；灰色的圆圈就是我们通过阻尼算出来实际人物应该运动的 trajectory_xxxx

程序中将 feature 存储在 query 中，如果是 LMM 就是使用 features_curr 作为 feature；如果使用 MM 就是使用 frame_index 在数据库中寻找 feature

array1d<float> query(db.nfeatures());

// feature 选择
// 如果使用 LMM 就是使用 features_curr 生成
// 如果使用  MM 就是使用 frame_index 在数据库中寻找 feature 
slice1d<float> query_features = lmm_enabled ? slice1d<float>(features_curr) : db.features(frame_index);

int offset = 0;
query_copy_denormalized_feature(query, offset, 3, query_features, db.features_offset, db.features_scale); // Left Foot Position
query_copy_denormalized_feature(query, offset, 3, query_features, db.features_offset, db.features_scale); // Right Foot Position
query_copy_denormalized_feature(query, offset, 3, query_features, db.features_offset, db.features_scale); // Left Foot Velocity
query_copy_denormalized_feature(query, offset, 3, query_features, db.features_offset, db.features_scale); // Right Foot Velocity
query_copy_denormalized_feature(query, offset, 3, query_features, db.features_offset, db.features_scale); // Hip Velocity
query_compute_trajectory_position_feature(query, offset, bone_positions(0), bone_rotations(0), trajectory_positions);
query_compute_trajectory_direction_feature(query, offset, bone_rotations(0), trajectory_rotations);

其中 query_copy_denormalized_feature() 函数是去归一化作用的，query_features 是归一化的 features ，features_scale 是放大的倍数，features_offset是偏移量，即：

$query = query_{features} * features_{scale}+ features_{offset}$

现在已经求得了当前状态的 feature 了，接下来如何生成下一时刻的 feature 呢？之前在 LMM 中有讲过，虽然说所有的下一时刻的状态都是靠 Stepper 生成的，但有两种情况需要重新匹配合适的 feature，对当前的 feature 进行修正：

输入的 feature 与当前动作的 feature 差异过大
靠 Stepper 生成下一步这个方式，连续使用太多步

// Check if we should force a search because input changed quickly
desired_velocity_change_prev = desired_velocity_change_curr;
desired_velocity_change_curr = (desired_velocity_curr - desired_velocity) / dt;
desired_velocity = desired_velocity_curr;

desired_rotation_change_prev = desired_rotation_change_curr;
desired_rotation_change_curr = quat_to_scaled_angle_axis(quat_abs(quat_mul_inv(desired_rotation_curr, desired_rotation))) / dt;
desired_rotation = desired_rotation_curr;

bool force_search = false;

if (force_search_timer <= 0.0f && (		// 第二种情况
	(length(desired_velocity_change_prev) >= desired_velocity_change_threshold && 
	length(desired_velocity_change_curr) < desired_velocity_change_threshold) || // 第一种情况
	(length(desired_rotation_change_prev) >= desired_rotation_change_threshold &&// 
	length(desired_rotation_change_curr) < desired_rotation_change_threshold)))// 第一种情况
{
	force_search = true;				// 需要重新匹配
	force_search_timer = search_time;	// 重新开始计算连续使用步数
}
else if (force_search_timer > 0){
	force_search_timer -= dt;			// S 网络生成下一步 且连续使用步数--
}

在代码中， force_search 用来判断是否需要 Projector 重新匹配，还是直接 Stepper 生成下一动作

// 是否需要匹配
if (force_search || search_timer <= 0.0f || end_of_anim) { 
	if (lmm_enabled) { 		// LMM 方法
		// 初始化 cost
		float best_cost = FLT_MAX;
		bool transition = false;
		
		// Projector 生成 pose
		projector_evaluate(...); 	

		// Projector 生成前后差距大，需要修正当前 feature
		if (transition) {
			// 生成 pose
			decompressor_evaluate(...);

			// 姿态更新
			inertialize_pose_transition(...);

			// 更新 current features and latents
			features_curr = features_proj;
			latent_curr = latent_proj;
		}
	}
	else{ ... }		 // MM 方法 

	// Reset search timer
	search_timer = search_time;
}

这一段代码的做法：

通过 Projector 生成一个 feature，如果生成的 feature 与当前 feature 距离大过阈值，则需要 transition = true去修正当前的 feature，否则认为生成的 feature 与当前 feature 差不多，没必要浪费时间去修正
将生成的 feature 传入至 Decompressor 产生一个 pose
利用这个 pose，在 inertialize_pose_transition()，中对姿态进行更新

PS： 只有差异过大（输入相差甚远或动作变化）这个时候 Projector 就发挥作用了，此时 LMM 做法对应于 MM 中即是换一个 dataset 源，接着那个dataset 源继续运动。

之前 search 只是用来判断模态有着较大的改变需要使用 Projector 对 feature 进行修正，但无论修正与否，其结果是停留在 t 时刻的，但是在这个循环中，我们需要的是 t+1 时刻的 feature ，所以说： “ Projector 都已经生成了 feature 了，就不用 Stepper 生成啦” 这是一个误区。

接着就是 Stepper 生成下一时刻的 feature ：

// Tick down search timer
search_timer -= dt;
if (lmm_enabled){	// LMM 方法
	// Stepper 更新下一时刻的 feature 和 latent
	stepper_evaluate(...);

	// Decompress 生成 pose
	decompressor_evaluate(...);
}
else{ ... }			// MM 方法

// 更新姿态
inertialize_pose_update(...);

到这里主循环就全部结束了，满足了最简单的 LMM ，当前此时还存在着许多问题，比如 IK 进行锁足功能之类的

五、LMM Train

由于原代码过分的冗余，重复的东西太多了，所以这里就以伪代码形式展示

Decompressor

输入：Feature 、Pose、Compressor 参数、Decompressor 参数
方法：类比为一个生成网络，Compressor 作为 encode；Decompressor 作为 decoder
评价标准：
- latent：稀疏性（ sparsity）& 平滑性（ smoothness）
- pose：生成的精度；平滑性（ smoothness）（位置 & 速度）
输出：Decompressor 网络参数

Train_Decompressor(Feature, Pose, θC, θD):
	Pose_Glob ← Forward_Kinematics(Pose)		# 生成全局 pose
	Latent ← C([Pose Pose_Glob]; θC)			# Compressor 生成 Latent 变量
	Pose_D ← D([X Z]T; θD)						# Decompressor 生成重建的 pose
	Pose_D_Glob ← Forward_Kinematics(Pose_D)	# 生成全局重建的 pose
	
	# 计算 Loss
	L_lreg ← w_1 * square(Latent)				# 平方				
	L_sreg ← w_2 * abs(Latent)					# 绝对值之和
	L_vreg ← w_3 * abs(△Latent / dt)			# 变化速度的绝对值之和
	
	L_loc ← w_4 * abs(Pose - Pose_D)			# 包括 joint 4个参数 root 的速度2个参数				
	L_chr ← w_5 * abs(Pose_Glob - Pose_D_Glob)	
	L_lvel ← w_6 * abs(△Pose / dt - △Pose_D/ dt)	
	L_cvel ← w_7 * abs(△Pose_Glob / dt - △Pose_D_Glob/ dt)	

	# loss 的总和
	Loss ← SUM(L_lreg + ... + L_cvel)

	# 更新参数
	θC, θD ← RAdam(θC θD, Loss)

Stepper

输入：Feature_In、Latent_In、生成次数 s 、Stepper 参数
方法：Stepper 生成 s 步，使生成的 Feature 和 Latent 尽可能与真实值吻合
评价标准：平滑性（位置 & 速度）
输出：Stepper 网络参数

需要注意的是这里的变化量是训练网络中的帧数决定的，所以使用 Stepper 的时候需要注意一下使用的 dt 即帧数的问题。

Train_Stepper(Feature_In, Latent_In, s, θS):
	# 备份
	Feature_S[0], Latent_S[0] ← Feature_In[0], Latent_In[0]

	# Stepper 生成 s步
	for i ← 1 to s
		# Stepper 生成变化量
		Feature_delta, Latent_delta ← S([Feature_In[i-1], Latent_In[i-1]]; θS)
		# 更新
		Feature_S[i], Latent_S[i] ← [Feature_S[i-1], Latent_S[i-1]] + [Feature_delta, Latent_delta]
	
	#计算loss
	L_xval ← w_1 * abs(Feature_In - Feature_S)	
	L_zval ← w_2 * abs(Latent_In- Latent_S)	
	L_xvel ← w_3 * abs(△Feature_In / dt - △Feature_S / dt)	
	L_zvel ← w_4 * abs(△Latent_In/ dt - △Latent_S/ dt)	
	
	# loss 的总和
	Loss ← SUM(L_xval + ... + L_zvel)
	
	# 更新参数
	θS ← RAdam(θS; Loss)

Projector

输入：Feature_1farme、Feature、Latent、Projector参数
方法：对 Feature 中选定一帧 Feature_1farme 进行加噪，输入 Projector 生成的Feature_P 和 Latent_P 尽可能和最近邻的 Feature_k 和 Latent_k 吻合
评价标准：数值接近（最重要），噪声的消去项
输出：Projector网络参数

Train_Projector(Feature_1farme, Feature, Latent, θP):
	# 生成一个噪声（均匀噪声 & 高斯噪声）
	N_U, N_G ← U(0, 1), N(0, 1)
	
	# 添加噪声 (方差 * 均匀噪声 * 高斯噪声)
	Feature_Noise ← Feature_1farme + (Feature_std * N_U * N_G)

	# 寻找最近邻的 Feature 
	Feature_k, Latent_k ← Nearest(Feature_Noise, [Feature Latent])

	# Projector 生成 Feature 和 Latent
	Feature_P, Latent_P ← P(Feature_Noise; θP)

	# 计算 Loss
	L_xval ← w_1 * abs(Feature_k - Feature_P)	# Feature 真实值匹配
	L_zval ← w_2 * abs(Latent_k - Latent_P)		# Latent 真实值匹配
	L_dist ← w_3 * abs(dist(Feature_Noise, Feature_k) - dist(Feature_Noise, Feature_P)) # 防止噪声太大的选项

	# loss 的总和
	Loss ← SUM(L_xval + L_zval + L_dist )

	# 更新参数
	θP ← RAdam(θP; Loss)

5.20 更新

六、LOCK FOOT

两关节的 Inverse Kinematics

这里借助大佬的文章：

链接: Simple Two Joint IK

两个关节的 IK 可以分为两步来做：

扩展 / 收缩关节链，使得 AC = AT
旋转 A 节点，使得末端滑入目标位置

这两步操作可以表示为：
在这里插入图片描述
Step 1: 扩展 / 收缩关节链

首先先分析当前各个关节点的角度：

$\angle CAB = acos(\frac{AB \cdot AC}{||AB|| \; ||AC||})\\ \angle ABC = acos(\frac{BA \cdot BC}{||BA|| \; ||BC||})$

再根据余弦定理推导旋转后各个关节点的角度：

$\angle CAB' = acos(\frac{BC^2-AB^2 - AT^2}{-2*AB*AT|})\\ \angle ABC' = acos(\frac{AT^2-BA^2 - BC^2}{-2*BA*BC})$

旋转角这不就出来了？

$\angle r0 =\angle CAB' - \angle CAB \\ \angle r1 = \angle ABC' - \angle ABC$

Step 2: 旋转 A 节点

这一步就更简单了，旋转角就摆在脸上就是∠TAC

$\angle r2 = \angle TAC = acos(\frac{AT \cdot AC}{||AT|| \; ||AC||})$

只不过程序中使用的是四元数，所以后续还需要将其转化为四元数表示旋转，并且由于浮点运算的误差，所以程序中还使用了clamp() 来保证余弦结果在 (-1, 1) 之间

下面列出了两个关节 IK 的完整代码。

void two_joint_ik(
    vec3 a, vec3 b, vec3 c, vec3 t, float eps,
    quat a_gr, quat b_gr, 
    quat &a_lr, quat &b_lr) {
  
    float lab = length(b - a);
    float lcb = length(b - c);
    float lat = clamp(length(t - a), eps, lab + lcb - eps);

    float ac_ab_0 = acos(clamp(dot(normalize(c - a), normalize(b - a)), -1, 1));
    float ba_bc_0 = acos(clamp(dot(normalize(a - b), normalize(c - b)), -1, 1));
    float ac_at_0 = acos(clamp(dot(normalize(c - a), normalize(t - a)), -1, 1));

    float ac_ab_1 = acos(clamp((lcb*lcb-lab*lab-lat*lat) / (-2*lab*lat), -1, 1));
    float ba_bc_1 = acos(clamp((lat*lat-lab*lab-lcb*lcb) / (-2*lab*lcb), -1, 1));

    vec3 axis0 = normalize(cross(c - a, b - a));
    vec3 axis1 = normalize(cross(c - a, t - a));    
    
    quat r0 = quat_angle_axis(ac_ab_1 - ac_ab_0, quat_mul(quat_inv(a_gr), axis0));
    quat r1 = quat_angle_axis(ba_bc_1 - ba_bc_0, quat_mul(quat_inv(b_gr), axis0));
    quat r2 = quat_angle_axis(ac_at_0, quat_mul(quat_inv(a_gr), axis1));

    a_lr = quat_mul(a_lr, quat_mul(r0, r2));
    b_lr = quat_mul(b_lr, r1);
}

锁足代码

// Contact fixup with foot locking and IK

adjusted_bone_positions = bone_positions;
adjusted_bone_rotations = bone_rotations;

if (ik_enabled) {
	for (int i = 0; i < contact_bones.size; i++) {
		int toe_bone = contact_bones(i);				// 脚趾   	index
		int heel_bone = db.bone_parents(toe_bone);		// 脚踝   	index
		int knee_bone = db.bone_parents(heel_bone);		// 膝关节 	index
		int hip_bone = db.bone_parents(knee_bone);		// 髋关节 	index
		int root_bone = db.bone_parents(hip_bone);		// root		index

		// 通过 FK 计算 toe 至 root 所有关节的全局坐标
		forward_kinematics_partial(...);

		// 判断脚是否接地？ 
		// 如果是则将锁足标志置1 并且给一个锁足点
		contact_update(...);

		// 保证锁足点在地上
		vec3 contact_position_clamp = contact_positions(i);
		contact_position_clamp.y = maxf(contact_position_clamp.y, ik_foot_height);

		// 两个关节的 IK 
		// target 为接触点减去一个 "脚向量"
		ik_two_bone(...);

	 	// 根据 IK 调整姿态
		for (int bone : {toe_bone, heel_bone, knee_bone}) {
			forward_kinematics_partial(...);
		}

		// 保证 toe 保持在地面上
		ik_look_at(...);	// heel
		for (int bone : {toe_bone, heel_bone}) {
			forward_kinematics_partial(...);
		}
		ik_look_at(...);	// toe
	}
}

// 计算所有关节点的全局坐标
forward_kinematics_full(...);