基于支持向量機的Laplacian網格曲面孔洞修補算法

許 斌，李忠科，宋大虎

（第二炮兵工程大學 計算機教研室，陜西 西安710025）

0 引 言

使用掃描儀獲取模型表面的三維數據時，由于物體表面反射性、掃描設備的缺陷、模型表面的自遮擋等原因，點云數據通常含有孔洞，進而導致重構后的網格曲面也存在著孔洞?？锥吹拇嬖跇O大地影響了模型的完整性和顯示效果，在工業應用場合則會導致后續操作無法進行，因此，對這些孔洞進行修復是曲面網格模型數據處理的一個關鍵步驟。

從實際應用角度考慮，一個好的曲面修復算法應盡可能地滿足下面兩個條件：①修復曲面盡量接近原始真實形狀，與周圍網格有很好的光順連接；②算法有較強的魯棒性，即能夠處理各種類型孔洞?？傮w來說，孔洞的修復問題可以歸結為空間多邊形的三角剖分問題，很多學者都對這個問題進行了深入研究，下面對其中一些具有代表性的算法進行簡要介紹。文獻 ［1］提出了一種采用體數據場融合進行孔洞修補的算法，該算法首先迭代計算擴展體數據場的描述范圍，然后在此基礎上建立一個能夠描述整個孔洞區域及其周圍曲面的場函數來完成孔洞缺失數據的修復。文獻 ［2］使用RBF隱式曲面來對孔洞曲面進行修復，算法中首先計算孔洞邊界頂點的鄰域頂點，然后根據其幾何屬性插值生成孔洞區域的RBF隱式曲面，最后對該曲面進行網格化處理并與原始網格曲面進行融合得到孔洞修復后的網格曲面。上述幾種算法的優點是能夠修復較大的孔洞，但是對于缺失部分包含豐富幾何細的模型來說上述幾種修復方法就存在局限性，因為它們不能很好地恢復模型缺失部分的幾何細節。

文獻 ［3］提出了一種基于PDE （partial differential equations）方法的孔洞修復方法，該方法首先分析模型的幾何特征，然后采用正交梯度的方法來修復孔洞模型。Remacle等［4］提出了一種在孔洞邊界采樣點和缺失部分曲面定點之間建立調和映射的孔洞修補算法，最后通過最小化一個帶約束的全局能量函數來實現對缺失部分集合特征的修復。Park等［5］提出一種基于局部參數化的孔洞模型幾何和紋理的修復算法，該算法首先通過孔洞檢測算法檢測出孔洞，然后對待修復區域進行局部參數化，再通過自選擇或者人機交互在模型上選擇一塊相似的區域來填補孔洞，最后通過解Poisson方程將選擇區域變形移植到孔洞區域，從而完成孔洞模型幾何和紋理的修復。

近年來，神經網絡及支持向量機等智能化的機器學習方法被大量應用到網格曲面的孔洞修補中，用以推斷缺失部分曲面的幾何特征。文獻 ［6］使用混合徑向基函數模擬網格曲面缺失部分頂點坐標與原始網格曲面頂點坐標之間的函數關系，該算法能夠較好地恢復修補曲面的幾何特征；文獻 ［7］在此基礎上將遺傳算法、模擬退火法和BP網絡方法相結合，在一定程度上提高了算法的效率，但是生成網格的質量不高；文獻 ［8］中的方法也存在類似問題。

綜合上述各種算法可以看出，由于孔洞部分的幾何信息已經完全缺失，怎樣根據模型已有部分的幾何特點推斷出缺失部分的幾何信息，是網格曲面孔洞修復的根本難點所在。以神經網絡和支持向量機等機器學習系統為基礎的修補算法正是針對這一核心難點提出的。本文采用在統計學習理論的基礎上發展起來的新一代機器學習方法［9］支持向量機 （support vector machines，SVM）作為幾何信息回歸手段，結合當前非常流行的網格曲面Laplacian坐標，提出了一種全新的，針對三角網格曲面模型的孔洞修復算法，該算法突出的優點就是能最大限度地使修補后的曲面在幾何風格上與模型曲面原有部分保持一致。

1 主要數學工具

本文中使用的兩種核心數學工具是網格曲面的Laplacian坐標和最小二乘支持向量機，其作用分別是描述網格曲面的幾何細節，和依據孔洞周圍部分曲面的幾何信息對孔洞缺失部分曲面的幾何信息進行回歸推斷。

1.1 網格曲面上的Laplacian坐標

網格上的Laplacian坐標是一種描述網格模型幾何細節的數學工具，頂點Vi處的Laplacian坐標的計算公式如式（1），Vj是頂點Vi的局部片中的任意一點

Laplacian坐標根據權重的不同分別是同一Laplacian坐標和余切Laplacian坐標，本文采用余切Laplacian坐標

