[图形学]拉普拉斯网格变形（Laplace Deformation）原理及复现

本文原理参考论文： Laplacian Mesh Processing

本文复现参考：Mesh模型的Laplace Deformation（网格形变 - 拉普拉斯形变） - C++代码实现

文章目录

拉普拉斯网格变形原理
- 拉普拉斯坐标
- 拉普拉斯矩阵
- 原理分析
- 添加锚点
- 求解形变后的坐标

拉普拉斯网格变形是一个相当经典的工作，本工作复现了一下基础的部分,这里记录下原理部分和结果。
复现结果演示如下：

拉普拉斯网格变形原理

可以用用一句话来解释拉普拉斯网格变形的原理： 形变前的点的拉普拉斯坐标与形变后的点的拉普拉斯坐标尽可能相等。这句话里面其实只有点的拉普拉斯坐标是我们所不清楚的，下面阐述拉普拉斯坐标为何物。

拉普拉斯坐标

拉普拉斯坐标是一个相当简单的概念，对于一个三角网格模型M=(V,E,F)，V为顶点集，E为边集，F为三角面片集。V（v1,…vn）中每一个点的坐标表示都是笛卡尔坐标，下图(源自论文Laplacian Mesh Processing)和公式定义了拉普拉斯坐标的表示，点vi的笛卡尔坐标减去所有vi的相邻点vj的笛卡尔坐标的平均值，di表示vi的相邻点的个数。

δi=(δi(x),δi(y),δi(z))=vi−1di∑j∈N(i)vj\delta_{i}=\left(\delta_{i}^{(x)}, \delta_{i}^{(y)}, \delta_{i}^{(z)}\right)=\mathbf{v}_{i}-\frac{1}{d_{i}} \sum_{j \in N(i)} \mathbf{v}_{j} δi=(δi(x),δi(y),δi(z))=vi−di1j∈N(i)∑vj
拉普拉斯坐标蕴含了网格模型中的细节信息，而细节信息正是我们在网格形变后不希望改变的，拉普拉斯网格变形正是基于这一约束。

拉普拉斯矩阵

拉普拉斯坐标的计算我们可以按照上面的公式分别求解，但是为了简化计算，通常用矩阵乘法计算拉普拉斯坐标。只需要定义拉普拉斯矩阵L，便可求解拉普拉斯坐标 δ = L*V ( V表示所有点构成的笛卡尔坐标矩阵)。下面给出L的定义：

