SIFT和surf

来源：互联网发布：佳能wifi软件app 编辑：程序博客网时间：2024/05/19 09:40

高斯拉普拉斯（Laplace of Gaussian）

kezunhai@gmail.com

http://blog.csdn.net/kezunhai

Laplace算子作为一种优秀的边缘检测算子，在边缘检测中得到了广泛的应用。该方法通过对图像求图像的二阶倒数的零交叉点来实现边缘的检测，公式表示如下：

由于Laplace算子是通过对图像进行微分操作实现边缘检测的，所以对离散点和噪声比较敏感。于是，首先对图像进行高斯卷积滤波进行降噪处理，再采用Laplace算子进行边缘检测，就可以提高算子对噪声和离散点的鲁棒性，如此，拉普拉斯高斯算子Log（Laplace of Gaussian）就诞生了。

高斯卷积（Gaussian convolution ），高斯函数的表达式如下：

原图像与高斯卷积的表达式如下：

因为：

所以Log可以通过先对高斯函数进行偏导操作，然后进行卷积求解，公式表示如下：

2D高斯拉普拉斯算子可以通过任何一个方形核进行逼近，只要保证该核的所有元素的和或均值为0，如下一个5×5的核进行逼近：

高斯拉普拉斯边缘检测算法的步骤：

1）对原图像进行Log卷积。

2）检测图像中的过零点（ Zero Crossings，也即从负到正或从正到负）。

3）对过零点进行阈值化。

尺度不变特征变换匹配算法详解
Scale Invariant Feature Transform(SIFT)
Just For Fun

zdd zddmail@gmail.com

对于初学者，从David G.Lowe的论文到实现，有许多鸿沟，本文帮你跨越。

1、SIFT综述

尺度不变特征转换(Scale-invariant feature transform或SIFT)是一种电脑视觉的算法用来侦测与描述影像中的局部性特征，它在空间尺度中寻找极值点，并提取出其位置、尺度、旋转不变量，此算法由 David Lowe在1999年所发表，2004年完善总结。

其应用范围包含物体辨识、机器人地图感知与导航、影像缝合、3D模型建立、手势辨识、影像追踪和动作比对。

此算法有其专利，专利拥有者为英属哥伦比亚大学。

局部影像特征的描述与侦测可以帮助辨识物体，SIFT 特征是基于物体上的一些局部外观的兴趣点而与影像的大小和旋转无关。对于光线、噪声、些微视角改变的容忍度也相当高。基于这些特性，它们是高度显著而且相对容易撷取，在母数庞大的特征数据库中，很容易辨识物体而且鲜有误认。使用 SIFT特征描述对于部分物体遮蔽的侦测率也相当高，甚至只需要3个以上的SIFT物体特征就足以计算出位置与方位。在现今的电脑硬件速度下和小型的特征数据库条件下，辨识速度可接近即时运算。SIFT特征的信息量大，适合在海量数据库中快速准确匹配。

SIFT算法的特点有：

1. SIFT特征是图像的局部特征，其对旋转、尺度缩放、亮度变化保持不变性，对视角变化、仿射变换、噪声也保持一定程度的稳定性；

2. 独特性（Distinctiveness）好，信息量丰富，适用于在海量特征数据库中进行快速、准确的匹配；

3. 多量性，即使少数的几个物体也可以产生大量的SIFT特征向量；

4. 高速性，经优化的SIFT匹配算法甚至可以达到实时的要求；

5. 可扩展性，可以很方便的与其他形式的特征向量进行联合。

SIFT算法可以解决的问题：

目标的自身状态、场景所处的环境和成像器材的成像特性等因素影响图像配准/目标识别跟踪的性能。而SIFT算法在一定程度上可解决：

1. 目标的旋转、缩放、平移（RST）

2. 图像仿射/投影变换（视点viewpoint）

3. 光照影响（illumination）

4. 目标遮挡（occlusion）

5. 杂物场景（clutter）

6. 噪声

SIFT算法的实质是在不同的尺度空间上查找关键点(特征点)，并计算出关键点的方向。SIFT所查找到的关键点是一些十分突出，不会因光照，仿射变换和噪音等因素而变化的点，如角点、边缘点、暗区的亮点及亮区的暗点等。

Lowe将SIFT算法分解为如下四步：

1. 尺度空间极值检测：搜索所有尺度上的图像位置。通过高斯微分函数来识别潜在的对于尺度和旋转不变的兴趣点。

2. 关键点定位：在每个候选的位置上，通过一个拟合精细的模型来确定位置和尺度。关键点的选择依据于它们的稳定程度。

3. 方向确定：基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向。所有后面的对图像数据的操作都相对于关键点的方向、尺度和位置进行变换，从而提供对于这些变换的不变性。

4. 关键点描述：在每个关键点周围的邻域内，在选定的尺度上测量图像局部的梯度。这些梯度被变换成一种表示，这种表示允许比较大的局部形状的变形和光照变化。

本文沿着Lowe的步骤，参考Rob Hess及Andrea Vedaldi源码，详解SIFT算法的实现过程。

2、高斯模糊

SIFT算法是在不同的尺度空间上查找关键点，而尺度空间的获取需要使用高斯模糊来实现，Lindeberg等人已证明高斯卷积核是实现尺度变换的唯一变换核，并且是唯一的线性核。本节先介绍高斯模糊算法。

2.1二维高斯函数

高斯模糊是一种图像滤波器，它使用正态分布(高斯函数)计算模糊模板，并使用该模板与原图像做卷积运算，达到模糊图像的目的。

N维空间正态分布方程为：

（1-1）

其中，是正态分布的标准差，值越大，图像越模糊(平滑)。r为模糊半径，模糊半径是指模板元素到模板中心的距离。如二维模板大小为m*n，则模板上的元素(x,y)对应的高斯计算公式为：

（1-2）

在二维空间中，这个公式生成的曲面的等高线是从中心开始呈正态分布的同心圆，如图2.1所示。分布不为零的像素组成的卷积矩阵与原始图像做变换。每个像素的值都是周围相邻像素值的加权平均。原始像素的值有最大的高斯分布值，所以有最大的权重，相邻像素随着距离原始像素越来越远，其权重也越来越小。这样进行模糊处理比其它的均衡模糊滤波器更高地保留了边缘效果。

理论上来讲，图像中每点的分布都不为零，这也就是说每个像素的计算都需要包含整幅图像。在实际应用中，在计算高斯函数的离散近似时，在大概3σ距离之外的像素都可以看作不起作用，这些像素的计算也就可以忽略。通常，图像处理程序只需要计算的矩阵就可以保证相关像素影响。

2.2 图像的二维高斯模糊

根据σ的值，计算出高斯模板矩阵的大小()，使用公式(1-2)计算高斯模板矩阵的值，与原图像做卷积，即可获得原图像的平滑(高斯模糊)图像。为了确保模板矩阵中的元素在[0,1]之间，需将模板矩阵归一化。5*5的高斯模板如表2.1所示。

下图是5*5的高斯模板卷积计算示意图。高斯模板是中心对称的。

2.3分离高斯模糊

如图2.3所示，使用二维的高斯模板达到了模糊图像的目的，但是会因模板矩阵的关系而造成边缘图像缺失(2.3 b,c)，越大，缺失像素越多,丢弃模板会造成黑边(2.3 d)。更重要的是当变大时，高斯模板(高斯核)和卷积运算量将大幅度提高。根据高斯函数的可分离性，可对二维高斯模糊函数进行改进。

高斯函数的可分离性是指使用二维矩阵变换得到的效果也可以通过在水平方向进行一维高斯矩阵变换加上竖直方向的一维高斯矩阵变换得到。从计算的角度来看，这是一项有用的特性，因为这样只需要次计算，而二维不可分的矩阵则需要次计算，其中，m,n为高斯矩阵的维数，M,N为二维图像的维数。

另外，两次一维的高斯卷积将消除二维高斯矩阵所产生的边缘。(关于消除边缘的论述如下图2.4所示，对用模板矩阵超出边界的部分——虚线框，将不做卷积计算。如图2.4中x方向的第一个模板1*5，将退化成1*3的模板，只在图像之内的部分做卷积。)

附录1是用opencv2.2实现的二维高斯模糊和分离高斯模糊。表2.2为上述两种方法和opencv2.3开源库实现的高斯模糊程序的比较。

3、尺度空间极值检测

尺度空间使用高斯金字塔表示。Tony Lindeberg指出尺度规范化的LoG(Laplacion of Gaussian)算子具有真正的尺度不变性，Lowe使用高斯差分金字塔近似LoG算子，在尺度空间检测稳定的关键点。

3.1 尺度空间理论

尺度空间(scale space)思想最早是由Iijima于1962年提出的，后经witkin和Koenderink等人的推广逐渐得到关注，在计算机视觉邻域使用广泛。

尺度空间理论的基本思想是：在图像信息处理模型中引入一个被视为尺度的参数，通过连续变化尺度参数获得多尺度下的尺度空间表示序列，对这些序列进行尺度空间主轮廓的提取，并以该主轮廓作为一种特征向量，实现边缘、角点检测和不同分辨率上的特征提取等。

尺度空间方法将传统的单尺度图像信息处理技术纳入尺度不断变化的动态分析框架中，更容易获取图像的本质特征。尺度空间中各尺度图像的模糊程度逐渐变大，能够模拟人在距离目标由近到远时目标在视网膜上的形成过程。

尺度空间满足视觉不变性。该不变性的视觉解释如下：当我们用眼睛观察物体时，一方面当物体所处背景的光照条件变化时，视网膜感知图像的亮度水平和对比度是不同的，因此要求尺度空间算子对图像的分析不受图像的灰度水平和对比度变化的影响，即满足灰度不变性和对比度不变性。另一方面，相对于某一固定坐标系，当观察者和物体之间的相对位置变化时，视网膜所感知的图像的位置、大小、角度和形状是不同的，因此要求尺度空间算子对图像的分析和图像的位置、大小、角度以及仿射变换无关，即满足平移不变性、尺度不变性、欧几里德不变性以及仿射不变性。

3.2 尺度空间的表示

一个图像的尺度空间，定义为一个变化尺度的高斯函数与原图像的卷积。

(3-1)

其中，*表示卷积运算，

(3-2)

与公式(1-2)相同，m，n表示高斯模板的维度(由确定)。(x, y)代表图像的像素位置。是尺度空间因子，值越小表示图像被平滑的越少，相应的尺度也就越小。大尺度对应于图像的概貌特征，小尺度对应于图像的细节特征。

