图像处理——SURF算法

首先感谢以下3位的渊博知识，帮助我理解SURF算法：

（1）☆Ronny丶  https://www.cnblogs.com/ronny/p/4045979.html

 

 

 

（2）-牧野-   http://blog.youkuaiyun.com/dcrmg/article/details/52601010

 

（2）-牧野-  http://blog.youkuaiyun.com/dcrmg/article/details/52601010

（3）PeOS http://blog.youkuaiyun.com/cxp2205455256/article/details/41311013

 

Speeded Up Robust Features（SURF，加速稳健特征），是一种稳健的局部特征点检测和描述算法。最初由Herbert Bay发表在2006年的欧洲计算机视觉国际会议（Europen Conference on Computer Vision，ECCV）上，并在2008年正式发表在Computer Vision and Image Understanding期刊上。

Surf是对David Lowe在1999年提出的Sift算法的改进，提升了算法的执行效率，为算法在实时计算机视觉系统中应用提供了可能。与Sift算法一样，Surf算法的基本路程可以分为三大部分：局部特征点的提取、特征点的描述、特征点的匹配。

但Surf在执行效率上有两大制胜法宝——一个是积分图在Hessian（黑塞矩阵）上的使用，一个是降维的特征描述子的使用。了解这两大法宝是如何出奇制胜前，先回顾一下传统Sift算法的基本思路及其优缺点。

一、Sift算法回顾

 

Sift算法的三大工序为，(1)提取关键点；(2)对关键点附加详细的信息（局部特征）也就是所谓的描述器；(3)通过两方特征点（附带上特征向量的关键点）的两两比较找出相互匹配的若干对特征点，也就建立了景物间的对应关系。提取关键点和对关键点附加详细的信息（局部特征）也就是所谓的描述器可以称做是Sift特征的生成，即从多幅图像中提取对尺度缩放、旋转、亮度变化无关的特征向量，Sift特征的生成一般包括以下几个步骤：

(1)、构建尺度空间，检测极值点，获得尺度不变性；

 

(2)、特征点过滤并进行精确定位；

 

(3)、为特征点分配方向值；

 

(4)、生成特征描述子。

以特征点为中心取16*16的邻域作为采样窗口，将采样点与特征点的相对方向通过高斯加权后归入包含8个bin的方向直方图，最后获得4*4*8的128维特征描述子。示意图如下：

 

当两幅图像的Sift特征向量生成以后，下一步就可以采用关键点特征向量的欧式距离来作为两幅图像中关键点的相似性判定度量。取图1的某个关键点，通过遍历找到图像2中的距离最近的两个关键点。在这两个关键点中，如果次近距离除以最近距离小于某个阙值，则判定为一对匹配点。

二、SURF算法原理

（1）积分图像

积分图像的概念是由Viola和Jones提出的。积分图像中任意一点（i，j）的值为原图像左上角到任意点（i，j）相应的对焦区域的灰度值的总和，其数学公式如图1所示:

那么，当我们想要计算图片一个区域的积分，就只需计算这个区域的四个顶点在积分图像里的值，便可以通过2步加法和2步减法计算得出，其数学公式如下：

 

 

（2）斑点检测

    斑点：与周围有着颜色和灰度差别的区域。
    在一个一维信号中，让它和高斯二阶导数进行卷积，也就是拉普拉斯变换，那么在信号的边缘处就会出现过零点，如下图所示：

    高斯拉普拉斯Log探测器的响应值就是在衡量图像的相似性，如下图是一个图像的高斯拉普拉斯变换的三维图和灰度图显示，在图像中的斑点尺寸与高斯拉普拉斯函数
的形状趋于一致时，图像的拉普拉斯响应抵达最大。

    Hession矩阵就是利用二阶微分来进行斑点检测，其矩阵如下：

（3）hession矩阵行列式的简化

 

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

 

H(x,σ)=[Lxx(x,σ)Lxy(x,σ)Lxy(x,σ)Lyy(x,σ)]H(x,σ)=[Lxx(x,σ)Lxy(x,σ)Lxy(x,σ)Lyy(x,σ)]

 

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

 

当Hessian矩阵的判别式取得局部极大值时，判定当前点是比周围邻域内其他点更亮或更暗的点，由此来定位关键点的位置。

我们知道在离散数字图像中，一阶导数是相邻像素的灰度差：

 

 

二阶导数是对一阶导数的再次求导：

 

 

反过来看Hessian矩阵的判别式，其实就是当前点对瓶方向二阶偏导乘以垂直方向的二阶偏导再减去当前点水平、垂直二阶偏导的二次方：

 

 

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

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

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

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

 

Det(H)=LxxLyy–LxyLxy=DxxLxxDxxDyyLyyDyy−DxyLxyDxyDxyLxyDxy≈DxxDyy–(wDxy)2Det(H)=LxxLyy–LxyLxy=DxxLxxDxxDyyLyyDyy−DxyLxyDxyDxyLxyDxy≈DxxDyy–(wDxy)2

 

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

 

w=|Lxy(σ)|F|Dxx(σ)F||Lxx(σ)|F|Dxy(σ)F|=0.912w=|Lxy(σ)|F|Dxx(σ)F||Lxx(σ)|F|Dxy(σ)F|=0.912

 

