賀 強，張樹生，白曉亮

(西北工業(yè)大學 現(xiàn)代設計與集成制造教育部重點實驗室，710072 西安，hqcq@mail.nwpu.edu.cn)

曲面特征恢復的三角網(wǎng)格模型孔洞修補算法

賀 強，張樹生，白曉亮

(西北工業(yè)大學 現(xiàn)代設計與集成制造教育部重點實驗室，710072 西安，hqcq@mail.nwpu.edu.cn)

為了恢復三角網(wǎng)格模型中的孔洞處的真實形狀，提出一種曲面特征恢復的孔洞修補算法.首先對模型中檢測出的孔洞進行三角化并細分，完成孔洞的粗修復.然后利用高斯球確定孔洞的鄰域及其曲面類型，對二次曲面類型的孔洞鄰域進行非線性最小二乘擬合以獲得曲面代數(shù)方程，對自由曲面類型的孔洞鄰域則進行B樣條曲面擬合.最后根據(jù)孔洞鄰域的曲面方程獲得孔洞內(nèi)新增頂點在曲面上的精確位置，完成孔洞修復.實驗結(jié)果表明，該孔洞修補算法不僅能完成孔洞區(qū)域的三角形填充，還能準確恢復出孔洞區(qū)域的曲面特征.

高斯球;孔洞鄰域;曲面特征;孔洞修補

由于物體表面的反射屬性、結(jié)構(gòu)閉塞等因素的影響，測量獲得的三角網(wǎng)格模型中不可避免地存在孔洞.孔洞不僅影響模型的外觀，而且不利于后續(xù)的幾何處理.為了獲得完整的網(wǎng)格模型，孔洞修補尤為關鍵.經(jīng)過大量研究，出現(xiàn)了基于網(wǎng)格和基于體數(shù)據(jù)2類具有代表性的孔洞修補算法.基于網(wǎng)格的修補算法［1－8］首先尋找網(wǎng)格模型中的孔洞，然后對孔洞進行三角化并細分，最后對新增的網(wǎng)格頂點的位置進行光順調(diào)整.基于網(wǎng)格的修補僅僅作用在孔洞及其鄰域，可以保持遠離孔洞區(qū)域的網(wǎng)格結(jié)構(gòu)不變，修補后孔洞區(qū)域能與周圍區(qū)域光滑融合，但這類算法對輸入的模型有一定的質(zhì)量要求且沒有或較少考慮孔洞處的曲面信息，常常丟失孔洞處的曲面特征.基于體數(shù)據(jù)的網(wǎng)格修補算法［9－12］首先將網(wǎng)格模型轉(zhuǎn)換為一種中間的體數(shù)據(jù)來表示，然后在體數(shù)據(jù)上進行修補操作.這類算法在處理三角面片自相交和重合的情況時比基于網(wǎng)格的算法更具優(yōu)勢，并能保證修復后輸出的網(wǎng)格是封閉的，但輸出的網(wǎng)格完全改變了原有模型中的連接關系，從而導致模型特征的丟失并會產(chǎn)生大量的、狹長的三角片.

上述算法較少或基本沒有考慮孔洞處的曲面特征，因而在對含有大量二次曲面的工業(yè)產(chǎn)品外形上的孔洞進行修補時，雖然可以得到完整且光順的模型，但同時會丟失孔洞處的曲面特征.為了恢復孔洞處的真實形狀，本文提出一種曲面的特征恢復的孔洞修補算法.首先對模型中檢測出的孔洞進行粗修復，然后利用高斯球確定孔洞鄰域及其曲面類型，最后根據(jù)孔洞鄰域的曲面信息對孔洞內(nèi)新增頂點的空間位置進行調(diào)整，獲得精確的修復結(jié)果.

1 孔洞的粗修復

對三角網(wǎng)格模型中檢測出的孔洞，直接進行三角化處理，并根據(jù)孔洞周圍網(wǎng)格的平均邊長對孔洞區(qū)域內(nèi)新增的三角形進行細分，這個過程稱為孔洞的粗修復.

1.1 孔洞的識別

