籍 冉 冉， 郑 晓 朋， 雷 娜， 罗 钟 铉

( 大连理工大学 软件学院, 辽宁 大连 116620 )

0 引 言

数字几何用于描述空间中的几何物体，在工业设计、影视、游戏等领域都有广泛的应用.数字几何首先需要在计算机中对空间几何模型进行描述，从而在此基础上进一步进行建模、分析、仿真等.三角形网格是表述中最常用的一种网格形式，生成技术也比较成熟，但在许多应用领域，四边形网格却具有无可替代的优势.在数值模拟仿真中，四边形网格的稳定性、精度、收敛速度等都优于三角形网格；四边形网格也可以用于生成样条曲面；四边形网格在模型参数化、纹理贴图等领域也有着广泛的应用.在四边形网格生成中保持模型的特征至关重要.

曲面四边形网格生成的方法很多，主要有基于方向场的方法[1-2]、参数化方法[3-4]、“分而治之”的方法[5-6]和基于Morse理论的方法.Dong等[7]基于Morse理论，分析拉普拉斯矩阵的特征向量得到四边形网格，算法首先计算输入网格的拉普拉斯矩阵，求解拉普拉斯矩阵的特征向量，选择合适的特征向量作为特征函数；然后利用该函数生成Morse-Smale复形完成对三角网格的四边单元剖分；最后在此剖分的基础上进行全局分片参数化，并在参数域上进行规则采样获得最终的四边形网格.在之后的工作中，Huang等[8]、Zhang等[9]和Ling等[10]加入了方向、特征对齐和各向异性的约束，通过优化求解带有约束的广义特征问题来求解特征函数.Huang等[8]将谱分析与离散余弦变换的关系应用在特征函数的优化上，建立了特征函数与其对应Morse复形的内在联系，通过改变拉普拉斯算子中的面积项来控制单元密度；通过在广义拉普拉斯特征问题中引入额外的能量项来控制Morse复形边的方向与位置，方向由用户输入的方向场控制，特征对齐对应于网格边的具体位置.Zhang等[9]提出一种基于驻波构造的各向异性四边形网格生成技术，通过挖掘四边形网格和驻波的内在关系，将所有的用户需求同时转化为驻波的各种物理属性约束，通过输入密度场实现了各向异性和密度控制，同样也通过输入方向场来实现方向和特征对齐约束.Ling等[10]通过引入边界条件来实现特征对齐控制，通过求解Helmhöltz方程来得到特征函数，在Helmhöltz方程中加入了密度场的控制.这些算法均需将所有用户控制和质量需求集中在一个全局优化问题中统一求解，算法复杂，随着模型规模的增大，计算代价将急速上升，噪声增多，也需要后期复杂的优化.2018年，Fang等[11]将参数化和Morse理论结合起来，算法通过用户输入标架场来进行参数化，在参数化的过程中也输入特征线来进行特征约束，在参数化发生退化的区域用Morse理论来生成四边形网格，通过将两种方法优势互补，解决了参数化的退化问题和Morse方法在规模较大时运算量大的缺陷.

本文将曲率信息加入到拉普拉斯矩阵中，通过迭代的方法计算特征函数，在迭代的过程中将特征线的信息加入，使得临界点能够准确定位到特征线上，并通过临界点来构造保特征的Morse-Smale复形，进而参数化生成四边形网格.

1 技术路线

基于Morse理论生成满足特征的四边形网格，算法步骤概括为5步，图1是原始算法与本文算法的对比，本文接下来将对每一步进行更加详细的介绍.

(a) 原始算法步骤

(b) 本文算法步骤

2 Morse理论与Morse-Smale复形

2.1 Morse理论

(a) 极小值点

(b) 鞍点

(c) 极大值点

Banchoff[13]用分段线性函数替代光滑函数f将Morse理论推广到三角网格上，f的值定义在三角网格顶点上，通过判断顶点的一环领域内的点来判断临界点.三角网格上光滑的Morse函数通常是求解一个拉普拉斯方程得到，即

(1)

式中：Ni为顶点i的一环领域的所有顶点的集合，wij表示与顶点i相连接的边(i，j)的权重，在三维流形曲面网格上，通过用拉普拉斯矩阵的形式来表示拉普拉斯算子，即

(2)

式中：αij、βij为边(i，j)的两个三角形的两个对角.-Lf=λf，特征函数即为计算特征值所对应的特征向量，特征值0=λ1≤λ2≤…≤λn所对应的特征向量e1，e2，…，en，随着特征值的增大，频率越来越高，临界点的个数越来越多，在模型上的可视化如图3所示.拉普拉斯方程的解称为调和函数，因此拉普拉斯矩阵的特征向量也是离散情形下球面调和函数的基.

2.2 Morse-Smale 复形

Morse函数f，如果它的稳定流形和不稳定流形之间只有横向相交，则f是Morse-Smale函数[12]，稳定流形与不稳定流形相交形成的胞腔称为Morse-Smale复形.Morse函数f定义在三角

