一种基于局部平均法向变形的网格参数化方法

焦 冲, 苏科华, 吴博文, 任术波, 辛 宁

(1. 武汉大学 计算机学院, 武汉 430072; 2. 中国空间技术研究院 通信与导航卫星总体部, 北京100094)

0 引 言

三维网格模型在工业制造、 辅助医疗、 3D打印、 文化娱乐和建筑设计等领域应用广泛. 由于物体的几何形状通常较复杂, 而这些复杂的几何形状使后续的网格处理很困难, 因此, 通常需要将三维网格参数化到一个简单的参数域, 例如平面参数域、 球面参数域等. 作为计算机图形学、 计算机辅助几何设计和数字几何处理中的一个重要工具, 网格参数化在纹理映射[1]、 网格变形[2-3]、 形状建模[4]及网格优化[5-6]等网格处理中都具有重要作用.

网格参数化问题一般可描述为: 给定一个二维流形的三维网格和参数域, 寻找从参数域的点到三维网格点的一一映射, 使参数域的网格与原网格拓扑同构, 且保证三角形之间互不重叠. 理想状态下, 三维网格到参数域之间的参数化映射是等距的, 即原始网格的边长与夹角在映射后应保持不变. 但除可展曲面外, 一般曲面都达不到这一理想条件. 一种参数化方法只能尽可能保持角度或者面积不产生扭曲, 但不能同时消除角度扭曲和面积扭曲. 因此, 为提高算法的适用性, 参数化的核心任务之一就是尽可能降低某种类型的扭曲, 例如保角参数化能最小化三角网格的角度扭曲, 从而有效地保留网格的形状信息.

早期的平面参数化算法一般采用凸组合的方式, 通过固定边界求解线性方程得到平面参数化[7-9]. 这类方法虽然能保证一一映射, 但都具有很大的扭曲. 因此, 一些方法采用自由边界或分割展平的方式直接控制扭曲[10-12]. 优化能量的方法使用无翻转的参数化作为初始, 并优化某种能量以降低扭曲实现参数化[13-17]. 由于扭曲的上界很难预先确定, 因此该类方法需要不断尝试寻找一个最佳上界, 同时由于扭曲度量通常是高度非线性的, 这些基于能量优化的方法计算代价较高.

能量优化方法可应用于球面拓扑的模型, 但对于高曲率的区域目标能量可能会导致翻转[18]. 一种解决方式是通过一个动态调整的参数生成双射的球面参数化[19]; 另一种方式是通过最小化调和能量将0亏格的流形变分为全局共形参数化[20-21]. 如文献[22]使用一些中间的参考三角形定义能量, 并生成具有低等距扭曲的无折叠参数化.

基于几何流的方法成功地实现了在同一框架下实现多种参数域的参数化. 平均曲率流(MCF)[23]是用于演化曲面几何的最基本流之一, 其可等价地表示为最小化曲面嵌入的梯度或最小化曲面面积的流. 但奇点的存在导致MCF计算过程不稳定. Bobenko等[24]提出了Willmore流, 但Willmore流需要曲面自身接近于球面, 并且依赖于高阶导数； Zhao等[25]提出了单位法向流的定义, 并将其应用于网格参数化, 但其计算效率较低并且对于某些复杂的模型效果不佳. 基于Monge-Brenier理论, 最优传输理论(OMT)也可用于计算保面积参数化[26]. 对于poly-annulus曲面, Su等[27]结合OMT和Ricci流计算保面积的参数化. 基于离散Calabi流, Zhao等[28]提出了一种保角的参数化算法. 为提高算法性能, Su等[29]为离散Calabi流引进了边的翻转、 拟牛顿法、 最优步长、 优先级嵌入和边界自定义等操作； Su等[30]利用离散的李导数流计算圆盘和球面的保面积参数化.

