張 琦， 藺素珍， 白佳璐， 王棟娟

(中北大學(xué) 大數(shù)據(jù)學(xué)院， 山西 太原 030051)

基于回溯雙向波前法虛擬修補點云孔洞

張 琦， 藺素珍， 白佳璐， 王棟娟

(中北大學(xué) 大數(shù)據(jù)學(xué)院， 山西 太原 030051)

針對虛擬修補邊界點凹凸不平且分布不均的點云孔洞效果欠佳的問題， 提出了回溯雙向波前法虛擬修補包含復(fù)雜邊界的孔洞. 以具有復(fù)雜邊界孔洞的青銅器模型為例： 首先， 三角網(wǎng)格化青銅器點云數(shù)據(jù)并依據(jù)網(wǎng)格化結(jié)果提取出孔洞邊界， 通過比較邊界點集的曲率波動幅度以去除偽孔洞邊界； 其次， 以孔洞邊界點集的回溯結(jié)果為初始點集， 并結(jié)合向量叉積對初始點集進(jìn)行凹凸性分類， 再對凹、 凸點分別采用正、 反向波前法逐圈新增點集直至補全； 最后， 利用最小二乘法擬合曲面平滑新增點集， 獲得最終修補結(jié)果. 實驗結(jié)果表明該方法與劃分子洞波前法、 曲線流收縮法相比， 結(jié)構(gòu)相似性的平均值分別提高了81.14%和93.8%， 且曲率差異性更低， 修補網(wǎng)格的頂點密度與原始網(wǎng)格更相近且過渡自然， 能有效修補復(fù)雜邊界孔洞.

虛擬修補； 點云孔洞； 回溯雙向波前法； 三角網(wǎng)格

0 引 言

利用計算機虛擬復(fù)原技術(shù)輔助點云三維模型修補， 是提高逆向工程中三維重建效果和效率的重要手段. 文物這類特殊的三維模型因自然腐蝕、 坍塌等原因所致， 其出土?xí)r往往存在大小不一、 形狀各異的復(fù)雜邊界孔洞. 對這些孔洞進(jìn)行虛擬精準(zhǔn)修補， 既是后續(xù)紋理映射和三維效果重建的基礎(chǔ)[1]， 也可為人工修補提供詳細(xì)的操作依據(jù)和評價標(biāo)準(zhǔn)， 因此文物點云孔洞修補已成為其修補的關(guān)鍵技術(shù)之一.

目前， 關(guān)于文物點云孔洞修補的研究較少. 與之相關(guān)的主要有隧道[2]、 工業(yè)零件[3]以及動物模型[4]等的虛擬修補， 這些研究均證明將點云三角形網(wǎng)格化后再進(jìn)行修補是行之有效的手段. 而基于三角形網(wǎng)格進(jìn)行孔洞修補包括提取孔洞邊界和新增點集. 依據(jù)邊界邊僅存在于一個三角面片的特性提取孔洞邊界被認(rèn)為是最為簡單有效的方法[5,6]， 但鑒于其僅適用于提取全封閉模型的孔洞邊界， 新近有學(xué)者根據(jù)三維點云二維投影后的8連通域?qū)Π敕忾]對象提取邊界， 雖成功提取到孔洞邊界， 卻對存在投影重疊的情況無能為力[7]. 至于新增點集的方法， 依據(jù)待修補點云孔洞邊界的特點大致可分為兩類： 一是針對邊界點全凹的孔洞， 多采用波前法[8]和基曲面網(wǎng)格添加法[9]等， 這些方法不僅能完整修補孔洞還能實現(xiàn)光滑連續(xù)； 二是對邊界點凹凸共存且分布均勻的孔洞， 采用劃分子洞波前法[10]和曲線流收縮法[11]等獲得較好的修補效果. 與這些修補對象相比， 文物不僅大多呈半封閉形狀而且在不同方向上投影大多存在點云重疊現(xiàn)象， 另外， 其孔洞多因腐蝕所致， 孔洞邊界點往往凹凸不平且分布不均勻， 用以上方法進(jìn)行修補難以獲得理想效果.

