OpenCV图像特征提取学习四，SIFT特征检测算法

一、SIFT特征检测概述

SIFT的全称是Scale Invariant Feature Transform，尺度不变特征变换，由加拿大教授David G.Lowe提出的。SIFT特征具有对旋转、尺度缩放、亮度变化等保持不变性，是一种非常稳定的局部特征。

1.1 SIFT算法具的特点

图像的局部特征，对旋转、尺度缩放、亮度变化保持不变，对视角变化、仿射变换、噪声也保持一定程度的稳定性。
独特性好，信息量丰富，适用于海量特征库进行快速、准确的匹配。
多量性，即使是很少几个物体也可以产生大量的SIFT特征
高速性，经优化的SIFT匹配算法甚至可以达到实时性
扩招性，可以很方便的与其他的特征向量进行联合。

1.2 SIFT特征检测的步骤

尺度空间的极值检测搜索所有尺度空间上的图像，通过高斯微分函数来识别潜在的对尺度和选择不变的兴趣点。
特征点定位在每个候选的位置上，通过一个拟合精细模型来确定位置尺度，关键点的选取依据他们的稳定程度。
特征方向赋值基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向，后续的所有操作都是对于关键点的方向、尺度和位置进行变换，从而提供这些特征的不变性。
特征点描述在每个特征点周围的邻域内，在选定的尺度上测量图像的局部梯度，这些梯度被变换成一种表示，这种表示允许比较大的局部形状的变形和光照变换。

================================================================

二、高斯核构建尺度空间

通过图像的模糊程度来模拟人在距离物体由远到近时物体在视网膜上成像过程，距离物体越近其尺寸越大图像也越模糊，这就是高斯尺度空间，使用不同的参数模糊图像（分辨率不变），是尺度空间的另一种表现形式。
我们知道图像和高斯函数进行卷积运算能够对图像进行模糊，使用不同的“高斯核”可得到不同模糊程度的图像。一副图像其高斯尺度空间可由其和不同的高斯卷积得到：

$L\left ( x,y,\sigma \right )=G\left ( x,y,\sigma \right )\cdot I\left ( x,y \right )$

其中， $G\left ( x,y,\sigma \right )$ 是高斯核函数：

$G\left ( x,y,\sigma \right )=\frac{1}{2\pi \sigma ^{2}}e^{\frac{x^{2}+y^{2}}{2\sigma ^{2}}}$

$\sigma$ 称为尺度空间因子，它是高斯正态分布的标准差，反映了图像被模糊的程度，其值越大图像越模糊，对应的尺度也就越大。 $L\left ( x,y,\sigma \right )$ 代表着图像的高斯尺度空间。

构建尺度空间的目的是为了检测出在不同的尺度下都存在的特征点，而检测特征点较好的算子是高斯拉普拉斯算子 $\Delta ^{2}G$ （LoG）

$\Delta ^{2}=\frac{\partial ^{2}}{\partial x^{2}}+\frac{\partial^{2} }{\partial y^{2}}$

设kk为相邻两个高斯尺度空间的比例因子，则DoG的定义:

$G\left ( x,y,\sigma \right )=\left [ G\left ( x,y,k\sigma \right ) -G\left ( x,y,\sigma \right )\right ]\cdot I\left ( x,y \right )\\ =L\left ( x,y,k\sigma \right )-L\left ( x,y,\sigma \right )$

将相邻的两个高斯空间的图像相减就得到了DoG的响应图像。为了得到DoG图像，先要构建高斯尺度空间，而高斯的尺度空间可以在图像金字塔降采样的基础上加上高斯滤波得到，也就是对图像金字塔的每层图像使用不同的参数σσ进行高斯模糊，使每层金字塔有多张高斯模糊过的图像。降采样时，金字塔上边一组图像的第一张是由其下面一组图像倒数第三张降采样得到。

高斯金字塔有多组，每组又有多层。一组中的多个层之间的尺度是不一样的（也就是使用的高斯参数σ是不同的），相邻两层之间的尺度相差一个比例因子k。如果每组有S层，则 $k=2^{\frac{1}{S}}$ 。上一组图像的最底层图像是由下一组中尺度为2σ的图像进行因子为2的降采样得到的（高斯金字塔先从底层建立）。高斯金字塔构建完成后，将相邻的高斯金字塔相减就得到了DoG金字塔。
高斯金字塔的组数一般是

