devbox
菜鸟追梦旅行
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

拉普拉斯变形的原理解析和python代码_拉普拉斯形变

作者：菜鸟追梦旅行 | 2024-02-15 20:35:45

赞

踩

拉普拉斯形变

背景

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

拉普拉斯网格变形原理

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

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

求：

  • 变形后的顶点位置 V ′ V' V′

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

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76
  • 77
  • 78
  • 79
  • 80
  • 81
  • 82
  • 83
  • 84
  • 85
  • 86
  • 87
  • 88

在这里插入图片描述

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

带有旋转的拉普拉斯变形

上一节展示的拉普拉斯矩阵有一个问题，当顶点变形涉及到旋转变化或者近似旋转的变化，即便形状一样，但是拉普拉斯距离已经发生了变化。
在这里插入图片描述
因此，在Laplacian Surface Editing原文论文中还提到了一种更准确的拉普拉斯变形方法。
首先假设每个顶点的旋转都是比较小，将旋转矩阵近似为一个简单形式，
在这里插入图片描述
然后要求旋转前后的拉普拉斯距离相同，则约束变成：
arg max ⁡ V ′ ∣ ∣ L V ′ − T L V ∣ ∣ + ∣ ∣ V 指定 ′ − A ∣ ∣ \operatorname{arg\,max}_V' || LV' - TLV|| + ||V'_{指定} - A|| argmaxV′​∣∣LV′−TLV∣∣+∣∣V指定′​−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

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76
  • 77
  • 78
  • 79
  • 80
  • 81
  • 82
  • 83

参考

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

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

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/菜鸟追梦旅行/article/detail/86760
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号