三角形網(wǎng)格模型是由一系列頂點和三角形構(gòu)成.2個相鄰三角形通常有1條公共邊，稱這樣的邊為網(wǎng)格的內(nèi)部邊，而網(wǎng)格模型中的邊界邊和孔洞邊通常只屬于1個三角形.根據(jù)這一性質(zhì)，對封閉的網(wǎng)格模型，如果存在某條邊僅屬于1個三角形，則此邊就是構(gòu)成孔洞多邊形的1條邊.由這樣的邊首尾相連就形成了1個孔洞.對于非封閉模型，提取出的孔洞中還包括了模型的邊界輪廓.一般而言，孔洞的邊的數(shù)量都小于邊界輪廓的邊的數(shù)量.根據(jù)這一經(jīng)驗，對檢測出的所有孔洞按照邊數(shù)量排序，不修補邊數(shù)量最大的孔洞.對于特殊情況，則需要用戶交互選擇模型輪廓，并將其從孔洞集合中剔除.

1.2 孔洞的三角化和細分

對孔洞多邊形使用Delaunay方法直接三角化，完成孔洞區(qū)域的拓撲填充.對面積較大的孔洞，獲得的修復網(wǎng)格與原始孔洞周圍的網(wǎng)格的采樣密度差異較大，這不僅影響了視覺效果而且破壞了孔洞區(qū)域及其鄰域的拓撲一致性.為了獲得與孔洞鄰域網(wǎng)格密度相近的網(wǎng)格分布，本文對新增的三角形的邊進行細分.細分的原理是首先計算孔洞周圍(選擇孔洞周圍網(wǎng)格6～10層)三角形邊長的平均值，然后檢測孔洞內(nèi)新增三角形的三邊長，如果某邊邊長超過平均值，則將該邊從中點處分裂，邊分裂過程如圖1所示.

圖1 邊分裂

圖2表示了孔洞的粗修復過程.其中圖2(a)是含有孔洞的三角網(wǎng)格模型，圖2(b)是采用Delaunay三角化的結(jié)果，圖2(c)是細分光順的結(jié)果.

圖2 孔洞的粗修復

2 孔洞鄰域及其曲面特征

為了使孔洞修補的結(jié)果能較好逼近孔洞處的真實曲面，利用孔洞周圍的網(wǎng)格(即確定孔洞的鄰域)尤為關鍵.物體的外形通常都由大量的二次曲面和少量的自由曲面構(gòu)成.在工業(yè)領域，這種情況尤為突出.鑒于高斯球［13－14］可以很好地識別二次曲面，本文利用三角網(wǎng)格模型良好的拓撲連接關系，向孔洞的周圍逐層擴展，同時利用高斯球識別出擴展部分的曲面類型，并將曲面類型不一致的頂點剔除.孔洞鄰域就是剩余的與孔洞處曲面類型一致的頂點的集合.

2.1 網(wǎng)格頂點的高斯映射

在完成孔洞的粗修復后，向孔洞周圍擴展6～10層三角形，獲得初步的孔洞鄰域點集.高斯映射是將曲面上任意一點的單位法矢的起點移動到坐標原點的過程.曲面上的點進行高斯映射后所得的結(jié)果是單位球內(nèi)的一個點集.這個單位球就是高斯球.由高斯映射的定義可知，網(wǎng)格頂點的法矢量即為其在高斯球上的坐標.

2.1.1 平面的高斯球

平面點的法矢量均相等，因此平面的高斯球應該是一個點，但實際應用中，由于法矢估算不準確或噪聲等原因，平面網(wǎng)格在高斯映射后的點會聚集在一個點區(qū)域中，如圖3(a)所示.

2.1.2 圓柱面的高斯球

圓柱面上任意點的法向量都應該與其軸線的單位矢量n垂直，因此圓柱面上任意一點的高斯球坐標都在過原點且法向為n的平面上，該平面可以表示為(n，0).實際應用中，圓柱面類型的網(wǎng)格頂點在高斯映射后的點集近似地聚集在一個環(huán)狀區(qū)域而不是嚴格位于同一平面，如圖3(b)所示.該區(qū)域可以被最小二乘法擬合成一張平面(n*，d)，d近似為0.

2.1.3 圓錐面的高斯球

