Laplacian_Surface_Editiing

概述: 为了解决传统的网格自由变形方法在曲面编辑操作通常会导致曲面细节严重失真的问题，本论文通过拉普拉斯坐标对网格进行编辑，并使用仿射变换矩阵调整拉普拉斯坐标，从而能在保留网格的细节和结构的前提下，对曲面进行平移，旋转和缩放等操作。本论文的方法能够用在网格编辑，涂层迁移和曲面块移植的场景下。
下图展示了把龙的嘴巴向上抬的过程。流程是首先固定红色区域，然后选择控制点，最后对局部区域进行网格变形。算法能够保持模型原有的特征。
在这里插入图片描述

1.预备知识

1.1 拉普拉斯坐标

拉普拉斯坐标是基于结点与其领域结点的相对位置关系进行构建的，所以拉普拉斯坐标是一种相对坐标。
好处在于，无论在空间中如何使网格变形，结点的相对坐标不会改变。
拉普拉斯坐标公式为：
$\delta_{i}=\mathscr{L}\left(\mathbf{v}_{i}\right) =\mathbf{v}_{i}-\frac{1}{d_{i}} \sum_{j \in \mathscr{N}_{i}} \mathbf{v}_{j} \tag{1}$
公式1可以理解为先求出领域结点绝对坐标的平均值，再用当前结点的绝对坐标减去这个均值，得到拉普拉斯坐标。

1.2 拉普拉斯坐标的性质

拉普拉斯坐标有三个性质。

1)线性变换
若 $T$ 是点集 $V$ 的线性变化，则 $T\left(L\left(v_{i}\right)\right) \equiv L\left(T\left(v_{i}\right)\right)$
2)平移不变性

若 $T$ 是平移变换，则
$\sum_{j \in N(i)} \omega_{i j}\left(T v_{i}-T v_{j}\right)=\sum_{j \in N(i)} \omega_{i j} T\left(v_{i}-v_{j}\right)=\sum_{j \in N(i)} \omega_{i j} T\left(\begin{array}{c} x_{v_{i}}-x_{v_{j}} \\ y_{v_{i}}-y_{v_{j}} \\ z_{v_{i}}-z_{v_{j}} \\ 0 \end{array}\right)=\sum_{j \in N(i)} \omega_{i j}\left(\begin{array}{c} x_{v_{i}}-x_{v_{j}} \\ y_{v_{i}}-y_{v_{j}} \\ z_{v_{i}}-z_{v_{j}} \\ 0 \end{array}\right)=\delta_{i}$

3)对旋转变换敏感,
若 $R$ 是旋转变换，则 $\sum_{j \in N(i)} \omega_{i j}\left(R v_{i}-R v_{j}\right)=R\delta_{i}$

2.算法原理

算法输入: 待变形的网格模型，控制网格变形前后的点序号和坐标，作为约束的固定的点序号和坐标；
算法输出: 变形后的网格模型。
算法目标：给定网格的绝对坐标 $V$ ，和拉普拉斯坐标，求解网格变形后的绝对坐标 $V^{\prime}$ 。
算法基本思路是：
1）将网格坐标（绝对坐标）编码为拉普拉斯坐标（相对坐标）
2）在拉普拉斯坐标中对网格进行处理
3）将拉普拉斯坐标解码回网格

可以用矩阵运算来表示网格中所有顶点的拉普拉斯坐标：
$\quad \Delta = L V \tag{2}$

其中， $L$ 是网格的拉普拉斯矩阵，定义为 $L=I-D^{-1} A$ ， $D$ 是这张网格的度矩阵， $A$ 是网格的邻接矩阵。可以发现， $R ank (L) = n - 1$ ，意味着无法直接求解线性方程 $L^{-1} \Delta$ 来还原顶点坐标。不过可以通过给这个线性方程添加约束，固定住其中的一些点来求解其他点，这些被固定的点称为控制点 (handle)。