3.3 高斯金字塔的构建

尺度空间在实现时使用高斯金字塔表示，高斯金字塔的构建分为两部分：

1. 对图像做不同尺度的高斯模糊；

2. 对图像做降采样(隔点采样)。

图像的金字塔模型是指，将原始图像不断降阶采样，得到一系列大小不一的图像，由大到小，从下到上构成的塔状模型。原图像为金子塔的第一层，每次降采样所得到的新图像为金字塔的一层(每层一张图像)，每个金字塔共n层。金字塔的层数根据图像的原始大小和塔顶图像的大小共同决定，其计算公式如下：

(3-3)

其中M，N为原图像的大小,t为塔顶图像的最小维数的对数值。如，对于大小为512*512的图像，金字塔上各层图像的大小如表3.1所示，当塔顶图像为4*4时，n=7，当塔顶图像为2*2时，n=8。

为了让尺度体现其连续性，高斯金字塔在简单降采样的基础上加上了高斯滤波。如图3.1所示，将图像金字塔每层的一张图像使用不同参数做高斯模糊，使得金字塔的每层含有多张高斯模糊图像，将金字塔每层多张图像合称为一组(Octave)，金字塔每层只有一组图像，组数和金字塔层数相等，使用公式(3-3)计算，每组含有多张(也叫层Interval)图像。另外，降采样时，高斯金字塔上一组图像的初始图像(底层图像)是由前一组图像的倒数第三张图像隔点采样得到的。

注：由于组内的多张图像按层次叠放，因此组内的多张图像也称做多层，为避免与金字塔层的概念混淆，本文以下内容中，若不特别说明是金字塔层数，层一般指组内各层图像。

注：如3.4节所示，为了在每组中检测S个尺度的极值点，则DOG金字塔每组需S+2层图像，而DOG金字塔由高斯金字塔相邻两层相减得到，则高斯金字塔每组需S+3层图像，实际计算时S在3到5之间。取S=3时，假定高斯金字塔存储索引如下：

第0组(即第-1组)： 0 1 2 3 4 5

第1组： 6 7 8 9 10 11

第2组：？

则第2组第一张图片根据第一组中索引为9的图片降采样得到，其它类似。

3.4 高斯差分金字塔

2002年Mikolajczyk在详细的实验比较中发现尺度归一化的高斯拉普拉斯函数的极大值和极小值同其它的特征提取函数，例如：梯度，Hessian或Harris角特征比较，能够产生最稳定的图像特征。

而Lindeberg早在1994年就发现高斯差分函数（Difference of Gaussian ，简称DOG算子）与尺度归一化的高斯拉普拉斯函数非常近似。其中和的关系可以从如下公式推导得到：

利用差分近似代替微分，则有：

因此有

其中k-1是个常数，并不影响极值点位置的求取。

如图3.2所示，红色曲线表示的是高斯差分算子，而蓝色曲线表示的是高斯拉普拉斯算子。Lowe使用更高效的高斯差分算子代替拉普拉斯算子进行极值检测，如下：

(3-4)

在实际计算时，使用高斯金字塔每组中相邻上下两层图像相减，得到高斯差分图像，如图3.3所示，进行极值检测。

3.5 空间极值点检测(关键点的初步探查)

关键点是由DOG空间的局部极值点组成的，关键点的初步探查是通过同一组内各DoG相邻两层图像之间比较完成的。为了寻找DoG函数的极值点，每一个像素点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或者小。如图3.4所示，中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较，以确保在尺度空间和二维图像空间都检测到极值点。

由于要在相邻尺度进行比较，如图3.3右侧每组含4层的高斯差分金子塔，只能在中间两层中进行两个尺度的极值点检测，其它尺度则只能在不同组中进行。为了在每组中检测S个尺度的极值点，则DOG金字塔每组需S+2层图像，而DOG金字塔由高斯金字塔相邻两层相减得到，则高斯金字塔每组需S+3层图像，实际计算时S在3到5之间。

当然这样产生的极值点并不全都是稳定的特征点，因为某些极值点响应较弱，而且DOG算子会产生较强的边缘响应。

3.6 构建尺度空间需确定的参数

—尺度空间坐标

O—组(octave)数

S— 组内层数

在上述尺度空间中，O和S，的关系如下：

(3-5)

其中是基准层尺度，o为组octave的索引，s为组内层的索引。关键点的尺度坐标就是按关键点所在的组和组内的层，利用公式(3-5)计算而来。

在最开始建立高斯金字塔时，要预先模糊输入图像来作为第0个组的第0层的图像，这时相当于丢弃了最高的空域的采样率。因此通常的做法是先将图像的尺度扩大一倍来生成第-1组。我们假定初始的输入图像为了抗击混淆现象，已经对其进行的高斯模糊，如果输入图像的尺寸用双线性插值扩大一倍，那么相当于。

取式(3-4)中的k为组内总层数的倒数，即

(3-6)

在构建高斯金字塔时，组内每层的尺度坐标按如下公式计算：

(3-7)

其中初始尺度，lowe取，s为组内的层索引，不同组相同层的组内尺度坐标相同。组内下一层图像是由前一层图像按进行高斯模糊所得。式(3-7)用于一次生成组内不同尺度的高斯图像，而在计算组内某一层图像的尺度时，直接使用如下公式进行计算：

(3-8)

该组内尺度在方向分配和特征描述时确定采样窗口的大小。

由上，式(3-4)可记为

(3-9)

图3.5为构建DOG金字塔的示意图，原图采用128*128的jobs图像，扩大一倍后构建金字塔。

4、关键点定位

以上方法检测到的极值点是离散空间的极值点，以下通过拟合三维二次函数来精确确定关键点的位置和尺度，同时去除低对比度的关键点和不稳定的边缘响应点(因为DoG算子会产生较强的边缘响应)，以增强匹配稳定性、提高抗噪声能力。

4.1关键点的精确定位

离散空间的极值点并不是真正的极值点，图4.1显示了二维函数离散空间得到的极值点与连续空间极值点的差别。利用已知的离散空间点插值得到的连续空间极值点的方法叫做子像素插值（Sub-pixel Interpolation）。

为了提高关键点的稳定性，需要对尺度空间DoG函数进行曲线拟合。利用DoG函数在尺度空间的Taylor展开式(拟合函数)为：

(4-1)

其中，。求导并让方程等于零，可以得到极值点的偏移量为：

(4-2)

对应极值点，方程的值为：

(4-3)

其中,代表相对插值中心的偏移量，当它在任一维度上的偏移量大于0.5时（即x或y或），意味着插值中心已经偏移到它的邻近点上，所以必须改变当前关键点的位置。同时在新的位置上反复插值直到收敛；也有可能超出所设定的迭代次数或者超出图像边界的范围，此时这样的点应该删除，在Lowe中进行了5次迭代。另外，过小的点易受噪声的干扰而变得不稳定，所以将小于某个经验值(Lowe论文中使用0.03，Rob Hess等人实现时使用0.04/S)的极值点删除。同时，在此过程中获取特征点的精确位置(原位置加上拟合的偏移量)以及尺度()。

4.2消除边缘响应

一个定义不好的高斯差分算子的极值在横跨边缘的地方有较大的主曲率，而在垂直边缘的方向有较小的主曲率。

DOG算子会产生较强的边缘响应，需要剔除不稳定的边缘响应点。获取特征点处的Hessian矩阵，主曲率通过一个2x2 的Hessian矩阵H求出：

(4-4)

H的特征值α和β代表x和y方向的梯度，

(4-5)

表示矩阵H对角线元素之和，表示矩阵H的行列式。假设是α较大的特征值，而是β较小的特征值，令，则

(4-6)

导数由采样点相邻差估计得到,在下一节中说明。

D的主曲率和H的特征值成正比，令为α最大特征值，β为最小的特征值，则公式的值在两个特征值相等时最小，随着的增大而增大。值越大，说明两个特征值的比值越大，即在某一个方向的梯度值越大，而在另一个方向的梯度值越小，而边缘恰恰就是这种情况。所以为了剔除边缘响应点，需要让该比值小于一定的阈值，因此，为了检测主曲率是否在某域值r下，只需检测

(4-7)

式(4-7)成立时将关键点保留，反之剔除。

在Lowe的文章中，取r＝10。图4.2右侧为消除边缘响应后的关键点分布图。

4.3有限差分法求导

有限差分法以变量离散取值后对应的函数值来近似微分方程中独立变量的连续取值。在有限差分方法中，我们放弃了微分方程中独立变量可以取连续值的特征，而关注独立变量离散取值后对应的函数值。但是从原则上说，这种方法仍然可以达到任意满意的计算精度。因为方程的连续数值解可以通过减小独立变量离散取值的间格，或者通过离散点上的函数值插值计算来近似得到。这种方法是随着计算机的诞生和应用而发展起来的。其计算格式和程序的设计都比较直观和简单，因而，它在计算数学中使用广泛。

有限差分法的具体操作分为两个部分：

1. 用差分代替微分方程中的微分，将连续变化的变量离散化，从而得到差分方程组的数学形式；

2. 求解差分方程组。

一个函数在x点上的一阶和二阶微商，可以近似地用它所临近的两点上的函数值的差分来表示。如对一个单变量函数f(x)，x为定义在区间[a,b]上的连续变量，以步长将区间[a,b]离散化，我们会得到一系列节点，

然后求出f(x)在这些点上的近似值。显然步长h越小，近似解的精度就越好。与节点相邻的节点有和，所以在节点处可构造如下形式的差值：

节点的一阶向前差分

节点的一阶向后差分

节点的一阶中心差分

本文使用中心差分法利用泰勒展开式求解第四节所使用的导数，现做如下推导。

函数f(x)在处的泰勒展开式为：

(4-8)

则，

(4-9)

(4-10)

忽略h平方之后的项，联立式(4-9)，(4-10)解方程组得：

(4-11)

(4-12)

二元函数的泰勒展开式如下：

将展开后忽略次要项联立解方程得二维混合偏导如下：

(4-13)

综上，推导了4.1,4.2遇到的所有导数计算。同理，利用多元泰勒展开式，可得任意偏导的近似差分表示。

在图像处理中，取h=1，在图4.2所示的图像中，将像素0的基本中点导数公式整理如下：

4.4 三阶矩阵求逆公式

高阶矩阵的求逆算法主要有归一法和消元法两种，现将三阶矩阵求逆公式总结如下：

若矩阵

可逆，即时，

(4-14)

5、关键点方向分配

