1. SURF比于SIFT
SURF（Speeded Up Robust Features）是对SIFT的改进版本，改进后的主要优点是速度更快，更适合做实时的特征检查。对于需要实时运算的场合，如基于特征点匹配的实时目标跟踪系统，每秒要处理8-24帧的图像，需要在毫秒级内完成特征点的搜索、特征矢量生成、特征矢量匹配、目标锁定等工作，这样SIFT算法就很难适应这种需求了。实验证明，SURF算法较SIFT在运算速度上要快3倍左右。
SURF算法是SIFT算法的加速版，opencv的SURF算法在适中的条件下完成两幅图像中物体的匹配基本实现了实时处理，其快速的基础实际上只有一个——积分图像haar求导。
不论科研还是应用上都希望可以和人类的视觉一样通过程序自动找出两幅图像里面相同的景物，并且建立它们之间的对应，前几年(David Lowe 于1999年首次提出，与2004年又予于完善)才被提出的SIFT（尺度不变特征）算法提供了一种解决方法，通过这个算法可以使得满足一定条件下两幅图像中相同景物的某些点（后面提到的关键点）可以匹配起来，为什么不是每一点都匹配呢？下面的论述将会提到。
SIFT算法实现物体识别主要有三大工序，1、提取关键点；2、对关键点附加详细的信息（局部特征）也就是所谓的描述器；3、通过两方特征点（附带上特征向量的关键点）的两两比较找出相互匹配的若干对特征点，也就建立了景物间的对应关系。
日常的应用中，多数情况是给出一幅包含物体的参考图像，然后在另外一幅同样含有该物体的图像中实现它们的匹配。两幅图像中的物体一般只是旋转和缩放的关系，加上图像的亮度及对比度的不同，这些就是最常见的情形。基于这些条件下要实现物体之间的匹配，SIFT算法的先驱及其发明者想到只要找到多于三对物体间的匹配点就可以通过射影几何的理论建立它们的一一对应。首先在形状上物体既有旋转又有缩小放大的变化，如何找到这样的对应点呢？于是他们的想法是首先找到图像中的一些“稳定点”，这些点是一些十分突出的点不会因光照条件的改变而消失，比如角点、边缘点、暗区域的亮点以及亮区域的暗点，既然两幅图像中有相同的景物，那么使用某种方法分别提取各自的稳定点，这些点之间会有相互对应的匹配点，正是基于这样合理的假设，SIFT算法的基础是稳定点。SIFT算法找稳定点的方法是找灰度图的局部最值，由于数字图像是离散的，想求导和求最值这些操作都是使用滤波器，而滤波器是有尺寸大小的，使用同一尺寸的滤波器对两幅包含有不同尺寸的同一物体的图像求局部最值将有可能出现一方求得最值而另一方却没有的情况，但是容易知道假如物体的尺寸都一致的话它们的局部最值将会相同。SIFT的精妙之处在于采用图像金字塔的方法解决这一问题，我们可以把两幅图像想象成是连续的，分别以它们作为底面作四棱锥，就像金字塔，那么每一个截面与原图像相似，那么两个金字塔中必然会有包含大小一致的物体的无穷个截面，但应用只能是离散的，所以我们只能构造有限层，层数越多当然越好，但处理时间会相应增加，层数太少不行，因为向下采样的截面中可能找不到尺寸大小一致的两个物体的图像。有了图像金字塔就可以对每一层求出局部最值，但是这样的稳定点数目将会十分可观，所以需要使用某种方法抑制去除一部分点，但又使得同一尺度下的稳定点得以保存。有了稳定点之后如何去让程序明白它们之间是物体的同一位置？研究者想到以该点为中心挖出一小块区域，然后找出区域内的某些特征，让这些特征附件在稳定点上，SIFT的又一个精妙之处在于稳定点附加上特征向量之后就像一个根系发达的树根一样牢牢的抓住它的“土地”，使之成为更稳固的特征点，但是问题又来了，遇到旋转的情况怎么办？发明者的解决方法是找一个“主方向”然后以它看齐，就可以知道两个物体的旋转夹角了。下面就讨论一下SIFT算法的缺陷。
SIFT/SURT采用henssian矩阵获取图像局部最值还是十分稳定的，但是在求主方向阶段太过于依赖局部区域像素的梯度方向，有可能使得找到的主方向不准确，后面的特征向量提取以及匹配都严重依赖于主方向，即使不大偏差角度也可以造成后面特征匹配的放大误差，从而匹配不成功；另外图像金字塔的层取得不足够紧密也会使得尺度有误差，后面的特征向量提取同样依赖相应的尺度，发明者在这个问题上的折中解决方法是取适量的层然后进行插值。SIFT是一种只利用到灰度性质的算法，忽略了色彩信息，后面又出现了几种据说比SIFT更稳定的描述器其中一些利用到了色彩信息，让我们拭目以待。
最后要提一下，我们知道同样的景物在不同的照片中可能出现不同的形状、大小、角度、亮度，甚至扭曲；计算机视觉的知识表明通过光学镜头获取的图像，对于平面形状的两个物体它们之间可以建立射影对应，对于像人脸这种曲面物体，在不同角度、距离，不同相机参数下获取的两幅图像，它们之间不是一个线性对应关系，就是说我们即使获得两张图像中的脸上若干匹配好的点对，还是无法从中推导出其他点的对应。
2.SURF特征提取的原理详解
这一部分将对SURF的实现原理进行详细的讲解。
上文做了说明，我们用像素点的海塞矩阵来描述一副图像的关键特征点。下面就来讲解图像的海塞矩阵在图像处理中如何计算，海塞矩阵就是求X，Y方向的二阶偏导。
公式2.1
以上公式就是对每一个像素点求出一个海赛矩阵，对于X方向和Y方向的二阶偏导。特征点要具备尺度无关性（这样做的目的就是进一步筛选那些关键点，在不同尺度下都达到要求的点，才列为关键点），所以在计算海赛矩阵以前，要对图像进行高斯滤波，得到不同尺度的图像。并且计算不同尺度下的图像特征点。那什么是图像的尺度？通俗的讲，一副图像在不同远近的情况下，清晰程度不同，越远越模糊，越近越清楚。在图像处理领域我们可以通过高斯滤波来模拟不同远近条件下图像的清晰程度不同，来模拟人眼在不同远近对物体的成像。
不同尺度的图像得到，通过高斯滤波，高斯滤波在二维空间的定义为：
公式2.2
图2.1 matlab中图像表示
在二维空间中，这个公式生成的曲面的等高线是从中心开始呈正态分布的同心圆。分布不为零的像素组成的卷积矩阵与原始图像做变换。每个像素的值都是周围相邻像素值的加权平均。原始像素的值有最大的高斯分布值，所以有最大的权重，相邻像素随着距离原始像素越来越远，其权重也越来越小。这样进行模糊处理比其它的均衡模糊滤波器更高地保留了边缘效果。其中，Sigma是正态分布的标准差，值越大，图像越模糊(平滑)。如二维模板大小为m*n，则模板上的元素(x,y)对应的高斯计算公式为：
公式2.3
通过以上公式我们可以计算出不同σ值下，不同大小的高斯模板。下图是5*5的高斯模板卷积计算示意图。高斯模板是中心对称的。
图2.2
图2.3
以上图片代表不同尺度因子σ取不同值的结果。由此得到不同尺度的图像，如图2.3
图像尺度不同实际就是图像的模糊程度不同，就像一副图像在人眼中表示远近不同，越近越清楚，越远越模糊。所以用高斯滤波可以很好的模拟不同尺度的同一幅图像。
图2.4
如图2.4所示，和SIFT算法不同的是SURF算法采用不同的高斯模板大小得到不同尺度的图像，一共分为四组，每组四个不同尺度的图像。就像最上面一层第一副图像的高斯模板为9X9,第二幅为15X15….以此类推。下面说明其中一副图像如何计算得到，其他图层就是采用不同大小的高斯模板，也可一样得到。
经过高斯滤波的海赛矩阵公式如下：
公式2.5
由于求海赛矩阵时先对图片进行高斯平滑，而后求二阶偏导，在图像处理中一般使用模板卷积形成。模板就是数学公式求值在图像上针对离散像素点的简化，模板在程序中就是一个矩阵。矩阵的成员数值通过数学计算得到。为了便于计算海赛矩阵行列式，还引入了积分图像的概念。对于一幅灰度的图像，积分图像中的任意一点(x,y)的值是指从图像的左上角到这个点的所构成的矩形区域内所有的点的灰度值之和。下面介绍积分图如何计算：
图2.4(灰度原图) 图2.5(积分图)
用公式表示就是：
图2.6
区域1 : = sum(A);
区域2 : = sum(A + B);
区域3 : = sum(A + C);
区域4 : = sum(A + B + C + D);
所以，如果需要计算D区域中的灰度和，则
sum(D) = 区域4 - 区域2 - 区域3 + 区域1 （都是灰度值）。
很明显，这里仅仅只需要通过查表得到 1、2、3、4点的积分图像的值即可得到。
下面给出积分图像的C语言程序片段：
void integral_matrix_Old(unsigned char* srcImage,unsigned int* sumMatrix,int width,int height,int line_length )
{int y,x,k;
int offset = line_length - width;
int srcstep = width;
int s = 0;//sum的第0列的值都是0
int sumstep = width+1;
int cn = 1;
memset( sumMatrix, 0,(width+1)*sizeof(sumMatrix[0]));/*将sum的第一行设为0*/
sumMatrix += sumstep + cn;
for( y = 0; y < height; y++, srcImage += srcstep - cn+offset, sumMatrix += sumstep - cn )
{
for( k = 0; k < cn; k++, srcImage++, sumMatrix++ )
{
s = sumMatrix[-cn] = 0;//sum的第0列的值都是0
for( x = 0; x <width; x += cn )
{
s += srcImage[x];
sumMatrix[x] = sumMatrix[x - sumstep] + s;//积分图像
}
}
}
}在计算图像特征点的时候，surf算法是采用计算原图中每个像素的的海赛矩阵行列式的近似值构成，其行列式计算公式如下：
公式2.4
其中0.9是原算法作者给出的经验值，有一套理论计算，具体请看原文论文。算法是在不同尺度图像下进行求海赛矩阵的，不同尺度图像通过高斯平滑得到，然后求不同方向的二阶偏导。在离散图像像素中，用一个图像处理模板代替就可以了。比如X方向，Y方向，XY方向上的模板如图2.7，我们通过X方向的模板求出Dxx，通过Y方向模板求出Dyy,通过XY方向模板求出Dxy. 最后分别求出之后带入公式2.4就可以得到每个像素海塞矩阵行列式的值。
图2.7
图2.7分别表示在X方向，Y方向，XY方向的运算模板。
下面就Y方向的滤波器详细说明这些箱式滤波器的设计原理
图2.8
上图的左边即用高斯平滑然后在y方向上求二阶导数的模板，为了加快运算用了近似处理，其处理结果如右图所示，这样就简化了很多。请看简化的区域有两块2X9的灰色区域，两块3X5的白色区域，一块3X5的黑色区域组成。这个9X9的模板被划分成5块不同区域。对于高斯滤波我们知道，对于一个像素的高斯滤波计算，是对周围9X9区域的像素的加权平均得到。区域内不同位置的像素，权值不一样，权值的意义是当前点像素的值对最终滤波值的影响程度。且看简化图像五个不同区域的权值为白色区域为1,黑色区域为-2，灰色区域为0。我们规定了不同区域权值大小。那么这个区域的像素灰度值乘以权值，就是这个区域对最终高斯滤波值的贡献。那么这个区域的灰度值如何得到，前面计算了积分图像就很容易得到这个区域的灰度值的总大小。如此可以看出极大的简化了计算。其他方向的原理类似，这里不做解释。注意1,-2,1的设计，这个权值的设计让我们完成了二阶偏导和高斯滤波的计算。图像中离散数据的求导运算，简化为减法运算。
其滤波器初始数据如下：
const int dx_s[3][5] = { {0, 2, 3, 7, 1}, {3, 2, 6, 7, -2}, {6, 2, 9, 7, 1} };
const int dy_s[3][5] = { {2, 0, 7, 3, 1}, {2, 3, 7, 6, -2}, {2, 6, 7, 9, 1} };
const int dxy_s[4][5] = { {1, 1, 4, 4, 1}, {5, 1, 8, 4, -1}, {1, 5, 4, 8, -1}, {5, 5, 8, 8, 1} };
void ResizeHaarPattern( const int* src, SURF_HARRI_FILTER* dst, unsigned int n, unsigned int oldSize, unsigned int newSize, unsigned int widthStep )
{
float ratio = (float)newSize/oldSize;
unsigned int k=0;
for( k = 0; k < n; k++ )
{
int dx1 =(int) ( ratio*(src+k*5)[0]+0.5f);/*+0.5主要是为了四舍五入*/
int dy1 = (int)( ratio*(src+k*5)[1]+0.5f);
int dx2 =(int) ( ratio*(src+k*5)[2]+0.5f );
int dy2 = (int)( ratio*(src+k*5)[3]+0.5f);
dst[k].p0 = dy1*widthStep + dx1;//第一个矩形相对于图像原点的坐标
dst[k].p1 = dy2*widthStep + dx1;//第二个矩形相对于图像原点的坐标
dst[k].p2 = dy1*widthStep + dx2;//第三个矩形相对于图像原点的坐标
dst[k].p3 = dy2*widthStep + dx2;//第四个矩形相对于图像原点的坐标
dst[k].w = (src+k*5)[4]/((float)(dx2-dx1)*(dy2-dy1));//这个矩形内的数据所占的权重
}
}上面的代码计算得到那些2X9，3X5的矩形区域相对当前像素点的坐标。知道这些子区域的位置坐标乘以他们所代表的权重，就是这块区域对当前像素点的贡献。下面举例计算在Y方向上的矩形区域，得到Dyy的过程：
分别设d1_y ,d2_y,d3_y为三个不同区域的像素灰度权值和，具体如图所示：
图2.9
d1_y=(s_ptr[pDy_0->p0]+s_ptr[pDy_0->p3] - s_ptr[pDy_0->p1]-s_ptr[pDy_0->p2])*pDy_0->w;
d2_y=(s_ptr[pDy_1->p0]+s_ptr[pDy_1->p3] - s_ptr[pDy_1->p1] - s_ptr[pDy_1->p2])*pDy_1->w;
d3_y=(s_ptr[pDy_2->p0]+s_ptr[pDy_2->p3] - s_ptr[pDy_2->p1] - s_ptr[pDy_2->p2])*pDy_2->w;
看上面一段程序s_ptr[pDy_0->p0]+s_ptr[pDy_0->p3] - s_ptr[pDy_0->p1]-s_ptr[pDy_0->p2]这就代d1_y这个矩形区域里面的灰度值的和是多少，就是根据前面计算的积分图得到的。那么DYY= d1_y + d2_y + d3_y,同理我们可以计算DXX，DXY，那么根据公式2.4，我们就可以计算出经过高斯滤波的每个像素点的海塞矩阵行列式的值大小。我们设定一个阈值来挑选大于某一行列式值的像素。把大于某一值的像素作为一个关键点的候选点。
图2.10
拿出这些关键点的候选点，在两个相邻尺度的的海赛矩阵行列式集合上的26个点做比较。注意因为我们设定一组由四幅尺度不同的图像，那么注定第一幅和最后一幅图像向上或者向下缺少一个邻域，这里我们不去计算即可。如果这个候选关键点的海赛矩阵行列式的值确实大于其他26个点的海赛矩阵行列式的值，那么这个关键点就被列为关键点。从这里我们得出一个这个算法的原理，对于关键点就是一副图像上的角点，边缘点，暗区域的亮点，亮区域的暗点，这些点其实都对应一条数学原理，就是这些点的二阶偏导较大。说白了就是图像在这个区域灰度值变化比较激烈。至于为何要进行高斯模糊，其实得到的是不同尺度下的图像，这样相邻的尺度比较得到的关键点，我想更为稳定一些，有效的去除那些干扰点。以上就是针对一副图像怎么计算关键点。两幅一样的图像或者一副图像包涵另一幅图像的情况，总会有位置一样的关键件。如果只在考虑图像的水平方向或者垂直方向，从这些关键点中匹配到目标图像的位置，也是可以了。但是我们要考虑，目标位置会发生旋转，角度偏移的情况。那么我们就给这些关键点，找到他们的主方向，让这些关键点有大小有方向，那么就可以应对目标位置发生角度偏移的情况了。
三 SURF算法关键点的描述符
3.1 什么是关键特征点的主方向
图3.1
原算法中规定，以前面得到关键点为中心，计算半径为6
（
为特征点所在的尺度值）的邻域内的点在x,y方向的Haar小波（Haar小波边长取4）响应，并给这些响应值赋予高斯权重系数，使得靠近特征点的响应贡献大，而远离特征点的响应贡献小。在360度范围内，每次扫描一个60度的扇区，把这个扇区内的Haar小波响应相加形成新的矢量，选择最长的方向为该特征点的主方向（如图3.1）。
注：Harr小波特征：Haar特征分为三类：边缘特征、线性特征、中心特征和对角线特征，组合成特征模板。特征模板内有白色和黑色两种矩形，并定义该模板的特征值为白色矩形像素的和减去黑色矩形像素的和。Haar特征值反映了图像的灰度变化情况。
图3.2
所示为哈尔小波滤波器在x方向和y方向上描述，黑色部分的权重是-1，白色部分的权重是1。计算出图像在哈尔小波的x和y方向上的响应值之后，对两个值进行因子为的高斯加权，加权后的值分别表示在水平和垂直方向上的方向分量。Harr特征值反应了图像灰度变化的情况，那么这个主方向就是描述那些灰度变化特别剧烈的区域方向。就像图3.1第三幅图像的样子。
3.2 构造关键点的主方向的描述子
我们已经得到了关键点的主方向了，然后以这个关键点为中心，画一个以边长为20S（S为该特征点所在的尺度），以该关键点的主方向为转角的正方形。如图：
图 3.3
V =
图3.4
然后把这些方框分为16个子区域，每个区域统计25个像素的水平方向和垂直方向的Harr小波特征，当然这里的水平和垂直方向都是相对主方向而言。该haar小波特征为水平方向值之和，水平方向绝对值之和，垂直方向之和，垂直方向绝对值之和以及分别的负方向之和。所以每个特征点就是16*8=128维的向量.
3.3 关键点的匹配
当搜索区域和搜索目标的关键点的主方向和描述子都计算出来之后，下一步进行关键点匹配，就是说对于目标图片上的任意关键点A，在搜索区域中找到最相似的点，与之对应。每一个关键点都有一个描述向量，我们采用关键点描述向量的欧式距离来作为两幅图像中关键点的相似性判定度量。取图像1中的某个关键点，并找出其与图像2中欧式距离最近的前两个关键点，在这两个关键点中，如果最近的距离除以次近的距离少于某个比例阈值，则接受这一对匹配点。调整这一比例阈值可以应对不同的匹配情况。这样一来我们可以得到，很多个关键点。
3.4区域匹配和旋转角度计算
上面的计算我们就得到了，各个关键点的特征向量。在工业检测里面，我们必须精确跟踪目标区域的位置，下面我讲解如何通过得到的关键点计算出目标区域的矩形坐标。
图3.5
如图3.5所示我们要在search的区域，跟踪object区域。方便理解我们构造了坐标系。
设object’为目标区域，search为搜索区域。经过前面的计算，我们得到了在object’里面的i,j,k和serach里面的is,js,ks这几个对应的关键点。那么我们要通过这些个对应的关键点，来在search计算出object的目标位置。通过这三个点，我们实际上构造出了，一对三角形。如果在object’里面的i,j,k关键点serach里面的is,js,ks对应的位置就是图像中对应的几何位置的话（就像i点对应object中M字母的一个角点，is点对应search中M字母的一个角点一样）.那么这六个点所构造的两个三角形，就是相似三角形。我们可以通过几何位置关系求出他们各个对应边的比值。在意图图像中的关键点有很多这样的对应的三角形。我们设定一个阈值，如果是完全相似的话，各个边对应的比值是一个固定的值，如果三角形一模一样比值就是1.在图像处理中考虑计算误差等图像干扰。不可能比值完全相等的。我们寻找多个这样的相对三角形，做计算，如果比例对应关系越好说明，我们寻找的这三个关键点的精确度越高。我们就把这三个点作为最佳匹配点。然后我们取这两个三角形的两条对应边，我们计算这两个边的角度，就可以得到目标区域的转角，就可以跟踪如果目标区域带角度偏移的情况。
以上就是整个算法的介绍，如有描述错误欢迎批评指正。
3. opencv API调用
在opencv中SURF的API调用和SIFT的调用时类似的，可以参考 OpenCV3中的SURF特征点的寻找和匹配
如果用的是opencv2.xx那么需要用不同的API，可以参考SIFT特征提取分析及Opencv API调用。
匹配效果如下：
#include<opencv.hpp>
#include<iostream>
#include<xfeatures2d.hpp>
using namespace cv;
using namespace std;
using namespace xfeatures2d;
int main(void)
{
Mat object = imread("d:/Opencv Picture/Match/box.png", 0);
Mat object_scene = imread("d:/Opencv Picture/Match/box_in_scene.png", 0);
//extract feature
int minHessian = 400;
Ptr<SURF> surfPtr = SURF::create(minHessian);
vector<KeyPoint> keypoint_object, keypoint_scene;
Mat descriptor_object, descriptor_scene;
surfPtr->detectAndCompute(object, Mat(), keypoint_object, descriptor_object);
surfPtr->detectAndCompute(object_scene, Mat(), keypoint_scene, descriptor_scene);
//matching
FlannBasedMatcher matcher;
vector<DMatch> matches;
matcher.match(descriptor_object, descriptor_scene, matches);
//find good matchers
double minDist = 1000;
double maxDist = 0;
for (int i = 0; i < matches.size(); ++i)
{
if (matches[i].distance > maxDist)
{
maxDist = matches[i].distance;
}
if (matches[i].distance < minDist)
{
minDist = matches[i].distance;
}
}
vector<DMatch> goodMatches;
for (int i = 0; i < matches.size(); ++i)
{
if (matches[i].distance < max(3 * minDist, 0.02))
{
goodMatches.push_back(matches[i]);
}
}
//draw mathes
Mat matchImage;
drawMatches(object, keypoint_object, object_scene, keypoint_scene, goodMatches, matchImage, Scalar::all(-1), Scalar::all(-1), vector<char>(), DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS);
vector<Point2f> obj;
vector<Point2f> obj_scene;
for (size_t t = 0; t < goodMatches.size(); ++t)
{
obj.push_back(keypoint_object[goodMatches[t].queryIdx].pt);
obj_scene.push_back(keypoint_scene[goodMatches[t].trainIdx].pt);
}
Mat H = findHomography(obj, obj_scene, RANSAC);
vector<Point2f> obj_corner(4);
vector<Point2f> scene_corner(4);
obj_corner[0] = Point(0, 0);
obj_corner[1] = Point(object.cols,0);
obj_corner[2] = Point(object.cols, object.rows);
obj_corner[3] = Point(0, object.rows);
perspectiveTransform(obj_corner, scene_corner, H);
//draw line
line(matchImage, scene_corner[0] + Point2f(object.cols, 0), scene_corner[1] + Point2f(object.cols, 0),Scalar(0,0,255),2,8,0);
line(matchImage, scene_corner[1] + Point2f(object.cols, 0), scene_corner[2] + Point2f(object.cols, 0), Scalar(0, 0, 255), 2, 8, 0);
line(matchImage, scene_corner[2] + Point2f(object.cols, 0), scene_corner[3] + Point2f(object.cols, 0), Scalar(0, 0, 255), 2, 8, 0);
line(matchImage, scene_corner[3] + Point2f(object.cols, 0), scene_corner[0] + Point2f(object.cols, 0), Scalar(0, 0, 255), 2, 8, 0);
imshow("matchImage", matchImage);
waitKey();
return 0;
}
