基于mediapipe的动作捕捉和人物模型的驱动

我实习做这个的时候，当时也找了很多资料，但很多资料都没有讲透，下面我就总结一下自己结合看代码后的感悟，还有一些我认为很重要的点。

qq_58484580

3378人浏览 · 2023-09-05 17:15:37

qq_58484580 · 2023-09-05 17:15:37 发布

我实习做这个的时候，当时也找了很多资料，但很多资料都没有讲透，下面我就总结一下自己结合看代码后的感悟，还有一些我认为很重要的点。

一.任务难点：

1.关键点如何检测？

2.检测出的关键点数据如何处理？

3.骨骼关键点数据和人物模型如何对应？

4.人物模型的关节如何约束？

二.解决方法：

1.关键点检测：

这里用到的是mediapipe，其中pose_landmark函数既能输出2d的关键点坐标，又能输出3d的关键点的坐标，2d关键点的坐标是以图像的左上角为原点，而3d关键点的坐标是以hips为原点，也就是图1中的left_hip和right_hip中间为原点。如果直接以3d关键点为输入的话，是没有全局的平移和旋转的。现在已知3d空间中的关键点和2d投影的点，可以用PNP算法求解相机的位姿，进而可以算出3d关键点相对于相机的位置变化。

图1

我第一次接触到3维视觉中的世界坐标系，相机坐标系，图像坐标系，像素坐标系也一脸懵，这个PNP算法其实就用到坐标系的转换。

1.1.3D视觉中的坐标系:

1.1.1世界坐标系:

世界坐标系就是图2中的Xw-Yw-Zw坐标系，简单理解就是我们生活的真实3维世界。

图2

1.1.2相机坐标系:

图2中的Xc-Yc-Zc坐标系，相当于世界坐标系任意旋转和平移，也可以理解为真实世界中有一个相机。

1.1.3图像坐标系:

图2中的U-V坐标系，去掉了Z轴，就相当于把相机坐标系的点投影在一个平行于相机坐标系的Xc-和Yc的平面。

1.1.4像素坐标系:

图2中的X-Y的坐标系，是以像素为单位，算相对于原点的偏移。

1.2.坐标系之间的转换:

世界坐标系-相机坐标系：

相机坐标系-图像坐标系:

这里面要用到一个f，表示相机的焦距，也就是相机的内参。

图像坐标系-相机坐标系:

PNP算法只需额外提供一个夹角，就可以求解相机坐标系相对于世界坐标系的位姿，也就是R旋转矩阵和T平移矩阵

相机坐标系中的点投影到平面上和2d点的x，y对不上，为了两者能保持一致，就保持相机坐标的z坐标不变，直接调整x,y使得投影结果一致。

2.关键点的数据处理:

首先介绍一下为什么要处理关键点数据，通过观察我们发现mediapipe检测人体关键点时，有些关键点会发生高频抖动，为了消除抖动，以及关键点位移更加平滑，我们采用滤波方法。这里的滤波方法有很多比如滑动均值滤波法，卡尔曼滤波法，euro滤波法，移动最小二乘法等等。这里我就简单的介绍一下:

2.1.滑动均值滤波法:

设定采样间隔，取间隔内的均值即可。例如：每隔5个采样点取一次均值，此时输出值就是前五个的3d坐标点的均值。代码如下：

# 可根据需要调整窗口大小
window_size = 5  
# 初始化滑动均值滤波的队列
keypoint_queue = []
# 处理每一帧图像中的3D关键点
def process_3d_keypoints(frame_3d_keypoints):
    # 将新的3D关键点添加到队列
    keypoint_queue.append(frame_3d_keypoints)
    # 如果队列大小超过窗口大小，则移除最旧的元素
    if len(keypoint_queue) > window_size:
        keypoint_queue.pop(0)
    # 计算滑动均值
    smoothed_keypoints = np.mean(keypoint_queue, axis=0)
    # 在这里执行对平滑后的关键点的操作，例如可视化或其他后续处理
    return smoothed_keypoints

2.2.卡尔曼滤波法:

卡尔曼滤波法中，有两种变量，一个是预测的状态变量(公式，定理算出的值，表示着理想值)，一个是实际测量的观察变量(有误差)。卡尔曼滤波法得思想就是：预测的状态变量和实际测量的观察变量进行数据融合得到接近真实的数据，然后作为下一个过程的状态变量，再与观察变量融合，如此往复。

