多视角立体视觉(multiple view stereo,MVS)重建是利用多张同一个场景的不同视角图像恢复场景3维模型的方法,在多媒体社区、虚拟现实、机器人技术、医学图像等领域有着广泛而迫切的需求。由文献[1]可知3维重建还可以为许多重要应用提供核心底层技术支持,如:全自主移动机器人的即时定位与地图构建(concurrent mapping and localization,CML)[2],通过2维视频构建3维多媒体数据[3],通过重建 3D组织表面进行计算机辅助腹腔镜手术[4]。基于多视图的3维重建一直是该领域研究的基本问题。
根据文献[5]可知,MVS算法可以分为4类:基于特征点扩展的方法[6–7]、基于表面演化的方法[8]、基于体素的方法[9]和基于深度图融合的方法[10–12]。由文献[12]分析可得:基于特征点扩展的方法先在有纹理区重建出稀疏点,再将这些点扩展到弱纹理区,在进行扩展时需对所有图像进行匹配搜索,但其计算复杂度高;基于表面演化的方法通过迭代演化初始曲面模型最小化能量函数,使其逐步收敛到目标的真实3维模型,但是在执行时易收敛于局部最优;基于体素的3维重建方法,需先建立包围目标模型的立方盒,再用优化方法将其逼近真实模型,但这类方法一般仅适用于小场景重建;基于深度图融合的方法先单独计算每幅图像所对应的深度图,再根据每幅图像对应的摄像机参数及其他条件还原完整的3D模型,这类方法灵活、方便,更适合大多数场景的3维重建。
作为许多深度图融合型3D建系统的核心组成部分,稠密立体匹配受到了研究者的高度关注[7–8,10,13]。在稠密立体匹配中,优化方式可以分为局部优化、半全局优化和全局优化。其中:基于局部的代价优化算法具有计算量小的优点,但仅是局部最优的,以致匹配结果精确度很低而不可靠;基于全局的优化算法是使用复杂计算的优化(如置信度传播算法、遗传算法),虽然可以获得比较理想的匹配结果,但是这类方法耗时、耗内存。基于半全局优化[14]的立体匹配方法表现出诸多优良特性,基于动态规划策略避免了NP问题,且匹配效果良好,在实际场合下最常被使用;但是,其计算开销比基于局部匹配优化算法更大;另外,由于因特网数据集规模的急剧增大及可用廉价计算资源的快速增加,将该方法应用于大场景图像集重建时,计算任务会更加繁重,使用单核CPU处理比较耗时。为了更加方便地在通用图形处理单元(GPGPU)中进行多视图立体匹配,作者提出一种并行半全局优化深度图计算方法。先通过平面扫描[15]方法和半全局优化方法计算深度图,再设计滤波为深度图降噪,利用3维点在图中的一致性剔除深度异常值并进行场景重建。本文方法充分考虑算法中各步骤的依赖关系和GPU处理海量数据时的优化规则,并基于此设计了一种新的并行处理架构和数据布局方式以克服计算效率瓶颈。
1 算法描述作者提出的基于GPU的半全局优化深度图计算方法包含以下步骤:1)由平面扫描方法计算初始匹配代价;2)通过并行半全局优化进行代价聚合,生成深度图;3)对深度图进行滤波,去除噪声点;4)通过一致性优化检测和剔除图中异常值。其中:第2)步并行半全局优化架构设计、实现将在第2节详细阐述。
1.1 平面扫描plane sweep匹配 1.1.1 平面扫描plane sweep框架文献[15]所提出的平面扫描(plane sweep)模型能够抽象为如图1所示的变换过程。给定参考图和相关图两幅图像,分别记为P1和P2,对应的摄像机记为C1和C2;扫描平面是C1坐标系下深度为d且平行于P1的成像仪平面;可定义单应矩阵
![]() |
图1 扫描平面变换过程 Fig. 1 Process of sweep plane transforming |
1.1.2 多视图匹配中的代价计算
将待计算深度的图像作为参考图像记为
${G_{{\rm{gvs}}}} = \{ {P_{{\rm{gv}}{{\rm{s}}_i}}}|i = 1,2, \cdots ,M\} $ | (1) |
将图像
$\mathit{\boldsymbol{H}} = {\mathit{\boldsymbol{K}}_{{\rm{gvs}}}}[{\mathit{\boldsymbol{R}}_{{\rm{gvs}}}}\mathit{\boldsymbol{R}}_{{\rm{ref}}}^{ - 1} - \frac{{({\mathit{\boldsymbol{C}}_{{\rm{gvs}}}} - {\mathit{\boldsymbol{R}}_{{\rm{gvs}}}}{\mathit{\boldsymbol{C}}_{{\rm{ref}}}}){\mathit{\boldsymbol{n}}^{\rm{T}}}}}{d}]\mathit{\boldsymbol{K}}_{{\rm{ref}}}^{ - 1}$ | (2) |
将图
${\mathit{\boldsymbol{x}}_{i,d}} = {\mathit{\boldsymbol{H}}_{i,d}}{\mathit{\boldsymbol{x}}_i}$ | (3) |
将
${\rho _i} = 1 - NCC({\mathit{\boldsymbol{x}}_{i,d}})$ | (4) |
当两像素相似度越大,对应
$Cost = \left(\sum\limits_{i = 1}^M {{\rho _i}} \right)/M$ | (5) |
由于单幅图像的匹配代价在估算过程中是在每个像素上独立进行的,没有考虑到像素间的制约关系,且邻居图和参考图之间可能会有尺度、光照上的差异,这些不确定因素都会引起初始相邻像素间的匹配代价相差过大,进而导致所得到的深度图不具有连续性。因此,使用优化方法处理第1.1节所得匹配代价异常值是十分自然的。
本算法构造类似文献[14]中的半全局匹配方法进行后续阶段的优化,立体匹配的本质为寻找使得匹配代价最小处对应的深度值。先在待匹配像素pos处沿多个(通常为8或4个)方向作动态规划。由于是在不同的方向上进行搜索,与传统的动态规划方法相比,顺序性约束一般难以满足。故必须使用扫描线最优算法,并将聚合方向设计为如图2所示。
![]() |
图2 聚合方向 Fig. 2 All paths’ directions for aggregation |
由文献[14]知,对于像素pos,在
$\begin{aligned}{L_r}(pos,d) = & Cost(pos,d) + \min ({L_r}(pos - r,d),\\& {L_r}(pos - r,d + 1) + {\delta _1},{L_r}(pos - r,d - 1) + {\delta _1},\\& \mathop {\min }\limits_i {L_r}(pos - r,i) + {\delta _2}) - \mathop {\min }\limits_k {L_r}(pos - r,k)\end{aligned}$ | (6) |
式中:
$S(pos,d) = \sum\limits_r {{L_r}(pos,d)} $ | (7) |
在计算
$depth(pos) = \arg \mathop {\min }\limits_d S(pos,d)$ | (8) |
半全局优化算法在能量优化过程中没有考虑噪声和遮挡问题,同时整个匹配优化过程难免会产生极值,这些因素都会严重影响深度图的质量。针对上述问题,先通过中值滤波剔除突变点;再将3维空间点在各深度图像上的一致性作为异常值检测的约束条件并作优化,核心步骤为:
1)根据所得深度图构造3D点对应关系,将所有像素反投影到全局坐标系下,其中,
2)每个3D点
${e_p} = {d_p}/{f_i}$ | (9) |
式中,
3)在3D空间中搜索与点
4)利用上一步建立起来的对应关系,为每张深度图像构建基于MAP-MRF优化模型,以过滤掉不可靠深度点。目标函数定义为:
${E^{{\rm{out}}}} = E_{\rm{d}}^{{\rm{out}}} + \gamma E_{\rm{s}}^{{\rm{out}}}$ | (10) |
式中:
$E_{\rm{d}}^{{\rm{out}}} = \sum\limits_p {[\exp ( - \beta {n_p})} - \alpha ]l(p)$ | (11) |
$E_{\rm{s}}^{{\rm{out}}} = \sum\limits_{q \in N(p)} {B(l(p),l(q))} $ | (12) |
式(11)~(12)中:
$B(x,y) = \left\{ \begin{aligned}& 1,\; {x \ne y};\\& 0,\; {x = y}\end{aligned} \right.$ | (13) |
目标函数(10)倾向于选择具有更多匹配点的空间3维点,同时维持点选择的平滑性,因为相邻点更可能具有相同的噪声水平,而异常点通常出现在局部地区。目标函数的优化可以通过文献[17]中所使用的图割方法实现。
2 并行计算架构半全局优化虽然较全局优化计算代价小,但是在处理单张图像时该算法的复杂度为
根据文献[14]可知,通过对多个方向进行代价聚合可提高匹配精度。要对多个方向进行代价聚合,最直接的方法是先将初始匹配代价传输到显存,依次在各方向进行代价聚合并将结果顺序传回内存,最后对聚合结果处理并得到深度值。这种模式存在两种缺点:一是,会引起GPU显存和CPU内存数据交换频繁,导致时间消耗严重;二是,各方向的代价聚合任务相互独立,当仅处理一个方向时,开启的线程数量有限,无法发挥GPU处理海量数据的优势。
为了改进这些缺陷,本文提出一种更合理的全局并行策略。首先,将待处理的匹配代价传入并始终存储于GPU显存;然后,同时开启在视差空间中每个深度平面上由左到右、由上到下、由左上到右下、由右上到左下以及它们的反方向共8个方向上的工作流,整个semiglobal优化过程和后续处理都在显存上完成,以便交给GPU足够多的任务;最后,处理结束后将结果传回内存,这样仅需要两次显存和内存的数据交换就能完成优化任务,并且能够充分利用GPU上的多核流处理器。
2.2 局部块并行在完成全局并行设计后,对各计算任务进行仔细划分。具体实施细节如图3所示,分别为由上到下(沿y正方向)和由左上到右下(沿x和y正方向)两个方向聚合,其余方向优化过程类似。值得注意的是如图3(b)所示当沿对角线聚合时应该在水平和竖直两个方向上同时进行。由于篇幅所限,仅介绍沿y正方向的代价聚合,算法1为在该方向上的伪码。
算法1 在y正方向上的代价聚合过程
$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!{\rm{Input:}} Cost(H,W,D),{H},{W},{D},{\delta _1},{\delta _2}$ |
$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!{\rm{Output}}:L_r(H,W,D)$ |
1. Parallel for
2. Block parallel for
3. Initialize aggr, min,
4. for
5.
6.
7.
8.
9.
10.
include the synchronization
图3(a)中,根据垂直于y平面上的待优化任务划分线程,同一x的不同d坐标被划分进同一线程块,每一个线程处理平面上相应位置的数据,并沿着y方向循环H次完成整个优化过程,如算法1中语句4所示;在循环过程中,进行平面上同一位置代价优化时,线程寄存器中原匹配代价能够直接被复用,邻居匹配代价需在共享内存中进行数据交换,如算法1中语句5、6、9所示;在同一线程块中进行并行规约,利用所有线程共同参与寻找最小值,如算法1中语句10所示。该设计从3个方面提升了计算效率。首先,充分考虑了该算法执行流程和各计算资源的依赖关系,使各线程在处理数据时能尽可能保持并行;其次,注重GPU上数据重用,以避免带宽限制:尽可能直接利用处理器周边数据(如寄存器、共享内存)避开全局内存瓶颈;最后,在规约时能够根据并行算法机制充分利用线程,并使用模板参数和展开循环优化代码,在执行该部分时理论时间复杂度由单线程的
![]() |
图3 并行半全局优化模式 Fig. 3 Scheme of parallel semiglobal optimization |
算法2为各方向代价求和深度赋值的伪码。在求和时由于每个像素都是独立进行的,故应按照像素坐标划分线程块,如语句1、2所示。为了减少全局内存的访问次数,使用语句4、5方式执行:首先,每个线程将和该像素相关的各方向上的代价求和;然后,同一个线程块内各线程合作寻找最小值。
算法2 代价求和与深度赋值
${\rm{Input:}}L_r = \{ L_r(H,W,D,r)|r = 1,2, \cdots ,8\} ,H,W,D$ |
$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!{\rm{Output:}}{depth}(H,W)$ |
1. Parallel for x=0 to W–1 do
2. Parallel for y=0 to H–1 do
3. Block parallel for
4.
5.
6. if
7.
所提出的方法已经使用C++和CUDA实现,并且在配有Intel Xeon E3 1230 CPU,NVidia GTX770 GPU和32G RAM的PC上测试。实验中各参数设置分别为:半全局优化聚合方向数量8,扫描平面数量128, 两个惩罚系数
为了展现提出方法在计算性能上的优势,将本文架构实现方法(其中,semiglobal优化和滤波两部分在GPU上实现,其他部分都在CPU单线程实现)分别与相同条件下本文算法在多核CPU中使用2线程和8线程实现方法进行定量比较。在4组数据集上,比较本文架构与CPU中2线程、CPU中8线程实现架构的运行时间,结果分别如表1、2所示。其中,所统计的时间是整个重建过程的时间耗费。由表1、2可知,用本文架构进行重建时效率更高。与CPU中2线程架构相比,本文架构实现了13.35倍到22.41倍的加速;与CPU中8线程架构相比,本文架构实现了5.71倍到9.13倍的加速。
表1 本文架构和CPU2线程架构的运行时间对比 Tab. 1 Comparison of runtime by architecture proposed in this paper and the CPU-2-thread implementation |
![]() |
表2 本文架构和CPU8线程架构的运行时间对比 Tab. 2 Comparison of runtime by architecture proposed in this paper and the CPU-8-thread implementation |
![]() |
将本文方法在相同条件下与文献[12]中的深度图计算方法进行比较,结果如图4所示。由图4可知,本文方法在各数据集上计算速度分别为文献[12]中的深度图计算方法5.41倍、6.52倍、5.40倍以及8.23倍,所提方法的加速效果都在5倍以上,在Lion123数据集上测试结果甚至可以超过8倍,这说明该方法在重建速度上表现出较大优势。
![]() |
图4 本文方法和文献[12]的重建时间对比 Fig. 4 Comparison of reconstruction time costs by the method proposed in reference[12] and this paper |
3.2 有效性 3.2.1 滤波效果
图5为滤波前后效果图。由图5可以看出,在GPU上实现的滤波起平滑作用,可过滤掉图像上的噪声点。
![]() |
图5 深度图滤波前后效果对比 Fig. 5 Comparison of results before and after filtering on depth map |
3.2.2 重构效果比较
使用本文方法进行深度图计算和优化,将所得深度图通过摄像机参数(该参数可通过Seitz等提供的Bundler[18]对原始图像标定获得)反投影回空间形成3维点云作为重建结果。图6为两组数据集Lion-p30和Lion123通过本文方法得到的重建结果。由图6可以看出重建出的点云能够较好地维持场景的主体结构。
![]() |
图6 输入图像与对应的输出点云 Fig. 6 Sample input images and corresponding reconstructed result in the experiments |
为了验证本文方法的有效性,将本文方法与Bailer等在文献[13]中提出的方法、Furukawa等研制的PMVS方法[7]进行比较。图7为该方法和文献[13]方法在ET、Kermit和Stone 3组数据集上的重建效果对比。
由图7可以看出,本文方法能够过滤掉大部分异常值,并且能得到更加密集和优质的3维点云。值得注意的是本文方法在Stone数据集重建出模型的汉字能够清晰识别,而使用文献[13]中方法得不到这种效果,说明所提方法在细节纹理丰富的区域,能够重建保留更多的信息。
![]() |
图7 文献[13]和本文方法的重建结果比较 Fig. 7 Comparision of reconstruction results by the method proposed in reference [13] and this paper |
使用本文方法与文献[13]中方法在Herziesu数据集上进行重建结果的比较,如图8所示。由图8可以看到,本文所提出的方法具有更高的重构一致性,例如,在整体效果上,本文方法可以成功重建出连续的模型(模型中较少出现空洞),而文献[13]中方法重建效果较差。本文方法不仅实现了较高的重构一致性,而且还能保持更好的局部细节信息,从局部放大细节可知:1)本文方法能够将该数据集中比较完整的浮雕重建出来,而文献[13]方法则不能;2)本文方法重建出的浮雕手臂包含更多的几何细节,而文献[13]无法捕获这些细节。
![]() |
图8 文献[13]和本文方法在Herziesu重建结果对比 Fig. 8 Comparision of reconstruction results on Herziesu by the method proposed in refenrence [13] and this paper |
3.2.3 精确度和完整度比较
通过计算各3D点云和真实值(测试数据集对应的基准模型,ground truth)之间的距离偏差,评估算法的精确度和完整度指标。本文采用基于文献[5]的开源点云比较软件对本文方法、文献[13]方法、PMVS方法[7]这3种方法在Herziesu数据集上进行精确度和完整度比较,结果如图9所示。由图9可以看出,在设定相同比率条件下,本文重建结果与真实值之间的平均匹配误差距离总是最小的,这说明使用本文所提方法可以有效提高重建的精度和完整性。
![]() |
图9 3种方法重建精度和完整度对比 Fig. 9 Comparision of accuracy and integrity on reconstruction results by the three methods |
4 结 论
针对大场景3维重建中深度图计算存在的瓶颈,提出了一种基于GPU的快速深度图计算方法。该方法充分考虑深度图计算时各步骤的依赖关系和GPU处理海量数据时的优化规则,并以此为基础设计出一种新的并行处理架构和数据布局方式;同时,设计滤波及优化技术对异常值进行剔除。从实验结果可知,该方法可以高效完成深度图计算且能够得到更好的重建效果。另外,由于半全局优化算法特点和GPU显存容量限制,因此,在下一步研究工作中,将尝试改进该优化算法和并行架构设计以适应对高分辨率图像做深度计算,并进一步提高重建方法的鲁棒性。
[1] |
Wu P, Liu Y, Ye M. Fast and adaptive 3D reconstruction with extensively high completeness[J]. IEEE Transactions on Multimedia, 2017, 19(2): 266-278. DOI:10.1109/TMM.2016.2612761 |
[2] |
Davison A J, Reid I D, Molton N D. MonoSLAM:Real-time single camera SLAM[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2007, 29(6): 1052-1067. DOI:10.1109/TPAMI.2007.1049 |
[3] |
Park M, Luo J, Gallagher A C. Learning to produce 3D media from a captured 2D video[J]. IEEE Transactions on Multimedia, 2013, 15(7): 1569-1578. DOI:10.1109/TMM.2013.2264926 |
[4] |
Maier-Hein L, Mountney P, Bartoli A. Optical techniques for 3D surface reconstruction in computer-assisted laparoscopic surgery[J]. Medical Image Analysis, 2013, 17(8): 974-996. DOI:10.1016/j.media.2013.04.003 |
[5] |
Seitz S M,Curless B,Diebel J,et al.A comparison and evaluation of multi-view stereo reconstruction algorithms[C]//Proceedings of the 2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition.New York:IEEE,2006,1:519–528.
|
[6] |
Goesele M,Snavely N,Curless B,et al.Multi-view stereo for community photo collections[C]//Proceedings of the IEEE 11th International Conference on Computer Vision.Rio de Janeiro:IEEE,2007:1–8.
|
[7] |
Furukawa Y, Ponce J. Accurate,dense,and robust multiview stereopsis[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2010, 32(8): 1362-1376. DOI:10.1109/TPAMI.2009.161 |
[8] |
Cremers D, Kolev K. Multiview stereo and silhouette consistency via convex functionals over convex domains[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2011, 33(6): 1161-1174. DOI:10.1109/TPAMI.2010.174 |
[9] |
Tabb A.Shape from silhouette probability maps:Reconstruction of thin objects in the presence of silhouette extraction and calibration error[C]//Proceedings of the 2013 IEEE Conference on Computer Vision and Pattern Recognition.Portland:IEEE,2013:161–168.
|
[10] |
Strecha C,von Hansen W,van Gool L,et al.On benchmarking camera calibration and multi-view stereo for high resolution imagery[C]//Proceedings of the 2008 IEEE Conference on Computer Vision and Pattern Recognition.Anchorage:IEEE,2008:1–8.
|
[11] |
Shen Shuhan. Accurate multiple view 3D reconstruction using patch-based stereo for large-scale scenes[J]. IEEE Transactions on Image Processing, 2013, 22(5): 1901-1914. DOI:10.1109/TIP.2013.2237921 |
[12] |
Liu Yiguang, Yi Shoulin, Wu Pengfei. A novel 3D reconstruction algorithm for large-scale scenes[J]. Journal of Sichuan University(Engineering Science Edition), 2015, 47(6): 91-96. [刘怡光, 易守林, 吴鹏飞. 一种新的大场景三维重建算法[J]. 四川大学学报(工程科学版), 2015, 47(6): 91-96.] |
[13] |
Bailer C,Finckh M,Lensch H P A.Scale robust multi view stereo[C]//Proceedings of the 2012 European Conference on Computer Vision.Berlin:Springer-Verlag,2012:398–411.
|
[14] |
Hirschmuller H. Stereo processing by semiglobal matching and mutual information[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2008, 30(2): 328-341. DOI:10.1109/TPAMI.2007.1166 |
[15] |
Collins R T.A space-sweep approach to true multi-image matching[C]//Proceedings of the 1996 IEEE Computer Society Conference on Computer Vision and Pattern Recognition.San Francisco:IEEE,1996:358–363.
|
[16] |
Hartley R,Zisserman A.Multiple view geometry in computer vision[M].Cambridge:Cambridge University Press,2003.
|
[17] |
Wu P, Liu Y, Ye M. Geometry guided multi-scale depth map fusion via graph optimization[J]. IEEE Transactions on Image Processing, 2017, 26(3): 1315-1329. DOI:10.1109/TIP.2017.2651383 |
[18] |
Seitz S M, Szeliski R, Snavely N. Photo tourism:Exploring photo collections in 3D[J]. ACM Transactions on Graphics, 2006, 25(3): 835-846. DOI:10.1145/1141911 |