基于回溯双向波前法虚拟修补点云孔洞

张 琦， 蔺素珍， 白佳璐， 王栋娟

(中北大学 大数据学院， 山西 太原 030051)

基于回溯双向波前法虚拟修补点云孔洞

张 琦， 蔺素珍， 白佳璐， 王栋娟

(中北大学 大数据学院， 山西 太原 030051)

针对虚拟修补边界点凹凸不平且分布不均的点云孔洞效果欠佳的问题， 提出了回溯双向波前法虚拟修补包含复杂边界的孔洞. 以具有复杂边界孔洞的青铜器模型为例： 首先， 三角网格化青铜器点云数据并依据网格化结果提取出孔洞边界， 通过比较边界点集的曲率波动幅度以去除伪孔洞边界； 其次， 以孔洞边界点集的回溯结果为初始点集， 并结合向量叉积对初始点集进行凹凸性分类， 再对凹、 凸点分别采用正、 反向波前法逐圈新增点集直至补全； 最后， 利用最小二乘法拟合曲面平滑新增点集， 获得最终修补结果. 实验结果表明该方法与划分子洞波前法、 曲线流收缩法相比， 结构相似性的平均值分别提高了81.14%和93.8%， 且曲率差异性更低， 修补网格的顶点密度与原始网格更相近且过渡自然， 能有效修补复杂边界孔洞.

虚拟修补； 点云孔洞； 回溯双向波前法； 三角网格

0 引 言

利用计算机虚拟复原技术辅助点云三维模型修补， 是提高逆向工程中三维重建效果和效率的重要手段. 文物这类特殊的三维模型因自然腐蚀、 坍塌等原因所致， 其出土时往往存在大小不一、 形状各异的复杂边界孔洞. 对这些孔洞进行虚拟精准修补， 既是后续纹理映射和三维效果重建的基础[1]， 也可为人工修补提供详细的操作依据和评价标准， 因此文物点云孔洞修补已成为其修补的关键技术之一.

目前， 关于文物点云孔洞修补的研究较少. 与之相关的主要有隧道[2]、 工业零件[3]以及动物模型[4]等的虚拟修补， 这些研究均证明将点云三角形网格化后再进行修补是行之有效的手段. 而基于三角形网格进行孔洞修补包括提取孔洞边界和新增点集. 依据边界边仅存在于一个三角面片的特性提取孔洞边界被认为是最为简单有效的方法[5,6]， 但鉴于其仅适用于提取全封闭模型的孔洞边界， 新近有学者根据三维点云二维投影后的8连通域对半封闭对象提取边界， 虽成功提取到孔洞边界， 却对存在投影重叠的情况无能为力[7]. 至于新增点集的方法， 依据待修补点云孔洞边界的特点大致可分为两类： 一是针对边界点全凹的孔洞， 多采用波前法[8]和基曲面网格添加法[9]等， 这些方法不仅能完整修补孔洞还能实现光滑连续； 二是对边界点凹凸共存且分布均匀的孔洞， 采用划分子洞波前法[10]和曲线流收缩法[11]等获得较好的修补效果. 与这些修补对象相比， 文物不仅大多呈半封闭形状而且在不同方向上投影大多存在点云重叠现象， 另外， 其孔洞多因腐蚀所致， 孔洞边界点往往凹凸不平且分布不均匀， 用以上方法进行修补难以获得理想效果.

为此， 本文以包含孔洞的青铜器点云为对象， 在提取边界后， 根据曲率波动幅度去除口部边缘等伪孔洞边界， 通过回溯到边界点集外圈建立新增点集以解决因边界点分布不均匀所致的修补误差问题， 然后对凹点采用正向波前法、 对凸点采用反向波前法逐圈新增点集， 在补全孔洞之后用原始的孔洞边界替换新增的第一圈点集， 再利用最小二乘法拟合曲面对剩余的新增点集进行微调， 即得到最终修补结果.

1 复杂边界点云孔洞修补总体思路

针对边界点凹凸不平且分布不均匀的点云孔洞， 本文提出的回溯双向波前法虚拟修补孔洞思路如图 1 所示， 其主要步骤为：

1) 对采集到的青铜器点云三角形网格化.

2) 根据边界边仅存在于一个三角形面片中的特点提取出孔洞边界.

3) 对孔洞边界和伪孔洞边界的曲率波动幅度进行比较， 去除伪孔洞边界.

4) 回溯到孔洞边界的外圈点作为新增点集的初始点集， 对凹、 凸点分别用正、 反向波前法建立新增点集， 反复进行直至孔洞补全.

5) 用孔洞原边界替换新增点集第一圈.

6) 用最小二乘法拟合曲面微调剩余新增点集， 得到最终修补结果.

图 1 孔洞虚拟修补思路图Fig.1 The schematic diagram of holes virtual repaired

2 点云孔洞修补的回溯双向波前法

2.1 孔洞边界自动提取

