基于几何矩的CAD模型形状匹配算法及应用

冷珏琳，张 哲，刘田田，郑 澎

(1.中物院高性能数值模拟软件中心，北京 100088；2.北京应用物理与计算数学研究所，北京 100088；3.中国工程物理研究院计算机应用研究所，四川 绵阳 621900)

形状匹配是计算机图形学、计算机视觉、模式识别等领域的一个重要研究课题。形状匹配通常分为特征提取和特征匹配2 个步骤[1]，即根据形状特征的相似程度来确定形状的匹配程度。现有的形状特征提取方法主要分为4 类[2]：基于几何结构的形状特征提取、基于拓扑关系的形状特征提取、基于函数投影的形状特征提取和基于统计特征的形状特征提取。目前还没有一种通用的方法能够对各类物体的形状特征进行描述。根据应用场景的不同，选用的方法也不尽相同。例如，基于几何结构的特征提取方法适用于严格定义的，具有一定刚度或形变较小的模型；基于拓扑关系的特征提取方法适用于描述有关节或分支的模型，且对噪声的敏感度高；基于函数投影的特征提取方法通常要求物体是封闭的；基于统计特性的特征提取方法计算简单，但对形状的描述不够充分，比较适用于粗粒度的匹配。根据已有的经验总结[3-4]，对于形状特征的描述应该具备的优点包括：信息的描述和分辨能力强、计算速度快、易于存储和查找、具有几何变换无关性、对噪声具有鲁棒性、描述方法独立于具体的物体对象等。

针对CAD模型的几何体形状匹配问题，国内外学者已经提出了大量的形状匹配算法[5-8]以提取和匹配CAD模型的形状特征，常用的方法包括基于骨架提取和基于表面的匹配算法等。基于骨架提取的匹配算法的主要思想是通过提取几何体的曲线骨架来描述三维形状特征，且对骨架图进行匹配[5-6]。骨架的提取通常是这类算法的难点，不仅运算量较大，而且对噪声异常敏感。基于表面的匹配算法采用表面的曲率、法向等特征的分布对三维形状进行描述[7]。例如，将曲面上各采样点的2 个主曲率的统计直方图作为三维形状的特征描述，通过比对曲率分布情况来匹配三维体。这类匹配算法的准确性需要通过足够的采样密度来保证。

对于CAD模型在刚体变换下的形状匹配问题，最理想的方式是找到独立于平移、旋转和缩放的形状特征不变量，然后根据特征不变量来定义物体形状的相似度。采用满足相似变换不变性的几何矩作为形状特征的匹配算法[9-14]已经在图像及三维形状检索中得到了广泛地应用。本文针对刚体变换下的CAD模型形状匹配问题，提出一个基于几何矩的CAD模型形状匹配算法。

面向CAD模型的形状匹配算法有着非常广泛地应用。如，用户在CAE 软件中对CAD模型设置材料、计算区域和边界条件等属性时，需要批量选取几何外形相似的几何体。采用自动化的几何体相似性匹配算法进行筛选，能够减轻用户逐个选取几何体的工作量，同时可减少手工操作时错选、漏选情况的发生。本文提出的CAD模型形状匹配算法将应用于CAE 软件的相似几何体拾取，算法的有效性将在实际应用中得到验证。

1 几何矩方法及相关理论

1.1 几何矩

矩在力学中用于表征物质的空间分布。第一篇关于矩方法的研究发表于1962 年[9]，其提出了基于二维情况下直角坐标系的几何矩的概念，并在理论上证明了每个物体对象都可由无穷阶的几何矩唯一确定。之后，矩的定义被推广至三维[10]，并给出了矩不变量的定义，即矩不变量是一种统计特征，满足平移、旋转和缩放不变性。

矩是形状密度函数在核函数下的积分。通过选取不同的核函数可定义不同类型的矩，如几何矩[9-12]、径向矩[13]、Fourier-Mellin 矩[14]等。

几何矩的核函数为直角坐标系下的基本集

其中，p1+p2+………+pn=p。几何矩实质上是形状密度函数的极数展开式的系数。由几何矩的存在唯一性定理可知[9]，几何矩序列可唯一确定形状密度函数。

1.2 三维几何体的几何矩

设Ω为三维几何体所在区域，且几何体的密度是均匀分布的(当(x,y,z)∈Ω时，f(x,y,z)＝1，(x,y,z)∉Ω时，f(x,y,z)＝0)，则该几何体的几何矩为

特别地，三维几何体的体积可由0 阶几何矩m000表示，重心坐标可由0 阶几何矩和1阶几何矩决定，即

1.3 三维几何体的几何矩不变量

如果一个几何体与目标几何体之间只存在位置、朝向、大小的区别，则称该几何体与目标几何体之间存在相似变换(或刚体变换)。由几何矩演变而来的具有平移、旋转、缩放不变性的不变量被统称为几何矩不变量。下面给出本文在第2 节算法中采用的几何矩不变量的定义[9]。

1.3.1 中心矩

