基于GPU并行计算的图像快速匹配

宋 骥,周松涛

(武汉大学 国际软件学院，湖北 武汉 430079)

图像匹配在数字摄影测量、计算机视觉、遥感影像处理、医学图像、文件夹和数据库的图像搜索等领域有广泛应用.在现在文字检索已然成熟的时代，以“图”搜“图”已经成为下一步的需求.在军事领域搜寻目标时，甚至需要一张图中迅速找到众多目标的精确位置.

目前用于图像匹配的方法主要有基于图像灰度相关方法、基于图像特征方法、基于神经网络和遗传算法等人工智能方法.基于图像灰度的匹配算法简单，匹配准确率高，适合硬件实现，但是计算量大，不利于实时处理，对灰度变化、旋转、形变以及遮挡等比较敏感；基于图像特征的方法计算量相对较小，对灰度变化、形变及遮挡有较好的适应性，但是匹配精度不高；基于神经网络、遗传算法等人工智能的方法发展较晚，算法还不成熟[1].最近几年在GPU基础上发展起来的CUDA并行计算技术，在高性能计算方面得到了广泛应用，在图像匹配上也用于图像灰度相关算法的加速[2]，虽得到了一定的加速比，但是目前这种算法设计的效率还不够高，当数据规模增大时对问题的适应性下降.本文基于图像灰度相关算法，提出一种更加优秀、高效的CUDA并行算法，使得在原来的基础上进一步加速，并且对问题规模的继续加大有很好的支持和适应.

图1 图像匹配过程Fig.1 Process of image matching

图2 相邻矩阵重复情况Fig.2 The repeat in neighborhood matrices

1 问题描述

假设需要在待搜索图像中匹配给定目标图像.一般来说，目标图像较小(设为m×n)，待搜索图像尺寸较大(设为(M+m)×(N+n)).图像规格是0-255的8位灰度图像，如果是彩色图像，则先灰度化.

如图1，在待搜索图像中枚举矩阵左上角(x,y)，逐一检查待匹配的m×n矩阵，将其与目标图像进行匹配运算，若相关系数大于设定的阈值，则认为该矩阵匹配成功.

记小图像为g，大图像为S，用Sx, y表示S中以(x,y)为左上角点与g大小相同的子块.将两个相同尺寸的矩阵的相关系数记为ρ(x,y).

将协方差和方差公式代入并化简[1]可得公式1：

(1)

2 算法分析与CPU优化

通过优化，只需要一步运算(O(1)的时间)便可以完成相邻像素点的递推，式(1)中求平均值也可一步运算解决.如果搜索整个图像，则式(1)的分母部分，以及分子中平均值相乘部分，基于CPU的串行计算的整体算法复杂度在O(M×N)，这已经完全可以被CPU胜任，不再成为算法的瓶颈.

而式(1)的分子中，S与g相乘并求和的部分的算法复杂度依然是O(M×N×m×n)，这将导致该算法在CPU上的计算很慢，该部分已经成为整个算法继续加速和优化的瓶颈所在.

为描述简便，将Sx, y与g相乘并求和的结果记为Sum(x,y)，用来表示以(x,y)为左上角的m×n大小的像素矩阵与目标矩阵对应像素乘积之和.Sum(x,y)的计算成为问题的关键，因此下面的GPU算法主要针对Sum(x,y)进行计算.

(2)

3 利用GPU的初步并行算法

现代的显示芯片已经具有高度的可程序化能力、相当高的内存带宽、以及大量的执行单元，因此可利用显示芯片来帮助进行一些计算工作，即 GPGPU[3].虽然GPU并没有CPU所具备的强大的复杂逻辑指令的处理能力，却非常适合胜任大量、简单、重复的运算工作，以及高度并行化的并行计算.图形硬件有强大的并行性，支撑高密度的运算，具有很高的浮点运算速度，因此采用图形硬件主要是为了加速运算.

根据GPU的特点可以看出，式(2)是一个简单、重复而大量的运算工作，十分适合并行处理.基于CUDA用GPU的并行计算来解决图像匹配问题，已经得到了初步的应用[2].

3.1 CUDA模型简介

CUDA是NVIDIA的GPGPU模型，它使用C语言为基础.在CUDA架构下，显示芯片执行时的最小单位是thread，数个thread可以组成一个block[3].

从CUDA程序逻辑上看，采用两层并行，每个block内的thread细粒度并行处理，并且用共享寄存器进行通信，而各个block之间粗粒度并行，并行度视硬件情况而定.两层并行的好处之一，在于实现了可扩展的程序模型.当硬件的运算内核数量增加时，CUDA可以自适应的方式，自动增大粗粒度并行度(即同时运算的block数量)[4].好处之二在于细粒度并行中，不同线程(thread)之间可以通过共享寄存器进行通信，达到优化数据传输速度和加速计算的目的.

