- 算法介绍
- SIFT算法步骤
- 图像金字塔
- 空间极值点检测
- 关键点方向分配
- 特征点描述符
- API即openCV应用
1.算法介绍
SIFT--尺度不变特征转换算法(Scale-invariant feature transform)，通常用来侦测与描述影像中的局部性特征，在空间尺度中寻找极值点，并提取出其位置，尺度，旋转不变量。
SIFT的应用范围包括物体识别，机器人感知与导航，3D模型建立，影像追踪等。
局部影像特征的描述与侦测可以帮助识别物体，**SIFT特征是基于物体上的一些局部外观的兴趣点而与影像的大小和旋转无关**
**SIFT算法的实质是在不同的尺度空间上查找关键点**，并计算出关键点的方向。SIFT所查找到的关键点十分突出，不会因光照，噪音而变化2.SIFT算法步骤
3.图像金字塔
3.1、高斯金字塔
高斯金字塔是采用高斯函数对图像进行模糊以及降采样处理得到的。
高斯模糊系数计算公式：
高斯函数与图像卷积
分离高斯卷积3.2的图像卷积方法，速度比较慢，同时图像边缘信息也会损失严重。可以使用分离的高斯卷积（即先用1XN的模板沿着X方向对图像卷积一次，然后用NX1的模板沿着Y方向对图像再卷积一次，其中N=[(6σ+1)]且向上取最邻近奇数），这样即省时，也见笑了直接卷积对图像边缘信息的算是
高斯金字塔源码
for (o = 0; o < octvs; o++)//金字塔组数为octvs，
for (i = 0; i < intvls + 3; i++)//每一组有intvls + 3 层，intvls一般为3
{
if (o == 0 && i == 0)//如果是第一组第1层
gauss_pyr[o][i] = cvCloneImage(base);//base 为原始灰度图像经过升采样或降采样得到的图像
/* base of new octvave is halved image from end of previous octave */
else if (i == 0)//建立非第一组的第1层
gauss_pyr[o][i] = downsample(gauss_pyr[o - 1][intvls]);//降采样图像
/* blur the current octave's last image to create the next one */
else//建立非第一组的非第1层
{
gauss_pyr[o][i] = cvCreateImage(cvGetSize(gauss_pyr[o][i - 1]),IPL_DEPTH_32F, 1);
cvSmooth(gauss_pyr[o][i - 1], gauss_pyr[o][i],CV_GAUSSIAN, 0, 0, sig[i], sig[i]);// sig[i]为模糊系数
}//cvSmooth 为平滑处理函数，也即模糊处理。CV_GAUSSIAN 为选用高斯函数对图像模糊
return gauss_pyr;//返回建好的金字塔3.2 高斯差分金字塔
高斯差分函数（DOG算子）的极大值和极小值能够产生比较稳定的图像特征，与尺度归一化的高斯拉普拉斯函数比较近似。
其中k-1是个常熟，不影响极值点位置的求取差分金字塔的建立
差分金字塔实在高斯金字塔的基础上操作的，其建立过程就是：在高斯金字塔中的每组中相邻两层相减（下一层减上一层，就生成了高斯差分金字塔）
for (o = 0; o < octvs; o++)//octvs为高斯金字塔组数
for (i = 0; i < intvls + 2; i++)//因为相减，故高斯金字塔中每组有（intvls + 2）层图像
{
dog_pyr[o][i] = cvCreateImage(cvGetSize(gauss_pyr[o][i]),IPL_DEPTH_32F, 1);
cvSub(gauss_pyr[o][i + 1], gauss_pyr[o][i], dog_pyr[o][i], NULL);//cvSub为opencv内置相减函数
}
return dog_pyr;//返回高斯差分金字塔4.空间极值点（关键点）检测
关键点是由DOG空间的局部极值点组成，关键点的初步探索是通过同一组内各DOG相邻两层图像之间比较完成的。为了寻找DOG函数的极值点，每一个像素点要和它相邻点比较，看是否比它的图像域和尺度域的相邻点大或者小。如下图所示，中间检测点要与其8个相邻点和上下相邻尺度的18个点共26个点比较。
源码：
if (val > 0)//极大值检测{
for (i = -1; i <= 1; i++)
for (j = -1; j <= 1; j++)
for (k = -1; k <= 1; k++)
if (val < pixval32f(dog_pyr[octv][intvl + i], r + j, c + k))//pixval32f为提取图像像素位置上的灰度值
return 0;}
else /* check for minimum */
{
for (i = -1; i <= 1; i++)
for (j = -1; j <= 1; j++)
for (k = -1; k <= 1; k++)
if (val > pixval32f(dog_pyr[octv][intvl + i], r + j, c + k))//r c为图像的行数和列数，dog_pyr为高斯差分图
return 0;4.2关键点定位
以上方法检测到的极值点是离散空间的极值点，接下来通过拟合三维二次函数来精确确定关键点的位置和尺度，同时去除低对比度和不稳定的边缘相应点（DOG算子会产生较强的边缘响应），以增强匹配的稳定性。
关键点精确定位
从上图可以看到，离散空间的极值点并不是真正的极值点。利用已知的离散空间点插值得到的连续空间极值点的方法叫做子像素插值为了提高关键点的稳定性，需要对尺度空间DOG函数进行 曲线插值。利用DOG函数在尺度空间的Taylor展开式如下：
上式算式的矩阵表示如下：
其中，x求导并让方程等于零，可以得到极值点的偏移量为：
对应极值点，方程的值为：
其中x^代表对应插值中心的偏移量，当它在任一维度上的偏移量大于0.5时，意味着插值中心已经偏移到它的临近点上了，所以必须改变当前关键点的位置。同时在新的位置上反复插值直到收敛；也有可能是迭代次数已超标，或超出图像边界的范围，这样的点则删除。
过小的点容易受到噪声的干扰而不稳定，所以将小于某个经验值的极值点删除。
同时，在此过程中获取特征点的精确位置（原位置加上拟合的偏移量）以及尺度(σ)
精确定位中的泰勒插值源码
while (i < SIFT_MAX_INTERP_STEPS)//SIFT_MAX_INTERP_STEPS=5为最大迭代次数，避免长时迭代
{
interp_step(dog_pyr, octv, intvl, r, c, &xi, &xr, &xc);// 泰勒展开拟合，xi,xr,xc依次为x、y、σ方向偏移量,
if (ABS(xi) < 0.5 && ABS(xr) < 0.5 && ABS(xc) < 0.5)//如果当前偏移量绝对值中的每个值均小于0.5，退出迭代
break;
c += cvRound(xc);//计算行坐标，cvRound 为四舍五入。
r += cvRound(xr);
intvl += cvRound(xi);
if (intvl < 1 ||//不在计算的图像层中
intvl > intvls ||//高斯差分每组的层数为intvls
c < SIFT_IMG_BORDER ||//靠近图像边缘5个像素的区域不做检测，SIFT_IMG_BORDER=5,
r < SIFT_IMG_BORDER ||
c >= dog_pyr[octv][0]->width - SIFT_IMG_BORDER ||//靠近图像边缘5个像素的区域不做检测
r >= dog_pyr[octv][0]->height - SIFT_IMG_BORDER)
{
return NULL;
}
i++;//迭代计数
}static void interp_step(IplImage*** dog_pyr, int octv, int intvl, int r, int c,double* xi, double* xr, double* xc)
{
CvMat* dD, *H, *H_inv, X;
double x[3] = { 0 };
dD = deriv_3D(dog_pyr, octv, intvl, r, c);//一阶偏导数
H = hessian_3D(dog_pyr, octv, intvl, r, c);//Hessian 矩阵即二阶导数组成的矩阵
H_inv = cvCreateMat(3, 3, CV_64FC1);
cvInvert(H, H_inv, CV_SVD);//求Hessian矩阵的逆矩阵
cvInitMatHeader(&X, 3, 1, CV_64FC1, x, CV_AUTOSTEP);
cvGEMM(H_inv, dD, -1, NULL, 0, &X, 0); //cvGEMM为矩阵乘法，//第一个矩阵的系数；//H_inv、dD第一二个矩阵//-1矩阵前的常数//X结果矩阵
cvReleaseMat(&dD);
cvReleaseMat(&H);
cvReleaseMat(&H_inv);
*xi = x[2];
*xr = x[1];
*xc = x[0];
}5.关键点方向分配
为了让描述符具有旋转不变性，需要利用图像的局部特征为每一个关键点分配一个基准方向。使用图像梯度的方法求取局部结构的稳定方向
5.1、特征点梯度
对于在DOG金字塔中检测出的关键点，采集器所在高斯金字塔图像3σ领域窗口内像素的梯度和方向分布特征
L为关键点所在的尺度空间值，按Lowe的建议，梯度的模值m(x,y)按 σ=1.5σ_oct 的高斯分布加成，按尺度采样的3σ原则，领域窗口半径为 3x1.5σ_oct。5.2 梯度直方图
在完成关键点的梯度计算后，使用直方图统计领域内像素的梯度和方向。梯度直方图将0-360°的方向范围分为36个柱（bins），每个柱10°。其中，直方图的峰值代表了关键点的主方向
5.3 特征点主方向的确定
方向直方图的峰值则代表了该特征点处邻域梯度的方向，以直方图中最大值作为该关键点的主方向。为了增强匹配的鲁棒性，只保留峰值大于主方向峰值80％的方向作为该关键点的辅方向。因此，对于同一梯度值的多个峰值的关键点位置，在相同位置和尺度将会有多个关键点被创建但方向不同。仅有15％的关键点被赋予多个方向，但可以明显的提高关键点匹配的稳定性。实际编程实现中，就是把该关键点复制成多份关键点，并将方向值分别赋给这些复制后的关键点，并且，离散的梯度方向直方图要进行插值拟合处理，来求得更精确的方向角度值。
5.4 梯度图像的平滑处理
为了防止某个梯度方向角度因受到噪声干扰而突变，可以对梯度方向直方图进行平滑处理
梯度直方图抛物线插值
#define interp_hist_peak( l, c, r ) ( 0.5 * ((l)-(r)) / ((l) - 2.0*(c) + (r)) )//插值计算式，l为左侧柱子值，r为左侧柱子值
static void add_good_ori_features(CvSeq* features, double* hist, int n,
double mag_thr, struct feature* feat)//精确主方向及辅方向
{
struct feature* new_feat;
double bin, PI2 = CV_PI * 2.0;//CV_PI=pi
int l, r, i;
for (i = 0; i < n; i++)// 直方图有n=36个小柱子
{
l = (i == 0) ? n - 1 : i - 1;//把小柱子看成是循环的，角度的取值为0-360即一个圆周
r = (i + 1) % n;
//只对小柱子的值大于等于主峰80%且此小柱子比左右两边小柱子都高的柱子进行抛物线插值
if (hist[i] > hist[l] && hist[i] > hist[r] && hist[i] >= mag_thr)// mag_thr为>=80%的最高峰值
{
bin = i + interp_hist_peak(hist[l], hist[i], hist[r]);//interp_hist_peak 插值函数
bin = (bin < 0) ? n + bin : (bin >= n) ? bin - n : bin;//角度取值约束在0-360之间，且是连续循环的
new_feat = clone_feature(feat);//幅值特征点
new_feat->ori = ((PI2 * bin) / n) - CV_PI;//？
cvSeqPush(features, new_feat);
free(new_feat);
}6.特征点描述符
通过以上步骤，对于每一个关键点，拥有三个信息：位置，尺度以及方向。接下来为每个关键点建立一个描述符，使其不随各种变化而变化，提高特征点正确匹配的概率。
将关键点附近的区域划分为d*d(Lowe建议d=4)个子区域，每个子区域作为一个种子点，每个种子点有8个方向。考虑到实际计算时，需要采用三线性插值，所需图像窗口边长为3x3xσ_oct x(d+1) 。在考虑到旋转因素(方便下一步将坐标轴旋转到关键点的方向)，如下图6.1所示，实际计算所需的图像区域半径为：
坐标轴旋转至主方向将坐标轴旋转为关键点的方向，以确保旋转不变性
生成梯度直方图将邻域内的采样点分配到对应的子区域内，将子区域的梯度值分配到8个方向上，计算其权值
6.2、三线性插值
插值计算每个种子点八个方向的梯度
采样点在子区域中的下标(x’’,y’’) (图中蓝色窗口内红色点)线性插值，计算其对每个种子点的贡献。如图中的红色点，落在第0行和第1行之间，对这两行都有贡献。对第0行第3列种子点的贡献因子为dr，对第1行第3列的贡献因子为1-dr，同理，对邻近两列的贡献因子为dc和1-dc，对邻近两个方向的贡献因子为do和1-do。则最终累加在每个方向上的梯度大小为：
其中k，m，n为0（像素点超出了对要插值区间的四个邻近子区间所在范围）或为1（像素点处在对要插值区间的四个邻近子区间之一所在范围）。
6.3、特征描述子
如上统计的4X4X8=128个梯度信息即为关键点的特征向量
特征向量形成后，为了去除光照变化的影响，需要对它们进行归一化处理，对于图像灰度值整体漂移，图像各点的梯度是邻域像素相减得到，所以也能去除。得到的描述子向量为H=(h1,h2,…,h128)，归一化后的特征向量为L=(L1,L2,…,L128)，则
描述子生成总括
三线性插值源码
static void interp_hist_entry(double*** hist, double rbin, double cbin,
double obin, double mag, int d, int n)
{
double d_r, d_c, d_o, v_r, v_c, v_o;
double** row, *h;
int r0, c0, o0, rb, cb, ob, r, c, o;
r0 = cvFloor(rbin);//向下取整
c0 = cvFloor(cbin);
o0 = cvFloor(obin);
d_r = rbin - r0;//小数余项
d_c = cbin - c0;
d_o = obin - o0;
for (r = 0; r <= 1; r++)//沿着行方向不超过1个单位长度
{
rb = r0 + r;
if (rb >= 0 && rb < d)//如果此刻还在真正的描述子区间内
{
v_r = mag * ((r == 0) ? 1.0 - d_r : d_r);//d_r = rbin - r0;
row = hist[rb];
for (c = 0; c <= 1; c++)//沿着行方向不超过1个单位长度
{
cb = c0 + c;
if (cb >= 0 && cb < d)
{
v_c = v_r * ((c == 0) ? 1.0 - d_c : d_c);
h = row[cb];
for (o = 0; o <= 1; o++)//沿着直方图方向不超过1个单位角度
{
ob = (o0 + o) % n；//n=8，8个小柱子
v_o = v_c * ((o == 0) ? 1.0 - d_o : d_o);
h[ob] += v_o;
}
}
}
}7.API及open CV应用
参数说明：
1.nfeatures：保留的最佳特征数量，特征按 其得分排序
2.高斯金字塔最小层数，由图像自动计算出
3.constrastThrehold：对比度阈值，用于过滤区域中的弱特征，域值越大，检测器产生的特征越少
4.edgeThreshold:用于过滤掉类似边缘特征的阈值。域值越大，产生的特征越少
5.sigma:高斯输入层级，如果图像分辨率越低，数值越小
#include<opencv2/opencv.hpp>
#include<opencv2/xfeatures2d>
#include<math.h>
#include<iostream>
using namespace cv;
using namespace std;
using namespace cv::xfeatures2d;
Mat src, gray_src;
const char* output_title = "output_win";
int corners = 20;
int max_corners = 50;
void Tomasi_Demo(int, void*);
void SubPixel_Demo(int, void*);
int main(int argc, char** argv) {
src = imread("C:/Users/18929/Desktop/博客项目/项目图片/18.jpg");
if (src.empty()) {
printf("could not load image");
return -1;
}
namedWindow("input image", CV_WINDOW_AUTOSIZE);
imshow("input_img", 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", CV_WINDOW_AUTOSIZE);
imshow("SIFT KeyPoints", keypoint_img);
waitKey(0);
return 0;
}