L=I−D−1AL=I-D^{-1} A L=I−D−1A
其中，I为单位矩阵，D矩阵是一个对角矩阵，对角线上的值Dii = di，di是第i个点的邻接点数目，A 矩阵是邻接矩阵。若点vi和vj相邻，Aij = 1，若不相邻，则 Aij = 0。
Aij={1(i,j)∈E0otherwise A_{i j}=\left\{\begin{array}{ll}{1} & {(i, j) \in E} \\ {0} & {\text { otherwise }}\end{array}\right. Aij={10(i,j)∈E otherwise
在计算机中求矩阵的逆是比较麻烦的，因此再使用上面公式的时候我们通常两边同乘D将公式变成如下形式。之后便用Ls代替L。
LS=DL=D−AL_{S}=D L=D-A LS=DL=D−A

(LS)ij={dii=j−1(i,j)∈E0otherwise \left(L_{S}\right)_{i j}=\left\{\begin{array}{ll}{d_{i}} & {i=j} \\ {-1} & {(i, j) \in E} \\ {0} & {\text { otherwise }}\end{array}\right. (LS)ij=⎩⎨⎧di−10i=j(i,j)∈E otherwise

原理分析

回到一开始说的拉普拉斯原理：形变前的点的拉普拉斯坐标与形变后的点的拉普拉斯坐标尽可能相等。现在将这句话用公式表达出来：

(L)nn * Xn3 = (L)nn * Cn3 = Bn3

用Ls矩阵代替L矩阵，等式两边同乘D，变为下式：

(Ls)nn * Xn3 = (Ls)nn * Cn3 = D*Bn3

n表示网格中点的个数，Lnn表示拉普拉斯坐标矩阵大小为n*n，(Ls)nn=DL( L为拉普拉斯矩阵，D为上述描述的对角阵)，X表示形变后的点的笛卡尔坐标（需要求解的目标），C表示形变前的点的笛卡尔坐标，B表示形变前的拉普拉斯坐标。

写到上面的矩阵的形式，是不是非常熟悉？这就是矩阵中最为常见的Ax=b的形式啊，对应到这个问题就是LX=B。

这里简单回顾一下矩阵里求解Ax=b的方法，当A是可逆的： x= A~ b（这里用A~表示A的逆），当A不是方阵的时候，通常我们用最小二乘求解最优解，等式两边同乘以AT ：ATAx=ATb−−>x=(ATA)−1ATbA^{T}A x=A^{T} b --> x = (A^{T}A)^{-1}A^{T} bATAx=ATb−−>x=(ATA)−1ATb（符号-1表示逆，ATA^{T}AT表示矩阵A的转置）

在拉普拉斯变换中如果Ls是可逆的我们就可以直接计算形变后的点的坐标矩阵 X = (Ls)~ * D * Bn3

但是很遗憾，通常Ls是不可逆的，对于使得Ls可逆的解释原文是那么说的

为了唯一地还原全局坐标，需要解决一个完整的线性系统。假设模型是相连的的，我们需要指定一个顶点的笛卡尔坐标来解决平移自由度。替换顶点i的坐标等效于从L删除第i行和第i列，这使矩阵可逆。

为什么需要增加锚点？我的理解是这样会增加Ls矩阵的秩，从而使得可以求解LsX = DB。

锚点是什么呢？现在可以滑动鼠标看看最开始的gif动态图，图中标识的红色点和绿色点即为锚点。锚点分成两种：固定锚点（形变过程中，该点是雷打不动的，如图中的绿色点），移动锚点（形变过程中，该点是控制移动的，如图中的红色点，这些红色点绕着y轴顺时针旋转，可以说是移动锚点牵引着形变的进行）

有一个博客说的也很有趣：

因为形变需要有控制点这种东西，所以在矩阵中增加锚点信息。也就是说，一个mesh，肯定有点移动了，才发生的laplace形变，不然形变个鬼啊。所以这些移动了的点和形变过程中保证不动的点，成为锚点，因为它们属于我们不需要计算坐标的点，因为我们已经知道了。

添加锚点

直接给出添加锚点的公式：

(Ln×nAa×n)x=(Bn×3Aba×3)\left(\begin{array}{c}{L_{n \times n}} \\ { A_{a \times n}}\end{array}\right)\mathbf{x}=\left(\begin{array}{c}{B_{n \times 3}} \\ { Ab_{a \times 3}}\end{array}\right)(Ln×nAa×n)x=(Bn×3Aba×3)

其中L和B为已知的，A和Ab为矩阵需要增加的锚点信息，Aa×3{ A_{a \times 3}}Aa×3是一个标识矩阵，其中的a表示添加的锚点个数（固定锚点和移动锚点个数之和），n为网格总点数。Aba×3{ Ab_{a \times 3}}Aba×3 表示锚点的笛卡坐标。

这样看上去非常抽象，举个栗子：

设n = 5，a = 2（一个移动锚点，一个固定锚点），固定锚点的坐标为（Fx,Fy,Fz），下标为1, 移动锚点的坐标为（Mx,My,Mz）下标为4。则：

Aa×3=(0100000001){ A_{a \times 3}} = ( \begin{matrix} 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 \end{matrix}) Aa×3=(0010000001) Aba×3=(FxFyFzMxMyMz){ Ab_{a \times 3}} = ( \begin{matrix} Fx & Fy & Fz \\ Mx & My & Mz \end{matrix}) Aba×3=(FxMxFyMyFzMz)

求解形变后的坐标

令L∗=(Ln×nAa×n),Ls∗=DL∗L^{*} = \left(\begin{array}{c}{L_{n \times n}} \\ { A_{a \times n}}\end{array}\right) , L_{s}^{*} = DL^{*} L∗=(Ln×nAa×n),Ls∗=DL∗

则最终求解只需求解：Ls−X=DBL_{s}^{-}X = DBLs−X=DB ,注意这时Ls−L_{s}^{-}Ls−不是方阵，需要两边同乘其转置，将X左边部分转为方阵。

(Ls∗)TLs∗X=(Ls∗)TDB(L_{s}^{*})^{T}L_{s}^{*} X = (L_{s}^{*})^{T}DB(Ls∗)TLs∗X=(Ls∗)TDB

最终求得X：

X=((Ls∗)TLs∗)−1(Ls∗)TDBX = ((L_{s}^{*})^{T}L_{s}^{*} )^{-1}(L_{s}^{*})^{T}DBX=((Ls∗)TLs∗)−1(Ls∗)TDB

接下来就只需要更新点的坐标即可，注意只改变非锚点的点的坐标。

写在最后：
该方法经过测试对于点数少于100的模型具有较好的实时性，对于上千的点的模型程序处理就会变得很慢。