$o=\left [ log_{2} min\left ( m,n \right )\right ]-a$

o表示高斯金字塔的层数，m，n分别是图像的行和列。减去的系数a可以在 $0-log_{2}\left ( m,n \right )$ 之间的任意值，和具体需要的金字塔的顶层图像的大小有关。
高斯模糊参数 $\sigma$ （尺度空间），可由下面关系式得到:

$\sigma \left ( o,s \right )= \sigma _{0}\cdot 2^{\frac{o+s}{S}}$

其中o为所在的组，s为所在的层， $\sigma _{0}$ 为初始的尺度，S为每组的层数。算法实现中 $\sigma _{0}=1.6$ , $o_{min}=-1$ 就是首先将原图像的长和宽各扩展一倍。S=3,从这里可以得知同一组内相邻层的图像尺度关系：

$\sigma _{s+1}=k\cdot\sigma _{s}=2^{\frac{1}{S}}\cdot \sigma _{s}$

相邻组之间的尺度关系:

$\sigma _{0+1}=2\sigma _{0}$

=================================================================

三、DoG空间极值检测

为了寻找尺度空间的极值点，每个像素点要和其图像域（同一尺度空间）和尺度域（相邻的尺度空间）的所有相邻点进行比较，当其大于（或者小于）所有相邻点时，改点就是极值点。如图所示，中间的检测点要和其所在图像的3×3邻域8个像素点，以及其相邻的上下两层的3×3领域18个像素点，共26个像素点进行比较。

从上面的描述中可以知道，每组图像的第一层和最后一层是无法进行比较取得极值的。为了满足尺度变换的连续性，在每一组图像的顶层继续使用高斯模糊生成3幅图像，高斯金字塔每组有S+3层图像，DoG金字塔的每组有S+2组图像。

设S=3，也就是每组有3层，则 $k=2^{\frac{1}{S}}=2^{\frac{1}{3}}$ ，也就是有高斯金字塔每组有(S−1)3层图像，DoG金字塔每组有(S−1)3层图像，DoG金字塔每组有(S-2)2层图像。在DoG金字塔的第一组有两层尺度分别是σ,kσ，第二组有两层的尺度分别是2σ,2kσ，由于只有两项是无法比较取得极值的（只有左右两边都有值才能有极值）。

由于无法比较取得极值，那么我们就需要继续对每组的图像进行高斯模糊，使得尺度形成 $\sigma ,k\sigma ,k^{2}\sigma ,k^{3}\sigma ,k^{4}\sigma$ ，这样就可以选择中间的三项 $\sigma ,k\sigma ,k^{2}\sigma ,k^{3}\sigma$ 。对应的下一组由上一组降采样得到的三项是 $2k\sigma ,2k^{2}\sigma ,2k^{3}\sigma$ ，其首项 $2k\sigma =2^{\frac{4}{3}}\sigma$ ，刚好与上一组的最后一项 $k^{3}\sigma =2^{\frac{2}{2}}\sigma$ 的尺度连续起来。

=====================================================

四、删除不好的极值点（特征点）

通过比较检测得到的DoG的局部极值点实在离散的空间搜索得到的，由于离散空间是对连续空间采样得到的结果，因此在离散空间找到的极值点不一定是真正意义上的极值点，因此要设法将不满足条件的点剔除掉。可以通过尺度空间DoG函数进行曲线拟合寻找极值点，这一步的本质是去掉DoG局部曲率非常不对称的点。

要剔除掉的不符合要求的点主要有两种：

低对比度的特征点
不稳定的边缘响应点

1.剔除低对比度的特征点

候选特征点x，其偏移量定义为 $\widetilde{x}$ ，其对比度为H(x),其绝对值∣H(x)∣，对H(x)应用泰勒展开式：

$H\left ( x \right )=H+\frac{\partial H^{T}}{\partial x}x+\frac{1}{2}x^{T}\frac{\partial ^{2}H}{\partial x^{2}}x$