基于重心坐标式(3)，可定义三维几何体的中心矩

中心矩具有平移不变的性质。

1.3.2 缩放不变矩

1.3.3 旋转不变矩

基于旋转不变核函数定义的矩，称为旋转不变矩。文献[9-11]中定义了一系列常用的旋转不变矩。本文将采用3 个旋转不变矩，即

其中，D(O,i)，A(O,i,j)均为旋转不变核函数。D(O,i)为坐标点pi=(xi,yi,zi)到坐标原点O的距离；A(O,i,j)为坐标点pi、坐标点pj与坐标原点组成的三角形的面积。

类似地，可以定义平移缩放旋转不变矩为

2 基于几何矩的CAD模型形状匹配算法

基于几何矩方法及其相关理论，本文针对CAD模型的形状匹配问题，提出一个基于几何矩的形状匹配算法，用于识别CAD模型中具有相似形状特征的几何体。算法采用几何矩不变量来描述几何体的形状特征，然后根据几何矩不变量的相似程度来确定几何体形状的相似度。

2.1 算法整体流程

基于几何矩的CAD模型形状匹配算法的主要步骤为：

输入：CAD模型，目标几何体，相似度阈值T。

输出：CAD模型中所有与目标几何体形状匹配的几何体。

步骤1.导入CAD模型。

步骤2.计算CAD模型中每个三维几何体的一组几何矩不变量

步骤3.计算CAD模型中各几何体的形状特征向量与目标几何体的形状特征向量的相似度，根据给定的相似度阈值T和2.3 节给出的判定准则，筛选出所有与目标几何体形状匹配的几何体，并返回结果。

在步骤2 中，本文选择了13 个低阶几何矩不变量组成几何体的形状特征向量。对于具有良好光滑性的CAD模型，采用低阶几何矩不变量就能很好地将几何形状存在差异的几何体进行区分。与高阶矩相比，低阶矩的计算量更小，也更稳定，具有计算简单、易存储、易查找等优点。关于几何矩的计算方法将在2.2 节进行介绍。

在步骤3 中，需根据几何体形状特征向量的相似程度判断几何体之间是否匹配。其中，相似度阈值T用于界定2 个几何体是否相似，取值范围为0～1。给定的阈值T越接近于1，则要求形状的匹配度越高。相似度判定准则将在2.3 节中详细介绍。

2.2 三维几何体的几何矩计算

CAD模型普遍采用边界表示法(B-Rep)表示，几乎所有主流格式的CAD模型都可转化为三角面片网格(STL)表示。因此，本文的算法采用CAD模型的三角面片网格来计算几何矩，具有较强的通用性。

2.2.1 三维几何矩的快速计算

形状匹配算法的运算量主要集中在几何矩的计算。为提升算法的效率，采用递归算法实现三维几何矩的计算[2]。通过高斯公式，将积分区域从三维退化为二维，再进一步从二维退化为一维。

设Ω为三维几何体所在的区域，S为几何体的边界面。将高斯公式应用于几何矩

2.2.2 几何体三角面片数据的预处理

为了使几何矩的计算更准确，需要对CAD模型的三角面片网格进行简单的预处理：将几何体的三角面片网格平移到以重心为坐标原点的位置，并进行标准化缩放。否则，当几何体离坐标原点距离很远时，在计算几何矩时容易出现2 个大数相减的情况。对三角面片网格进行标准化缩放则是为了避免在计算平移缩放不变矩时出现2 个大数或2 个小数相除的情况。这里，采用的策略是将每个几何体缩放至体积为1。

2.3 几何体相似性判定准则

根据上述判别条件，可确定2 个几何体之间的相似性：

(1) 若同时满足判别条件I～IV，则判定几何体A 和B 完全相似，不存在平移、旋转和缩放变换；

(2) 若仅满足判别条件I 和III，则判定几何体A 和B 相似，仅存在平移变换；

(3) 若仅满足判别条件II 和III，则判定几何体A 和B 相似，仅存在缩放变换；

(4) 若仅满足判别条件I，II 和IV，则判定几何体A 和B 相似，仅存在旋转变换；

(5) 若仅满足判别条件III，则判定几何体A 和B 相似，且存在平移缩放变换；

(6) 若仅满足判别条件I 和IV，则判定几何体A 和B 相似，且存在平移旋转变换；

(7) 若仅满足判别条件II 和IV，则判定几何体A 和B 相似，且存在旋转缩放变换；

(8) 若仅满足判别条件IV，则判定几何体A 和B 相似，且存在旋转平移缩放变换。

2.4 算法效率和误差分析

形状匹配算法主要分为2 步：几何矩不变量的计算和几何体的相似性判定。几何矩不变量的计算采用了一种快速递归算法，将积分区域从三维退化为二维，再进一步从二维退化为一维。每个几何体的形状特征向量由13 个几何矩不变量组成，算法的计算复杂度为O(N)，其中N为几何体表面三角面片的数量。在判断2 个几何体是否相似时，需计算特征向量的相似度，计算复杂度为O(1)。

