冷玨琳，張 哲，劉田田，鄭 澎

(1.中物院高性能數(shù)值模擬軟件中心，北京 100088；2.北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京 100088；3.中國工程物理研究院計算機應(yīng)用研究所，四川 綿陽 621900)

形狀匹配是計算機圖形學(xué)、計算機視覺、模式識別等領(lǐng)域的一個重要研究課題。形狀匹配通常分為特征提取和特征匹配2 個步驟[1]，即根據(jù)形狀特征的相似程度來確定形狀的匹配程度?，F(xiàn)有的形狀特征提取方法主要分為4 類[2]：基于幾何結(jié)構(gòu)的形狀特征提取、基于拓?fù)潢P(guān)系的形狀特征提取、基于函數(shù)投影的形狀特征提取和基于統(tǒng)計特征的形狀特征提取。目前還沒有一種通用的方法能夠?qū)Ω黝愇矬w的形狀特征進(jìn)行描述。根據(jù)應(yīng)用場景的不同，選用的方法也不盡相同。例如，基于幾何結(jié)構(gòu)的特征提取方法適用于嚴(yán)格定義的，具有一定剛度或形變較小的模型；基于拓?fù)潢P(guān)系的特征提取方法適用于描述有關(guān)節(jié)或分支的模型，且對噪聲的敏感度高；基于函數(shù)投影的特征提取方法通常要求物體是封閉的；基于統(tǒng)計特性的特征提取方法計算簡單，但對形狀的描述不夠充分，比較適用于粗粒度的匹配。根據(jù)已有的經(jīng)驗總結(jié)[3-4]，對于形狀特征的描述應(yīng)該具備的優(yōu)點包括：信息的描述和分辨能力強、計算速度快、易于存儲和查找、具有幾何變換無關(guān)性、對噪聲具有魯棒性、描述方法獨立于具體的物體對象等。

針對CAD模型的幾何體形狀匹配問題，國內(nèi)外學(xué)者已經(jīng)提出了大量的形狀匹配算法[5-8]以提取和匹配CAD模型的形狀特征，常用的方法包括基于骨架提取和基于表面的匹配算法等?；诠羌芴崛〉钠ヅ渌惴ǖ闹饕枷胧峭ㄟ^提取幾何體的曲線骨架來描述三維形狀特征，且對骨架圖進(jìn)行匹配[5-6]。骨架的提取通常是這類算法的難點，不僅運算量較大，而且對噪聲異常敏感。基于表面的匹配算法采用表面的曲率、法向等特征的分布對三維形狀進(jìn)行描述[7]。例如，將曲面上各采樣點的2 個主曲率的統(tǒng)計直方圖作為三維形狀的特征描述，通過比對曲率分布情況來匹配三維體。這類匹配算法的準(zhǔn)確性需要通過足夠的采樣密度來保證。

對于CAD模型在剛體變換下的形狀匹配問題，最理想的方式是找到獨立于平移、旋轉(zhuǎn)和縮放的形狀特征不變量，然后根據(jù)特征不變量來定義物體形狀的相似度。采用滿足相似變換不變性的幾何矩作為形狀特征的匹配算法[9-14]已經(jīng)在圖像及三維形狀檢索中得到了廣泛地應(yīng)用。本文針對剛體變換下的CAD模型形狀匹配問題，提出一個基于幾何矩的CAD模型形狀匹配算法。

面向CAD模型的形狀匹配算法有著非常廣泛地應(yīng)用。如，用戶在CAE 軟件中對CAD模型設(shè)置材料、計算區(qū)域和邊界條件等屬性時，需要批量選取幾何外形相似的幾何體。采用自動化的幾何體相似性匹配算法進(jìn)行篩選，能夠減輕用戶逐個選取幾何體的工作量，同時可減少手工操作時錯選、漏選情況的發(fā)生。本文提出的CAD模型形狀匹配算法將應(yīng)用于CAE 軟件的相似幾何體拾取，算法的有效性將在實際應(yīng)用中得到驗證。

1 幾何矩方法及相關(guān)理論

1.1 幾何矩

矩在力學(xué)中用于表征物質(zhì)的空間分布。第一篇關(guān)于矩方法的研究發(fā)表于1962 年[9]，其提出了基于二維情況下直角坐標(biāo)系的幾何矩的概念，并在理論上證明了每個物體對象都可由無窮階的幾何矩唯一確定。之后，矩的定義被推廣至三維[10]，并給出了矩不變量的定義，即矩不變量是一種統(tǒng)計特征，滿足平移、旋轉(zhuǎn)和縮放不變性。

矩是形狀密度函數(shù)在核函數(shù)下的積分。通過選取不同的核函數(shù)可定義不同類型的矩，如幾何矩[9-12]、徑向矩[13]、Fourier-Mellin 矩[14]等。