為此， 本文以包含孔洞的青銅器點云為對象， 在提取邊界后， 根據(jù)曲率波動幅度去除口部邊緣等偽孔洞邊界， 通過回溯到邊界點集外圈建立新增點集以解決因邊界點分布不均勻所致的修補誤差問題， 然后對凹點采用正向波前法、 對凸點采用反向波前法逐圈新增點集， 在補全孔洞之后用原始的孔洞邊界替換新增的第一圈點集， 再利用最小二乘法擬合曲面對剩余的新增點集進(jìn)行微調(diào)， 即得到最終修補結(jié)果.

1 復(fù)雜邊界點云孔洞修補總體思路

針對邊界點凹凸不平且分布不均勻的點云孔洞， 本文提出的回溯雙向波前法虛擬修補孔洞思路如圖 1 所示， 其主要步驟為：

1) 對采集到的青銅器點云三角形網(wǎng)格化.

2) 根據(jù)邊界邊僅存在于一個三角形面片中的特點提取出孔洞邊界.

3) 對孔洞邊界和偽孔洞邊界的曲率波動幅度進(jìn)行比較， 去除偽孔洞邊界.

4) 回溯到孔洞邊界的外圈點作為新增點集的初始點集， 對凹、 凸點分別用正、 反向波前法建立新增點集， 反復(fù)進(jìn)行直至孔洞補全.

5) 用孔洞原邊界替換新增點集第一圈.

6) 用最小二乘法擬合曲面微調(diào)剩余新增點集， 得到最終修補結(jié)果.

圖 1 孔洞虛擬修補思路圖Fig.1 The schematic diagram of holes virtual repaired

2 點云孔洞修補的回溯雙向波前法

2.1 孔洞邊界自動提取

Delaunay三角剖分的三角網(wǎng)格具有空圓特性和最小內(nèi)角最大化特性， 廣泛應(yīng)用于平面和空間的三角剖分[12]. 本文采用Delaunay三角形網(wǎng)格化法建立網(wǎng)格， 根據(jù)其邊界邊僅存在于一個三角形中的特性[5]提取孔洞邊界. 由于青銅器多屬半封閉模型， 所以提取出的邊界還包含偽孔洞邊界， 即原模型中的口部邊緣等非閉合處. 鑒于偽孔洞邊界多為器皿自身邊緣， 其曲率連續(xù)性較好， 而腐蝕所致孔洞邊界曲率連續(xù)性往往較差， 所以利用曲率分布特性可區(qū)分孔洞邊界和偽孔洞邊界點集. 去除偽孔洞邊界的具體過程為：

2.1.1 計算各點曲率并建立其分布圖.

圖 2 曲率分布折線圖Fig.2 The polyline diagram of the curvature distribution

不同孔洞邊界的曲率分布如圖 2 所示. 其中， 虛線與實線分別是偽孔洞邊界和實際孔洞邊界的曲率分布， 可以看出前者的波動幅度遠(yuǎn)小于后者.

2.1.2 計算曲率波動幅度D.

2.1.3 去除偽孔洞邊界.

去除曲率波動幅度最小的偽孔洞邊界， 即得到孔洞邊界點集L=(p1,p2,…,pi,…,pn). 需要說明的是此方法僅能去除半封閉模型中邊界曲率連續(xù)性較好的偽孔洞邊界.

2.2 回溯雙向波前法建立新增點集

2.2.1 回溯孔洞邊界

青銅器點云孔洞邊界點分布不均勻， 直接從此處開始建立新增點集會造成點集過密， 不利于保持與原網(wǎng)格的相似性， 進(jìn)而會影響后續(xù)的紋理映射和三維效果重建. 本文利用網(wǎng)格的拓?fù)浣Y(jié)構(gòu)從孔洞邊界進(jìn)行回溯， 以其外圈點作為新增點集的初始點集.

