基于多向波前法的島嶼孔洞修補

李 松 馬聰聰 陸 帆 曹菁菁 周 勇

武漢理工大學物流工程學院，武漢，430063

0 引言

文物受自然環境、人為活動的影響，表面極易出現孔洞、殘缺。為減少傳統接觸式修復方法對文物的二次破壞，人們采用三維掃描系統對文物模型進行數據采集，以指導后面的修復。受測量技術與精度的限制、文物自身形狀與缺陷等因素的影響，被掃描物體的點云模型易出現孔洞?？锥雌茐牧宋奈锏耐暾裕瑢罄m數據的處理也產生較大影響，因此，對模型孔洞進行修補是逆向工程中十分重要的環節。

針對網格模型孔洞的修補，研究人員進行了大量研究。PERNOT等[1]提出曲率最小化原則，在孔洞周邊區域曲率最小化的前提下通過新增三角面修補孔洞。GIRSHICK等[2]采用非均勻有理B樣條(NURBS)曲面生成補丁、修補孔洞。ARGUDO等[3]提出雙諧波場的網格修補算法，可以較好地修補復雜孔洞和島嶼孔洞。謝倩茹等[4]采用波前法使孔洞區域不斷縮小，根據曲率特征細分網格。劉云華等[5]基于區域生長的方法對尖銳特征丟失的孔洞進行修補，獲得了較好的效果。晏海平等[6]通過建立徑向基函數的隱式曲面完成孔洞修復。袁天然等[7]通過構建拉格朗日方程對復雜孔洞進行修補。張立國等[8]基于徑向基函數對島嶼孔洞進行修補，提升修補速度。

對于島嶼孔洞的修補，上述方法均難以得到理想的修補效果：直接刪除島嶼面片將其變為普通封閉孔洞會造成特征信息丟失；在曲面擬合過程中將島嶼三角面作為控制點會引起網格相交。為了有效實現島嶼孔洞的修補，本文提出一種基于改進波前法的多向波前法，以有效提升島嶼類孔洞的修補效果。

1 多向波前法概述

三角網格模型表面的孔洞可以分為封閉孔洞、半封閉孔洞和島嶼孔洞[9]。本文算法的實驗模型為封閉的網格模型。封閉模型中不可能存在半封閉孔洞，因此下面僅對封閉孔洞與島嶼孔洞進行研究。

封閉孔洞是最常見的孔洞，其邊界連接成一個封閉環。島嶼孔洞是指在一個較大孔洞內有少量孤立的三角面孔洞。將島嶼孔洞中范圍最大的邊界定義為孔洞邊界，將孔洞邊界內的三角面集合稱為島嶼面片，將島嶼面片對應的邊界定義為島嶼邊界。圖1為三角網格模型封閉孔洞與島嶼孔洞示意圖。

圖1 網格模型孔洞類型Fig.1 Hole type of grid model

孔洞修補就是從孔洞邊界開始不斷向孔洞區域插入離散點與三角面，從而使孔洞不斷縮小。波前法可控性高、網格生成速度快[4]，是孔洞修補中常用的方法。該方法首先要檢測孔洞邊界，并以此作為初始波前，按照一定的規則在孔洞區域插入新增點與三角面片，每完成一次插入，孔洞邊界就會發生變化，需調整波前，循環上述過程直至完成孔洞修補。

圖2 多向波前法原理Fig.2 Principle of multi-directional advancing method

封閉孔洞與島嶼孔洞的區別在于波前的組成。傳統波前法針對只有一條邊界的孔洞。多向波前法原理如圖2所示，將孔洞邊界與島嶼邊界同時作為初始波前，插入離散點。以波前與插入的離散點為基礎，不斷向孔洞區域內新增三角面，逐漸縮小孔洞區域。孔洞邊界發生變化時，調整更新波前。重復上述過程直至修補結束。

2 算法詳細步驟

2.1 孔洞邊界識別

