工程科学与技术   2017, Vol. 49 Issue (6): 114-121
一种基于GPU的快速半全局优化深度图计算方法
刘怡光, 赵洪田, 吴鹏飞, 徐振宇, 都双丽, 李杰     
四川大学 计算机学院,四川 成都 610065
基金项目: 国家自然科学基金资助项目(61571313);四川省科技创新苗子工程资助项目(2016017)
摘要: 由于图像集规模巨大、匹配信息丰富,快速精准多视图立体匹配受计算效率严重制约。针对该问题,提出一种基于GPU的快速半全局优化深度图计算方法。首先,在CPU上通过平面扫描方法计算单张图像初始匹配代价。然后,提出GPU半全局优化并行计算架构,对匹配代价进行聚合,其核心算法为:在全局进行各方向聚合任务流并行以提升众核处理器的利用率;在局部通过将各像素计算任务准确分配到各线程块内实现并行处理,且注重GPU上数据重用以避免带宽限制。再通过GPU滤波剔除突变点进行图像增强。最后,将3维空间点在各深度图像上的一致性作为异常值检测和优化的约束条件。在多组数据集上测试结果显示,该方法计算速度最高为多核CPU系统中开启2线程实现方法的22.41倍,为开启8线程实现方法的9.13倍,且与两者精度相当;与同类深度图计算方法比较结果表明,该方法在重建过程中加速效果均为其他算法的5倍及以上;通过使用开源点云比较软件在标准测试数据集上与其他算法比较,验证了该方法能有效提高重建结果的精度和完整度。
关键词: 3维重建    平面扫描    深度图    并行半全局优化    GPU    滤波    
A Fast Semiglobal Optimization Algorithm for Computing Depth Map with GPU
LIU Yiguang, ZHAO Hongtian, WU Pengfei, XU Zhenyu, DU Shuangli, LI Jie     
College of Computer Sci.,Sichuan Univ.,Chengdu 610065,China
Abstract: The accuracy and efficiency of multi-view stereo matching are seriously hindered by the scale of the image datasets and the wide variety of matching information.To address this issue,in this paper,a fast semiglobal optimizing method based on GPU was proposed for depth map calculation with parallel computing technology applied.First,the matching cost of single image was computed by the plane sweep method on CPU. Second,the matching costs were aggregated by using semiglobal optimizing method executed on GPU in parallel.The core module included two strategies:1) by aggregating the cost in each direction globally and parallelly,the utilization rate of the multi-cores of the GPU was increased;2) The computation on each single pixel was assigned to different threads and computed in parallel.To alleviate the constraint of limited bandwidth,the data reuse of GPU had also been taken into consideration.Then,the filtering on GPU was proposed to remove the outliers and enhance the image.Finally,the consistency of the 3D point projections on every image was used as a constraint for outlier detection and optimization.The proposed method had been tested on multi-group datasets.Experiments showed that the proposed method achieved as high as 22.41 times and 9.13 times speedup compared to two-threaded and eight-threaded implementation on a multi-core CPU system respectively,while also preserving similar quality.Compared with other algorithms of depth map calculation on tested datasets,the proposed method always achieved as fast as more than 5 times speedup in the reconstruction process.Compared with other algorithms on several datasets via open source pointcloud comparison software (CloudCompare),the proposed method could effectively improve the accuracy and integrity of the reconstruction results.
Key words: 3D reconstruction    plane sweep    depth map    parallel semiglobal optimization    GPU    filtering    

多视角立体视觉(multiple view stereo,MVS)重建是利用多张同一个场景的不同视角图像恢复场景3维模型的方法,在多媒体社区、虚拟现实、机器人技术、医学图像等领域有着广泛而迫切的需求。由文献[1]可知3维重建还可以为许多重要应用提供核心底层技术支持,如:全自主移动机器人的即时定位与地图构建(concurrent mapping and localization,CML)[2],通过2维视频构建3维多媒体数据[3],通过重建 3D组织表面进行计算机辅助腹腔镜手术[4]。基于多视图的3维重建一直是该领域研究的基本问题。

