guodong's blog

PhD@zhejiang university
   

三维重建(4):SMPL模型源码阅读

SMPL论文阅读部分:点此

SMPL源码下载:http://smpl.is.tue.mpg.de/downloads   需要注册

本文主要讨论SMPL源码核心代码。需要储备的知识还有:三维重建基础知识,以及 chumpy阅读以及opendr阅读LBS DQS

SMPL文件主要包括 vert.py serialization.py lbs.py。下面将逐一说说明。

0、关节位置

SMPL模型关节点名称

self.j_names = {

0: ‘Pelvis’,
1: ‘L_Hip’,
2: ‘R_Hip’,
3:  ‘Spine1’,
4: ‘L_Knee’,
5:  ‘R_Knee’,
6:  ‘Spine2’,
7:  ‘L_Ankle’,
8:  ‘R_Ankle’,
9:  ‘Spine3’,
10: ‘L_Foot’,
11: ‘R_Foot’,
12: ‘Neck’,
13: ‘L_Collar’,
14: ‘R_Collar’,
15: ‘Head’,
16: ‘L_Shoulder’,
17: ‘R_Shoulder’,
18: ‘L_Elbow’,
19: ‘R_Elbow’,
20: ‘L_Wrist’,
21: ‘R_Wrist’,
22: ‘L_Hand’,
23: ‘R_Hand’, }

关节树图形:

 

1、posemapper .py

主要是实现关节角度到姿势混合形状的映射。

import chumpy as ch
import numpy as np
import cv2

class Rodrigues(ch.Ch):  #创建函数类
    dterms = 'rt'  #可微变量旋度
    
    def compute_r(self):  #计算旋转矩阵
        return cv2.Rodrigues(self.rt.r)[0] 
    
    def compute_dr_wrt(self, wrt): # 计算旋度
        if wrt is self.rt:
            return cv2.Rodrigues(self.rt.r)[1].T #返回旋度的转置


def lrotmin(p): 
    if isinstance(p, np.ndarray): # 判断p的类型是不是np,基本上都不是,跳过
        p = p.ravel()[3:] # 转化成3*n矩阵
        return np.concatenate([(cv2.Rodrigues(np.array(pp))[0]-np.eye(3)).ravel() for pp in p.reshape((-1,3))]).ravel()        
    if p.ndim != 2 or p.shape[1] != 3: # 判断姿势的维度,如果不是2维或第二维shape不为3
        p = p.reshape((-1,3)) # 将pose转换成(24,3)的形式
    p = p[1:]  # 去掉root出,提取23个关节处
    return ch.concatenate([(Rodrigues(pp)-ch.eye(3)).ravel() for pp in p]).ravel()  # 将旋转向量转化成旋转矩阵,并拼接在一起

def posemap(s):
    if s == 'lrotmin':
        return lrotmin
    else:
        raise Exception('Unknown posemapping: %s' % (str(s),))

 

2、LBS.py

就是实现LBS

from posemapper import posemap
import chumpy
import numpy as np

