一种适应无人机平台的快速立体匹配方法
2014-06-07杨靖宇张永生
于 英, 杨靖宇, 张永生, 薛 武
(1.信息工程大学测绘学院,河南郑州450052; 2.中国人民解放军61733部队,河南郑州 450052;3.江西省数字国土重点实验室,江西抚州 344000)
一种适应无人机平台的快速立体匹配方法
于 英1,3, 杨靖宇2, 张永生1, 薛 武1,3
(1.信息工程大学测绘学院,河南郑州450052; 2.中国人民解放军61733部队,河南郑州 450052;3.江西省数字国土重点实验室,江西抚州 344000)
设计了一种高质量和快速的立体匹配方法,首先在粗分辨率的影像上获取立体影像的重叠区域和初始视差值,而后采用互信息作为匹配测度并采用一种通过一维路径聚合的全局优化方法得到精确的视差值。匹配方法采用图形处理器进行了加速,得到了很好的加速比。实验结果表明了该匹配方法的正确性和高效率性。
计量学;无人机;立体匹配;互信息;半全局;图形处理器并行
1 引 言
无人机低空摄影具有快速、灵活、低成本、高影像分辨率等特点,可在云下遥感的能力弥补了卫星光学遥感和普通航空摄影经常受云层遮挡获取不到影像的缺陷,相应的图像处理技术也得到了广泛的关注和重视。无人机图像处理中核心的技术之一是影像的快速立体匹配技术,因此研究一种适合无人机平台的快速立体匹配方法很有必要[1,2]。
本文研究的匹配方法要求可以生成稠密的视差图,因此特指基于灰度的匹配算法。基于灰度的匹配算法又可分为局部优化算法和全局优化算法。局部优化算法是依据匹配测度,计算过程比较简单,但匹配的可靠性比较差。而全局优化算法则是通过选用一个全局能量函数,并最小化该函数来得到更为准确的视差图[3]。图割法算法和信任传播算法是目前公认的效果最好的全局优化算法,可以获取高精度的稠密视差图,但是这两种算法的时间复杂度比较高,计算效率低[5~9]。德国宇航中心提出并采用的半全局匹配算法,不仅可以得到与图割法、置信传播法相媲美的处理结果,执行效率远高于这些算法,且半全局匹配方法具有较高的边缘保持性能[10]。全局匹配算法保证了匹配的稳定性和匹配的精度,但匹配的速度就难以保证。多核CPU和多核图形处理器(GPU:Graphic Processing Unit)的出现意味着并行系统已逐渐成为主流的处理器芯片,特别是GPU在处理能力和存储器带宽上相对CPU有明显的优势,使得GPU-CPU协同处理架构成为优异的高速计算平台[11],本文的匹配算法将采用并行的处理方式提高匹配的实时性以适应无人机平台对匹配速度的要求。
2 匹配算法
由于无人机上固定基线的立体相机的外方位元素精确已知,因此可以很容易地得到核线影像,这样本文匹配算法输入的影像均是核线影像。如图1所示,本文匹配方法采用降采样的方法生成粗分辨率的影像,在粗分辨率的影像上采用灰度相关匹配的方法计算影像的重叠度和初始视差值。后续的匹配只在重叠区域中进行,这样不仅使计算量成倍地较少,还有效地防止了图像非重叠区域中的信息对算法的干扰,提高了算法的精度;在后面精确的同名点搜索的过程中,搜索范围仅需要在初始视差值一定范围内进行,可以提高匹配的效率。然后,在原始分辨率的影像上按行计算匹配代价计算并采用半全局匹配的方法在利用多个方向的一维平滑约束来近似一个二维平滑约束得到最优的视差值,其中在匹配代价计算和匹配代价聚合采用了GPU并行计算进行了加速。最后,进行左右一致性检查和中值滤波,确保匹配结果的唯一性和误差的剔除。
图1 本文匹配算法计算流程图
2.1 匹配代价计算
2幅图像的互信息可以定义为:
式中:H(A)和H(B)分别为图像A和B的平均信息量,而H(A,B)则是它们的相关平均信息量。图像A和B的平均信息量和相关平均信息量为:
式中:ρA(a)和ρB(b)分别为图像A中具有灰度值a的概率密度函数和图像B中具有灰度值b的概率密度函数;ρA,B(a,b)为图像A和B的联合概率密度函数,它们可以由2幅图像A和B的联合直方图计算得到。
互信息并不直接依赖于灰度值来衡量不同图像的一致程度,而是依赖于它们在每幅图像中各自发生的概率和两幅图像组合产生的联合发生概率.因此它对灰度改变或一对一的灰度变换不敏感,能同时处理积极的和消极的图像灰度相互关系[3]。
2.2 匹配代价聚合
为由于防止噪声引起的匹配代价计算误差以及由此引起的深度污染扩散,将引进一个额外的约束加入到能量函数中:
式中:第1项表示所有像素点的匹配代价之和;第2项和第3项分别利用系数P1和P2对像素点p与其邻域内的像素点深度差存在的较小变化和较大变化2种情况进行了惩罚,即平滑约束,显然P1〈P2。T为判断函数,当且仅当其参数为真时函数值为1,否则为0。
对于二维图像,寻找式(3)的全局最小值已被证明是NP完全问题(Non-deterministic Polynomial:NP,多项式复杂程度的非确定性问题),而一维路径上的能量最小化则可以使用动态规划算法高效实现,但经典的动态规划算法只能沿着扫描行进行一维优化,使得匹配结果产生拖尾效应,因此在匹配算法中利用8个(或16个)方向上的一维平滑约束来近似拟合一个二维平滑约束。
在每一条路径L上依据动态规划的思想按照后面式(4)、式(5)进行计算,如图2所示。
式(4)中的第1项表示对像素点p赋予深度d的匹配代价;第2项是当前路径上p的上一个点p-r包含了惩罚系数的最小匹配代价;第3项则对最优路径的产生没有施加影响,加入这一项的目的仅仅是为了防止L过大,使得L≤Cmax+P2。其中r为路径;Cmax为最大匹配代价值;P2为视呈跳变惩罚值。
将各个方向上的匹配代价相加形成为总的匹配代价,如式(5)、图3所示。
图2 最小匹配代价路径Lr(p,d)
图3 8个方向上的匹配代价聚合
在通过匹配代价聚合更新S(p,d)得到所有像点对的匹配代价后,视差图的确定就是一个简单的选择过程:基准图像上的每个像点p的视差dp=min[S(p,d)],即对应的总匹配代价最小的视差值;而参考图像Im中的每个像素点q对应的视差dm=min[S(emb(q,d),d]。通过比较Db和Dm(Db和Dm分别为判断遮挡和错误匹配的依据),即进行一致性检查,式(7)。若二者差别在阈值之内,则认为匹配正确;否则若两者差别过大则不接受算法输出结果,并将该点标识为误匹配点Dinv。综上所述,匹配算法通过对视差值的不同变化施以相应的惩罚系数保证了平滑性约束,通过8或16个方向的一维路径动态规划来拟合二维平滑约束使得结果更加可靠,对噪声表现出一定的鲁棒性[10]。
2.3 GPU并行处理
在全局优化匹配算法的匹配代价立方体生成过程中,每个立方体元素的计算过程完全独立,具有很高的并行性,且每一层的视差空间图像的生成又存在着大量的局部性重复访问,为提高GPU的计算效率,需要对影像进行纹理存储器优化和共享存储器优化。
而在匹配代价的聚合过程中,每个像素的计算过程不再是相互独立的,而是需要利用当前路径上前一像点的匹配代价,这使得传统的按照对图像矩阵进行分块的模式不再适用,而应该按照扫描行(或列)对匹配代价立方体进行分块。如图4所示,在计算从上到下3个方向的Lr(p,d)时,图像的每一行被分割为N=Iw/Ap段,每一段由一个包含Ap×Rd个线程的线程块来进行计算,Iw为扫描宽度,Ap为计算过程中每个线程要循环执行Ll(扫描行数)次来扫描整幅图像。
图4 匹配代价聚合的GPU并行化
此外为进一步加快处理速度,在一次扫描过程中,同时计算多个方向的匹配代价聚合,如图5所示仅需要扫描两遍(从上而下和从下而上)即可对多个方向的匹配代价进行聚合。同时匹配代价立方体的中间计算子块被直接写入到全局存储空间,并在共享存储空间内的保留其副本数据,避免在下一行聚合过程中需要这些数据时,再从访问缓慢的全局存储空间读取数据。
在通过匹配代价聚合更新S(p,d)后,就可以简单的最小值选择来生成视差图。
3 实 验
图5 匹配代价聚合的快速扫描方法
为验证本文匹配算法的正确性,本文以UCD航空数字相机获取的湖北宝应地区的大重叠度面阵影像进行实验,如图6所示,影像大小为3000× 3000像素,像素大小为9.0μm,焦距为101.4000 mm。
匹配质量实验:分别采用匹配窗口大小为5× 5的相关系数法和本文方法对上面的核线影像进行了匹配,匹配结果如图7所示,可以看出本文方法得到的匹配视差图的质量明显高于相关系数法。
图6 本文所用核线影像
图7 匹配结果视差图
GPU匹配速度实验:实验中采用的GPU为NVIDIA Quadro Fx580显卡,CPU为Inter(R)Core--(TM)i5 2.8GHz×4核,系统内存为4GB。分别设置左右视差的搜索范围为128、256像素,在左影像上选择不同大小的多个待匹配的区域,统计相应CPU与GPU处理时间(ms)和加速效率(表1)。
4 结 论
针对无人机立体测绘的实际需求,设计了一种高效的立体影像匹配方法,实验结果表明该匹配方法在匹配质量和匹配速度上均取得了满意的结果。实现了2张像片的全局并行快速的匹配,后续将会研究实现序列图像的快速匹配方法(多视匹配),以更好地适应无人机大批量影像数据快速处理的需求。
表1 CPU与GPU处理时间及加速比
[1] 吴正鹏.无人机载双相机低空遥感系统应用初探[J].城市勘测,2011,(1):76-80.
[2] 胡堃.基于无人机遥感平台的震后灾情监测系统[J].科学探索与知识创新,2009,(1):100-101.
[3] 王珺.计算机立体视觉算法研究与实现[D].大连:大连理工大学,2009.
[4] 周文晖.一种鲁棒的基于互信息的实时立体匹配算法[J].遥感技术学报,2006,19(4):1243-1249.
[5] Boykow Y,Veksler O,Zabih R.Fastapproximate energy minimization via graph cuts[J].IEEETransactionson PatternAnalysisandMachineIntelligence,2001,23(11):1222-1239.
[6] Sun J,Zheng N N,Shum H Y.Stereo matching using belief propagation[J].IEEETranscationsonPattern AnalysisandMachineIntelligence,2003,25:787-800.
[7] Felzenszwalb P F,Huttenlocher D P.Efficient Belief Propagation for Early Vision[J].InternationalJournalof ComputerVision,2006.
[8] Forstmann S,Ohya J,Kanou Y,etal.Real-time stereo by using dynamic programming[C]//Proc of CVPR Workshop on Real-time 3D Sensors and Their Use,2004.
[9] Gong M L,Yang Y H.Fast stereo matching using reliability-based dynamic programming and consistency constraints[J].ProcedingsoftheNithIEEEInternational ConferenceonComputerVision,2003,(1):610-617.
[10] Egnal G.Mutual information as a stereo correspondence measure[R].Comp and Inf Science,U of Pennsylvania,2000.
[11] 张舒,褶艳利.GPU高性能运算之CUDA[M].北京:中国水利水电出版社,2009.
A Quick Stereo Matching Algorithm for Helicopter Platform
YU Ying1,3, YANG Jing-yu2, ZHANG Yong-sheng3, XUEWu1,3
(1.Institute of Surveying and Mapping,Information Engineering University,Zhengzhou,Henan 450052,China;
2.Troops 61733,The Chinese People's Liberation Army,Zhengzhou,Henan 450052,China;
3.Jiangxi Province Key Lab for Digital Land,Fuzhou,Jiangxi344000,China)
A stereo matching algorithm which is good matching quality and quick is designed.Firstly,an initial disparity and overlap based on coarse-resolution images is got.Then amutual information is used asmatching-measure and through one dimension path aggregation the accurate disparity will be got.Thematching algorithm uses GPU to accelerate the process,and a good acceleration ratio will be achieved.Finally,the experiment prove it right and efficient.
Metrology;Helicopter;Stereo matching;Mutual information;Semi-global;Graphic processing unit parallel
TB96
A
1000-1158(2014)02-0143-04
10.3969/j.issn.1000-1158.2014.02.10
2012-03-05;
2012-08-13
江西省数字国土重点实验室开放基金(DLLJ201303)
于英(1985-),男,河南郑州人,信息工程大学在读博士研究生,主要研究方向为基于移动平台的视频图像快速处理方法。yuying5559104@163.com