为了使描述符具有旋转不变性，需要利用图像的局部特征为给每一个关键点分配一个基准方向。使用图像梯度的方法求取局部结构的稳定方向。对于在DOG金字塔中检测出的关键点点，采集其所在高斯金字塔图像3σ邻域窗口内像素的梯度和方向分布特征。梯度的模值和方向如下：

(5-1)

L为关键点所在的尺度空间值，按Lowe的建议，梯度的模值m(x,y)按的高斯分布加成，按尺度采样的3σ原则，邻域窗口半径为。

在完成关键点的梯度计算后，使用直方图统计邻域内像素的梯度和方向。梯度直方图将0~360度的方向范围分为36个柱(bins)，其中每柱10度。如图5.1所示，直方图的峰值方向代表了关键点的主方向，(为简化，图中只画了八个方向的直方图)。

方向直方图的峰值则代表了该特征点处邻域梯度的方向，以直方图中最大值作为该关键点的主方向。为了增强匹配的鲁棒性，只保留峰值大于主方向峰值80％的方向作为该关键点的辅方向。因此，对于同一梯度值的多个峰值的关键点位置，在相同位置和尺度将会有多个关键点被创建但方向不同。仅有15％的关键点被赋予多个方向，但可以明显的提高关键点匹配的稳定性。实际编程实现中，就是把该关键点复制成多份关键点，并将方向值分别赋给这些复制后的关键点，并且，离散的梯度方向直方图要进行插值拟合处理，来求得更精确的方向角度值，检测结果如图5.2所示。

至此，将检测出的含有位置、尺度和方向的关键点即是该图像的SIFT特征点。

6、关键点特征描述

通过以上步骤，对于每一个关键点，拥有三个信息：位置、尺度以及方向。接下来就是为每个关键点建立一个描述符，用一组向量将这个关键点描述出来，使其不随各种变化而改变，比如光照变化、视角变化等等。这个描述子不但包括关键点，也包含关键点周围对其有贡献的像素点，并且描述符应该有较高的独特性，以便于提高特征点正确匹配的概率。

SIFT描述子是关键点邻域高斯图像梯度统计结果的一种表示。通过对关键点周围图像区域分块，计算块内梯度直方图，生成具有独特性的向量，这个向量是该区域图像信息的一种抽象，具有唯一性。

Lowe建议描述子使用在关键点尺度空间内4*4的窗口中计算的8个方向的梯度信息，共4*4*8=128维向量表征。表示步骤如下：

1. 确定计算描述子所需的图像区域

特征描述子与特征点所在的尺度有关，因此，对梯度的求取应在特征点对应的高斯图像上进行。将关键点附近的邻域划分为d*d(Lowe建议d=4)个子区域，每个子区域做为一个种子点，每个种子点有8个方向。每个子区域的大小与关键点方向分配时相同，即每个区域有个子像素，为每个子区域分配边长为的矩形区域进行采样(个子像素实际用边长为的矩形区域即可包含，但由式(3-8)，不大，为了简化计算取其边长为，并且采样点宜多不宜少)。考虑到实际计算时，需要采用双线性插值，所需图像窗口边长为。在考虑到旋转因素(方便下一步将坐标轴旋转到关键点的方向)，如下图6.1所示，实际计算所需的图像区域半径为：

(6-1)

计算结果四舍五入取整。

2. 将坐标轴旋转为关键点的方向，以确保旋转不变性，如6.2所示。

旋转后邻域内采样点的新坐标为：

(6-2)

3. 将邻域内的采样点分配到对应的子区域内，将子区域内的梯度值分配到8个方向上，计算其权值。

旋转后的采样点坐标在半径为radius的圆内被分配到的子区域，计算影响子区域的采样点的梯度和方向，分配到8个方向上。

旋转后的采样点落在子区域的下标为

(6-3)

Lowe建议子区域的像素的梯度大小按的高斯加权计算，即

(6-4)

其中a，b为关键点在高斯金字塔图像中的位置坐标。

4. 插值计算每个种子点八个方向的梯度。

如图6.3所示，将由式(6-3)所得采样点在子区域中的下标(图中蓝色窗口内红色点)线性插值，计算其对每个种子点的贡献。如图中的红色点，落在第0行和第1行之间，对这两行都有贡献。对第0行第3列种子点的贡献因子为dr，对第1行第3列的贡献因子为1-dr，同理，对邻近两列的贡献因子为dc和1-dc，对邻近两个方向的贡献因子为do和1-do。则最终累加在每个方向上的梯度大小为：

(6-5)

其中k，m，n为0或为1。

5. 如上统计的4*4*8=128个梯度信息即为该关键点的特征向量。特征向量形成后，为了去除光照变化的影响，需要对它们进行归一化处理，对于图像灰度值整体漂移，图像各点的梯度是邻域像素相减得到，所以也能去除。得到的描述子向量为，归一化后的特征向量为则

(6-7)

6. 描述子向量门限。非线性光照，相机饱和度变化对造成某些方向的梯度值过大，而对方向的影响微弱。因此设置门限值(向量归一化后，一般取0.2)截断较大的梯度值。然后，再进行一次归一化处理，提高特征的鉴别性。

7. 按特征点的尺度对特征描述向量进行排序。

至此，SIFT特征描述向量生成。

描述向量这块不好理解，我画了个草图，供参考：

7、SIFT的缺点

SIFT在图像的不变特征提取方面拥有无与伦比的优势，但并不完美，仍然存在：

1. 实时性不高。

2. 有时特征点较少。

3. 对边缘光滑的目标无法准确提取特征点。

等缺点，如下图7.1所示，对模糊的图像和边缘平滑的图像，检测出的特征点过少，对圆更是无能为力。近来不断有人改进，其中最著名的有SURF和CSIFT。

8、总结

本人研究SIFT算法一月有余，鉴于相关知识的缺失，尺度空间技术和差分近似求导曾困我良久。Lowe在论文中对细节提之甚少，甚至只字未提，给实现带来了很大困难。经过多方查阅，实现，总结成此文。自认为是到目前为止，关于SIFT算法最为详尽的资料，现分享给你，望批评指正。

一同分享给你的还有同时实现的高斯模糊源码，sift算法源码，见附录1，2。源码使用vs2010+opencv2.2实现。

zdd

2012年4月28日于北师大

2012年5月17日15:33:23第一次修正

修正内容：第3.3部分内容，图3.1，图3.5。

修正后代码：http://download.csdn.net/detail/zddmail/4309418

参考资料

1、David G.Lowe Distinctive Image Features from Scale-Invariant Keypoints. January 5, 2004.

2、David G.Lowe Object Recognition from Local Scale-Invariant Features. 1999

3、Matthew Brown and David Lowe Invariant Features from Interest Point Groups. In British Machine Vision Conference, Cardiff, Wales, pp. 656-665.

4、PETER J. BURT, MEMBER, IEEE, AND EDWARD H. ADELSON, The Laplacian Pyramid as a Compact Image Code. IEEE TRANSACTIONS ON COMMUNICATIONS, VOL. COM-3l, NO. 4, APRIL 1983

5、宋丹 10905056 尺度不变特征变换匹配算法Scale Invariant Feature Transform （SIFT）(PPT)

6、RaySaint 的博客SIFT算法研究http://underthehood.blog.51cto.com/2531780/658350

7、Jason Clemons SIFT: SCALE INVARIANT FEATURE TRANSFORM BY DAVID LOWE(ppt)

8、Tony Lindeberg Scale-space theory: A basic tool for analysing structures at different scales.1994

9、SIFT官网的Rob Hess <hess@eecs.oregonstate.edu> SIFT源码

10、Opencv2.2 Andrea Vedaldi(UCLA VisionLab)实现的SIFT源码 http://www.vlfeat.org/~vedaldi/code/siftpp.html, opencv2.3改用Rob Hess的源码

11、科学计算中的偏微分方程有限差分法杨乐主编

12、维基百科SIFT词条：http://zh.wikipedia.org/zh-cn/Scale-invariant_feature_transform

13、百度百科SIFT词条：http://baike.baidu.com/view/2832304.htm

14、其它互联网资料

附录1 高斯模糊源码

http://blog.csdn.net/zddmail/article/details/7450033

http://download.csdn.net/detail/zddmail/4217704

附录2 SIFT算法源码

http://download.csdn.net/detail/zddmail/4309418

资助：如果此渣文对大家有帮助，敬请资助, 一毛两毛的都行呀！不给媳妇筹钱的程序员不是好程序员！

1. 什么是斑点

斑点通常是指与周围有着颜色和灰度差别的区域。在实际地图中，往往存在着大量这样的斑点，如一颗树是一个斑点，一块草地是一个斑点，一栋房子也可以是一个斑点。由于斑点代表的是一个区域，相比单纯的角点，它的稳定性要好，抗噪声能力要强，所以它在图像配准上扮演了很重要的角色。

同时有时图像中的斑点也是我们关心的区域，比如在医学与生物领域，我们需要从一些X光照片或细胞显微照片中提取一些具有特殊意义的斑点的位置或数量。

比如下图中天空的飞机、向日葵的花盘、X线断层图像中的两个斑点。

在视觉领域，斑点检测的主要思路都是检测出图像中比它周围像素灰度值大或比周围灰度值小的区域。一般有两种方法来实现这一目标：

基于求导的微分方法，这类的方法称为微分检测器；
基于局部极值的分水岭算法。

这里我们重点介绍第一种方法，主要检测LOG斑点。而OpenCV中SimpleBlobDetector斑点检测算子就实现了第二种方法，我们这里也会介绍它的接口使用方法。

2. LOG斑点检测

2.1 基本原理

利用高斯拉普通拉斯（Laplace of Gaussian,LOG）算子检测图像斑点是一种十分常用的方法，对于二维高斯函数：

G (x, y; σ) = 1 2 π σ 2 e x p (- x 2 + y 2 2 σ 2)

它的拉普拉斯变换为:

\nabla 2 g = \partial 2 g \partial x 2 + \partial 2 g \partial y 2

规范化的高斯拉普变换为：

\nabla 2 n o r m = σ 2 \nabla 2 g = σ 2 (\partial 2 g \partial x 2 + \partial 2 g \partial y 2) = - 1 2 π σ 2 [1 - x 2 + y 2 σ 2] \cdot e x p (- x 2 + y 2 2 σ 2)

规范化算法子在二维图像上显示是一个圆对称函数，如下图所示。我们可以用这个算子来检测图像中的斑点，并且可以通过改变σ的值，可以检测不同尺寸的二维斑点。

2.2 LOG原理解释