封閉的三角網格包含邊界邊與內部邊。每條邊都有鄰接三角形，只有1個鄰接三角形的為邊界邊，有2個鄰接三角形的為內部邊。邊界識別原理如圖3所示，遍歷模型的所有邊，找到邊界邊AB并令其為種子邊，從種子邊的2個頂點開始尋找鄰接邊，可以快速找到邊界邊。當邊界邊首尾相連、構成封閉的空間多邊形ABCDEFGHIJK時，完成一個完整的孔洞邊界搜索。

圖3 邊界邊識別原理Fig.3 Principle of boundary recognition

對于島嶼孔洞，由于在大的孔洞內部存在島嶼面片，因此會搜索出多個邊界，利用孔洞邊界頂點、根據最小二乘法[10]擬合最小二乘平面，將孔洞邊界投影到擬合的平面上，組成目標區域。將孔洞邊界頂點的一階鄰域也投影到擬合的最小二乘平面上，若投影點位于目標區域內，則為島嶼邊界，反之，則為孔洞邊界。

2.2 孔洞邊界預處理

孔洞邊界邊的長度相差較大，會影響離散點的插入，而且后續網格化的過程會產生狹長的三角面片，降低網格質量。因此在新增離散點之前，需要對孔洞邊界進行預處理，使其均勻化。

計算所有孔洞邊長度及其平均值lavg[11]，令閾值為1.5lavg，孔洞頂點以Vi(i=1，2，…)表示，當某條邊ViVi+1的長度l>1.5lavg時，計算邊ViVi+1的中點Vnew。刪除原有孔洞邊和對應鄰接三角形，孔洞邊更新為ViVnew和Vi+1Vnew，孔洞三角形更新為ViVnewVj以及VnewVi+1Vj。

2.3 離散點三角面插入

對島嶼孔洞進行修補時，島嶼面片的邊界最小角θmin往往會大于180°，傳統波前法并未對此進行考慮，會生成狹長三角面，降低網格質量。因此，本文對傳統波前法中新增點插入規則進行改進，改進的插入規則如下：

(1)以孔洞邊界和島嶼邊界為初始波前。

(2)計算兩兩孔洞邊的夾角，得到最小夾角∠Vi-1ViVi+1。

(3)將夾角最小的孔洞頂點作為初始填充位置，根據θmin的大小設置不同的頂點與三角面插入準則，具體準則如下：①0°<θmin<30°時(圖4a)，θmin屬于狹小角，直接連接兩頂點會產生狹長的三角面單元，此時應合并該三角面，刪除原有頂點Vi+1、Vi-1，計算Vi+1Vi-1的中點Vnew，如圖4b所示，以Vnew為新頂點，重新建立三角面，如圖4c所示。②30°<θmin<80°時(圖5a)，直接連接點Vi-1、Vi+1，新增一個三角面。③80°<θmin<150°時(圖5b)，在θmin的角平分線方向上，在孔洞邊Vi-1Vi和ViVi+1邊長的均值處添加新點Vnew，再連接Vi-1Vnew以及Vi+1Vnew，新增2個三角面。④150°<θmin<210°時(圖5c)，在θmin的三等分線方向上，在孔洞邊Vi-1Vi和ViVi+1邊長的均值處添加新點Vnew1、Vnew2，連接Vi-1Vnew1、Vnew1Vnew2、Vi+1Vnew2，新增3個三角面。⑤210°<θmin<270°以及270°<θmin<360°時，按上述方式分別新增4個和5個三角面。

(4)按步驟(3)向孔洞區域內插入頂點與三角面后，更新波前。

(5)重復步驟(2)～步驟(4)，當邊界內頂點數為3時，結束網格生長，完成孔洞修補。

圖4 夾角過小時剖分準則Fig.4 Split criterion under relatively small angle

圖5 最小角度原則新增點規則Fig.5 New point rule of minimum angle principle

2.4 新增頂點和新增面檢驗

使用波前法常會出現狹長三角形以及由不合理新增點引起的錯誤三角形，因此需要對新增頂點和新增三角面進行檢驗。

