本文只实现最简单Laplace形变，按照论文《Laplacian Mesh Processing》进行展开。
注：内容相当容易理解、本文所有讨论都是基于三角形网格模型（ triangular mesh）进行的
Laplace Deformation中最基础的两个概念就是， global Cartesian coordinates（全局笛卡尔坐标系） 和 differential representations（微分表示）。笛卡尔坐标系不用多讲，就是常见的（X, Y, Z）三个轴表示的坐标信息，而微分坐标（也叫 δ-coordinates）是和 Laplacian operator（Laplace算子） 相关的概念。
首先，M = (V,E,F) 表示一个网格模型，M 代表 Model，V 代表 Vertices，E 代表 Edges，F 代表 Faces，模型上的每一个点，我们都能得到一个世界坐标 
，接下来定义微分坐标 δ-coordinates，
点的微分坐标表示为 
，它表示点的世界坐标和这个点一圈邻接点的世界坐标的加权平均（ the center of mass of its immediate neighbors in the mesh），用公式表示更容易理解：
可以从公式看出来，三个轴上的微分坐标是可以单独计算的。 
就是点 的邻接点，也就是公用一条边的点，
是这个点的邻接点的个数，比如周围有5个点，那么这5个点的坐标都对中心点有1/5 的贡献，微分坐标就是点的坐标和这个中心点的坐标的差（其实就是一个向量）。
这样每个点都这么计算就能计算出每个点的微分坐标了，计算微分坐标的意义在于，Laplace形变认为这个微分坐标能够记录或者说表示模型的局部细节信息，比如鼻子，嘴巴，凸起，凹陷这样的信息，当模型形变后，我们希望模型的局部信息不要被破坏，那么换成坐标的说法就是，模型上点的世界坐标变了，但是我们希望模型上所有点的微分坐标没变，也就是点和点的相对位置关系没变。
一个一个点计算当然可以，但是不好记录和表示，所以Laplace形变用一个Laplace矩阵来记录，我们叫做L，然后用L * V = 
，来得到Laplace坐标，也就是微分坐标。通过上述公式，我们能够得到L：
其中，D矩阵是一个对角矩阵，
，也就是第 i 个点周围有几个邻居。A 矩阵是邻接矩阵（adjacency）：
这样，左乘上点的坐标的矩阵就可以得到所有点的微分坐标，值得注意的是，这里的V矩阵，可以使 n * 3 大小的（n是点的个数，3表示xyz三个维度），也可以单独计算某一个轴的值，也就是 n * 1 大小。
上述矩阵就是Laplace矩阵，但是实际使用过程中，我们常常使用等式两边同时左乘D矩阵的结果，
，也就是：
可以看到 
矩阵是一个对称矩阵。用19个点的模型进行Laplace形变，矩阵的样子如下所示：
这样我们可以得到syz三个分量上的等式，
，
，
。
这个 矩阵被称为 topological (or graph) Laplacian of the mesh，不知道怎么翻译比较好，就直接放英文吧。
论文原文中有专门证明为什么这个向量可以近似表示模型的局部信息，不过和实现没有什么关系这里就不放了。
上述公式中的权重是可以变得，每个点的权重按理说肯定应该不一样的， 最常用的是用cotangent加权，比如下图中 
点的权重可以用两个角的cot的平均 
来计算。数学上应该可以证明，这样的计算方式下，所有邻接点的权重的和还是1。当然也有更多种权重设置方式，比如 
。
这里只不过是权重的计算方式更加复杂，实现的流程上并没有什么区别，所以方便起见本文只描述用统一权重的加权方式来计算的流程。
论文中和网上的教程中一直讲，因为这个矩阵不可逆，不能通过 
来计算坐标，所以需要增加锚点（anchors）来让矩阵可逆。这里我看了两天才明白，但还是不认为是这个原因才需要增加锚点，我觉得就是因为形变需要有控制点这种东西，所以在矩阵中增加锚点信息。也就是说，一个mesh，肯定有点移动了，才发生的laplace形变，不然形变个鬼啊。所以这些移动了的点和形变过程中保证不动的点，成为锚点，因为它们属于我们不需要计算坐标的点，因为我们已经知道了。
这样一下就好解释了，我们知道mesh模型形变前的所有点的坐标，以及形变后的一部分点（锚点）的坐标，通过所有点的Laplace坐标不变这一个限定条件，来计算线性方程组，得到所有点新坐标，对除锚点以外的点的坐标进行更新，就完成了Laplace 形变，就这么简单。网上很多教程讲得非常麻烦还没逻辑，搞得我第一次看都看懵了。
在矩阵中添加锚点的方式是。在等式 左边的 下方添加 anc 行（anc是锚点数量），每一行都是只有一个位置是1（或者非0），其他位置全是0，这个点的索引是几，那个位置上的值就是非0，这个非0值代表这个锚点的权重。
用论文中的例子来说明，下图中的mesh模型，红色点为锚点，第一个矩阵就是没有添加锚点信息的矩阵，第二个矩阵就是添加了锚点信息的矩阵，每个锚点权重都是1。（看懂下边三个图，Laplace Deformation 就完事了）
左边的矩阵添加了 anc 行，右边的矩阵肯定也要加 anc 行。这 anc 的目的就是告诉等式，这些锚点在形变后的坐标是多少，这样能够计算弄出其他点的位移。也就是等式变成了如下形式：
是m行的单位矩阵，m就是上文说道的anc，右边的c就是锚点坐标。
简写一下就是 
（b可以是向量，也可以是(m + n) * 3大小的矩阵）。要求的就是 
（是所有点的坐标，但是我们只需要更新不是锚点的点的坐标）。
补充说明一点，就是锚点的数量的决定是有最低限要求的，一个完全联通的mesh是最少需要一个锚点的，数学上可以证明。不过锚点的数目多，对形变的控制性更强，不会出现很奇怪的情况。
由于左边的 
是m*n的，不可逆，所以需要乘上转置，转成方阵再求逆，也就是 
。
我是用VS2015上对MeshLab2016的源码进行编写的，添加了Laplace形变的功能，使用的是C++，以及VCG库里边的Eigen矩阵库。使用Eigen库很容易实现，用 MatrixXf 声明矩阵，用 VectorXf 声明向量，然后用三个 VectorXf 向量接收三个轴的所有点的坐标即可。10k+个点（共3个锚点）的话9分钟完成，实在是有点慢。
可以看出来我这种实现方式，没有保持原有的体积，或者说这种laplacian deformation是对旋转不可靠的。有一种比较可靠的方式就是上图的Volumetric Laplacian Defroamation。
论文里对这个过程提到了加速的方式。因为数学上可以证明，矩阵 
是一个稀疏的正定矩阵（sparse and positive definite），所以可以计算它的cholesky分解 sparse Cholesky factorization，
，R是上三角或者下三角的系数矩阵，这样计算能够省略很多重复的计算。论文中使用的是 TAUCS 这个库，研究了一下发现和Eigen好像差不多，而且Eigen接口更好用，MeshLab2016也是用的VCG的库，所以直接就用了Eigen库的里边的cholesky分解来求线性方程组。
也就是求 Ax = b 这个线性方程组的时候，使用cholesky分解（LU三角分解法的变形）进行加速。
比如我们要得到的所有点坐标是 Eigen::VectorXf v_new，A 是 MatrixXf 类型的矩阵（就是上文的M），b 就是上文添加了锚点信息的b。那么直接 v_new = M.llt().solve(b); 就可以了。v_new 的大小都不需要手动设定。这种加速方式下，10k个点2分钟多一点完成形变，但是论文中提到的是10s不到就可以结束，可能还是TAUCS库更快？
Eigen库中还有其他的解线性方程组的方法，如下图（来源 https://www.cnblogs.com/python27/p/EigenQuickRef.html ）：
总结来说，对一个mesh模型进行Laplace Deformation的流程就是：
1、计算所有点邻接点的权重；
2、通过权重得到Laplace矩阵
3、通过Laplace矩阵（方阵）计算所有点的Laplace 坐标；
4、选择锚点，记录所有锚点的索引值，以及形变后的坐标值；
5、对Laplace矩阵添加锚点信息，得到矩阵 A。添加的每行都是（0 0 0 0 
0 0 0）这样的， 为这个锚点的权重， 所在的列数为锚点的索引数，比如第3个点（从0开始），就应该第2列为非0，其余位置为0；
6、对右侧坐标的矩阵或向量添加锚点信息，得到矩阵 b。如果是每个轴单独计算，则就是按照上一步添加锚点信息的顺序一次添加，填写的是锚点移动后的坐标，比如移动后为（1, 2, 3），如果指计算x轴，那么在向量后边加上1即可，如果三个轴一起计算，那就再b矩阵下边添加一行（1, 2, 3）即可；
7、用Cholesky分解法加速求解 Ax = b 这个线性方程组；
8、由于记录了锚点的索引值，所以遍历所有点，对不是锚点的点，用上一步得到的坐标进行位置的更新即可。
所以说，比较重要的就是，邻接点的权重的计算方式以及计算，锚点的索引和形变后的坐标值。
这只是简单的Laplace形变，不过更高级的Laplace形变也都是在这个基础上的，区别不大。
接下来还会在写一份完整的代码流程以及结果展示的总结。