其实从更直观的角度去解释为什么LOG算子可以检测图像中的斑点是：

图像与某一个二维函数进行卷积运算实际就是求取图像与这一函数的相似性。同理，图像与高斯拉普拉斯函数的卷积实际就是求取图像与高斯拉普拉斯函数的相似性。当图像中的斑点尺寸与高斯拉普拉斯函数的形状趋近一致时，图像的拉普拉斯响应达到最大。

从概率的角度解释为：假设原图像是一个与位置有关的随机变量X的密度函数，而LOG为随机变量Y的密度函数，则随机变量X+Y的密度分布函数即为两个函数的卷积形式(这一部分的理论，可以参见本博客概率与统计相关文章)。如果想让X+Y能取到最大值，则X与Y能保持步调一致最好，即X上升时，Y也上升，X最大时，Y也最大。

那么LOG算子是怎么被构想出来的呢？

事实上我们知道Laplace可以用来检测图像中的局部极值点，但是对噪声敏感，所以在我们对图像进行Laplace卷积之前，我们用一个高斯低通滤波对图像进行卷积，目标是去除图像中的噪声点。这一过程可以描述为:

先对图像f(x,y)用方差为σ的高斯核进行高斯滤波，去除图像中的噪点。

L (x, y; σ) = f (x, y) * G (x, y; σ)

然后对图像的拉普拉斯图像则为：

\nabla 2 = \partial 2 L \partial x 2 + \partial 2 L \partial y 2

而实际上有下面等式：

\nabla 2 [G (x, y) * f (x, y)] = \nabla 2 [G (x, y)] * f (x, y)

所以，我们可以先求高斯核的拉普拉斯算子，再对图像进行卷积。也就是一开始描述的步骤。

2.3 LOG算子的实现

