本发明涉及图像配准技术领域，尤其涉及一种基于gpu的图像快速配准方法。
背景技术：
由于gpu的灵活性和低成本，gpu的可用范围已经扩展到信号处理，计算机视觉，计算几何和科学计算。然而，当gpu用于实现通用并行计算时，存在一些限制。首先，因为gpu最初是为绘图而设计的，所以用户必须具备图形编程方面的知识，并将通用计算转换为图形计算。其次，由于并行性是通过在不同的数据元素上执行相同的程序来实现的。并行计算领域涉及的同步、互斥和原子操作不能有效地实现。幸运的是，在cuda架构出现之后，上面提到的这些限制已经得到了很好的解决。cuda架构是专门为通用gpu并行计算而设计的。与其他gpu编程系统(如开放图形库、microsoft’sdirectx11等)相比，最明显的优点是cuda使用c语言，因此原来为cpu编写的c函数通常通过很少修改就可以移植到cuda内核上运行。在cuda架构下，一个完整的程序分为两部分：一部分运行在cpu(所谓的主机)上。为了实现并行性，cuda允许数千线程(最小执行单元)在设备上同时运行。在cuda中，一定数量的线程被分组到一个块中，并且一定数量的块按照规则的网格模式被逻辑地排列，如图1所示。
每个块被映射到单个多处理器，并且一个多处理器可以同时运行多个线程块。由于本地资源(寄存器和共享存储器)在块之间分割，在相同块中的线程可以访问相同的共享存储器并直接实现同步。除了本地寄存器和共享内存之外，所有的线程都可以访问全局内存、常量内存和纹理内存。这就是cuda专门为通用计算而设计的原因。另外不同的多处理器之间的通信是通过使用位于板上的全局存储器实现的，必须使用全局存储器来在主机和设备之间交换数据。成千上万的线程可以通过一个单一的执行单元来启动，但并不是所有的线程都是同时执行的。实际上，通常一个大小为32的线程是gpu作为一个整体可以同时执行的最大的一组线程。在确定块的大小时，应确保块中的线程数必须是32的倍数。
目前，图像配准主要有如下几种方案：
1)基于图像灰度的配准方法。通常直接利用整幅图像的灰度信息，建立两幅图像之间的相似性度量，然后采用某种搜索方法，寻找使相似性度量值最大或最小的变换模型的参数值。但是，该方法计算量大，难以达到实时性要求，而且一旦进入信息贫乏的区域，会导致误匹配率的上升。
2)基于图像特征的配准方法。该方法需要对图像进行预处理，然后提取图像中保持不变的特征，如边缘点、闭区域的中心、线特征、面特征、矩特征等，作为两幅图像配准的参考信息。但是，由于其不依赖于图像的灰度信息，这种方法对特征提取和特征匹配的错误十分敏感，匹配性能依赖于特征提取的质量，需要可靠的特征提取和鲁棒的特征一致性。
3)基于模型的匹配方法。该方法在计算机视觉领域中的应用非常广泛，它可以分成刚体形状匹配和变形模板匹配两大类。kass提出的snake主动轮廓模型是比较典型的自由式变形模板。由于不受全局结构的限制，所以snake模型能表示任意的形状，但是，该模板对于模板的初始位置和噪声比较敏感，对于凹边缘的收敛性比较差，而且容易陷入局部最小值。
技术实现要素：
本发明的目的是提供一种基于gpu的图像快速配准方法，可以实现图像快速配准，并确保配准时的精度与准确度。
本发明的目的是通过以下技术方案实现的：
一种基于gpu的图像快速配准方法，包括：
对源图像进行仿射变换，获得变换图像，并分别计算源图像与变换图像的特征向量集合，分别记为da和db，特征向量集合da和db中的特征向量所对应的关键点记为集合a和b；
利用两个cuda内核在gpu上来分别实现关键点的粗匹配与精细匹配：
粗匹配时，块的数量由特征向量集合da的数量n确定，每个块中线程数量由特征向量集合db的数量m确定，利用上述方式来组织内核，每个块实现关键点集合a中每一个关键点与关键点集合b中所有关键点之间的匹配，匹配时计算相应特征向量之间的距离，并配合第一阈值，进行匹配，所有块的匹配结果存储在全局存储器中，共有n个匹配点对；
精细匹配时，将内核设计为执行ransac算法的每个元素，且只启动一个块；选择三个非线性的匹配点对由cpu计算仿射变换时所使用的变换矩阵h；所启动的一个块中包含n-3个线程，利用这n-3个线程来结合第二阈值、变换矩阵h与其余的n-3个匹配点对，来区分匹配点对属于内点匹配或外点匹配，从而筛选出精确的匹配点对。
由上述本发明提供的技术方案可以看出，在进行并行匹配中，启动两个内核来实现特征匹配的并行化，减少了运行时间，提高了配准时的精度与准确度。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。
图1为本发明背景技术提供的cuda中线程的组织和内存访问示意图；
图2为本发明实施例提供的一种基于gpu的图像快速配准方法的流程图。
具体实施方式
下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。
本发明实施例提供一种基于gpu的图像快速配准方法，如图2所示，其主要包括如下步骤：
步骤1、对源图像进行仿射变换，获得变换图像，并分别计算源图像与变换图像的特征向量集合。
本发明实施例中，仿射变换主要包括：平移、旋转、缩放和剪切变换。
本发明实施例中，设p＝(x,y)t，p'＝(x',y')t分别是源图像与变换图像的坐标，二者关系为：
p'＝hp+t；
写成矩阵形式为：
其中，t是平移分量，tx、ty分别为在x、y轴上的平移分量；a11、a12、a21与a22为源图像与变换图像的坐标经过变换时的映射关系；变换矩阵h控制缩放、旋转、剪切效应；
通过使用齐次坐标，上述矩阵形式变为：
再分别计算源图像与变换图像的特征向量集合，分别记为da和db，特征向量集合da和db中的特征向量所对应的关键点记为集合a和b。假设：da＝{d1,d2,...,dn-1,dn}，db＝{d1',d'2...,d'm}，a＝{p1,p2,....,pn}，b＝{p'1,p'2,....,p'm}；其中，di、d'j对应的表示特征向量集合da、db中的特征向量，pi、p'j为di、d'j对应的关键点，下标为关键点的序号，i＝1,2,...,n，j＝1,2,...,m。
步骤2、利用两个cuda内核在gpu上来分别实现关键点的粗匹配与精细匹配。
在cpu上运行算法的每个部分都存在一个并行的实现版本，且这些部分是隔离的，因此设计几个内核来实现完整的拼接。本发明实施例中最耗时的部分是特征匹配，因此将注意力集中在并行特征匹配和相应的优化策略上。
1、粗匹配。
粗匹配时，对关键点集合a中每一个关键点，计算与关键点集合所有关键点匹配关系，计算时带入相应关键点对应的特征向量，从而计算特征向量之间的距离；采用合并排序的并行排序算法，来按照升序对已经计算了距离值使用同步机制进行排序；计算排序后的序列中第一和第二距离的比值，如果比值小于第一阈值，则认为相应关键点为匹配点，否则不是匹配点。假设粗匹配阶段，共获得n个匹配点对。
粗匹配时，首先计算集合a中的一个关键点pk(k＝1,2,3...n)与集合b中所有关键点之间的欧式距离，然后将所获得的欧式距离值按升序排列。计算按照升序排列的第一个值与第二个值的比值，如果比值小于预定的阈值t1,则关键点pk将被识别为匹配点，否则它不是匹配点。
为了实现粗匹配的并行化，粗匹配时，块的数量由特征向量集合da的数量n确定，每个块中线程数量由特征向量集合db的数量m确定，利用上述方式来组织内核，每个块实现关键点集合a中每一个关键点与关键点集合b中所有关键点之间的匹配，匹配时计算相应特征向量之间的距离，并配合第一阈值，进行匹配，所有块的匹配结果存储在全局存储器中，共有n个匹配点对。
为了提高并行实现的性能，优化策略是通过以下方式进行的：(1)代替使用全局内存，将第一阈值放入常量内存来使用它的缓存，一方面缓存的访问速度高于全局内存，另一方面gpu上运行的所有线程都可以；(2)由于共享内存的访问速度为快速寄存器，将特征向量集合da放入共享内存中可以加快数据访问的速度；(3)为了有效实现合并排序，共享内存也用于优化策略。
2、精细匹配
采用随机样本一致性(ransac)算法对匹配对做进一步的精细匹配，由于变换矩阵h只有6个参数，因此精细匹配时，随机选择三个不在一行的匹配点对(即三个非线性的匹配点)，来计算变换矩阵h；再使用变换矩阵h来计算其余的n-3个匹配点对的误差，如果误差小于第二阈值，则相应匹配点对属于内点匹配，即内点对；否则，属于外点匹配，即外点对。
当进行大量的迭代计算(k>200，k为迭代的次数)每次迭代记录内点对和外点对。选择包含最大数量的内点对的迭代，并将这些内点对视为精细的匹配点对。
为了实现精细匹配的并行化，将内核设计为执行ransac算法的每次迭代。因为ransac算法是一种迭代算法，因此只启动一个块，并且线程数是n-3。对于ransac算法来说，并行的主要过程遵循以下原则：选择三个非线性的匹配点对由cpu计算仿射变换时所使用的变换矩阵h；利用n-3个线程来结合第二阈值、变换矩阵h与其余的n-3个匹配点对，来区分匹配点对属于内点匹配或外点匹配，从而筛选出精确的匹配点对，以此提高并行计算算法的效率。
本发明实施例中，将第二阈值和得到的变换矩阵h放入访问速度比全局内存高得多的常量内存中，其余的n-3个匹配放入共享内存，这样做提高了并行执行的性能。这样的内存分配可以实现对并行精细匹配的优化。
本发明实施例上述方案是利用cuda架构在gpu上运行的并行图像特征点匹配的方法，通过多个cuda内核将顺序算法转换为并行版本，使用不同类型的内存进而实现了并行优化算法，将gpu的性能与cpu的性能进行了比较，实现了27.6的加速比。
另一方面，在获得精确的匹配点对后，还利用筛选出精确的匹配点对，结合最小二乘法估计变换矩阵h中的6个参数。
由于之前，在精细匹配的时候，在一定的误差范围内，我们得到了一系列的变换矩阵h的参数，利用最小二乘法可以更好的估计我们需要的变换矩阵h的6个参数。
本发明实施例中，由于精细匹配是一个迭代过程，迭代计算时随机选择很多次(每次三对点)，每次计算生成一个变换矩阵h，这样就得到了一系列的h矩阵。一定误差是由于前面粗匹配时可能产生误匹配点对，在精细匹配的第一步是随机选择的三对点，这可能会选到误匹配点对，故存在误差。
再使用计算出6个参数的变换矩阵h将变换图像变换到参考坐标系中。变换过程通常称为翘曲或映射。可以使用两种基本类型的翘曲(正向映射和反向映射)来校正由于所产生的马赛克中既没有出现孔洞，也没有发生重叠，本发明采用反向映射模式。
考虑到精度和计算复杂度之间的折衷，在实际中，一些映射的像素坐标可能不是非正整数，因此这些原始像素将没有强度值。在这种情况下，还需要进行插值操作来确定原始像素的强度值。同时，插值可以产生更平滑、更精确的亮度值，并进一步减小由变形产生的误差。最常用的插值方法是最近邻插值、双线性插值和双三次插值。考虑精度和计算复杂度之间的关系，采用双线性插值算法。
本发明实施例中，将变换图像变换到参考坐标系时，所涉及的变形和插值操作通过一个内核来实现并行化处理；内核的线程数量等于重叠图像中包含的像素(像素指的是图像的灰度值)，这样可以保证每个像素都由一个线程处理。必须注意，每个块中的最佳线程数应该是32的倍数。在开始时，且将计算到的变换矩阵h的逆矩阵h-1存入常量内存，并被所有线程使用。(这里的重叠中的意思是源图像与变换图像像素的重叠区域)
在开始时，计算变换矩阵h的逆矩阵h-1并将其放入常量内存中，因为变换矩阵h的逆矩阵h-1将被所有线程频繁使用。
由于插值操作为局部操作，而纹理存储器是专门为本地存取模式设计的，因此将相关的变换图像存入纹理存储器中；为了进一步提高性能，在插值操作时，如果两个坐标值的小数部分小于0.2，则将属于整数坐标的强度值分配给校正像素，否则采用双线性插值计算新的强度值。
本领域技术人员可以理解，变形与插值操作的具体方案可以参见现有技术，也即，本发明实施例关于变形与插值操作中未详细介绍的流程都可以参照现有技术。
通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例可以通过软件实现，也可以借助软件加必要的通用硬件平台的方式来实现。基于这样的理解，上述实施例的技术方案可以以软件产品的形式体现出来，该软件产品可以存储在一个非易失性存储介质(可以是cd-rom，u盘，移动硬盘等)中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述的方法。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。