#计算模型的全局刚体变换,包括旋度和平移量。实现论文里的公式4
def global_rigid_transformation(pose, J, kintree_table, xp):
    results = {}
    # 看了后面parents 终于知道这个kintree_table是干啥的了:关节树,可以理解为各个关节相互依赖的关系。
    pose = pose.reshape((-1,3)) #静止的姿势参数:shape=(24,3) 可能都是0向量
    # 下面两行主要是变换列表,并剔除[0,0]坐标的异常值(可以认为是root的方向异常)。kintree_table的shape=(2,24)
    id_to_col = {kintree_table[1,i] : i for i in range(kintree_table.shape[1])}
    parent = {i : id_to_col[kintree_table[0,i]] for i in range(1, kintree_table.shape[1])}
    #如果xp是chumpy类型
    if xp == chumpy:
        from posemapper import Rodrigues # 则引入posemapper的函数 
        rodrigues = lambda x : Rodrigues(x) # 创建匿名类,ch对象
    else:
        import cv2
        rodrigues = lambda x : cv2.Rodrigues(x)[0]  #否则的话使用opencv的函数
    # stack纵向拼接,改变的是行数,类似于concat。xp在这里就相当于chumpy类
    with_zeros = lambda x : xp.vstack((x, xp.array([[0.0, 0.0, 0.0, 1.0]]))) 
    # 先确定根关节的结果:目前还不知道这个result有什么用。
    # hstack横向拼接,改变的是列数。输入的是pose的root节点和关节的root节点,最后输出是主对角线为1的4*4矩阵 
    # pose[0] 表示静止姿势的根关节位置
    results[0] = with_zeros(xp.hstack((rodrigues(pose[0,:]), J[0,:].reshape((3,1))))) 
    # 利用关节的依赖关系求得其他关节的result    
    for i in range(1, kintree_table.shape[1]): # i从1到23,代表每个关节的编号
        results[i] = results[parent[i]].dot(with_zeros(xp.hstack((  # 父关节的轴角乘以关节相对旋度等于当前关节的轴角
            rodrigues(pose[i,:]),
            ((J[i,:] - J[parent[i],:]).reshape((3,1)))  # 计算关节的相对位置
            ))))
    
    pack = lambda x : xp.hstack([np.zeros((4, 3)), x.reshape((4,1))])
    
    results = [results[i] for i in sorted(results.keys())]
    results_global = results

    if True:
        results2 = [results[i] - (pack(
            results[i].dot(xp.concatenate( ( (J[i,:]), 0 ) )))
            ) for i in range(len(results))]
        results = results2
    result = xp.dstack(results)
    return result, results_global

# 计算每个顶点的关节影响度的混合 
def verts_core(pose, v, J, weights, kintree_table, want_Jtr=False, xp=chumpy):
    A, A_global = global_rigid_transformation(pose, J, kintree_table, xp) #调用函数,生成全局和局部的旋转矩阵。A.shape=[4,4,24] A_global.shape=[24,4,4]
    T = A.dot(weights.T)  # weight.shape=[6890,24],生成结果的shape=[4,4,6890]
    
    rest_shape_h = xp.vstack((v.T, np.ones((1, v.shape[0]))))  # 静止姿势的shape
        
    v =(T[:,0,:] * rest_shape_h[0, :].reshape((1, -1)) + 
        T[:,1,:] * rest_shape_h[1, :].reshape((1, -1)) + 
        T[:,2,:] * rest_shape_h[2, :].reshape((1, -1)) + 
        T[:,3,:] * rest_shape_h[3, :].reshape((1, -1))).T # 应用公式2

    v = v[:,:3] 
    
    if not want_Jtr:
        return v
    Jtr = xp.vstack([g[:3,3] for g in A_global])
    return (v, Jtr)

3、vert.py

verts_decorated函数没有被用到就不过多注释了


import chumpy
import lbs
from posemapper import posemap
import scipy.sparse as sp
from chumpy.ch import MatVecMult

def ischumpy(x): return hasattr(x, 'dterms') #如果没有dterm属性就返回chumpy

