【笔记】《Laplacian Surface Editing》的思路

近来在做三维网格编辑相关的工作，于是看了04年的这篇高引用的经典论文，这篇文章在三维中使用拉普拉斯坐标配合多个限制方法实现了效果不错的网格编辑。因为最近太忙了所以现在才抽空写好总结发出来

这篇文章对于数学的要求可能比较高尤其是构造和求解线性方程组那里，需要好好阅读，里面也许还有很多理解不透彻或理解错误的地方，请谅解

这篇文章我发现网络上几乎没有开源的实现，于是我花了一段时间对这篇文章进行实现，目前已经实现了前两个部分，第三个部分可能不会继续实现了，代码等整理好后会开源到Github上到时会通知大家

PartA 总览

交互式的自由曲面变形在以前是一个难以解决的问题, 因为传统的网格自由变形方法会导致表面细节的严重失真. 这篇04年的文章使用拉普拉斯坐标来对网格进行编辑(Mesh editing), 从而能在保留表面细节和结构的情况下对表面进行变形. 除此之外还利用拉普拉斯坐标对涂层迁移(Coating transfer)和曲面块移植(Transplanting surface patches)问题进行了应用.

思路:

1. 将网格坐标编码为拉普拉斯坐标
2. 在拉普拉斯坐标中对网格进行处理
3. 将拉普拉斯坐标解码回网格

PartB 拉普拉斯坐标

这一节是全篇文章的重点

包含了邻近点信息的坐标表示方式已经有很多了, 这被称为差分坐标系, 最常见的应用就是二维上的泊松图像融合. 而之所以要将这样的坐标应用到三维中就是为了找到一种能够在相对坐标中表达出绝对坐标的方法, 这样的表示能够让我们在对网格进行处理时一定程度上忽略掉网格本身的绝对关系, 忽略掉网格在编辑时发生的平移, 旋转, 缩放等情况

拉普拉斯坐标在之前的传统定义中, 虽然有平移不变性, 但是缺少旋转和缩放的不变性. 这篇文章中就通过优化来解决了拉普拉斯坐标的旋转和缩放不变性

首先基本的拉普拉斯坐标的定义如下公式, 下面这个式子有多种不同的权重表示, 文中所采用的是最最基础的直接平均权重(1/di):

在这个式子中Ni代表顶点vi的邻接点, 我们可以想到这是把一个顶点的绝对坐标表示为顶点与其邻接点平均坐标的差值, 这也称之为质心差值. 这样表示的好处就是能够让绝对坐标以一种平移不变性的相对坐标表示出来

指导了一个顶点如何生成拉普拉斯坐标后, 假如我们现在有一整个网格面那么多的顶点, 最好的处理方法就是采用矩阵. 我们可以用 L=I-(D^-1)A 来构造一个拉普拉斯算子矩阵. 这里D是表示每个顶点的度的对角矩阵, A是包括了网格所有顶点的邻接稀疏矩阵. 得到这个表达式后我们就可以通过 Δ=LV 来一口气计算出所有顶点的拉普拉斯坐标了.这一步如果不好理解的话可以自己手动用三个顶点推导一下

得到拉普拉斯坐标Δ后, 我们在拉普拉斯坐标上进行处理, 处理后想要还原坐标回到V时会发现拉普拉斯坐标的秩是n-1, 这意味着我们不能直接解V=(L^-1)Δ这个线性方程组来还原顶点坐标. 解决方法是给这个线性方程组增加条件,固定住其中的一些点融合来求解其他的点. 这些被固定的点就称为控制点handle point

在实践中文章发现如果控制点的目标位置能够满足线性方程组的最小二乘解而不是精确相等的话, 可以得到更好的效果. 因此此时还原顶点的求解就转变为了最小化下面这个能量函数:

在这个能量函数中, 前半部分是要最小化还原的顶点拉普拉斯坐标与变形前的拉普拉斯坐标的差值, 也就是为了保证还原的网格的表面纹理能够保持和变形前一致, 后半部分是最小化新顶点中与控制点有关联的顶点的坐标误差