图3 算法2的分块方式及一个block内部线程分配Fig.3 Designs of blocks and threads in algorithm 2

图4 使用共享存储器的初步设想Fig.4 Abstract design of shared memory

3.2 初步的并行策略

通过对CUDA的并行机制的分析，做出以下分块策略[2].首先，将一个Sum(x,y)的计算作为一个block，并使block(x,y)计算Sum(x,y)，此时共有M×N个block，根据硬件条件若干个block并行计算，直到M×N个Sum(x,y)计算完毕.而后将每一block分配m×n个线程，每个线程在S和g中各取一个像素点的灰度并将其相乘累加，计算一个Sum(x,y).

初步的并行方法的速度已经较CPU有很大的提升，但是该并行方法存在以下诸多缺点：

1)数据重复调度：在GPU中计算，必须先把内存中的数据拷贝到显存中，然后访问显存进行计算，而访问显存需要消耗大量的时钟周期，所以设计时应该尽量减少对显存的访问次数.CUDA在显卡端(设备端)提供了一些可以高速缓存的寄存器，包括register和共享寄存器(shared memory)，可以将显存中的数据放到共享寄存器中，以此减少数据传输时间.

在上述分块方法中，每一个像素对应的block之间粗粒度并行，并不能用共享寄存器共享数据，所以此时每一个block都必须把m×n的矩阵数据从显存中读取出来.然而相邻像素之间(也即相邻block之间)访问的显存数据大量重复(目标矩阵重复m×n，即完全重复，待搜索矩阵重复(m-1)×n).

2)数据尺寸适应性不强：在目前的主流GPU中，每一个block所能同时管理的线程数量是有上限的(512)，当目标矩阵的尺寸m×n大于线程上限时，线程数将不够，无法使得每一个像素对应一个线程.此时必须修改，让每一个线程计算多个像素.而同时，当目标矩阵的尺寸变大时，在缺点1)中所阐述的数据的重复量将变大，目标矩阵尺寸越大，数据重复问题越严重.在这样的情况下，该算法对于数据规模的适应性大打折扣.

4 高效的GPU并行算法

4.1 改变分块方式

针对上述缺点对分块方式进行改进，将连续的BLOCKSIZE*BLOCKSIZE(16×16)数量个Sum(x,y)的计算放到一个block中(形成一个方阵)，同时使一个block包含16×16个线程(参见图3)，从而使每一个线程负责计算一个Sum(x,y).此时整个待搜索矩阵分块为(M/16)×(N/16)个block，每个block计算16*16个Sum(x,y)(无法整除情况可填补像素).

这样的好处在于可以同时解决数据重复调度和尺寸适应性问题.

首先，连续的16×16个线程处于同一个block中，他们可以通过共享寄存器共享数据.在原来的方法中，16×16个Sum(x,y)的计算需要调度显存共(m×n×2)×(16×16)次，而通过共享数据的方式，在所有线程开始阶段，可以将他们共同使用的待搜索图像的(16+m)×(16+n)灰度数据以及目标图像(m×n)读入，栅栏同步(barrier synchronization)后，各个线程开始计算自己负责的m×n的矩阵，算出一个Sum(x,y)，但是计算过程中可以共享地访问由前面读入到共享寄存器内的数据.数据调度效率大约是原来的16×16倍.(参考图4)

其次，每一个线程处理一个m×n的矩阵的计算，得到一个Sum(x,y)，当m×n增大时，不会受到线程数量限制的影响.

4.2 共享寄存器的处理方法

虽然理论上可以通过共享寄存器一次性读入数据，但是由于目前硬件限制，每一个SM中的共享寄存器空间被限制为一个上限(16 KB)，所以不能一次性读入和计算.该限制可以用人工分块的方法(blocking)解决.通过分块的方法，分块、分次读入，同时也分块、分次计算.

假设目标图像的尺寸m×n= 48×48，则将其分解成16×16的小块，一共3×3=9个小块.在一个block内同步的所有线程中，枚举这9个小块，处理完所有的小块，则累加计算完成，Sum(x,y)也就计算完毕了.处理一个小块的方法和1中所述的类似，对于一个16×16的小块，先将目标矩阵g对应的16×16个像素灰度值读入共享寄存器，而对于待搜索图像的数据此时则需要读入32×32个像素的数据(如图(5)，对于16×16个线程，Sum(x,y)对应的16×16的小块是从(x+i,y+j)到(x+i+15,y+j+15)，Sum(x+15,y+15)对应的16×16的小块是从(x+i+15,y+j+15)到(x+i+30,y+j+30)，其他Sum所用数据均在这31×31之中，但是为了满足合并访问，此处读入32×32).栅栏同步(barrier synchronization)后，每个线程并行开始各自的16×16次的相乘计算，并且累加到用于加和的变量中(该变量是每个线程私有的，定义为register寄存器类型).计算完毕后再次栅栏同步，而后进入下一个16×16的小块的计算，重复上述过程(参见图5).