# verts_decorated函数:从另一个SMPL模型中参数来创建新的SMPL模型实例
def verts_decorated(trans, pose,
    v_template, J, weights, kintree_table, bs_style, f,
    bs_type=None, posedirs=None, betas=None, shapedirs=None, want_Jtr=False):
    #各个参数的含义 

    for which in [trans, pose, v_template, weights, posedirs, betas, shapedirs]:
        if which is not None:
           assert ischumpy(which)  #判断是否是chumpy类型,如果不是则报错

    v = v_template

    if shapedirs is not None:
       if betas is None:
          betas = chumpy.zeros(shapedirs.shape[-1])
       v_shaped = v + shapedirs.dot(betas)
    else:
       v_shaped = v

    if posedirs is not None:
       v_posed = v_shaped + posedirs.dot(posemap(bs_type)(pose))
    else:
       v_posed = v_shaped

    v = v_posed

    if sp.issparse(J):
       regressor = J
       J_tmpx = MatVecMult(regressor, v_shaped[:,0])
       J_tmpy = MatVecMult(regressor, v_shaped[:,1])
       J_tmpz = MatVecMult(regressor, v_shaped[:,2])
       J = chumpy.vstack((J_tmpx, J_tmpy, J_tmpz)).T
    else:
       assert(ischumpy(J))

    assert(bs_style=='lbs')
    result, Jtr = lbs.verts_core(pose, v, J, weights, kintree_table, want_Jtr=True, xp=chumpy)

    tr = trans.reshape((1,3))
    result = result + tr
    Jtr = Jtr + tr

    result.trans = trans
    result.f = f
    result.pose = pose
    result.v_template = v_template
    result.J = J
    result.weights = weights
    result.kintree_table = kintree_table
    result.bs_style = bs_style
    result.bs_type =bs_type
    if posedirs is not None:
       result.posedirs = posedirs
       result.v_posed = v_posed
    if shapedirs is not None:
       result.shapedirs = shapedirs
       result.betas = betas
       result.v_shaped = v_shaped
    if want_Jtr:
       result.J_transformed = Jtr
    return result

# verts_core :overloaded function inherited by lbs.verts_core
def verts_core(pose, v, J, weights, kintree_table, bs_style, want_Jtr=False, xp=chumpy):

    if xp == chumpy:
       assert(hasattr(pose, 'dterms'))
       assert(hasattr(v, 'dterms'))
       assert(hasattr(J, 'dterms'))
       assert(hasattr(weights, 'dterms'))

    assert(bs_style=='lbs')
    result = lbs.verts_core(pose, v, J, weights, kintree_table, want_Jtr, xp)

    return result

4、serialization.py

SMPL模型的序列化函数。

def save_model(model, fname): # 保存模型model到文件名fname中
    m0 = model
    # 将模型参数存入字典。各个模型参数的含义:
    # v_template:顶点模板 J:关节坐标,weight:关节权重,kintree_table:关节树,f:相机参数,bs_type:蒙皮方法LBS还是DQBS
    trainer_dict = {'v_template': np.asarray(m0.v_template),'J': np.asarray(m0.J),'weights': np.asarray(m0.weights),'kintree_table': m0.kintree_table,'f': m0.f, 'bs_type': m0.bs_type, 'posedirs': np.asarray(m0.posedirs)}    
    if hasattr(model, 'J_regressor'):
        trainer_dict['J_regressor'] = m0.J_regressor
    if hasattr(model, 'J_regressor_prior'):
        trainer_dict['J_regressor_prior'] = m0.J_regressor_prior
    if hasattr(model, 'weights_prior'):
        trainer_dict['weights_prior'] = m0.weights_prior
    if hasattr(model, 'shapedirs'):
        trainer_dict['shapedirs'] = m0.shapedirs
    if hasattr(model, 'vert_sym_idxs'):
        trainer_dict['vert_sym_idxs'] = m0.vert_sym_idxs
    if hasattr(model, 'bs_style'):
        trainer_dict['bs_style'] = model.bs_style
    else:
        trainer_dict['bs_style'] = 'lbs'
    pickle.dump(trainer_dict, open(fname, 'w'), -1) #写入文件

# 加载model时更改参数名称
def backwards_compatibility_replacements(dd):

    # replacements
    if 'default_v' in dd:
        dd['v_template'] = dd['default_v']
        del dd['default_v']
    if 'template_v' in dd:
        dd['v_template'] = dd['template_v']
        del dd['template_v']
    if 'joint_regressor' in dd:
        dd['J_regressor'] = dd['joint_regressor']
        del dd['joint_regressor']
    if 'blendshapes' in dd:
        dd['posedirs'] = dd['blendshapes']
        del dd['blendshapes']
    if 'J' not in dd:
        dd['J'] = dd['joints']
        del dd['joints']

    # defaults
    if 'bs_style' not in dd:
        dd['bs_style'] = 'lbs'