三維空間中的合理性檢驗比較復雜，可將空間三角形投影至二維平面完成檢驗。創建過點Vi且以n0=(n1+n2)/2(n1、n2為點Vi的兩鄰接三角形的法向量)為法矢的平面，將孔洞三角形投影至該平面，形成目標區域。若新增點的投影不在目標區域內，說明新增點出現在角平分線的反向延長線上，則該點為不合理點，應將其沿角平分線方向移動到頂點Vi的對稱位置。

Delaunay三角剖分算法將三角網格中的最小角最大化，以預防狹長三角面的產生，提高網格質量[12]。平面Delaunay檢驗利用“空圓準則”即四點不共圓準則[13]，如果2個三角形不滿足Delaunay準則，則根據“最大化最小角”原則將2個三角形的對角線進行交換，使其滿足Delaunay準則，從而避免產生狹長三角形。

2.5 基于法矢與曲率的頂點調整

2.5.1法矢與曲率計算

為獲得更好的曲面孔洞修補效果，利用頂點法矢和曲率調整新增頂點位置，使新增三角面最大程度地擬合孔洞原有表面特征。

三角網格模型頂點Vi的鄰接三角形如圖6所示，其中，Ni為Vi的單位法矢，fik為Vi的一個鄰接三角面，Nik為fik的單位法矢，αk是fik在Vi處的頂角，ei,j和ei,j+1是頂角αk的2條邊，gk為fik的質心，Gk為Vi到gk的距離。

圖6 頂點Vi的鄰接三角形Fig.6 Adjacent triangles of vertex Vi

三角網格模型中某三角面片fik的單位法矢為

(1)

頂點法矢為該點所有鄰接三角面片法矢的加權平均，若頂點有n個鄰接三角面，令三角面fik的權重為wk，則頂點Vi的單位法矢為

(2)

權重wk由三角形形狀因子λ和質心距Gk計算得到，λ的計算公式為

(3)

其中，a、b、c為三角形的邊長。λ(0,1]，λ越大，三角形的規整程度越高[14]。λ與質心距Gk可以較好地評價三角形幾何特性，將式(2)變形為

(4)

利用式(4)對點Vi的單位法矢進行計算，為了估算結果的準確性，在孔洞頂點的法矢計算之前，先將缺失部分補齊。

(5)

將新增點Vnew替換點Vj，點Vi處沿方向T的法曲率為

(6)

2.5.2頂點調整

由Vi指向新增點Vnew構成向量nnew，nnew的單位向量為nb，nb與頂點法矢Ni共同確定平面Π。將新增點在平面Π上進行調整，得到V′new，V′new與Vj組成的向量的單位向量為nB，如圖7所示。

圖7 新增頂點調整示意圖Fig.7 New vertex adjustment diagrammatic sketch

若不進行頂點調整，則利用波前法求得新增點坐標：

Vnew=Vi+(|Vi-1Vi|+|ViVi+1|)nb/2

(7)

根據曲率ki(T)對頂點進行調整。ki(T)>0時，曲面在點Vi沿T方向處為凸，應將nb沿法向方向旋轉θ角；ki(T)<0時，曲面在點Vi沿T方向處為凹，應將nb沿法向的反方向旋轉θ角；ki(T)=0時，無需調整。聯立式(4)、式(5)得調整角度θ以及nB：

(8)

θ=arccos(nB·nb)

(9)

調整后的坐標Vnew=Vi+(|Vi-1Vi|+|ViVi+1|)nB/2，將所有新增點按該過程完成調整。

本文算法在修補孔洞過程中基于波前法，以孔洞邊界與島嶼邊界為基礎增加頂點與三角面，并基于法矢與曲率進行自適應調整，使新增區域盡可能擬合缺失部分的原有幾何形狀，使孔洞區域與非孔洞區域能光滑銜接，不斷更新波前信息，有效減小修補痕跡。

3 實驗結果分析

