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