为了实现局部编辑，我们只需固定一个点就能求解 $V$ ，但是实验发现求解速度太慢。如果固定一组点，求解速度加快，但求解 $V$ 会产生误差，我们希望最小化误差，并得到一个解，这个解就是对网格编辑后求的绝对坐标的结果。
也就是通过最小化能量函数 $E\left(V^{\prime}\right)$ 来求解新的顶点坐标 $V^{\prime}$ ：
$E\left(V^{\prime}\right)=\sum_{i=1}^{n}\left\|\delta_{i}-\mathscr{L}\left(\mathbf{v}_{i}^{\prime}\right)\right\|^{2}+\sum_{i=m}^{n}\left\|\mathbf{v}_{i}^{\prime}-\mathbf{u}_{i}\right\|^{2} \tag{3}$
该公式的前半部分表示 $V^{\prime}$ 的相对坐标与原⽹格的相对坐标 $\delta_{i}$ 尽可能保持⼀致，后半部分表示要让控制点在变形后处于目标位置 $\mathbf{u}_{i}$ 。

由于拉普拉斯坐标仅具有平移不变性，但对旋转和缩放敏感。为了解决这⼀问题，⽂章提出的⽅法是为每个顶点计
算⼀个基于新坐标 $V^{\prime}$ 的仿射变换矩阵 $T_i$ ，它记录了每个顶点及其邻域在变换过程中发⽣的旋转和缩放。公式3转换为下式：
$E\left(V^{\prime}\right)=\sum_{i=1}^{n}\left\|T_{i}\left(V^{\prime}\right) \delta_{i}-\mathscr{L}\left(\mathbf{v}_{i}^{\prime}\right)\right\|^{2}+\sum_{i=m}^{n}\left\|\mathbf{v}_{i}^{\prime}-\mathbf{u}_{i}\right\|^{2} \tag{4}$
这个式子中最明显的改变就是在变形前的拉普拉斯坐标 $δ_i$ 前乘上了仿射变换矩阵 $T_i$ 。仿射变换矩阵记录了每个顶点 $v_i$ 和其领域在原网格转换为新网格过程中发生的缩放和旋转变换。

注意这⾥ $T_i$ 和 $V^{\prime}$ 都是未知的，但由于 $T_i$ 是关于 $V^{\prime}$ 的线性函数，因此求解出 $V^{\prime}$ ⾃然可以推出 $T_i$ ：
$\min _{T_{i}}\left(\left\|T_{i} \mathbf{v}_{i}-\mathbf{v}_{i}^{\prime}\right\|^{2}+\sum_{j \in \mathscr{N}_{i}}\left\|T_{i} \mathbf{v}_{j}-\mathbf{v}_{j}^{\prime}\right\|^{2}\right) \tag{5}$
这个能量函数的前半部分代表了我们要让原顶点 $v_i$ 经过 $T_i$ 变换后能接近于新求解出来的坐标，后半部分表示对于顶点的邻接顶点也要保持有一样的效果。

拥有了这样的约束之后，我们的问题就是如何求解这个变换矩阵 $T_i$ 。 $T_i$ 可以写成：
$T_{i}=\left(\begin{array}{cccc} s & -h_{3} & h_{2} & t_{x} \\ h_{3} & s & -h_{1} & t_{y} \\ -h_{2} & h_{1} & s & t_{z} \\ 0 & 0 & 0 & 1 \end{array}\right)$
对于上面这个矩阵， $s$ 是缩放比率参数， $h$ 是三个轴的旋转参数, $t$ 是位移参数。我们通过求解应用这个矩阵来解决平移，旋转，缩放的不变性的问题。

我们可以将这个公式5转写为下面的矩阵形式来优化运算，将目标改为最小化下面的这个式子，即变换后的点和理
想点尽可能地相等。
$\left\|A_{i}\left(s_{i}, \mathbf{h}_{i}, \mathbf{t}_{i}\right)^{T}-\mathbf{b}_{i}\right\|^{2}$
其中 $A_i$ 是经过调整的邻接矩阵，即原顶点及其邻接点坐标的矩阵， $b_i$ 就是新求解出来的点和其邻接点的绝对坐标，向量 $s_i,h_i,t_i)$ 就是变换矩阵 $T_i$ 转写的线性组合形式： $\left(s_{i}, \overrightarrow{h_{i}}, \overrightarrow{t_{i}}\right)^{T}=\left(s_{i}, h_{i 1}, h_{i 2}, h_{i 3}, t_{i 1}, t_{i 2}, t_{i 3}\right)^{T}$
更具体地，
$A_{i}=\left(\begin{array}{ccccccc} v_{k_{x}} & 0 & v_{k_{z}} & -v_{k_{y}} & 1 & 0 & 0 \\ v_{k_{y}} & -v_{k_{z}} & 0 & v_{k_{x}} & 0 & 1 & 0 \\ v_{k_{z}} & v_{k_{y}} & -v_{k_{x}} & 0 & 0 & 0 & 1 \\ \vdots & & & & & & \end{array}\right), k \in\{i\} \cup \mathscr{N}_{i},$

