opencv:特征检测与提取——SIFT and SURF

一：前言

特征检测是计算机视觉和图像处理中的一个概念。它指的是使用计算机提取图像信息，决定每个图像的点是否属于一个图像特征。特征检测的结果是把图像上的点分为不同的子集，这些子集往往属于孤立的点、连续的曲线或者连续的区域。在opencv中，我们常用的特征检测算法有SIFT，SURF以及HOG，LBP，Haar特征检测等等，下面我们将分别介绍这几个算法。篇幅有点长，我尽量每个地方都能说到，有错误的地方还请指正！

二：SIFT特征检测

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

我们从SIFT算法的基本步骤一步一步的说起吧。

1. 主要是建立尺度空间，高斯金字塔以及高斯差分金字塔(DOG)。
2. 尺度空间中寻找关键点（keypoints），并对其进行精确定位。
3. 求解关键点的梯度幅值和幅角，对关键点进行方向赋值。
4. 对关键点进行描述，形成128维的空间向量。

第一步：尺度空间的构造

先来了解一下什么是尺度空间：自然界中的物体随着观测尺度不同有不同的表现形态。例如我们形容建筑物用“米”，观测分子、原子等用“纳米”。更形象的例子比如高德地图，滑动鼠标轮可以改变观测地图的尺度，看到的地图绘制也不同；还有电影中的拉伸镜头等等……

尺度空间中各尺度图像的模糊程度逐渐变大，能够模拟人在距离目标由近到远时目标在视网膜上的形成过程。因此描述图像的尺度越大，看到的图像也就越模糊。

我们为什么要讨论尺度空间呢，主要是因为我们在用计算机视觉分析未知场景时，计算机它并不知道图像中物体的尺度。我们需要同时考虑图像在多尺度下的描述，获知感兴趣物体的最佳尺度。另外如果不同的尺度下都有同样的关键点，那么在不同的尺度的输入图像下就都可以检测出来关键点匹配，也就是尺度不变性。

了解了尺度空间的概念后，我们接下来的任务就是构建高斯金字塔以及DOG金子塔。

对高斯金子塔和DOG金子塔概念不是很清晰的朋友可以参考一下这篇博客(点击这里)，在我看来这位大佬已经说的很清楚了。最主要的图就是下面这张图：

建立高斯金子塔主要包括两个步骤，首先通过高斯掩膜对图像卷积，进行高斯模糊，然后对高斯模糊之后的图像进行降采样，也就是对前一组图像的行列像素取一半。

上图中，左边的就是高斯金字塔，示意图中每五层为一组，层之间的顺序按照从下到上进行索引，如上图Gaussian第一组（最下面的那一组）*有5层，标记为0,1,2,3,4层 ，第二组（中间那一组）*有5层，同样标记为0,1,2,3,4  。第一组的第0层是原始图像放大一倍后得到的，第1层是由第0层经过高斯平滑后得到的，... ，第4层是由第3层经过高斯平滑得到的；而第二组的第0层是由第1组的倒数第三张图像降采样得到的，我们后面会介绍为什么是这样获取的。

如何建立DoG金字塔，大家可以从上面的示意图看出，右边的DOG金字塔是通过左边的Gaussian金字塔同一组当中相邻两层相减而得到的，如：DOG金字塔的第1组第0层是通过第Gaussian金字塔的第1组第0层减去第1组的第1层得到的，依次类推。

 第二步：关键点搜索与定位

在进行关键点的搜索与定位之前，我们需要弄明白尺度空间连续这一概念，如下图所示：

我们需要认识到以下几点：(a)尺度空间中同一组相邻2层之间的尺度倍数为k（）,S为尺度空间的组数，默认的S为3）；（b）尺度空间中相邻两组之间的尺度倍数为2（如第1组第0层尺度为sigma，第2组第0层尺度为2*sigma）；（3）Gaussian金字塔的的组数为S+3，DoG金字塔的组数为S+2，而最后寻找极值点的尺度空间为S(上图中红色框标注的区域)。

                                                                                    