在改进算法的人工分块(blocking)方法中，相邻共享寄存器覆盖的灰度数据还会出现重复的情况(如图6).所以在共享寄存器的读取数据时，可以采取将上一次共享寄存器的右边16×32个灰度读出来，直接移动到这一次共享寄存器的左边16×32的位置中，而不用重新调度显存中的这16×32个数据，只需要调度右边16×32个位置(首次从显存中读入)对应的灰度(如图6).这样进一步提高共享寄存器的利用率，也进一步减少对显存的重复访问.

图5 人工分块中共享存储器的处理流程与方法 图6 重复访问的共享存储器的处理Fig.5 The method of blocking using shared memory Fig.6 The method to solve repeat problem

4.3 改进后的并行算法的优点

1)充分利用共享寄存器的作用，降低显存访问次数，速度大幅提升.

2)不会受到目标图像尺寸的影响，目标图像尺寸增大时，不受限制，适应性好.

3)不会降低粗粒度层面(block)的并行度，虽然block数量下降(下降16×16倍，但是由于待搜索图像的尺寸巨大和海量的特点，block数量依然足够，而且不妨碍多GPU并行.

4)优化计算指令，使用CUDA内置整数乘法函数，达到最高运算速度.

5 实验

实验所使用的硬件：CPU为AMD Turion 64×2，主频2.2GHz.显卡GeForce 8600M GS，多处理器数目4，CUDA核心数目32，显存512M[4].

实验数据如表1所示，记初步GPU并行算法为算法1；改进后的高效算法为算法2；用时均指Kernel时间；加速比均相较CPU用时.

CPU优化中，已经将除开式(2)的其他部分做了优化，该部分用时在实验数据中均在10ms左右，已经不对最后算法时间造成过大影响，所以不参与比较.该表格中CPU用时是指用CPU的方法运算式(2)的所用时间.

表1 三种算法的匹配时间比较

实验数据分析：

初步的GPU并行计算速度大概是CPU的5倍左右，数据尺寸越大，由于并行度提高，加速比有所提升.在数据较小时，计算量不大，算法1的重复访问内存的影响凸显，在数据增大时，数据延迟的不足虽然被大计算量和大并行度稍微掩盖，但是依然有不少的影响.

由于初步的GPU算法受限于传输数据，所以没有达到最佳效率.以第1个数据为例，根据硬件条件和计算公式[5]，理论带宽为：Theoretical Bandwidth=(内存时钟×106×(内存接口宽度/8) ×2)/109=12.8GB/s(GeForce 8600M GS内存时钟400 MHz，内存接口宽度128 bit).实际宽度为：Effective bandwidth = ((Br+Bw)/109)/time经过计算，算法1的实际带宽为：(16×16×512×512×4)/10^9/time=0.67 GB/s，远远小于理论带宽，需要改进.而用改进后的高效并行算法2的实际带宽为：(16×16×512×512×4)/10^9/time=9.2GB/s，已经几乎达到理论带宽，此时改进后的算法的瓶颈已经不会出现在数据传输上，而是计算指令上.

从数据可以看出，在目前主流软硬件条件下，改进的算法是原有算法速度的7～12倍，大大提高了效率.

6 结语

本文提出一种高效GPU并行算法，避免了原有算法过多使用低速存储器所带来的问题，利用高速共享寄存器大幅提高效率，进一步加快图像匹配问题的求解速度.

现今，GPU并行计算的高速性为各个领域解决复杂而海量的计算问题提供了一种通用方法.这种方法门槛低、见效快、效果好.但是，在设计CUDA并行计算算法时，应综合考虑计算指令优化和数据调度传输优化，充分利用粗粒度并行带来的扩展性和细粒度并行中的数据共享通讯特性，将数据调度传输问题解决，不让其成为瓶颈.

现在若需要将问题进一步加速，关键则在于以下三个方面：新型的问题求解方法和计算公式，计算机硬件本身的发展和加速，以及运用多设备和设备集群等大规模联合计算，甚至是云计算.

[1] 李卓，邱慧娟.基于相关系数的快速图像匹配研究[J].北京理工大学学报，2007，27 (11)：998-1000.

[2] 肖汉，张祖勋.基于GPGPU的并行影像匹配算法[J].测绘学报，2010，39(1):46-51.

[3] Hotball.深入浅出谈CUDA[OL].[2008-06-04].http://www.kimicat.com/cuda简介.

[4] NVIDIA.CUDA C Programming Guide 3.2[OL].[2010-10-22].http://developer.download.nvidia.com/compute/cuda/3_2_prod/toolkit/docs/CUDA_C_Programming_Guide.pdf.

[5] NVIDIA.CUDA C Best Practices Guide 3.2[OL].[2010-8-20].http://developer.download.nvidia.com/compute/cuda/3_2_prod/toolkit/docs/CUDA_C_Best_Practices_Guide.pdf.