Delaunay三角剖分的三角网格具有空圆特性和最小内角最大化特性， 广泛应用于平面和空间的三角剖分[12]. 本文采用Delaunay三角形网格化法建立网格， 根据其边界边仅存在于一个三角形中的特性[5]提取孔洞边界. 由于青铜器多属半封闭模型， 所以提取出的边界还包含伪孔洞边界， 即原模型中的口部边缘等非闭合处. 鉴于伪孔洞边界多为器皿自身边缘， 其曲率连续性较好， 而腐蚀所致孔洞边界曲率连续性往往较差， 所以利用曲率分布特性可区分孔洞边界和伪孔洞边界点集. 去除伪孔洞边界的具体过程为：

2.1.1 计算各点曲率并建立其分布图.

图 2 曲率分布折线图Fig.2 The polyline diagram of the curvature distribution

不同孔洞边界的曲率分布如图 2 所示. 其中， 虚线与实线分别是伪孔洞边界和实际孔洞边界的曲率分布， 可以看出前者的波动幅度远小于后者.

2.1.2 计算曲率波动幅度D.

2.1.3 去除伪孔洞边界.

去除曲率波动幅度最小的伪孔洞边界， 即得到孔洞边界点集L=(p1,p2,…,pi,…,pn). 需要说明的是此方法仅能去除半封闭模型中边界曲率连续性较好的伪孔洞边界.

2.2 回溯双向波前法建立新增点集

2.2.1 回溯孔洞边界

青铜器点云孔洞边界点分布不均匀， 直接从此处开始建立新增点集会造成点集过密， 不利于保持与原网格的相似性， 进而会影响后续的纹理映射和三维效果重建. 本文利用网格的拓扑结构从孔洞边界进行回溯， 以其外圈点作为新增点集的初始点集.

2.2.2 基于双向波前法建立新增点集

传统的波前法虽是一种简单有效的修补方法， 但因青铜器孔洞边界点凹凸共存， 若对凸区用该方法建立新增点集， 则会增加到原始点云区域. 所以， 本文分别处理孔洞边界的凹凸区域， 即凹区采用传统的正向波前法、 凸区采用反向波前法修补. 正向波前法如文献[8]所述， 包括建立和更新波前点集两部分. 反向波前法的具体实现方法如下：

1) 以已提取的边界点集为波前点集， 计算波前点集每一个顶点pi相邻的边界边〈ei,ei+1〉间夹角θi.

2) 以最小夹角(minθi)对应点为初始点， 根据式(3)判断其凹凸性， 其中·为点积. 若d<0为凸点， 做反向处理； 反之则为凹点， 正向处理.

3) 基于pi对应的夹角θ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 利用最小二乘法拟合曲面微调新增点集

为保证所修补的孔洞能光滑连续， 在补全后还需要利用孔洞边界外圈的局部区域基于最小二乘法拟合曲面微调新增点集[13]. 首先， 用孔洞边界替换掉新增的第一圈点， 然后对剩余的每一圈都利用其外圈点集来近似拟合一个光滑的完整曲面， 实际上就是调整每一个点的z坐标， 计算方法为

式中：x和y为横纵坐标， 系数a1…a6通过式(8)计算获得.

式中：n为拟合曲面所选点的个数；ω(xi)是权函数， 取值为1；p(x)是基数，pT(x)=[1,x,y,x2,y2,xy].

3 实验结果及评价

3.1 定性评价

图 3(a) 为从山西省博物馆采集到的包含孔洞的青铜器数据(No.1)， 该数据有247 521个顶点， 包含两个孔洞， 图 3(b)～图 3(e)为采用本文方法修补前后的整体效果图和最终的模型三维图. 总体看， 本文方法不但恢复了缺失的孔洞部分， 而且修补区域能与原始区域光滑连接.

图 3 No.1青铜器孔洞修补效果图Fig.3 The result diagram of No.1 bronze holes repaired

文献[10]曾采用划分子孔洞波前法并通过最小二乘法拟合曲面微调获得比传统波前法更好的孔洞修补效果,文献[11]则利用曲线流收缩来修补孔洞. 为便于比较， 现将利用文献[10]、 文献[11]和本文方法修补孔洞1和孔洞2的结果分别放大， 如图 4 所示， 其中图4(a)～图4(c)为孔洞1的修补结果网格图， 图4(d) 和图4(e)为本文方法修补前后模型图； 图4(f)～图4(h)为孔洞2的修补结果网格图， 图4(i)和图4(j)为本文方法修补前后模型图. 图4(a), 图4(f), 图4(b), 图4(g)和图4(c), 图4(h)分别是用文献[10]、 文献[11]和本文方法的修补结果图， 从中可以看出本文方法修补网格更相似于原始网格， 说明其光滑连续性更好.

图 4 No.1青铜器孔洞修补结果细节图Fig.4 The results detail diagram of No.1 bronze holes repaired

图 5 是另外两组青铜器点云孔洞及其修补结果， 其中图5(a)， 图5(g)为原始数据的三角形网格图， 图5(b)， 图5(h)是文献[10]的修补结果， 图5(c)， 图5(i)是文献[11]的修补结果， 图5(d)， 图5(j)是本文方法的修补结果， 图5(e)， 图5(k)是修补前模型图， 图5(f)， 图5(l)是本文方法修补的模型结果. 尽管两组数据中孔洞所在面的曲率起伏不同， 但采用本文方法修补结果与原网格的相似性都优于文献[10]和文献[11].