由于x是H(x)的极值点，所以对上式求导并令其为0，得到

$\widetilde{x}=\frac{\partial ^{2}H^{-1}}{\partial H^{2}}\cdot \frac{\partial H}{\partial x}$

然后再把求得的 $\widetilde{x}$ 代入到H(x)的泰勒展开式中:

$H\left ( \widetilde{x} \right )= H+\frac{1}{2}\frac{\partial H^{T}}{\partial x}\widetilde{x}$

设对比度的阈值为T，若 $H\left ( \widetilde{x} \right )>T$ ，则该特征点保留，否则剔除掉。

2.剔除不稳定的边缘响应点

在边缘梯度的方向上主曲率值比较大，而沿着边缘方向则主曲率值较小。候选特征点的DoG函数H(x)的主曲率与3×3的Hessian矩阵HH的特征值成正比。

$H=\begin{vmatrix} d_{xx},d_{xy}\\ d_{yy},d_{yx}\\ \end{vmatrix}$

其中， $d_{xx},d_{xy},d_{yy}$ 是候选点邻域对应位置的差分求得的。为了避免求具体的值，可以使用H特征值得比例。设 $a=\lambda _{max}$ 为H的最大特征值， $\beta =\lambda _{min}$ 为H的最小特征值。可以得到：

$T_{r}\left ( H \right )=d^{xx}+d_{yy}=\alpha +\beta$

$Det\left ( H \right )=d_{xx}+d_{yy}-d_{xy}^{2}=\alpha \cdot \beta$

其中， $Tr\left ( H \right )$ 为矩阵H的迹， $Det\left ( H \right )$ 为矩阵H的行列式。设 $\gamma =\frac{\alpha }{\beta }$ 表示最大特征值和最小特征值的比值，则：

$\frac{Tr\left ( H \right )^{2}}{Det\left ( H \right )}=\frac{\left ( \alpha +\beta \right )^{2}}{\alpha \cdot \beta }=\frac{ \left ( \gamma \beta +\beta \right )^{2}}{\gamma \beta ^{2}}=\frac{\left ( \gamma +1 \right )^{2}}{\gamma }$

上式的结果与两个特征值的比例有关，和具体的大小无关，当两个特征值想等时其值最小，并且随着γγ的增大而增大。因此为了检测主曲率是否在某个阈值 $T\gamma$ 下，只需检测:

$\frac{Tr\left ( H \right )^{2}}{Det\left ( H \right )}> \frac{\left ( T\gamma +1 \right )^{2}}{T\gamma }$

如果上式成立，则剔除该特征点，否则保留。

=========================================================================

五、求取特征点的主方向

经过上面的步骤已经找到了在不同尺度下都存在的特征点，为了实现图像旋转不变性，需要给特征点的方向进行赋值。利用特征点邻域像素的梯度分布特性来确定其方向参数，再利用图像的梯度直方图求取关键点局部结构的稳定方向。找到了特征点，也就可以得到该特征点的尺度 $\sigma$ ，也就可以得到特征点所在的尺度图像.

计算以特征点为中心、以3×1.5 $\sigma$ 为半径的区域图像的幅角和幅值，为了防御噪声的干扰，计算得到圆形区域内所有点的梯度大小与方向之后，以极值点为中心，以1.5σ为高斯分布的标准差，对圆形区域内所有点的梯度大小进行高斯加权（越靠近中心点权重越大）：