基于此, 本文针对单/多边界的开网格以及亏格为0的封闭网格, 提出一种基于局部平均法向变形的网格参数化方法. 首先通过计算局部平均面法向, 以该平均法向为引导, 交替地将三角形旋转到新位置, 然后利用Poisson系统解算旋转拉伸能量方程, 得到分散的三角形顶点新位置, 在视觉上如同将分散的三角形重新缝合. 通过不断迭代, 将顶点不断向邻居位置推移, 使重建网格达到目标曲率. 该方法将局部平均法向的思想引入到MCF的框架中, 并基于该思想用动态控制平均法向权重的方式避免三角形翻转, 极大减少了算法的迭代次数, 提高了求解速度, 解决了传统MCF算法求解不稳定、 迭代速度过慢的问题. 同时, 优化策略也尽量保证了最小化参数化扭曲. Torse模型的原始网格通过迭代变形为最终球面参数的参数化过程如图1所示.

图1 Torse模型的参数化过程Fig.1 Parameterization process of Torse model

1 算法描述

类似于MCF方法, 本文的变形方法通过将顶点移向其邻居顶点实现, 但不直接计算顶点的邻域平均位置, 而是通过面法向量场平均实现顶点位置的计算. 因此, 本文算法共分为两个阶段.

1) 局部平均法向旋转： 计算每个三角形k-阶邻居域的平均法向量, 并以该法向为目标法向, 以笛卡尔坐标系下的原点作为旋转中心, 为每个三角形执行旋转操作；

2) Poisson表面形变： 结合Poisson系统“缝合”网格, 通过优化一种拉伸能量计算顶点的新位置.

交替执行上述两个阶段, 直到网格收敛, 即网格上每个点的平均曲率都相等. 本文用M={V,F}表示三角网格, 其中V={vi,i=1,2,…,Nv}为顶点集合,F={fi,i=1,2,…,Nf}为三角形面集合.算法的实现过程中采用半边结构存储网格信息.

1.1 局部平均法向旋转

本文以k-阶邻居域的平均法向作为目标法向量, 对每个三角形应用旋转变换, 并记录旋转矩阵.

1.1.1 邻居面法向平均

(1)

1.1.2k-阶邻居面选取

图2 不同k值下不同算法的收敛速度对比Fig.2 Comparison of convergence speed of different algorithms under different k values

1.2 Poisson表面形变

通过局部的平均法向分别对每个三角形执行旋转操作会导致整个网格的拓扑发生改变, 顶点被分裂. 因此, 算法基于Poisson系统重建变形的三角网格, 保证每个面都满足当前指定的目标法向量. 由于Poisson系统只能逼近当前目标法向, 采用迭代式的过程可引导面法向不断靠近最终的目标法向. 本质上, 迭代过程可视为是网格表面不断变形的过程, 常见网格变形算法的目的包括两方面: 一是保持原始网格的局部特征; 二是尽可能保持网格的度量, 即保证最小化边长变化. 由于参数化算法需要使扭曲尽可能小, 因此本文在重建过程中应尽量保持网格的度量.

由文献[33]可知, 网格表面形变能量可显式地分为两部分： 拉伸能量和弯曲能量. 拉伸能量与度量相关, 而弯曲能量与曲率相关. 本文算法以拉伸能量为基础, 并以每个三角形的局部平均法向为限制, 建立拉伸能量方程：

(2)

对于每次迭代,Rij可通过三角形的当前法向量和目标法向量计算, 因此Es可视为是只关于顶点坐标的二次函数.通过对Es求导可知, 在梯度为0时, 其解为最优解.通过化简, 可建立关于顶点坐标的一个线性方程组为

(3)

其解为变形后网格顶点的新坐标, 其中vj～vi表示顶点vj与顶点vi相邻接.在线性方程组(3)中, 需已知网格中每个三角形面的旋转矩阵Rji.虽然预先未知迭代变形后顶点的确切位置, 但通过上述计算可知每个三角形在此次迭代过程中的目标法向量, 因此可为每个面计算旋转矩阵.由于式(3)右边为已知量, 因此可将其简写为三维向量bi.对于方程组(3)的系数, 定义为

(4)

因此, 可将方程组(3)简写为

Lv′=b.

(5)