圓錐上任意一點的法矢與軸線方向n的夾角與圓錐角的1/2 α互余.因此圓錐上任意點的高斯球坐標構(gòu)成平面(n，sin α).由于網(wǎng)格噪聲等原因，實際的圓錐網(wǎng)格高斯映射后的點集近似地聚集在一個環(huán)狀區(qū)域而不是嚴格位于同一平面，如圖3(c)所示.該區(qū)域可以被擬合成平面(n*，d).

2.1.4 球面的高斯球

球面上任意點的法向各不相同，因此球面上的點的高斯球坐標各不相同，如圖3(d)所示.球面上各個點及其法矢量構(gòu)成的直線之間的交點理論上應該都在球心處，而對實際應用中的球面網(wǎng)格則是在包含球心的一個點區(qū)域中，與平面的高斯映射類似.本文根據(jù)這個性質(zhì)來識別球面.

圖3 二次曲面的高斯球

2.2 孔洞鄰域的確定

對初步的孔洞鄰域點集進行高斯映射后，利用高斯球就可以對其進行精煉，進而得到孔洞鄰域的曲面類型和鄰域點集，孔洞鄰域及其曲面類型的確定過程為:

HN為初步的孔洞鄰域點集，其高斯球坐標形成的點集為C，建立兩者之間一一對應的索引關系.首先利用通用點聚類法判斷HN是否為平面，該過程為:計算C的中心點Pc，然后計算C中每一個點與Pc的距離，如果該距離小于給定的閾值，則保留該點，否則拋棄該點.完成上述操作后，減小閾值，遞歸調(diào)用上述過程，直至C中的所有點到Pc的距離都小于給定的全局閾值或迭代次數(shù)達到設定的上限.此時若C中點的數(shù)量與HN中點的數(shù)量之比大于設定的有效點集率，則HN是1個平面區(qū)域.

如果平面判定失敗，則依次進行圓柱面、圓錐面的判定.此時需要判斷C是否處于一塊有效的平面區(qū)域，利用通用平面聚類法完成平面區(qū)域的判定.該過程為:對C中的點行最小二乘平面擬合，得到一個中心平面Pp.計算C中的點到平面Pp的距離，如果該距離小于給定的閾值，則保留該點，反之刪除該點.完成上述操作后，減小閾值，遞歸調(diào)用上述過程，直至C中的所有點到Pp的距離都小于給定的全局閾值或迭代次數(shù)達到設定的上限.根據(jù)C中剩下的點，利用建立的索引關系，就可以得到HN中與孔洞區(qū)域曲面類型一致的頂點集合，即孔洞的鄰域.

以上曲面判定都失敗后，還需對HN進行球面的判定.利用HN中的頂點和C中的法矢量確定的直線的交點來判定球面.若是球面類型的網(wǎng)格，這些直線會相交于球心附近，形成一個點區(qū)域，對該點區(qū)域進行通用點集聚類即可確定孔洞的鄰域.

若以上曲面類型的判定都失敗，本文將這樣的孔洞鄰域歸為自由曲面類型.孔洞鄰域則為擴展部分的所有頂點.

2.3 孔洞鄰域的曲面擬合

對平面孔洞，無需擬合平面(因為三角化后細分新增的頂點還是在平面上).對球面、圓柱面、圓錐面孔洞，則根據(jù)孔洞的鄰域，使用非線性最小二乘擬合獲得孔洞鄰域的代數(shù)方程.該方程作為孔洞內(nèi)新增頂點位置調(diào)整的依據(jù);若孔洞鄰域是自由曲面類型，本文采用文獻［15］中的雙三次B樣條曲面擬合算法來確定孔洞鄰域的潛在曲面.

球面的代數(shù)方為

式中:r為半徑;(x0，y0，z0)為球心，是待求的參數(shù).

圓柱面的代數(shù)方程為

式中:(x0，y0，z0)為待求的圓柱軸線上任意一點;(nx，ny，nz)為待求的圓柱的軸線方向;r為待求的截面圓的半徑.

圓錐面的代數(shù)方程為

式中:α為待求的圓錐角的1/2;(x0，y0，z0)為待求的圓錐的頂點;(nx，ny，nz)為待求軸線方向.

