前言

之前介绍了使用传统算法求解BS系数的表情驱动方法，其中提到过的三种方法之一是基于网格形变迁移做的，那么这篇文章就是对《Deformation Transfer for Triangle Meshes》做表情驱动的解析。

国际惯例，参考博客：

论文原文《Deformation Transfer for Triangle Meshes》
大佬的代码

本博文实现几乎照搬大佬代码，但是大佬代码实在是太多太复杂了，搞了很多库和函数，还有很多重复的实现，所以本博客一边解析论文，一边按照大佬的代码简化实现，把几千行代码缩短到几百行。

python的预备知识

因为3D模型的网格数量一般比较大，所以论文的实现需要稀疏矩阵

稀疏矩阵相关库说明scipy.sparse
由于涉及到临近点查找，所以需要建立cKDTree，同时需要熟悉其查找函数query，专门用于在KD树中查找指定点的临近点。
从一个网格到另一个网格的形变量，涉及到大型矩阵的求解级Ax=B ，其中A和B是超大的矩阵建立的稀疏矩阵，所以求解x需要使用稀疏矩阵求解，代码使用的是LU分解

理论与实现

概述

论文主要包括两部分：

形变迁移：原始模型的表情变化量(网格形变情况)迁移到目标模型中，比如角色1从正常到笑脸的表情变化量迁移到角色2的正常表情上，使角色2变成笑脸
网格对齐：角色1和角色2需要找到网格面片的对应关系，告诉算法原模型的面片100发生变换时，目标模型应该是哪个网格面片发生对应变换。

论文第5章介绍的是网格对齐，第3、4章分别介绍面片和顶点做形变迁移；我们按照正常顺序解析，先找到源模型和目标模型的网格对应关系，再去做形变迁移。

网格对应关系

找网格面片对应关系的思想是将源模型尽量形变到目标模型的同时，要保证：

源网格形变后能够保证自身的大概结构，不能本来是猫，形变以后变成狗了；这一步要求每个网格形变矩阵尽量接近单位阵
形变过程中，具有相同顶点的网格要具有接近甚至是相同的形变矩阵，不然同一个顶点左边的网格让他往上形变，而右边的网格让他往下形变，这样就造成了这个顶点不连续的现象；这一步要求相邻网格的形变矩阵的差值为零矩阵
源模型的所有顶点形变后，目标模型中与源模型指定顶点的最临近点要满足人工指定的顶点约束；即形变后，最开始的人为指定的约束顶点对的空间位置要尽量接近；这一步要求形变后，在目标模型中查找形变后源模型指定顶点的最近点的位置要与人为指定的顶点接近。

上面三个约束就对应论文第五章介绍的三种损失函数，这样我们计算两个模型对应点的目标就是
$min w_sE_s+w_IE_I+w_CEc$
式子约束为将源模型的顶点映射到目标模型的某个对应顶点上去

网格信息的建立

为计算上述三个约束的损失函数值，我们需要对原始模型做一些额外信息的提取：

首先是提取对应源模型和目标模型各自的相邻面片索引，如与第0个面片与 $[4, 5, 6]$ 面片相邻等，按照边来计算，如果两个面片具有同样的边，那么他俩就是相邻面，是实现如下：

def compute_adjacent_by_edges(objFaces):# 每条边涉及到哪些面candidates = defaultdict(set)for i in range(objFaces.shape[0]):f0, f1, f2 = sorted(objFaces[i])candidates[(f0, f1)].add(i) # 注意i是面索引candidates[(f0, f2)].add(i)candidates[(f1, f2)].add(i)# 每个面与哪些面邻接；candidates的value存的就是共享边的面faces_adjacent = defaultdict(set)  # Face -> Facesfor faces in candidates.values():for f in faces:faces_adjacent[f].update(faces)# 按面的顺序排列所有邻接面faces_sort = []for f, adj in faces_adjacent.items():exclude_f = []for a in adj :if a != f:exclude_f.append(a)faces_sort.append([f,exclude_f])faces_sort = sorted(faces_sort, key=lambda e: e[0])# 只返回邻接面faces_adj = []for _,ff in faces_sort:faces_adj.append(ff)return faces_adj

可以看到流程大概是：先记录每条边属于那些面，然后按照每个边记录的面索引集合，提取每个面的相邻面，最后按照顺序将自己去除，按索引顺序记录他的相邻面，比如 $[0, 1, 2, 3]$ 是相邻面，那么第0个面的相邻面索引记录为 $[1, 2, 3]$ ，而第2个面的相邻面索引记录为 $[0, 1, 3]$