表示图像的初始尺度，o是组数的索引值，r是每一组层数的索引值，s即尺度空间的组数。

进入正题，为什么说尺度空间是连续的，看下表

之前我们说过，第2组的第0层是由第1组的倒数第三张图像降采样得到的，这个时候可以告诉大家原因了，Gaussian1组的倒数第三张图像的尺度是，假设Gaussian2组的第一张图像是由此图片降采样得到的，可以得到Gaussian2组的第一张图像的尺度也是；经过这样的采样可以保证Finalscale space是连续的。读者仔细观察表1就可以发现Final scale space 中的尺度是连续的。

                                                                     

尺度空间的极值点代表这样一类点，如：角点、暗区域的亮点以及亮区域的黑点，这些点是图像中十分突出的点，而且这些点相对比较稳定，在对两幅图像当中具有相同的物体进行SIFT特征点检测的时候，可以分别提取到这些稳定点，然后对这些进行特征点进行匹配。极值点的搜索与定位，如下图所示：

寻找DoG极值点时，每一个像素点和它所有的相邻点比较，当其大于（或小于）它的图像域和尺度域的所有相邻点时，即为极值点。如下图所示，比较的范围是个3×3的立方体：中间的检测点和它同尺度的8个相邻点，以及和上下相邻尺度对应的9×2个点——共26个点比较，以确保在尺度空间和二维图像空间都检测到极值点。

但是以上极值点的搜索是在离散空间进行搜索的，由下图可以看到，在离散空间找到的极值点不一定是真正意义上的极值点，因此需要对极值点进行精确定位。可以通过对尺度空间DoG函数进行曲线拟合寻找极值点来减小这种误差。

利用DoG函数在尺度空间的Taylor展开式：

则极值点为：

 第三步：方向赋值

我们已经找到了关键点。为了实现图像旋转不变性，需要根据检测到的关键点局部图像结构为特征点方向赋值。

在上面，精确定位关键点后也找到该特征点的尺度值σ，根据这一尺度值，得到最接近这一尺度值的高斯图像：

使用有限差分，计算以关键点为中心，以3×1.5σ为半径的区域内图像梯度的幅角和幅值，这里的σ是指关键点所在层相对于所在组的基准层的高斯尺度图像的尺度，即,r为所在层的索引值，具体公式如下：

x,y是关键点的位置坐标，L(x,y)为x,y这一点的像素值。

在完成关键点邻域内高斯图像梯度计算后，使用直方图统计邻域内像素对应的梯度方向和幅值。直方图可以看做是离散点的概率表示形式。此处方向直方图的核心是统计以关键点为原点，一定区域内的图像像素点对关键点方向生成所作的贡献。

梯度方向直方图的横轴是梯度方向角，纵轴是剃度方向角对应的梯度幅值累加值。梯度方向直方图将0°~360°的范围分为36个柱，每10°为一个柱。下图是从高斯图像上求取梯度，再由梯度得到梯度方向直方图的例图。

在计算直方图时，每个加入直方图的采样点都使用圆形高斯函数函数进行了加权处理，也就是进行高斯平滑。这主要是因为SIFT算法只考虑了尺度和旋转不变形，没有考虑仿射不变性。通过高斯平滑，可以使关键点附近的梯度幅值有较大权重，从而部分弥补没考虑仿射不变形产生的特征点不稳定。直方图的峰值就代表该关键点处的主方向，而辅方向是代表峰值大于主方向峰值80%以上的方向。在opencv中梯度直方图共有36个bins，上图只是代表示意图。所以一个关键点可能检测得到多个方向，这可以增强匹配的鲁棒性。发明这个算法的论文指出大概有15%关键点具有多方向，但这些点对匹配的稳定性至为关键。

获得图像关键点主方向后，每个关键点有三个信息(x,y,σ,θ)：位置、尺度、方向。由此我们可以确定一个SIFT特征区域。通常使用一个带箭头的圆或直接使用箭头表示SIFT区域的三个值：中心表示特征点位置，半径表示关键点尺度（r=2.5σ）,箭头表示主方向。具有多个方向的关键点可以复制成多份，然后将方向值分别赋给复制后的关键点。如下图：