本文采用 Levenberg-Marquardt［16］算法求解曲面方程中的待求參數(shù).設:F(x)為二次曲面的代數(shù)方程;x={x1，x2，…，xn}為待求參數(shù)組成的向量.對曲面方程中待求參數(shù)進行泰勒展開得

設xk為第k次迭代獲得的待求參數(shù)向量.該求解過程為

式中:δ為x中參數(shù)的偏差;J為m×n階雅可比矩陣;m為孔洞鄰域點的數(shù)量;μ＞0為用戶設定的參數(shù)，初始設為 10－6.

上述過程首先需要設定待求參數(shù)即x0的初始估計值.對球面擬合，球心的初值設為孔洞鄰域內(nèi)任意2點和其法矢量構(gòu)成的直線的交點，半徑初值為其中1點到該交點的距離.對圓柱面擬合，利用任意2點及其法矢量，根據(jù)法矢量與軸線垂直的性質(zhì)就可估算軸線，其中一點到軸線的距離即為半徑.對圓錐面擬合，選擇不在同一平面的3點及其法矢量，利用法矢量與軸線的夾角為圓錐角的1/2這個性質(zhì)就可計算半徑、軸線方向、圓錐角、圓錐頂點.

3 孔洞內(nèi)新增頂點的位置調(diào)整

由孔洞粗修復產(chǎn)生的新增頂點，如果孔洞鄰域曲面類型是平面，則無需進一步處理，修補完成.若是球面、圓柱面和錐面中的任何一種，則需要利用過程2中獲得的代數(shù)方程計算新增頂點在曲面上的真實空間位置.對新增頂點pi=(xi，yi，zi)，其法矢量n=(nx，ny，nz)，得到過該點且方向為法矢量的直線方程為

式(6)與孔洞鄰域曲面代數(shù)方程聯(lián)立求解就得到孔洞內(nèi)的新增點的空間坐標.若孔洞鄰域是自由曲面類型，本文采用文獻［15］的網(wǎng)格參數(shù)化算法對粗修復的網(wǎng)格進行參數(shù)化，獲得孔洞處填充的頂點的參數(shù)值，再將參數(shù)值帶入由孔洞鄰域擬合得到的B樣條曲面方程中就獲得了孔洞處新增頂點的空間位置.

4 實驗結(jié)果

在VS2005的環(huán)境下實現(xiàn)了算法，在主頻1.8 G，內(nèi)存512 M的PC機上運行程序，且采用VTK(Visualization Toolkit一個面向?qū)ο蟮目梢暬悗?顯示孔洞模型及其修補結(jié)果.為了更好地說明孔洞修補過程，本文以圖4的球面孔洞為例詳述其修補過程.圖4(a)是原始的球面孔洞模型，其真實球心坐標為(0，0，0)，半徑為 8 mm.圖4(b)是孔洞粗修復的結(jié)果即完成了孔洞多邊形的三角化和細分.對孔洞鄰域進行擬合得到的球心為(－0.000 422，0.000 116，－0.000 208)，半徑為7.999 997 mm，與真實情況的差異非常小，因此利用擬合的結(jié)果能精確恢復孔洞處的曲面特征，如圖4(c)所示.圖5是圓錐面孔洞的修復，圖6是自由曲面孔洞的修復，圖7是多種曲面類型孔洞的修復.從多種曲面類型孔洞的修補實驗可知，除了平面類型的孔洞，孔洞粗修復后產(chǎn)生的新頂點都偏離了其在曲面上的真實位置，而采用本文算法處理后，能精確地恢復孔洞處的曲面特征.

圖4 球面孔洞及其修復

圖5 圓錐面孔洞及其修復

圖6 自由曲面孔洞及其修復

圖7 機械零件模型中的孔洞及其修復

5 結(jié)論

1)通過提取孔洞多邊形，孔洞三角化和細分完成了孔洞區(qū)域的拓撲填充.

2)利用高斯映射確定了孔洞鄰域及其曲面類型，并利用非線性最小二乘擬合得到相關的曲面參數(shù).

3)利用孔洞鄰域的曲面信息，完成了孔洞區(qū)域曲面特征恢復的修補，實驗結(jié)果證明了算法的有效性.