几何矩是基于几何体的离散表面三角形网格计算得到的，存在一定的离散误差。本文采用的离散参数为：①相邻三角面片法向夹角小于20°；②三角面片到真实曲面的最大距离小于几何体包围盒对角线长度的0.01 倍。为了保证形状匹配的准确性，需要选择一个合适的阈值，使得匹配结果不受离散误差的影响。经过测试，连续几何体与离散几何体之间的相似度通常大于0.95，因此，本文将0.95 作为默认的相似度阈值。

3 测试结果和应用效果

3.1 算法的正确性测试

为了验证算法的有效性，本文对图1 中12 个CAD模型进行了测试。分别对每个原始几何体进行10 次随机的相似变换，从而得到10 个相似几何体。平移量(x,y,z)和缩放量λ为0.0～10.0 之间的随机数，旋转角(θ,φ)为0～π之间的随机数。由此，共得到120 个测试几何体。依次将12 个原始几何体作为目标几何体，与120 个测试几何体进行相似度比对，共计1 440 次。针对不同的相似度阈值，几何体匹配结果见表1。其中，存伪率表示将不相似的2 个几何体判定为匹配的百分比；弃真率表示将相似的2 个几何体判定为不匹配的百分比。从结果可以看出，阈值设置过高会导致弃真率的增加，反之会导致存伪率的增加。当阈值设定为0.9 以上时，能够保证95%以上的成功率，初步验证了算法的有效性。在计算时间方面，测试模型中单个几何体的几何矩特征向量的计算时间不超过0.50 s，每对几何体的相似度匹配判定时间不超过0.01 s。

图1 正确性测试模型 Fig.1 Test models for correctness test

表1 不同阈值下的形状匹配率统计表 (%) Table 1 Statistical table of shape matching rates with different thresholds (%)

3.2 应用效果

基于CAD模型形状匹配算法，本文研发了一个相似几何体识别模块，集成到了自主研发的前处理引擎SuperMesh (http://www.caep-scns.ac.cn/Super Mesh.php)中。同时，在图形用户界面(GUI)上定制了相应的功能，为CAE 软件用户提供特征相似几何体自动化拾取功能。

本文通过几个应用示例来展示该算法的应用效果。图2～4 分别为大坝模型、链条模型和电子学器件3 个测试模型。其中，大坝模型由1 773 个几何体组成，各坝段存在大量相似的几何体。链条模型由200 个几何体组成，包含3 组几何形状相似的几何体，其中第1 组有40 个，第2 组和第3组分别有80 个。电子学器件模型由865 个几何体组成，含有大量相似的柱状结构。经过测试，测试模型中形状相似的几何体都能被准确的筛选出来(图2～4)，验证了方法的有效性。计算时间见表2。模型中所有几何体的形状特征向量均需预先计算，对于数十万规模的三角面片网格，单CPU核的计算时间可控制在30 s 之内。在执行拾取操作时只需完成形状特征向量的相似性比对，对于包含数百上千几何体的CAD模型，仍然能够做到即时响应。

图2 大坝模型的测试结果 Fig.2 Test results of the dam model

表2 时间统计表 Table 2 Statistical table of computing time

在GUI 中执行相似几何体拾取的操作流程分为4 个步骤：①导入CAD模型；②选取一个几何体作为目标几何体；③设置参数；④点击拾取键，执行相似几何体拾取操作。形状匹配算法可自动识别CAD模型中所有与目标几何体相似的几何体，识别出的几何体将在图形交互区高亮显示，并在几何体列表中被选中。图5 为相似几何体识别功能的使用流程示意图。用户可根据需要调整相似度阈值(默认为0.95)、平移、旋转和缩放开关。例如，如果只考虑平移变换，可关闭旋转和缩放开关，由此选出的相似几何体与目标几何体之间只存在平移变换，即只有与目标几何体存在平移变换的几何体能够被筛选出来。

图5 相似几何体识别功能的GUI 操作流程 Fig.5 Procedure of picking similar geometry entities using graphical user interface

4 结 论

本文面向相似变换下的CAD模型匹配问题，基于矩方法及其理论，提出了一个基于几何矩的CAD模型形状匹配算法，用于识别CAD模型中具有相似形状特征的几何体。该算法采用一组几何矩不变量对三维几何体的形状特征进行描述，并根据形状特征向量的相似程度评估几何体之间的相似性，具有易于实现、计算速度快、通用性强等优点。测试结果表明，本算法具有较高的匹配率，并实现了形状特征的快速提取和匹配。最后，本文算法被应用于CAE 软件的相似几何体拾取中，能够通过GUI 交互的方式实时拾取与目标几何体形状特征相似的几何体，取得了良好地应用效果。

在本文算法中，基于几何矩的形状特征向量计算仍有提速的空间，下一步需考虑实现算法的并行化。目前，只考虑了相似变形情况下的形状匹配问题。后续将对算法进行改进，将其推广到非相似变形情况下的形状匹配中。