第四步：关键点描述

之前找到的关键点即SIFT特征点，包含位置、尺度和方向的信息。接下来的步骤是关键点描述，即用用一组向量将这个关键点描述出来，这个描述子不但包括关键点，也包括关键点周围对其有贡献的像素点。用来作为目标匹配的依据（所以描述子应该有较高的独特性，以保证匹配率），也可使关键点具有更多的不变特性，如光照变化、3D视点变化等。

特征描述子与关键点所在尺度有关，因此对梯度的求取应在特征点对应的高斯图像上进行。将关键点附近划分成d×d个子区域，每个子区域尺寸为mσ个像元（d=4，m=3，σ为尺特征点的尺度值）。考虑到实际计算时需要双线性插值，故计算的图像区域为mσ(d+1)，再考虑旋转，则实际计算的图像区域为，如下图所示：

为了保证特征矢量具有旋转不变性，要以特征点为中心，在附近邻域内旋转θ角，即旋转为特征点的方向。

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

将旋转后区域划分为d×d个子区域（每个区域间隔为mσ像元），在子区域内计算8个方向的梯度直方图，绘制每个方向梯度方向的累加值，形成一个种子点。

与求主方向不同的是，此时，每个子区域梯度方向直方图将0°~360°划分为8个方向区间，每个区间为45°。即每个种子点有8个方向区间的梯度强度信息。由于存在d×d，即4×4个子区域，所以最终共有4×4×8=128个数据，形成128维SIFT特征矢量。

对特征矢量需要加权处理，加权采用mσd/2的标准高斯函数。为了除去光照变化影响，还有一步归一化处理。

五：SIFT算法实践

首先让我们先看一下SIFT的构造函数。

SIFT::SIFT(int nfeatures=0, int nOctaveLayers=3, double contrastThreshold=0.04, double edgeThreshold=
10, double sigma=1.6)

代码奉上。

#include<opencv2/xfeatures2d.hpp>
#include<opencv2/opencv.hpp>
#include<iostream>

using namespace std;
using namespace cv;
using namespace cv::xfeatures2d;

int main1(int argc, char** argv) {
	Mat src = imread("D:/opencv/yuner.jpg", IMREAD_GRAYSCALE);
	if (src.empty()) {
		cout << "could not install image" << endl;
		return -1;
	}
	namedWindow("input image", CV_WINDOW_AUTOSIZE);
	Mat dst;
	//resize(src, dst, Size(), 0.5, 0.5);
	imshow("input image", src);

	int numFeatures = 1000;
	//创建能够检测numFeatures个特征的特征检测器指针对象
	//SiftFeatureDetector  siftdtc;
	Ptr<SIFT> detector = SIFT::create(numFeatures);
	vector<KeyPoint> keyPoints;
	//特征点检测
	detector->detect(src, keyPoints, Mat());
	printf("the total keyPoints is:%d", keyPoints.size());

	//绘制特征点
	Mat keyPoints_image;
	drawKeypoints(src, keyPoints, keyPoints_image, Scalar::all(-1), DrawMatchesFlags::DEFAULT);//颜色随机，带有方向
	namedWindow("SIFT Image", CV_WINDOW_AUTOSIZE);
	imshow("SIFT Image", keyPoints_image);

	waitKey(0);
	return 0;
}

程序执行的效果如下所示：

 

三：SURF特征检测算法

SURF 算法，全称是 Speeded-Up Robust Features。该算子在保持 SIFT 算子优良性能特点的基础上，同时解决了 SIFT 计算复杂度高、耗时长的缺点，对兴趣点提取及其特征向量描述方面进行了改进，且计算速度得到提高。

SURF算法之所以比SIFT算法更快更强大，主要是运用了一个黑塞矩阵(Hessian),SIFT采用的是DOG图像，而SURF采用的是Hessian矩阵行列式近似值图像。图像中某个像素点的Hessian矩阵表示如下(即每一个点都可以求一个黑塞矩阵)：

                                                                        