此时我们开始需要考虑一开始提到的问题: 拉普拉斯坐标有平移不变性但是没有旋转和缩放的不变性, 因为求拉普拉斯坐标时本质上是与邻接顶点的绝对坐标差, 这会带来麻烦. 文章所给出的优化是改变上面的能量函数为下面的式子:

这个式子中最明显的改变就是在变形前的拉普拉斯坐标δi前乘上了Ti, Ti表示的是顶点变换矩阵. 顶点变换矩阵记录了每个顶点vi和其领域在原网格转换为新网格过程中发生的缩放和旋转变换, 是一个图形学中的仿射变换矩阵

要求解这个变换矩阵Ti并不容易, 我们首先想到Ti实际上可以求解下面的能量函数来得到:

这个能量函数的前半部分代表了我们要让原顶点vi经过Ti变换后能接近于新求解出来的坐标, 后半部分表示对于顶点的邻接顶点也要保持有一样的效果

拥有了这样的约束之后, 我们的问题就是如何求解这个变换矩阵Ti, 首先我们将Ti列出来如下式, 这是一个图形学中很经典的仿射变换矩阵, 注意到这里三个轴的缩放参数都是s, 这代表我们只允许小面片变形时发生各向同性缩放:

对于上面这个矩阵, s是缩放比率参数, h是三个轴的旋转参数, t是位移参数, 我们通过求解应用这个矩阵来解决平移, 旋转, 缩放的不变性限制. 回到前面求解Ti的能量函数, 我们可以将这个能量函数转写为下面的矩阵形式来优化运算, 将目标改为最小化下面的这个式子, 其中Ai是原顶点本身和其邻接点的组合, bi就是新求解出来的坐标和其邻接点, 向量(si,hi,ti)其实就是变换矩阵Ti转写的线性组合形式:

进一步写出这个Ai的排列如下式, 同样只要手动用三个点测试一下就会明白这个矩阵乘上向量(si,hi,ti)的效果就和每个点单独与矩阵Ti计算相同:

而bi自然是原顶点的展开形式排列得到的向量:

我们可以采用下面常用的最小二乘法来求解这个线性方程组, 但是在第一次计算的时候我们没有bi的具体信息, 此时我们可以采用常见的迭代求解法, 假定bi为任意一个坐标向量例如最简单的处理方式是选择全1向量处理, 然后对此求解出一个有误差的旋转矩阵Ti, 进行后续计算得到真正的变换顶点v后再迭代求解Ti, 实践中第一次求解尽管有较大误差但是结果已经可以初见雏形, 后续再经过不超过十次的迭代就可以非常精准:

因此根据以上步骤我们可以得到Ti并将其代回求解新顶点vi并迭代, 但是在有些情况下简单求解出来的Ti并不适用. 一种可能就是顶点发生的旋转角度过大, 另一种可能就是目标实际缩放是各向异性的缩放. 对于这两个特殊情况, 解决方法很类似:

1. 旋转角度过大的情况, 在实际进行旋转矩阵计算前先使用作者另一篇文章Differential coordinates for interactive mesh editing的方法对表面提前进行旋转, 这样就可以将大角度的旋转提前转为小角度的旋转
2. 各向异性缩放的情况, 在实际应用前先估算出大致的缩放比率然后对原顶点提前进行缩放从而归一化矩阵所需的缩放尺度

PartC 网格编辑

知道了如何求解拉普拉斯坐标和如何将拉普拉斯坐标还原为普通坐标后, 后续的操作也就简单了起来. 对于三维网格编辑, 我们所需的操作就是先选择感兴趣的变形区域ROI, 得到ROI边界的顶点, 这一方面是为了保证变形不要影响到整个网格区域造成不良的效果, 另一方面是为了减少需要计算的点从而加快计算的速度. 选完ROI后, 我们在网格中选择几个想要的控制点, 然后输入顶点移动到想要的目标位置, 这一步就是控制要编辑的网格需要得到的目标位置

