CN111968218A

CN111968218A - 基于gpu集群的三维重建算法并行化方法

Info

Publication number: CN111968218A
Application number: CN202010704955.0A
Authority: CN
Inventors: 黄方; 彭书颖; 杨浩; 铁博; 陈胤杰
Original assignee: University of Electronic Science and Technology of China
Current assignee: University of Electronic Science and Technology of China
Priority date: 2020-07-21
Filing date: 2020-07-21
Publication date: 2020-11-20

Abstract

本发明公开了一种基于GPU集群的三维重建算法并行化方法，涉及计算机视觉技术领域。该方法以SfM算法为基础，研究无人机影像三维重建技术流程，同时针对无人机三维稠密重建处理耗时这一问题，采用GPU集群作为处理平台，具体地，基于SFM_MVS相关理论，掌握基于图片序列的实景三维重建相关流程，同时分别采用了MPI并行编程技术和GPU并行编程技术对三维重建流程部分环节进行优化加速研究工作，其中利用集群进行稀疏重建算法算子替换有效解决了无人机航拍影像数据量大、计算耗时的问题，稠密重建算法粗粒度数据并行与稠密匹配算法特征提取环节的细粒度并行及优化有效加速了三维重建后期稠密点云重建流程。

Description

基于GPU集群的三维重建算法并行化方法

技术领域

本发明属于GPU集群与并行计算领域，涉及一种基于GPU集群的三维重建算法并行化方法。

背景技术

近年来，无人机技术已经发展成为一种广泛使用的较为新颖的遥感数据获取平台。和其他遥感平台相比，无人机平台具有成本低廉、使用方便，结合地面站的航线规划功能上手快、外出实验操作简单等方面的优势。因无人机飞速发展产生的无人机三维重建技术，通过无人机搭载多角度相机拍摄地物，较为轻易地获取物体侧面纹理信息，从而能够获取更为全面的三维信息，在城市规划、建设、管理和应急响应中起着极其重要的作用。

目前无人机三维重建技术主要是在Windows环境下依托专业的自动化程度较高的无人机测绘和摄影测量软件，能够进行全自动联合空中三角测量、三维模型构建、贴纹理等操作，形成最终三维真实场景。相比于传统的建模方法，利用无人机进行三维实景重建，在数据获取、处理方面效率都有着显著的提升。但是，在使用这些专业的商业软件时，依旧会存在以下问题：

(1)相对于传统手工建模的方式，使用商业软件进行三维重建有着不少的效率提升，但进行大规模作业时，处理周期依然较长；

(2)商业软件的硬件环境要求也随着软件版本的更新换代也在不断提高，普通办公所用的台式机往往无法满足，在一定程度上也带来了诸多不便；

(3)某些商业软件虽然也能支持在局域网下进行并行化处理，但是其通信模块不能够很好承担各处理平台之间的通信问题，经常出现通信死锁、任务缺失问题；且商业软件闭源，无法进行进一步针对性的修改与开发。

计算机视觉和传统摄影测量技术的融合发展，带来了运动恢复结构(SfM)这一新型的地形测量技术，它可以让人们以较低的成本获得关于地貌高质量的、密集的三维点云数据。采用开源的计算机视觉算法进行实景三维重建拥有着源码易获取、移植性高等特点，但是随着人们对点云质量要求的提高，SfM得到的稀疏点云无法直接进行应用，需要进一步三维重建得到三维稠密点云，此时的数据处理耗时变成为一大问题。

采用高性能计算技术进行数据处理加速是一个高效的解决方案，其中的一个重要分支，CPU+GPU异构计算，结合了CPU和GPU两者的优势，被不少研究人员用于处理图形图像解决效率低下的首选方式。然而，单个文件/任务处理速度的降低依旧无法掩盖数据量的加大导致的处理时间漫长的问题，面对此类问题，相比于其他处理方式，采用集群处理方式使用多个节点同时进行处理则有着明显的优势。GPU集群系统和传统的集群有所不同，GPU集群中的每个计算节点在配备常规CPU计算核心的前提下，每个节点必备一个或者多个GPU处理器，这样能充分发掘CPU的逻辑运算能力和GPU的浮点计算能力，相较于传统计算平台能够发挥出难以匹敌的优势。

因此，本发明基于计算机视觉相关算法，以SfM算法为基础，研究无人机影像三维重建技术流程，同时针对无人机三维稠密重建处理耗时这一问题，采用GPU集群作为处理平台，具有重要的理论研究和实践应用价值。

发明内容

本发明公开了一种基于GPU集群的三维重建算法并行化方法，该方法包括以下步骤：

S1)给出基于图片特征的三维重建算法的具体流程；

S2)分析三维重建算法的具体流程中最耗时的部分；

S3)对三维重建流程中最耗时的部分进行并行优化；

S4)对三维重建流程中的其余部分进行优化；

其中，步骤S1)中基于图片特征的三维重建算法的具体流程包括：利用无人机获取到地物图像序列后，首先进行特征提取，在一系列不同的图像中识别同名关键点，正是这些关键点实现了不同图像间的匹配，以及场景几何的重构。一旦定位了每个图像中的关键点，就需要确定不同图像中关键点之间的对应关系，为了确保仅仅保留正确的对应关系，有必要采取进一步的步骤过滤掉剩余的任何错误的匹配，由于关键点仅限于具有几何一致性匹配的关键点，所以，每个图像对之间的链接可以被识别，并且组织为轨迹；然后运动恢复结构理论基于这些轨迹和关键点以及其他相关的图片信息，进行几何解算，可以粗略描述场景几何信息的稀疏点云，SfM产生的稀疏点云通常只是使用MVS生成更加密集点云数据的中间步骤；MVS的目标是从已知的相机参数的图像集合中提供完整的三维场景重建；选用开源的Bundler处理稀疏重建流程，选用PMVS算法软件处理稠密重建流程。

步骤S2)中分析三维重建算法的具体流程中最耗时的部分包括：在Linux环境下通过处理不同规模的数据集，记录各环节耗时情况，采用多次实验求平均的思想，分析三维重建算法的具体流程中的瓶颈环节，得到在基于图像序列的三维重建算法的具体流程中，稠密重建流程耗时占比最大，至少占全部流程的60％以上。

步骤S3)在深入分析算法原理与流程的基础上，分析并确定算法的可并行性及并行优化方法，实现在集群环境下基于MPI并行编程技术的CMVS稠密点云生成算法；为进一步加快三维重建流程，充分发挥GPU处理图形图像优势，本研究细化分析三维重建流程，根据流程热点分析结果结合GPU硬件特性，在GPU单节点环境下设计优化稠密重建面片优化环节，实现MPI+GPU的二级稠密重建并行流程；使用GPU并行技术处理稠密重建特征提取环节，同时采用合并核函数等方式，进一步优化对应的GPU并行算法，建立GPU集群环境下二级并行模型；

步骤S4)将稀疏重建环节的特征提取阶段替换为基于GPU平台的优化算法，形成基于GPU环境下的稀疏重建流程，实现GPU环境下稀疏重建加速处理。

本发明基于SFM_MVS相关理论，掌握了基于图片序列的实景三维重建相关流程。同时分别采用了MPI并行编程技术和GPU并行编程技术对三维重建流程部分环节进行优化加速研究工作，其中利用集群进行稀疏重建算法算子替换有效解决了无人机航拍影像数据量大、计算耗时的问题，稠密重建算法粗粒度数据并行与稠密匹配算法特征提取环节的细粒度并行及优化有效加速了三维重建后期稠密点云重建流程。

附图说明

图1为本发明的技术路线图

图2为本发明的基于图像序列的三维重建算法流程示意图

图3为本发明的CMVS并行优化设计图

图4为本发明的调用串行PMVS过程程序流程图

图5为本发明的面片及光度差异函数示意图

图6为本发明的MPI+GPU二级并行优化稠密重建处理示意图

图7为本发明的Harris特征提取原理图

图8为本发明的梯度变化图

图9位本发明的Harris特征提取环节并行化思路示意图

图10为本发明的梯度计算并行算法的处理过程示意图

图11为本发明的矩阵计算并行算法的处理过程示意图

图12为本发明的Harris特征提取环节优化后并行示意图

图13为本发明的SIFT-GPU加速示意图

具体实施方式

下面结合附图对本发明的具体实施方式以及工作原理作进一步详细说明。

本发明技术路线如图1所示，首先掌握三维重建方法整体流程，了解各流程耗时情况，明确流程中耗时主体部分，即为三维稠密重建部分的前提下；之后依据原理流程分析可并行性，采用MPI并行编程技术，在多节点集群环境下实现基于CMVS算法的粗粒度并行优化策略；在明确稠密重建内部基本流程后，进一步分析稠密三维重建流程内部的算法热点，在不同环节设计并行算法：在GPU单节点环境下实现稠密点云面片优化阶段并行设计和稠密点云特征处理阶段并行设计，同时建立GPU集群环境下二级并行模型；最后在稀疏重建环节，替换特征提取算子，改进相关流程，实现GPU环境下稀疏重建加速处理。

本发明提出了一种基于GPU集群的三维重建算法并行化方法，该方法包括以下步骤：

S1)给出基于图片特征的三维重建算法的具体流程；

S2)分析三维重建算法的具体流程中最耗时的部分；

S3)对三维重建流程中最耗时的部分进行并行优化；

S4)对三维重建流程中的其余部分进行优化。

步骤S3)中对三维重建流程中最耗时的部分，即稠密重建流程进行并行优化包括：

S31)稠密重建任务划分：当图像数量增加时，处理稠密重建流程PMVS算法软件的计算负担迅速增加，进而出现可扩展性的问题，对于内存的需求，也随着稠密重建过程中使用图像数量的增加而增加，并且实际限制了同时匹配的图像数量，图像聚类是这类内存问题的解决方法，也就是把大型项目拆分成为不同的块，MVS聚类视图(Clustering Viewsfor MVS，CMVS)的预处理步骤将图像集分解为重叠视图群集，以使稠密重建任务被划分为一个个占用内存较小的重建任务图像集，由SfM生成的稀疏点云被用于生成大小可控的重叠图像集群，使得目标生成点云中的每个三维点至少由一个图像任务群重构，分解完毕之后，再利用PMVS算法进行稠密点云生成，即基于CMVS/PMVS稠密点云重建的流程一般分为三步，如图3所示：

(1)任务拆分，确定子任务图像集；

(2)依据子任务图像集生成配置文件；

(3)根据配置文件进行PMVS稠密重建；

S32)CMVS并行算法实现：使用MPI并行编程技术，并采用数据的并行，将CMVS依据需求分割生成多个子图像任务集，并生成对应的任务信息文件option 0，option 1，…….option N，任务信息文件由Master按情况分发给各个Worker，Worker收到任务信息文件后开始各自调用PMVS稠密重建算法进行稠密点云生成，进行工作得到结果，其具体的实现过程如下：

(1)主节点和计算节点初始化以后，主节点开始预处理操作，这些操作包括配置错误日志文件，处理和校正输入参数，访问shell文件中所有的‘option’记录；

(2)根据用户选择将记录分解成一系列小的‘shell’文件，如果有需要的话，还需要检查‘option’文件的正确性；

(3)主节点将随后把必要的参数广播给其他计算节点；

(4)当所有的进程接收到相应的参数后，主节点和计算节点进入到处理和沟通环节，即每个计算节点在处理之前都先创建适合自己使用的正确的PMVS运行环境；

(5)当所有的子任务完成后，整个程序将终止。

1.主节点的运行原理

主节点的通讯部分用一个while循环实现，通过循环来给各个子节点分发任务，当没有子任务发送给计算节点时，循环终止。同时，主节点能够接收到任一计算节点的状态信息。如果主节点接收到一个计算节点的“空闲”标志信息，并且仍然有未处理的子任务，这时主节点将在新的子任务信息和“工作”这一更新状态信息一起发送给该计算节点。否则，主节点将发送“终止”命令，告诉计算节点退出循环并等待程序终止。如果发送标志信息中包括“结果发送”信息，主节点将收到此计算节点返回的结果。主节点收到的信息是Work节点处理重建结果的状态，子任务成功执行则状态为0，否则状态是不为0的数字。当所有子任务完成后，主节点退出循环。

注意以下几点：(1)因为主节点在初始化过程中为每一个子任务准备了任务文件(或者多个任务的批处理文件)，并且每个计算节点都有自己的运行环境，所以在1个子任务发送给一个计算节点之前，主节点没有必要为了避免不同的计算节点同时处理相同的子任务等到特定的时间段才开始另一项工作；(2)主节点不运行任何的子任务，因为它主要负责与计算节点的通讯从而获得整个集群的最佳负载平衡。也就是说在实际计算中共有N-1个进程；(3)不像没有负载均衡的方法，这个方法能使主节点控制每一个计算节点的状态并且向计算节点发送相应的指令。这样，无论一个分布式任务中配置文件的数量是1，2，……或n/N，主节点将同时处理它们，例如拆分子任务，收到处理结果和指定一个新任务给计算节点。

2.计算节点的运行原理

为了与主节点保持协调，计算节点的通讯部分也是一个while循环。一个计算节点首先向主节点发送当前状态信息，例如空闲、工作或发送结果。然后，计算节点等待主节点的指令，按照主节点的指令，计算节点将终止工作或继续工作。如果它处于工作状态，将收到一个任务文件并且开始运行稠密重建(PMVS)模块，除非计算节点正在处理任务。如果PMVS模块运行结束，最终获得了正确的结果或者由于某些错误而终止，结果将被发送给主节点，主节点将处理这些结果并且将相应的提示发送给用户。

3.调用串行PMVS的过程

在计算节点上，接收到的子任务可以是单独一个待处理的option文件，也可以是一种shell批处理配置文件的形式。接收完成后，如果是单独的option文件，则直接进行PMVS算法模块的调用，如果是shell批处理文件，首先检查shell文件形式，然后通过shell文件调用PMVS模块。

在落实到MPI算法实现过程中，通常需要以下几步：初始化MPI运行环境；获取每一个进程的序列号；MPI程序发送消息；MPI程序接收消息；MPI程序结束。开发人员可以根据以上各处理步骤之间加入自己需要的功能。程序流程图见图4。

如图4所示，其流程具体描述如下：

(1)中心节点(Master)初始化运行环境，获取数据目录，解析将要分配给子节点的具体任务；

(2)中心节点向计算节点广播(分配)数据路径，计算节点按照各自的路径在存储节点中找到各自的任务数据；

(3)计算节点处理数据，数据处理完成后，将结果写入存储数据，并将自身节点的运行状况反馈给中心节点；

(4)中心节点分析从计算节点反馈回来的消息，按照反馈消息和剩余任务数量决定是否继续向计算节点分配任务数据路径；

(5)依次重复(2)、(3)、(4)流程，直到所有任务分配；

(6)中心节点结束运行环境，结束程序运行。

S33)剖析PMVS稠密重建内部流程，依照实际算法运行顺序进行耗时热点分析，在此基础上，针对耗时热点尝试采用GPU加速处理，进一步加速PMVS稠密重建流程：PMVS算法具有三个步骤：图像特征处理，即匹配步骤；面片生成与扩散，即扩散步骤；过滤不正确匹配，即过滤步骤；在匹配步骤中，使用Harris算子等特征提取算法检测角点和团块特征，然后在多个图像中进行匹配，所得的特征匹配点作为下一步面片生成的基础；在扩散步骤中，从初始匹配，即稀疏面片开始，考虑对投影面片图像中的相邻像素进行扩散，这样可以通过扩散重建来创建密集的面片，如果相邻图像单元已经被重建或发生深度不连续的情况，那么不再进行扩散；在过滤步骤中，通过处理模型中的遮挡，可见性约束被随后用于过滤不正确的匹配，通过滤除异常块强制执行全局可见一致性，然后通过进一步的过滤删除仅在少数深度图上可见的面片；

S331)PMVS稠密重建的算法耗时热点在稠密重建的面片处理阶段，对其实现基于GPU环境下的加速处理：

经过前期特征处理后，对于每张图像中的任意特征点f，都能通过特征匹配找出其在其他图像中的对应同名点f′，从而构成一个集合F；将集合中的每个(f,f′)进行三角测量原理构建即可生成三维空间点，生成三维点之后将这个三维空间点作为一个可能的面片中心，尝试生成空间面片。面片(patch)是近似于物体表面的局部正切平面，包括中心c(p)，法向量n(p)，参考影像R(p)。面片的一条边与参考相机(即拍摄参考影像的相机)的x轴平行。patch上有μ×μ的格网(μ一般取5或7)。

构建过程如下：

其中，参考影像R(p)要求面片在参考影像中可见；O(I)为参考图像光心，I_i表示一个任务集合中的影像。

由于生成的面片可能会含有错误的点，所以还需要考虑面片向量与图像光心向量的相互关系，即这两者的夹角需要小于给定阈值τ，只有符合该条件，才认为面片对于相片，即影像所在的投影平面是可见的，即：

其中，V(p)为经过筛选后的参考影像集，I_m表示一个面片任务的影像集合。

对侯选面片的参数进行初始化后，首先利用面片优化方法对n(p)和c(p)进行优化，得到更为准确的面片及对应的图像集信息。在进行面片优化的过程中，主要依赖于光度差异函数，现详细说明如下：

假设物点X对应的面片为p(X)＝p，R(p)为面片p的参考图像，V(p)代表所有可以在照片中看见空间面片p的照片集合。面片p的像素点的个数为n，在面片p的可见图像I₁,I₂∈V(p)上的投影分别为pp₁和pp₂。如图5所示，面片p在I₁,I₂上的光度差异性被定义为h(p,I₁,I₂)：

h(p,I₁,I₂)＝1-NCC(pp₁,pp₂) (3)

其中，NCC(Normalized Cross Correlation)是两个面片投影之间的归一化互相关关系。pp₁(j)和pp₂(j)分别代表了面片p的第j个像素在图像I₁,I₂上的像素值。

和

分别代表了pp₁和pp₂的像素平均值。

利用h(p,I₁,I₂)，便可定义物点X对应面片p的总的光度差异函数：

这儿的\表示V(p)中除了R(p)剩下的集合。

为了降低干扰，算法还需要筛选出那些与参考图像R(p)的光度差异值低于给定阈值的图像，即

V^*(p)＝{I|I∈V(p),h(p,I,R(p))≤α} (6)

优化过程中更新的|V^*(p)|≥γ，若面片生成成功，则把p存储在对应的cell里，并同时更新以(i,j)为中心的生成面片Q(i,j)和其中符合面片可见条件的面片Q^*(i,j)。为了加快计算速度，规定：只要存在一个面片被重建成功，则这个存储面片的cell终止更新。

初始化面片生成之后，开始试图在面片领域内生成新的面片。给定面片p，首先确定一个满足一定规则的邻居图像cel1集C(p)，然后对每个cel1执行面片扩展策略。其中：

C(p)＝{C_i(x′,y′)|p∈Q_i(x,y),|x-x′|+|y-y′|＝1} (7)

其中，C_i(x′,y′)是由(x,y)的邻域(x′,y′)组成的cell集。

当面片p与片p′满足以下公式时，不需要扩散。并从C_i(x′,y′)中删除C(p)。

|(c(p)-c(p′))^*n(p)|+|(c(p)-c(p′))^*n(p′)|＜2ρ₁ (8)

其中，ρ₁是人为预先设定的误差。

如果不满足上述情况，那么对于每一个属于c(p)的图像块均尝试生成新面片，然后重复进行面片优化工作，最后判定V^*(p)。如果|V^*(p)|≥γ，γ为预先设定的阈值，就接受新扩散的面片，并相应的更新Q(x,y)和Q^*(x,y)。

在稠密重建面片生成与优化阶段需要进行大量的向量操作和图像矩阵间运算，而GPU设备因为拥有诸多图像处理核心(CU)，能够快速处理矩阵(向量)计算操作，因此将面片优化阶段的NCC函数运算和光度差异值计算流程迁移至GPU环境下加速处理；采用GPU并行基于OpenCL编程环境，OpenCL并行程序分为宿主机程序和内核程序两个部分，一个并行算法可能会包含一个或多个内核程序，采用的refinePatch并行算法的设计也是按宿主机程序和内核程序两个部分完成；运行在主机端的宿主机程序，除了要负责程序的初始化、参数设定和平台设备的获取等工作外，还要在设备端计算的发起与结束时承担主机端与设备端的内存空间映射和数据传输任务，待计算设备处理完成之后进行计算资源的释放工作，并将结果数据输出，该并行程序只设定了一个运行在GPU设备端的kernel程序，当宿主机向OpenCL计算设备发起计算任务时，kernel函数被映射到OpenCL计算设备的各个计算单元分别处理，各个计算单元运行kernel函数一次只处理一个像素/面片，处理完成再读取未处理的像素/面片进行处理，直至所有数据均得到处理，CPU在refinePatch并行算法中在完成计算任务的发起外，还承担一定量的计算工作，而GPU则承担了refinePatch并行算法相关函数计算及优化的主要计算任务，如图6所示，由此实现MPI+GPU二级稠密点云重建流程。

S332)Harris特征提取环节使用GPU进行加速处理：

Harris角点检测器，该技术主要基于图像的偏导数构造2×2矩阵，然后分析特征值。以图像中的一小块区域为例，我们的目标是确定该区域是否有个角点。因此，考虑所有相邻的区域并计算区域和所有相邻区域之间的强度差异。如果在所有方向上的差异都很大，那么我们就知道该区域中有一个角点。角是沿两个方向具有强烈强度差异的点，如图7所示。

在灰度图像G中我们使用一个窗口w(x,y)扫描该图像：在左右方向上位移u，在垂直方向上位移v)，(x,y)为窗口中心坐标，并将计算窗口滑动前后强度的变化。

处理化简可得：

其中，

在窗口大小在滑动期间固定时，我们只需要关注M：一个椭圆函数。其扁率和尺寸是由矩阵的特征值决定的，如图8所示。

可通过M矩阵的特征值来判断检测区域的属性：

(1)当M的两个特征值相差巨大，则说明检测区域内存在直线。

(2)当M的两个特征值近似相等，且都比较少，则检测区内为平面，没有关联特征；

(3)当M的两个特征值都比较大，且两个特征值数值相近，则可以判断为检测范围内存在角点。

根据二次项函数特征值的计算公式，我们可以求M(x,y)矩阵的特征值。但是Harris给出的角点检测方法并不需要计算具体的特征值，而是计算一个角点响应值R来判断角点。R的计算公式为：

R＝detM-α(trace M)² (11)

式中，detM为矩阵M的行列式；trace M为矩阵M的直迹；α为人为常数，取值范围为0.04～0.06。

在串行算法实际实现过程中算法的具体实现步骤可描述为：

(1)输入数据转换：根据输入图像I将图像数据从RGB色彩空间转化到灰度图像；

(2)计算图像G(x,y)在X和Y两个方向的梯度；

(3)计算图像两个方向梯度的乘积；

(4)使用高斯函数对

和G_xy进行高斯加权；

(5)计算每个像素的Harris响应值R，并对小于某一预设阈值t的R置为零；

R＝{R:detM-α(traceM)²＜t} (13)

Harris算法的耗时热点根据算法原理为：计算梯度与高斯加权等矩阵计算操作。因此，将针对这些地方进行并行优化设计改造，将算法这两步全部集中到GPU计算设备中进行处理处理，CPU负责算法初始化、数据读取、逻辑控制及最终的结果数据写入输出；

图9展示了整体并行化策略，主机端主要负责图像数据读取、数据跨设备间传输和核函数选定等初始化工作，还有负责后期GPU设备处理完成的GPU内存数据传输回主机端进行结果输出，GPU主要负责算法的主体计算任务。

如图10-12所示，下面开始介绍基于GPU的图像梯度计算并行化策略：可以拆分为三大块：

(1)在实际运行时，待处理图片的每个像素需要进行窗口滑动计算变化梯度，可以参考数据并行的原理进行并行化设计，即设计kernel核函数，将不同像素点的梯度计算任务分配给GPU中的不同线程，每个线程使用一个计算节点(CU)处理多个像素点，GPU多线程同时计算每个像素点的梯度；

(2)在得到图像的梯度计算矩阵之后，根据Harris计算原理及设计流程，需要进行梯度矩阵的几何计算，此环节也可依据GPU进行并行化设计，每个线程单独获取一个像素点的梯度计算，同时进行；

(3)高斯加权的并行化处理：高斯加权的并行化处理思路与(1)类似，即将待处理的矩阵拷贝进入GPU内存，所有线程共享矩阵数据，每个矩阵元素分配一个线程，进行单点高斯加权计算，实现方式仅在核函数上与(1)不同。

进行上述三次GPU加速处理的串行流程实际上是上下联接的，中间没有参杂其他操作，因此，我们考虑尝试将三次GPU加速合并到一次设备间传输当中，减少设备间数据拷贝次数，从而进行优化。

主机端在CPU环境下依然负责图像数据读取、数据跨设备间传输和核函数选定等初始化工作，还有负责后期GPU设备处理完成的GPU内存数据传输回主机端进行结果输出。GPU设备分别调用kernel核函数，计算梯度矩阵，几何计算与高斯加权并传回结果。

步骤S4)中对三维重建流程中的其余部分进行优化为对开源的Bundler处理稀疏重建流程进行优化：基于Bundler的稀疏重建大致可以分为3步：初始化环境，并进行特征提取；基于特征提取进行特征匹配；运行主程序，进行基于运动恢复结构原理的稀疏三维重建，Bundler程序中调用的特征提取算法是SIFT算法，其全称叫做尺度不变特征变换匹配算法。

如图13所示，主机端负责图像灰度化，高斯滤波卷积核生成等预处理过程。GPU设备负责将图像各像素点分配到计算节点中，图像以不同的尺度递增地与高斯函数进行卷积，并进行差分运算；然后，通过对比每一采样点及其在当前图像中的8个邻域，以及上、下尺度中的9个邻域来检测局部极值作为关键点；通过分析局部强度梯度的主要方向，为每个关键点分配一致的方向；最后对每个关键点周围的梯度幅值和方向进行取样，应用高斯加权函数生成特征描述子。SIFT算法高计算复杂度限制了它在大规模数据和实时系统中的进一步应用，在GPU上实现高性能SIFT，在SIFT各流程环节中，高斯滤波需要在设定的范围半径内进行所有像素点的卷积操作，卷积核固定的情况下，每个点的最终结果看起来会受到周围像素点的影响，但实际上在图像固定的情况下，每个像素点的值在运行前均已固定，各点的响应值仅与卷积核有关，相互独立。同样如此，极值点的检测环节、特征主方向确定环节、直方图统计环节与高斯滤波类似，看起来各点的结果与周围像素紧密相连，一旦周围的值产生变化，响应值就可能发生变化，实际上在运算过程中各点的依赖值仅仅是程序运算前的图像值，在程序开始运行时便已经固定不再变化，依旧可以互不干扰，适合并行化处理，在得到经过GPU优化后的SIFT特征处理之后，可将其带入Bundler处理流程，实现稀疏点云获取的GPU优化版本。

以上所述，仅为本发明的具体实施方式，本说明书中所公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换；所公开的所有特征、或所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以任何方式组合；本领域的技术人员根据本发明技术方案的技术特征所做出的任何非本质的添加、替换，均属于本发明的保护范围。

Claims

1.一种基于GPU集群的三维重建算法并行化方法，其特征在于，该方法包括以下步骤：

S1)给出基于图片特征的三维重建算法的具体流程；

S2)分析三维重建算法的具体流程中最耗时的部分；

S3)对三维重建流程中最耗时的部分进行并行优化；

S4)对三维重建流程中的其余部分进行优化；

其中，步骤S1)中基于图片特征的三维重建算法的具体流程包括：利用无人机获取到地物图像序列后，首先进行特征提取，在一系列不同的图像中识别同名关键点，正是这些关键点实现了不同图像间的匹配，以及场景几何的重构，一旦定位了每个图像中的关键点，就需要确定不同图像中关键点之间的对应关系，为了确保仅仅保留正确的对应关系，有必要采取进一步的步骤过滤掉剩余的任何错误的匹配，由于关键点仅限于具有几何一致性匹配的关键点，所以，每个图像对之间的链接可以被识别，并且组织为轨迹；然后运动恢复结构理论基于这些轨迹和关键点以及其他相关的图片信息，进行几何解算，可以粗略描述场景几何信息的稀疏点云，SfM产生的稀疏点云通常只是使用MVS生成更加密集点云数据的中间步骤；MVS的目标是从已知的相机参数的图像集合中提供完整的三维场景重建；选用开源的Bundler处理稀疏重建流程，选用PMVS算法软件处理稠密重建流程；

步骤S2)中分析三维重建算法的具体流程中最耗时的部分包括：在Linux环境下通过处理不同规模的数据集，记录各环节耗时情况，采用多次实验求平均的思想，分析三维重建算法的具体流程中的瓶颈环节，得到在基于图像序列的三维重建算法的具体流程中，稠密重建流程耗时占比最大，占三维重建算法的具体流程的60％以上。

S31)稠密重建任务划分：当图像数量增加时，处理稠密重建流程PMVS算法软件的计算负担迅速增加，进而出现可扩展性的问题，对于内存的需求，也随着稠密重建过程中使用图像数量的增加而增加，并且实际限制了同时匹配的图像数量，图像聚类是这类内存问题的解决方法，也就是把大型项目拆分成为不同的块，MVS聚类视图(Clustering Views forMVS，CMVS)的预处理步骤将图像集分解为重叠视图群集，以使稠密重建任务被划分为一个个占用内存较小的重建任务图像集，由SfM生成的稀疏点云被用于生成大小可控的重叠图像集群，使得目标生成点云中的每个三维点至少由一个图像任务群重构，分解完毕之后，再利用PMVS算法进行稠密点云生成，即基于CMVS和PMVS稠密点云重建的流程一般分为三步：

(1)任务拆分，确定子任务图像集；

(2)依据子任务图像集生成配置文件；

(3)根据配置文件进行PMVS稠密重建；

S32)CMVS并行算法实现：使用MPI并行编程技术，并采用数据的并行，将CMVS依据需求分割生成多个子图像任务集，并生成对应的任务信息文件option 0，option 1，…….option N，任务信息文件由Master按情况分发给各个Worker，Worker收到任务信息文件后开始各自调用稠密重建算法进行稠密点云生成，进行工作得到结果，其具体的实现过程如下：

(3)主节点将随后把必要的参数广播给其他计算节点；

(5)当所有的子任务完成后，整个程序将终止；

S33)剖析PMVS稠密重建内部流程，依照实际算法运行顺序进行耗时热点分析，在此基础上，针对耗时热点尝试采用GPU加速处理，进一步加速稠密重建流程：PMVS算法具有三个步骤：图像特征处理，即匹配步骤；面片生成与扩散，即扩散步骤；过滤不正确匹配，即过滤步骤；在匹配步骤中，使用Harris算子等特征提取算法检测角点和团块特征，然后在多个图像中进行匹配，所得的特征匹配点作为下一步面片生成的基础；在扩散步骤中，从初始匹配，即稀疏面片开始，考虑对投影面片图像中的相邻像素进行扩散，这样可以通过扩散重建来创建密集的面片，如果相邻图像单元已经被重建或发生深度不连续的情况，那么不再进行扩散；在过滤步骤中，通过处理模型中的遮挡，可见性约束被随后用于过滤不正确的匹配，通过滤除异常块强制执行全局可见一致性，然后通过进一步的过滤删除仅在少数深度图上可见的面片；

S331)稠密重建的算法耗时热点在稠密重建的面片处理阶段，对其实现基于GPU环境下的加速处理：

经过前期特征处理后，对于每张图像中的任意特征点f，都能通过特征匹配找出其在其他图像中的对应同名点f′，从而构成一个集合F；将集合中的每个(f,f′)进行三角测量原理构建即可生成三维空间点，生成三维点之后将这个三维空间点作为一个可能的面片中心，尝试生成空间面片，面片(patch)是近似于物体表面的局部正切平面，包括中心c(p)，法向量n(p)，参考影像R(p)。面片的一条边与参考相机，即拍摄参考影像的相机的x轴平行，patch上有μ×μ的格网；

构建过程如下：

其中，参考影像R(p)要求面片在参考影像中可见；O(I)为参考图像光心，I_i表示一个任务集合中的影像；

对侯选面片的参数进行初始化后，首先利用面片优化方法对n(p)和c(p)进行优化，得到更为准确的面片及对应的图像集信息，在进行面片优化的过程中，主要依赖于光度差异函数，现详细说明如下：

假设物点X对应的片面为p(X)＝p，R(p)为面片p的参考图像，V(p)代表所有可以在照片中看见空间面片p的照片集合，面片p的像素点的个数为n，在面片p的可见图像I₁,I₂∈V(p)上的投影分别为pp₁和pp₂，那么，面片p在I₁,I₂上的光度差异性被定义为h(p,I₁,I₂)：

h(p,I₁,I₂)＝1-NCC(pp₁,pp₂) (公式3)

和

分别代表了pp₁和pp₂的像素平均值；

V^*(p)＝{I|I∈V(p),h(p,I,R(p))≤α} (公式6)

优化过程中更新的|V^*(p)|≥γ，γ为预先设定的阈值，若面片生成成功，则把p存储在对应的cell里并同时更新以(i,j)为中心的生成面片Q(i,j)和其中符合面片可见条件的面片Q^*(i,j)，为了加快计算速度，只要存在一个面片被重建成功，则这个存储面片的cell终止更新；

初始化面片生成之后，开始试图在面片领域内生成新的面片，给定面片p，首先确定一个满足一定规则的邻居图像cel1集C(p)，然后对每个cel1执行面片扩展策略，其中：

C(p)＝{C_i(x′,y′)|p∈Q_i(x,y),|x-x′|+|y-y′|＝1} (公式7)

其中，C_i(x′,y′)是由(x,y)的邻域(x′,y′)组成的cell集；

当面片p与片p′满足以下公式时，不需要扩散，并从C_i(x′,y′)中删除C(p)，

|(c(p)-c(p′))^*n(p)|+|(c(p)-c(p′))^*n(p′|＜2ρ₁ (公式8)

其中，ρ₁是人为预先设定的误差；

如果不满足上述情况，那么对于每一个属于c(p)的图像块均尝试生成新面片，然后重复进行面片优化工作，最后判定V^*(p)，如果|V^*(p)|≥γ，γ为预先设定的阈值，就接受新扩散的面片，并相应的更新Q(x,y)和Q^*(x,y)；

在稠密重建面片生成与优化阶段需要进行大量的向量操作和图像矩阵间运算，而GPU设备因为拥有诸多图像处理核心(CU)，能够快速处理矩阵计算操作，因此将面片优化阶段的NCC函数运算和光度差异值计算流程迁移至GPU环境下加速处理；采用GPU并行基于OpenCL编程环境，OpenCL并行程序分为宿主机程序和内核程序两个部分，一个并行算法可能会包含一个或多个内核程序，采用的refinePatch并行算法的设计也是按宿主机程序和内核程序两个部分完成；运行在主机端的宿主机程序，除了要负责程序的初始化、参数设定和平台设备的获取等工作外，还要在设备端计算的发起与结束时承担主机端与设备端的内存空间映射和数据传输任务，待计算设备处理完成之后进行计算资源的释放工作，并将结果数据输出，该并行程序只设定了一个运行在GPU设备端的kernel程序，当宿主机向OpenCL计算设备发起计算任务时，kernel函数被映射到OpenCL计算设备的各个计算单元分别处理，各个计算单元运行kernel函数一次只处理一个像素或面片，处理完成再读取未处理的像素或面片进行处理，直至所有数据均得到处理，CPU在refinePatch并行算法中在完成计算任务的发起外，还承担一定量的计算工作，而GPU则承担了refinePatch并行算法相关函数计算及优化的主要计算任务，由此实现MPI+GPU二级稠密点云重建流程；

S332)Harris特征提取环节使用GPU进行加速处理：

Harris角点检测器，该技术主要基于图像的偏导数构造2×2矩阵，然后分析特征值，以图像中的一小块区域为例，目标是确定该区域是否有个角点，因此，考虑所有相邻的区域并计算区域和所有相邻区域之间的强度差异，如果在所有方向上的差异都很大，那么我们就知道该区域中有一个角点，角是沿两个方向具有强烈强度差异的点；

在灰度图像G中我们使用一个窗口w(x,y)扫描该图像：在左右方向上位移u，在垂直方向上位移v)，并将计算窗口滑动前后强度的变化：

处理化简可得：

其中，

在窗口大小在滑动期间固定时，我们只需要关注M：一个椭圆函数，其扁率和尺寸是由矩阵的特征值决定的。

可通过M矩阵的特征值来判断检测区域的属性：

(1)当M的两个特征值相差巨大，则说明检测区域内存在直线；

(3)当M的两个特征值都比较大，且两个特征值数值相近，则可以判断为检测范围内存在角点；

根据二次项函数特征值的计算公式，求M(x,y)矩阵的特征值，但是Harris给出的角点差别方法并不需要计算具体的特征值，而是计算一个角点响应值R来判断角点，R的计算公式为：

R＝det M-α(trace M)² (公式11)

式中，detM为矩阵M的行列式；trace M为矩阵M的直迹；α为预设常数；

在串行算法实际实现过程中算法的具体实现步骤可描述为：

(2)计算图像G(x,y)在X和Y两个方向的梯度；

(3)计算图像两个方向梯度的乘积；

(4)使用高斯函数对

和G_xy进行高斯加权；

(5)计算每个像素的Harris响应值R，并对小于某一阈值t的R置为零；

R＝{R:detM-α(traceM)²＜t} (公式13)

Harris算法的耗时热点根据算法原理为：计算梯度与高斯加权等矩阵计算操作；因此，将针对这些地方进行并行优化设计改造，将算法这两步全部集中到GPU计算设备中进行处理处理，CPU负责算法初始化、数据读取、逻辑控制及最终的结果数据写入输出；

主机端主要负责图像数据读取、数据跨设备间传输和核函数选定等初始化工作，还有负责后期GPU设备处理完成的GPU内存数据传输回主机端进行结果输出，GPU主要负责算法的主体计算任务；

下面开始介绍基于GPU的图像梯度计算并行化策略：可以拆分为三大块：

(1)在实际运行时，待处理图片的每个像素需要进行窗口滑动计算变化梯度，可以参考数据并行的原理进行并行化设计，即可考虑将不同像素点的梯度计算任务分配给GPU中的不同线程，GPU多线程同时计算每个像素点的梯度；

进行上述三次GPU加速处理的串行流程实际上是上下联接的，中间没有参杂其他操作，因此，将三次GPU加速合并到一次设备间传输当中，减少设备间数据拷贝次数，从而进行优化；

主机端在CPU环境下依然负责图像数据读取、数据跨设备间传输和核函数选定的初始化工作，还有负责后期GPU设备处理完成的GPU内存数据传输回主机端进行结果输出；GPU设备分别调用kernel核函数，计算梯度矩阵，几何计算与高斯加权并传回结果。

2.根据权利要求1所述的基于GPU集群的三维重建算法并行化方法，其特征在于，所述步骤S4)中对三维重建流程中的其余部分进行优化为对开源的Bundler处理稀疏重建流程进行优化：基于Bundler的稀疏重建大致可以分为3步：初始化环境，并进行特征提取；基于特征提取进行特征匹配；运行主程序，进行基于运动恢复结构原理的稀疏三维重建，Bundler程序中调用的特征提取算法是SIFT算法，其全称叫做尺度不变特征变换匹配算法；

SIFT算法高计算复杂度限制了它在大规模数据和实时系统中的进一步应用，在GPU上实现高性能SIFT，在SIFT各流程环节中，高斯滤波需要在设定的范围半径内进行所有像素点的卷积操作，卷积核固定的情况下，每个点的最终结果看起来会受到周围像素点的影响，但实际上在图像固定的情况下，每个像素点的值在运行前均已固定，各点的响应值仅与卷积核有关，相互独立。同样如此，极值点的检测环节、特征主方向确定环节、直方图统计环节与高斯滤波类似，看起来各点的结果与周围像素紧密相连，一旦周围的值产生变化，响应值就可能发生变化，实际上在运算过程中各点的依赖值仅仅是程序运算前的图像值，在程序开始运行时便已经固定不再变化，依旧可以互不干扰，适合并行化处理，在得到经过GPU优化后的SIFT特征处理之后，可将其带入Bundler处理流程，实现稀疏点云获取的GPU优化版本；主机端负责图像灰度化，高斯滤波卷积核生成的预处理过程；GPU设备负责将图像各像素点分配到计算节点中，图像以不同的尺度递增地与高斯函数进行卷积，并进行差分运算；然后，通过对比每一采样点及其在当前图像中的8个邻域，以及上、下尺度中的9个邻域来检测局部极值作为关键点；通过分析局部强度梯度的主要方向，为每个关键点分配一致的方向；最后对每个关键点周围的梯度幅值和方向进行取样，应用高斯加权函数生成特征描述子。

3.根据权利要求2所述的基于GPU集群的三维重建算法并行化方法，其特征在于，所述步骤S331)中μ取5或7。

4.根据权利要求3所述的基于GPU集群的三维重建算法并行化方法，其特征在于，所述步骤S332)中α取值范围为0.04-0.06。