為驗證本文方法的有效性，選擇3組網格數據進行實驗。實驗環境為VS2013，采用OpenGL庫驗證算法。首先驗證該算法對普通封閉孔洞的修補效果，采用帶有孔洞的兔子模型進行實驗。圖8a所示模型的孔洞為狹長孔洞，利用文獻[1]算法和文獻[4]算法對其進行修補，得到的效果圖分別為圖8b、圖8c，利用本文算法得到的結果如圖8d所示，本文算法明顯優于前兩種算法。

手模型1中的孔洞為島嶼孔洞，圖9a所示為孔洞模型，圖9b為文獻[1]算法的修補效果圖，利用文獻[4]算法得到的效果如圖9c所示，本文算法對島嶼孔洞修補后的效果如圖9d所示。對比發現，利用本文所提算法得到的新增面片與非孔洞區域過渡自然，保留了原有曲面特征，修補效果明顯優于前兩種算法。

圖9 手模型1孔洞不同算法修補效果Fig.9 Repair effect of different algorithms for Hand 1 model

手模型2孔洞的實驗結果如圖10所示，圖10a所示為島嶼孔洞模型，可以看出該孔洞既包含大規模島嶼面片，也有較小的島嶼面片。利用文獻[1]算法和文獻[4]算法得到的修補效果分別如圖10b、圖10c所示。圖10d所示為本文算法的修補結果，多向波前法在較大程度上保留了原有特征，過渡自然，無明顯修補痕跡，效果較前兩種算法更理想。

圖10 手模型2孔洞不同算法修補效果Fig.10 Repair effect of different algorithms for Hand 2 model

采用形狀因子定量評價修補效果，由2.5.1節可知，形狀因子越接近于1，三角形面片越規整。表1所示為運用不同算法進行修補得到的新增三角網格平均形狀因子，利用本文算法進行修補得到的平均因子超過0.85，在3種算法中效果最好，網格質量更高。

表1 新增網格平均形狀因子

由上述3組實驗可知，本文算法修補效果明顯優于其他算法，采用本文所提多向波前法對青銅古劍模型(圖11)進行修復。使用手持式三維掃描儀Artec Space Spider獲得文物表面點云數據，對數據進行濾波、精簡等預處理后，采用本文提出的多向波前法對孔洞進行修補?？锥粗饕霈F在劍柄與尾端部分，圖12所示為劍柄部分孔洞修補前后的效果，劍尾部分孔洞修補前后的效果如圖13所示，可以看出，利用本文算法進行孔洞修補后，文物表面能夠與非孔洞區域自然過渡，很好地還原文物表面特征，并且較為復雜的島嶼與大范圍孔洞也能得到較為理想的修補效果。

圖12 劍柄部分孔洞修復前后對比Fig.12 Comparison of the sword hilt before and after repair

圖13 劍尾部分孔洞修復前后對比Fig.13 Comparison of the sword tail before and after repair

使用聯泰Lite600HD工業級3D打印機，采用光固化成形技術，以白色光敏樹脂為材料，按1∶1的比例打印修補后的古劍模型，如圖14所示?？梢钥吹?，經本文算法修補后進行打印的3D模型表面光滑，孔洞區域與非孔洞區域銜接光順，也從側面反映本文算法的優勢。

圖14 3D打印實物模型Fig.14 Physical model of 3D printing

4 結論

本文針對島嶼孔洞問題提出了多向波前法。首先識別孔洞邊界并對其進行預處理，改進傳統波前規則，完成離散點與三角面片的插入。對新增頂點與三角面片進行檢驗與優化，減少狹長三角面片。采用法矢與曲率信息完成對新增頂點的調整，使新增三角面最大程度地擬合原有表面特征，實現島嶼孔洞的修補。采用多組網格數據進行實驗并對文物模型進行修復，結果表明本文算法生成網格質量較高，可以有效完成島嶼孔洞修補?？锥葱扪a算法的適用范圍受原有模型的曲面復雜程度影響較大，當孔洞原有區域形狀不規整、不均勻且孔洞部分過大時，利用孔洞修補算法難以得到理想的修補效果。對于曲面復雜度不同的孔洞，本文算法的適用范圍將作為下一步的研究重點。