网格上，f的值定义在三角网格顶点上，通过比较顶点的值与一环领域内点的值来判断临界点的类型，如图4所示.

Morse-Smale复形是在一个流形上标量函数的胞腔剖分.如图5所示，从鞍点出发沿梯度上升或者下降最快的方向，直到达到极大值点或者极小值点，一般来说会从鞍点引出两条上升梯度线和两条下降梯度线，这些路径把流形分成若干四边形区域，所有的四边形区域都有两个鞍点、一个极大值点、一个极小值点.

Morse-Smale复形的生成直接关乎到最终生成的四边形网格的质量，而生成Morse-Smale复形的关键在于求解Morse函数，在计算Morse函数时通常会充分考虑特征对齐约束的要求，但由于离散化中的精度损失以及噪声点的影响，临界点不能完全准确地定位在特征线上，导致Morse-Smale复形中的边没能与特征对齐，因此需要对临界点进行局部扰动，使其准确定位到特征线上.本文针对这一点提出了相应的改进算法，下面将详细说明.

3 加入特征约束的Morse-Smale复形

生成满足模型特征的网格在曲面四边形化中是非常重要的一个约束要求.而生成满足模型特征的Morse-Smale复形，关键在于求解拉普拉斯矩阵的特征函数，即Morse函数.本文给出如何求解特征函数，使得临界点准确地定位在特征线上.

3.1 加入曲率信息的拉普拉斯矩阵

在计算三角网格上标量函数时，通过拉普拉斯矩阵的某个特征值的特征向量来定义，即

-Lf=λf

(3)

通过借鉴之前的研究[8，14]，对拉普拉斯矩阵的对角线进行扰动.在拉普拉斯方程中加入了模型曲率的信息，使用三角网格曲面模型顶点上的高斯曲率，并归一化：

(4)

其中θi为点v的一环领域的三角形中所对应的顶点为v的角的大小.曲率反映了模型表面的弯曲程度，刻画了模型的几何特征信息，求解特征函数的方程定义如下：

-Lf=λ(I+K)f⟹-(I+K)-1Lf=λf

(5)

其中矩阵K为对角矩阵，对角线上的元素ki等于第i个顶点vi处的高斯曲率.

拉普拉斯矩阵在加入曲率信息后，临界点会因为加入的曲率值而进行相应的位置扰动，临界点能够更加精确地定位到模型特征的地方，生成的Morse-Smale复形更加贴合模型的特征，如图6所示.

3.2 特征函数的迭代算法

通过计算拉普拉斯矩阵L的某个特征值λk的特征向量来作为特征函数f，在这里一般使用的是ARPACK稀疏矩阵特征问题求解器来求解，求解出来的只能是特定的数值λk所对应的特征函数，特征函数的选取不具有灵活性，具体体现在以下两个方面：

(a) 加入曲率信息之前的结果

(b) 加入曲率信息之后的结果

(1)2.1节分析的临界点个数随着特征值的数值增大而增多，如果在求解中λk所对应的特征函数fk的临界点个数偏少，而λk+1所对应的特征函数fk+1的临界点个数偏多，根据原始的方法不能求得fk和fk+1之间的特征函数；

(2)想要对特征值λk所对应的特征函数fk做少许的扰动，或者是对局部顶点上的特征函数做扰动，使得临界点移动到模型特征的地方，根据原始的方法也没有办法实现.

从以上两点可以看出，特征函数的选取受到局限，因此本文提出了一种迭代算法来计算特征函数.

(6)

对于任意的顶点vi，式(6)的第i行可以写为

(7)

这里λ不仅限于特征值，可以是任意数值，则相对应求解出来的f就更具有灵活性，则迭代公式即为式(7).根据拉普拉斯算子的离散形式(1)，定义全局M上的调和能量[15]

(8)

当迭代求解趋于稳定时，调和能量逐渐减小并趋于0.想要对特征值λk所对应的特征函数fk做相应扰动，就将fk作为迭代初值.算法流程如下：

输入： 拉普拉斯矩阵L.

输出： 特征函数f.

Step 1计算矩阵L的特征值λk所对应的特征函数fk，将fk作为迭代的初值f0.

Step 2计算调和能量E0.

Step 3对所有网格上的点进行更新：

Step 4计算调和能量Em+1，当|Em+1-Em|<ε时，停止迭代，否则返回Step3.

3.3 加入特征信息

在模型表面生成四边形网格时，希望生成的网格的边能与模型的特征线对齐，在基于Morse理论生成四边形网格的过程中，Morse-Smale复形的边一定是四边形网格的边，所以生成的Morse-Smale复形的边与模型的特征线对齐.在本文算法中，关键是让Morse-Smale复形的临界点准确定位到模型的特征线上，步骤如下：