论文写了三角面并不能准确的表现出网格的变换，需要为每个面增加额外的一个顶点，代表的方向为当前面的法线方向，所以每个面需要被重新调整一下，增加面的法线顶点

def compute_face_norms(verts,faces):fnorms = []# 计算每个面片的三组方向for f in faces:v1 = verts[f[0]]v2 = verts[f[1]]v3 = verts[f[2]]a = v2 - v1b = v3 - v1tmp = np.cross(a, b)c = tmp.T / np.linalg.norm(tmp)fnorms.append([a,b,c])fnorms = np.array(fnorms) # 更新顶点，添加第四个顶点v4 = v1 + fnorms[...,-1]    new_verts = np.concatenate((verts,v4),axis=0)# 更新面片顶点索引v4_indices = np.arange(verts.shape[0],verts.shape[0]+v4.shape[0])new_faces = np.concatenate((faces,v4_indices.reshape((-1,1))),axis=1)return np.transpose(fnorms, (0, 2, 1)),new_verts,new_faces

找最近点的时候，我们不仅要看这个点的位置是否与目标点接近，还得考虑他的邻接面的法线方向是否一致，比如我找到模型两个顶点接近，但是一个发现向上一个法线向下，他俩也不能作为对应点，但是判断顶点附近面的法线很麻烦，所以为当前顶点利用临近面法线的平均计算了顶点法线。

'''
计算顶点法线：面法线的平均
'''
def compute_vert_norms(faceNorms,faces):vf = defaultdict(set)# 每个顶点与哪些面片有关for n, (f0, f1, f2) in enumerate(faces[:, :3]):vf[f0].add(n)vf[f1].add(n)vf[f2].add(n)# 计算每个顶点的法线，用面法线的平均vn = np.zeros((len(vf),3))for i in range(len(vf)):vn[i] = np.mean(faceNorms[list(vf[i])],axis=0)vn[i] = vn[i] / np.linalg.norm(vn[i])return vn

求解稀疏矩阵的前提是要建立稀疏矩阵，所以需要将每个面片的信息保留在稀疏矩阵中
```
'''
创建稀疏矩阵
'''
row = np.array([0, 1, 2] * 4)
def expand(f, inv, size):i0, i1, i2, i3 = fcol = np.array([i0, i0, i0, i1, i1, i1, i2, i2, i2, i3, i3, i3])data = np.concatenate([-inv.sum(axis=0), *inv])return sparse.coo_matrix((data, (row, col)), shape=(3, size), dtype=float)
def construct(faces, invVs, size):assert len(faces) == len(invVs)return sparse.vstack([expand(f, inv, size) for f, inv in zip(faces, invVs)], dtype=float)
```
上面函数的意思是模型所有面片faces在invVs中记录了相关信息，需要把每个面片分别构建稀疏矩阵，然后组合起来；具体每个面片存储的信息就是论文公式(3)中的有面片边向量和法线向量组成的矩阵，然后按照公式(4)求逆，即
$V = [v_2-v_1,v_3-v_1,v_4-v_1]$
因为v4代表的顶点，由 $v 1$ 加上法线方向组成，所以 $v 4 - v 1$ 就是法线方向。

因为上面仅仅是计算了 $V^{-1}$ ，而我们的目标是计算损失 $B-V^{-1}A$ ，其中B是目标比如约束1的目标是单位阵，约束2的目标是零矩阵，所以还得有一个计算损失的函数

def apply_markers(A, b, target, markers):"""Ei的形变量接近单位阵Es的形变量是邻接面的形变量接近:param A: Matrix (NxM):param b: Result vector (Nx3):param target: Target mesh:param markers: Marker (Qx2) with first column the source indices and the second the target indices.:return: Matrix (Nx(M-Q)), result vector (Nx3)"""assert markers.ndim == 2 and markers.shape[1] == 2invmarker = np.setdiff1d(np.arange(A.shape[1]), markers[:, 0])# 不在marker中的点zb = b - A[:, markers.T[0]] * target[markers.T[1]] # 使得形变量接近breturn A[:, invmarker].tocsc(), zb

上述函数的功能就是计算指定顶点的形变矩阵与目标矩阵的损失值。也即变换矩阵应该与单位阵或者零矩阵差多少。

构建保证结构的损失

