python Deformation Transfer for Triangle Meshes-CSDN博客

本文链接：https://blog.csdn.net/hongmaodaxia/article/details/103800736

这个主要是给自己看的，这个脚本是自己第一次思考如何构建如此复杂稀疏矩阵的方法，并进行加速，里面有一些思考的要素，值得以后记住。

@20-5-20
补充之前缺少的代码【在最后】

算法要点：稀疏矩阵是 $(nF*3*3)\times(3*(nV+NF))$
构造的时候把握每个面片的构造，应该是 $A(3*3)=B(3*4)\times C(4*3)$ . 其中B是个稀疏矩阵，本意就是构造 $A=((v2-v1)^T, (v3-v1)^T, (v4-v1)^T)$ ，也可以看出来其实对于每个面片其实最后组成的系数矩阵每行也只有四个元素起作用。 $C=(v1^T, v2^T, v3^T, v4^T)$ 。

$B=\begin{bmatrix} -1 & 1 & & \\ -1 & & 1 & \\ -1 & & & 1 \end{bmatrix}$

注意 v4是每个面有一个，共nF 个，因此优化的顶点数为 nV+nF 个。
然后就需要构造最终的大矩阵，设法把B提前乘进去，这样变量 x 就只包括顶点位置了。然后求解类似 Ax=b 这种的方程就行（和上面那个A不一样）。

事后感觉, 这里似乎就是为了避免求法向量, 直接把v4 当作变量优化出来.

S: V1->V2
T: V3->V4

def run():
    fobjs=["../tempdir/scans_hd_smooth2/36/1.obj", 
        "../tempdir/scans_hd_smooth2/36/15.obj",
        "../tempdir/scans_hd_smooth/36/1.obj"
            ]
    #nV*3, nF*3
    v1, F = get_obj_v_f(fobjs[0])
    v2 = get_obj_v(fobjs[1])
    v3 = get_obj_v(fobjs[2])
    
    f, v1, v2, v3 = [ np2th(e) for e in [F.astype(np.int64), v1, v2, v3]  ] 

    tm_bg= time.time()

    V1 = get_init_V_(v1, f)
    V2 = get_init_V_(v2, f)
    V3 = get_init_V_(v3, f)

    # V1->V2
    # 可能出现奇异值的情况
    Q=V1.pinverse().bmm( V2 ) # 果然不能死抄公式,要与自己的分析一致
    
    # 注意 待求解的是 3*(nV+nF)
    V3_inv = V3.pinverse()

    nF = f.size(0)
    nV = v1.size(0)
    
    tri=[]
    s_h = nF*9
    s_w = (nF+nV) * 3
    sparse_info=get_sparse_info( nF*3*3*4 )
    cM=torch.zeros((3, 4))
    if nptype==np.float64: cM=cM.double()
    cM[:, 0] = -1
    cM[:, 1:] = torch.eye(3)

    # A_full : (nF*9) * ((nV+nF)*3)
    _As = V3_inv.bmm(cM.unsqueeze(0).expand( nF, 3, 4 ).clone()).numpy() # 其实这个可以用einsum 

    sparse_info.i[:] = np.arange(nF*3*3).repeat(4)

    sparse_info.d[:] = _As.reshape(nF, 3, 1, 4).repeat(3, axis=2).reshape(-1)

    # 上面的比下面的快 100倍多, 简直不敢想象
    # 每行有4个, 作用nF, 9 行
    #j1 = F.reshape(nF,3, 1 ).repeat(1, 1, 9).transpose(0, 2, 1).reshape(nF, 3, 3, 3) * 3
    j1 = F.reshape(nF, 1, 3 ).repeat(9, axis=1).reshape(nF, 3, 3, 3) * 3
    j1[:, :, 1] += 1
    j1[:, :, 2] += 2
    j1 = j1.reshape(-1, 3)
    j2 = (np.arange(nF).astype(np.int32) + nV).repeat(9).reshape(nF, 3, 3) *3
    j2[:, :, 1] += 1
    j2[:, :, 2] += 2
    j2 = j2.reshape(-1, 1)
    sparse_info.j[:] = np.hstack( (j1, j2) ).reshape(-1)

    # 写的时候主要看下面到底是哪一列需要重复, 然后就 repeat

    #for i in range(nF):
    #    A = _As[i]
    #    for j in range(3):
    #        for k in range(3):
    #            bid = i*9+j*3+k
    #            for kk in range(4):
    #                sparse_info.i[bid*4+kk]=bid
    #                sparse_info.d[bid*4+kk]=A[j, kk]
    #                if kk<3:
    #                    sparse_info.j[bid*4+kk]=f[i, kk]*3+k
    #                else:
    #                    sparse_info.j[bid*4+kk]=(nV+i)*3+k

    p("time use:", time.time()-tm_bg)
    s = sparse_info
    
    A = scipy.sparse.coo_matrix(( s.d, (s.i, s.j)), shape=(s_h, s_w)).tocsr()

    b = Q.view(-1).numpy()

    AtA=A.T.dot(A)
    Atb=A.T.dot(b)
    x = spsolve(AtA, Atb)
    
    p("time use:", time.time()-tm_bg)

    v =x.reshape(-1, 3)[:nV]
    save_obj("../tempdir/ck.obj", v, F)
    
def get_sparse_info(n):
    class info:
        i=np.zeros(n , dtype=np.int32)
        j=np.zeros(n , dtype=np.int32)
        d=np.zeros(n, dtype=nptype)
    return info

# 论文中需要的V
def get_init_V_(v, f, isnp=True):
    if nptype==np.float64: v= v.double()
    vs = v[f.view(-1)].view(-1, 3, 3)
    
    a = vs[:, 0]
    b = vs[:, 1]
    c = vs[:, 2]
    
    v2_1 = b-a
    v3_1 = c-a
    v4_1 = thF.normalize(v2_1.cross(v3_1), eps=1e-20, dim=1)
    V = torch.cat( (v2_1, v3_1, v4_1), 1 )
    # 注意这里 我是按照行的和论文中不一样,但是这样求解就免去转置了

    return V.view(-1, 3, 3) 

def kronecker(A, B):
    return torch.einsum("ab,cd->acbd", A, B).view(A.size(0)*B.size(0),  A.size(1)*B.size(1))

之前缺少的代码

读取obj 有关的在下面找到
https://github.com/changhongjian/SimRender/blob/master/lib/comm/geometry.py

其他的

def np2th(e): return torch.from_numpy(e)
def np2th_gpu(e): return torch.from_numpy(e).cuda()
def th2np(e): return e.detach().cpu().numpy()

# 其他就是Import了，缺什么自己补吧
import torch
import numpy as np
import time
import scipy
import torch.nn.functional as thF
from scipy.sparse.linalg import spsolve