2.2.2 基于雙向波前法建立新增點集

傳統(tǒng)的波前法雖是一種簡單有效的修補方法， 但因青銅器孔洞邊界點凹凸共存， 若對凸區(qū)用該方法建立新增點集， 則會增加到原始點云區(qū)域. 所以， 本文分別處理孔洞邊界的凹凸區(qū)域， 即凹區(qū)采用傳統(tǒng)的正向波前法、 凸區(qū)采用反向波前法修補. 正向波前法如文獻(xiàn)[8]所述， 包括建立和更新波前點集兩部分. 反向波前法的具體實現(xiàn)方法如下：

1) 以已提取的邊界點集為波前點集， 計算波前點集每一個頂點pi相鄰的邊界邊〈ei,ei+1〉間夾角θi.

2) 以最小夾角(minθi)對應(yīng)點為初始點， 根據(jù)式(3)判斷其凹凸性， 其中·為點積. 若d<0為凸點， 做反向處理； 反之則為凹點， 正向處理.

3) 基于pi對應(yīng)的夾角θi的大小， 分為3類分別在ei和ei+1所在平面上新增點.θi≤90°時， 不添加新的點； 90°<θi≤135°時， 在角平分線的反向方向上以式(4)新增點pnew；θi>135°時， 則在角3等分線的反方向上以式(5)， 式(6)新增點pnew1,pnew2. 然后按順時針方向依次添加新的點， 并更新波前點集.

式中： |pi+1-pi|是向量(pi+1,pi)的模， 其余同理.

2.2.3 利用最小二乘法擬合曲面微調(diào)新增點集

為保證所修補的孔洞能光滑連續(xù)， 在補全后還需要利用孔洞邊界外圈的局部區(qū)域基于最小二乘法擬合曲面微調(diào)新增點集[13]. 首先， 用孔洞邊界替換掉新增的第一圈點， 然后對剩余的每一圈都利用其外圈點集來近似擬合一個光滑的完整曲面， 實際上就是調(diào)整每一個點的z坐標(biāo)， 計算方法為

式中：x和y為橫縱坐標(biāo)， 系數(shù)a1…a6通過式(8)計算獲得.

式中：n為擬合曲面所選點的個數(shù)；ω(xi)是權(quán)函數(shù)， 取值為1；p(x)是基數(shù)，pT(x)=[1,x,y,x2,y2,xy].

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

3.1 定性評價

圖 3(a) 為從山西省博物館采集到的包含孔洞的青銅器數(shù)據(jù)(No.1)， 該數(shù)據(jù)有247 521個頂點， 包含兩個孔洞， 圖 3(b)～圖 3(e)為采用本文方法修補前后的整體效果圖和最終的模型三維圖. 總體看， 本文方法不但恢復(fù)了缺失的孔洞部分， 而且修補區(qū)域能與原始區(qū)域光滑連接.

圖 3 No.1青銅器孔洞修補效果圖Fig.3 The result diagram of No.1 bronze holes repaired

文獻(xiàn)[10]曾采用劃分子孔洞波前法并通過最小二乘法擬合曲面微調(diào)獲得比傳統(tǒng)波前法更好的孔洞修補效果,文獻(xiàn)[11]則利用曲線流收縮來修補孔洞. 為便于比較， 現(xiàn)將利用文獻(xiàn)[10]、 文獻(xiàn)[11]和本文方法修補孔洞1和孔洞2的結(jié)果分別放大， 如圖 4 所示， 其中圖4(a)～圖4(c)為孔洞1的修補結(jié)果網(wǎng)格圖， 圖4(d) 和圖4(e)為本文方法修補前后模型圖； 圖4(f)～圖4(h)為孔洞2的修補結(jié)果網(wǎng)格圖， 圖4(i)和圖4(j)為本文方法修補前后模型圖. 圖4(a), 圖4(f), 圖4(b), 圖4(g)和圖4(c), 圖4(h)分別是用文獻(xiàn)[10]、 文獻(xiàn)[11]和本文方法的修補結(jié)果圖， 從中可以看出本文方法修補網(wǎng)格更相似于原始網(wǎng)格， 說明其光滑連續(xù)性更好.