求解方程(5)即可得到网格顶点在本次迭代中的新坐标.迭代操作的终止条件是所有顶点的平均曲率相等, 由于计算机精度的原因, 可将此条件放宽至顶点的平均曲率相差在一个较小的阈值内.由于L是对称正定的, 因此, 可利用Cholesky因式分解求解线性方程组, 比能量优化方法更简单高效.在每次迭代中为每个三角形重新计算旋转矩阵, 用新的旋转矩阵与顶点坐标更新L和b.对于不同拓扑以及规模大小不同的网格, 所需的迭代次数一般不同.对于小型网格, 一般仅需几十次迭代即可达到收敛.

1.3 算法的加速策略

实验发现, 仅使用k-阶邻居平均法向量时, 邻居域范围越大, 顶点移动范围越大, 收敛速度越快.但在某些区域, 尤其是高曲率区域, 过大的邻居域可能会导致三角形翻转.本文通过动态调整邻居域的平均法向量权重, 在保证三角形不翻转的同时加快算法收敛.即在引进全局平均法向量的同时, 设置惩罚函数调整其对目标法向量的影响, 保证三角形不发生翻转.

(6)

由于在形变过程中, 网格的整体形状会发生变化,n*也是不断变化的, 固定权重μ通常无法取得理想结果, 因此本文定义惩罚函数动态地改变μ值.对于平面参数化, 惩罚函数的选择需满足以下条件: 对于面fi, 如果n*与ni之间为钝角, 则应减小权重μ; 由于翻转的非边界点高斯曲率绝对值收敛于2π, 所以在高曲率区域使用n*加速收敛也会导致翻转, 因此当顶点的曲率较大时, 应减小权重μ.基于上述条件, 可设计惩罚函数为

(7)

(8)

对于平面参数域, 本文的基本思想是通过预测每个面的法向方向加速收敛.对于球面参数域, 该策略同样可行, 其惩罚函数需满足以下条件: 对于面fi, 如果n*与ni之间为钝角, 则应减小权重μ; 如果网格M是非凸的, 则应减小权重μ.基于上述条件, 可设计惩罚函数为

(9)

(10)

2 实验与评估

本文以C++实现该算法, 并在AHSP[22]提供的模型上进行测试, 同时与其他算法进行比较. 首先在一个模型集合中对本文的参数化算法进行测试, 包括各种拓扑模型, 验证其有效性和高效性; 然后与目前其他先进的网格参数化方法进行对比; 最后利用具有不同剖分的网格证明其对网格剖分不敏感. 实验平台为配置AMD Ryzen-7-1700处理器, GTX 1060显卡, 16 GB内存, Windows 10操作系统的台式电脑； 图像处理库为OpenMesh 6.3, 矩阵处理库为Eigen 3.3.3; 输入的网格文件为obj格式.

2.1 扭曲度量

为评价参数化效果, 目前已有很多种类参数化扭曲的度量方式. 假设σ1和σ2分别是从原三角形到参数化三角形转换的Jacobi矩阵的最大和最小特征值, 用ρi表示fi的面积与网格总面积之比, 扭曲的度量方式如下：

文献[34]研究表明,Darea和Dangle的值越接近于2, 其面积扭曲或角度扭曲越小.此外, 为进一步比较算法性能, 用等距能量EARAP、 保角能量EASAP和Green-Lagrange能量(EGL)[1]评估参数化效果.

2.2 数据集测试

为验证本文提出的基于局部平均法向变形参数化方法的有效性和鲁棒性, 在2 063个模型上测试该算法, 所有模型均来自文献[22]中数据集.对每个模型的参数化结果, 先计算其扭曲并记录运行时间, 再将所有结果绘制成频率直方图, 如图3所示, 其中max,min,ave和std分别表示对应扭曲或时间的最大值、 最小值、 平均值和标准偏差.由图3可见, 本文算法能高效生成低扭曲的参数化.

图3 2 063个模型的参数化结果直方图Fig.3 Histogram of parameterization results with 2 063 models

对具有多边界的开网格及高曲率模型, 本文方法也能计算出高质量的参数化, 图4和图5分别为这两种类型网格的参数化结果.

