1、利用SIFT算法实现并优化基于GPU的通用计算摘要：现代图形硬件中加入了可编程序组件后，利用图形处理器（GPUs）进行通用计算就变得越来越有趣，尤其是在利用了并行缓冲处理机制后，这一功能得到了更广泛的应用。在本文中，我们介绍了利用现代图形硬件进行实时跟踪和识别特征点的方法。重点是对于一些特征向量的研究，这些特征向量是从不同阶段的输入图像中提取出来的。尺度不变特征变换(SIFT)方法对图像的旋转、尺度、和光照条件的改变具有较高稳定性，因而被广泛使用于对图像特征向量的处理。我们提供了不同阶段中利用图像处理器对特征向量进行处理的结果，并且将它与基于SIFT算法的中央处理器的处理结果进行比较。这个方法
2、可以在Geforce6图形版系列中很好的实现，并且能够充分利用新的硬件特点，比如碎片处理器中的动态分支和多重着色目标(MRT)。我们可以在图形处理器上实现上述实时帧率特征跟踪方法，与此同时中央处理器还可以执行其它任务。关键字：图形处理器、尺度不变特征变换、特征提取、跟踪。1 引言随着现代图形硬件中可编程元件的加入，比如顶点和片段处理器，开发人员开始使用图形处理器(GPUs)进行通用计算，而不再仅仅利用它制作精美的图片和高端游戏引擎。这项研究是从GPGPU (GPU实现通用计算) 中发展起来的。在研究利用GPU实现通用计算的过程中，当务之急就是变换正在研究的算法，使它能更充分地利用现代图形硬件上
3、的并行缓冲处理模型。大量的算法和应用已经在GPGPU 上得以实现，例如，kd-树搜索FS05，分类算法加速GHLM05和数据搜索GLW*04。针对GPU 的架构进行变换算法的概述在Har05和OLG*05中给出。本文中，我们会介绍特征提取算法如何利用现代图形硬件实现功能，以及优化各阶段算法达到加速目的。确切地说，就是在内含256MB 音频存储器的NVIDIA QuadroFX 3400 GPU上利用 SIFT算法 Low04a实现新的功能，例如，动态分支、片段上多个渲染目标(MRTs)以及顶点处理器。在了解了相关工作之后，我们先简单介绍了一下算法，接着是GPU执行步骤以及如何在标准CPU 上实
4、现算法全面加速。2 相关工作计算机视觉的一个主要研究问题即是自然图像稳定特征点的提取。这些特征点主要用于对应匹配，以便找到已知对象，并且获得它们在其它相应图像中的效果、位置、大小或者旋转度的信息BL02。在Fun05中给出了一种特征跟踪的方法，但是该方法中提取的特征点不具有尺度不变性，并且该方法主要应用在多并行显卡上。DavidLowLow04提出了一种从图像中提取高度稳定特征向量的方法，该方法叫做尺度不变特征变换(SIFT) 。 SIFT特征具有旋转不变性，尺度和照明度/对比度不变性，因此SIFT方法能很好的应用在场景模型，识别，跟踪GL04以及全景制作中。3 SIFT概述SIFT方法介绍了
5、获取相关特征点的各个步骤。分析调整各步骤以便最大化地使用图形处理器(GPU) 的并行处理机制，并且占用较少的 CPU进程。单个SIFT特征的提取步骤是：1.寻找兴趣点，兴趣点由高斯差异(DOG)尺度空间金字塔模型识别，在尺度空间中图像中由极值表示。2.对第一步中得到的点进行进一步的过滤和筛选，以得到具有高匹配度的稳定点。3.通过确定特征方向，对每个点的主方向进行指定。4.由第三步中所确定的方向进行特征向量的计算，确保其具有旋转不变性。5.整个过程中，只有一部分元素从一个步骤的开始传递到下一步骤，并由此展开。为了实现尺度不变性，有必要创建一种表示图像频率的模型，如WitkinWit83所介绍的一
6、样，该模型利用尺度空间金字塔可以实现。金字塔里面的每张图像都对应不同的图像频率。通过搜索尺度空间金字塔中所有图像，获取潜在的特征点，并使其具有尺度不变性。尺度空间金字塔由灰度版本的原始图像和尺度不断增大的高斯卷积核反复卷积而组成。因此大量清晰度不同的图像由此生成，就像在图1左上方所显示的四个640*480的图像堆叠一样。接下来，这堆叠的图像中，最模糊的图像缩减两倍尺度。随后与一套同样的高斯卷积核进行卷积，就像之前一样生成第二组4张320*240图像。反复进行此过程直到达到指定的大小，也就是图1右上方例子中的80*60大小的那张图。图1. 高斯金字塔( 顶部)和高斯的不同( 底部)。原始图像经过
7、持续不断的运算生成大小和清晰度不同的图像版本，这个金字塔即由这些不同版本图像组成。为了计算一个图像的单一频率，在金字塔中相同大小的邻近图像或图层必须进行相减运算，生成高斯差异模型，如图1下方所示。最终，通过筛选这个DOG金字塔，从而找到极值。将每一个采样点和它所有的相邻点进行比较，中间的检测点和它同尺度的8个相邻点进行比较，并和上下相邻尺度对应的各9个点进行比较，观察其是否比它的图像域和尺度域的相邻点大或者小，从而确定该采样点是否是极值点。在下一步中，对上述检测到的可能的极值点进行进一步筛选，以去除那些不能够进行匹配的特征点。这里，需要考虑的应删除点主要有两类。首先，一些非极值点之所以被检测为
8、极值点是因为输入图像噪声的影响，其次，该点可能位于图像的边缘。对于可能的特征点及其邻近点是通过引进亮度门限值进行消除的。然而，如果该点的值超过门限值，这个点将被当作候选的特征点进入下一步的处理及运算。边缘点需要排除，因为它们不适合跟踪和对应匹配。对于边缘点检测，在位置(x，y)的候选点D(x，y),其表面曲率将用二阶矩阵Hessian 矩阵H进行分析：就像Harris和StephensHS88所介绍的那样，D(x,y) 的主曲率和H是成正比的。因为我们只关注该点是否是边缘点，对于两个特征值e1和 e2，只有比例在 e1e2时才符合要求。设r = e1/e2，然后：和 （2）因此，非边缘点的判别
9、标准可由下式表示：在Lowe中，为了排除边缘点，最好的结果取 r=10。在第三步中，为特征点分配主方向，以便使其具有旋转不变性。因此先将每个特征点邻域像素的梯度方向转换为极坐标：然后，从相位值(x,y) 和加权极值m(x,y) 中，构建出直方图，以获得特征点周围主要的梯度方向。图2. 由关键点邻域梯度信息生成特征向量最后，由梯度方向确定特征向量。在图2左边，显示出了图像梯度。对于每个4*4的区域，其主要方向均在新的坐标系中显示，这也意味着新的纹理坐标也要插在相应的中间位置。在新的坐标系统中，得到特定区域的16个梯度。原坐标中的16个梯度方向投影到纹理坐标中的8个方向，以此得到特征点的描述子，如
10、图2右边所示。按照这种方法投影到所有4*4区域，因而就生成了一个128维的特征向量。最后，对特征向量进行长度归一化处理，可以进一步去除由光照变化造成的影响。4 GPU的实现由于在GPU上应用SIFT 算法，那么原始的CPU 算法就需要配合该算法进入到图形管道Zel05中去并且要充分利用 GPU并行处理能力。因此，研究重点就是如何重组 SIFT各个步骤，以适应GPU 纹理格式。在一个系统中即可完成此实验，该系统应采用英特尔Xeon 3.2兆赫的 CPU，2GB的RAM和NVIDIA QuadroFX 3400 GPU，内带256MB音频RAM和PCI Express x16标准总线。4.1 Do
11、G 金字塔的生成SIFT特征通常取自图像的灰度层，而GPU纹理缓存器则设计为一个三色加阿尔法通道。GPU不仅具有基于片段处理器对每个像素进行并行处理的能力，而且还可以作为一个载体进行GPU 四色值的并行计算。只有灰度图像才可以参与与另一张图像进行卷积的计算，这会浪费GPU75%的处理能力。为了充分利用GPU的矢量处理能力，对输入图像的灰度级要进行改进。这里，我们将灰度图像数据重新排列到一个四通道RGBA图像中去，如图3所示。图3. 纹理打包到RGBA16 的 GPU表中在RGBA图像中的一种颜色值，代表灰度级图像中的2*2像素，因此，RGBA图像将图像区域减小为原来的4四分之一。这样卷积就不会
12、浪费GPU的计算能力。在一个卷积中，对打包数据的处理是很简单的，因为大部分是线性操作，例如对像素进行相加或者相乘。在具体操作的情况下，对像素的处理也取决于相邻的像素，因为所有相邻参考点都需要被重新定向，所以对于打包的数据来说，它们对算法的适应也变得复杂起来。对输入图像的重组会造成一些计算开销，但这些是比较小的，因为在整个尺度空间和DOG金字塔生成的全过程中，数据一致保持为压缩格式。数据打包是采用一个简单的片段着色来实现，它是将一块2*2相邻像素块排列到一个RGBA像素中去。高斯卷积是直接应用到 RGBA图像打包格式中，如图4中9个高斯卷积核函数。这里，高斯核被分成偶数和奇数值，从而进行独立的半
13、卷积，最后相加得到最终结果。在一个计算中，每个高斯核的像素都与四色组件相乘，因此只需要一个访问进程。这个奇数和偶数高斯核像素的计算，应用相同的片段着色器，因此这个相同纹理访问可以同时用于两步骤。图4. 在RGBA16 GPU 图表中高斯卷积，(a)奇数和(b)偶数的样例。利用不同大小的高斯卷积核进行垂直和水平的筛选已经得到成功应用，不同清晰度的图像也可通过做差生成DOG金字塔，如上所述。通过这个技术，我们可以将一张彩色图像转换成一张灰度级图像，并且按描述方式在一个渲染通道上打包这些像素。4.2 关键点筛选和定向如前所述，为了特征点的检测，动态分支必须保证整个选择处理过程都在GPU中进行。因此，
14、对潜在特征点按照亮度差异门限值进行重新排列，按照新的判别标准进行筛选后，得到的新的潜在特征点不会超过原潜在特征点的50% 。然后进行全局极值搜索，首先在同一个缓存器内的一个点与它的8个相邻点进行比较，只有0.6%的潜在特征点会与DOG金字塔中临近缓存区的9个相邻点进行比较。潜在的特征点如图5(b)所示。随后，排除了噪声和边缘点后，留下稳定的特征点，如图5(c) 所示。图5. 特征的提取和过滤在筛选和定位潜在特征点之后，就应计算相应的特征向量。为了计算梯度方向和大小，就需要使用MRT 功能。对于大小和方向这两个值，只需考虑四个直接相邻的像素，它们可以相互参照使计算相应简单一些。这些像素访问和操作
15、如图6所示。图6. 红色通道的梯度级(top)和方向计算 (底部)红色通道在这里，对中央纹素的“RGBA”模型和它的四个原始像素进行处理，这四个原始像素已经被打包成四个颜色组件。对中间组件“r”的大小和方向的计算需要四个像纯色一样突出的相邻组件，如图6中显示。每个中间组件都由中央纹素中的两个组件和相邻纹素的两个组件构成。这个计算大小和方向的请求处理同样如图6所示。这两个计算都需要用到相同的输入数据。因为输入图像数据已经打包，计算大小和方向时用到的两个颜色组件不会独立作用。相反，由于两个计算可以同时进行，并且可以写入两个独立的渲染目标中，因此，MRT的使用大大加快了处理进程。所以，因图形程序接口
16、中内容切换所造成的时间损失就可以避免，同时，由于大小和方向的计算使用了同样的中间过渡计算方式(即水平和垂直减法) ，所以当两个动作同时执行时只会占用一个纹素通道。4.3 特征描述子的生成现在两个渲染目标都用来生成特征点周围4*4区域的加权直方图，如图2理论部分所示。每个区域都包括八个方向，加一起生成了一个128维的输出向量。因为每个单一的输入点（或极值点）都会生成一个128维的输出向量，所以现在的操作不同于之前的实现方式。这个操作不能实时执行，即使使用MRTS技术也做不到。由于针对直方图本身的计算不可能分解，因此每个区域都会在一个片段着色调用模块中进行处理。为了实现这个计算，我们必须找到一个合
17、理的结构，以便每个片段可以访问不同的数据。在其它计算中用到的简单的矩形结构在这里是不能满足要求的，因为通过差值算法会将四个角落的属性插入到整个区域，而我们需要的是每个点都具有独立的属性。这种独立属性可以用顶点栅格来表示，顶点栅格是由几何点通过glVertex2f()函数和片段处理生成。为了便于处理，大小和方向值从两个渲染目标中重新排列到一个纹理中去，以便对它们进行纹理方向的预计算处理，如图7所示。图7.渐变映射解压成一个纹理这样，大小和方向打包的数据值经过解压缩和交错处理后，就可以实现一个输出值包含一个大小值和与之关联的方向值。这种形式下，一个纹理就包含了两个值，对纹理的进一步处理就不需要改变
18、格式了。为了生成最终特征，我们需要对每个极值点的4*4区域创建一个梯度直方图。在没有进行复杂计算的情况下，针对直方图本身的计算不会分解开来，所以每个区域必须在一个片段着色调用器中进行处理。为了使直方图计算在一个着色周期内实现，8个输出向量必须同时进行计算。在先进的显卡上使用MRTs技术，我们就可以实现上述要求。在图8（a）中，我们可以看到特征生成的数据结构，图8（b）中我们可以看到一个示例的特征向量，它包含16个定点，并且被映射到了两个渲染目标中。图8. 特征生成的渲染目标。（a）两个渲染目标中的帧缓存区（b）图像中的特征点位置和对应的纹理坐标每一个特征向量的顶点都会在CPU中具有相应的属性。
19、这些属性包括对应的纹理坐标，梯度区域的大小和方向值。相应的计算可以独立的进行，所有必要的参数都会在特征向量中进行编码。因此，在图8（a）中，每个渲染目标都包含由一半的特征向量设置而成的一个完整数据。这些SIFT特征向量可以在不同的图像中进行相应的特征匹配，例如对图像序列进行跟踪。为此，我们需要计算两个长度为N的特征向量 V1和V2之间的欧式距离D，如式（5）：在我们的SIFT 方法实现中，每个向量都有 N=128个元素。对应的（5）式可以进行并行计算。在测试中，CPU 实现通用计算和GPU 实现通用计算会以处理时间作为标准进行比较，在这个测试中，我们会采用一系列的特征向量K ，其中 K500到
20、K3000逐渐变化。每个特征向量都会进行相互比较，因此会进行K 2次比较。当处理 3000*3000次比较时，CPU需要用时13s ，而GPU 只需0.5s 。4.4 结果为了实现高效的GPU处理而进行SIFT算法步骤优化后，整个算法都与原始的CPU处理进行了测试和比较，同时也和手动SSE（SIMD 流指令扩展）的优化版本进行了比较，结果如图9所示。图9.基于SIFT 描述的GPU 实现与标准CPU处理的比较结果这里，手动的SSE优化版本需要 0.312s而原始的CPU 实现要0.406s。从对比中发现，SIFT算法可以大大加快利用大规模并行处理的GPU的运算速度，因而实现了处理时间只有0.058s 。5 结论在本文中，我们已经介绍，SIFT算法可以大大加快利用GPU 的并行处理速度。图像文件输入到GPU 的RGBA格式的纹理缓冲后，重新改变其亮度值，所有的后续步骤也都会适应这种纹理格式，并使用新的GPU技术，即检测相关的特征点的动态分支技术、对梯度方向及大小进行并行计算的MRTs技术。因此，SIFT算法可以应用于每秒 20帧640*480像素的图像序列中。今后的工作将主要集中在使用SIFT 特征实现实时应用，例如图像序列的三维场景重建的校准估计。