由于我们的特征点需要具备尺度无关性，所以在进行Hessian矩阵构造前，需要对其进行高斯滤波。其中，为高斯滤波后图像g（σ）在各个方向的二阶导数。对图像的二阶导数不熟悉的同学可以参考这篇博客(点击这里)，二阶导数计算公式如下：

                                      

其中L(x)=g(h(x))（假设，I(x)为原始图像的灰度值，L(x)是将I(x)高斯滤波处理后的图像）。 为了找出图像中的特征点，需要对原图进行变换，在SIFT算法中，是在DOG图像中进行，在surf算法中，该变换图就是原图每个像素的Hessian矩阵行列式的近似值构成的。公式为 ：

                                                                       

0.9为经验值，高斯函数的高阶微分与离散的图像函数做卷积运算时相当于使用高斯滤波模板对图像做滤波处理。但在实际运用中，高斯二阶微分进行离散化和裁剪处理得到盒子滤波器近似代替高斯滤波板进行卷积运算：

图中灰色像素代表0；上图1为y方向的先高斯滤波然后二阶求导的处理，近似处理为图3，图2为x和y方向上的先高斯滤波然后二阶混合偏导，近似为图4。 有了这个近似的模板以后，计算高斯滤波和二阶导数两个步骤就可以一个步骤完成，同时，为了提高计算效率，还引入了积分图像的概念，提高了速度。积分图像的概念可以参考这篇文章(点击这里)。

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

SURF算法不会对图片进行下采样，SURF算法会先从9 x 9尺寸的盒子滤波器开始，对盒子滤波器的尺寸进行扩展，9 x 9尺寸的盒子滤波器是为为1.2时的高斯二阶微分函数经过离散和裁剪后的滤波模板。在SURF中，我们保持图像不变，仅仅改变高斯滤波窗口的大小来获得不同尺度的图像，即构成了尺度空间。

每一层对应的与滤波模板尺寸之间的关系式为,建议将尺度空间分为四组，每组中包括四层.

为了保持尺度空间的连续性，SURF算法尺度空间相邻组中有部分层重叠，同时每组中的盒子滤波器的尺寸都是逐渐增大的。

在SURF算法的尺度空间中，每一组中任意一层包括三种盒子滤波器。对一幅输入图像进行滤波后通过Hessian行列式计算公式可以得到对于尺度坐标下的Hessian行列式的值，所有Hessian行列式值构成一幅Hessian行列式图像。

(等同于上述的)一幅灰度图像经过尺度空间中不同尺寸盒子滤波器的滤波处理，可以生成多幅Hessian行列式图像，从而构成了图像金字塔。

在每一组中选取相邻的三层Hessian行列式图像，对于中间层的每一个Hessian行列式值都可以做为待比较的点，在空间中选取该点周围的26个点进行比较大小，若该点大于其他26个点，则该点为特征点。从上诉过程可以知道，当尺度空间每组由四层构成时，非极大值抑制只会在中间两层进行，相邻的组之间不进行比较。

低于Hessian行列式阀值的点不能作为最终的特征点。在实际选择阀值时，根据实际应用中对特征点数量和精确度的要求改变阀值。阀值越大，得到的特征点的鲁棒性越好。在处理场景简单的图像时，其阀值可以适当的调低。在复杂的图像中，图像经旋转或者模糊后特征点变化的数量较大，测试需要适当提高阀值。

为了保证旋转不变性，在SURF中，不统计其梯度直方图，而是统计特征点领域内的Harr小波特征。即以特征点为中心，计算半径为6s(S为特征点所在的尺度值)的邻域内，统计60度扇形内所有点在x(水平)和y(垂直)方向的Haar小波响应总和(Haar小波边长取4s)，并给这些响应值赋高斯权重系数，使得靠近特征点的响应贡献大，而远离特征点的响应贡献小，然后60度范围内的响应相加以形成新的矢量，遍历整个圆形区域，选择最长矢量的方向为该特征点的主方向。这样，通过特征点逐个进行计算，得到每一个特征点的主方向。该过程的示意图如下：