幾何矩的核函數(shù)為直角坐標(biāo)系下的基本集

其中，p1+p2+………+pn=p。幾何矩實質(zhì)上是形狀密度函數(shù)的極數(shù)展開式的系數(shù)。由幾何矩的存在唯一性定理可知[9]，幾何矩序列可唯一確定形狀密度函數(shù)。

1.2 三維幾何體的幾何矩

設(shè)Ω為三維幾何體所在區(qū)域，且?guī)缀误w的密度是均勻分布的(當(dāng)(x,y,z)∈Ω時，f(x,y,z)＝1，(x,y,z)?Ω時，f(x,y,z)＝0)，則該幾何體的幾何矩為

特別地，三維幾何體的體積可由0 階幾何矩m000表示，重心坐標(biāo)可由0 階幾何矩和1階幾何矩決定，即

1.3 三維幾何體的幾何矩不變量

如果一個幾何體與目標(biāo)幾何體之間只存在位置、朝向、大小的區(qū)別，則稱該幾何體與目標(biāo)幾何體之間存在相似變換(或剛體變換)。由幾何矩演變而來的具有平移、旋轉(zhuǎn)、縮放不變性的不變量被統(tǒng)稱為幾何矩不變量。下面給出本文在第2 節(jié)算法中采用的幾何矩不變量的定義[9]。

1.3.1 中心矩

基于重心坐標(biāo)式(3)，可定義三維幾何體的中心矩

中心矩具有平移不變的性質(zhì)。

1.3.2 縮放不變矩

1.3.3 旋轉(zhuǎn)不變矩

基于旋轉(zhuǎn)不變核函數(shù)定義的矩，稱為旋轉(zhuǎn)不變矩。文獻(xiàn)[9-11]中定義了一系列常用的旋轉(zhuǎn)不變矩。本文將采用3 個旋轉(zhuǎn)不變矩，即

其中，D(O,i)，A(O,i,j)均為旋轉(zhuǎn)不變核函數(shù)。D(O,i)為坐標(biāo)點pi=(xi,yi,zi)到坐標(biāo)原點O的距離；A(O,i,j)為坐標(biāo)點pi、坐標(biāo)點pj與坐標(biāo)原點組成的三角形的面積。

類似地，可以定義平移縮放旋轉(zhuǎn)不變矩為

2 基于幾何矩的CAD模型形狀匹配算法

基于幾何矩方法及其相關(guān)理論，本文針對CAD模型的形狀匹配問題，提出一個基于幾何矩的形狀匹配算法，用于識別CAD模型中具有相似形狀特征的幾何體。算法采用幾何矩不變量來描述幾何體的形狀特征，然后根據(jù)幾何矩不變量的相似程度來確定幾何體形狀的相似度。

2.1 算法整體流程

基于幾何矩的CAD模型形狀匹配算法的主要步驟為：

輸入：CAD模型，目標(biāo)幾何體，相似度閾值T。

輸出：CAD模型中所有與目標(biāo)幾何體形狀匹配的幾何體。

步驟1.導(dǎo)入CAD模型。

步驟2.計算CAD模型中每個三維幾何體的一組幾何矩不變量

步驟3.計算CAD模型中各幾何體的形狀特征向量與目標(biāo)幾何體的形狀特征向量的相似度，根據(jù)給定的相似度閾值T和2.3 節(jié)給出的判定準(zhǔn)則，篩選出所有與目標(biāo)幾何體形狀匹配的幾何體，并返回結(jié)果。

在步驟2 中，本文選擇了13 個低階幾何矩不變量組成幾何體的形狀特征向量。對于具有良好光滑性的CAD模型，采用低階幾何矩不變量就能很好地將幾何形狀存在差異的幾何體進(jìn)行區(qū)分。與高階矩相比，低階矩的計算量更小，也更穩(wěn)定，具有計算簡單、易存儲、易查找等優(yōu)點。關(guān)于幾何矩的計算方法將在2.2 節(jié)進(jìn)行介紹。

在步驟3 中，需根據(jù)幾何體形狀特征向量的相似程度判斷幾何體之間是否匹配。其中，相似度閾值T用于界定2 個幾何體是否相似，取值范圍為0～1。給定的閾值T越接近于1，則要求形狀的匹配度越高。相似度判定準(zhǔn)則將在2.3 節(jié)中詳細(xì)介紹。

2.2 三維幾何體的幾何矩計算

CAD模型普遍采用邊界表示法(B-Rep)表示，幾乎所有主流格式的CAD模型都可轉(zhuǎn)化為三角面片網(wǎng)格(STL)表示。因此，本文的算法采用CAD模型的三角面片網(wǎng)格來計算幾何矩，具有較強的通用性。