为了保证变形后，猫还是猫，而不是变成狗，我们需要变形尽量是单位阵

def construct_identity_cost(faces,invVs,verts):AEi = construct(faces,invVs,verts.shape[0])AEi.eliminate_zeros()Bi = np.tile(np.identity(3, dtype=float), (len(faces), 1))assert AEi.shape[0] == Bi.shape[0]return AEi.tocsr(), Bi

上述返回的就是论文中的 $V^{-1}$ 和单位阵，然后调用时候，计算损失的方法为：

Ei,Bi = construct_identity_cost(source_f,invVs,source_v)
EiSelect,EiLoss = apply_markers(Ei,Bi,target_v,marks)

这样就能返回指定顶点的 $V^{-1}$ 和与保持原形状目标(单位阵)的损失。

构建保证平滑性的损失

上面提到过，相邻面片的形变量最好接近甚至是一致，不然左边面片让顶点往下偏移，右边面片让顶点往上偏移，这样就断裂不连续了，如果这俩面片的形变矩阵相同，那么他们都会往同一个位置偏移了，毕竟面片的形变本质上就是把形变矩阵用到面片的每个顶点上。

所以我们要对相邻的面片对构建他们的形变矩阵差值，使其目标接近零矩阵

def construct_smoothness_cost(faces,invVs,verts,adjacent):# 预先计算并保存lhs = []rhs = []face_idx = 0for f,inv in zip(faces,invVs):for adjIndex in adjacent[face_idx]:            lhs.append(expand(f,inv,verts.shape[0]).tocsc())rhs.append(expand(faces[adjIndex],invVs[adjIndex],verts.shape[0]).tocsc())face_idx = face_idx + 1AEs = sparse.vstack(lhs) - sparse.vstack(rhs)AEs.eliminate_zeros()    count_adjacent = sum(len(a) for a in adjacent)Bs = np.zeros((count_adjacent * 3, 3))assert AEs.shape[0] == Bs.shape[0]return AEs, Bs

上面有个lhs和rhs就是对当前邻接面片的遍历，构建成对的信息，计算他俩的差值，因此后续出现了 $l h s - r h s$ ，而目标Bs就是零矩阵。

【注】当面片数量比较大的时候，这一步的for循环和vstack大概率会把内存怼满，就无法进行后续的计算了，所以可以采用暂存的方法，把 $A E s$ 存起来，后面复用

sparse.save_npz("cat-lion.npz",AEs)
AEs = sparse.load_npz("head.npz")

接下来与结构损失的计算同理

Es,Bs = construct_smoothness_cost(source_f,invVs,source_v,source_face_adj)
EsSelect,EsLoss = apply_markers(Es,Bs,target_v,marks)

这样平滑损失也构建完毕

构建距离损失

这个稍微麻烦点，因为涉及到两个模型的最近点索引得查找，所以应用到了上面说过的cKDTree去快速查询临近点

def get_closet_points(kd_tree,verts,verts_norms,target_normal,max_angle=np.radians(90),ks=200):'''寻找最近点，同时满足法线角度差小于max_angle,寻找最近点的数量不大于ks'''assert len(verts) == len(verts_norms),"source verts and norms is not same length"closest_points = []dists,indices = kd_tree.query(verts,min(len(target_normal),ks)) # 找到在最近点的距离和索引for i,(dist,ind) in enumerate(zip(dists,indices)):angles = np.arccos(np.dot(target_normal[ind],verts_norms[i]))if((np.abs(angles)<max_angle).any()):cind = ind[np.abs(angles)<max_angle][0]closest_points.append([i,cind])return np.array(closest_points)

上面的函数就是利用query查找临近点，由于要计算法线夹角约束，可能最近的点并不一定满足法线约束，所以提取了多个临近点，然后选取满足法线约束的第一个顶点即为满足要求的临近点。

那么这一点就是要求查找的最近点与人工指定的对应点的距离最短，那么：

先查点：

closest_points = get_closet_points(kd_tree_target,vertices_clipped,compute_vert_norms(vertices_fnorm_copy[...,-1],source_org_f),target_vn)

再去计算距离误差：

 AEc = AEc[closest_points[:,0]]Bc = target_v[closest_points[:,1]]mAEc,mBc = apply_markers(AEc,Bc,target_v,marks)

其中marks就是人工挑选的对应点，AEc为单位阵，target_v为目标模型的点，Bc为计算得到的最近点的位置

最小化损失函数

上面提到损失函数就是三个约束的损失值得和，那么我们可以把它堆叠起来：

