基于支持向量機的Laplacian網(wǎng)格曲面孔洞修補算法

許 斌，李忠科，宋大虎

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

0 引 言

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

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

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

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

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

1 主要數(shù)學工具

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

1.1 網(wǎng)格曲面上的Laplacian坐標

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

Laplacian坐標根據(jù)權(quán)重的不同分別是同一Laplacian坐標和余切Laplacian坐標，本文采用余切Laplacian坐標

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

圖1 頂點Vi的局部片

在計算網(wǎng)格頂點的Laplacian坐標時，本文使用n×n階的Laplacian矩陣L，n為需要計算的網(wǎng)格區(qū)域內(nèi)的頂點數(shù)

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

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

1.2 最小二乘支持向量機

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

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

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

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

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

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

2 算法具體流程

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

2.1 孔洞邊緣檢測

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

步驟1 計算每個三角網(wǎng)格的中心以及它和三角形頂點間的聯(lián)系；

步驟2 通過三角網(wǎng)格的中心建立八叉樹；

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

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

2.2 孔洞邊界預處理

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

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

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

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

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

圖2 孔洞邊界預處理

2.3 建立孔洞特征平面

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

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

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

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

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

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

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

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

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

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

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

圖4 孔洞平面網(wǎng)格填充

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

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

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

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

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

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

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

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

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

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

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

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

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

3 算法評價及實驗結(jié)果

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

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

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

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

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

4 結(jié)束語

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

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

［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，F(xiàn)lyer 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）.［劉德平，于水晶.基于最小二乘支持向量機的三角網(wǎng)格修補算法 ［J］.計算機集成制造系統(tǒng)，2009，15 （9）：1868－187.］