根据文献[5]可知,MVS算法可以分为4类:基于特征点扩展的方法[67]、基于表面演化的方法[8]、基于体素的方法[9]和基于深度图融合的方法[1012]。由文献[12]分析可得:基于特征点扩展的方法先在有纹理区重建出稀疏点,再将这些点扩展到弱纹理区,在进行扩展时需对所有图像进行匹配搜索,但其计算复杂度高;基于表面演化的方法通过迭代演化初始曲面模型最小化能量函数,使其逐步收敛到目标的真实3维模型,但是在执行时易收敛于局部最优;基于体素的3维重建方法,需先建立包围目标模型的立方盒,再用优化方法将其逼近真实模型,但这类方法一般仅适用于小场景重建;基于深度图融合的方法先单独计算每幅图像所对应的深度图,再根据每幅图像对应的摄像机参数及其他条件还原完整的3D模型,这类方法灵活、方便,更适合大多数场景的3维重建。

作为许多深度图融合型3D建系统的核心组成部分,稠密立体匹配受到了研究者的高度关注[78,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所示的变换过程。给定参考图和相关图两幅图像,分别记为P1P2,对应的摄像机记为C1C2;扫描平面是C1坐标系下深度为d且平行于P1的成像仪平面;可定义单应矩阵 ${\mathit{\boldsymbol{H}}_i}\left( {i = 1,2, \cdots ,n} \right) $ 将图像P2平面上的点集位置与成像仪平面Di的点集位置相关联,并将相应变换记为H变换(或H透视变换)。

图1 扫描平面变换过程 Fig. 1 Process of sweep plane transforming

1.1.2 多视图匹配中的代价计算

将待计算深度的图像作为参考图像记为 ${P_{{\rm{ref}}}}$ ,相关邻居图记为 ${P_{{\rm{gv}}{{\rm{s}}_i}}}\left( {i = 1,2, \cdots ,V} \right)$ 。根据文献[13]中提出的选图方法,在输入图像序列中选出与 ${P_{{\rm{ref}}}}$ 相关联的前M幅图像参与立体匹配,将所选图像集合记为 ${G_{{\rm{gvs}}}}$ ,表达式为:

${G_{{\rm{gvs}}}} = \{ {P_{{\rm{gv}}{{\rm{s}}_i}}}|i = 1,2, \cdots ,M\} $ (1)

将图像 ${P_{{\rm{ref}}}}$ ${G_{{\rm{gvs}}}}$ 中一个子集 ${P_{{\rm{gvs}}}}$ 对应的摄像机参数分别记为 $\{ {\mathit{\boldsymbol{K}}_{{\rm{ref}}}},{\mathit{\boldsymbol{R}}_{{\rm{ref}}}},{\mathit{\boldsymbol{C}}_{{\rm{ref}}}}\} $ $\{ {\mathit{\boldsymbol{K}}_{{\rm{gvs}}}},{\mathit{\boldsymbol{R}}_{{\rm{gvs}}}},{\mathit{\boldsymbol{C}}_{{\rm{gvs}}}}\} $ ,其中, ${\mathit{\boldsymbol{K}}_{{\rm{ref}}}}$ ${\mathit{\boldsymbol{K}}_{{\rm{gvs}}}}$ 分别为对应摄像机内部参数, ${\mathit{\boldsymbol{R}}_{{\rm{ref}}}}$ ${\mathit{\boldsymbol{R}}_{{\rm{gvs}}}}$ 分别为对应旋转矩阵, ${\mathit{\boldsymbol{C}}_{{\rm{ref}}}}$ ${\mathit{\boldsymbol{C}}_{{\rm{gvs}}}}$ 分别为对应摄像机中心。定义 $\textit{π}{(d)} =( {\mathit{\boldsymbol{n}}^{\rm{T}}},d)$ 为深度是 $d$ 的扫描平面。在本算法中将 ${\mathit{\boldsymbol{n}}^{\rm{T}}} = (0,0,1)$ 作为Z轴方向(以参考摄像机为原点)。根据文献[16]可知,这两幅图像关于扫描平面 $\textit{π} {(d)}$ 的单应矩阵为:

$\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)