其中|X|F|X|F为Frobenius范数。另外，响应值还要根据滤波器大小进行归一化处理，以保证任意大小滤波器的F范数是统一的。0.9^2是滤波器响应的相关权重w是为了平衡Hessian行列式的表示式。这是为了保持高斯核与近似高斯核的一致性。理论上来说对于不同的σ的值和对应尺寸的模板尺寸，w值是不同的，但为了简化起见，可以认为它是同一个常数。 使用近似的Hessian矩阵行列式来表示图像中某一点x处的斑点响应值，遍历图像中所有的像元点，便形成了在某一尺度下琉璃点检测的响应图像。使用不同的模板尺寸，便形成了多尺度斑点响应的金字塔图像，利用这一金字塔图像，就可以进行斑点响应极值点的搜索。

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

 

（4）尺度金字塔构造

    相比于sift算法的高斯金字塔构造过程，sift算法速度有所提高。在sift算法中，每一组（octave）的图像大小是不一样的，下一组是上一组图像的降采样（1/4大小）；在每一组里面的几幅图像中，他们的大小是一样的，不同的是他们采用的尺度σ不同。而且在模糊的过程中，他们的高斯模板大小总是不变的，只是尺度σ改变。对于surf算法，图像的大小总是不变的，改变的只是高斯模糊模板的尺寸，当然，尺度σ也是在改变的。 

 

在SURF中，采用不断增大盒子滤波器模板尺寸与积分图像求取Hession矩阵响应，然后在响应图像上采用3D非极大值抑制，求取各种不同尺度

的斑点，以下是两种不同的金字塔，SURF的金字塔属于第二种：

上图中a为高斯模板保持不变，图像大小改变的情况，适用于sift算法，图b是高斯模板改变，图像大小保持不变的情况，适用于surf算法。因为surf算法没有了降采样的过程，因此处理速度得到提高。 

 

    SURF中采用9X9尺寸的滤波器作为起始滤波器，之后的滤波器尺寸可由以下公式计算得出：

 

octave、interval在公式中都是从1开始，也就是当第0组第0层时，在公式中octave = 1， interval = 1。

   

滤波器响应长度l（小写L）、滤波器尺寸L、组索引O、层索引S、尺度sigma之间的关系如下：

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

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

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

 

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

        

采用类似的方法来处理其他几组的模板序列。其方法是将滤波器尺寸增加量翻倍（6，12，24，38）（这句话不太明白，不过将本段绿色数字代入公式推导，得出其意思是：第0组相邻层之间像素增加2l0，两边各增加l0；第1组相邻层之间像素增加4l0，两边各增加2l0；第2组相邻层之间像素增加8l0，两边各增加4l0）。这样，可以得到第二组的滤波器尺寸，它们分别为15，27，39，51。第三组的滤波器尺寸为27，51，75，99。如果原始图像的尺寸仍然大于对应的滤波器尺寸，尺度空间的分析还可以进行第四组，其对应的模板尺寸分别为51，99，147和195。下图显示了第一组至第三组的滤波器尺寸变化。

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

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

（5）非最大值抑制

所有小于预设极值的取值都被丢弃，增加极值使检测到的特征点数量减少，最终只有几个特征最强点会被检测出来。检测过程中使用与该尺度层图像解析度相对应大小的滤波器进行检测，以3×3的滤波器为例，该尺度层图像中9个像素点之一图2检测特征点与自身尺度层中其余8个点和在其之上及之下的两个尺度层9个点进行比较，共26个点，图中标记‘x’的像素点的特征值若大于周围像素则可确定该点为该区域的特征点。

 

 （6）局部极大值精确定位

 

    跟sift算法类似，采用3维线性插值法得到亚像素级的特征点，同时也去掉那些值小于一定阈值的点，增加极值使检测到的特征点数量减少，最终只有几个特征最强点会被检测出来。

（7）特征点主方向分配

Sift特征点主方向分配是采用在特征点邻域内统计其梯度直方图，取直方图bin值最大的以及超过最大bin值80%的那些方向作为特征点的主方向。而在Surf中，采用的是统计特征点圆形邻域内的harr小波特征。

以特征点为中心，计算半径为6s(S为特征点所在的尺度值)的邻域内的点在x、y方向的Haar小波(Haar小波边长取4s)响应，Harr小波

模板如图所示：（ps:关于Haar小波，要完全理解的话要看不少资料，我个人在这里把它姑且视为x、y方向的梯度变化，不知是否正确）

 

    计算出图像在哈尔小波的x和y方向上的响应值之后，对两个值进行因子为2S的高斯加权，加权后的值分别表示在水平和垂直方向上的方向分量。

Harr特征值反应了图像灰度变化的情况，那么这个主方向就是描述那些灰度变化特别剧烈的区域方向。

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

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

 

mw=∑wdx+∑wdymw=∑wdx+∑wdy

 

 

θw=arctan(∑wdx/∑wdy)θw=arctan(∑wdx/∑wdy)

 

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

 

θ=θw|max{mw}θ=θw|max{mw}

 

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

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

（8）特征点特征矢量生成

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

图2 特征描述子表示

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

 

V子块=[∑dx,∑|dx|,∑dy,∑|dy|]V子块=[∑dx,∑|dx|,∑dy,∑|dy|]

 

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

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

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

 

x=x0–j×scale×sin(θ)+i×scale×cos(θ)x=x0–j×scale×sin(θ)+i×scale×cos(θ)

 

 

y=y0–j×scale×cos(θ)+i×scale×sin(θ)y=y0–j×scale×cos(θ)+i×scale×sin(θ)

 

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

 

dx′=w(−dx×sin(θ)+dy×cos(θ))dx′=w(−dx×sin(θ)+dy×cos(θ))

 

 

dy′=w(−dx×cos(θ)+dy×sin(θ))dy′=w(−dx×cos(θ)+dy×sin(θ))

 

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

 

（9）匹配

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

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

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

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