class KalmanFilterWrapper:
    def __init__(self, input_dim, init_error, init_process_var, init_measure_var):
        self.input_dim = input_dim
        self.init_error = init_error
        self.init_process_var = init_process_var
        self.init_measure_var = init_measure_var
        self.kf,self.statapose= self._create_kalman_filter()

    def _create_kalman_filter(self):
        kf = cv2.KalmanFilter(self.input_dim, self.input_dim)
        kf.transitionMatrix = np.eye(self.input_dim,dtype=np.float32)
        kf.measurementMatrix = np.eye(self.input_dim,dtype=np.float32)
        kf.processNoiseCov = self.init_process_var * np.eye(self.input_dim,dtype=np.float32)
        kf.measurementNoiseCov = self.init_measure_var * np.eye(self.input_dim,dtype=np.float32)
        kf.errorCovPost = self.init_error * np.eye(self.input_dim,dtype=np.float32)
        return kf,kf.statePost

    def filter(self, observation):
        return self.kf.correct(observation)

    def predict(self):
        return self.kf.predict()

以上的程序展示的是1D卡尔曼滤波器，系统的初始化状态的设置是第一次检测时关键点的坐标，系统的预测方差:表示了在卡尔曼滤波器中考虑的系统模型的不确定性.系统的测量方差:表示了传感器(测量工具)测量的不确定性。系统的不确定性大，那么预测方差可以调高一些。测量工具存在噪声，测量方差可以调大一些。初始化误差和但两个都要根据实际情况来调整。Input_dim表示要打散成一维。这里的卡尔曼滤波器，使用了三个，分别对躯体还有左手和右手创建一个卡尔曼滤波器，有些部分的相关性强，有的弱，这样分开滤波效果会好一些。

2.3.euro滤波法:

euro是一种能够实时滤除噪声的低通滤波器，而且只需要调整两个参数（最小截止频率和速度系数），一般需要调节这两个参数达到动态平衡的效果，降低最小截至频率可以大幅度的降低抖动，但是会出现延迟，增加速度系数可以降低延迟，但是去抖效果会大大降低。

下面就是相应的代码:

class OneEuroFilter:

    def __init__(self, t0, x0, dx0=0.0, min_cutoff=1.0, beta=0.0, d_cutoff=1.0):
        """Initialize the one euro filter."""
        # The parameters.
        self.min_cutoff = float(min_cutoff)
        self.beta = float(beta)
        self.d_cutoff = float(d_cutoff)

        # Previous values.
        self.x_prev = x0
        self.dx_prev = float(dx0)
        self.t_prev = t0

    def smoothing_factor(self, t_e, cutoff):
        r = 2 * math.pi * cutoff * t_e
        return r / (r + 1)

    def exponential_smoothing(self, a, x, x_prev):
        return a * x + (1 - a) * x_prev

    def filter_signal(self, t, x):
        """Compute the filtered signal."""
        t_e = t - self.t_prev

        # The filtered derivative of the signal.
        a_d = self.smoothing_factor(t_e, self.d_cutoff)
        dx = (x - self.x_prev) / t_e
        dx_hat = self.exponential_smoothing(a_d, dx, self.dx_prev)

        # The filtered signal.
        cutoff = self.min_cutoff + self.beta * abs(dx_hat)
        a = self.smoothing_factor(t_e, cutoff)
        x_hat = self.exponential_smoothing(a, x, self.x_prev)

        # Memorize the previous values.
        self.x_prev = x_hat
        self.dx_prev = dx_hat
        self.t_prev = t

        return x_hat

下面表示整个代码的逻辑：

2.4.移动最小二乘法(MLS):

在MLS法当中，需要在一组不同位置的节点（node）附近建立拟合曲线，每个节点都有自己的一组系数值aj，每个节点的系数(aj)取值只考虑其临近采样点，且距离节点越近的采样点的贡献就越大。

def mls_smooth_numpy(input_t: List[float], input_y: List[np.ndarray], query_t: float, smooth_range: float):
    # 1-D MLS: input_t: (N), input_y: (..., N), query_t: scalar
    if len(input_y) == 1:
        return input_y[0]
    # input_t:和当前帧的时间差
    input_t = np.array(input_t) - query_t
    # input_y:
    input_y = np.stack(input_y, axis=-1)
    broadcaster = (None,)*(len(input_y.shape) - 1)
    # w: 每个节点对采样点的权值
    w = np.maximum(smooth_range - np.abs(input_t), 0)
    # input_t[broadcaster]:(1,len(input_t))--2维
    # w[broadcaster]:(1,len(input_t))---2维
    coef = moving_least_square_numpy(input_t[broadcaster], input_y, w[broadcaster])
    return coef[..., 0]

