Laplacian_Surface_Editiing

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

1.预备知识

1.1 拉普拉斯坐标

拉普拉斯坐标是基于结点与其领域结点的相对位置关系进行构建的，所以拉普拉斯坐标是一种相对坐标。
好处在于，无论在空间中如何使网格变形，结点的相对坐标不会改变。
拉普拉斯坐标公式为：
δi=L(vi)=vi−1di∑j∈Nivj(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} δi=L(vi)=vi−di1j∈Ni∑vj(1)
公式1可以理解为先求出领域结点绝对坐标的平均值，再用当前结点的绝对坐标减去这个均值，得到拉普拉斯坐标。

1.2 拉普拉斯坐标的性质

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

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

若TTT是平移变换，则
∑j∈N(i)ωij(Tvi−Tvj)=∑j∈N(i)ωijT(vi−vj)=∑j∈N(i)ωijT(xvi−xvjyvi−yvjzvi−zvj0)=∑j∈N(i)ωij(xvi−xvjyvi−yvjzvi−zvj0)=δi\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} j∈N(i)∑ωij(Tvi−Tvj)=j∈N(i)∑ωijT(vi−vj)=j∈N(i)∑ωijT⎝⎛xvi−xvjyvi−yvjzvi−zvj0⎠⎞=j∈N(i)∑ωij⎝⎛xvi−xvjyvi−yvjzvi−zvj0⎠⎞=δi

3)对旋转变换敏感,
若RRR是旋转变换，则 ∑j∈N(i)ωij(Rvi−Rvj)=Rδi\sum_{j \in N(i)} \omega_{i j}\left(R v_{i}-R v_{j}\right)=R\delta_{i}∑j∈N(i)ωij(Rvi−Rvj)=Rδi

2.算法原理

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

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

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

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

由于拉普拉斯坐标仅具有平移不变性，但对旋转和缩放敏感。为了解决这⼀问题，⽂章提出的⽅法是为每个顶点计
算⼀个基于新坐标V′V^{\prime}V′的仿射变换矩阵TiT_iTi，它记录了每个顶点及其邻域在变换过程中发⽣的旋转和缩放。公式3转换为下式：
E(V′)=∑i=1n∥Ti(V′)δi−L(vi′)∥2+∑i=mn∥vi′−ui∥2(4)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} E(V′)=i=1∑n∥Ti(V′)δi−L(vi′)∥2+i=m∑n∥vi′−ui∥2(4)
这个式子中最明显的改变就是在变形前的拉普拉斯坐标δiδ_iδi前乘上了仿射变换矩阵TiT_iTi。仿射变换矩阵记录了每个顶点viv_ivi和其领域在原网格转换为新网格过程中发生的缩放和旋转变换。

注意这⾥TiT_iTi和V′V^{\prime}V′都是未知的，但由于TiT_iTi是关于V′V^{\prime}V′的线性函数，因此求解出V′V^{\prime}V′⾃然可以推出TiT_iTi：
min⁡Ti(∥Tivi−vi′∥2+∑j∈Ni∥Tivj−vj′∥2)(5)\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} Timin⎝⎛∥Tivi−vi′∥2+j∈Ni∑∥∥Tivj−vj′∥∥2⎠⎞(5)
这个能量函数的前半部分代表了我们要让原顶点viv_ivi经过TiT_iTi变换后能接近于新求解出来的坐标，后半部分表示对于顶点的邻接顶点也要保持有一样的效果。

拥有了这样的约束之后，我们的问题就是如何求解这个变换矩阵TiT_iTi。TiT_iTi可以写成：
Ti=(s−h3h2txh3s−h1ty−h2h1stz0001)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) Ti=⎝⎛sh3−h20−h3sh10h2−h1s0txtytz1⎠⎞
对于上面这个矩阵， sss是缩放比率参数， hhh是三个轴的旋转参数, ttt是位移参数。我们通过求解应用这个矩阵来解决平移，旋转，缩放的不变性的问题。

我们可以将这个公式5转写为下面的矩阵形式来优化运算，将目标改为最小化下面的这个式子，即变换后的点和理
想点尽可能地相等。
∥Ai(si,hi,ti)T−bi∥2\left\|A_{i}\left(s_{i}, \mathbf{h}_{i}, \mathbf{t}_{i}\right)^{T}-\mathbf{b}_{i}\right\|^{2}∥∥Ai(si,hi,ti)T−bi∥∥2
其中AiA_iAi是经过调整的邻接矩阵，即原顶点及其邻接点坐标的矩阵， bib_ibi就是新求解出来的点和其邻接点的绝对坐标，向量(si,hi,ti)(s_i,h_i,t_i)(si,hi,ti)就是变换矩阵TiT_iTi转写的线性组合形式：(si,hi→,ti→)T=(si,hi1,hi2,hi3,ti1,ti2,ti3)T\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}(si,hi,ti)T=(si,hi1,hi2,hi3,ti1,ti2,ti3)T
更具体地，
Ai=(vkx0vkz−vky100vky−vkz0vkx010vkzvky−vkx0001⋮),k∈{i}∪Ni,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}, Ai=⎝⎛vkxvkyvkz⋮0−vkzvkyvkz0−vkx−vkyvkx0100010001⎠⎞,k∈{i}∪Ni,

