張 煒， 金 濤

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

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

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

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

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

1 準(zhǔn)備知識

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

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

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

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

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

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

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

1.3 特征點(diǎn)判別指標(biāo)

文獻(xiàn)[12]給出了一個(gè)判斷沿網(wǎng)格頂點(diǎn)某個(gè)方向的光滑程度的指標(biāo)：

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

圖2 頂點(diǎn)沿某個(gè)方向的光滑程度

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

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

2 特征邊識別的新算法

2.1 基本幾何體分析

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

2.2 閾值的確定

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

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

2.3 識別算法描述

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

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

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

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

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

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

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

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

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

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

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

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

3 數(shù)值例子

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

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

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

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

4 總 結(jié)

本文改進(jìn)了判別網(wǎng)格特征點(diǎn)的指標(biāo)，進(jìn)而給出了一種識別三角網(wǎng)格特征邊的新算法。該算法能夠很好地識別出用文獻(xiàn)[11]中的算法不能識別的二面法向夾角比較小的網(wǎng)格邊，從而在時(shí)間復(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é)報(bào)(自然科學(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.