def moving_least_square_numpy(x: np.ndarray, y: np.ndarray, w: np.ndarray):
    # 1-D MLS: x: (..., N), y: (..., N), w: (..., N)
    # p:(1,2,len(input_t))
    p = np.stack([np.ones_like(x), x], axis=-2)             # (..., 2, N)
    M = p @ (w[..., :, None] * p.swapaxes(-2, -1))
    a = np.linalg.solve(M, (p @ (w * y)[..., :, None]))
    a = a.squeeze(-1)
    return a

MLS的效果和处理数据的长度有关，也和每个采样点的权值大小有关。

2.5.引入质心:

虽然运用上边的各种滤波方法对关键点进行滤波平滑处理，但是人体还是会有整体的晃动，特别是有大幅动作时前后晃动，看起来很不自然。看到有一个大佬的解决方法，就是对各个关键点加以权值，然后通过加权求和，求出人体近似的一个质心，对质心再利用上述的平滑方法，对质心的轨迹做强有力的平滑，对全身的大范围动作具有稳定的作用。那各个关键点的权值的确定，可以根据距离人体的实际质心的远近来确定，但最后所有的权值相加必须等于1。也就是说距离人体中心越近的地方的权值越大(例如：left_hip，right_hip这两点距离重心近可以加大的权值)，相反越小。

3.骨骼关键点和人物模型的对应:

3.1.关键点对应:

图3

图3显示的就是unity中的骨骼关键点，但是可能和mediapipe中的关键点的数目无法对应起来，也就是说unity中定义了的mediapipe中没定义，mediapipe中定义了的mediapipe中可能没有定义。所以我们有些没有值的就要自己根据比例来算：

图4

就比如，mediapipe中没有hips，spine，chest，upper chest，neck这些等等，那这些关键点坐标的计算都要依靠已有骨骼的位置。比如我想计算hips，那么需要粗略测量模型的hips相对于left_hip和right_hip的位置，然后根据实际的left_hip和right_hip再进行缩放和偏移。其他的骨骼计算都是如此。还有mediapipe中的3d坐标都是根据hips而移动的，所以在unity中只需要计算hips的偏移就保证人体的移动，其它关节就绕父关节旋转就行，这样就实现了人体的各种动作。

3.2.骨骼长度的对应:

这个也很好想，每个人高矮胖瘦不同，那么骨骼的长度也不同，mediapipe中检测的也就不同。同样虚拟人物模型的比例也不相同（尤其是一些卡通人物），这时为了更好的驱动人物模型，要身体各部位的关键点和人物模型对应起来，也要考虑要目标关键点之间相互耦合，人体的骨骼也有耦合，不能分若干部位独立求解，要整体优化。

3.3.骨骼旋转：

在unity中要自己一个个计算关键点的坐标以及关节与关节之间形成的夹角。主要用到了四元数矩阵，来表示骨骼相对于父骨骼的旋转，为什么不用欧拉角来旋转？是因为会出现万向锁的问题，我这里就是简单解释一下：例如：xyz三个轴，绕y转90度，那么就会导致另外两轴x，z重合，使得不知道最后是绕x轴旋转，还是z轴。当然解决的方法也有很多，四元数就是其中之一。具体的可以看下面的链接：

https://blog.csdn.net/euphorias/article/details/123612227

unity中的旋转主要用到四元数矩阵。

Quaternion Inverse(Quaternion rotation)

旋转方式：

Quaternion LookRotation(Vector3 forward, [DefaultValue("Vector3.up")] Vector3 upwards);

下面是求中间矩阵：

中间矩阵就相当于求骨骼转到和父骨骼相同角度所需要的角度

hip.Inverse = Quaternion.Inverse(Quaternion.LookRotation(forward));
hip.InverseRotation = hip.Inverse * hip.InitRotation;

然后再计算骨骼的旋转，最后在通过中间矩阵，得到骨骼最后旋转的位置。

root = animator.GetBoneTransform(HumanBodyBones.Hips);
midRoot = Quaternion.Inverse(root.rotation) * Quaternion.LookRotation(forward);

上面是需要一个个求角度，比较麻烦，网上有个大佬直接拿pytorch中的优化器来写的一个角度求解的优化器，关节与关节之间不独立求解角度，而是整体求解。这也是我觉得大佬整个工程中的精髓所在，我着重讲解。

3.4.整体求解：

大佬用了L-BFGS算法，在有限内存下进行BFGS算法求解，比较适合在大规模的数值计算中，具备牛顿法收敛速度快的特点，但不需要牛顿法那样存储Hesse矩阵，因此节省了大量的空间以及计算资源。

3.4.1牛顿法求解：

牛顿法求根:

a.在x轴任意选一点x1，经过x1做x的垂线，得到图像的交点f(x1)