图4 具有多边界的汽车模型平面参数化结果Fig.4 Plane parameterization results of car model with multiple boundaries

图5 具有高曲率区域的网格模型及其球面参数化结果Fig.5 Spherical parameterization results on mesh model with high curvature area

2.3 算法对比测试

为验证本文算法的实用性和高效性, 将本文算法与其他参数化方法进行比较, 包括单位法向流(UNF)[25]、 可伸缩局部内射映射(SLIM)[17]和多层次球面参数化(AHSP)[22]. 对于UNF, 遵循其默认参数, 其他算法均采用其源代码, 所有实验均在同一台电脑上运行. 由于篇幅所限, 本文仅选取部分网格的参数化度量结果列于表1和表2. 图6和图7分别为Bear和VaseLion模型的参数化视觉效果.

表1 不同算法平面参数化方法的对比结果

表2 不同算法球面参数化方法的对比结果

由表1、 表2及图6、 图7可见： 基于局部平均法向思想的加速策略减少了迭代时间, 相比其他算法有更高的计算效率； 在平面参数化中, 虽然保面积的效果稍逊于其他方法, 但却拥有较好的保形效果, 并能很好地降低拉伸扭曲; 在球面参数化中, 本文方法尽管在降低面积扭曲方面略差于AHSP方法, 但在其他度量中都取得了较好结果. 因此, 本文方法可在同一个框架下高效地实现高质量的平面参数化和球面参数化.

图6 不同算法对Bear模型平面参数化方法的比较Fig.6 Comparison of different algorithms for Bear model by plane parameterization methods

图7 不同算法对VaseLion模型球面参数化方法的比较Fig.7 Comparison of different algorithms for VaseLion model by spherical parameterization methods

2.4 算法分析

2.4.1 算法效率

通过统计算法在数据集上的运行效率, 可得出以下结论: 算法的收敛速度与网格的规格大小密切相关, 这是因为对于规模较大的网格, 针对每个顶点与三角形, 都需要进行一系列的计算, 因此在每次迭代过程中, 计算规模随着网格规模增大而增加, 总运算时间也会随之延长. 但无论在平面参数化还是在球面参数化中, 本文算法的效率均具有明显优势.

2.4.2 对不同网格剖分的鲁棒性

为测试本文算法对不同网格剖分的敏感性, 对部分网格模型进行重新处理, 例如随机增加顶点数, 随机减少顶点数, 或使顶点分布均匀等. 将本文算法应用于不同剖分的网格模型并统计其运行时间及度量扭曲. 图8为本文算法在不同剖分网格上的测试结果. 通过将棋盘格重新贴回模型中发现, 同一网格不同剖分下的参数化扭曲基本一致; 对比不同剖分模型的运行时间, 即使模型剖分变大, 运行时间却仅有很小的变化. 实验结果表明, 本文算法在不同剖分的网格中都能获得高质量的参数化, 具有较强的鲁棒性.

图8 Buddha模型在4种网格剖分下的参数化扭曲鲁棒性测试Fig.8 Robustness test of parametric distortion under four mesh subdivisions of Buddha model

综上所述, 本文针对具有单/多边界的开网格以及亏格为0的封闭网格设计了一种基于局部平均法向变形的网格参数化方法. 通过交替执行两阶段的操作将顶点不断向邻居域移动, 每次迭代通过求解稀疏线性系统使网格变形, 使网格逐步收敛为一个常平均曲率曲面, 从而生成低扭曲的参数化结果. 一方面, 使用k-阶邻居面平均法向的思想极大加快了算法的运行效率； 另一方面, 通过惩罚函数动态调整每次迭代的目标法向, 有效保证了网格无翻转情况的发生. 通过在2 063个模型上的测试, 验证了本文算法的有效性和鲁棒性. 由于局部平均法向思想仅受曲率限制, 因此本文算法对于三角剖分不敏感, 实用性较强. 此外, 与其他参数化算法相比, 本文算法能高效生成低扭曲的参数化结果.