張 煒， 金 濤

(1.浙江機電職業(yè)技術(shù)學(xué)院機械工程學(xué)院，浙江 杭州 310053；2.浙江大學(xué)化工機械研究所，浙江 杭州 310027)

三角網(wǎng)格在求解科學(xué)與工程計算的問題中向來扮演著重要的角色，而近年來，有關(guān)網(wǎng)格特征邊識別的課題正在成為一個研究熱點。這是因為，隨著三維掃描儀的成本下降和質(zhì)量提升，三角網(wǎng)格的生成變得輕車熟路，點云或網(wǎng)格曲面在計算機輔助設(shè)計(computer aided design，CAD)與計算機圖形(computer graphics，CG)中使用更加普遍，作用日益凸顯，從而大大促進了網(wǎng)格重建[1]、網(wǎng)格編輯[2]、網(wǎng)格參數(shù)化[3]、網(wǎng)格光順[4]等一系列以網(wǎng)格計算和網(wǎng)格變換為主題的研究。而上述研究都是以網(wǎng)格特征邊識別為先決條件的。另一方面，借助掃描儀與CAD軟件(如AutoCAD、Pro/E，、Unigraphics、Solidworks等)，人們能夠快速獲得逼近三維幾何物體表面的三角網(wǎng)格，這些三角網(wǎng)格的數(shù)據(jù)主要保存在一種業(yè)已成為行業(yè)標準的STL文件當中；而基于STL文件的計算機輔助制造(computer aided manufacturing，CAM)，必須以三角網(wǎng)格的特征邊識別為首要前提。綜上可知，給出一種魯棒的網(wǎng)格特征邊識別算法，對CAD/CAM/CG而言，實為當務(wù)之急。

關(guān)于網(wǎng)格特征邊已經(jīng)有大量識別算法。文獻[5]首先通過判斷曲率的局部極值來識別網(wǎng)格的特征點，然后根據(jù)主曲率分析方法(principal curvature analysis，PCA)用3次B-樣條把這些特征點連接成特征線。文獻[6]提出了利用網(wǎng)格邊的二面法向夾角(2個網(wǎng)格三角形平面的法向夾角)和邊角(兩條相鄰的邊的外轉(zhuǎn)角)之間關(guān)系的一種算法來識別特征邊。此后文獻[7]又改進了這種算法，并進一步提出了一種識別C2不連續(xù)網(wǎng)格邊的算法。文獻[8]以張量投票(tensor voting)為工具對網(wǎng)格頂點進行分類，然后用區(qū)域擴張合并的遞歸算法對網(wǎng)格進行分塊。假設(shè)其中n是頂點p的1-ring的三角形ii的法向，分別是矩陣A的特征值和相應(yīng)的特征向量。所謂張量投票是指用分別逼近曲面法向(e1)和主曲率(e2,e3)。但是文獻[9]指出，這種逼近方式有以下缺點：①法向?qū)?quán)重十分敏感而且符號未定，可能是內(nèi)法向也可能是外法向；②當頂點所在的邊的二面法向夾角大于90°時，用逼近法向的效果比用為優(yōu)；③當頂點所在的邊的二面法向夾角等于90°時，三個特征向量都不能很好地逼近法向。文獻[10]通過網(wǎng)格重建得到一個各向同性的對特征敏感的度量，然后利用局部鄰域的積分不變量來識別多尺度網(wǎng)格上的特征。以上算法在精度上都不盡如人意。文獻[11]通過統(tǒng)計網(wǎng)格邊的二面法向夾角信息，計算出二面法向夾角的平均值和標準差來確定特征邊的閾值，并通過遞歸的算法來識別特征邊。但該文有一個明顯的缺點：不能識別二面法向夾角比較小的特征邊，如圖1所示。最近，文獻[12]利用網(wǎng)格頂點處的法向和曲率半徑等離散微分幾何工具給出了一個判別特征點的指標，并給出了識別網(wǎng)格特征點的魯棒算法。但該算法存在以下問題：①文中給出的指標十分不合理，影響了識別的準確度；②文中只是給出了特征點的識別算法，并沒有給出特征邊的識別算法，難以在CAM中加以使用。