式 （2）、式 （3）是權重系數計算公式，采用式 （2）和式 （3）計算出來的就是頂點Vi處的余切Laplacian坐標，記做δi。如式 （4）所示余切Laplacian在數值上線性逼近網格頂點處的平局曲率，方向上逼近該點的主曲率法向，可以非常準確地描述網格頂點的局部片的幾何特性

圖1 頂點Vi的局部片

在計算網格頂點的Laplacian坐標時，本文使用n×n階的Laplacian矩陣L，n為需要計算的網格區域內的頂點數

式 （5）中的Lij的形式如式 （6）所示

式 （8）中的Vd是頂點坐標在某一坐標軸方向的分量，式 （7）中的Δd是Laplacian坐標在某一坐標方向的分量。

1.2 最小二乘支持向量機

最小二乘支持向量機 （least－squares support vector machines）是由Suykens提出的，這種方法采用最小二乘線性系統作為損失函數，突出優點是求解速度快，魯棒性強。最小二乘支持向量機在本文中的應用是一個典型的回歸問題，回歸的目標是最優化函數f（x）。本文給定訓練樣本集（如式 （9）所示）為孔洞周圍區域網格頂點的某一空間坐標和Laplacian坐標在同一坐標軸方向的分量

式中：d∈ ｛x ，y ，z｝，l——訓練樣本集中的樣本數量。平衡考慮回歸函數的復雜度和擬合誤差，依照結構風險最小化原理，將本文中的回歸問題描述為下面的約束優化問題

式 （11）中的w是權矢量，ei是松弛變量，γ為正則化參數，b為偏差量，為求解上述優化問題，需要建立Lagrange函數 （如式 （12）所示）

式中：a＝ （a1，a2，…，al）為 Lagrange乘子，在 KKT 最優化條件 （Karush－kuhn－Tuchker）下優化函數各變量的偏導數，并在鞍點處得到零值［10］（如式 （13）所示）

將式 （13）寫成矩形形式，并消去w和e可以得到式（14）

式中：Ωij＝K（xi，xj），求解得到式 （15）中的回歸函數

2 算法具體流程

本文算法分為以下4個部分：①孔洞邊界處理：包括孔洞邊界檢測與邊界邊預處理；②孔洞平面網格填充：包括建立孔洞特征平面、孔洞平面缺失采樣點估計與孔洞平面的網格生成；③三維缺失數據獲?。喊ㄓ柧殬颖炯牟杉妥钚《酥С窒蛄繖C回歸函數求解；④空間網格生成：包括三維空間中孔洞區域修補網格的生成。

2.1 孔洞邊緣檢測

在三角網格曲面上孔洞通常是由那些只包含于一個三角形的邊組成的，這些邊都屬于模型的邊界邊。為了能夠加快邊界檢測速度，本文采用層次空間剖分技術來對點云數據進行剖分。目前層次空間剖分技術主要包括八叉樹、四叉樹、BSP樹和K－D樹等，其中八叉樹表現出了優良的特性。因此，采用八叉樹來加快邊界檢測的速度。邊界檢測算法的具體步驟如下：

步驟1 計算每個三角網格的中心以及它和三角形頂點間的聯系；

步驟2 通過三角網格的中心建立八叉樹；

步驟3 搜索八叉樹，檢測與它相關的邊的拓撲連接性，如果某條邊只屬于一個三角形，則是邊界邊。

步驟4 通過頂點的不同幾何連通性分辨出不同的邊界回路，表示出不同的孔洞。

2.2 孔洞邊界預處理

由于在實際測量中孔洞區域的數據信息非常少，因此，重構后孔洞邊界處的三角片形狀可能會不理想，會出現一些狹長的三角片。這些三角片上的邊同其它邊界邊不協調，在孔洞填充時會造成不好的影響。為了降低這種影響，本文先對孔洞邊界進行預處理，將一些比較長的邊進行剖分，然后再對孔洞進行修補。

設vi（i＝1，2，…，n）為孔洞邊界上的n個頂點，ε為孔洞邊界上的邊，則孔洞邊界預處理可以分為三步。

步驟1 計算孔洞邊界上所有邊的平均長度Lave。

步驟2 遍歷所有的邊，如果存在某一條邊εij＝ （vi，vj）的長度Lεij＞k·Lave，k為比例系數 （本文中取k＝2），則計算邊εij的中點vnew。

步驟3 將vnew增加為孔洞的邊界點，同時刪除原來的邊界邊εij與三角片vvivj，增加邊界邊vivnew、vivnew和三角片vivvnew、vnewvvj，如圖2所示。最后迭代執行上一步操作，直到邊界邊沒有變化為止。

圖2 孔洞邊界預處理

2.3 建立孔洞特征平面