(1)首先在模型的表面画出特征线，在本文中，为了验证算法的效果，只选取了部分特征线，特征线的集合为Se，定义这些边的特征属性作为输入，如图7所示

图7 加入特征线信息

(2)在3.2节的迭代算法流程中，针对Step 3，如果[vi，vj]∈Se，μij>1，否则μij=1.针对不同的模型μij的取值不同，在本文的模型中μij∈[1，2].特征函数的迭代公式为

(9)

即增大特征线的权重.通过增大特征线的权重，就可以使得原本在特征线附近的临界点，准确定位到特征线上，并且采用迭代的算法，对不是特征线上的点没有产生太大的影响，如图8给出了增加特征线权重前后的结果对比.

图8 加入特征线的临界点前后对比图

从图中可以看出，临界点都准确地定位到了特征线上，并且其他位置的特征点没有发生太大的变化.从之前的分析也可以知道，Morse-Smale复形的边都是从鞍点出发，连向极大值点或极小值点，因此，只要临界点在特征线上，Morse-Smale复形的边大多就会沿着特征线，少数情况下在特征线上存在两个相邻的临界点为极值点，两个点之间没有复形的边，此时可以通过对偶操作或者补齐操作来优化，补齐操作的优化过程在下一节会详细介绍，优化使得生成的复形的边可以沿着模型的特征线，如图9所示，蓝线是补齐算法补齐的特征线.

图9 沿特征线的复形结果图

4 Morse-Smale复形的优化

4.1 噪声点的消除

初步生成的Morse-Smale复形会产生很多的噪声点，如图10所示.

通过持久性简化[16]来消除不稳定的点对.在生成的Morse-Smale复形中，对于两个相邻的临界点vi、vj，它们对应于特征函数所对应的标量场的一组拓扑特征，持久性值定义为两个点的函数值之差的绝对值，即|f(vi)-f(vj)|，其是描述移除这个拓扑特征时，特征函数的变化量，持久性值较小的点对属于不稳定点对，要进行消除.一对相邻的临界点，必然一个是鞍点，另一个是极值点，定义在它们上的消除操作是指将鞍点和极值点消除，同时将所有连向该极值点的梯度线连向鞍点的另一个同属性的极值点.如图11所示.

(a) 存在大量噪声的结果图

(b) 消除噪声的结果图

图11 持久性简化操作

4.2 几何优化

希望生成的Morse-Smale复形是近似于正方形或者矩形，但是由于根据特征函数计算梯度线过程中的不可控制性，初始复形的形状可能会出现菱形、梯形等，如果不对其优化，会对后续参数化生成四边形网格的质量产生影响，因此对全局的复形进行形状的优化.

(1)文献[14]中提到通过参数化和松弛迭代的方法来重新定位复形的边缘点，文献[8]中由于加入了特征线，为了不使复形的边缘偏离特征线，对文献[14]的方法进行了扩充，固定特征线上的参数化值，借鉴文献[8]的算法，可以对复形的边缘进行优化；

(2)本文迭代算法中，需要对特征线上相邻的两个极值点之间进行补齐操作，由此一个复形区域被分割成了两个三角形，针对因为这种情况生成的三角形区域，可以对这些区域采用“分而治之”的思想，将1个三角形剖分成3个四边形.

5 四边形网格生成

优化完成后的Morse-Smale复形可以保证生成高质量的四边形网格.将每一个复形区域参数化到平面上，这一步可以通过调和映射完成[15]，根据密度要求d，在二维平面上生成规则的d×d维的四边形网格，再将其映回到模型表面，输出四边形网格.

本文给出了几种模型迭代结果，如图12(c)，临界点可以很准确地与特征线对齐，通过补齐操作可以实现最终的保特征四边形网格生成；图12(d)给出了补齐操作的结果.几何优化部分和最后的分区域参数化四边形网格生成都是非常成熟的算法，在本文中没有显示，可以在以后的优化中实现.

图12 4个模型的结果

6 结 语

本文在Morse理论的基础上，提出了一种特征函数的新计算方法.通过迭代算法计算特征函数，使得特征函数的选取更加具有灵活性；在迭代算法中加入特征约束，可以将临界点准确定位在特征线上，生成Morse-Smale复形；最后对初始形成的Morse-Smale复形进行拓扑优化和几何优化，进而生成保特征的四边形网格.本文所展示的结果与文献[8-11]的结果相比，存在不足，但是本文提出了计算特征函数的一种特殊的新算法，并且这个算法实现了临界点对齐这一好的结果，对后续的保特征的四边形网格生成起到了关键作用，在优化完成后也可以生成质量较好的四边形网格，并且本文提出的算法简单，易于实现，输入信息较少，计算效率高.

后续的研究工作包括以下内容：(1)实现对目前四边形网格的优化，并将本文算法应用到实际问题中，验证其有效性；(2)将本文的算法推广到高维，生成三维的Morse-Smale复形，进而生成六面体网格.