图像特征提取（三）——SURF特征提取算法详解

矫情吗；* 2022-05-20 00:17 1260阅读 0赞

本文大量参考了:http://blog.csdn.net/cxp2205455256/article/details/41311013;并且在此博客的基础之上，增加了一些新的内容。

一、SURF算法的概述

SURF，英文的全称为Speed Up Robust Features,直译为:加速版的具有鲁棒特性的特征算法，是由Bay在2006年首次提出的。该算法对经典的尺度不变特征变换算法(SIFT算法)进行了改进，SIFT算法最大的缺点就是如果不借助硬件或者专门的图像处理器进行加速的话，SIFT算法很难达到实时处理的效果。SURF算法最大的特征在于采用了Haar特征以及积分图像的概念，这大大的加速了程序的运行时间。SURF算法不仅保持了SIFT算法的尺度不变和旋转不变的特性，而且对光照变化和放射变化同样具有很强的鲁棒性。SURF算法一般应用于计算机视觉中的物体识别、图像拼接、图像配准以及3D重建中。

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

![Center][]

![Center 1][]

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

![Center 2][]

三、SURF算法的原理

1、斑点检测

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

![Center 3][]

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

![Center 4][]

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

![Center 5][]

2、构建Hessian矩阵构建高斯图像金字塔

其实SURF构造的金字塔图像与SIFT构造的金字塔图像有很大的不同，正是因为这些不同才加快了SURF算法关键点的检测速度。SIFT采用的是DOG(高斯差分图像)图像，而SURF采用的是Hessian矩阵行列式近似值图像。Hessian矩阵式SURF算法的核心。在数学中，Hessian矩阵是一个自变量为向量的实值函数的二阶偏导数组成的方块矩阵。为了方便运算，假设函数f(x,y)，Hessian矩阵H是函数的二阶偏导数组成的方块矩阵。首先来看看图像中某个像素点的Hessian矩阵，如下图所示:

![Center 5][]

即，我们在图像中的每一个像素点，都可以求出一个Hessian矩阵，然后，我们再根据Hessian矩阵的判别式det(H)的正负号，来判别该点是不是极值点，Hessian矩阵的判别式如下所示:

![Center 6][]

但是，由于特征点需要具备尺度无关性，所在我们在Hessian矩阵的构造前，需要首先对图像进行高斯平滑滤波。这样，在经过高斯平滑滤波之后，在对图像中的每一个像素点进行Hessian矩阵的求解，具体的计算公式图下所示:

L(x,y) = G(t)\*I(x,t)

其中L(x,t)表示的是一幅图像在不同解析度下的表示，其本质表示的就是原始图像在不同尺度高斯核平滑之后的高斯平滑图像，可以利用高斯核G(t)和图像函数I(x)在点x的卷积来实现，其中的高斯核G(t)的计算公式为：高斯函数g(x)对变量x的二阶导数。

为了方便应用，Herbert Bay提出用近似值代替L(x,t)。为了平衡准确值与近似值之间的误差引入了权值，权值随着尺度的变化而变化，此时，Hessian矩阵的判别公式如下所示:

![Center 7][]

由于求Hessian时，首先要先进行高斯平滑，然后求二阶导数，这在离散的像素点中应用模板的卷积就可以形成，那么，我们就在想能不能用一个模板就可以将高斯平滑过程中的高斯平滑模板和在求图像二阶导数过程中的拉普拉斯算子的二阶微分模板代替，答案是：这是可以的！然后，作者就提出了一种新的模板代替了这两宗模板的组合。比如说Y方向上的模板就可以用下图中的第二行的第二个模板代替；X方向上的模板就可以用第二行的第一个模板代替；X和Y方向上的二阶混合偏导就可以用第一行中的第三个模板代替。其实，这就是Haar特征中的:积分图+边缘模板、对角模板、中心环绕模板。至此，我们就用一个积分图+模板的方式代替了高斯平滑和图像的二阶求导的过程。具体的模板可以参考如下所示：

![Center 8][]

上面讲了这么多，其实就只为了得到一张近似Hessian矩阵的行列式图像，SURF算法中的Hessian矩阵行列式的图像就相当于SIFT算法中的高斯差分图像，即DOG图像。但是，在金字塔图像中分为很多层，每一层叫一个Octave，每一个Octave中，又有几张尺度不同的图片(解析度不同，即高斯核函数的sigma不同)。在SIFT算法中，同一个Octave层中的图片尺寸相同，但是图片的解析度(即图片的模糊程度)不同。而不同的Octave层中的图片尺寸大小也不相同，因为它是由上一层图片降采样得到的。在进行高斯模糊时，SIFT的高斯模板的大小是始终不变的，只是在不同的Octave之间，图片的大小不一样。而在SURF中，图片的大小是一直不变的，不同Octave层待检测图片是改变了高斯模糊尺寸的大小得到的，当然了，同一个Octave中不同图片用到的高斯模板的尺度(模糊程度)也是不同的。算法允许尺度空间多层图像同时被处理，不需要对图像进行二次抽样，从而提高了算法的性能。当然在求解Hessian矩阵行列式图像的过程中，积分图像的使用，也是SURF算法性能突飞猛进的原因之一。