綜上所述，迄今為止，各種網(wǎng)格特征邊識別算法都存在著明顯缺陷，嚴重地阻礙了以網(wǎng)格計算和網(wǎng)格變換為主題的研究，也影響了基于STL文件的CAM的正常開展。因此，本文深入剖析了文獻[11]和文獻[12]的理論機理中的軟肋，首先改造文獻[12]所提出的識別特征點的指標，進而利用改造所得的新指標來改進文獻[11]中的算法，得到了一種識別網(wǎng)格特征邊的新的魯棒的算法，并予編程實現(xiàn)。在此基礎(chǔ)上，本文把新算法與文獻[11]中的算法作了深入的數(shù)值實例比較，并給出了圖形顯示與繪制。這些例子表明，文獻[11]的算法易導(dǎo)致失效，而本文的算法可挽救這種失效。盡管文獻[11]是2004年發(fā)表的會議論文，但是本文與之比較還是有意義的，原因如下：①近年來，在網(wǎng)格特征邊識別并沒有取得突破性的發(fā)展，沒有一個魯棒的算法可以在機械制造領(lǐng)域得到廣泛的應(yīng)用；②文獻[11]在特征邊識別領(lǐng)域有著重要的意義，近年來許多相關(guān)的文獻都與之做出比較，如文獻[6,12]。

圖1 用文獻[11]中的算法所得到的網(wǎng)格特征邊

1 準備知識

1.1 網(wǎng)格頂點的法向

由于網(wǎng)格曲率的估計依賴于網(wǎng)格頂點法向，這里先給出網(wǎng)格法向的計算方法。必須指出，網(wǎng)格頂點的法向是不能精確計算的，已經(jīng)有很多關(guān)于網(wǎng)格頂點的法向估計的文章。其中最常用的兩種方法是加權(quán)平均和張量投票。正如文獻[9]所指出的，張量投票方法存在著一些缺點，所以本文采用加權(quán)平均來估計網(wǎng)格頂點的法向。假設(shè)頂點p的1-ring三角形的法向為vi,i=1,...,m，m是頂點p的1-ring的三角形數(shù)量，則頂點p處的法向可計算如下：這里wi是權(quán)重，本文取wi為三角形的面積。

1.2 網(wǎng)格頂點的曲率

由于網(wǎng)格頂點的曲率也是不能精確計算的，已經(jīng)有很多關(guān)于網(wǎng)格頂點的曲率估計的文章。 常用的網(wǎng)格頂點曲率估計的方法有離散算法[13]、曲面局部擬合法[14-15]、曲率張量投票法[8]等。文獻[15]指出，用離散算法估計網(wǎng)格頂點的曲率在某些情況下是不準確的，而曲率張量投票的逼近算法也存在著很多缺點[9]，因此本文采用曲面局部擬合的方法來估計網(wǎng)格頂點的曲率。文獻[14]采用二次曲面去擬合頂點與其1-ring的頂點，這種擬合方法是一定插值該頂點的。文獻[15]也采用二次Bézier曲面去擬合頂點與其1-ring的頂點，不過這種擬合方法是不一定插值該頂點的。盡管文獻[14]中的曲率估計方法不是最準確的，但是用它來估計網(wǎng)格的局部微分性質(zhì)已經(jīng)足夠了。為了論文的完整性，下面簡單介紹一下這種方法。

在網(wǎng)格頂點p處建立一個局部坐標系，其中坐標系的z軸為該點處的外法向。采用二次曲面s(u,v)=(u,v,au2+buv+cv2)來擬合該頂點及其1-ring的頂點，這里u,v都是局部坐標系下的自變量。假設(shè)該頂點的1-ring頂點在局部坐標系下的坐標值為(ui,vi,hi),i=1,...,m，可以得到下面的線性方程組：

用最小二乘法求解式(2)可得二次曲面的系數(shù)a,b,c。由此可得該頂點處的主曲率和相應(yīng)的主曲率方向如下：

假設(shè)向量m在主曲率方向m1,m2張成的平面上，且m與m1的夾角為θ，則曲面s(u,v)沿m方向的法曲率為：

1.3 特征點判別指標

文獻[12]給出了一個判斷沿網(wǎng)格頂點某個方向的光滑程度的指標：

這里iT是頂點p鄰域的三角形，bi是該三角形的重心，vi是該三角形的法向，ki是在向量pbi和法向n張成的平面上頂點p處的法曲率(見圖2)。由于文獻[12]中并沒證明式(3)中的變量滿足：

圖2 頂點沿某個方向的光滑程度

而且在實驗中發(fā)現(xiàn)式(4)在某些情況下是不成立的，這將干擾識別的準確性。在網(wǎng)格大小比較均勻處，式(4)是成立的，因為從圖2中可以看出，若在該處擬合出來的曲面是比較接近三角形，則在bi處的曲率半徑會小于，即式(4)成立。但是，在網(wǎng)格大小不均勻處，擬合的二次曲面曲率比較大，則在bi處的曲率半徑會大于。進一步地，可以發(fā)現(xiàn)當這種情況出現(xiàn)時，點p處是不光滑的。因此，本文把式(3)改成如下：