Astack = [EiSelect * Wi, EsSelect * Ws]
Bstack = [EiLoss * Wi, EsLoss * Ws]
Astack.append(mAEc * Wc[iteration])
Bstack.append(mBc * Wc[iteration])

但是按照论文说明，第一次迭代没必要把距离损失接进来，因为刚开始如果两个模型大小相差一万倍，就没必要去算这个最近点了，所以可以先利用其它两个损失，将源模型形变到与目标模型差不多大小和方向，然后第二次迭代开始再加入顶点距离损失。

接下来怎么利用这个方程找到两个模型得对应面片索引，按照代码中说的：

直接进行LU分解:

 A = sparse.vstack(Astack,format="csc")A.eliminate_zeros()b = np.concatenate(Bstack)LU = sparse_linalg.splu((A.T @ A).tocsc())x = LU.solve(A.T @ b)

从上式可以看出，本来我们要解 $A x = b$ 的，但是代码直接转换成了 $A^TAx= A^Tb$ ，那是因为 $A^TA$ 可以把左矩阵变成方阵，便于求解。

随后，我们就可以在迭代过程中，利用损失进一步重构源模型的顶点：

def revert_marks(A,x,target_v,markers,out):if out is None:out = np.zeros((A.shape[1] + len(markers), 3))else:assert out.shape == (A.shape[1] + len(markers), 3)invmarker = np.setdiff1d(np.arange(out.shape[0]),markers[:,0])out[invmarker] = xout[markers[:,0]] = target_v[markers[:,1]]return out
revert_marks(A,x,target_v,marks,vertices_copy) 
vertices_fnorm_copy,vertices_copy,_ = compute_face_norms(vertices_copy[:source_org_v.shape[0]],source_org_f)

以上函数的意思就是，按照论文公式(9)(10)的描述

求解的x就是未知顶点的形变位置

所以把源模型中非人为指定的顶点更新为x，而已指定的顶点位置更新为目标模型中指定顶点的位置。

所以按照论文中的迭代求解方法，代码如下：

vertices_copy = source_v.copy()
vertices_fnorm_copy = source_f_norms.copy()
for iteration in range(iterations):Astack = [EiSelect * Wi, EsSelect * Ws]Bstack = [EiLoss * Wi, EsLoss * Ws]if(iteration>0 and Wc[iteration]!=0):AEc = sparse.identity(source_v.shape[0], dtype=float, format="csc")[:source_org_v.shape[0]]vertices_clipped = vertices_copy[:source_org_v.shape[0]]# 从source中获取到与target顶点的最近点，且满足法线夹角约束closest_points = get_closet_points(kd_tree_target,vertices_clipped,compute_vert_norms(vertices_fnorm_copy[...,-1],source_org_f),target_vn)AEc = AEc[closest_points[:,0]]Bc = target_v[closest_points[:,1]]mAEc,mBc = apply_markers(AEc,Bc,target_v,marks)Astack.append(mAEc * Wc[iteration])Bstack.append(mBc * Wc[iteration])A = sparse.vstack(Astack,format="csc")A.eliminate_zeros()b = np.concatenate(Bstack)LU = sparse_linalg.splu((A.T @ A).tocsc())x = LU.solve(A.T @ b)# 重建顶点vertices_copy = revert_marks(A,x,target_v,marks,vertices_copy) vertices_fnorm_copy,vertices_copy,_ = compute_face_norms(vertices_copy[:source_org_v.shape[0]],source_org_f)

意思就是第一次迭代不适用距离损失，后续的迭代再把距离损失加入，同时不断更新模型顶点

查找对应面片

多次迭代优化以后，就可以按照查找最近点的方法找到两个模型的对应面片了；

可以用最大的三角形的边长作为寻找半径，以每个三角面的中心点坐标为寻找目标：

计算三角形最大边长：

def max_triangle_length(v,faces):fnorms = []max_len = 0.0# 计算每个面片的三组方向for f in faces:v1 = v[f[0]]v2 = v[f[1]]v3 = v[f[2]]a = v2 - v1b = v3 - v1if(max_len<max(max_len,np.linalg.norm(a),np.linalg.norm(b))):max_len = max(max_len,np.linalg.norm(a),np.linalg.norm(b))return max_len

计算每个三角形的中心：

def get_centroids(v,f):return v[f[:,:3]].mean(axis=1)