Mat Feat::getHOGKernel(Size& ksize, double sigma){    Mat kernel(ksize, CV_64F);    Point centPoint = Point((ksize.width -1)/2, ((ksize.height -1)/2));    // first calculate Gaussian    for (int i=0; i < kernel.rows; i++)    {        double* pData = kernel.ptr<double>(i);        for (int j = 0; j < kernel.cols; j++)        {            double param = -((i - centPoint.y) * (i - centPoint.y) + (j - centPoint.x) * (j - centPoint.x)) / (2*sigma*sigma);            pData[j] = exp(param);        }    }    double maxValue;    minMaxLoc(kernel, NULL, &maxValue);    for (int i=0; i < kernel.rows; i++)    {        double* pData = kernel.ptr<double>(i);        for (int j = 0; j < kernel.cols; j++)        {            if (pData[j] < EPS* maxValue)            {                pData[j] = 0;            }        }    }    double sumKernel = sum(kernel)[0];    if (sumKernel != 0)    {        kernel = kernel / sumKernel;    }    // now calculate Laplacian    for (int i=0; i < kernel.rows; i++)    {        double* pData = kernel.ptr<double>(i);        for (int j = 0; j < kernel.cols; j++)        {            double addition = ((i - centPoint.y) * (i - centPoint.y) + (j - centPoint.x) * (j - centPoint.x) - 2*sigma*sigma)/(sigma*sigma*sigma*sigma);            pData[j] *= addition;        }    }    // make the filter sum to zero    sumKernel = sum(kernel)[0];    kernel -= (sumKernel/(ksize.width  * ksize.height));        return kernel;}

2.4 多尺度检测

我们注意到当σ尺度一定时，只能检测对应半径的斑点，那么检测的是多大半径的斑点呢，我们可以通过对规范化的二维拉普拉斯高斯算子求导：

规范化的高斯拉普拉斯函数为：

\nabla 2 n o r m = - 1 2 π σ 2 [1 - x 2 + y 2 σ 2] \cdot e x p (- x 2 + y 2 2 σ 2)

求∇2norm的极点值等价于求取下式：

\partial ( \nabla 2 n o r m ) \partial σ = 0

得到：

(x 2 + y 2 - 2 σ 2) \cdot e x p (- ( x 2 + y 2 ) 2 σ 2)

r 2 - 2 σ 2 = 0

对于图像中的斑点，在尺度σ=r/2√时，高斯拉普拉斯响应值达到最大。同理，如果图像中的圆形斑点黑白反向，那么，它的高斯拉普拉斯响应值在σ=r/2√时达到最小。将高斯拉普拉斯响应达到峰值时的尺度σ值，称为特征尺度。

那么在多尺度的情况下，同时在空间和尺度上达到最大值（或最小值）的点就是我们所期望的斑点。对于二维图像I(x,y)，计算图像在不同尺度下的离散拉普拉斯响应值，然后检查位置空间中的每个点；如果该点的拉普拉斯响应值都大小于或小于其他26个立方空间领域(9+8+9)的值，那么该点就是被检测到的图像斑点。

3 OpenCV进行斑点检测

opencv中检测Blobs的类为SimpleBlobDetector，这个类在opencv中的定义如下：

class SimpleBlobDetector : public FeatureDetector{public:struct Params{    Params();    float thresholdStep;    float minThreshold;    float maxThreshold;    size_t minRepeatability;    float minDistBetweenBlobs;    bool filterByColor;    uchar blobColor;    bool filterByArea;    float minArea, maxArea;    bool filterByCircularity;    float minCircularity, maxCircularity;    bool filterByInertia;    float minInertiaRatio, maxInertiaRatio;    bool filterByConvexity;    float minConvexity, maxConvexity;};SimpleBlobDetector(const SimpleBlobDetector::Params &parameters = SimpleBlobDetector::Params());protected:    ...};

算法的大致步骤如下：

对[minThreshold,maxThreshold)区间，以thresholdStep为间隔，做多次二值化。
对每张二值图片，使用findContours()提取连通域并计算每一个连通域的中心。
根据2得到的中心，全部放在一起。一些很接近的点［由theminDistBetweenBlobs控制多少才算接近］被归为一个group,对应一个bolb特征..
从3得到的那些点,估计最后的blob特征和相应半径，并以key points返回。

同时该支持提取特征的方法，一共有5个选项，这里就不多加描述了，默认是提取黑色圆形的Blob特征。下面是一个示例

int main(int argc, char** argv) {     Mat image = imread(argv[1]);     vector<KeyPoint> keyPoints;     SimpleBlobDetector::Params params;    SimpleBlobDetector blobDetect(params);     blobDetect.create("SimpleBlob");     blobDetect.detect(image, keyPoints);     cout << keyPoints.size() << endl;     drawKeypoints(image, keyPoints, image, Scalar(255,0,0));    namedWindow("blobs");     imshow("blobs", image);     waitKey();     return 0; }

总体来说，OpenCV的斑点检测效果还算不错，但是在有些图像的效果上明显不如LOG算子检测的检测效果。

4. 扩展阅读

一个与LOG滤波核近似的是高斯差分DOG滤波核，它的定义为：

D (x, y, σ) = (G (x, y, k σ) - G (x, y, σ)) * I (x, y) = L (x, y, k σ) - L (x, y, σ)

其中k为两个相邻尺度间的比例因子。

DOG可以看作为LOG的一个近似，但是它比LOG的效率更高。

前面介绍的微分算子在近圆的斑点检测方面效果很好，但是这些检测算子被限定于只能检测圆形斑点，而且不能估计斑点的方向，因为LOG算子等都是中心对称的。如果我们定义一种二维高斯核的变形，记它在X方向与Y方向上具有不同的方差，则这种算子可以用来检测带有方向的斑点。

G (x, y) = A \cdot e x p (- [(a x 2 + 2 b x y + c y 2)])

a = c o s 2 θ 2 σ 2 x + s i n 2 θ 2 σ 2 y, b = - s i n 2 θ 2 σ 2 x + s i n 2 θ 4 σ 2 y, c = s i n 2 θ 2 σ 2 x + c o s 2 θ 2 σ 2 y

其中A是规一性因子。

5. 参考资料

1. 《现代数字图像 -- 处理技术提高与应用案例详解》

2. 《图像局部不变性特征与描述》

3. Lindeberg, T. Feature Detection with Automatic Scale Selection

4. Hui Kong. A Generalized Laplacian Of Gaussian Filter for Blob Detection and Its Applications.

如果说SIFT算法中使用DOG对LOG进行了简化，提高了搜索特征点的速度，那么SURF算法则是对DoH的简化与近似。虽然SIFT算法已经被认为是最有效的，也是最常用的特征点提取的算法，但如果不借助于硬件的加速和专用图像处理器的配合，SIFT算法以现有的计算机仍然很难达到实时的程度。对于需要实时运算的场合，如基于特征点匹配的实时目标跟踪系统，每秒要处理8-24帧的图像，需要在毫秒级内完成特征点的搜索、特征矢量生成、特征矢量匹配、目标锁定等工作，这样SIFT算法就很难适应这种需求了。SURF借鉴了SIFT中简化近似的思想，把DoH中的高斯二阶微分模板进行了简化，使得模板对图像的滤波只需要进行几个简单的加减法运算，并且，这种运算与滤波器的尺度无关。实验证明，SURF算法较SIFT在运算速度上要快3倍左右。

1. 积分图像

SURF算法中要用到积分图像的概念。借助积分图像，图像与高斯二阶微分模板的滤波转化为对积分图像的加减运算。

积分图像中任意一点(i,j)的值ii(i,j)，为原图像左上角到点(i,j)相应的对角线区域灰度值的总和，即

i i (i, j) = \sum r \leq i, c \leq j p (r, c)

式中，p(r,c)表示图像中点(r,c)的灰度值，ii(i,j)可以用下面两式迭代计算得到

S (i, j) = S (i, j - 1) + p (i, j)

i i (i, j) = i i (i - 1, j) + S (i, j)

式中，S(i,j)表示一列的积分，且S(i,−1)=0,ii(−1,j)=0。求积分图像，只需要对原图像所有像素进行一遍扫描。

OpenCV中提供了用于计算积分图像的接口

/** src :输入图像，大小为M*N* sum: 输出的积分图像，大小为（M+1）*（N+1）* sdepth：用于指定sum的类型，-1表示与src类型一致*/void integral(InputArray src, OutputArray sum, int sdepth = -1);

值得注意的是OpenCV里的积分图大小比原图像多一行一列，那是因为OpenCV中积分图的计算公式为：

i i (i, j) = \sum r < i, c < j p (r, c)

一旦积分图计算好了，计算图像内任何矩形区域的像素值的和只需要三个加法，如上图所示。

2. DoH近似

在斑点检测这篇文章中已经提到过，我们可以利用Hessian矩阵行列式的极大值检测斑点。下面我们给出Hessian矩阵的定义。

给定图像I中的一个点x(i,j)，在点x处，尺度为σ的Hessian矩阵H(x,σ)定义如下：

H (x, σ) = [L x x (x, σ) L x y (x, σ) L x y (x, σ) L y y (x, σ)]

式中，Lxx(x,σ)是高斯二阶微分∂2g(σ)∂x2在点x处与图像I的卷积，Lx,y(x,σ)和Lyy(x,σ)具有类似的含义。

下面显示的是上面三种高斯微分算子的图形。

但是利用Hessian行列式进行图像斑点检测时，有一个缺点。由于二阶高斯微分被离散化和裁剪的原因，导致了图像在旋转奇数倍的π/4时，即转换到模板的对角线方向时，特征点检测的重复性降低（也就是说，原来特征点的地方，可能检测不到特征点了）。而在π/2时，特征点检测的重现率真最高。但这一小小的不足不影响我们使用Hessian矩阵进行特征点的检测。

为了将模板与图产像的卷积转换为盒子滤波运算，我们需要对高斯二阶微分模板进行简化，使得简化后的模板只是由几个矩形区域组成，矩形区域内填充同一值，如下图所示，在简化模板中白色区域的值为正数，黑色区域的值为负数，灰度区域的值为0。

对于σ=1.2的高斯二阶微分滤波器，我们设定模板的尺寸为9×9的大小，并用它作为最小尺度空间值对图像进行滤波和斑点检测。我们使用Dxx、Dxy和Dyy表示模板与图像进行卷积的结果。这样，便可以将Hessian矩阵的行列式作如下的简化。

D e t (H) = L x x L y y - L x y L x y = D x x L x x D x x D y y L y y D y y - D x y L x y D x y D x y L x y D x y \approx D x x D y y - (w D x y) 2

滤波器响应的相关权重w是为了平衡Hessian行列式的表示式。这是为了保持高斯核与近似高斯核的一致性。

w = | L x y ( σ ) | F | D x x ( σ ) F | | L x x ( σ ) | F | D x y ( σ ) F | = 0.912

其中|X|F为Frobenius范数。理论上来说对于不同的σ的值和对应尺寸的模板尺寸，w值是不同的，但为了简化起见，可以认为它是同一个常数。

使用近似的Hessian矩阵行列式来表示图像中某一点x处的斑点响应值，遍历图像中所有的像元点，便形成了在某一尺度下琉璃点检测的响应图像。使用不同的模板尺寸，便形成了多尺度斑点响应的金字塔图像，利用这一金字塔图像，就可以进行斑点响应极值点的搜索，其过程完全与SIFT算法类同。

3. 尺度空间表示

通常想要获取不同尺度的斑点，必须建立图像的尺度空间金字塔。一般的方法是通过不同σ的高斯函数，对图像进行平滑滤波，然后重采样图像以获得更高一层的金字塔图像。SIFT特征检测算法中就是通过相邻两层图像金字塔相减得到DoG图像，然后再在DoG图像上进行斑点和边缘检测工作的。

由于采用了盒子滤波和积分图像，所以，我们并不需要像SIFT算法那样去直接建立图像金字塔，而是采用不断增大盒子滤波模板的尺寸的间接方法。通过不同尺寸盒子滤波模板与积分图像求取Hessian矩阵行列式的响应图像。然后在响应图像上采用3D非最大值抑制，求取各种不同尺度的斑点。

如前所述，我们使用9×9的模板对图像进行滤波，其结果作为最初始的尺度空间层（此时，尺度值为s=1.2，近似σ=1.2的高斯微分），后续的层将通过逐步放大滤波模板尺寸，以及放大后的模板不断与图像进行滤波得到。由于采用盒子滤波和积分图像，滤波过程并不随着滤波模板尺寸的增加而使运算工作量增加。

与SIFT算法类似，我们需要将尺度空间划分为若干组（Octaves）。一个组代表了逐步放大的滤波模板对同一输入图像进行滤波的一系列响应图。每个组又由若干固定的层组成。由于积分图像离散化的原因，两个层之间的最小尺度变化量是由高斯二阶微分滤波器在微分方向上对正负斑点响应长度l0决定的，它是盒子滤波器模板尺寸的1/3。对于9×9的模板，它的l0=3。一下层的响应长度至少应该在l0的基础上增加2个像素，以保证一边一个像素，即l0=5。这样模板的尺寸就为15×15。以此类推，我们可以得到一个尺寸增大模板序列，它们的尺寸分别为：9×9，15×15，21×21，27×27，黑色、白色区域的长度增加偶数个像素，以保证一个中心像素的存在。

采用类似的方法来处理其他几组的模板序列。其方法是将滤波器尺寸增加量翻倍（6，12，24，38）。这样，可以得到第二组的滤波器尺寸，它们分别为15，27，39，51。第三组的滤波器尺寸为27，51，75，99。如果原始图像的尺寸仍然大于对应的滤波器尺寸，尺度空间的分析还可以进行第四组，其对应的模板尺寸分别为51，99，147和195。下图显示了第一组至第三组的滤波器尺寸变化。

在通常尺度分析情况下，随着尺度的增大，被检测到的斑点数量迅速衰减。所以一般进行3-4组就可以了，与此同时，为了减少运算量，提高计算的速度，可以考虑在滤波时，将采样间隔设为2。

对于尺寸为L的模板，当用它与积分图运算来近似二维高斯核的滤波时，对应的二维高斯核的参数σ=1.2×L9，这一点至关重要，尤其是在后面计算描述子时，用于计算邻域的半径时。

4. 兴趣点的定位

为了在图像及不同尺寸中定位兴趣点，我们用了3×3×3邻域非最大值抑制。具体的步骤基本与SIFT一致，而且Hessian矩阵行列式的最大值在尺度和图像空间被插值。

下面显示了我们用的快速Hessian检测子检测到的兴趣点。

5. SURF源码解析

这份源码来自OpenCV nonfree模块。

这里先介绍SURF特征点定位这一块，关于特征点的描述下一篇文章再介绍。

5.1 主干函数 fastHessianDetector

特征点定位的主干函数为fastHessianDetector，该函数接受一个积分图像，以及尺寸相关的参数，组数与每组的层数，检测到的特征点保存在vector<KeyPoint>类型的结构中。

static void fastHessianDetector(const Mat& sum, const Mat& msum, vector<KeyPoint>& keypoints,    int nOctaves, int nOctaveLayers, float hessianThreshold){    /*first Octave图像采样的步长，第二组的时候加倍，以此内推    增加这个值，将会加快特征点检测的速度，但是会让特征点的提取变得不稳定*/    const int SAMPLE_STEP0 = 1;    int nTotalLayers = (nOctaveLayers + 2)*nOctaves; // 尺度空间的总图像数    int nMiddleLayers = nOctaveLayers*nOctaves; // 用于检测特征点的层的 总数，也就是中间层的总数    vector<Mat> dets(nTotalLayers); // 每一层图像 对应的 Hessian行列式的值    vector<Mat> traces(nTotalLayers); // 每一层图像 对应的 Hessian矩阵的迹的值    vector<int> sizes(nTotalLayers); // 每一层用的 Harr模板的大小    vector<int> sampleSteps(nTotalLayers); // 每一层用的采样步长     vector<int> middleIndices(nMiddleLayers); // 中间层的索引值    keypoints.clear();    // 为上面的对象分配空间，并赋予合适的值    int index = 0, middleIndex = 0, step = SAMPLE_STEP0;    for (int octave = 0; octave < nOctaves; octave++)    {        for (int layer = 0; layer < nOctaveLayers + 2; layer++)        {            /*这里sum.rows - 1是因为 sum是积分图，它的大小是原图像大小加1*/            dets[index].create((sum.rows - 1) / step, (sum.cols - 1) / step, CV_32F); // 这里面有除以遍历图像用的步长            traces[index].create((sum.rows - 1) / step, (sum.cols - 1) / step, CV_32F);            sizes[index] = (SURF_HAAR_SIZE0 + SURF_HAAR_SIZE_INC*layer) << octave;            sampleSteps[index] = step;            if (0 < layer && layer <= nOctaveLayers)                middleIndices[middleIndex++] = index;            index++;        }        step *= 2;    }    // Calculate hessian determinant and trace samples in each layer    for (int i = 0; i < nTotalLayers; i++)    {        calcLayerDetAndTrace(sum, sizes[i], sampleSteps[i], dets[i], traces[i]);    }    // Find maxima in the determinant of the hessian    for (int i = 0; i < nMiddleLayers; i++)    {        int layer = middleIndices[i];        int octave = i / nOctaveLayers;        findMaximaInLayer(sum, msum, dets, traces, sizes, keypoints, octave, layer, hessianThreshold, sampleSteps[layer]);    }    std::sort(keypoints.begin(), keypoints.end(), KeypointGreater());}

5.2 计算Hessian矩阵的行列式与迹calcLayerDetAndTrace

这个函数首先定义了尺寸为9的第一层图像的三个模板。模板分别为一个3×5、3×5、4×5的二维数组表示，数组的每一行表示一个黑白块的位置参数。函数里只初始化了第一层图像的模板参数，后面其他组其他层的Harr模板参数都是用resizeHaarPattern这个函数来计算的。这个函数返回的是一个SurfHF的结构体，这个结构体由两个点及一个权重构成。

struct SurfHF{    int p0, p1, p2, p3;    float w;    SurfHF() : p0(0), p1(0), p2(0), p3(0), w(0) {}};

resizeHaarPattern这个函数非常的巧妙，它把模板中的点坐标。转换到在积分图中的相对（模板左上角点）坐标。

static voidresizeHaarPattern(const int src[][5], SurfHF* dst, int n, int oldSize, int newSize, int widthStep){    float ratio = (float)newSize / oldSize;    for (int k = 0; k < n; k++)    {        int dx1 = cvRound(ratio*src[k][0]);        int dy1 = cvRound(ratio*src[k][1]);        int dx2 = cvRound(ratio*src[k][2]);        int dy2 = cvRound(ratio*src[k][3]);        /*巧妙的坐标转换*/        dst[k].p0 = dy1*widthStep + dx1; // 转换为一个相对距离，距离模板左上角点的  在积分图中的距离 !!important!!        dst[k].p1 = dy2*widthStep + dx1;         dst[k].p2 = dy1*widthStep + dx2;        dst[k].p3 = dy2*widthStep + dx2;        dst[k].w = src[k][4] / ((float)(dx2 - dx1)*(dy2 - dy1));// 原来的+1,+2用 覆盖的所有像素点平均。    }}

在用积分图计算近似卷积时，用的是calcHaarPattern函数。这个函数比较简单，只用知道左上与右下角坐标即可。

inline float calcHaarPattern(const int* origin, const SurfHF* f, int n){    /*orgin即为积分图，n为模板中 黑白 块的个数 */    double d = 0;    for (int k = 0; k < n; k++)        d += (origin[f[k].p0] + origin[f[k].p3] - origin[f[k].p1] - origin[f[k].p2])*f[k].w;    return (float)d;}

最终我们可以看到了整个calcLayerDetAndTrack的代码

static void calcLayerDetAndTrace(const Mat& sum, int size, int sampleStep,    Mat& det, Mat& trace){    const int NX = 3, NY = 3, NXY = 4;    const int dx_s[NX][5] = { { 0, 2, 3, 7, 1 }, { 3, 2, 6, 7, -2 }, { 6, 2, 9, 7, 1 } };    const int dy_s[NY][5] = { { 2, 0, 7, 3, 1 }, { 2, 3, 7, 6, -2 }, { 2, 6, 7, 9, 1 } };    const int dxy_s[NXY][5] = { { 1, 1, 4, 4, 1 }, { 5, 1, 8, 4, -1 }, { 1, 5, 4, 8, -1 }, { 5, 5, 8, 8, 1 } };    SurfHF Dx[NX], Dy[NY], Dxy[NXY];    if (size > sum.rows - 1 || size > sum.cols - 1)        return;    resizeHaarPattern(dx_s, Dx, NX, 9, size, sum.cols);    resizeHaarPattern(dy_s, Dy, NY, 9, size, sum.cols);    resizeHaarPattern(dxy_s, Dxy, NXY, 9, size, sum.cols);    /* The integral image 'sum' is one pixel bigger than the source image */    int samples_i = 1 + (sum.rows - 1 - size) / sampleStep; // 最大能遍历到的 行坐标，因为要减掉一个模板的尺寸    int samples_j = 1 + (sum.cols - 1 - size) / sampleStep; // 最大能遍历到的 列坐标    /* Ignore pixels where some of the kernel is outside the image */    int margin = (size / 2) / sampleStep;    for (int i = 0; i < samples_i; i++)    {        /*坐标为(i,j)的点是模板左上角的点，所以实际现在模板分析是的i+margin,j+margin点处的响应*/        const int* sum_ptr = sum.ptr<int>(i*sampleStep);        float* det_ptr = &det.at<float>(i + margin, margin); // 左边空隙为 margin        float* trace_ptr = &trace.at<float>(i + margin, margin);        for (int j = 0; j < samples_j; j++)        {            float dx = calcHaarPattern(sum_ptr, Dx, 3);            float dy = calcHaarPattern(sum_ptr, Dy, 3);            float dxy = calcHaarPattern(sum_ptr, Dxy, 4);            sum_ptr += sampleStep;            det_ptr[j] = dx*dy - 0.81f*dxy*dxy;            trace_ptr[j] = dx + dy;        }    }}

5.3 局部最大值搜索findMaximaInLayer

这里算法思路很简单，值得注意的是里面的一些坐标的转换很巧妙，里面比较重的函数就是interpolateKeypoint函数，通过插值计算最大值点。

/** Maxima location interpolation as described in "Invariant Features from* Interest Point Groups" by Matthew Brown and David Lowe. This is performed by* fitting a 3D quadratic to a set of neighbouring samples.** The gradient vector and Hessian matrix at the initial keypoint location are* approximated using central differences. The linear system Ax = b is then* solved, where A is the Hessian, b is the negative gradient, and x is the* offset of the interpolated maxima coordinates from the initial estimate.* This is equivalent to an iteration of Netwon's optimisation algorithm.** N9 contains the samples in the 3x3x3 neighbourhood of the maxima* dx is the sampling step in x* dy is the sampling step in y* ds is the sampling step in size* point contains the keypoint coordinates and scale to be modified** Return value is 1 if interpolation was successful, 0 on failure.*/static intinterpolateKeypoint(float N9[3][9], int dx, int dy, int ds, KeyPoint& kpt){    Vec3f b(-(N9[1][5] - N9[1][3]) / 2,  // Negative 1st deriv with respect to x        -(N9[1][7] - N9[1][1]) / 2,  // Negative 1st deriv with respect to y        -(N9[2][4] - N9[0][4]) / 2); // Negative 1st deriv with respect to s    Matx33f A(        N9[1][3] - 2 * N9[1][4] + N9[1][5],            // 2nd deriv x, x        (N9[1][8] - N9[1][6] - N9[1][2] + N9[1][0]) / 4, // 2nd deriv x, y        (N9[2][5] - N9[2][3] - N9[0][5] + N9[0][3]) / 4, // 2nd deriv x, s        (N9[1][8] - N9[1][6] - N9[1][2] + N9[1][0]) / 4, // 2nd deriv x, y        N9[1][1] - 2 * N9[1][4] + N9[1][7],            // 2nd deriv y, y        (N9[2][7] - N9[2][1] - N9[0][7] + N9[0][1]) / 4, // 2nd deriv y, s        (N9[2][5] - N9[2][3] - N9[0][5] + N9[0][3]) / 4, // 2nd deriv x, s        (N9[2][7] - N9[2][1] - N9[0][7] + N9[0][1]) / 4, // 2nd deriv y, s        N9[0][4] - 2 * N9[1][4] + N9[2][4]);           // 2nd deriv s, s    Vec3f x = A.solve(b, DECOMP_LU);    bool ok = (x[0] != 0 || x[1] != 0 || x[2] != 0) &&        std::abs(x[0]) <= 1 && std::abs(x[1]) <= 1 && std::abs(x[2]) <= 1;    if (ok)    {        kpt.pt.x += x[0] * dx;        kpt.pt.y += x[1] * dy;        kpt.size = (float)cvRound(kpt.size + x[2] * ds);    }    return ok;}static void findMaximaInLayer(const Mat& sum, const Mat& mask_sum,    const vector<Mat>& dets, const vector<Mat>& traces,    const vector<int>& sizes, vector<KeyPoint>& keypoints,    int octave, int layer, float hessianThreshold, int sampleStep){    // Wavelet Data    const int NM = 1;    const int dm[NM][5] = { { 0, 0, 9, 9, 1 } };    SurfHF Dm;    int size = sizes[layer];    // 当前层图像的大小    int layer_rows = (sum.rows - 1) / sampleStep;    int layer_cols = (sum.cols - 1) / sampleStep;    // 边界区域大小，考虑的下一层的模板大小    int margin = (sizes[layer + 1] / 2) / sampleStep + 1;    if (!mask_sum.empty())        resizeHaarPattern(dm, &Dm, NM, 9, size, mask_sum.cols);    int step = (int)(dets[layer].step / dets[layer].elemSize());    for (int i = margin; i < layer_rows - margin; i++)    {        const float* det_ptr = dets[layer].ptr<float>(i);        const float* trace_ptr = traces[layer].ptr<float>(i);        for (int j = margin; j < layer_cols - margin; j++)        {            float val0 = det_ptr[j]; // 中心点的值            if (val0 > hessianThreshold)            {                // 模板左上角的坐标                int sum_i = sampleStep*(i - (size / 2) / sampleStep);                int sum_j = sampleStep*(j - (size / 2) / sampleStep);                /* The 3x3x3 neighbouring samples around the maxima.                The maxima is included at N9[1][4] */                const float *det1 = &dets[layer - 1].at<float>(i, j);                const float *det2 = &dets[layer].at<float>(i, j);                const float *det3 = &dets[layer + 1].at<float>(i, j);                float N9[3][9] = { { det1[-step - 1], det1[-step], det1[-step + 1],                    det1[-1], det1[0], det1[1],                    det1[step - 1], det1[step], det1[step + 1] },                    { det2[-step - 1], det2[-step], det2[-step + 1],                    det2[-1], det2[0], det2[1],                    det2[step - 1], det2[step], det2[step + 1] },                    { det3[-step - 1], det3[-step], det3[-step + 1],                    det3[-1], det3[0], det3[1],                    det3[step - 1], det3[step], det3[step + 1] } };                /* Check the mask - why not just check the mask at the center of the wavelet? */                if (!mask_sum.empty())                {                    const int* mask_ptr = &mask_sum.at<int>(sum_i, sum_j);                    float mval = calcHaarPattern(mask_ptr, &Dm, 1);                    if (mval < 0.5)                        continue;                }                /* 检测val0，是否在N9里极大值，？？为什么不检测极小值呢*/                if (val0 > N9[0][0] && val0 > N9[0][1] && val0 > N9[0][2] &&                    val0 > N9[0][3] && val0 > N9[0][4] && val0 > N9[0][5] &&                    val0 > N9[0][6] && val0 > N9[0][7] && val0 > N9[0][8] &&                    val0 > N9[1][0] && val0 > N9[1][1] && val0 > N9[1][2] &&                    val0 > N9[1][3] && val0 > N9[1][5] &&                    val0 > N9[1][6] && val0 > N9[1][7] && val0 > N9[1][8] &&                    val0 > N9[2][0] && val0 > N9[2][1] && val0 > N9[2][2] &&                    val0 > N9[2][3] && val0 > N9[2][4] && val0 > N9[2][5] &&                    val0 > N9[2][6] && val0 > N9[2][7] && val0 > N9[2][8])                {                    /* Calculate the wavelet center coordinates for the maxima */                    float center_i = sum_i + (size - 1)*0.5f;                    float center_j = sum_j + (size - 1)*0.5f;                    KeyPoint kpt(center_j, center_i, (float)sizes[layer],                        -1, val0, octave, CV_SIGN(trace_ptr[j]));                    /* 局部极大值插值，用Hessian，类似于SIFT里的插值，里面没有迭代5次，只进行了一次查找，why?  */                    int ds = size - sizes[layer - 1];                    int interp_ok = interpolateKeypoint(N9, sampleStep, sampleStep, ds, kpt);                    /* Sometimes the interpolation step gives a negative size etc. */                    if (interp_ok)                    {                        /*printf( "KeyPoint %f %f %d\n", point.pt.x, point.pt.y, point.size );*/                        keypoints.push_back(kpt);                    }                }            }        }    }}

6. 总结

总体来说，如果理解了SIFT算法，再来看SURF算法会发现思路非常简单。尤其是局部最大值查找方面，基本一致。关键还是一个用积分图来简化卷积的思路，以及怎么用不同的模板来近似原来尺度空间中的高斯滤波器。

这一篇主要讨论分析的是SURF的定位问题，下面还有SURF特征点的方向计算与描述子的生成，将在下一篇文章中详细描述。

上一篇文章 SURF算法与源码分析、上中主要分析的是SURF特征点定位的算法原理与相关OpenCV中的源码分析，这篇文章接着上篇文章对已经定位到的SURF特征点进行特征描述。这一步至关重要，这是SURF特征点匹配的基础。总体来说算法思路和SIFT相似，只是每一步都做了不同程度的近似与简化，提高了效率。

1. SURF特征点方向分配

为了保证特征矢量具有旋转不变性，与SIFT特征一样，需要对每个特征点分配一个主方向。为些，我们需要以特征点为中心，以6s（s=1.2∗L/9为特征点的尺度）为半径的圆形区域，对图像进行Haar小波响应运算。这样做实际就是对图像进行梯度运算只不过是我们需要利用积分图像，提高计算图像梯度的效率。在SIFT特征描述子中我们在求取特征点主方向时，以是特征点为中心，在以4.5σ为半径的邻域内计算梯度方向直方图。事实上，两种方法在求取特征点主方向时，考虑到Haar小波的模板带宽，实际计算梯度的图像区域是相同的。用于计算梯度的Harr小波的尺度为4s。

与SIFT类似，使用σ=2s的高斯加权函数对Harr小波的响应值进行高斯加权。为了求取主方向值，需要设计一个以特征点为中心，张角为π/3的扇形滑动窗口。以步长为0.2弧度左右，旋转这个滑动窗口，并对滑动窗口内的图像Harr小波响应值dx、dy进行累加，得到一个矢量(mw,θw)：

m w = \sum w d x + \sum w d y

θ w = a r c t a n (\sum w d x / \sum w d y)

主方向为最大Harr响应累加值所对应的方向，也就是最长矢量所对应的方向，即

θ = θ w | m a x {m w}

可以依照SIFT求方方向时策略，当存在另一个相当于主峰值80%能量的峰值时，则将这个方向认为是该特征点的辅方向。一个特征点可能会被指定具有多个方向（一个主方向，一个以上辅方向），这可以增强匹配的鲁棒性。和SIFT的描述子类似，如果在mw中出现另一个大于主峰能量max{mw}80时的次峰，可以将该特征点复制成两个特征点。一个主的方向为最大响应能量所对应的方向，另一个主方向为次大响应能量所对应的方向。

图 1 求取主方向时扇形滑动窗口围绕特征点转动，统计Haar小波响应值，并计算方向角

2. 特征点特征矢量生成

生成特征点描述子与确定特征点方向有些类似，它需要计算图像的Haar小波响应。不过，与主方向的确定不同的是，这次我们不是使用一个圆形区域，而是在一个矩形区域来计算Haar小波响应。以特征点为中心，沿上一节讨论得到的主方向，沿主方向将s20s×20s的图像划分为4×4个子块，每个子块利用尺寸2s的Harr模板进行响应值进行响应值计算，然后对响应值进行统计∑dx、∑|dx|、∑dy、∑|dy|形成特征矢量。如下图2所示。图中，以特征点为中心，以20s为边长的矩形窗口为特征描述子计算使用的窗口，特征点到矩形边框的线段表示特征点的主方向。

图2 特征描述子表示

将20s的窗口划分成4×4子窗口，每个子窗口有5s×5s个像素。使用尺寸为2s的Harr小波对子窗口图像进行其响应值计算，共进行25次采样，分别得到沿主方向的dy和垂直于主方向的dx。然后，以特征点为中心，对dy和dx进行高斯加权计算，高斯核的参数为σ=3.3s(即20s/6)。最后，分别对每个子块的响应值进行统计，得到每个子块的矢量：

V 子 块 = [\sum d x, \sum | d x |, \sum d y, \sum | d y |]

由于共有4×4个子块，因此，特征描述子共由4×4×4=64维特征矢量组成。SURF描述子不仅具有尺度和旋转不变性，而且对光照的变化也具有不变性。使小波响应本身就具有亮度不变性，而对比度的不变性则是通过将特征矢量进行归一化来实现。图3 给出了三种不同图像模式的子块得到的不同结果。对于实际图像的描述子，我们可以认为它们是由这三种不同模式图像的描述子组合而成的。

图3 不同的图像密度模式得到的不同的描述子结果

为了充分利用积分图像进行Haar小波的响应计算，我们并不直接旋转Haar小波模板求得其响应值，而是在积图像上先使用水平和垂直的Haar模板求得响应值dy和dx，然后根据主方向旋转dx和dy与主方向操持一致，如下图4所示。为了求得旋转后Haar小波响应值，首先要得到旋转前图像的位置。旋转前后图偈的位置关系，可以通过点的旋转公式得到：

x = x 0 - j \times s c a l e \times s i n (θ) + i \times s c a l e \times c o s (θ)

y = y 0 - j \times s c a l e \times c o s (θ) + i \times s c a l e \times s i n (θ)

在得到点(j,i)在旋转前对应积分图像的位置(x,y)后，利用积分图像与水平、垂直Harr小波，求得水平与垂直两个方向的响应值dx和dy。对dx和dy进行高斯加权处理，并根据主方向的角度，对dx和dy进行旋转变换，从而，得到旋转后的dx’和dy’。其计算公式如下：

d x' = w (- d x \times s i n (θ) + d y \times c o s (θ))

d y' = w (- d x \times c o s (θ) + d y \times s i n (θ))

图4 利用积分图像进行Haar小波响应计算示意图，左边为旋转后的图像，右边为旋转前的图像

3. 特征描述子的维数

一般而言，特征矢量的长度越长，特征矢量所承载的信息量就越大，特征描述子的独特性就越好，但匹配时所付出的时间代价就越大。对于SURF描述子，可以将它扩展到用128维矢量来表示。具体方法是在求∑dx、∑|dx|时区分dy<0和dy≥0情况。同时，在求取∑dy、∑|dy|时区分dx<0和dx≥0情况。这样，每个子块就产生了8个梯度统计值，从而使描述子特征矢量的长度增加到8×4×4=128维。

为了实现快速匹配，SURF在特征矢量中增加了一个新的变量，即特征点的拉普拉斯响应正负号。在特征点检测时，将Hessian矩阵的迹的正负号记录下来，作为特征矢量中的一个变量。这样做并不增加运算量，因为特征点检测进已经对Hessian矩阵的迹进行了计算。在特征匹配时，这个变量可以有效地节省搜索的时间，因为只有两个具有相同正负号的特征点才有可能匹配，对于正负号不同的特征点就不进行相似性计算。

简单地说，我们可以根据特征点的响应值符号，将特征点分成两组，一组是具有拉普拉斯正响应的特征点，一组是具有拉普拉斯负响应的特征点，匹配时，只有符号相同组中的特征点才能进行相互匹配。显然，这样可以节省特征点匹配的时间。如下图5所示。

图5 黑背景下的亮斑和白背景下的黑斑因为它们的拉普拉斯响应正负号不同，不会对它们进行匹配

4. 源码解析

特征点描述子的生成这一部分的代码主要是通过SURFInvoker这个类来实现。在主流程中，通过一个parallel_for_()函数来并发计算。

struct SURFInvoker{    enum{ORI_RADIUS = 6, ORI_WIN = 60, PATCH_SZ = 20};    // Parameters    const Mat* img;    const Mat* sum;    vector<KeyPoint>* keypoints;    Mat* descriptors;    bool extended;    bool upright;    // Pre-calculated values    int nOriSamples;    vector<Point> apt; // 特征点周围用于描述方向的邻域的点    vector<float> aptw; // 描述 方向时的 高斯 权    vector<float> DW;    SURFInvoker(const Mat& _img, const Mat& _sum,        vector<KeyPoint>& _keypoints, Mat& _descriptors,        bool _extended, bool _upright)    {        keypoints = &_keypoints;        descriptors = &_descriptors;        img = &_img;        sum = &_sum;        extended = _extended;        upright = _upright;        // 用于描述特征点的 方向的 邻域大小： 12*sigma+1 (sigma =1.2) 因为高斯加权的核的参数为2sigma        // nOriSampleBound为 矩形框内点的个数        const int nOriSampleBound = (2 * ORI_RADIUS + 1)*(2 * ORI_RADIUS + 1); // 这里把s近似为1 ORI_DADIUS = 6s        // 分配大小         apt.resize(nOriSampleBound);        aptw.resize(nOriSampleBound);        DW.resize(PATCH_SZ*PATCH_SZ); // PATHC_SZ为特征描述子的 区域大小 20s(s 这里初始为1了)        /* 计算特征点方向用的 高斯分布 权值与坐标 */        Mat G_ori = getGaussianKernel(2 * ORI_RADIUS + 1, SURF_ORI_SIGMA, CV_32F); // SURF_ORI_SIGMA = 1.2 *2 =2.5        nOriSamples = 0;        for (int i = -ORI_RADIUS; i <= ORI_RADIUS; i++)        {            for (int j = -ORI_RADIUS; j <= ORI_RADIUS; j++)            {                if (i*i + j*j <= ORI_RADIUS*ORI_RADIUS) // 限制在圆形区域内                {                    apt[nOriSamples] = cvPoint(i, j);                    // 下面这里有个坐标转换，因为i,j都是从-ORI_RADIUS开始的。                    aptw[nOriSamples++] = G_ori.at<float>(i + ORI_RADIUS, 0) * G_ori.at<float>(j + ORI_RADIUS, 0);                }            }        }        CV_Assert(nOriSamples <= nOriSampleBound); // nOriSamples为圆形区域内的点，nOriSampleBound是正方形区域的点        /* 用于特征点描述子的高斯 权值 */        Mat G_desc = getGaussianKernel(PATCH_SZ, SURF_DESC_SIGMA, CV_32F); // 用于生成特征描述子的 高斯加权 sigma = 3.3s (s初取1)        for (int i = 0; i < PATCH_SZ; i++)        {            for (int j = 0; j < PATCH_SZ; j++)                DW[i*PATCH_SZ + j] = G_desc.at<float>(i, 0) * G_desc.at<float>(j, 0);        }        /* x与y方向上的 Harr小波,参数为4s */        const int NX = 2, NY = 2;        const int dx_s[NX][5] = { { 0, 0, 2, 4, -1 }, { 2, 0, 4, 4, 1 } };        const int dy_s[NY][5] = { { 0, 0, 4, 2, 1 }, { 0, 2, 4, 4, -1 } };        float X[nOriSampleBound], Y[nOriSampleBound], angle[nOriSampleBound]; // 用于计算特生点主方向        uchar PATCH[PATCH_SZ + 1][PATCH_SZ + 1];        float DX[PATCH_SZ][PATCH_SZ], DY[PATCH_SZ][PATCH_SZ]; // 20s * 20s区域的 梯度值        CvMat matX = cvMat(1, nOriSampleBound, CV_32F, X);        CvMat matY = cvMat(1, nOriSampleBound, CV_32F, Y);        CvMat _angle = cvMat(1, nOriSampleBound, CV_32F, angle);        Mat _patch(PATCH_SZ + 1, PATCH_SZ + 1, CV_8U, PATCH);        int dsize = extended ? 128 : 64;        int k, k1 = 0, k2 = (int)(*keypoints).size();// k2为Harr小波的 模板尺寸        float maxSize = 0;        for (k = k1; k < k2; k++)        {            maxSize = std::max(maxSize, (*keypoints)[k].size);        }        // maxSize*1.2/9 表示最大的尺度 s        int imaxSize = std::max(cvCeil((PATCH_SZ + 1)*maxSize*1.2f / 9.0f), 1);        Ptr<CvMat> winbuf = cvCreateMat(1, imaxSize*imaxSize, CV_8U);        for (k = k1; k < k2; k++)        {            int i, j, kk, nangle;            float* vec;            SurfHF dx_t[NX], dy_t[NY];            KeyPoint& kp = (*keypoints)[k];            float size = kp.size;            Point2f center = kp.pt;            /* s是当前层的尺度参数 1.2是第一层的参数，9是第一层的模板大小*/            float s = size*1.2f / 9.0f;            /* grad_wav_size是 harr梯度模板的大小 边长为 4s */            int grad_wav_size = 2 * cvRound(2 * s);            if (sum->rows < grad_wav_size || sum->cols < grad_wav_size)            {                /* when grad_wav_size is too big,                * the sampling of gradient will be meaningless                * mark keypoint for deletion. */                kp.size = -1;                continue;            }            float descriptor_dir = 360.f - 90.f;            if (upright == 0)            {                // 这一步 是计算梯度值，先将harr模板放大，再根据积分图计算，与前面求D_x，D_y一致类似                resizeHaarPattern(dx_s, dx_t, NX, 4, grad_wav_size, sum->cols);                resizeHaarPattern(dy_s, dy_t, NY, 4, grad_wav_size, sum->cols);                for (kk = 0, nangle = 0; kk < nOriSamples; kk++)                {                    int x = cvRound(center.x + apt[kk].x*s - (float)(grad_wav_size - 1) / 2);                    int y = cvRound(center.y + apt[kk].y*s - (float)(grad_wav_size - 1) / 2);                    if (y < 0 || y >= sum->rows - grad_wav_size ||                        x < 0 || x >= sum->cols - grad_wav_size)                        continue;                    const int* ptr = &sum->at<int>(y, x);                    float vx = calcHaarPattern(ptr, dx_t, 2);                    float vy = calcHaarPattern(ptr, dy_t, 2);                    X[nangle] = vx*aptw[kk];                    Y[nangle] = vy*aptw[kk];                    nangle++;                }                if (nangle == 0)                {                    // No gradient could be sampled because the keypoint is too                    // near too one or more of the sides of the image. As we                    // therefore cannot find a dominant direction, we skip this                    // keypoint and mark it for later deletion from the sequence.                    kp.size = -1;                    continue;                }                matX.cols = matY.cols = _angle.cols = nangle;                // 计算邻域内每个点的 梯度角度                cvCartToPolar(&matX, &matY, 0, &_angle, 1);                float bestx = 0, besty = 0, descriptor_mod = 0;                for (i = 0; i < 360; i += SURF_ORI_SEARCH_INC) // SURF_ORI_SEARCH_INC 为扇形区域扫描的步长                {                    float sumx = 0, sumy = 0, temp_mod;                    for (j = 0; j < nangle; j++)                    {                        // d是 分析到的那个点与 现在主方向的偏度                        int d = std::abs(cvRound(angle[j]) - i);                        if (d < ORI_WIN / 2 || d > 360 - ORI_WIN / 2)                        {                            sumx += X[j];                            sumy += Y[j];                        }                    }                    temp_mod = sumx*sumx + sumy*sumy;                    // descriptor_mod 是最大峰值                    if (temp_mod > descriptor_mod)                    {                        descriptor_mod = temp_mod;                        bestx = sumx;                        besty = sumy;                    }                }                descriptor_dir = fastAtan2(-besty, bestx);            }            kp.angle = descriptor_dir;            if (!descriptors || !descriptors->data)                continue;            /* 用特征点周围20*s为边长的邻域 计算特征描述子 */            int win_size = (int)((PATCH_SZ + 1)*s);            CV_Assert(winbuf->cols >= win_size*win_size);            Mat win(win_size, win_size, CV_8U, winbuf->data.ptr);            if (!upright)            {                descriptor_dir *= (float)(CV_PI / 180); // 特征点的主方向 弧度值                float sin_dir = -std::sin(descriptor_dir); //  - sin dir                float cos_dir = std::cos(descriptor_dir);                float win_offset = -(float)(win_size - 1) / 2;                float start_x = center.x + win_offset*cos_dir + win_offset*sin_dir;                float start_y = center.y - win_offset*sin_dir + win_offset*cos_dir;                uchar* WIN = win.data;                int ncols1 = img->cols - 1, nrows1 = img->rows - 1;                size_t imgstep = img->step;                for (i = 0; i < win_size; i++, start_x += sin_dir, start_y += cos_dir)                {                    double pixel_x = start_x;                    double pixel_y = start_y;                    for (j = 0; j < win_size; j++, pixel_x += cos_dir, pixel_y -= sin_dir)                    {                        int ix = cvFloor(pixel_x), iy = cvFloor(pixel_y);                        if ((unsigned)ix < (unsigned)ncols1 &&                            (unsigned)iy < (unsigned)nrows1)                        {                            float a = (float)(pixel_x - ix), b = (float)(pixel_y - iy);                            const uchar* imgptr = &img->at<uchar>(iy, ix);                            WIN[i*win_size + j] = (uchar)                                cvRound(imgptr[0] * (1.f - a)*(1.f - b) +                                imgptr[1] * a*(1.f - b) +                                imgptr[imgstep] * (1.f - a)*b +                                imgptr[imgstep + 1] * a*b);                        }                        else                        {                            int x = std::min(std::max(cvRound(pixel_x), 0), ncols1);                            int y = std::min(std::max(cvRound(pixel_y), 0), nrows1);                            WIN[i*win_size + j] = img->at<uchar>(y, x);                        }                    }                }            }            else            {                float win_offset = -(float)(win_size - 1) / 2;                int start_x = cvRound(center.x + win_offset);                int start_y = cvRound(center.y - win_offset);                uchar* WIN = win.data;                for (i = 0; i < win_size; i++, start_x++)                {                    int pixel_x = start_x;                    int pixel_y = start_y;                    for (j = 0; j < win_size; j++, pixel_y--)                    {                        int x = MAX(pixel_x, 0);                        int y = MAX(pixel_y, 0);                        x = MIN(x, img->cols - 1);                        y = MIN(y, img->rows - 1);                        WIN[i*win_size + j] = img->at<uchar>(y, x);                    }                }            }            // Scale the window to size PATCH_SZ so each pixel's size is s. This            // makes calculating the gradients with wavelets of size 2s easy            resize(win, _patch, _patch.size(), 0, 0, INTER_AREA);            // Calculate gradients in x and y with wavelets of size 2s            for (i = 0; i < PATCH_SZ; i++)            for (j = 0; j < PATCH_SZ; j++)            {                float dw = DW[i*PATCH_SZ + j]; // 高斯加权系数                float vx = (PATCH[i][j + 1] - PATCH[i][j] + PATCH[i + 1][j + 1] - PATCH[i + 1][j])*dw;                float vy = (PATCH[i + 1][j] - PATCH[i][j] + PATCH[i + 1][j + 1] - PATCH[i][j + 1])*dw;                DX[i][j] = vx;                DY[i][j] = vy;            }            // Construct the descriptor            vec = descriptors->ptr<float>(k);            for (kk = 0; kk < dsize; kk++)                vec[kk] = 0;            double square_mag = 0;            if (extended)            {                // 128维描述子，考虑dx与dy的正负号                for (i = 0; i < 4; i++)                for (j = 0; j < 4; j++)                {                    // 每个方块内是一个5s * 5s的区域，每个方法由8个特征描述                    for (int y = i * 5; y < i * 5 + 5; y++)                    {                        for (int x = j * 5; x < j * 5 + 5; x++)                        {                            float tx = DX[y][x], ty = DY[y][x];                            if (ty >= 0)                            {                                vec[0] += tx;                                vec[1] += (float)fabs(tx);                            }                            else {                                vec[2] += tx;                                vec[3] += (float)fabs(tx);                            }                            if (tx >= 0)                            {                                vec[4] += ty;                                vec[5] += (float)fabs(ty);                            }                            else {                                vec[6] += ty;                                vec[7] += (float)fabs(ty);                            }                        }                    }                    for (kk = 0; kk < 8; kk++)                        square_mag += vec[kk] * vec[kk];                    vec += 8;                }            }            else            {                // 64位描述子                for (i = 0; i < 4; i++)                for (j = 0; j < 4; j++)                {                    for (int y = i * 5; y < i * 5 + 5; y++)                    {                        for (int x = j * 5; x < j * 5 + 5; x++)                        {                            float tx = DX[y][x], ty = DY[y][x];                            vec[0] += tx; vec[1] += ty;                            vec[2] += (float)fabs(tx); vec[3] += (float)fabs(ty);                        }                    }                    for (kk = 0; kk < 4; kk++)                        square_mag += vec[kk] * vec[kk];                    vec += 4;                }            }            // 归一化 描述子 以满足 光照不变性            vec = descriptors->ptr<float>(k);            float scale = (float)(1. / (sqrt(square_mag) + DBL_EPSILON));            for (kk = 0; kk < dsize; kk++)                vec[kk] *= scale;        }    }};

5. 总结

实际上有文献指出，SURF比SIFT工作更出色。他们认为主要是因为SURF在求取描述子特征矢量时，是对一个子块的梯度信息进行求和，而SIFT则是依靠单个像素梯度的方向

0 0