拥有以上数据后, 我们就可以构建线性方程组, 很容易可以想到前面所说的方程组系数矩阵L就是由原顶点, ROI边界顶点和需要改变的控制点组成,其中原顶点满足的是能量函数的前半部分约束, 边界和控制点统称为控制点满足的是后半部分能量函数的约束. 最小化约束就可以还原出绝对坐标也就是重建出网格编辑后的新顶点, 将这些点应用到原网格上就完成了对网格的修改

在实际计算中, 我们会发现构建稀疏矩阵来得到线性方程组的过程运行速度很慢, 如果想要达到文章所说的交互式曲面变形的话我们需要对代码流程进行一些调整. 首先在第一次变形时构建一次方程组, 然后求解得到变形效果, 后续变形中我们只修改方程组最后几项也就是对应控制点的部分, 这样我们就无需重复建立方程组可以大大提升执行效率

PartD 涂层迁移

涂层迁移环节要比网格编辑来得更复杂一点, 但是整体的逻辑也是一样的. 涂层编辑中, 首先我们需要定义出哪些部分是网格的涂层. 网格的涂层, 我们一般理解为网格表面的高频细节, 也可以叫做网格的三维纹理, 如果我们想要得到网格的表面涂层, 这里提到了高频细节, 方法也就呼之欲出了

首先对于需要提取涂层的表面, 我们对表面求出整个表面的拉普拉斯坐标δ, 然后我们对这个表面使用例如拉普拉斯过滤器来进行平滑化, 得到的平滑表面我们也求出拉普拉斯坐标δ_hat, 就这样我们定义这两个坐标相减的差值就是表面的高频细节ξ:

得到高频细节之后, 我们求出要迁移上去的目标表面的拉普拉斯坐标, 然后如果我们将这个高频细节加到目标表面的拉普拉斯坐标上, 我们就可以得到涂层迁移后的拉普拉斯坐标. 对这个新得到的拉普拉斯坐标组合出新的线性方程组然后用PartB的方法进行表面重建就可以得到迁移后的表面

涂层迁移的整个思路都比较简单, 但是实现时却会遇到一些其他的问题, 首先要迁移的目标表面可能本身就存在很多的纹理, 此时如果我们直接应用高频细节的话可能会得到不理想的变形结果, 因此我们可以先将目标表面进行类似的滤波平滑化再处理

然后是当我们对方向不同的表面进行迁移时, 由于方向的改变, 我们需要类似PartB中计算出方向来应用, 但是这次不是算出T来应用在求解中, 而是对原表面的平滑模型的每个顶点计算出额外的旋转矩阵然后应用在目标表面的拉普拉斯坐标中,通过这样的简单对齐来使得拉普拉斯高频细节尽量符合逻辑. 而我们计算旋转矩阵的方法是首先对齐两个对应顶点处的法线, 然后选择两个顶点所连接着的一个对应邻接点, 用邻接点的边投影到顶点的切平面上, 以此作为另一个向量对齐旋转, 这样就对齐了顶点平面的坐标系. 由于拉普拉斯坐标的特性, 我们也可以直接对高频细节差值进行旋转

而一直说到了两个表面的对应顶点, 在实际使用中由于我们要迁移的表面的拓扑结构各不相同, 我们必须对表面进行处理并一一查找两个表面之间的对应关系. 对于拓扑结构的不同, 首先可以用重采样的方式让两个表面的尺度和顶点数目接近, 然后通过参数化的方式得到点与点的对应关系. 文章中使用了流行的均值坐标参数化方法(mean-value coordinate parameterization), 我们为求简单使用普通的平面参数化, 球面参数化等等方法也是可以的. 有些特殊的迁移应用可能会需要和PartC的网格编辑一样要求某些点的坐标对应, 此时需要对参数化方法进行一些额外的调整