圖 4 No.1青銅器孔洞修補結(jié)果細(xì)節(jié)圖Fig.4 The results detail diagram of No.1 bronze holes repaired

圖 5 是另外兩組青銅器點云孔洞及其修補結(jié)果， 其中圖5(a)， 圖5(g)為原始數(shù)據(jù)的三角形網(wǎng)格圖， 圖5(b)， 圖5(h)是文獻(xiàn)[10]的修補結(jié)果， 圖5(c)， 圖5(i)是文獻(xiàn)[11]的修補結(jié)果， 圖5(d)， 圖5(j)是本文方法的修補結(jié)果， 圖5(e)， 圖5(k)是修補前模型圖， 圖5(f)， 圖5(l)是本文方法修補的模型結(jié)果. 盡管兩組數(shù)據(jù)中孔洞所在面的曲率起伏不同， 但采用本文方法修補結(jié)果與原網(wǎng)格的相似性都優(yōu)于文獻(xiàn)[10]和文獻(xiàn)[11].

圖 5 No.2和No.3青銅器孔洞修補結(jié)果圖Fig.5 The result diagram of No.2 and No.3 bronze holes repaired

3.2 定量評價

修補區(qū)域與原始區(qū)域網(wǎng)格越相似越有利于后續(xù)的紋理映射等三維重建； 修補區(qū)域與原始區(qū)域曲率連續(xù)性越好表明曲面越光滑， 綜合兩者可衡量點云孔洞修補質(zhì)量的好壞.

3.2.1 結(jié)構(gòu)相似性(structural similarity, SSIM)

結(jié)構(gòu)相似性是常用的衡量三維模型空間相似性的指標(biāo)， 其值越大表示結(jié)構(gòu)越相似. 表 1 為利用文獻(xiàn)[10]、 文獻(xiàn)[11]方法和本文方法對修補區(qū)域和原始區(qū)域的結(jié)構(gòu)相似性統(tǒng)計結(jié)果， 從表 1 中可以看出， 以上3組孔洞修補結(jié)果的SSIM指標(biāo)， 本文方法都要高于文獻(xiàn)[10]和文獻(xiàn)[11]， 平均值分別高出81.14% 和93.8%， 說明本文方法具有更準(zhǔn)確的修補結(jié)果.

表 1 3組孔洞數(shù)據(jù)的SSIM統(tǒng)計結(jié)果

3.2.2 曲率連續(xù)

曲率連續(xù)可以測量多個曲面在相交處的平滑程度. 曲率分布圖可以直觀體現(xiàn)曲率的連續(xù)性， 本文借助通用的Imageware中三角形網(wǎng)格曲率分布圖對曲面的連續(xù)性進(jìn)行評價. 圖 6 為Imageware中3組數(shù)據(jù)的三角形網(wǎng)格曲率分布圖， 圖中各深淺色分別表示凹、 凸和較平坦區(qū)域， 顏色越接近說明連續(xù)性越好. 其中圖6(a)， 圖 6(d)， 圖6(g)和圖6(j)是用文獻(xiàn)[10]方法的修補結(jié)果曲率分布圖， 圖6(b)，圖6(e)， 圖6(h) 和圖6(k)是用文獻(xiàn)[11]方法的修補結(jié)果曲率分布圖， 圖6(c)， 圖6(f)， 圖6(i)和圖6(l)是本文方法的修補結(jié)果曲率分布， 畫圈區(qū)域為原孔洞所在區(qū)域. 直觀上看， 孔洞周圍區(qū)域以凹凸為主， 文獻(xiàn)[10] 修補結(jié)果較平坦， 文獻(xiàn)[11]修補結(jié)果凸區(qū)較多， 而本文修補結(jié)果與周圍區(qū)域顏色分布更為相似.