图 5 No.2和No.3青铜器孔洞修补结果图Fig.5 The result diagram of No.2 and No.3 bronze holes repaired

3.2 定量评价

修补区域与原始区域网格越相似越有利于后续的纹理映射等三维重建； 修补区域与原始区域曲率连续性越好表明曲面越光滑， 综合两者可衡量点云孔洞修补质量的好坏.

3.2.1 结构相似性(structural similarity, SSIM)

结构相似性是常用的衡量三维模型空间相似性的指标， 其值越大表示结构越相似. 表 1 为利用文献[10]、 文献[11]方法和本文方法对修补区域和原始区域的结构相似性统计结果， 从表 1 中可以看出， 以上3组孔洞修补结果的SSIM指标， 本文方法都要高于文献[10]和文献[11]， 平均值分别高出81.14% 和93.8%， 说明本文方法具有更准确的修补结果.

表 1 3组孔洞数据的SSIM统计结果

3.2.2 曲率连续

曲率连续可以测量多个曲面在相交处的平滑程度. 曲率分布图可以直观体现曲率的连续性， 本文借助通用的Imageware中三角形网格曲率分布图对曲面的连续性进行评价. 图 6 为Imageware中3组数据的三角形网格曲率分布图， 图中各深浅色分别表示凹、 凸和较平坦区域， 颜色越接近说明连续性越好. 其中图6(a)， 图 6(d)， 图6(g)和图6(j)是用文献[10]方法的修补结果曲率分布图， 图6(b)，图6(e)， 图6(h) 和图6(k)是用文献[11]方法的修补结果曲率分布图， 图6(c)， 图6(f)， 图6(i)和图6(l)是本文方法的修补结果曲率分布， 画圈区域为原孔洞所在区域. 直观上看， 孔洞周围区域以凹凸为主， 文献[10] 修补结果较平坦， 文献[11]修补结果凸区较多， 而本文修补结果与周围区域颜色分布更为相似.

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

表 2 3种区域所占百分比%

为便于定量描述， 采用所修补孔洞与其周围10圈(实物长度约0.5 cm)区域中凸、 凹、 平3种区域的分布百分比进行比较， 差异范围越小与原始区域凹凸性越相似， 如表 2 所示： 在3组对象中， 本文方法与修补前的百分比差异范围在0～4.4之间， 而文献[10]的差异在0.8～34.1之间， 文献[11]的差异在2.6～15.3之间， 表明本文方法修补后能更好地保持孔洞与周围区域的凹凸一致性. 综上， 在用本文方法修补后的曲率分布图中孔洞均与孔洞周围区域的曲率图呈光滑连接， 而用文献[10]、 文献[11]修补后的孔洞区域边界处可以保持光滑连接但越靠近孔洞中心位置曲率分布越不连续. 由此可见， 本文方法修补结果曲率连续性更好， 保证了孔洞与周围区域的光滑连接.

4 结 论

本文的回溯双向波前法基于半封闭模型， 在三角形网格基础上有效地提取出孔洞边界并实现了复杂边界孔洞的修补， 具有较高的结构相似性和曲率连续性， 其主要特点体现在： ① 利用回溯双向波前法可有效修补青铜器等半封闭模型中边界点凹凸不平且分布不均匀的点云孔洞； ② 依据孔洞边界点曲率波动幅度能有效去除伪孔洞边界； ③ 将孔洞边界回溯到其外圈， 有利于减小因孔洞边界点分布不均匀给新增点集带来的误差； ④ 使用双向波前法对兼具凹凸点的边界建立新增点集简单实用； 用最小二乘法微调修补部分， 可使得修补区域与周围连接更光滑.

值得指出的是本文方法的局限性是修补尖锐部分及非闭合状的孔洞效果欠佳， 下一步将以此作为重点进行研究.

[1] 姜翰青, 王博胜, 章国锋， 等. 面向复杂三维场景的高质量纹理映射[J]. 计算机学报, 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] 钱伯至, 蓝秋萍. 隧道三维点云孔洞修复方法[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] 王钦瑞. 三角网格模型特征线提取与孔洞修补方法研究[D]. 大连： 大连理工大学, 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] 刘震, 王艳宾, 白丽丽, 等. 曲面细节特征保持的三维模型孔洞修复算法[J]. 计算机辅助设计与图形学学报, 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] 魏潇然, 耿国华, 张雨禾. 几何信息预测的三角网格模型拓扑压缩[J]. 西安电子科技大学学报, 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

山西省重点研发计划(指南)资助项目(201603D321128)； 山西省研究生教育创新资助项目(2016SY051)； 山西省应用基础研究项目(201701D121062)

张 琦(1992-)， 女， 硕士生， 主要从事计算机图形图像处理研究.

TP391.9

A

10.3969/j.issn.1671-7449.2017.06.008