四、利用非极大值抑制初步确定关键点的位置(特征点)

1、尺度金字塔构造

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

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

![Center 9][]

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

![Center 10][]

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

采用这种方式来定义滤波器尺寸的理由如下;

![Center 11][]

滤波器响应长度、滤波器尺寸、组索引O、层索引S

尺度sigma之间的关系如下：

![Center 12][]

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

![Center 13][]

在通常尺度分析情况下，随着尺度的增大，被检测到的斑点数量迅速衰减。所以一般进行3-4组就可以了，与此同时，为了减少运算量，提高计算的速度，可以考虑在滤波时，将采样间隔设为2。  
    2、为了在图像及不同尺寸中定位兴趣点，我们用了3×3×3邻域非最大值抑制：

此步骤和SIFT算法类似，将经过Hessian矩阵处理过的每个像素点与其三维领域的26个点的大小进行比较，如果它是这26个点中的最大值或者最小值，则保留下来，当做初步确定的关键点(特征点)。检测过程中，使用与该尺度图像解析度相对应大小的滤波器进行检测。以3\*3的滤波器为例，该尺度层图像中的9个像素点之一的检测特征点与自身的8个特征点和其上下两个尺度层中的各9个点进行比较，共计与26个点比较。

![Center 14][]

3、局部极大值精确定位

这里也和SIFT算法类似，采用三维线性插值算法得到亚像素级别的特征点(关键点)，同时也去掉小于一定阈值的点，增加极值是检测到的特征点数量减少，最终只有几个特征最强的点被检测了出来。

五、特征点描述符

1、特征点方向分配

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

模板如图所示：

![Center 15][]

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

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

接着，以特征点为中心，张角为π/3的扇形滑动，计算窗口内的Harr小波响应值dx、dy的累加：

![Center 16][]

扇形窗口的滑动如图所示：

![Center 17][]

在OpenSURF的C++代码实现如下：

\[cpp\]  [view plain][]   [copy][view plain]

1.  **for**(**int** i = -6; i <= 6; ++i)   
2.  \{  
3.    **for**(**int** j = -6; j <= 6; ++j)   
4.    \{  
5.      **if**(i\*i + j\*j < 36)   
6.      \{  
7.        gauss = **static\_cast**<**float**>(gauss25\[id\[i+6\]\]\[id\[j+6\]\]);  
8.        resX\[idx\] = gauss \* haarX(r+j\*s, c+i\*s, 4\*s);  
9.        resY\[idx\] = gauss \* haarY(r+j\*s, c+i\*s, 4\*s);  
10.       Ang\[idx\] = getAngle(resX\[idx\], resY\[idx\]);  
11.       ++idx;  
12.     \}  
13.   \}  
14. \}

通过i，j来控制以特征点为中心的6X6的范围，(i\*i + j\*j < 36)则筛选落在以特征点为中心，半径为6s的圆形区域内的点，然后计算HaarX与HarrY并与通过事先计算好的Gauss

滤波。并计算出每个点的角度。

最后将最大值那个扇形的方向作为该特征点的主方向。在OpenSURF寻找扇形中具有最大值得方向代码如下：

\[cpp\]  [view plain][]   [copy][view plain]

1.  **for**(ang1 = 0; ang1 < 2\*pi;  ang1+=0.15f) \{  
2.     ang2 = ( ang1+pi/3.0f > 2\*pi ? ang1-5.0f\*pi/3.0f : ang1+pi/3.0f);  
3.     sumX = sumY = 0.f;   
4.     **for**(unsigned **int** k = 0; k < Ang.size(); ++k)   
5.     \{  
6.       // get angle from the x-axis of the sample point  
7.       **const** **float** & ang = Ang\[k\];  
8.    
9.       // determine whether the point is within the window  
10.      **if** (ang1 < ang2 && ang1 < ang && ang < ang2)   
11.      \{  
12.        sumX+=resX\[k\];    
13.        sumY+=resY\[k\];  
14.      \}   
15.      **else** **if** (ang2 < ang1 &&   
16.        ((ang > 0 && ang < ang2) || (ang > ang1 && ang < 2\*pi) ))   
17.      \{  
18.        sumX+=resX\[k\];    
19.        sumY+=resY\[k\];  
20.      \}  
21.    \}  
22.   
23.    // if the vector produced from this window is longer than all   
24.    // previous vectors then this forms the new dominant direction  
25.    **if** (sumX\*sumX + sumY\*sumY > max)   
26.    \{  
27.      // store largest orientation  
28.      max = sumX\*sumX + sumY\*sumY;  
29.      orientation = getAngle(sumX, sumY);  
30.    \}  
31.  \}