将图 ${P_{{\rm{gv}}{{\rm{s}}_i}}}$ 通过式(2)计算关于 $\textit{π} {(d)}$ 的单应矩阵 ${\mathit{\boldsymbol{H}}_{i,d}}$ 。图 ${P_{{\rm{gv}}{{\rm{s}}_i}}}$ 上的点记为 ${\mathit{\boldsymbol{x}}_i}$ ,将该点经式(3)变换后记作 ${\mathit{\boldsymbol{x}}_{i,d}}$ ,变换过程如图1所示,在 ${P_{{\rm{ref}}}}$ 中相应位置的点记为 ${\mathit{\boldsymbol{x}}_{i,{\rm {ref}}}}$

${\mathit{\boldsymbol{x}}_{i,d}} = {\mathit{\boldsymbol{H}}_{i,d}}{\mathit{\boldsymbol{x}}_i}$ (3)

${G_{{\rm{gvs}}}}$ 中不同图像像素分别经相应H透视变换到各 $\textit{π} (d)$ 平面上,对于参考图像中每个像素 ${x}_{i,{\rm ref}}$ ,以 $ {x}_{i,{\rm ref}}$ 为中心构建局部窗口,并与 $ {{x}_{i,d}}$ 一同计算归一化互相关系数(normalized correlation coefficient,NCC),并通过式(4)计算匹配代价:

${\rho _i} = 1 - NCC({\mathit{\boldsymbol{x}}_{i,d}})$ (4)

当两像素相似度越大,对应 ${\rho _i}$ 越小。最后,由式(5)计算得到 $\textit{π} {(d)}$ 关于 ${P_{{\rm{ref}}}}$ 中各像素的初始匹配代价:

$Cost = \left(\sum\limits_{i = 1}^M {{\rho _i}} \right)/M$ (5)
1.2 半全局(semiglobal)优化

由于单幅图像的匹配代价在估算过程中是在每个像素上独立进行的,没有考虑到像素间的制约关系,且邻居图和参考图之间可能会有尺度、光照上的差异,这些不确定因素都会引起初始相邻像素间的匹配代价相差过大,进而导致所得到的深度图不具有连续性。因此,使用优化方法处理第1.1节所得匹配代价异常值是十分自然的。

本算法构造类似文献[14]中的半全局匹配方法进行后续阶段的优化,立体匹配的本质为寻找使得匹配代价最小处对应的深度值。先在待匹配像素pos处沿多个(通常为8或4个)方向作动态规划。由于是在不同的方向上进行搜索,与传统的动态规划方法相比,顺序性约束一般难以满足。故必须使用扫描线最优算法,并将聚合方向设计为如图2所示。

图2 聚合方向 Fig. 2 All paths’ directions for aggregation

由文献[14]知,对于像素pos,在 $r$ 方向上的匹配代价可定义为:

$\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)

式中: $Cost(pos,d)$ $\textit{π} {(d)}$ 在像素pos点由式(5)计算得到的初始匹配代价; ${\delta _1}$ ${\delta _2}$ 分别为惩罚系数,保持深度图的平滑性和防止边缘出现锯齿状;最后一项为防止累积匹配代价值过大而减去上一像素点的最小优化代价。式(6)所表示的传递优化过程可使邻居像素代价按照顺序依次取最小,并提升相邻像素深度值变化的一致性。将各个方向上的优化代价相加,作为该点的最终代价:

$S(pos,d) = \sum\limits_r {{L_r}(pos,d)} $ (7)

在计算 $S(pos,d)$ 后,为了减少优化复杂度,直接采用获胜者全拿(winner takes all,WTA)策略选择最小深度值:

$depth(pos) = \arg \mathop {\min }\limits_d S(pos,d)$ (8)
1.3 滤 波

半全局优化算法在能量优化过程中没有考虑噪声和遮挡问题,同时整个匹配优化过程难免会产生极值,这些因素都会严重影响深度图的质量。针对上述问题,先通过中值滤波剔除突变点;再将3维空间点在各深度图像上的一致性作为异常值检测的约束条件并作优化,核心步骤为:

1)根据所得深度图构造3D点对应关系,将所有像素反投影到全局坐标系下,其中, $i \in \varPhi $ 表示图像的索引, $\varPhi $ 为所有图像的集合。

2)每个3D点 $p$ 的影响范围 $e$ 可以定义为:

${e_p} = {d_p}/{f_i}$ (9)

式中, ${d_p}$ $p$ 的深度, ${f_i}$ 为包含点 $p$ 对应图像的焦距。

3)在3D空间中搜索与点 $p$ 距离小于 ${e_p}$ 的点,并将这些点设置为 $p$ 的对应点。值得注意的是, $p$ 的对应点和 $p$ 不可以来自于同一幅深度图。将点 $p$ 的对应点个数设置为 ${n_p}$ ,以建立图像间可视3维点的一致性。

4)利用上一步建立起来的对应关系,为每张深度图像构建基于MAP-MRF优化模型,以过滤掉不可靠深度点。目标函数定义为:

${E^{{\rm{out}}}} = E_{\rm{d}}^{{\rm{out}}} + \gamma E_{\rm{s}}^{{\rm{out}}}$ (10)

式中: $\gamma $ 为权重参数; $E_{\rm{d}}^{{\rm{out}}}$ $E_{\rm{s}}^{{\rm{out}}}$ 分别为数据项和平滑项,分别被定义为:

$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)中: $l( * )$ 是给定像素的标签,若要过滤掉该像素则置为0,否则置为1; $\alpha $ 为不选择某个像素的惩罚参数,当 $\alpha $ =0时倾向于过滤掉某个像素,但是并不想总是过滤掉这些点,因此将其设为 $\alpha > 0$ $\beta $ 是一个权重参数 $N(p)$ 表示 $p$ 的临近点集合,通常定义为 $p$ 投影在图上点的4邻域; $B( * )$ 为分段平滑函数,主要起到插值扩散和图像增强的作用,