$m^{'}\left ( x,y \right )= m\left ( x,y \right )\cdot W\left ( x,y \right )$

$W\left ( x,y \right )=\frac{w\left ( x,y \right )}{\sum w\left ( x,y \right )}$

$w\left ( x,y \right )=e^{\frac{d^{2}}{2\cdot \left ( 1.5\sigma \right )^{2}}}$

每个点 $L\left ( x,y \right )$ 的梯度的模 $m\left ( x,y \right )$ 以及方向 $\theta \left ( x,y \right )$ 可以通过以下公式求得：

计算得到梯度方向后，就要使用直方图统计特征点邻域内像素对应的梯度方向和幅值。梯度方向的直方图的横轴是梯度方向的角度（梯度方向的范围是0到360度，直方图每36度一个柱共10个柱，或者每45度一个柱共8个柱），纵轴是梯度方向对应梯度赋值的累加。

为了到更精确的方向，通常还可以对离散的梯度直方图进行插值拟合。具体而言，关键点的方向可以由和主峰值最近的三个柱值通过抛物线插值得到。在梯度直方图中，当存在一个相当于主峰值80%能量的柱值时，则可以将这个方向认为是该特征点辅助方向，一个特征值可以检测到多个方向。

得到特征点的主方向后，对于每个特征点可以得到三个信息 $\left ( x,y,\sigma ,\theta \right )$ ，即位置、尺度和方向。由此可以确定一个SIFT特征区域，一个SIFT特征区域由三个值表示，中心表示特征点位置，半径表示关键点的尺度，箭头表示主方向。具有多个方向的关键点可以被复制成多份，然后将方向值分别赋给复制后的特征点，一个特征点就产生了多个坐标、尺度相等，但是方向不同的特征点。

=========================================================================

六、生成特征描述

通过以上的步骤已经找到了SIFT特征点位置、尺度和方向信息，下面就需要使用一组向量来描述关键点也就是生成特征点描述子，这个描述符不只包含特征点，也含有特征点周围对其有贡献的像素点。描述子应具有较高的独立性，以保证匹配率。
特征描述符的生成大致有三个步骤：

校正旋转主方向，确保旋转不变性。
生成描述子，最终形成一个128维的特征向量
归一化处理，将特征向量长度进行归一化处理，进一步去除光照的影响。

为了保证特征矢量的旋转不变性，要以特征点为中心，在附近邻域内将坐标轴旋转θθ（特征点的主方向）角度，即将坐标轴旋转为特征点的主方向。旋转后邻域内像素的新坐标为：

$\begin{vmatrix} x^{'}\\ y^{'}\\ \end{vmatrix}=\begin{bmatrix} cos\theta ,-sin\theta \\ sin\theta ,-cos\theta \\ \end{bmatrix}\begin{vmatrix} x\\ y\\ \end{vmatrix}$

旋转后以主方向为中心取 8×8的窗口。下图所示，左图的中央为当前关键点的位置，每个小格代表为关键点邻域所在尺度空间的一个像素，求取每个像素的梯度幅值与梯度方向，箭头方向代表该像素的梯度方向，长度代表梯度幅值，然后利用高斯窗口对其进行加权运算。最后在每个4×4的小块上绘制8个方向的梯度直方图，计算每个梯度方向的累加值，即可形成一个种子点，如右图所示。每个特征点由4个种子点组成，每个种子点有8个方向的向量信息。这种邻域方向性信息联合增强了算法的抗噪声能力，同时对于含有定位误差的特征匹配也提供了比较理性的容错性。

与求主方向不同，此时每个种子区域的梯度直方图在0-360之间划分为8个方向区间，每个区间为45度，即每个种子点有8个方向的梯度强度信息。对每个关键点使用4×4共16个种子点来描述，这样一个关键点就可以产生128维的SIFT特征向量。

=========================================================================

代码实现：

#include"stdafx.h"
#include <opencv2/opencv.hpp>
#include <iostream>
#include "math.h"#include <opencv2/features2d.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>using namespace cv;
using namespace std;
//using namespace cv::features2d;int main(int argc, char** argv) {Mat src = imread("F:/photo/i.jpg", IMREAD_GRAYSCALE);if (src.empty()) {printf("could not load image...\n");return -1;}namedWindow("input image", WINDOW_AUTOSIZE);imshow("input image", src);int numFeatures = 400;Ptr<SIFT> detector = SIFT::create(numFeatures);vector<KeyPoint> keypoints;detector->detect(src, keypoints, Mat());printf("Total KeyPoints : %d\n", keypoints.size());Mat keypoint_img;drawKeypoints(src, keypoints, keypoint_img, Scalar::all(-1), DrawMatchesFlags::DEFAULT);namedWindow("SIFT KeyPoints", WINDOW_AUTOSIZE);imshow("SIFT KeyPoints", keypoint_img);waitKey(0);return 0;
}

效果实现: