[Python知识库]拉普拉斯变形的原理解析和python代码

背景

拉普拉斯变形是图形学处理mesh的常用方法，它假定mesh的顶点，在变化前后，顶点的拉普拉斯距离应该是一致的。
最常见的拉普拉斯矩阵的定义如下：
$$\begin{gathered} L=D?A \\ =D(I?D^{?1}A) \end{gathered}$$
D是每个顶点的度，A是邻接矩阵。
假设有变形前的顶点为V，有L*V，将其拆解，可以发现就是求每个顶点$i$的顶点位置减去相邻的顶点位置*(1/$d_i$)。
$$L{V}_i=V_i?\sum _{j\in N_i}\frac{1}{d_i}{V}_j$$
j是顶点i的邻接顶点索引。
因此，拉谱拉斯矩阵中保存了顶点的局部信息。顶点周围的关系通过这种方式记录下来。
另外，通常我们使用归一化的拉普拉斯矩阵，即$I?D^{?1}A$，每一行之和为1.

拉普拉斯网格变形原理

为了简单描述，本文全部以2d数据作为实验对象。
有：

• 变形前的顶点$V\in R^{n\times 2}$
• 锚点A：描述了几个顶点的新位置。

求：

• 变形后的顶点位置$V^{\prime }$

因此，拉普拉斯变形问题可以看成，已知网格形状，并指定一些顶点落在新的位置上，求其他未知顶点的位置，并要求所有顶点在变化前后的拉普拉斯距离一致。因此有两个优化目标：
$${\mathrm{argmin}}_V^{\prime }||L{V}^{\prime }?LV||+||V_{指定}^{\prime }?A||$$
将LV看成b1，A看成b2，则化简为
$${\mathrm{argmin}}_V^{\prime }||L{V}^{\prime }?b1||+||V_{指定}^{\prime }?b2||$$
形式为AX-b的优化问题，即最小二乘法。最小二乘就有很多解法了。
继续化简，因为两个式子的优化对象都是V’，因此我们可以联立方程：
∣ ∣ { L ∣ B } V ′ ? { L V ∣ A } ∣ ∣ ||\{L | B\} V' - \{LV | A\} || ∣∣{L∣B}V′?{LV∣A}∣∣
|是行拼接符号，即联立方程。B是一个特殊矩阵，他为了实现锚点的已知条件，即$V^{\prime }=A$，因此B的每一行只有一个1,1的位置对应一个未知数，该未知数是锚点对应的位置。

Python实战

假设我们有一组2d点集，如蓝色的线，anchor点（黄色的点）分别对应了蓝色线在变形后的起始位置和终止位置。

求其他点的位置。

import numpy as np
from scipy.sparse import coo_matrix, block_diag
from scipy.sparse.linalg import lsqr
from collections import defaultdict

# https://github.com/SecretMG/Laplacian-Mesh-Deformation
class LaplacianDeformation:
    def __init__(self, vps, faces) -> None:
        self.vps = vps
        self.faces = faces
        self.adj_info = self.get_adj_info()
        self.mode = 'mean'
        self.dim = vps.shape[1]
    
    def get_adj_info(self):
        adj_dic = defaultdict(set)
        for idx, face in enumerate(self.faces):
            for f in face:
                adj_dic[idx].add(f)
        return adj_dic


    def get_Ls_matrix(self, anchor_ids, anchor_weight=None):
        k = anchor_ids.shape[0]
        n = self.vps.shape[0]
        # 生成稀疏矩阵，不建议使用np.zeros依次赋值
        data = []
        I = [] # 行序号
        J = []  # 列序号
        for i in range(n):
            neighbors = [v for v in self.adj_info[i]]  # ids
            degree = len(neighbors)
            # data += [degree] + [-1] * degree  # D-A
            data += [1] + [-1 / degree] * degree  # (I - D^{-1}A)
            I += [i] * (degree + 1)
            J += [i] + neighbors
        # 为了anchor增加的方程组，组成了超定方程 
        for i in range(k):
            if anchor_weight is None:
                data += [1]
            else:
                data += [anchor_weight[i]]
            I += [n+i]
            J += [anchor_ids[i]] 
        # [coco_matrix] https://blog.csdn.net/kittyzc/article/details/126077002
        Ls = coo_matrix((data, (I, J)), shape=(n+k, n)).todense()
        return Ls

    def solve(self, anchors, anchor_ids):
        k = anchor_ids.shape[0]
        n = self.vps.shape[0]
        Ls = self.get_Ls_matrix(anchor_ids)
        print(Ls)
        delta = Ls.dot(self.vps) # n+k, dim
        # print(delta.shape);exit()
        # 为后k行的原始顶点坐标修改为anchor的坐标（即修改为人为指定的已知条件）
        for i in range(k):
            delta[n+i] = anchors[i]  
            # update mesh vertices with least-squares solution
        res = np.zeros((n, self.dim))
        delta = np.array(delta)
        for i in range(self.dim):
            # res[:, i] = lsqr(Ls, delta[:, i])[0] # n,
            res[:, i] = np.linalg.pinv(Ls).dot(delta[:, i])
        return res
   
if __name__ == '__main__':
    import matplotlib.pyplot as plt
    x = np.linspace(0, 40, 100)
    y = np.sin(0.5*x)
    anchor_ids = np.array([100, 1]) - 1
    anchors = np.array([x[-1], y[-1]+10, x[0], y[0]]).reshape(-1, 2)
    vps = np.stack((x, y), axis=1) # n, 2
    faces = []
    for idx in range(len(x)):
        if idx == 0:
            faces.append((idx+1,))
        elif idx == (len(x)-1):
            faces.append((idx-1,))           
        else:
            faces.append((idx-1, idx+1))

    model = LaplacianDeformation(vps, faces)
    new_pnts = model.solve(anchors, anchor_ids)
    plt.scatter(x, y)
    plt.scatter(new_pnts[:, 0], new_pnts[:, 1])
    plt.show()
    