孔洞特征多邊形的平面化是通過將特征多邊形投影到特征平面上實現的，所以特征平面的選取是一個關鍵步驟，設Ve＝ ｛vi，i＝1，2，…，m｝為孔洞的邊界點集合，本文選取與Ve中的所有點的距離平方和最小的平面作為特征平面，這樣的特征平面的法向量如式 （16）所示

式中：N——特征平面法向量，nsum——孔洞周邊所有三角面片法矢量的向量和，這樣的法平面選取可以保證N和最后修復得到的孔洞部分填充曲面的主法向量 （所有新生三角面片的矢量和）偏差最小，進而保證對平面三角化后得到的修補網格頂點進行重新定位時Laplacian坐標的旋轉失真最小。

2.4 特征面平面上的孔洞填充

為了提高最終生成的三角網格質量，本文提出一種基于步長均值采樣的孔洞平面缺失數據估計算法來估計孔洞平面內缺失的數據，其基本思想是使用兩組步長相等的平行線對孔洞平面均值采樣，孔洞平面內的交點即為估計的數據。算法的具體步驟如下：

步驟1 在特征平面內求取孔洞特征多邊形的最小包圍盒；

步驟2 計算投影后邊界邊的平均邊長珚Lb；

步驟3 以珚Lb為步長，用兩組互成60°的平行線對特征多邊形的包圍盒均勻采樣，如圖3所示。

圖3 孔洞特征多邊形均值采樣

步驟4 采用多邊形內外點判斷方法對步驟3中的交點進行判斷，保留特征多邊形的內點；

步驟5 分別計算步驟4中所得估計點與邊界邊和邊界點間的距離le和lp，如果le＜0.2珚Lb或是lp＜0.2珚Lb，則刪除對應的估計點。

步驟6 最終所剩下的點則為所要計算的缺失數據點，如圖4（a）所示。

圖4 孔洞平面網格填充

在特征平面內通過采樣得到缺失數據點的坐標值之后，需要生成三角網格以方便后續的處理。直接在局部坐標系下將上一步求得的數據點與孔洞的邊界點一起進行Delaunay三角剖分，得到三角網格，完成孔洞平面的網格填充，填充后的網格如圖4（b）所示。

2.5 為最小二乘支持向量機選取學習樣本

迄今為止所有的孔洞修復算法都是依據已有網格曲面的幾何特征去推測缺失部分曲面的幾何特征，本文算法就是選取孔洞周圍的點及改點上Laplacian坐標作為最小二乘支持向量機的訓練樣本集，在實現過程中以孔洞邊界為中心不斷向外擴展，逐層采集相鄰幾層的三角形頂點，直到采集到滿足需要的層數為止。

假設Ve為孔洞邊界B 上的頂點集合，SVj（j＝0，1，…k）為孔洞外面第j層的頂點集合，樣本采集的具體算法如下：

步驟1 在邊界上取出一點vi，然后尋找與它相關的三角面片頂點，如果該頂點不在Ve中，則將它加入到孔洞的下一層頂點集SV1中；

步驟2 遍歷數組Ve中的每一個頂點，重復步驟1，這樣就采集到了由孔洞多邊形向外的第一層三角片頂點，并存放在頂點集SV1中；

步驟3 對頂點集SVj（j＝1，2…k－1）中的頂點，尋找與它相關的三角形的頂點，并將不在SVj和SVj－1中的點加入到頂點集SVj＋1中；

步驟4 迭代執行步驟3，直至采集達到所設定的層數為止，則SVj（j＝0，1，…k）中所有的點及其Laplacian坐標即構成了訓練樣本集，采樣點分布如圖5所示。

圖5 為最小二乘為支持向量機選取的樣本點集

2.6 填充部分頂點的重新定位

將第2.5小節中得到的訓練樣本集代入最小二乘支持向量機，按照2.2小節中的方法可以得到3個回歸函數：fx（x），fy（x），fz（x）。

設定Vnew為2.4小節中填充的新頂點集合。式 （17）中的V′nd是對Vnew中的頂點重新定位后的頂點坐標在某一坐標軸方向上的分量，重新定位的過程是對網格頂點坐標的3個分量分別進行求解，求解的是式 （18）中所示的n×n的線性系統，n為Vnew中頂點數

式 （18）中的LVn是在2.4小節中剖分得到的修補三角網格上，按照2.1小節中的方法建立的余切Laplacian矩陣，Δ′dn是重新定位后的修補區域網格頂點余切Laplacian坐標的某一分量構成的向量，是按照式 （19）構造的。式（19）中的fd（v′nd）是由回歸函數推測出的缺失區域網格頂點的Laplacian坐標的某一分量。

3 算法評價及實驗結果