bi=(vkx′vky′vkz′⋮),k∈{i}∪Ni\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} bi=⎝⎛vkx′vky′vkz′⋮⎠⎞,k∈{i}∪Ni

那么，使⽤最⼩⼆乘法求解下面线性⽅程组就能得到sss，hhh和ttt的值，进⽽计算出TiT_iTi及V′V'V′的值。
(si,hi,ti)T=(AiTAi)−1AiTbi(6)\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} (si,hi,ti)T=(AiTAi)−1AiTbi(6)

AiA_iAi包含已知点ViV_iVi，bib_ibi包含未知点Vi′V_i'Vi′.

3.算法实现

3.1 实验环境

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

3.2 实现思路

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

拥有以上数据后，我们就可以构建线性方程组公式6，方程组系数矩阵LLL就是由原顶点， 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_verticessub_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']] = isub2glo[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] = 1L_prime[3 * edit_num + 1, 1 * edit_num + handle_vtx_sub_idx] = 1L_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)] = 1L_prime[(3 * (edit_num + 1)) + (i * 3) + 1, ROI_vtx_sub_idx + (1*edit_num)] = 1L_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.参考文献

laplacian surface editing
【笔记】《Laplacian Surface Editing》的思路
图形学高被引论文赏析系列2：Laplacian Surface Editing
Laplacian surface editing博客
拉普拉斯网格变形
图形处理（三）简单拉普拉斯网格变形-Siggraph 2004
Laplacian_Surface_Editing_reproduce
Laplacian-Surface-Editing
laplacian-surface-editing code2

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

c语言连通域算法递归,VC++ 6.0编写计算机图形学中的种子填充算法，想用递归的八向连通域，求助！...
VC++ 6.0编写计算机图形学中的种子填充算法,想用递归的八向连通域,求助!0 填充函数代码如下: void CComputerGraphicsView::PolygonFill2()//区域填充函 ...
计算机图形学裁剪算法源代码,OpenGL计算机图形学梁友栋裁剪算法实验代码及运行结果.doc...
OpenGL计算机图形学梁友栋裁剪算法实验代码及运行结果.doc (10页) 本资源提供全文预览,点击全文预览即可全文预览,如果喜欢文档就下载吧,查找使用更方便哦! 14.9 积分 .<计算 ...
计算机图形学直线裁剪原理,计算机图形学-3.2用Liang-Barsky算法实现直线段裁剪...
计算机图形学-3.2用Liang-Barsky算法实现直线段裁剪计算机图形学-3.2用Liang-Barsky算法实现直线段裁剪 (1)算法设计原理依次处理(p1,q1)(p2,q2)(p3,q3 ...
计算机图形学05：中点BH算法对任意斜率的直线扫描转换方法
作者:非妃是公主专栏:<计算机图形学> 博客地址:https://blog.csdn.net/myf_666 个性签:顺境不惰,逆境不馁,以心制境,万事可成.--曾国藩文章目录专栏推 ...
TIT 计算机图形学实验三使用重心坐标算法绘制颜色渐变的正六面体
TIT 计算机图形学实验三使用重心坐标算法绘制颜色渐变的正六面体前言参考视频计算机图形学全套算法讲解和C++编码实现(共23讲配套源码),计算机图形学案例视频讲解以及主页相关算法.孔老师是我的 ...
图形学画直线c语言,002计算机图形学之直线画线算法
002计算机图形学之直线画线算法我们知道直线方程的斜截式是如下的样子: y = kx +b 在显示器上显示直线的话,如果使用如上的方程,每描一个点需要进行一次浮点乘法,一次浮点加法,和取整操作. ...
计算机图形学空间曲线,部分计算机图形学参数曲线和曲面.ppt
部分计算机图形学参数曲线和曲面参数曲线和曲面基础 Ray ray@mail.buct.edu.cn 内容曲线曲面参数表示位置矢量.切矢量.法矢量.曲率和挠率插值.拟合.逼近和光顺参数曲线的代 ...
计算机图形学直线段的生成算法
计算机图形学直线段的生成算法C++实现,包括:DDA,中点画线,改进的Bresenham画线文章目录 1.实验目的和内容 1.1实验目的 1.2实验内容 2.算法原理 2.1 DDA(数值微分算法) ...
计算机图形学直线算法论文,《计算机图形学》中直线生成算法的教学心得
摘要:<计算机图形学>是计算机科学与技术专业一门重要的专业课,其中直线生成算法是教学重点之一.该文通过分析几种直线生成算法的特点,阐述了理论教学和实践教学的重点和难点,总结了教学的体会和心 ...
[OpenGL]计算机图形学：明暗处理的基本算法
下面来逐个介绍一下基本要点以便加深对于明暗处理的基本算法. 光照方程: 1.光的无穷次反射和吸收过程可以用光照方程(rendering equation)描述.一般是无法求解的,即使应用数值方法也是非 ...

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