泊松变形技术理论知识解读

一、简单介绍

微分域变形领域最重要的两篇论文分别是：

[1] Yu Y , Zhou K , Xu D , et al. Mesh editing with poisson-based gradient field manipulation[J]. ACM Transactions on Graphics, 2004, 23(3):644. [ciation:679]
[2] Olga Sorkine-Hornung, Daniel Cohen-Or, Yaron Lipman, Marc Alexa, Christian Rössl, Hans-Peter Seidel . Laplacian surface editing[C]. SGP// 2004. [citation:1231]

两篇文章分别以经典的泊松方程和拉普拉斯方程为出发点很好的解决网格编辑的问题，形式上非常漂亮。两篇论文的作者也都是图形学领域非常著名的学者，按Hu老师上课的说法，那是一个风云变幻的时代。图形学界公认最好的期刊是TOG，但你会发现发在SGP上的laplace几乎是发在TOG上poisson的2倍！为什么呢？那是因为poisson方法在理论上相对于laplace方法过于抽象，因此对于很多人难以理解和实现。

这篇博客重点讲述泊松变形技术，也希望通过开源代码可以帮助到更多的初学者。（如果对你有帮助，请给我star！）
代码地址：https://github.com/Zhengjun-Du/PoissonDeformation

首先我们看一下泊松变形的效果,其中白色部分为固定不动的部分，红色部分为手动交互的部分。

