Mesh Editing with Poisson-Based Gradient Field Manipulation

原文：http://blog.sina.com.cn/s/blog_73428e9a0101h4l9.html

参考论文：Mesh Editing with Poisson-Based Gradient Field Manipulation. Yi ZhouYu, Kun Zhou

需要源代码请联系我：duzjqhu@aliyun.com ,发邮件的时候请说清楚三件事： 你是谁？ 你来自哪里？ 用代码做什么？

实验效果：

  

          

                                           图1. 泊松变形效果图

  是时候重新写这篇文章了，去年我曾写过一篇泊松变形的文章，并且公布了源码，博客发表后收到很多小伙伴的来信，让我受宠若惊！我本人在研究生阶段的研究方向就是网格变形，我个人喜欢开源，也很喜欢跟大家一起讨论，我觉得分享是一件有趣的事情，大家可以相互交流、相互学习，而现在很多有名的论文尤其是关于几何处理的算法很少公布源码，比如MSRA,ZJU等虽然SIGGRAPH论文一大把，但是难得开源。而文章里的描述很多都假定你有较好的数学基础，而对于很多初学者而言，这是非常困难的事情，之前公布的源码实现较为繁琐，并且文章写得不够具体，故现在重新梳理之，希望对相关领域刚刚入门的小伙伴有一定的帮助。

     这篇文章是发在SIGGRAPH 2004上的，至今引用次数400+，可以说是几何处理领域一篇比较不错的文章，2003年P Pére等人发了一篇关于图像处理的文章：Poisson Image Editing，这篇文章的引用次数至今已经1000+！这篇文章就是将poisson图像处理的方法扩展到三维的几何处理。说来惭愧，花了一个月时间才实现了这篇论文，实现文章的时候给几位作者都发信问过一些细节，但是都没有回复。在实现论文的过程中，计算机系的Lin Gao博士给我提供了一些必要的帮助，特此感谢！

     文章的核心即是求解如下所示的泊松方程实现模型变形，该泊松方程实际上是微分几何中一个常见的公式，大意是说定义在三角网格上的分段线性函数 f 的拉普拉斯算子等于三角网格上常矢量场w的散度。

                       （1）

 其中，w表示定义在网格表面的常矢量场，▽w表示该矢量场的散度， f 是方程的未知数即待求解的变形后的网格模型， △f 表示网格的拉普拉斯算子，f*是该泊松方程的狄利克雷边界值，同时也是变形的位置约束。下面我们就对这个方程进行分解：

拉普拉斯算子

    首先，对于单个的网格顶点而言，其拉普拉斯算子定义为其与1-邻域顶点的质心偏差，如图2所示，这里我们采用常用的余切权拉普拉斯算子的形式，如公式（2）所示：

  

图2   余切权拉普拉斯算子

            （2）

 图示的拉普拉斯算子仅仅是针对Pi点，对于整个网格而言，其拉普拉斯算子可以写成3所示的矩阵的乘积的形式 :

                （3）

其中，

     实际上3中左端矩阵的第i行即是网格上第i的顶点的拉普拉斯算子，读者可以自行验证，这里的△P即是方程(1)中的△f，我们也将(3)中右端的第一个矩阵称为拉普拉斯矩阵。

三角网格的梯度场

     接下来，定义在网格上的矢量场 w 一般使用三角网格的梯度代替，那么三角网格的梯度是什么样的呢？

 三角网格梯度的描述基于分段线性函数的假设，即假设三角网格的任意顶点都定义了一个分段线性函数f, 同时该函数满足f (p_x )= f_x，并且越远离p_x，f (p_x )渐进的减小直到为0。

 

图3   三角网格上任意三角形

  对于图3所示的三角网格上的任意三角形△p_i p_j p_k，我们知道对于其中的任一点v都可以由其所在三角形的三个顶点线性插值得到：

             （4）

相应的，其梯度可以表示为：

由重心坐标的性质可知φ_i + φ_j + φ_k  = 1 ，因此 ▽φ_i+▽φ_j+▽φ_k= 0 ，于是用 ▽φ_j、▽ φ_k 表示 ▽φ_i 进一步得到：