代码借鉴了一些github代码，已经在注释给出。
如果有想应用在3d点上，则需要一些修改。

带有旋转的拉普拉斯变形

上一节展示的拉普拉斯矩阵有一个问题，当顶点变形涉及到旋转变化或者近似旋转的变化，即便形状一样，但是拉普拉斯距离已经发生了变化。

因此，在Laplacian Surface Editing原文论文中还提到了一种更准确的拉普拉斯变形方法。
首先假设每个顶点的旋转都是比较小，将旋转矩阵近似为一个简单形式，

然后要求旋转前后的拉普拉斯距离相同，则约束变成：
$${\mathrm{argmax}}_V^{\prime }||L{V}^{\prime }?TLV||+||V_{指定}^{\prime }?A||$$
T是旋转矩阵。为了求T，设立T的求解方程为

通过求解发现，T是V’是线性关系，可以用V’的线性方程表示，接着带入约束条件，使用最小二乘法求解即可。
得到的结果如下图所示，可以发现整体的曲线更加保持了原始的正弦曲线的弧度。

# https://github.com/luost26/laplacian-surface-editing
class LaplacianDeformationWithT(LaplacianDeformation):
    def __init__(self, vps, faces) -> None:
        super().__init__(vps, faces)
        self.extra_A = None

    def get_Ls_matrix(self):
        n = self.vps.shape[0]
        # 生成稀疏矩阵，不建议使用np.zeros依次赋值
        data = []
        I = [] # 行序号
        J = []  # 列序号
        for i in range(n):
            neighbors = [v for v in self.adj_info[i]]  # ids
            degree = len(neighbors)
            # data += [degree] + [-1] * degree  # D-A
            data += [1] + [-1 / degree] * degree  # (I - D^{-1}A)
            I += [i] * (degree + 1)
            J += [i] + neighbors
        Ls = coo_matrix((data, (I, J)), shape=(n, n)).todense()
        return Ls

    def solve(self, anchors, anchor_ids, anchor_weight=1):
        k = anchor_ids.shape[0]
        n = self.vps.shape[0]
        Ls = self.get_Ls_matrix() # n, n
        delta = Ls.dot(self.vps) # n, dim
        LS = np.zeros([self.dim*n, self.dim*n])
        for idx in range(self.dim):
            LS[idx*n:(idx+1)*n, idx*n:(idx+1)*n] = (-1) * Ls
        # [
        #   [s, -w],
        #   [w, s],
        # ]
        for idx in range(n):
            neighbors = [v for v in self.adj_info[idx]]  # ids
            ring = np.array([idx] + neighbors)
            V_ring = self.vps[ring]
            n_ring = V_ring.shape[0]
            A = np.zeros((n_ring * self.dim, 4))
            for j in range(n_ring):
                A[j]          = [V_ring[j,0], -V_ring[j,1], 1, 0]
                A[j+n_ring]   = [V_ring[j,1], V_ring[j,0],  0, 1]
                # A[j+2*n_ring] = [V_ring[j,2],  V_ring[j,1], -V_ring[j, 0],           0 , 0, 0, 1]
            inv_A = np.linalg.pinv(A)  # 4, n_ring*2
            s = inv_A[0]
            w = inv_A[1]
            T_delta = np.stack(
                [
                    delta[idx, 0] * s - delta[idx, 1] * w,
                    delta[idx, 0] * w + delta[idx, 1] * s
                ], axis=0
            ) # 2, 2*n_ring
            LS[idx, np.hstack([ring, ring+n,])] += T_delta[0]
            LS[idx+n, np.hstack([ring, ring+n])] += T_delta[1]

        constraint_coef = []
        constraint_b = []
        for idx in range(k):
            vps_idx = anchor_ids[idx]
            tmp_coeff_x = np.zeros((self.dim * n))
            tmp_coeff_x[vps_idx] = anchor_weight
            tmp_coeff_y = np.zeros((self.dim * n))
            tmp_coeff_y[vps_idx+n] = anchor_weight
            
            constraint_coef.append(tmp_coeff_x)
            constraint_coef.append(tmp_coeff_y)
            constraint_b.append(anchor_weight * anchors[idx, 0])
            constraint_b.append(anchor_weight * anchors[idx, 1])
        constraint_coef = np.matrix(constraint_coef)
        constraint_b = np.array(constraint_b)
        A = np.vstack([LS, constraint_coef])
        b = np.hstack([np.zeros(self.dim * n), constraint_b])
        spA = coo_matrix(A).todense()
        # 求解的时候保证shape为二维
        V_prime = lsqr(spA, b)[0]
        # V_prime = np.linalg.lstsq(spA, b.reshape(-1, 1), rcond=None)[0]
        # V_prime = np.linalg.pinv(spA).dot(b)
        new_pnts = []
        for idx in range(n):
            new_pnts.append(list(V_prime[idx + i*n] for i in range(self.dim)))
        new_pnts = np.array(new_pnts)
        return new_pnts

参考

一个国外大学课件
简单拉普拉斯网格变形-Siggraph 2004
用数学编辑3D模型（一）- Mesh Deformation with Laplacian Coordinates
github C++
github Python
图神经网络和图表征学习笔记（三）拉普拉斯矩阵与谱方法
Laplacian Deformation on 2D/3D Mesh Editing 有详细的推导

拉普拉斯变形后面有一些改进方法，论文地址为

• Laplacian Mesh Optimization