圖 6 3組孔洞曲率分布圖Fig.6 The curvature distributed diagram of three group holes

表 2 3種區(qū)域所占百分比%

為便于定量描述， 采用所修補孔洞與其周圍10圈(實物長度約0.5 cm)區(qū)域中凸、 凹、 平3種區(qū)域的分布百分比進(jìn)行比較， 差異范圍越小與原始區(qū)域凹凸性越相似， 如表 2 所示： 在3組對象中， 本文方法與修補前的百分比差異范圍在0～4.4之間， 而文獻(xiàn)[10]的差異在0.8～34.1之間， 文獻(xiàn)[11]的差異在2.6～15.3之間， 表明本文方法修補后能更好地保持孔洞與周圍區(qū)域的凹凸一致性. 綜上， 在用本文方法修補后的曲率分布圖中孔洞均與孔洞周圍區(qū)域的曲率圖呈光滑連接， 而用文獻(xiàn)[10]、 文獻(xiàn)[11]修補后的孔洞區(qū)域邊界處可以保持光滑連接但越靠近孔洞中心位置曲率分布越不連續(xù). 由此可見， 本文方法修補結(jié)果曲率連續(xù)性更好， 保證了孔洞與周圍區(qū)域的光滑連接.

4 結(jié) 論

本文的回溯雙向波前法基于半封閉模型， 在三角形網(wǎng)格基礎(chǔ)上有效地提取出孔洞邊界并實現(xiàn)了復(fù)雜邊界孔洞的修補， 具有較高的結(jié)構(gòu)相似性和曲率連續(xù)性， 其主要特點體現(xiàn)在： ① 利用回溯雙向波前法可有效修補青銅器等半封閉模型中邊界點凹凸不平且分布不均勻的點云孔洞； ② 依據(jù)孔洞邊界點曲率波動幅度能有效去除偽孔洞邊界； ③ 將孔洞邊界回溯到其外圈， 有利于減小因孔洞邊界點分布不均勻給新增點集帶來的誤差； ④ 使用雙向波前法對兼具凹凸點的邊界建立新增點集簡單實用； 用最小二乘法微調(diào)修補部分， 可使得修補區(qū)域與周圍連接更光滑.

值得指出的是本文方法的局限性是修補尖銳部分及非閉合狀的孔洞效果欠佳， 下一步將以此作為重點進(jìn)行研究.

[1] 姜翰青, 王博勝, 章國鋒， 等. 面向復(fù)雜三維場景的高質(zhì)量紋理映射[J]. 計算機學(xué)報, 2015, 38(12)： 2349-2360.

Jiang Hanqing, Wang Bosheng, Zhang Guofeng, et al. High-quality texture mapping for complex 3D scenes[J]. Chinese journal of Computers, 2015, 38(12)： 2349-2360. (in Chinese)

[2] 錢伯至, 藍(lán)秋萍. 隧道三維點云孔洞修復(fù)方法[J]. 測繪工程, 2017, 26(3)： 46-50， 55.

Qian Bozhi, Lan Qiuping. Repairing tunnel 3D point cloud hole[J]. Engineering of Surveying and Mapping, 2017, 26(3)： 46-50， 55. (in Chinese)

[3] 王欽瑞. 三角網(wǎng)格模型特征線提取與孔洞修補方法研究[D]. 大連： 大連理工大學(xué), 2016.

[4] Fortes M A, Gonzalez P, Palomares A, et al. Filling holes with shape preserving conditions[J]. Mathematics and Computers in Simulation, 2015, 118： 198-212.

[5] Nathan L F, Yaoguo L. Automatic boundary extraction from magnetic field data using triangular meshes[J]. Geophysics, 2016, 81(3)： 147-160.