2.2.1 三維幾何矩的快速計算

形狀匹配算法的運算量主要集中在幾何矩的計算。為提升算法的效率，采用遞歸算法實現(xiàn)三維幾何矩的計算[2]。通過高斯公式，將積分區(qū)域從三維退化為二維，再進(jìn)一步從二維退化為一維。

設(shè)Ω為三維幾何體所在的區(qū)域，S為幾何體的邊界面。將高斯公式應(yīng)用于幾何矩

2.2.2 幾何體三角面片數(shù)據(jù)的預(yù)處理

為了使幾何矩的計算更準(zhǔn)確，需要對CAD模型的三角面片網(wǎng)格進(jìn)行簡單的預(yù)處理：將幾何體的三角面片網(wǎng)格平移到以重心為坐標(biāo)原點的位置，并進(jìn)行標(biāo)準(zhǔn)化縮放。否則，當(dāng)幾何體離坐標(biāo)原點距離很遠(yuǎn)時，在計算幾何矩時容易出現(xiàn)2 個大數(shù)相減的情況。對三角面片網(wǎng)格進(jìn)行標(biāo)準(zhǔn)化縮放則是為了避免在計算平移縮放不變矩時出現(xiàn)2 個大數(shù)或2 個小數(shù)相除的情況。這里，采用的策略是將每個幾何體縮放至體積為1。

2.3 幾何體相似性判定準(zhǔn)則

根據(jù)上述判別條件，可確定2 個幾何體之間的相似性：

(1) 若同時滿足判別條件I～I(xiàn)V，則判定幾何體A 和B 完全相似，不存在平移、旋轉(zhuǎn)和縮放變換；

(2) 若僅滿足判別條件I 和III，則判定幾何體A 和B 相似，僅存在平移變換；

(3) 若僅滿足判別條件II 和III，則判定幾何體A 和B 相似，僅存在縮放變換；

(4) 若僅滿足判別條件I，II 和IV，則判定幾何體A 和B 相似，僅存在旋轉(zhuǎn)變換；

(5) 若僅滿足判別條件III，則判定幾何體A 和B 相似，且存在平移縮放變換；

(6) 若僅滿足判別條件I 和IV，則判定幾何體A 和B 相似，且存在平移旋轉(zhuǎn)變換；

(7) 若僅滿足判別條件II 和IV，則判定幾何體A 和B 相似，且存在旋轉(zhuǎn)縮放變換；

(8) 若僅滿足判別條件IV，則判定幾何體A 和B 相似，且存在旋轉(zhuǎn)平移縮放變換。

2.4 算法效率和誤差分析

形狀匹配算法主要分為2 步：幾何矩不變量的計算和幾何體的相似性判定。幾何矩不變量的計算采用了一種快速遞歸算法，將積分區(qū)域從三維退化為二維，再進(jìn)一步從二維退化為一維。每個幾何體的形狀特征向量由13 個幾何矩不變量組成，算法的計算復(fù)雜度為O(N)，其中N為幾何體表面三角面片的數(shù)量。在判斷2 個幾何體是否相似時，需計算特征向量的相似度，計算復(fù)雜度為O(1)。

幾何矩是基于幾何體的離散表面三角形網(wǎng)格計算得到的，存在一定的離散誤差。本文采用的離散參數(shù)為：①相鄰三角面片法向夾角小于20°；②三角面片到真實曲面的最大距離小于幾何體包圍盒對角線長度的0.01 倍。為了保證形狀匹配的準(zhǔn)確性，需要選擇一個合適的閾值，使得匹配結(jié)果不受離散誤差的影響。經(jīng)過測試，連續(xù)幾何體與離散幾何體之間的相似度通常大于0.95，因此，本文將0.95 作為默認(rèn)的相似度閾值。

3 測試結(jié)果和應(yīng)用效果

3.1 算法的正確性測試

為了驗證算法的有效性，本文對圖1 中12 個CAD模型進(jìn)行了測試。分別對每個原始幾何體進(jìn)行10 次隨機的相似變換，從而得到10 個相似幾何體。平移量(x,y,z)和縮放量λ為0.0～10.0 之間的隨機數(shù)，旋轉(zhuǎn)角(θ,φ)為0～π之間的隨機數(shù)。由此，共得到120 個測試幾何體。依次將12 個原始幾何體作為目標(biāo)幾何體，與120 個測試幾何體進(jìn)行相似度比對，共計1 440 次。針對不同的相似度閾值，幾何體匹配結(jié)果見表1。其中，存?zhèn)温时硎緦⒉幌嗨频? 個幾何體判定為匹配的百分比；棄真率表示將相似的2 個幾何體判定為不匹配的百分比。從結(jié)果可以看出，閾值設(shè)置過高會導(dǎo)致棄真率的增加，反之會導(dǎo)致存?zhèn)温实脑黾?。?dāng)閾值設(shè)定為0.9 以上時，能夠保證95%以上的成功率，初步驗證了算法的有效性。在計算時間方面，測試模型中單個幾何體的幾何矩特征向量的計算時間不超過0.50 s，每對幾何體的相似度匹配判定時間不超過0.01 s。