$B(x,y) = \left\{ \begin{aligned}& 1,\; {x \ne y};\\& 0,\; {x = y}\end{aligned} \right.$ (13)

目标函数(10)倾向于选择具有更多匹配点的空间3维点,同时维持点选择的平滑性,因为相邻点更可能具有相同的噪声水平,而异常点通常出现在局部地区。目标函数的优化可以通过文献[17]中所使用的图割方法实现。

2 并行计算架构

半全局优化虽然较全局优化计算代价小,但是在处理单张图像时该算法的复杂度为 $O(W \times H \times D)$ ,其中, $W$ $H$ 为图像的宽度和高度, $D$ 为深度范围(一般设置为64或128)。当将半全局优化运用在多视图立体匹配的核心优化中,整个算法执行过程会更加耗时。文献[14]的算法设计比较复杂,通过更改该方法的执行流程可以达到降低计算复杂度和内存消耗的目的,但是,最终所得图像配准精度也会随之降低。为了保留原算法的优良特性同时提高执行效率,作者在GPU上设计一种分层并行架构以适应第1节基于半全局(semiglobal)优化模块的实现。

2.1 全局数据流并行

根据文献[14]可知,通过对多个方向进行代价聚合可提高匹配精度。要对多个方向进行代价聚合,最直接的方法是先将初始匹配代价传输到显存,依次在各方向进行代价聚合并将结果顺序传回内存,最后对聚合结果处理并得到深度值。这种模式存在两种缺点:一是,会引起GPU显存和CPU内存数据交换频繁,导致时间消耗严重;二是,各方向的代价聚合任务相互独立,当仅处理一个方向时,开启的线程数量有限,无法发挥GPU处理海量数据的优势。

为了改进这些缺陷,本文提出一种更合理的全局并行策略。首先,将待处理的匹配代价传入并始终存储于GPU显存;然后,同时开启在视差空间中每个深度平面上由左到右、由上到下、由左上到右下、由右上到左下以及它们的反方向共8个方向上的工作流,整个semiglobal优化过程和后续处理都在显存上完成,以便交给GPU足够多的任务;最后,处理结束后将结果传回内存,这样仅需要两次显存和内存的数据交换就能完成优化任务,并且能够充分利用GPU上的多核流处理器。

2.2 局部块并行

在完成全局并行设计后,对各计算任务进行仔细划分。具体实施细节如图3所示,分别为由上到下(沿y正方向)和由左上到右下(沿xy正方向)两个方向聚合,其余方向优化过程类似。值得注意的是如图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 $x = 0$ to $W - 1$ do

2. Block parallel for $Thread = 0$ to $D - 1$ do

3. Initialize aggr, min, $share(D)$ with max;

4. for $y = 0$ to $H - 1$ do

5. $front = {\rm{ }}share(Thread - 1)$ ;

6. $back = {\rm{ }}share(Thread + 1)$ ;

7. $\begin{array}{l}aggr = Cost(y,x,Thread) + {\rm{minimum}}(aggr,\\ \quad\quad\quad front + {\delta _1},back + {\delta _1},min + {\delta _2}) - min;\end{array}$

8. $L_r(y,x,Thread){\rm{ = }}aggr$ ;

9. $share(Thread) = aggr$ ;

10. $min = min\_merging\_reduce\left( {aggr} \right)$ ; //

include the synchronization

图3(a)中,根据垂直于y平面上的待优化任务划分线程,同一x的不同d坐标被划分进同一线程块,每一个线程处理平面上相应位置的数据,并沿着y方向循环H次完成整个优化过程,如算法1中语句4所示;在循环过程中,进行平面上同一位置代价优化时,线程寄存器中原匹配代价能够直接被复用,邻居匹配代价需在共享内存中进行数据交换,如算法1中语句5、6、9所示;在同一线程块中进行并行规约,利用所有线程共同参与寻找最小值,如算法1中语句10所示。该设计从3个方面提升了计算效率。首先,充分考虑了该算法执行流程和各计算资源的依赖关系,使各线程在处理数据时能尽可能保持并行;其次,注重GPU上数据重用,以避免带宽限制:尽可能直接利用处理器周边数据(如寄存器、共享内存)避开全局内存瓶颈;最后,在规约时能够根据并行算法机制充分利用线程,并使用模板参数和展开循环优化代码,在执行该部分时理论时间复杂度由单线程的 $O(N)$ 减少为众线程的 $O({\rm lb}\;N)$

图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 $Thread = 0$ to $D - 1$ do

4. $tmp = \displaystyle\sum\limits_{r = 1}^8 {L_r(y,x,Thread,r)} $ ;

5. $\min Idx = min\_merging\_reduce(tmp,$ $Thread)$ ; // include the synchronization

6. if $Thread{\rm{ }} = = {\rm{ }}0$ then

7. ${depth}(y,x) = \min Idx$ ;

3 实 验

所提出的方法已经使用C++和CUDA实现,并且在配有Intel Xeon E3 1230 CPU,NVidia GTX770 GPU和32G RAM的PC上测试。实验中各参数设置分别为:半全局优化聚合方向数量8,扫描平面数量128, 两个惩罚系数 ${\delta _1}$ =0.04, ${\delta _2}$ =0.5,滤波窗口大小3 $ \times $ 3,优化参数分别设置为 $\alpha = 1$ $\beta = 0.5$ $\gamma = 3$ 。并在6组相关数据集中进行测试,它们是Seitz等[18]提供的ET、Kermit和Strecha等提供的Herziesu[10],以及3组户外采集的数据集Lion-p30、Lion123、Stone。

3.1 速 度

为了展现提出方法在计算性能上的优势,将本文架构实现方法(其中,semiglobal优化和滤波两部分在GPU上实现,其他部分都在CPU单线程实现)分别与相同条件下本文算法在多核CPU中使用2线程和8线程实现方法进行定量比较。在4组数据集上,比较本文架构与CPU中2线程、CPU中8线程实现架构的运行时间,结果分别如表12所示。其中,所统计的时间是整个重建过程的时间耗费。由表12可知,用本文架构进行重建时效率更高。与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