而定义在三角网格顶点的函数φ的梯度大小等于该点对边的高的倒数，方向从对面垂直指向该点。

     如图3所示，我们可以这样理解：三角形上任一点v可以由三个顶点p_i p_j p_k插值得到，如公式(4)所示，即v同时受到三个顶点的影响，对于φ_i 而言，我们知道，当v离p_i 越近，φ_i 越大，极端情况下，v和p_i 重合，此时φ_i = 1，φ_j = φ_k = 0，反过来，当接近 p_i 的对边时，φ_i = 0，因为这时v仅仅是p_j p_k线性插值的结果。因此，φ_i 在p_i处取得最大值1，而在对边p_j p_k取得最小值0，因此梯度方向为p_j p_k垂直指向p_i，大小等于1/|bipi|。 

     如图3所示，pj点的梯度方向即为蓝色线段b_jp_j的指向，大小等于1/|b_jp_j|。因此上述的等式进一步可以转化为：

而这里我们将f_i,、f_j、f_k分别取为p_i、p_j、p_k的坐标值，于是:

至此，我们得到了三角网格上的梯度场。

梯度场的散度

     接下来，有了上面的知识，梯度的散度计算就比较容易了,三角网格矢量场的散度定义为：

 

     其中▽Bik表示定义在p_i 点的分段线性函数在三角形T_k 上的梯度，如图3中b_iv_{i向量所示}，这个梯度是一个常数，大小等于T_k 中p_i 对边的高的倒数，方向由p_i 垂直指向其对边，A_k 是三角形T_k 的面积。

泊松方程的求解

 

                                                        (5)

    对泊松方程(1)进一步细化得到等式(5)所示的形式，泊松方程求解过程中，拉普拉斯矩阵只需计算一次，此后保持不变。等式左边的列矩阵即是方程的未知数，同时也是待求解的变形后的模型。变形的实际就是通过用户交互改变网格的梯度场，进而得到改变的散度，通过求解泊松方程获取变形后的模型。

    那么用户的交互如何影响网格的梯度场呢？对梯度场的操纵是泊松网格编辑中的一个关键步骤，它决定最终网格的形状，如果对网格模型中的每一个三角片都定义一个局部变换的话，几乎是一个不可能完成的任务，因此将局部的变换自动扩展到全局网格显得很有必要。论文中设计了Uniform，Linear， Gaussian三种传播局部变换的方法将控制区域的变化传播到兴趣区域的各个三角形面，从而改变兴趣区域的梯度场。

    变形过程大致描述为：用户交互地操作模型的控制点，通过线性或热辐射扩散的方法将控制点的变换扩散到模型的兴趣区域，进而模型表面的矢量场和散度相应的发生改变，即泊松方程的右侧部分发生改变，每次变换后，求解该方程得到最终的变形结果。

存在的问题

泊松变形技术通过改变矢量场求解原始网格的标量场，作为一个网格编辑框架不仅可以实现网格变形，同时在网格拼接、网格去噪方面都有很好的应用。但是泊松重建方法也有其缺点：

(1)  Laplacian坐标不具有旋转不变性，网格局部旋转的时候，泊松变形方法将原本非线性的能量优化为线性变形，因此对于一些变形显得不直观；

(2) 相对于拉普拉斯变形交互较为复杂；

(3) 平移不变性造成的变形畸形，因为平移不改变模型的梯度场，因此平移不敏感。如图4所示：锥形体的右侧平移，而变形后的模型，锥形体不再垂直于基面，而保持了最初的朝向。

              图4   梯度场的平移不变性造成变形不合理

     大概就这样了，感觉依然写得不太好，不过大致求解方法基本说清楚。当初脑残选择了新浪博客，根本不适合写技术博客，里面很多公式都是粘贴的图片。。。。以后还是搬去博客园吧。希望大家提问，我乐于跟大家讨论。

参考文献：

[1]  Yu Y, Zhou K, Xu D, et al. Mesh editing with poisson-based gradient field  manipulation[C]//ACM Transactions on Graphics. ACM, 2004, 23(3): 644-651.

[2]  Botsch M, Kobbelt L, Pauly M, et al. Polygon mesh processing [M]. CRC press, 2010.

[3]  许栋. 微分网格处理技术[D]. 杭州: 浙江大学, 2006.

[4]  Botsch M, Sorkine O. On linear variational surface deformation methods[J]. Visualization and Computer Graphics, IEEE Transactions on, 2008, 14(1): 213-230.