泊松变形交互展示，原地址：(https://www.bilibili.com/video/BV1H54y1Q7dF/)

二、背景工作

泊松变形的背景工作实际上是poisson image editing：

简单回顾一下，这篇论文主要实现图像克隆，把一张图像的某个区域拷贝到另一张图像中，实现完美的融合。方法分别利用待求解像素的（上图g区域）laplace算子和源图像梯度场（上图v区域）的散度构建了一个拉普拉斯方程，并以目标图像的边界为约束条件，通过求解线性方程组实现光滑的过渡。

具体过程可以参考：图像的泊松(Poisson)编辑、泊松融合，文章很详细，通俗易懂。

图像空间，无论是梯度、散度还是laplace算子的计算都非常直观。但拓展到三维空间却不是那么的直观和简明。
因此难点主要体现在三方面：
1、如何定义三角网格上的梯度场和散度？
2、如何定义三角网格的laplace算子?
3、如何通过交互修改散度

三、三角网格梯度场

首先要清楚的是三角网格本身是一个分段线性的表示，它本质上不是一个光滑的曲面，它是由很多三角形平面组合而成的。要定义网格梯度场，首先要定义实值函数。

在每个三角形内部，其实可以通过三角形三个顶点的属性插值得到。如下所示， $f_i, f_j,f_k$ 分别是三个顶点关联的属性值， $f_u$ 可以通过 $α,β,γ\alpha, \beta, \gamma$ 的线性插值得到，一般情况下， $α,β,γ\alpha, \beta, \gamma$ 分别可以表示为三个小三角形与大三角形的面积之比。

因为 $α,β,γ\alpha, \beta, \gamma$ 分别表示 $v_i, v_j,v_k$ 对 $v_u$ 的影响，因此这里我们将 $α,β,γ\alpha, \beta, \gamma$ 定义为三个函数： $B_i(u)$ ， $B_j(u)$ 和 $B_k(u)$ , 同时也满足： $B_i(u)+B_j(u)+B_k(u) = 1$ .

为此，三角网格上任意一点可以表示为： $f(u)=fiBi(u)+fjBj(u)+fkBk(u)(1)f(u)=f_iB_i(u)+f_jB_j(u)+f_kB_k(u)\tag1$

接下来，可以得到梯度： $∇f(u)=fi∇Bi(u)+fj∇Bj(u)+fk∇Bk(u)(2)\nabla f(u)=f_i\nabla B_i(u)+f_j\nabla B_j(u)+f_k\nabla B_k(u)\tag2$

因为： $B_i(u)+B_j(u)+B_k(u) = 1$

所以： $∇Bi(u)+∇Bj(u)+∇Bk(u)=0(3)\nabla B_i(u)+\nabla B_j(u)+\nabla B_k(u) = 0\tag3$

因此， $∇Bi(u)=−∇Bj(u)−∇Bk(u)\nabla B_i(u) =-\nabla B_j(u)-\nabla B_k(u)$

将 (5) 代入 (2) : $∇f(u)=(fj−fi)∇Bj(u)+(fk−fi)∇Bk(u)(4)\nabla f(u)=(f_j-f_i)\nabla B_j(u)+(f_k-f_i)\nabla B_k(u)\tag4$

那么，问题来了， $∇Bj(u)\nabla B_j(u)$ 又该怎么计算呢？ 下面以 $∇Bi(u)\nabla B_i(u)$ 为例讲解。

如上左图所示，我们先做一条平行于 $v_jv_k$ 的直线 $l$ , 在平行线上取两个点： $v_x,v_y$ . 那么有：
${fx=fiBi(x)+fjBj(x)+fkBk(x)fy=fiBi(y)+fjBj(y)+fkBk(y)(5)\left\{ \begin{aligned} f_x = f_iB_i(x) + f_jB_j(x) + f_kB_k(x)\\ f_y = f_iB_i(y) + f_jB_j(y) + f_kB_k(y) \end{aligned} \tag5 \right.$

我们知道：
${Bi(x)=s△xjk/s△ijkBi(y)=s△yjk/s△ijk(6)\left\{ \begin{aligned} B_i(x) = s\triangle xjk / s\triangle ijk\\ B_i(y) = s\triangle yjk / s\triangle ijk \end{aligned} \tag6 \right.$
其中， $s△s\triangle$ 表示三角形的面积。因为 $s△xjks\triangle xjk$ 和 $s△yjks\triangle yjk$ 同底等高, 所以： $s△xjk=s△yjks\triangle xjk =s\triangle yjk$ ，因此：
$B_i(x) = B_i(y)$

基于这个重要的结论，实际上所有平行于 $v_x,v_y$ 的直线上的 $B_i$ 均相等。为此，我们可以做出右侧的一个关于 $B_i$ 的等值线。那么我们就可以知道， $B_i$ 梯度方向是垂直于等值线且指向顶点 $v_i$ 的。在长度为底边高的线段上从0变化到1.那么，梯度的模长为高的倒数。因此：

$∇Bi(u)=(xk−xj)⊥/(2A)(7)\nabla B_i(u)=(x_k-x_j)^\perp/(2A)\tag7$

将（7）代入（4）可以得到三角面片的梯度场为：

$∇f(u)=(fj−fi)(xi−xk)⊥/(2A)+(fk−fi)(xj−xi)⊥/(2A)(8)\nabla f(u)=(f_j-f_i)(x_i-x_k)^\perp/(2A)+(f_k-f_i)(x_j-x_i)^\perp/(2A)\tag8$

其中， $⊥\perp$ 表示方向垂直底边指向另一个顶点， $A$ 表示三角形的面积。

四、梯度场的散度

梯度场在点 $v_i$ 处的散度定义如下，它实际上是1-邻域的求和。

$\boldsymbol w(v_i)=\sum_{T\in N(v_i)}\nabla B_i|_T\cdot \boldsymbol w_T\cdot A_T\tag9$

其中， $∇Bi\nabla B_i$ 如（7）所示， $wT\boldsymbol w_T$ 如（8）所示， $A_T$ 为三角形的面积， $T$ 为 $v_i$ 关联的三角面片。

五、顶点的laplace算子及线性系统

顶点的laplace算子定义如下左图所示，那么可以写成矩阵的形势，则权重矩阵A定义为：

$Aij={1/2⋅∑vj∈Ni(cotαj+cotβj),i=j−1/2⋅(cotαj+cotβj),i≠j(10)A_{ij} = \begin{cases} 1/2\cdot\sum_{v_j\in N_i}(cot\alpha_j + cot\beta_j) ,& i = j\\ -1/2\cdot(cot\alpha_j + cot\beta_j) ,& i \not= j\end{cases}\tag{10}$

可以证明：三角网格顶点的laplace算子等于散度。为此，可以构建如下右图所示的线性方程组。 $A x$ 表示三角网格在任意顶点的laplace算子， $b$ 为顶点的散度。

在实际的变形过程中，用户通过交互间接修改三角网格的散度，而左侧权重矩阵 $A$ 保持不变，通过求解上述线性方程组求解变形后顶点的坐标。我们注意到：无论是laplace算子还是散度，都只刻画了顶点与one-ring结点的局部关系，因此泊松变形可以保持局部细节。

六、如何修改散度

实际变形中，如下图所示，我们只需要对 $B C$ (局部控制点)进行旋转、平移等简单的操作来驱动网格变形。那么上面我们提到，变形过程通过修改散度来求解变形后网格顶点的坐标。那么，这里的核心问题就是如何通过操作局部的控制点来修改其它地方三角网格的散度。

其实方法有很多，我是对每个点设置了形变系数（0，1）之间，红色控制点处设置为1，然后求解调和方程将形变系数光滑过渡到网格的所有顶点。形变系数的核心意义就是，离控制点越近形变越大，离控制点越远，形变越小。然后利用四元数插值，得到所有点的变形矩阵。任意三角面片的三个顶点一旦改变，那么就可以求解新的散度。然后就可以求解线性方程组了。

参考文献

[1] Yu Y , Zhou K , Xu D , et al. Mesh editing with poisson-based gradient field manipulation[J]. ACM Transactions on Graphics, 2004, 23(3):644.
[2] Press C. Polygon Mesh Processing[J]. Crc Press, 2010.
[3] Pérez, Patrick, Gangnet, Michel, Blake, Andrew. Poisson image editing[J]. Acm Transactions on Graphics, 22(3):313
[4] Tong Y , Lombeyda S , Hirani A N , et al. Discrete Multiscale Vector Field Decomposition[J]. ACM Transactions on Graphics, 2003, 22(3):445-452.
[5] Cohenor D. Laplacian surface editing[C]// 2004.