若SHI(p)大于某個給定的閾值α，則認為頂點p為特征點，否則認為頂點p為光滑點。

2 特征邊識別的新算法

2.1 基本幾何體分析

我們對Unigraphics導(dǎo)出STL文件的基本幾何體(見圖3)進行分析，可得出如下結(jié)果：誤差取默認的輸出誤差，圓柱體的上下底面的邊界為特征邊，側(cè)面的邊(非底面的邊)的最大二面法向夾角為9°；圓錐體的底面邊界為特征邊，側(cè)面的邊(非底面的邊)的最大二面法向夾角為8°；球體本身沒有特征邊，球表面網(wǎng)格邊的最大二面法向夾角為4.5°；圓環(huán)體本身沒有特征邊，圓環(huán)表面網(wǎng)格邊的最大二面法向夾角為6.6°。這里網(wǎng)格邊的二面法向夾角定義為以該邊為公共邊的兩個三角形的外法向夾角。

2.2 閾值的確定

本文的算法需要用到2個閾值α1,α2，其中α1是網(wǎng)格邊的二面法向夾角的閾值，而α2是網(wǎng)格特征點的閾值。根據(jù)2.1節(jié)的分析，可將把二面法向夾角大于20°的網(wǎng)格邊定義為特征邊，而把二面法向夾角介于[α1,20°]之間的網(wǎng)格邊成為次特征邊。若式(6)中定義的指標SHI(p)大于α2，則認為頂點p為特征點。這里的閾值α1,α2都可以通過用戶交互確定。另外，若α1取得比較大，則可能會漏掉一些特征邊，若α1取得比較小，則可能會識別出一個錯誤的特征邊。α1的取值還跟模型的網(wǎng)格大小質(zhì)量等有關(guān)，所以無法確定出一種最優(yōu)通用的閾值方案。一種常用的方案是取所有二面法向夾角的均值作為閾值的初始值。對于閾值α2, 可以參考文獻[12]給出的一種通過求解某個最優(yōu)化問題的方法來確定該閾值。

圖3 基本幾何體的特征邊分析

2.3 識別算法描述

文獻[12]中的識別特征點指標比文獻[11]中的根據(jù)最小主曲率的主曲率方向來識別特征點的準確度要高，故用改造文獻[12]的識別特征點的指標來改進文獻[11]的算法，能提高識別的準確度。本節(jié)給出網(wǎng)格特征邊識別的新算法描述。假設(shè)輸入的網(wǎng)格模型是M=(V,E,F)，這里V,E,F分別是模型的頂點、邊、三角形的數(shù)組。

網(wǎng)格特征邊識別算法：

輸入.三角網(wǎng)格模型M=(V,E,F)，以及用戶的輸入?yún)?shù)μ,λ。

輸出.三角網(wǎng)格模型的特征邊數(shù)組L。

步驟1.初始化一個臨時數(shù)組tempL，計算網(wǎng)格模型的相關(guān)數(shù)據(jù)，如網(wǎng)格邊的二面法向夾角， 網(wǎng)格特征點的判斷指標SHI(p)。

步驟2.計算網(wǎng)格邊的二面法向夾角的平均值θavr和標準差θsd，計算閾值α1=θavr+μθsd，α2=λ。

步驟3.對網(wǎng)格模型作一次平面檢測。

步驟4.對每個網(wǎng)格頂點p進行檢測，判斷該頂點是否為特征點。若指標SHI(p)大于α2，則認為頂點p為特征點，否則認為它為光滑點。

步驟5.對每條網(wǎng)格邊進行檢測。 若該邊的二面法向夾角大于20°，則把它放入特征邊數(shù)組L；若該邊的二面法向夾角屬于[α1,20°]，則把它放入臨時數(shù)組tempL。

步驟6.對臨時數(shù)組tempL中的邊進行過濾。若邊的兩個端點中不都是特征點，則把該邊從tempL刪除。

步驟7.對臨時數(shù)組tempL中的邊按二面法向夾角從大到小的順序進行排序。這里采用二分法排序，時間復(fù)雜度為O(nlogn)。

步驟8.對臨時數(shù)組tempL中的邊進行特征邊判別。假設(shè)p1,p2為該邊的兩個端點，若與p1,p2相連的邊中已經(jīng)有兩條或兩條以上的特征邊，則認為該邊不是特征邊，否則認為該邊是特征邊并將之放入特征邊數(shù)組L。