為了檢驗本文算法的實際性能，以Visual C＋＋2010為開發平臺，結合OpenGL圖形庫實現了文中算法，并在head模型和happyvrip模型上進行了對比試驗，試驗結果如下：

圖6 中的模型就是head模型，在頭發部分人為制造了一個孔洞，作為待修補孔洞。試驗中選取的對比算法為傳統的直接修復算法和文獻 ［11］中的算法，之所以選擇文獻 ［11］算法是因為該算法和本文算法都是基于最小二乘支持向量機的，區別主要是文獻 ［11］算法直接將網格頂點的一個歐氏坐標分量做回歸函數的輸出，而本文引入了網格頂點的余切Laplacian坐標作為回歸函數的輸出，再以由回歸函數推斷得到的缺失部分網格曲面的Laplacian坐標為基礎建立線性系統求解得到缺失部分的網格頂點坐標。

圖6 原始head模型和制造孔洞后的head模型

圖7 中是使用原始的直接修復算法得到的修付效果，明顯比較粗糙。圖8中是使用文獻 ［11］中算法得到的修復效果，效果很好，修復區域的網格曲面質量也不錯。圖9是使用本文算法得到的修復效果，修復區域的網格質量也很好，同時可以看到相對于文獻 ［11］中的算法，使用本文算法的到的修補曲面與原始曲面的融合度明顯更高，效果與沒有制造孔洞前的原始模型曲面幾乎沒有差別。

圖10 是使用本文算法對happyvrip模型的修復效果，這個試驗的目的是檢驗本文算法在高曲率曲面上的修復效果，在具體實驗中，為了應付高曲率環境，加大了2.5小節中支持向量機的訓練樣本的選取范圍 （head試驗中的選取范圍是孔洞外的3層點，happyvrip試驗中為7層點），從試驗結果看本文算法對高曲率曲面的修復效果也是比較理想的。

4 結束語

本文算法的創新點在于引入了網格頂點的余切Laplacian坐標作為回歸函數的輸出，并通過求解線性系統的方式對修復區域的網格頂點進行二次定位，這樣做的優勢是發揮了Laplacian坐標對網格曲面幾何細節的刻畫能力，可以使修復后的網格曲面與原始曲面融合的更加自然。

本文算法的局限性在于對面積過大的孔洞的修復效果不理想，這也是目前所有的曲面修復算法面對的共同問題，對付這個問題的根本途徑是在構建修復區域網格時不僅要以孔洞周圍曲面的幾何特性作為依據，還要考慮模型的整體幾何特征，這也是本文下一步研究和改進的方向。

［1］Geuzaine C，Remacle J.Gmsh：A three－dimensional finite element mesh generator with built－in pre－and post－processing facilities ［J］.International Journal for Numerical Methods in Engineering，2009，79 （11）：1309－1331.

［2］Wu X，Wang M，Han B.An automatic hole－？lling algorithm for polygonal meshes ［J］.Computer－Aided Design and Applications，2008，5 （6）：889－899.

［3］Marchandise E，Piret C.CAD and mesh repair with radial basis functions［J］.Journal of Computational Physics，2012，231 （5）：2376－2387.

［4］Remacle J，Geuzaine C，Compère G，et al.High quality surface remeshing using harmonic maps ［J］.International Journal for Numerical Methods in Engineering，2010，83 （4）：403－425.

［5］Marchandise E，Carton de Wiart C，Vos W，et al.High quality surface remeshing using harmonic maps—Part II：Surfaces with high genus and of large aspect ratio ［J］.International Journal for Numerical Methods in Engineering，2011，86（11）：1303－1321.

［6］Wright G B，Flyer N，Yuen D.A hybrid radial basis function－pseudo－spectral method for thermal convection in a 3dspherical shell［J］.Chemistry Geophysics Geosystems，2010，11 （7）：1－18.

［7］CHEM Jie，GAO Chenghui，HE Bingwei.Automatic recognition of boundary features of non－closed triangulation model［J］.Machinery Design ＆ Manufacture，2011 （11）：89－94.

［8］Macdonald C B，Ruuth S J.The implicit closest point method for the numerical solution of partial differential equations on surfaces［J］.SIAM Journal on Scienti？c Computing，2009，31（6）：4330－4350.

［10］WU Dehui.Intelligent prediction model for surface roughness in milling ［J］.Computer Integrated Manufacturing Systems，2007，13 （6）：1137－1141 （in Chinese）.

［11］LIU Deping，YU Shuijing.Hole repairing in triangular meshes based on least－squares support vector machines ［J］.Computer Integrated Manufacturing Systems，2009，15 （9）：1868－187（in Chinese）.［劉德平，于水晶.基于最小二乘支持向量機的三角網格修補算法 ［J］.計算機集成制造系統，2009，15 （9）：1868－187.］