【网格形变】Laplacian形变的python实现

import numpy as np
class Laplacian():
    def __init__(self, mesh):
        self.mesh=mesh
        self.v = mesh.vertex # numpy array with shape of [n_v, 3], n_v is the number of  vertices
        self.f = mesh.face # numpy array with shape of [n_f, 3], n_f is the number of vertices
        self.nv = len(self.v)
        self.matLap=self.getMatLap()
        print('已获取Laplacian矩阵')

        self.posLap=self.getPosLap()
        print('已获取Laplacian坐标')

    def useAnchor(self,indexes,ps):#使用锚点对网格进行形变
        weight=0.5  # 0.5~5 #设置锚点在形变中的权重
        import copy
        matLap2 = copy.deepcopy(self.matLap)#复制list对象
        posLap2 = copy.deepcopy(self.posLap)

        for index in indexes:
            vector=self.initMat(1,self.nv)[0]
            vector[index]=weight
            matLap2.append(vector)
        for pos in ps:
            for k in range(len(pos)):
                pos[k]=pos[k]*weight
            posLap2.append(pos)
        print('开始求解超定齐次线性方程')
        return self.solveLinearEquation(matLap2,posLap2)
    def transform(self,indexes,ps):
        v2=self.useAnchor(indexes,ps)
        self.mesh.vertex=v2
        self.mesh.updateVertex()

        self.v = v2
        self.posLap=self.getPosLap()
    def solveLinearEquation(self,A,B):
        from numpy.linalg import lstsq # 解超定方程
        a=np.mat(A)#A=[[2, 3], [1, 3], [1, 1]]
        b=np.mat(B)#B=[5, 4, 2]
        result=lstsq(a, b, rcond=None)#ax=b #a'ax=a'b
        return result[0].tolist()
    def getPosLap(self):#拉普拉斯坐标 (n,n)*(n,3)
        p=np.array(self.v)
        l=np.array(self.matLap)
        posLap=np.dot(l,p)
        return posLap.tolist()
    @staticmethod
    def initMat(i,j):
        m=[]
        for k1 in range(i):
            m.append([])
            for k2 in range(j):
                m[k1].append(0)
        return m
    def getMatLap(self):#计算均值拉普拉斯矩阵
        #创建0矩阵
        matLap=self.initMat(self.nv,self.nv)

        #根据self.f给矩阵赋值
        for face0 in self.f:
            if len(face0)==3:
                a=face0[0]-1
                b=face0[1]-1
                c=face0[2]-1
            else:
                a=face0[0]-1
                b=face0[3]-1
                c=face0[6]-1
            #matLap[b][a]=matLap[a][b]=matLap[a][b]+1
            #matLap[c][a]=matLap[a][c]=matLap[a][c]+1
            #matLap[c][b]=matLap[b][c]=matLap[b][c]+1
            matLap[b][a]=matLap[a][b]=1
            matLap[c][a]=matLap[a][c]=1
            matLap[c][b]=matLap[b][c]=1

        for i in range(self.nv):
            k=0
            for j in range(self.nv):
                k=k+matLap[i][j]

            for j in range(self.nv):
                if not matLap[i][j]==0:
                    matLap[i][j]=-1/k
            matLap[i][i]=1
        return matLap
class Test():
      def __init__(self, v,f):
        self.vertex=v
        self.face=f
if __name__ == "__main__":
    m0 = Test([
                  [1,2,3],
                  [1,2,3],
                  [1,2,3],
                  [1,2,3]
              ],
              [
                  [0,1,3],
                  [1,2,3]
              ]
          )

    l=Laplacian(m0)
    v2=l.useAnchor(
                [1],
                [
                    [1.5,2.5,3.5]
                ])
    print(v2)