而由于参数化的原因, 有时候我们避免不了将大区域和小区域进行对齐, 此时对于中间没有被参数化到的点一般有两种处理方法, 一种是将邻接环映射到区域上然后按照邻接环的关系得到新的映射值, 这样的执行效果比较好但是编写困难, 另一种方法是直接取没有处理到的点的周围几个点来做平均, 这样得到的迁移会有种模糊感但是编写简单. 当然, 参数化最好的情况就是两个区域能尽量一对一配对, 尽量减少参数化中的信息损失

在涂层迁移时, 还有一个常常遇到的问题就是当目标表面也有涂层想要保留时应该怎么办, 文中给出了细节混合策略, 简单说就是将原表面的拉普拉斯高频差值和提取出来的目标表面的拉普拉斯高频差值进行加权求和,通过让两者按比例的混合来达到纹理自然融合的效果

PartE 曲面块移植

曲面块移植问题, 就是想要将一个模型的一部分区域连同整体的几何结构都移植到另一个表面上, 例如将鸟的翅膀连接到兔子的背上, 这个操作的关键一方面是两个部分的摆放问题, 另一方面是两个部分之间的接缝问题, 也就是如何生成合适的三角剖分来让两个部分顺畅地连接在一起

首先我们需要将两个所选区域放到想要连接的位置, 如何我们需要调整到合适的尺寸并在需要的时候进行重采样以让两者的顶点数量接近. 接着我们对边缘部分进行裁剪以将两个区域交叉出界部分的裁去

连接后接缝的处理问题接近于前面的细节融合部分, 此时的接缝由于前面的裁剪所以是断开的, 此时在各自原表面上找到对应接缝的区域, 然后对接缝部分进行参数化, 让接缝部分对两个表面都得到一对一的关系. 然后为了合理地融合起来, 文中按照下面图b所表示的给两个表面都计算了一个尺寸权重, 也就是红蓝线所指示的到接缝处的距离,接着按照这个权重来用类似纹理融合时进行加权插值得到新的顶点位置于是完成接缝处理. 在接缝重建过程中应该定义一些如下图的红色控制点, 它们起到边界的作用保证了区域形状在融合中不发生改变

PartF 总结

拉普拉斯处理这个方法的好处是确实能够对上面三个应用情况得到较好的执行效果, 但也有自己的缺陷. 一个是拉普拉斯的计算需要多次迭代才能得到比较好的执行效果, 另一个是由于拉普拉斯算子矩阵的构造需要邻接矩阵的介入, 因此一旦物体拥有太多的顶点数矩阵的大小就会成倍增大, 会使得内存不堪重负, 还有就是表面参数化这一步骤也会消耗不少的计算时间

PartG 相关工作

这篇文章毕竟已经很老了, 很多文章到了今天参考价值都不大了, 大概总结了相关的工作如下,还是可以用来了解领域扩展视野的:

1. 在参数表面上的编辑操作 HOSCHEK J., LASSER D.: Fundamentals of Computer Aided Geometric Design. A.K. Peters, 1993.
2. 用网格多分辨率方法保留细节下进行编辑, 细节会随着分辨率的提高而出现 GUSKOV I., SWELDENS W., SCHRÖDER P.: Multiresolution signal processing for meshes. In Proceedings of SIGGRAPH 99, pp. 325–334.
3. 多分辨率的差分思想引发了图像泊松融合的出现 PÉREZ P., GANGNET M., BLAKE A.: Poisson image editing. In Proceedings of SIGGRAPH 2003.
4. 三维上的差分坐标也被提出, 是这篇文章的核心来源 ALEXA M.: Differential coordinates for local mesh morphing and deformation. The Visual Computer 19, 2 (2003), 105–114.
5. 通过求解最小二乘解来约束表面的变形方法, 这篇文章的前身和一些细节方法的相关介绍 LIPMAN Y., SORKINE O., COHEN-OR D., LEVIN D., RÖSSL C., SEIDEL H.-P.: Differential coordinates for interactive mesh editing. In Proceedings of Shape Modeling International (2004).