圖1 正確性測試模型 Fig.1 Test models for correctness test

表1 不同閾值下的形狀匹配率統(tǒng)計表 (%) Table 1 Statistical table of shape matching rates with different thresholds (%)

3.2 應(yīng)用效果

基于CAD模型形狀匹配算法，本文研發(fā)了一個相似幾何體識別模塊，集成到了自主研發(fā)的前處理引擎SuperMesh (http://www.caep-scns.ac.cn/Super Mesh.php)中。同時，在圖形用戶界面(GUI)上定制了相應(yīng)的功能，為CAE 軟件用戶提供特征相似幾何體自動化拾取功能。

本文通過幾個應(yīng)用示例來展示該算法的應(yīng)用效果。圖2～4 分別為大壩模型、鏈條模型和電子學(xué)器件3 個測試模型。其中，大壩模型由1 773 個幾何體組成，各壩段存在大量相似的幾何體。鏈條模型由200 個幾何體組成，包含3 組幾何形狀相似的幾何體，其中第1 組有40 個，第2 組和第3組分別有80 個。電子學(xué)器件模型由865 個幾何體組成，含有大量相似的柱狀結(jié)構(gòu)。經(jīng)過測試，測試模型中形狀相似的幾何體都能被準(zhǔn)確的篩選出來(圖2～4)，驗證了方法的有效性。計算時間見表2。模型中所有幾何體的形狀特征向量均需預(yù)先計算，對于數(shù)十萬規(guī)模的三角面片網(wǎng)格，單CPU核的計算時間可控制在30 s 之內(nèi)。在執(zhí)行拾取操作時只需完成形狀特征向量的相似性比對，對于包含數(shù)百上千幾何體的CAD模型，仍然能夠做到即時響應(yīng)。

圖2 大壩模型的測試結(jié)果 Fig.2 Test results of the dam model

表2 時間統(tǒng)計表 Table 2 Statistical table of computing time

在GUI 中執(zhí)行相似幾何體拾取的操作流程分為4 個步驟：①導(dǎo)入CAD模型；②選取一個幾何體作為目標(biāo)幾何體；③設(shè)置參數(shù)；④點擊拾取鍵，執(zhí)行相似幾何體拾取操作。形狀匹配算法可自動識別CAD模型中所有與目標(biāo)幾何體相似的幾何體，識別出的幾何體將在圖形交互區(qū)高亮顯示，并在幾何體列表中被選中。圖5 為相似幾何體識別功能的使用流程示意圖。用戶可根據(jù)需要調(diào)整相似度閾值(默認(rèn)為0.95)、平移、旋轉(zhuǎn)和縮放開關(guān)。例如，如果只考慮平移變換，可關(guān)閉旋轉(zhuǎn)和縮放開關(guān)，由此選出的相似幾何體與目標(biāo)幾何體之間只存在平移變換，即只有與目標(biāo)幾何體存在平移變換的幾何體能夠被篩選出來。

圖5 相似幾何體識別功能的GUI 操作流程 Fig.5 Procedure of picking similar geometry entities using graphical user interface

4 結(jié) 論

本文面向相似變換下的CAD模型匹配問題，基于矩方法及其理論，提出了一個基于幾何矩的CAD模型形狀匹配算法，用于識別CAD模型中具有相似形狀特征的幾何體。該算法采用一組幾何矩不變量對三維幾何體的形狀特征進(jìn)行描述，并根據(jù)形狀特征向量的相似程度評估幾何體之間的相似性，具有易于實現(xiàn)、計算速度快、通用性強等優(yōu)點。測試結(jié)果表明，本算法具有較高的匹配率，并實現(xiàn)了形狀特征的快速提取和匹配。最后，本文算法被應(yīng)用于CAE 軟件的相似幾何體拾取中，能夠通過GUI 交互的方式實時拾取與目標(biāo)幾何體形狀特征相似的幾何體，取得了良好地應(yīng)用效果。

在本文算法中，基于幾何矩的形狀特征向量計算仍有提速的空間，下一步需考慮實現(xiàn)算法的并行化。目前，只考慮了相似變形情況下的形狀匹配問題。后續(xù)將對算法進(jìn)行改進(jìn)，將其推廣到非相似變形情況下的形狀匹配中。