摘要
在网格生成软件处理几何模型过程中，针对几何表面信息缺失带来的孔洞问题，提出一种基于B样条曲面的填充修复方法。根据拓扑关系在给定的B样条曲线集中提取孔洞边界，针对单个孔洞包含的曲线采用曲线逼近拟合与组合技术进行预处理得到相容曲线。由曲线构造单向插值直纹面与张量积曲面，然后将曲面通过布尔和操作生成双线性差值B样条曲面来修复孔洞。此外，为保证方法的健壮性，针对复杂的特殊孔洞，可进一步采用直纹面生成填充作为候补方法。实验结果表明，方法具有很好的通用性，能适用于真实工业数模中各类形态孔洞的脏几何修复，为后续的网格生成提供干净、封闭的几何模型。
Abstract
In the process of geometric model processing by mesh generation software, a filling method based on B-spline surface was proposed to solve the hole problem caused by missing geometric surface information. The hole boundaries were extracted from the given set of B-spline curves based on their topological relationship, and curve approximation fitting and combination techniques were employed to process the curves within individual holes to obtain compatible curves. Hole repair was achieved by first constructing unidirectional interpolating ruled surfaces and tensor product surfaces from the curves, these surfaces were then combined through interpolation and Boolean sum operations to generate bilinear difference B-spline surfaces for filling the holes. In addition, the ruled surface was applied as an alternate method in complex special holes to ensure the overall robustness of the method. Experimental results show that the method is highly general and can be applied to the dirty geometry repair of various types of morphological holes in real-world industrial geometric models, providing clean and closed geometric models for subsequent mesh generation.
Keywords
几何数模作为开展科学工程计算的原始对象，在前处理流程一般需要通过网格生成将其离散表示为计算网格。三维几何数据的完整性是网格生成最基本的要求，同时对网格生成质量及后续科学工程计算[1]的结果有重要影响。由于在设计制作时本身的数据缺陷和在不同标准之间转化数据时可能导致几何信息丢失等，用于数值模拟的几何模型大多情况下成了存在“噪声”的“脏几何”[2]，无法满足后续网格生成要求。此外，几何模型在设计过程中通常会包含完整的局部设计特征细节，而这些微小特征往往不会影响数值模拟的结果，但会显著增加网格生成算法处理的难度，导致局部非关键区域进行加密，扩大网格总体规模，浪费计算资源，甚至会造成网格质量降低从而影响后续科学工程计算的精度与效率。科学工程计算在网格生成前针对几何模型的处理主要有几何修复和特征简化两种方式：几何修复是在保证模型不失真的条件下对“脏几何”中瑕疵部分进行修复；特征简化主要是消除不影响分析结果的设计细节，得到更易于网格生成以及科学计算的模型。由于引起几何孔洞瑕疵的原因多样，无论是模型本身的原因还是由于特征简化所引起，针对孔洞的处理贯穿于几何处理整个流程。
自由曲面孔洞的修复是计算机辅助设计（computer aided design，CAD）常见问题，解决填充方案主要分为：约束求解裁剪曲面[3]、细分曲面[4-5]、插值边界生成曲面[6-8]和构造中心点多边域修复[9-11]。其中，后面两种方式采用工程应用中大量使用的非均匀有理B样条（nonuniform rational B-spline，NURBS）曲面，相比细分曲面和裁剪曲面更易于数据交换。Shi等[6]利用退化参数曲面的性质，通过插值解决具有不相容边界四边区域曲面生成的方法，但无法保证在不兼容角点的连续性；Mu[8]通过插值边界曲线及其边界导数，并采用能量优化控制点逼近给定曲面点构造B样条曲线曲面，但边界曲线数量有严格的控制；Piegl和Tiller[9]提出的方法对边界曲线的数量没有限制，并且可以独立指定跨边界导数，同时也要求切向相关条件；Yang等[10]对Piegl和Tiller文章中的方法在NURBS曲线曲面应用上进行了扩展，同时跨边界导数由边界的NURBS曲线插值近似，使该方法可以完成更为复杂的孔洞修补。综上所述，相关研究虽能给出比较好的效果但一般会伴随特殊的条件。
面向国产高性能平台[12]上实际科学工程计算应用中的前处理流程，针对几何模型处理中常见的孔洞问题，本文提出一种基于孔洞边界几何与拓扑数据生成B样条曲面的通用修复方法。
1 孔洞边界处理
1.1 B样条曲线曲面定义
待修复的孔洞边界由B样条曲线构成，采用控制顶点定义的曲线方程可写为：
(1)
式中：di（i=0，1，···，n）为控制顶点，其按顺序连接的折线为B样条控制多边形；Ni，k为k次B样条基函数，由被称之为节点矢量的非递减序列U0≤U1≤···≤Un+k+1所决定的k次分段多项式。
B样条曲面为曲线从二维到三维的推广，其标准方程为：
(2)
式中，di，j为（m+1）×（n+1）大小的控制点阵，k和l分别为u和v方向上的次数，Ni，k（u）与Nj，l（v）由u、v方向节点矢量U=[u0，u1，···，um+k+1]与V=[v0，v1，···，vn+l+1]决定。
1.2 提取孔洞边界
面对复杂多样孔洞空间形态，设计了一套通用的处理方法，具体过程如图1所示，主要包括提取孔洞边界、双线性插值B样条填充以及直纹面填充。
图1B样条曲面填充孔洞流程
Fig.1Completion process of filling hole with B-spline surface
针对孔洞的处理目前在工业界还无法做到全自动精准识别并修复孔洞，对于存在孔洞瑕疵的几何模型，往往需要借助用户交互等辅助手段[13]，在待修复的几何模型中标识出包含孔洞边界的曲线集。提取孔洞边界的功能是在标识出的曲线集合中筛选出封闭环，即B样条曲线组成的孔洞边界。如存在多个孔洞边界，每一个孔洞均单独按照后文介绍的双线性插值曲面或直纹面填充算法进行修复。
提取孔洞边界具体实现的主要步骤如下：
1）获取所有曲线的端点拓扑编号或三维坐标；
2）循环对比编号或坐标确定曲线之间连接关系，并按顺序进行排序，如形成闭环即确定一条完整的孔洞边界；
3）对于无法组成闭环的曲线，根据需求舍弃或者在距离最短且无连接关系的曲线端点之间创建B样条直线直至构成闭环。
2 双线性插值B样条曲面填充算法
该算法的主要核心思想是将两对相同次数、相同节点矢量的B样条曲线，即相容曲线通过双线性插值生成B样条曲面填充孔洞，每个曲面生成包含至少一对、最多两对相容曲线。
2.1 孔洞边界曲线预处理
提取出的孔洞边界所包含的曲线数量不定，内部曲线的次数与节点矢量也不尽相同。边界通过预处理生成相容曲线满足曲面生成条件。预处理包括拟合逼近曲线与组合曲线。
2.1.1 拟合逼近曲线
为了统一曲线的次数及节点矢量，在保证误差范围的前提下对边界曲线进行拟合，用相容的逼近曲线代替原始曲线进行曲面生成。B样条曲线逼近需要确定的未知量有三个：曲线次数、节点矢量、控制点坐标。
高次曲线往往会带来更多难给出的边界问题，在广泛实践中往往采用三次B样条曲线[14]。采用三次曲线进行拟合。
节点矢量根据控制多边形与B样条曲线的关系进行计算。控制多边形为特殊样条曲线，通过细化本身节点矢量使其与控制的原始样条曲线相近似[15]。通过创建临时B样条曲线，以曲线与原曲线之间若干个等参点的空间距离为误差控制参数，采用迭代对节点矢量逐步细化。临时B样条曲线的控制点位于原始曲线上，空间坐标由当前迭代步的节点矢量ti确定。控制点对应的参数坐标序列由Piegl和Tiller[16]所提到的节点平均来确定：
(3)
式中，k为曲线次数，m+1为临时B样条曲线控制点数量。
由误差控制参数计算出曲线之间的平均误差距离与最大误差距离，根据两者比重不同，对当前迭代的节点矢量进行全域或者最大误差区间细化。对于区间[ti，ti+1]，通过在区间1/4处和3/4处进行线性插值完成细化。根据端点插值条件与定义域要求，一般采用定义域两端重复度为k+1的端点条件，即定义初始节点矢量ti为t0=···=tk=0，tk+1=···=t2k+1=1。将初始节点矢量通过上述过程进行迭代，直至最大距离满足要求或者节点矢量序列长度达到允许的最大数量。此外，曲线拟合目的是同步两条曲线的节点矢量时，拟合曲线的节点矢量由两条曲线的自身的节点矢量取并集，并加入端点条件确定。
曲线中控制点di（i=0，1，···，m）采用Shene[17]介绍的一套标准的最小二乘算法加上端点约束来进行计算。对原曲线上的有序数据点D0，D1，···，Dh，逼近B样条曲线c（u′）满足端点约束D0=c（0）=d0，Dh=c（1）=dm，其余控制点由最小二乘算法所确定的函数方程确定，即：
(4)
其中，u′i为数据点Di在原曲线p（u）上对应的参数坐标，Di=p（u′i），由节点矢量ti在每个参数区间按曲线次数k均匀采样。
为了保证拟合精度，需要限定逼近曲线与原始曲线在某个误差界范围内。选取Hausdorff 距离[18]作为误差指标，通过增加曲线控制点数量进行迭代，直至满足两条曲线的误差小于常数容差ε或逼近曲线控制点数量达到限定数量的最大值，常数容差ε一般可设置为0.8。
2.1.2 同次数曲线组合
双线性插值曲面生成对孔洞边界所包含曲线的数量有严格限制，需要将边界中多条按顺序相接的独立曲线用统一的B样条表示以减少曲线数量。
曲线的组合顺序是所提方法的技术关键点。以曲线连接端点左右切矢量夹角的余弦值大小作为组合判断条件，按曲线弧长顺序，对余弦值大于给定临界值的曲线进行组合。如组合后曲线数量仍大于四条，则利用四个最小余弦值所属的连接端点作为分割点将孔洞曲线分为四组，将每一组中的曲线按相接顺序组合曲线即可。
组合前的曲线利用逼近的方式统一了曲线的次数k，即组合后曲线为三次B样条。
各曲线在组合曲线整体参数域中的定义区间根据自身弧长Li与累加弧长之比，按照连接顺序划分曲线占有的相应区间：
(5)
然后将各曲线局部参数域uj∈[0，1]变换至各自的区间Uj∈[Ui-1，Ui]上，Uj=Ui-1+uj（Ui-Ui-1）。在公共连接点的节点重复度取k，两端点节点重复度取k+1，得到组合后曲线节点矢量。
组合曲线控制顶点di由各曲线控制点按式（6）确定：
(6)
其中：mi+1为组合前各曲线控制点的数量，m+1为组合后曲线控制点数量，；d1i与dlj为按曲线连接关系排序后第1条与第l条曲线控制点。
2.2 双线性插值曲面生成
以相容曲线构造u方向与v方向单向插值直纹面与张量积曲面，通过对三张曲面进行布尔和运算得到双线性插值曲面。
以v方向边界曲线为准线、对应参数点连线的直线族为母线构造u向3×3次直纹面q（u，v）。其控制顶点通过相容边界线对应控制顶点d0，j和dm， j插值计算，（i=0，1，···，m; j=0，1，···，n）。其中，系数在u向边界节点矢量区间U=[u0，u1，···，um+k]上通过式（1）计算，m、n分别为u、v方向边界曲线的控制点数量。
同理，在v方向边界之间构造直纹面r（u，v），对应控制顶点为。
根据张量积曲面的性质，利用角点d0，0、dm，0、d0，n、dm，n生成张量积曲面，并对曲面进行升阶，构成3×3次双向插值曲面s（u，v），对应控制顶点为。
如图2所示，上述三个曲面的布尔和p（u，v）=q（u，v）+r（u，v）-s（u，v）就是所求的双线性插值曲面。曲面的次数和节点矢量与边界曲线一致，对应的控制顶点：
(7)
特别地，孔洞边界所包含的曲线数量为两条或三条时需要通过不同程度的特殊处理。其中由一对相容曲线构造的单向插值直纹面即可作为两条孔洞曲线填充的曲面；而三条曲线的情况下可将其中的一个顶点当作一条退化的曲线，以点与对边、点连接的临边作为两对相容曲线完成孔洞域内生成曲面填充。此外，孔洞曲线为一条闭环曲线时会采用下面介绍的直纹面填充修复方法。
图2B样条曲面布尔和运算
Fig.2Boolean sum operations of B-spline surface
3 直纹面填充孔洞方法
在实际应用中，由边界曲线构成的孔洞可能形成复杂的凹多边形，此时生成曲面会存在扭曲折叠或者无法沿边界生成贴合的曲面，导致修复效果不佳，需要采取鲁棒性更好的方法进行修复填充，此时则可采用直纹面填充修复方法进行处理。
曲面中u、v两方向上的等参线中有一族为直线的曲面为直纹面，该类型曲面能很好地保持孔洞边界特性，且易于后续网格生成，可用于填充较为复杂的凹形孔洞，作为插值曲面生成失败的补救方案。实现过程以构造中心点填充多边区域为基本思路，先确定边界孔洞的中心点，然后以中心点与边界划分区域，每个区域内生成直纹面作为孔洞填充曲面。
为了保证曲面整体性、减少区域划分、尽量消除对后续网格生成和科学工程计算的影响，需要将孔洞边界曲线进行组合。与前文同次数曲线组合中介绍的方式区别点在于：
1）保证组合前曲线定义域两端节点的重复度为ki+1。
2）为了避免组合前曲线节点矢量中存在过小区间，以及便于后续节点矢量组合处理，需要重新参数化组合前各曲线：针对单条曲线，以曲线弧长L为比例系数，对节点矢量等比缩放，使节点矢量中最大值为max（Li）。
3）需要确定组合曲线次数k，其中k=max（k1，k2，···，kn），并将次数低于k次曲线升阶到k次。升阶原理是通过改变节点矢量的节点重复度将k′+1次基函数Ni，k′+1用k′次基函数Ni，k′线性组合进行表达，然后通过原曲线控制点平均得到升阶后曲线的控制点坐标，通过迭代完成从ki到k的升阶，具体过程见文献[19]。
孔洞边界曲线经过曲线组合操作后，若边界曲线已合并成为一条完整封闭曲线，如图3（a）所示，则以中心点与合并后曲线连线的直线族为母线创建单张直纹面修复孔洞。若组合过程因曲线控制点数量超过限制等原因造成曲线合并失败，则按边界数量与中心点划分多个子区域，在子区域内生成直纹面共同修复整体孔洞，效果如图3（b）所示。
图3直纹面修复孔洞的两种不同方式
Fig.3Two different ways to repair holes in ruled surfaces
中心点坐标O通过各曲线pi（u）上的坐标点确定：
(8)
其中，N为孔洞边界曲线数量。
最后，以中心点与区域内边界曲线上对应点连线为母线族构成直纹面修复当前区域。
4 算法复杂度分析与应用实例
4.1 算法复杂度分析
由数模孔洞填充修复方法的算法流程可知，其计算复杂度包括曲线降阶拟合、拟合曲线的接收、曲线组合、双线性插值曲面创建、直纹面创建5个方面的计算。令孔洞边界的曲线数量为N，则有：
1）曲线降阶拟合：对孔洞每一条边界曲线降阶拟合，最好情况是一次拟合满足要求，最坏情况是直到拟合曲线控制点数量达到限制数量，复杂度为O（N）。
2）拟合曲线的接收：每次曲线拟合都需利用Hausdorff距离进行接收判断，复杂度为O（N）。
3）曲线组合：边界曲线数量大于4时需要执行边界组合，执行次数为N-4，复杂度为O（N）。
4）双线性插值曲面创建：只需创建一次，且由4条曲面创建，复杂度为O（1）。
5）直纹面创建：由中心点与边界曲线共同组合创建N片直纹面，复杂度为O（N）。
综合上述分析可知，算法的整体计算复杂度为O（N）。
4.2 应用实例
图4展示了在几何修复中，孔洞自身所包含边界为相容曲线的填充应用实例。图4（a）是待修复孔洞的模型，图中红色边线是存在几何缺失的孔洞边界，左侧曲面缺失区域边界以直线构成，右侧区域边界曲线为3条三次B样条曲线，且其中有一对曲线具有相容特性。图4（b）是填充修复后的模型，修复后曲面具有较好的逼近程度，且与模型相邻曲面在边界上G0连续。
图5展示了边界曲线孔洞多曲线构成修复实例。图5（a）中给出的是模型整体孔洞修补情况及其局部细节，其中导入模型包含两块缺失区域，一块区域的孔洞边界包含14条曲线，孔洞形态复杂，另一块区域曲线数量相对较少但曲线之间连接点位置角度变化大。如图5（a）右图所示，两块区域修补结果均体现了孔洞边界几何特征，与周围曲面连接紧密，具有良好的光顺性。但在修复过程中，拟合逼近过程会使得曲面边界与原始曲线存在偏差，且多曲线组合会加剧偏差现象，最大误差距离发生在图5（a）红框处，大小约为0.7。模型整体包围盒大小约为3 500×1 000×1 500，在最大网格尺寸为10的条件下生成的网格如图5（b）左图所示，整体约包含20万个网格单元面。以纵横比作为非结构曲面网格的质量评定标准，其对应的三角形曲面网格质量定义公式为：
图4导弹模型孔洞修复
Fig.4Hole repair of missile model
(9)
式中，R为网格单元外接圆半径，r为单元内接圆半径。QAR∈[0，1]，当其等于0时代表网格单元面积退化为0；等于1时表示网格单元质量最好，为正三角形。参考图5（a）右图所示的局部对应的网格细节如图5（b）右图，误差对网格生成无明显影响，区域内网格面质量系数均达到0.5及以上，满足计算要求。
图5汽车模型孔洞修复和网格生成
Fig.5Hole repair and mesh generation of car model
图6和图7均展示了孔洞修复在特征修复中的应用。如图6（a）所示，水杯模型整体包围盒大小约为1 000×1 000×700，图中蓝色曲面自身不属于小曲面，但在红框标识位置包含非常狭长的区域，该处最窄宽度大小为1。以最大网格尺寸为3的条件下生成的网格，网格整体数量约为36万，经网格质量检查，存在30个低于0.3质量系数的网格单元面，均在蓝色曲面最窄宽度处。利用曲线切分、投影等曲线处理算法以及孔洞填充方法在保持原模型几何特征的基础上对握把处曲面进行重新划分，如图6（b）所示。将简化后模型与原始模型在同一网格尺寸控制下进行网格划分，两者网格数量基本一致；同时相对原始模型，简化后模型中0.9质量系数的网格单元面数量提升了2%，并消除了0.4质量系数及其以下的网格单元面。
图6孔洞修复应用在处理设计缺陷的实例
Fig.6Example of hole repair application in dealing with design defects
图7（a）为原始几何模型，图7（b）是经过消除孔洞、凹槽等设计特征后通过直纹面填充的模型，图7（c）、（d）分别是原始模型与修复后模型在相同网格控制尺寸下的网格生成结果。原始模型中在特征小孔洞附近发生局部非必要加密，导致网格数量显著增加，最终网格数量达到65 244个单元，而简化特征后网格仅包含24 953个单元。
图7孔洞修复在特征简化中的应用实例
Fig.7Application example of hole repair in feature simplification
5 结论
针对用于网格生成的几何模型修复过程中遇到的孔洞修补问题，提出了一种有效的解决方法。该方法具有良好的通用性与鲁棒性，能广泛地应用于整个几何模型修复处理流程。孔洞边界的预处理是通过拟合逼近与组合，使得边界曲线满足后续曲面填充修复要求，并保证在逼近过程产生距离偏差不影响后续网格生成精度；基于布尔和操作生成双线性插值曲面，该方式可以保证整张曲面光滑，边界信息明确。候选的直纹面方法无法保证填充曲面整体光滑，但对于复杂凹形轮廓边界孔洞具有良好的适用性，能够较好地保证方法整体的健壮性。