六、构造SURF特征点的特征描述子

在SIFT算法中，是在特征点周围取16\*16的邻域，并且把该邻域化为4\*4的小区域，每个小区域统计8个方向的梯度，最后得到4\*4\*8=128维的向量，并且把该向量作为SIFT特征点的特征描述子。

在SURF算法中，也是在特征点周围取一个正方形框，框的边长为20s(s是所检测到该特征点所在的尺度)。该框带方向，方向当然就是上面一步检测出来的主方向。然后把该框分为16个子区域，每个子区域统计25个像素的水平方向和垂直方向的haar小波特征，这里的水平方向和垂直方向是相对于主方向而言的。该haar小波特征为水平方向值之和、水平方向绝对值之和和垂直方向之和、垂直方向绝对值之和。这样每个小区域就有4个值，所以每个特征点就是16\*4维的向量，相比SIFT而言，少了一半，这在特征匹配的过程中会大大的加速匹配的速度。

以特征点为中心，沿主方向将20SX20S的图像划分为4X4个子块，每个子块用尺寸2S的Harr模板进行响应值计算，并统计每个子块中

![Center 18][]  
这样就有4X4X4=64维的特征数据。如下图所示：  
![Center 19][]

在计算这个矩形区域时并不是先把它旋转到主方向，而是先计算出每一个点的Harr响应值dx、dy并高斯加权处理后，把dx、dy进行旋转变换，计算

公式如下：

![Center 20][]

在OpenSURF的实现源码中采用的是另外一种方式，通过点旋转公式，把点旋转到主方向上并进行最近邻插值的对应点，公式如下：

![Center 21][]

![Center 22][]

七、匹配

为了加速匹配过程，SURF借助Laplacian（在之前计算Hessian是可以顺便得出，不占用太多的时间）的符号使匹配过程索引加快。这样可以将下面的情况区分开，然后在进行描述符匹配：

![Center 23][]

[Center]: /images/20220520/5ffe4262090c4a53b65ad25360e24f7d.png
[Center 1]: /images/20220520/d1e159d1b350488c84782bba0f9410ea.png
[Center 2]: /images/20220520/ec1e1245a16c4ae998cfeb4f5ea942ff.png
[Center 3]: /images/20220520/46310d64fe464b3eaf9210efb73641cf.png
[Center 4]: /images/20220520/fbbf5683914941ea8750d3b3c2ca9467.png
[Center 5]: /images/20220520/03506add58ef4ee1bf74e0d8cb3fb8ec.png
[Center 6]: /images/20220520/5c3f6a620232432fae483f47e898b732.png
[Center 7]: /images/20220520/7613413f5d614fdf9868b98f81e0ecc0.png
[Center 8]: /images/20220520/9bfc3a461f9342fb8cafccaf94fde486.png
[Center 9]: /images/20220520/c66c414d49eb49c6a1fc46e3ef842cf8.png
[Center 10]: /images/20220520/f67448d67dab4329b25f9d80a51881dd.png
[Center 11]: /images/20220520/bfa54ed10df64b49871cfc52d2202410.png
[Center 12]: /images/20220520/173bc267f4f6458088b5af63b3a26a4b.png
[Center 13]: /images/20220520/41f7907e8bf54001bbd35ebe76143282.png
[Center 14]: /images/20220520/d64a5b45573b4cdebd2abd9315778d51.png
[Center 15]: /images/20220520/aaae6a47ee18481fa570903a69e14495.png
[Center 16]: /images/20220520/df5df4157e2f4a17be57f82b4a7c5668.png
[Center 17]: /images/20220520/cc9c819ba7f54ab7aacf66415ea96624.png
[view plain]: http://blog.csdn.net/cxp2205455256/article/details/41311013#
[Center 18]: /images/20220520/9728ac25672641d382f6e520253c0f93.png
[Center 19]: /images/20220520/97b645ee104f4710af478d208f9d69a2.png
[Center 20]: /images/20220520/d125f45fa8e946a8a2fd1acbf6c6b4ef.png
[Center 21]: /images/20220520/dfd5440e20e745dd9043979c234c5210.png
[Center 22]: /images/20220520/dc28171a0f8b43f7a7101dbdd24faa6b.png
[Center 23]: /images/20220520/99780c41077e479da80204d767c214d9.png