# 准备参数
def ready_arguments(fname_or_dict):

    if not isinstance(fname_or_dict, dict):
        dd = pickle.load(open(fname_or_dict))
    else:
        dd = fname_or_dict
        
    backwards_compatibility_replacements(dd)
        
    want_shapemodel = 'shapedirs' in dd
    nposeparms = dd['kintree_table'].shape[1]*3 # 这些不是pose的变量,是关节树

    if 'trans' not in dd:
        dd['trans'] = np.zeros(3) # 定义平移量为0 
    if 'pose' not in dd:
        dd['pose'] = np.zeros(nposeparms) # 定义姿势为 0
    if 'shapedirs' in dd and 'betas' not in dd:
        dd['betas'] = np.zeros(dd['shapedirs'].shape[-1])

    for s in ['v_template', 'weights', 'posedirs', 'pose', 'trans', 'shapedirs', 'betas', 'J']:
        if (s in dd) and not hasattr(dd[s], 'dterms'):
            dd[s] = ch.array(dd[s])

    if want_shapemodel: # 需要形状模型
        dd['v_shaped'] = dd['shapedirs'].dot(dd['betas'])+dd['v_template']
        v_shaped = dd['v_shaped']
        J_tmpx = MatVecMult(dd['J_regressor'], v_shaped[:,0])        
        J_tmpy = MatVecMult(dd['J_regressor'], v_shaped[:,1])        
        J_tmpz = MatVecMult(dd['J_regressor'], v_shaped[:,2])        
        dd['J'] = ch.vstack((J_tmpx, J_tmpy, J_tmpz)).T    
        dd['v_posed'] = v_shaped + dd['posedirs'].dot(posemap(dd['bs_type'])(dd['pose']))
    else:    
        dd['v_posed'] = dd['v_template'] + dd['posedirs'].dot(posemap(dd['bs_type'])(dd['pose']))
            
    return dd



def load_model(fname_or_dict):
    dd = ready_arguments(fname_or_dict)
    
    args = {
        'pose': dd['pose'],
        'v': dd['v_posed'],
        'J': dd['J'],
        'weights': dd['weights'],
        'kintree_table': dd['kintree_table'],
        'xp': ch,
        'want_Jtr': True,
        'bs_style': dd['bs_style']
    }
    
    result, Jtr = verts_core(**args)
    result = result + dd['trans'].reshape((1,3))
    result.J_transformed = Jtr + dd['trans'].reshape((1,3))

    for k, v in dd.items():
        setattr(result, k, v)
        
    return result

 




上一篇:
下一篇:

头像

guodong

14条评论

你好,问个问题,chumpy定义地的对象什么时候会自动调用dr_wrt,为什么每个自定义的对象都定义了
compute_dr_wrt

可以得到所有顶点和关节点3D坐标,只是是以臀部的那个关节点为原点的局部坐标系

生成的smpl模型不是只有表面顶点吗,内嵌的关节点如何得到呢?也可能是我理解很浅~

我结合你的文章看了一下hmr的源码,发现73行所得到的keypoints应该只直接从model文件中读取的,所以不太清楚具体计算过程。我在https://github.com/YinghaoHuang91/MuVS/blob/c106adc7bee6f14b0717088beaca7d74e20435d4/Code_TF/smpl_batch.py这里找到了由smpl参数计算3d joints的方法~谢谢您的回复,祝科研顺利~

你好,我看到了你给的计算 smpl 3d joints 坐标的链接。在读完代码后,我还是有不懂的地方,希望能咨询一下你,方便加我qq吗? qq:1546375589.

补充一下**LBS.py**,result:local axis-angle,results_global:global axis-angle(by root),个人理解,不对请指正。

您好, 我想请问一下,源码中 为什么直接更新pose和betas参数,模型的结果就会自动计算更新呢? 不好意思, 谢谢了!

请问一下,这里:ch.concatenate([(Rodrigues(pp) – ch.eye(3)).ravel() for pp in p]).ravel()
罗德里格斯公式算出来的旋转矩阵为什么还要减去一个对角矩阵?谢谢

xinqi进行回复 取消回复

电子邮件地址不会被公开。 必填项已用*标注