在特征点周围取一个正方形框，框的边长为20s(s是所检测到该特征点所在的尺度)。该框带方向，方向当然就是第4步检测出来的主方向了。然后把该框分为16个子区域，每个子区域统计25个像素的水平方向和垂直方向的haar小波特征，这里的水平和垂直方向都是相对主方向而言的。该haar小波特征为水平方向值之和，水平方向绝对值之和，垂直方向之和，垂直方向绝对值之和。该过程的示意图如下所示：

这样每个小区域就有4个值，所以每个特征点就是16*4=64维的向量，即每个关键点描述是64维，比SIFT描述少了一半，这在特征匹配过程中会大大加快匹配速度。

SURF算法实践

构造函数：

SURF::SURF(double hessianThreshold, int nOctaves=4, int nOctaveLayers=2, bool extended=true,bool upright=false )

1. hessianThreshold 为Hessian 矩阵行列式响应值的阈值
2. nOctaves 为图像堆的组数
3. nOctaveLayers 为图像堆中每组中的中间层数，该值加2 等于每组图像中所包含的层数
4. extended 表示是128 维描述符，还是64 维描述符，为true 时，表示128 维描述
5. upright 表示是否采用U-SURF 算法，为true 时，采用U-SURF 算法

特征检测实践：

#include<opencv2/xfeatures2d.hpp>
#include<opencv2/opencv.hpp>
#include<iostream>

using namespace std;
using namespace cv;
using namespace cv::xfeatures2d;

int main(int argc, char** argv) {
	Mat src = imread("D:/opencv/yuner.jpg", IMREAD_GRAYSCALE);
	if (src.empty()) {
		cout << "could not install image" << endl;
		return -1;
	}
	namedWindow("input image", CV_WINDOW_AUTOSIZE);
	Mat dst;
	//resize(src, dst, Size(), 0.5, 0.5);
	imshow("input image", src);

	int minHessian = 400;
	//创建能够检测numFeatures个特征的特征检测器指针对象
	Ptr<SURF> detector = SURF::create(minHessian);
	vector<KeyPoint> keyPoints;
	//特征点检测
	detector->detect(src, keyPoints, Mat());
	printf("the total keyPoints is:%d", keyPoints.size());

	//绘制特征点
	Mat keyPoints_image;
	drawKeypoints(src, keyPoints, keyPoints_image, Scalar::all(-1), DrawMatchesFlags::DEFAULT);//颜色随机，带有方向
	namedWindow("SIFT Image", CV_WINDOW_AUTOSIZE);
	imshow("SIFT Image", keyPoints_image);

	waitKey(0);
	return 0;
}

运行效果如下：
 

 

SURF算法与SIFT算法总结对比

1. 在生成尺度空间方面，SIFT算法利用的是差分高斯金字塔与不同层级的空间图像相互卷积生成。SURF算法采用的是不同尺度的box filters与原图像卷积
2. 在特征点检验时，SIFT算子是先对图像进行非极大值抑制，再去除对比度较低的点。然后通过Hessian矩阵去除边缘的点。而SURF算法是先通过Hessian矩阵来检测候选特征点，然后再对非极大值的点进行抑制
3. 在特征向量的方向确定上，SIFT算法是在正方形区域内统计梯度的幅值的直方图，找到最大梯度幅值所对应的方向。SIFT算子确定的特征点可以有一个或一个以上方向，其中包括一个主方向与多个辅方向。SURF算法则是在圆形邻域内，检测各个扇形范围内水平、垂直方向上的Haar小波响应，找到模值最大的扇形指向，且该算法的方向只有一个。
4. SIFT算法生成描述子时,是将16 x 16 的采样点划分为 4 x 4 的区域，从而计算每个分区种子点的幅值并确定其方向，共计128维度。SURF是将20 x 20s 的正方形区域划分为4 x 4的小方格，每个子区域采样25个点，计算小波相应，一共64维

综上，SURF算法在各个步骤上都简化了一些繁琐的工作，仅仅计算了特征点的一个主方向，生成的特征描述子也与前者相比降低了维数。