按照中心和半径查找最近点，同时注意法线约束，参考的大佬的开源代码这里写的有问题，所以稍作修正了：

def get_closet_triangles(vn1,vn2,center1,center2,radius):max_angle = np.radians(90)k = 500triangles = set()kkd_tree = cKDTree(center2)dists,indicies = kkd_tree.query(center1,min(center2.shape[0],k))for index_source,(dist,ind) in enumerate(zip(dists,indicies)):angles = np.arccos(np.dot(vn2[ind],vn1[index_source]))if((angles<max_angle).any()):index_target = ind[angles<max_angle][0]triangles.add((index_source,index_target))return triangles

为了找到更多的对应面，我们需要先从源模型找目标模型对应面，再反过来从目标模型找源模型对应面：

def match_triangles(v1,v2,f1,f2,vn1,vn2,factor=2):source_f_center = get_centroids(v1,f1)target_f_center = get_centroids(v2,f2)   radius = max(max_triangle_length(v1,f1),max_triangle_length(v2,f2))*factor    triangles = get_closet_triangles(vn1,vn2,source_f_center,target_f_center,radius)tmp_triangles = get_closet_triangles(vn2,vn1,target_f_center,source_f_center,radius)triangles.update((t[1],t[0]) for t in tmp_triangles)return list(triangles)

这样就找到两个模型的对应面片索引了。

形变迁移

形变迁移就是将源模型从形状1变换到形状2时候，每个面片的变化矩阵被应用到另一个模型上，此模型与模型1具有相似的语义(如都是正常表情),且知道模型间的面片对应关系。

计算面法线：

source_f_norms,source_v,source_f = compute_face_norms(source_org_v,source_org_f)
target_f_norms,target_v,target_f = compute_face_norms(target_org_v,target_org_f)
pose_f_norms,pose_v,pose_f = compute_face_norms(pose_org_v,pose_org_f)

还是先对目标模型构建稀疏矩阵：

Am = construct(target_f[mapping[:,1]],inv_target_span[mapping[:,1]],len(target_v))

随后构建源模型到形变后模型的对应顶点变化量，比如从一个模型正常表情到笑脸的网格变化量：

s = (pose_f_norms @ np.linalg.inv(source_f_norms)).transpose(0,2,1)
Bm = np.concatenate(s[mapping[:,0]])

考虑到有的面片可能并没有映射关系，所以需要一个平滑损失，去保证形变时不会出现断裂，和建立网格对应关系时候的那个约束一样：

# 没有对应面的部分需要平滑处理
adjacent = compute_adjacent_by_edges(target_f[:,:-1])
inv_target_span = np.linalg.inv(target_f_norms)
missing = np.setdiff1d(np.arange(len(target_f)),np.unique(mapping[:,1]))
count_adjacent = sum(len(adjacent[m]) for m in missing)
if(count_adjacent==0):AEs = sparse.csc_matrix((0,len(target_v)),dtype=float)Bs = np.zeros((0,3))
else:size = len(target_v)lhs = []rhs = []face_idx = 0for f,inv in zip(target_f,inv_target_span):for adjIndex in adjacent[face_idx]:            lhs.append(expand(f,inv,target_v.shape[0]).tocsc())rhs.append(expand(target_f[adjIndex],inv_target_span[adjIndex],target_v.shape[0]).tocsc())face_idx = face_idx + 1AEs = sparse.vstack(lhs) - sparse.vstack(rhs)Bs = np.zeros((AEs.shape[0],3))

接下来构建求解目标：

Wm = 1.0 
Ws = 1.0Astack = [Am*Wm, AEs*Ws]
Bstack = [Bm*Wm, Bs*Ws]
A = sparse.vstack(Astack,format='csc')
A.eliminate_zeros()
b = np.concatenate(Bstack)

求解:

LU = sparse.linalg.splu((A.T @ A).tocsc())
x = LU.solve(A.T @ b)
result = x[:source_v.shape[0]]

其中x就是目标模型做对应形变后的结果啦，保存为obj即可：

def saveObj(verts,fvidxs,name):f = open(name,'w')color = np.random.rand(3)for v in verts:f.write("v {0} {1} {2} {3} {4} {5} \n".format(v[0],v[1],v[2],color[0],color[1],color[2]))for fv in fvidxs:f.write("f ")for ii in range(len(fv)):f.write("{0} ".format(fv[ii]+1))f.write("\n")f.close()
saveObj(result,target_org_f,"result.obj")

可视化结果如下：

对于人脸：