[6] Yann Q, Claire L. Filling holes in digitized point cloud using a morphing-based approach to preserve volume characteristics[J]. The International Journal of Advanced Manufacturing Technology, 2015, 81： 411-421.

[7] Van S N, Ha T M, Thanh N T. Filling holes on the surface of 3D point clouds based on tangent plane of hole boundary points[J]. Symposium on Information and Communication Technology(ACM), 2016： 331-338.

[8] Yasushi I, Alan M S, Anil K E, et al. Parallel unstructured mesh generation by an advancing front method[J]. Mathematics and Computers in Simulation, 2007, 75(5-6)： 200-209.

[9] 劉震, 王艷賓, 白麗麗, 等. 曲面細(xì)節(jié)特征保持的三維模型孔洞修復(fù)算法[J]. 計算機輔助設(shè)計與圖形學(xué)學(xué)報, 2016, 28(12)： 2052-2059.

Liu Zhen, Wang Yanbin, Bai Lili, et al. Detail-preserving hole-filling for complex 3D models[J]. Journal of Computer-Aided Design & Computer Graphics, 2016, 28(12)： 2052-2059. (in Chinese)

[10] Jun Y. A piecewise hole filling algorithm in reverse engineering[J]. Computer-Aided Design, 2005, 37(2)： 263-270.

[11] Altantsetseg E, Khorloo O, Matsuyama K, et al. Complex hole-filling algorithm for 3D model[C]. Computer Graphics International Conference, 2017.

[12] 魏瀟然, 耿國華, 張雨禾. 幾何信息預(yù)測的三角網(wǎng)格模型拓?fù)鋲嚎s[J]. 西安電子科技大學(xué)學(xué)報, 2015, 45(5)： 194-199.

Wei Xiaoran, Geng Guohua, Zhang Yuhe. Connectivity compression of triangle meshes based on geometric parameter predict[J]. Journal of Xidian University, 2015, 45(5)： 194-199. (in Chinese)

[13] Wei Z, Feng F, Wei W. Non-linear partial least squares response surface method for structural reliability analysis[J]. Reliability Engineering and System Safety, 2017, 161： 69-77.

HoleRepairinginPointCloudBasedonBacktrackingBidirectionalAdvancingFrontMethod

ZHANG Qi， LIN Suzhen， BAI Jialu， WANG Dongjuan

(School of Computer and Control Engineering， North University of China， Taiyuan 030051， China)

A method utilized to virtual repair the complex boundary hole is proposed based on the backtracking bidirectional wavefront method which aim to solve the problem that the repair effect is poor due to the even and uneven distribution of the boundary of holes. The approach regard bronzes as experimental objects. Firstly, the bronze cloud data are Disposed of triangular mesh and the boundary of hole extracted in mesh. Then, the boundary of the pseudo-hole is removed by comparing the curvature fluctuation range of the boundary point set. Secondly, the point of backtracking result is the initial point set and combine the method of vector cross product to classify the concave and convex points, which, respectively, are operated by positive and negative wavefront method before adding a new set of points until the completion. Then the least squares method is used to fit the surface to smooth the added point set and then the final patching result is obtained. The experimental results show that the average data of our approach is increased respectively by 81.14%、93.8% in SSIM and has lower curvature differences compared with simple hole divided and curves flow shorted. The hole-filling meshes obtained by our method interpolate the shape and have consistent mesh distribution and smooth transition with the surrounding meshes. And it can effectively repair the cloud holes of complex boundary.

virtual repair; points cloud hole; the backtracking bidirectional wavefront method; triangular mesh

1671-7449(2017)06-0512-07

2017-03-18

山西省重點研發(fā)計劃(指南)資助項目(201603D321128)； 山西省研究生教育創(chuàng)新資助項目(2016SY051)； 山西省應(yīng)用基礎(chǔ)研究項目(201701D121062)

張 琦(1992-)， 女， 碩士生， 主要從事計算機圖形圖像處理研究.

TP391.9

A

10.3969/j.issn.1671-7449.2017.06.008