［1］LIEPA P.Filling holes in meshes［C］//Proceedings of Eurographics/ACM SIGGRAPH symposium on Geometry processing. Switzerland:Eurographics Association，2003:200－205.

［2］CHUI C，LAI M J.Filling polygonal holes using C1 cubic triangular spline patches［J］.Computer Aided Geometry Design，2000，17(4):297 －307.

［3］BRANCH J，PRIETO F，BOULANGER P.Automatic Hole-Filling of triangular mesh using local radial basis function［C］//Proceedings of the third International Symposium on 3D Data Processing，Visualization and Transmission.Washington，NC:IEEE，2006:727 －734.

［4］TEKUMALLA L S，COHEN E.A Hole-Filling Algorithm for Triangular Meshes［M］.USA:Technical Report UUCS-04-019，School of Computing University of Utah，2004.

［5］PERNOT J P，MORARU G，VERON P.Filling holes in meshes using a mechanical model to simulate the curvature variation minimization［J］.Computer ＆ Graphics，2006，30(6):892－902.

［6］JUN Y.A piecewise hole filling algorithm in reverse engineering［J］.Computer-Aided Design，2005，37(2):263－270.

［7］張麗艷，周儒榮，周來水.三角網(wǎng)格模型孔洞修補算法研究［J］.應用科學學報，2002，20(3):221－224.

［8］張潔，岳瑋寧，王楠，等.三角網(wǎng)格模型的各向異性孔洞修補算法［J］.計算機輔助設計與圖形學學報，2007，19(7):893 －897.

［9］DAVIS J，MARSCHNER S R，GARR M，et al.Filling holes in complex surfaces using volumetric diffusion［C］//Proceedings of First International Symposium on 3D Data Processing，Visualization and Transmission.Washington，NC:IEEE，2002:428－438.

［10］NOORUDDIN F S，TURK G.Simplification and repair of polygonal models using volumetric techniques［J］.IEEE Transactions on Visualization and Computer Graphics，2003，9(2):191 －205.

［11］JU T.Robust repair of polygonal models［J］.ACM Transaction on Graphics，2004，23(3):888 －895.

［12］BISCHOFF S，PAVIC D，KOBBELT L.Automatic restoration of polygon models［J］.ACM Transaction on Graphics，2005，24(4):1332－1352.

［14］丁展，陳志楊，張三元，等.基于Gauss Ball的二次曲面細分解與識別［J］.計算機輔助設計與圖形學學報，2007，1(19):31－36.

［15］白曉亮.逆向工程中混合CSG/B－rep模型重構(gòu)技術研究［D］.西安:西北工業(yè)大學，2005.

［16］MADSEN K，NIELSEN H B.Method for Non-Linear Least Squares Problems［M］.Copenhagen:Technical University of Denmark，1999.

A hole repairing algorithm based on surface feature recovery in triangular mesh model

HE Qiang，ZHANG Shu-sheng，BAI Xiao-liang

(The Key Laboratory of Contemporary Design and Integrated Manufacturing Technology，Ministry of Education China，Northwestern Polytechnical University，710072 Xi’an，China，hqcq@mail.nwpu.edu.cn)

In order to recover original shape of the holes in triangular meshes，a hole-repairing algorithm based on surface feature was proposed.First，the hole triangulation and subdivision were performed as coarse filling.Then，gauss ball was used to determine the hole neighbors and their surface types.Quadric surface equations were fixed by non linear least square.The hole neighbors of freeform were fitted by B spline surfaces.Finally，accurate positions of vertexes were acquired according to surface equations of the hole neighbors and the holerepairing was completed.The experimental results show that the presented hole repairing algorithm can not only fill the holes but also recover surface feature in the hole regions.

gauss ball;hole-neighbor;surface feature;hole-repairing

TP391

A

0367－6234(2011)11－0120－05

2010－01－10.

國家高技術研究發(fā)展計劃資助項目(2007AA04Z137);

國家自然科學基金資助項目(60573177).

賀 強(1985—)，男，博士研究生;

張樹生(1956—)，男，教授，博士生導師.

(編輯 張 紅)