b.过f(x1)做函数的切线，得到切线与x轴的交点x2。

c.重复a,b两个步骤。

d.经过不断迭代,xn就等于f(x)的根。

牛顿法是无限逼近根，当上一个Xk-1和Xk的差值小于阈值时，就可以认为Xk是f(x)的根

通过f(Xn)导数，求根Xn=Xn-1- ( f(Xn-1)/Xn-1 )

牛顿法求驻点的本质：二阶泰勒展开式，导函数为0。得到X = Xk - [ f '(Xk)/f ''(Xk) ]

在机器学习中，x不是一个数字而是一个向量，f '(Xk)也是一个向量，f ''(Xk)是一个矩阵也叫海塞矩阵。所以上边的式子就变成 $X_{k+1}=X_{k}-H_{k}^{-1}*g_{k}$ ，其中 $g_{k}$ 表示 $f_{}^{'}(x)$ ， $H_{k}^{-1}$ 表示二阶导函数的倒数。但是X维度多的时候，二阶导数难求，而且要多次迭代。

BFGS算法的出现，解决了 $H_{k}^{-1}$ 难求的问题。本质是：通过迭代来逼近 $H_{k}^{-1}$ ，进而近似代替 $H_{k}^{-1}$

逼近法：

$D_{k+1}=\left [ I-\frac{S_{k}*Y_{k}^{T}}{Y_{k}^{T}*S_{k}} \right ]*D_{k}*\left [ I-\frac{Y_{k}*S_{k}^{T}}{Y_{k}^{T}*S_{k}} \right ]+\frac{S_{k}*S_{k}^{T}}{Y_{k}^{T}*S_{k}}$ 其中I是单位矩阵， $D_{k}$ 对应 $H_{k}^{-1}$ ， $S_{k}=X_{k+1}-X_{k}$ , $Y_{k}=g_{k+1}-g_{k}$ , $g_{k}$ 表示原函数的导函数。

步骤：

a.随机取一个x1，那么损失函数在x1的导数g1可求，D1的初始值为单位矩阵I。

b.x2=x1- $D_{1}$ *g1，那么x2可求，同样g2可求。

c.s1 = x2-x1, y1=g2-g1,带入公式可求D2。

d.x3 = x2 -D2*g2,可求x3，同样g3可求，带入公式求D3，然后再算x4。

e.重复上述步骤进行迭代，直到Dn达到阈值时，停止迭代，此时可近似看作 $H_{k}^{-1}$ 的值。

L-BFGS算法：

BFGS算法每次要储存Dk矩阵，需要用很大的内存空间，可能使计算机无法运行。因此L-BFGS算法只保存最近的m次迭代信息，之前的信息都丢掉(可以看作是BFGS算法的又一次近似)，从而大大减少数据的存储空间。

上面讲述了L-BFGS算法由牛顿法一步步推导的由来，其实也就是在牛顿法的基础上，做了代替，进而在运算速度上和运算的内存上做了优化。那个大佬把全身的骨骼联系起来，也就是父骨骼的运动会影响子骨骼，子骨骼又会影响孙骨骼。以此类推会把影响传递下去。这里使用L-BFGS算法优化的是每块骨骼每次的旋转的欧拉角，损失函数是每块骨骼旋转欧拉角的约束+弱的角度的l2正则化+骨骼旋转前后的长度(不变)

大佬是用blender去做的，我用的unity，当时有个想法就是能不能让unity中的骨骼旋转也用L-BFGS算法去优化。但是后来自己各种尝试就没有成功。下面就是我当时想的，利用unity中的localword矩阵和worldtolocal矩阵进行变换，p：parent，c：child，cc：grandson

4.人物模型的关节旋转角度的约束：

在Blender中，上面的L-BFGS算法里面，就做了角度的约束，自己定义绕x，y，z轴旋转的角度限制，当超出这些角度就会给出相应的惩罚，来抑制下一次角度的变化。

在Unity中，可以先计算关节的角度，对角度进行限制后，在转换成四元数矩阵。

5.代码部分：

python部分的代码：https://github.com/ykyk000

如有错误，欢迎大家指正！！！

AtomGit 开源协作平台测评赛

瓜分20万奖金获得内推名额丰厚实物奖励易参与易上手

更多推荐

【Spring Boot 】Spring Boot + HikariCP 连接池使用示例

文章目录示例工具版本HikariCP 依赖HikariCP 配置1. connectionTimeout2. minimumIdle3. maximumPoolSize4. idleTimeout5. maxLifetime6. autoCommitSpring Boot Data + HikariCP + MySQL示例测试应用程序1. 使用 Maven 命令2. 使用 Eclipse3. 使用