3 數(shù)值例子

本文實現(xiàn)了文獻[11]中的算法，并將之與本文的新算法作比較。實現(xiàn)新算法的PC配置如下：Intel Core2 CPU T7250，RAM 2 GB。本文算法是用C++和OpenGL實現(xiàn)的。圖4給出的是一個經(jīng)過圓角處理的長方體，從圖1可以看出，文獻[11]的算法不能得到二面法向夾角很小(黑色圈出的部分)的那些特征邊，而圖4中用新算法卻得到了所有的特征邊。圖5和圖6分別給出了一個葉輪模具和一個零件模具的網(wǎng)格模型，圖5(b)和圖6(b)中用黑色圓圈圈出的部分皆為沒有識別到的地方，由之可以看出，文獻[11]的算法易導(dǎo)致失效，而本文的算法挽救了這種失效。

圖4 本文算法所得到的網(wǎng)格特征邊

圖5 對葉輪的網(wǎng)格特征邊作識別對比 (最右邊的一列圖是右2列圖中黑色矩形內(nèi)部分的放大圖)

圖6 對零件的網(wǎng)格特征邊作識別對比

4 總 結(jié)

本文改進了判別網(wǎng)格特征點的指標，進而給出了一種識別三角網(wǎng)格特征邊的新算法。該算法能夠很好地識別出用文獻[11]中的算法不能識別的二面法向夾角比較小的網(wǎng)格邊，從而在時間復(fù)雜度一樣的基礎(chǔ)上大大提高CAM中模具加工的效能。本文以后的工作是將本文的算法推廣到C2不連續(xù)網(wǎng)格邊的識別。

[1]Volpin O, Sheffer A, Bercovier M, Joskowicz L.Mesh simplification with smooth surface reconstruction [J].Computer Aided Design, 1998, 30(11): 875-882.

[2]Zhang Shaoting, Huang Junzhou, Metaxas D N.Robust mesh editing using Laplacian coordinates [J].Graphics Models, 2011, 73(1): 10-19.

[3]Morigi S.Feature-sensitive parameterization of polygonal meshes [J].Applied Mathematics and Computation, 2009, 215(4): 1561-1572.

[4]Li Zhong, Ma Lizhuang, Jin Xiaogang, Zheng Zuoyong.A new feature-preserving mesh-smoothing algorithm [J].The Visual Computer, 2009, 25(2): 139-148.

[5]寧上官, 劉 斌.三角網(wǎng)格模型特征線提取方法[J].華僑大學(xué)學(xué)報(自然科學(xué)版), 2010, 31(5): 487-490.

[6]Jiao Xiangming, Heath M T.Feature detection for surface meshes [C]//Proceedings of 8thInternational Conference on Numerical Grid Generation in Computational Field Simulations.Honolulu, 2002:705-714.

[7]Jiao Xiangming, Bayyana N R.Identification of C1and C2discontinuities for surface meshes in CAD [J].Computer Aided Design, 2008, 40(2): 160-175.

[8]Kim H S, Choi H K, Lee K H.Feature detection of triangular meshes based on tensor voting theory [J].Computer Aided Design, 2009, 41(1): 47-58.

[9]Jiao Xiaoming, Alexander P J.Parallel feature-preserving mesh smoothing [C]//Proceedings of International Conference on Computational Science and Its Applications.Singapore, 2005: 1180-1189.

[10]Lai Yukun, Zhou Qianyi, Hu Shimin, Wallner J,Pottmann H.Robust feature classification and editing [J].IEEE Transaction on Visualization and Computer Graphics, 2007, 13(1): 34-45.

[11]Baker T J.Identification and preservation of surface features [C]//Proceedings of 13thInternational Meshing Roundtable.Williamsburg, 2004: 299-310.

[12]Di Angelo L, Di Stefano P.C1continuities detection in triangular meshes [J].Computer Aided Design, 2010,42(9): 828-839.

[13]Myers M, Desbrun M, Schroder P, Barr A H.Discrete differential geometry operators for triangulated 2-manifolds [M].Berlin: Springer-Verlag Press, 2002:34-57.

[14]Haman B.Curvature approximation for triangulated surfaces [M].Berlin: Springer-Verlag Press, 1993:139-153.

[15]Razdan A, Baea M.Curvature estimation scheme for triangle meshes using biquadratic Bézier patches [J].Computer Aided Design, 2005, 37(14): 1481-1491.