$\mathbf{b}_{\mathbf{i}}=\left(\begin{array}{c} v_{k_{x}}^{\prime} \\ v_{k_{y}}^{\prime} \\ v_{k_{z}}^{\prime} \\ \vdots \end{array}\right), k \in\{i\} \cup \mathscr{N}_{i}$

那么，使⽤最⼩⼆乘法求解下面线性⽅程组就能得到 $s$ ， $h$ 和 $t$ 的值，进⽽计算出 $T_i$ 及 $V^{'}$ 的值。
$\left(s_{i}, \mathbf{h}_{i}, \mathbf{t}_{i}\right)^{T}=\left(A_{i}^{T} A_{i}\right)^{-1} A_{i}^{T} \mathbf{b}_{i} \tag{6}$

$A_i$ 包含已知点 $V_i$ ， $b_i$ 包含未知点 $V_i'$ .

3.算法实现

3.1 实验环境

编程环境：python
使用的库：igraph, networkx用于构建图，scipy.sparse.linalg.lsqr求解最小二乘法问题
电脑: ubuntu20.04

3.2 实现思路

知道了如何求解拉普拉斯坐标，和如何将拉普拉斯坐标还原为绝对坐标后，后续的操作也就简单了起来。
对于三维网格编辑，所需的操作就是先选择感兴趣的变形区域ROI（对局部区域进行编辑），得到ROI边界的顶点。选完ROI后，然后我们在网格中选择控制点，最后输入控制点移动到希望的目标位置。

拥有以上数据后，我们就可以构建线性方程组公式6，方程组系数矩阵 $L$ 就是由原顶点， ROI边界顶点和需要改变的控制点组成，其中原顶点满足的是能量函数的前半部分约束，边界和控制点统称为控制点，控制点满足的是能量函数后半部分的约束。最小化约束就可以还原出绝对坐标，也就是重建出网格编辑后的新顶点，将这些点应用到原网格上就完成了对网格的变形。

3.3 关键代码

def Laplacian_surface_editing():
    # 读取文件 并 对所有网格建图和相应坐标
    ply_data = PlyData.read("./ply/bunny.ply")
    glo_graph = mesh2graph(ply_data)
    glo_vtx_pos = np.vstack(list(vtx_attrs)[0:3] for vtx_attrs in ply_data.elements[0].data)
    
    #设置 handle 点（一个） 和 ROI区域
    # handle_vertex = {idx:[x, y, z]}
    handle_vertex = {14651: [-0.0123491, 0.153911, -0.0264322]}
    ROI_vertices = [15692, 7357, 9877, 28992]
    
    # 提取 ROI 子图
    border_vertices = get_border_vertices(glo_graph, ROI_vertices)
    edit_vertices   = get_edit_vertices(glo_graph, border_vertices, list(handle_vertex.keys())[0])
    edit_vertices   = border_vertices + edit_vertices
    sub_graph = glo_graph.subgraph(edit_vertices)
    
    # 全局/局部（待编辑区域） 序号映射
    edit_num  = len(edit_vertices)
    fixed_num = len(handle_vertex) + len(ROI_vertices)
    glo2sub = {}
    sub2glo = {}
    for i in range(edit_num):
        sub_vtx = sub_graph.vs[i]
        glo2sub[sub_vtx['index']] = i
        sub2glo[i] = sub_vtx['index']
    # 拉普拉斯算子
    L = calc_Lapla_oper(sub_graph, edit_num)
    V = np.vstack(glo_vtx_pos[sub2glo[i]] for i in range(edit_num))
    Delta = L.dot(V)
    
    L_prime = np.zeros([3 * (edit_num + fixed_num), 3 * edit_num])
    L_prime[0*edit_num : 1*edit_num, 0*edit_num : 1*edit_num] = L *(-1)
    L_prime[1*edit_num : 2*edit_num, 1*edit_num : 2*edit_num] = L *(-1)
    L_prime[2*edit_num : 3*edit_num, 2*edit_num : 3*edit_num] = L *(-1)
    
    for i in range(edit_num):
        # i在这里是点vertex
        # 邻接点
        nei_vtxs = sub_graph.neighbors(i)
        nei_vtxs.append(i)
        # 邻接矩阵
        Ai = np.zeros([3 * len(nei_vtxs), 7])
        for j in range(len(nei_vtxs)):
            pos = V[nei_vtxs[j]]
            Ai[j] = [pos[0],0,pos[2],-pos[1],1,0,0]
            Ai[j + len(nei_vtxs)] = [pos[1],-pos[2],0,pos[0],0,1,0]
            Ai[j + 2* len(nei_vtxs)] = [pos[2],pos[1],-pos[0],0,0,0,1]
    
        # 仿射变换矩阵
        Ti = np.linalg.inv(Ai.T.dot(Ai)).dot(Ai.T)
        si = Ti[0]
        hi = Ti[1:4]
        ti = Ti[4:7]
        
        # 拉普拉斯坐标 δ 分量
        Delta_ix = Delta[i, 0]
        Delta_iy = Delta[i, 1]
        Delta_iz = Delta[i, 2]
        # 公式5左半边部分Ti(V') * δ
        T_delta = [si*Delta_ix    - hi[2]*Delta_iy + hi[1]*Delta_iz,
                   hi[2]*Delta_ix + si*Delta_iy    - hi[0]*Delta_iz,
                  -hi[1]*Delta_ix + hi[0]*Delta_iy + si*Delta_iz]
        
        # 公式5左边Ti(V') * δ - L'
        nei_vtxs = np.array(nei_vtxs)
        row_idx = np.hstack([nei_vtxs, nei_vtxs + edit_num, nei_vtxs + 2*edit_num])
        L_prime[i + 0*edit_num, row_idx] += T_delta[0]
        L_prime[i + 1*edit_num, row_idx] += T_delta[1]
        L_prime[i + 2*edit_num, row_idx] += T_delta[2]
    
    # 设置 handle 点的约束方程
    b = np.array([])
    handle_vtxs = np.array(list(handle_vertex.values()))
    b = np.append(b, handle_vtxs) # handle 的 x y z 坐标
    
    handle_vtx_sub_idx = glo2sub[list(handle_vertex.keys())[0]]

    L_prime[3 * edit_num + 0, 0 * edit_num + handle_vtx_sub_idx] = 1
    L_prime[3 * edit_num + 1, 1 * edit_num + handle_vtx_sub_idx] = 1
    L_prime[3 * edit_num + 2, 2 * edit_num + handle_vtx_sub_idx] = 1
        
    # 设置 ROI 点的约束方程
    for i in range(len(ROI_vertices)):
        b = np.append(b, V[glo2sub[ROI_vertices[i]]])
        ROI_vtx_sub_idx = glo2sub[ROI_vertices[i]]
        L_prime[(3 * (edit_num + 1)) + (i * 3) + 0, ROI_vtx_sub_idx + (0*edit_num)] = 1
        L_prime[(3 * (edit_num + 1)) + (i * 3) + 1, ROI_vtx_sub_idx + (1*edit_num)] = 1
        L_prime[(3 * (edit_num + 1)) + (i * 3) + 2, ROI_vtx_sub_idx + (2*edit_num)] = 1
    
    # 最小二乘法求解
    b = np.hstack([np.zeros(3 * edit_num), b])
    A = scipy.sparse.coo_matrix(L_prime)
    V_prime = scipy.sparse.linalg.lsqr(A, b)[0]
    
    # 将编辑后的坐标重新代入全局坐标中
    for i in range(edit_num):
        glo_vtx_pos[sub2glo[i]] = [V_prime[i], V_prime[i + edit_num], V_prime[i + (2*edit_num)]]
    # 写入文件
    for i in range(glo_vtx_pos.shape[0]):
        ply_data.elements[0].data[i] = tuple(glo_vtx_pos[i]) + tuple(ply_data.elements[0].data[i])[3:]
    ply_data.write("./result/res.ply")

实现结果：展示兔子耳朵弯曲的过程。（左边为原图，右边为变形后的图）
在这里插入图片描述

4.心得体会

知识可以触类旁通。Laplacian Surface Editing 和 Poisson Image Editing 的思想很像，都是通过保持离散点的Laplacian Coordinates 进行模型/图像重建，从而保证细节完整。
将公式写成代码的时候，要仔细。比如列表序号应从0开始；矩阵相乘要注意维度。
学会用igraph从网格中构建图。
学会用工具Blender查找坐标点。

5.参考文献

本文含有隐藏内容，请开通VIP 后查看

【计算机图形学】Laplacian_Surface_Editiing拉普拉斯曲面编辑算法