SIFT算法详解
Posted 夜半罟霖
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了SIFT算法详解相关的知识，希望对你有一定的参考价值。
大纲
引言
SIFT算法是为了解决图片的匹配问题，想要从图像中提取一种对图像的大小和旋转变化保持鲁棒的特征，从而实现匹配。这一算法的灵感也十分的直观：人眼观测两张图片是否匹配时会注意到其中的典型区域（特征点部分），如果我们能够实现这一特征点区域提取过程，再对所提取到的区域进行描述就可以实现特征匹配了。于是问题就演变成了以下几个子问题：
- 应该选取什么样的点作为特征点呢？：人眼对图像中的高频区域更加的敏感，由此我们应该选择变化剧烈的边缘或者角点进行检测，这里我们选择检测角点。（这里的intuition不是很明白，但感觉来说可能是我们提取的是特征点进行匹配，而边缘往往是由多个点组成。）
- 怎样使得选取的特征点有尺度不变性呢？：使用高斯金字塔获取不同尺寸下的图像变体，由这些变体获得尺度不变特征。
- 怎样使得选取的特征点有缩放不变性呢？：使用高斯金字塔获取不同尺寸下的图像变体，每个变体里都提取特征点再缩放回原大小以获得图像不同尺寸下的特征点。
- 怎样使得选取的特征点有旋转不变性呢？：后续采取将特征点区域旋转到主方向的设定以获得旋转不变性。
- 如何描述特征点区域呢？： 使用特征点区域内每个方向的梯度赋值，类似HOG算子。
以上几个问题在SIFT算法里都用了很有意思的trick，后续会一一介绍。
一、高斯金字塔
引入高斯金字塔的目的在引言中已经介绍过了–解决图片缩放及尺度变化下特征提取的问题，高斯金字塔的Intuition有两个：1. 人看物体时近大远小，可以对图片下采样实现（金字塔->组）；2. 人看物体时近处清晰，远处模糊，可以对图像高斯平滑实现（高斯->层）；具体的推导可以参见我的另一篇博客：
Opencv学习笔记（六）图像金字塔
在SIFT里，高斯金字塔的层数和组数有着如下设定：
组数：
O
=
[
l
o
g
2
m
i
n
(
M
,
N
)
]
−
3
O=[log_2min(M,N)]-3
O=[log2min(M,N)]−3
层数：
S
=
n
+
3
S=n+3
S=n+3
组数的设定是来自于提出SIFT算法的原始论文给出的经验值，理论上来说知要
O
≤
[
l
o
g
2
m
i
n
(
M
,
N
)
]
O\\leq[log_2min(M,N)]
O≤[log2min(M,N)]即可，层数的设定则是有着理论依据的，这里的
n
n
n是我们想要提取特征点的图片层数，由于提取出高斯金字塔后需要计算层间差分以获得高斯差分金字塔(DOG, Difference of Gausssian)，高斯金字塔层数需要比DOG层数多1，而计算特征值时要求在尺度层面，即上下相邻层间计算，则DOG层数要比特征层数多2，则要求
S
=
n
+
3
S=n+3
S=n+3。
附：在SIFT中对于高斯金字塔有几个需要补充的知识点。
- SIFT算法中的S指的是"scale"即尺度的概念，这一概念有别于传统的图像尺寸，而是一种连续变化的参数，在高斯金字塔里的 σ \\sigma σ就是这样一个尺度空间的参数，更特殊的，高斯核是唯一可以产生连续多尺度变化的线性核，这就是为什么我们用高斯金字塔，这一部分内容涉及到尺度空间的概念，可以自行了解。
- 图像金子塔一文中提到了两个基础参数，层间尺度变化系数 k k k和初始尺度 σ 0 \\sigma_0 σ0，在原始算法中规定 k = 2 o + r n ， r = 0 , 1... , n + 2 k=2^o+\\fracrn，r=0,1...,n+2 k=2o+nr，r=0,1...,n+2，由于每组的起始层是由前一组倒数第三层降采样得到，实际上保证了每组开始的尺度为 σ 0 , 2 σ 0 , 3 σ 0 . . . \\sigma_0,2\\sigma_0,3\\sigma_0... σ0,2σ0,3σ0...。此外我们提到 σ 0 = 1.6 \\sigma_0=1.6 σ0=1.6，而实际拍摄图片时相机已经原始景物进行了一次尺度变换，此变换的系数原文设定经验值为 σ ′ = 0.5 \\sigma'=0.5 σ′=0.5，由此我们得到 σ 0 = 1. 6 2 − σ ′ 2 = 1.52 \\sigma_0=\\sqrt1.6^2-\\sigma'^2=1.52 σ0=1.62−σ′2=1.52。这一式子的来历为：高斯滤波体现为高斯核与原图进行卷积 f ( x ) ⨂ G ′ ( x ) f(x)\\bigotimes G'(x) f(x)⨂G′(x)，两次高斯滤波的级联相当于 f ( x ) ⨂ ( G ′ ( x ) ⨂ G ( x ) ) f(x)\\bigotimes (G'(x) \\bigotimes G(x)) f(x)⨂(G′(x)⨂G(x))，而两个高斯卷积的结果仍为高斯卷积，且方差满足平方和关系，可用时域卷积或者傅里叶变换证明，详见参考文献。
二、高斯差分金字塔
通过高斯金字塔，我们获取了不同尺度的图片，接下来的问题是如何获取高频区域呢，一个很简单的思路就是按照边缘检测的算法使用差分滤波器如拉普拉斯滤波器、sobel滤波器在图片上滑动找到灰度值变化剧烈的区域。而经前人研究，归一化的高斯拉普拉斯算子的极大值极小值相较于其他特征提取函数可以获得最稳定的图像特征，因此我们打算使用归一化的高斯拉普拉斯算子在多尺度图片上提取特征，但是使用这种方式提取复杂度会很高，又尺度归一化高斯拉普拉斯算子和DOG函数有着如下关系: 说明：本文旨在给出 SIFT 算法的具体实现，而在 SIFT 详解上只是做出简单介绍，在这里可以给大家推荐一篇好文：https://blog.csdn.net/zddblog/article/details/7521424；结合这篇文章和下文的具体代码实现，我相信你能很快掌握并运用 SIFT 算法，加油哦！！！如有疑问大家可以一起讨论学习！！！ SIFT，即尺度不变特征变换（Scale-invariant feature transform，SIFT），是用于图像处理领域的一种描述。这种描述具有尺度不变性，可在图像中检测出关键点，是一种局部特征描述子。该方法于1999年由David Lowe 首先发表于计算机视觉国际会议（International Conference on Computer Vision，ICCV），2004年再次经David Lowe整理完善后发表于International journal of computer vision（IJCV）。截止2014年8月，该论文单篇被引次数达25000余次。 （1） SIFT特征是图像的局部特征，其对旋转、尺度缩放、亮度变化保持不变性，对视角变化、仿射变换、噪声也保持一定程度的稳定性； （2）独特性（Distinctiveness）好，信息量丰富，适用于在海量特征数据库中进行快速、准确的匹配； （3）多量性，即使少数的几个物体也可以产生大量的SIFT特征向量； （4）高速性，经优化的SIFT匹配算法甚至可以达到实时的要求； （5）可扩展性，可以很方便的与其他形式的特征向量进行联合。 目标的自身状态、场景所处的环境和成像器材的成像特性等因素影响图像配准/目标识别跟踪的性能。而SIFT算法在一定程度上可解决： （1）目标的旋转、缩放、平移（RST） （2）图像仿射/投影变换（视点viewpoint） （3）光照影响（illumination） （4）目标遮挡（occlusion） （5）杂物场景（clutter） （6）噪声 SIFT算法的实质是在不同的尺度空间上查找关键点(特征点)，并计算出关键点的方向。SIFT所查找到的关键点是一些十分突出，不会因光照，仿射变换和噪音等因素而变化的点，如角点、边缘点、暗区的亮点及亮区的暗点等。 搜索所有尺度上的图像位置。通过高斯微分函数来识别潜在的对于尺度和旋转不变的兴趣点。 在每个候选的位置上，通过一个拟合精细的模型来确定位置和尺度。关键点的选择依据于它们的稳定程度。 基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向。所有后面的对图像数据的操作都相对于关键点的方向、尺度和位置进行变换，从而提供对于这些变换的不变性。 在每个关键点周围的邻域内，在选定的尺度上测量图像局部的梯度。这些梯度被变换成一种表示，这种表示允许比较大的局部形状的变形和光照变化。 说明：代码实现主要有七个源文件，我和大家一样有时候下载了好几个资源然而没有一个能用的，导致最后都不知道要下载哪一个、是否要下载；所以这里直接把代码贴出来了！ 以上是关于SIFT算法详解的主要内容，如果未能解决你的问题，请参考以下文章 OpenCV+Python特征提取算法与图像描述符之SIFT / SURF / ORB SIFT图像特征提取 python3.6 + opencv3.3代码 不能在 Python OpenCV v4.20 中使用 SIFT
G
(
x
,
y
,
k
σ
)
−
G
(
x
,
y
,
σ
)
≈
(
k
−
1
)
σ
2
∇
2
G
G(x,y,k\\sigma)-G(x,y,\\sigma)\\approx(k-1)\\sigma^2 \\nabla^2 G
G(x,y,kσ)−G(x,y,σ)≈(k−1)σ2∇2G
证明如下：
忽略高斯函数系数:
G
(
x
,
y
,
σ
)
=
1
σ
2
e
x
p
(
−
x
2
+
y
2
2
σ
2
)
∂
G
∂
x
=
−
x
σ
4
e
x
p
(
−
x
2
+
y
2
2
σ
2
)
∂
2
G
∂
x
2
=
−
σ
2
+
x
2
σ
6
e
x
p
(
−
x
2
+
y
2
2
σ
2
)
∇
2
G
(
x
,
y
)
=
∂
2
G
∂
x
2
+
∂
2
G
∂
y
2
=
−
2
σ
2
+
x
2
+
y
2
σ
6
e
x
p
(
−
x
2
+
y
2
2
σ
2
)
∂
G
∂
σ
=
−
2
σ
2
+
x
2
+
y
2
σ
5
e
x
p
(
−
x
2
+
y
2
2
σ
2
)
⇒
σ
∇
2
G
=
∂
G
∂
σ
\\beginaligned &G(x,y,\\sigma)=\\frac1\\sigma^2exp(-\\fracx^2+y^22\\sigma^2)\\\\ &\\frac\\partialG\\partial x=-\\fracx\\sigma^4exp(-\\fracx^2+y^22\\sigma^2)\\\\ &\\frac\\partial^2G\\partial x^2=\\frac-\\sigma^2+x^2\\sigma^6exp(-\\fracx^2+y^22\\sigma^2)\\\\ &\\nabla^2 G(x,y)=\\frac\\partial^2G\\partial x^2+\\frac\\partial^2G\\partial y^2\\\\ &\\qquad \\qquad \\ =\\frac-2\\sigma^2+x^2+y^2\\sigma^6exp(-\\fracx^2+y^22\\sigma^2)\\\\ &\\frac\\partialG\\partial \\sigma=\\frac-2\\sigma^2+x^2+y^2\\sigma^5exp(-\\fracx^2+y^22\\sigma^2)\\\\ &\\Rightarrow\\\\ &\\sigma\\nabla^2G=\\frac\\partial G\\partial \\sigma \\endaligned
G(x,y,σ)=σ21exp(−2σ2x2+y2)∂x∂G=−σ4xexp(−2σ2x2+y2)∂x2∂2G=σ6−σ2+x2exp(−2
一、SIFT算法简介
1、SIFT算法的特点
2、SIFT算法可以解决的问题
二、SIFT算法分为如下四步
1. 尺度空间极值检测
2. 关键点定位
3. 方向确定
4. 关键点描述
三、SIFT具体实现
1、mySift.h
#pragma once //防止头文件重复包含和下面作用一样
#include<iostream>
#include<vector>
#include<opencv2\\core\\core.hpp>
#include<opencv2\\features2d\\features2d.hpp>
using namespace std;
using namespace cv;
/*************************定义常量*****************************/
//高斯核大小和标准差关系，size=2*(GAUSS_KERNEL_RATIO*sigma)+1,经常设置GAUSS_KERNEL_RATIO=2-3之间
const double GAUSS_KERNEL_RATIO = 3;
const int MAX_OCTAVES = 8; //金字塔最大组数
const float CONTR_THR = 0.04f; //默认是的对比度阈值(D(x))
const float CURV_THR = 10.0f; //关键点主曲率阈值
const float INIT_SIGMA = 0.5f; //输入图像的初始尺度
const int IMG_BORDER = 2; //图像边界忽略的宽度，也可以改为 1
const int MAX_INTERP_STEPS = 5; //关键点精确插值次数
const int ORI_HIST_BINS = 36; //计算特征点方向直方图的BINS个数
const float ORI_SIG_FCTR = 1.5f; //计算特征点主方向时候，高斯窗口的标准差因子
const float ORI_RADIUS = 3 * ORI_SIG_FCTR; //计算特征点主方向时，窗口半径因子
const float ORI_PEAK_RATIO = 0.8f; //计算特征点主方向时，直方图的峰值比
const int DESCR_WIDTH = 4; //描述子直方图的网格大小(4x4)
const int DESCR_HIST_BINS = 8; //每个网格中直方图角度方向的维度
const float DESCR_MAG_THR = 0.2f; //描述子幅度阈值
const float DESCR_SCL_FCTR = 3.0f; //计算描述子时，每个网格的大小因子
const int SAR_SIFT_GLOH_ANG_GRID = 8; //GLOH网格沿角度方向等分区间个数
const int SAR_SIFT_DES_ANG_BINS = 8; //像素梯度方向在0-360度内等分区间个数
const float SAR_SIFT_RADIUS_DES = 12.0f; //描述子邻域半径
const int Mmax = 8; //像素梯度方向在0-360度内等分区间个数
const double T = 100.0; //sobel算子去除冗余特征点的阈值
const float SAR_SIFT_GLOH_RATIO_R1_R2 = 0.73f;//GLOH网格中间圆半径和外圆半径之比
const float SAR_SIFT_GLOH_RATIO_R1_R3 = 0.25f;//GLOH网格最内层圆半径和外圆半径之比
#define Feature_Point_Minimum 1500 //输入图像特征点最小个数
#define We 0.2
#define Wn 0.5
#define Row_num 3
#define Col_num 3
#define SIFT_FIXPT_SCALE 48 //不理解，后面可查原论文
/************************sift类*******************************/
class mySift
public:
//默认构造函数
mySift(int nfeatures = 0, int nOctaveLayers = 3, double contrastThreshold = 0.03,
double edgeThreshold = 10, double sigma = 1.6, bool double_size = true) :nfeatures(nfeatures),
nOctaveLayers(nOctaveLayers), contrastThreshold(contrastThreshold),
edgeThreshold(edgeThreshold), sigma(sigma), double_size(double_size)
//获得尺度空间每组中间层数
int get_nOctave_layers() return nOctaveLayers;
//获得图像尺度是否扩大一倍
bool get_double_size() return double_size;
//计算金字塔组数
int num_octaves(const Mat& image);
//生成高斯金字塔第一组，第一层图像
void create_initial_image(const Mat& image, Mat& init_image);
//使用 sobel 算子创建高斯金字塔第一层图像，以减少冗余特征点
void sobel_create_initial_image(const Mat& image, Mat& init_image);
//创建高斯金字塔
void build_gaussian_pyramid(const Mat& init_image, vector<vector<Mat>>& gauss_pyramid, int nOctaves);
//创建高斯差分金字塔
void build_dog_pyramid(vector<vector<Mat>>& dog_pyramid, const vector<vector<Mat>>& gauss_pyramid);
//该函数生成高斯差分金字塔
void amplit_orient(const Mat& image, vector<Mat>& amplit, vector<Mat>& orient, float scale, int nums);
//DOG金字塔特征点检测
void find_scale_space_extrema(const vector<vector<Mat>>& dog_pyr, const vector<vector<Mat>>& gauss_pyr,
vector<KeyPoint>& keypoints);
//DOG金字塔特征点检测，特征点方向细化版
void find_scale_space_extrema1(const vector<vector<Mat>>& dog_pyr, vector<vector<Mat>>& gauss_pyr,
vector<KeyPoint>& keypoints);
//计算特征点的描述子
void calc_sift_descriptors(const vector<vector<Mat>>& gauss_pyr, vector<KeyPoint>& keypoints,
Mat& descriptors, const vector<Mat>& amplit, const vector<Mat>& orient);
//构建空间尺度—主要是为了获取 amplit 和 orient 在使用 GLOH 描述子的时候使用
void build_sar_sift_space(const Mat& image, vector<Mat>& sar_harris_fun, vector<Mat>& amplit, vector<Mat>& orient);
//GLOH 计算一个特征描述子
void calc_gloh_descriptors(const vector<Mat>& amplit, const vector<Mat>& orient, const vector<KeyPoint>& keys, Mat& descriptors);
//特征点检测
void detect(const Mat& image, vector<vector<Mat>>& gauss_pyr, vector<vector<Mat>>& dog_pyr, vector<KeyPoint>& keypoints,
vector<vector<vector<float>>>& all_cell_contrasts,
vector<vector<float>>& average_contrast, vector<vector<int>>& n_cells, vector<int>& num_cell, vector<vector<int>>& available_n,
vector<int>& available_num, vector<KeyPoint>& final_keypoints,
vector<KeyPoint>& Final_keypoints, vector<KeyPoint>& Final_Keypoints);
//特征点描述
void comput_des(const vector<vector<Mat>>& gauss_pyr, vector<KeyPoint>& final_keypoints, const vector<Mat>& amplit,
const vector<Mat>& orient, Mat& descriptors);
private:
int nfeatures; //设定检测的特征点的个数值,如果此值设置为0，则不影响结果
int nOctaveLayers; //每组金字塔中间层数
double contrastThreshold; //对比度阈值（D(x)）
double edgeThreshold; //特征点边缘曲率阈值
double sigma; //高斯尺度空间初始层的尺度
bool double_size; //是否上采样原始图像
;//注意类结束的分号
2、mySift.cpp
#include"mySift.h"
#include<string>
#include<cmath>
#include<iostream> //输入输出
#include<vector> //vector
#include<algorithm>
#include<numeric> //用于容器元素求和
#include<opencv2\\core\\hal\\hal.hpp>
#include<opencv2\\core\\hal\\intrin.hpp>
#include<opencv2\\opencv.hpp>
#include<opencv2\\core\\core.hpp> //opencv基本数据结构
#include<opencv2\\highgui\\highgui.hpp> //图像界面
#include<opencv2\\imgproc\\imgproc.hpp> //基本图像处理函数
#include<opencv2\\features2d\\features2d.hpp> //特征提取
#include<opencv2\\imgproc\\types_c.h>
#include<opencv2\\videoio.hpp>
#include <highgui\\highgui.hpp>
#include <imgproc\\imgproc.hpp>
/******************根据输入图像大小计算高斯金字塔的组数****************************/
/*image表示原始输入灰度图像,inline函数必须在声明处定义
double_size_image表示是否在构建金字塔之前上采样原始图像
*/
int mySift::num_octaves(const Mat& image)
int temp;
float size_temp = (float)min(image.rows, image.cols);
temp = cvRound(log(size_temp) / log((float)2) - 2);
if (double_size)
temp += 1;
if (temp > MAX_OCTAVES) //尺度空间最大组数设置为MAX_OCTAVES
temp = MAX_OCTAVES;
return temp;
/************************sobel滤波器计算高斯尺度空间图像梯度大小*****************************/
void sobelfilter(Mat& image, Mat& G)
// image 是经过归一化后的数据 (0,1)
int rows = image.rows;
int cols = image.cols;
float dx = 0.0, dy = 0.0;
//cv::Mat Gx = cv::Mat::zeros(rows, cols, CV_32FC1); //包含了图像像素在水平方向上的导数的近似值的图像
//cv::Mat Gy = cv::Mat::zeros(rows, cols, CV_32FC1); //包含了图像像素在垂直方向上的导数的近似值的图像
G = Mat::zeros(rows, cols, CV_32FC1); //在每个像素点处的灰度大小由Gx和Gy共同决定
double v = 0.0, vx, vy;
//利用 sobel 算子求梯度幅度图像
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
v = 0.0;
if (i == 0 || j == 0 || i == rows - 1 || j == cols - 1)
G.at<float>(i, j) = 0.0;
else
float dx = image.at<float>(i - 1, j + 1) - image.at<float>(i - 1, j - 1)
+ 2 * image.at<float>(i, j + 1) - 2 * image.at<float>(i, j - 1)
+ image.at<float>(i + 1, j + 1) - image.at<float>(i + 1, j - 1);
float dy = image.at<float>(i + 1, j - 1) - image.at<float>(i - 1, j - 1)
+ 2 * image.at<float>(i + 1, j) - 2 * image.at<float>(i - 1, j) +
image.at<float>(i + 1, j + 1) - image.at<float>(i - 1, j + 1);
v = abs(dx) + abs(dy); //简化后 G = |Gx| + |Gy|
//保证像素值在有效范围内
v = fmax(v, 0); //返回浮点数中较大的一个
v = fmin(v, 255); //返回浮点数中较小的一个
if (v > T) //T为阈值等于50
G.at<float>(i, j) = (float)v;
else
G.at<float>(i, j) = 0.0;
//水平方向上的导数的近似值的图像
/*for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
vx = 0;
if (i == 0 || j == 0 || i == rows - 1 || j == cols - 1)
Gx.at<float>(i, j) = 0;
else
dx = image.at<float>(i - 1, j + 1) - image.at<float>(i - 1, j - 1)
+ 2 * image.at<float>(i, j + 1) - 2 * image.at<float>(i, j - 1)
+ image.at<float>(i + 1, j + 1) - image.at<float>(i + 1, j - 1);
vx = abs(dx);
vx = fmax(vx, 0); vx = fmin(vx, 255);
Gx.at<float>(i, j) = (float)vx;
*/
//垂直方向上的导数的近似值的图像
/*for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
vy = 0;
if (i == 0 || j == 0 || i == rows - 1 || j == cols - 1)
Gy.at<float>(i, j) = 0;
else
dy = image.at<float>(i + 1, j - 1) - image.at<float>(i - 1, j - 1)
+ 2 * image.at<float>(i + 1, j) - 2 * image.at<float>(i - 1, j) +
image.at<float>(i + 1, j + 1) - image.at<float>(i - 1, j + 1);
vy = abs(dy);
vy = fmax(vy, 0); vx = fmin(vy, 255);
Gy.at<float>(i, j) = (float)vy;
*/
//cv::imshow("Gx", Gx); // horizontal
//cv::imshow("Gy", Gy); // vertical
//cv::imshow("G", G); // gradient
/*********该函数根据尺度和窗口半径生成ROEWA滤波模板************/
/*size表示核半径，因此核宽度是2*size+1
scale表示指数权重参数
kernel表示生成的滤波核
*/
static void roewa_kernel(int size, float scale, Mat& kernel)
kernel.create(2 * size + 1, 2 * size + 1, CV_32FC1);
for (int i = -size; i <= size; ++i)
float* ptr_k = kernel.ptr<float>(i + size);
for (int j = -size; j <= size; ++j)
ptr_k[j + size] = exp(-1.f * (abs(i) + abs(j)) / scale);
/************************创建高斯金字塔第一组，第一层图像************************************/
/*image表示输入原始图像
init_image表示生成的高斯尺度空间的第一层图像
*/
void mySift::create_initial_image(const Mat& image, Mat& init_image)
Mat gray_image;
if (image.channels() != 1)
cvtColor(image, gray_image, CV_RGB2GRAY); //转换为灰度图像
else
gray_image = image.clone();
Mat floatImage; //转换到0-1之间的浮点类型数据归一化，方便接下来的处理
//float_image=(float)gray_image*(1.0/255.0)
gray_image.convertTo(floatImage, CV_32FC1, 1.0 / 255.0, 0);
double sig_diff = 0;
if (double_size)
Mat temp_image;
//通过插值的方法改变图像尺寸的大小
resize(floatImage, temp_image, Size(2 * floatImage.cols, 2 * floatImage.rows), 0.0, 0.0, INTER_LINEAR);
//高斯平滑的标准差，值较大时平滑效果比较明显
sig_diff = sqrt(sigma * sigma - 4.0 * INIT_SIGMA * INIT_SIGMA);
//高斯滤波窗口大小选择很重要，这里选择(4*sig_diff_1+1)-(6*sig_diff+1)之间，且四舍五入
int kernel_width = 2 * cvRound(GAUSS_KERNEL_RATIO * sig_diff) + 1;
Size kernel_size(kernel_width, kernel_width);
//对图像进行平滑处理(高斯模糊)，即降低图像的分辨率,高斯模糊是实现尺度变换的唯一变换核，并其实唯一的线性核
GaussianBlur(temp_image, init_image, kernel_size, sig_diff, sig_diff);
else
sig_diff = sqrt(sigma * sigma - 1.0 * INIT_SIGMA * INIT_SIGMA);
//高斯滤波窗口大小选择很重要，这里选择(4*sig_diff_1+1)-(6*sig_diff+1)之间
int kernel_width = 2 * cvRound(GAUSS_KERNEL_RATIO * sig_diff) + 1;
Size kernel_size(kernel_width, kernel_width);
GaussianBlur(floatImage, init_image, kernel_size, sig_diff, sig_diff);
/************************使用 sobel 算子创建高斯金字塔第一组，第一层图像****************************/
//目的是为了减少冗余特征点
void mySift::sobel_create_initial_image(const Mat& image, Mat& init_image)
Mat gray_image, gray_images; //gray_images用于存放经过sobel算子操作后的图像
if (image.channels() != 1)
cvtColor(image, gray_image, CV_RGB2GRAY); //转换为灰度图像
else
gray_image = image.clone();
sobelfilter(gray_image, gray_images);
Mat floatImage; //转换到0-1之间的浮点类型数据归一化，方便接下来的处理
//float_image=(float)gray_image*(1.0/255.0)
gray_images.convertTo(floatImage, CV_32FC1, 1.0 / 255.0, 0);
double sig_diff = 0;
if (double_size)
Mat temp_image;
//通过插值的方法改变图像尺寸的大小
resize(floatImage, temp_image, Size(2 * floatImage.cols, 2 * floatImage.rows), 0.0, 0.0, INTER_LINEAR);
//高斯平滑的标准差，值较大时平滑效果比较明显
sig_diff = sqrt(sigma * sigma - 4.0 * INIT_SIGMA * INIT_SIGMA);
//高斯滤波窗口大小选择很重要，这里选择(4*sig_diff_1+1)-(6*sig_diff+1)之间，且四舍五入
int kernel_width = 2 * cvRound(GAUSS_KERNEL_RATIO * sig_diff) + 1;
Size kernel_size(kernel_width, kernel_width);
//对图像进行平滑处理(高斯模糊)，即降低图像的分辨率,高斯模糊是实现尺度变换的唯一变换核，并其实唯一的线性核
GaussianBlur(temp_image, init_image, kernel_size, sig_diff, sig_diff);
else
sig_diff = sqrt(sigma * sigma - 1.0 * INIT_SIGMA * INIT_SIGMA);
//高斯滤波窗口大小选择很重要，这里选择(4*sig_diff_1+1)-(6*sig_diff+1)之间
int kernel_width = 2 * cvRound(GAUSS_KERNEL_RATIO * sig_diff) + 1;
Size kernel_size(kernel_width, kernel_width);
GaussianBlur(floatImage, init_image, kernel_size, sig_diff, sig_diff);
/**************************生成高斯金字塔*****************************************/
/*init_image表示已经生成的高斯金字塔第一层图像
gauss_pyramid表示生成的高斯金字塔
nOctaves表示高斯金字塔的组数
*/
void mySift::build_gaussian_pyramid(const Mat& init_image, vector<vector<Mat>>& gauss_pyramid, int nOctaves)
vector<double> sig;
sig.push_back(sigma);
double k = pow(2.0, 1.0 / nOctaveLayers); //高斯金字塔每一层的系数 k
for (int i = 1; i < nOctaveLayers + 3; ++i)
double prev_sig = pow(k, i - 1) * sigma; //每一个尺度层的尺度
double curr_sig = k * prev_sig;
//组内每层的尺度坐标计算公式
sig.push_back(sqrt(curr_sig * curr_sig - prev_sig * prev_sig));
gauss_pyramid.resize(nOctaves);
for (int i = 0; i < nOctaves; ++i)
gauss_pyramid[i].resize(nOctaveLayers + 3);
for (int i = 0; i < nOctaves; ++i) //对于每一组
for (int j = 0; j < nOctaveLayers + 3; ++j) //对于组内的每一层
if (i == 0 && j == 0) //第一组，第一层
gauss_pyramid[0][0] = init_image;
else if (j == 0)
resize(gauss_pyramid[i - 1][3], gauss_pyramid[i][0],
Size(gauss_pyramid[i - 1][3].cols / 2,
gauss_pyramid[i - 1][3].rows / 2), 0, 0, INTER_LINEAR);
else
//高斯滤波窗口大小选择很重要，这里选择(4*sig_diff_1+1)-(6*sig_diff+1)之间
int kernel_width = 2 * cvRound(GAUSS_KERNEL_RATIO * sig[j]) + 1;
Size kernel_size(kernel_width, kernel_width);
GaussianBlur(gauss_pyramid[i][j - 1], gauss_pyramid[i][j], kernel_size, sig[j], sig[j]);
/*******************生成高斯差分金字塔，即LOG金字塔*************************/
/*dog_pyramid表示DOG金字塔
gauss_pyramin表示高斯金字塔*/
void mySift::build_dog_pyramid(vector<vector<Mat>>& dog_pyramid, const vector<vector<Mat>>& gauss_pyramid)
vector<vector<Mat>>::size_type nOctaves = gauss_pyramid.size();
for (vector<vector<Mat>>::size_type i = 0; i < nOctaves; ++i)
//用于存放每一个梯度中的所有尺度层
vector<Mat> temp_vec;
for (auto j = 0; j < nOctaveLayers + 2; ++j)
Mat temp_img = gauss_pyramid[i][j + 1] - gauss_pyramid[i][j];
temp_vec.push_back(temp_img);
dog_pyramid.push_back(temp_vec);
temp_vec.clear();
/***********生成高斯差分金字塔当前层对应的梯度幅度图像和梯度方向图像***********/
/*image为高斯差分金字塔当前层图像
*amplit为当前层梯度幅度图像
*orient为当前层梯度方向图像
*scale当前层尺度
*nums为相对底层的层数
*/
void mySift::amplit_orient(const Mat& image, vector<Mat>& amplit, vector<Mat>& orient, float scale, int nums)
//分配内存
amplit.resize(Mmax * nOctaveLayers);
orient.resize(Mmax * nOctaveLayers);
int radius = cvRound(2 * scale);
Mat kernel; //kernel(2 * radius + 1, 2 * radius + 1, CV_32FC1);
roewa_kernel(radius, scale, kernel); //返回滤波核，也即指数部分,存放在矩阵的右下角
//四个滤波模板生成
Mat W34 = Mat::zeros(2 * radius + 1, 2 * radius + 1, CV_32FC1); //把kernel矩阵下半部分复制到对应部分
Mat W12 = Mat::zeros(2 * radius + 1, 2 * radius + 1, CV_32FC1); //把kernel矩阵上半部分复制到对应部分
Mat W14 = Mat::zeros(2 * radius + 1, 2 * radius + 1, CV_32FC1); //把kernel矩阵右半部分复制到对应部分
Mat W23 = Mat::zeros(2 * radius + 1, 2 * radius + 1, CV_32FC1); //把kernel矩阵左半部分复制到对应部分
kernel(Range(radius + 1, 2 * radius + 1), Range::all()).copyTo(W34(Range(radius + 1, 2 * radius + 1), Range::all()));
kernel(Range(0, radius), Range::all()).copyTo(W12(Range(0, radius), Range::all()));
kernel(Range::all(), Range(radius + 1, 2 * radius + 1)).copyTo(W14(Range::all(), Range(radius + 1, 2 * radius + 1)));
kernel(Range::all(), Range(0, radius)).copyTo(W23(Range::all(), Range(0, radius)));
//滤波
Mat M34, M12, M14, M23;
double eps = 0.00001;
//float_image为图像归一化后的图像数据，做卷积运算
filter2D(image, M34, CV_32FC1, W34, Point(-1, -1), eps);
filter2D(image, M12, CV_32FC1, W12, Point(-1, -1), eps);
filter2D(image, M14, CV_32FC1, W14, Point(-1, -1), eps);
filter2D(image, M23, CV_32FC1, W23, Point(-1, -1), eps);
//计算水平梯度和竖直梯度
Mat Gx, Gy;
log((M14) / (M23), Gx);
log((M34) / (M12), Gy);
//计算梯度幅度和梯度方向
magnitude(Gx, Gy, amplit[nums]); //梯度幅度图像，平方和开平方
phase(Gx, Gy, orient[nums], true); //梯度方向图像
/***********************该函数计算尺度空间特征点的主方向，用于后面特征点的检测***************************/
/*image表示特征点所在位置的高斯图像，后面可对着源码进行修改
pt表示特征点的位置坐标(x,y)
scale特征点的尺度
n表示直方图bin个数
hist表示计算得到的直方图
函数返回值是直方图hist中的最大数值*/
static float clac_orientation_hist(const Mat& image, Point pt, float scale, int n, float* hist)
int radius = cvRound(ORI_RADIUS * scale); //特征点邻域半径(3*1.5*scale)
int len = (2 * radius + 1) * (2 * radius + 1); //特征点邻域像素总个数（最大值）
float sigma = ORI_SIG_FCTR * scale; //特征点邻域高斯权重标准差(1.5*scale)
float exp_scale = -1.f / (2 * sigma * sigma); //权重的指数部分
//使用AutoBuffer分配一段内存，这里多出4个空间的目的是为了方便后面平滑直方图的需要
AutoBuffer<float> buffer((4 * len) + n + 4);
//X保存水平差分，Y保存数值差分，Mag保存梯度幅度，Ori保存梯度角度，W保存高斯权重
float* X = buffer, * Y = buffer + len, * Mag = Y, * Ori = Y + len, * W = Ori + len;
float* temp_hist = W + len + 2; //临时保存直方图数据
for (int i = 0; i < n; ++i)
temp_hist[i] = 0.f; //数据清零
//计算邻域像素的水平差分和竖直差分
int k = 0;
for (int i = -radius; i < radius; ++i)
int y = pt.y + i; //邻域点的纵坐标
if (y <= 0 || y >= image.rows - 1)
continue;
for (int j = -radius; j < radius; ++j)
int x = pt.x + j;
if (x <= 0 || x >= image.cols - 1)
continue;
float dx = image.at<float>(y, x + 1) - image.at<float>(y, x - 1); //水平差分
float dy = image.at<float>(y + 1, x) - image.at<float>(y - 1, x); //竖直差分
//保存水平差分和竖直差分及其对应的权重
X[k] = dx;
Y[k] = dy;
W[k] = (i * i + j * j) * exp_scale;
++k;
len = k; //邻域内特征点的个数
cv::hal::exp(W, W, len); //计算邻域内所有像素的高斯权重
cv::hal::fastAtan2(Y, X, Ori, len, true); //计算邻域内所有像素的梯度方向，角度范围0-360度
cv::hal::magnitude32f(X, Y, Mag, len); //计算邻域内所有像素的梯度幅度，计算的是数学意义上的梯度
//遍历邻域的像素
for (int i = 0; i < len; ++i)
int bin = cvRound((n / 360.f) * Ori[i]); //利用像素的梯度方向，约束bin的范围在[0,(n-1)]
//像素点梯度方向为360度时，和0°一样
if (bin >= n)
bin = bin - n;
if (bin < 0)
bin = bin + n;
temp_hist[bin] = temp_hist[bin] + Mag[i] * W[i]; //统计邻域内像素各个方向在梯度直方图的幅值(加权后的幅值)
//平滑直方图
temp_hist[-1] = temp_hist[n - 1];
temp_hist[-2] = temp_hist[n - 2];
temp_hist[n] = temp_hist[0];
temp_hist[n + 1] = temp_hist[1];
for (int i = 0; i < n; ++i)
hist[i] = (temp_hist[i - 2] + temp_hist[i + 2]) * (1.f / 16.f) +
(temp_hist[i - 1] + temp_hist[i + 1]) * (4.f / 16.f) +
temp_hist[i] * (6.f / 16.f);
//获得直方图中最大值
float max_value = hist[0];
for (int i = 1; i < n; ++i)
if (hist[i] > max_value)
max_value = hist[i];
return max_value;
/***********************使用 sobel 滤波器定义的新梯度计算尺度空间特征点的主方向**************************/
static float clac_orientation_hist_2(Mat& image, Point pt, float scale, int n, float* hist)
Mat output_image; //使用 sobel 滤波器计算的图像的梯度幅度图像
sobelfilter(image, output_image); //使用 sobel 滤波器求高斯差分图像的梯度幅度图像
int radius = cvRound(ORI_RADIUS * scale); //特征点邻域半径(3*1.5*scale)
int len = (2 * radius + 1) * (2 * radius + 1); //特征点邻域像素总个数（最大值）
float sigma = ORI_SIG_FCTR * scale; //特征点邻域高斯权重标准差(1.5*scale)
float exp_scale = -1.f / (2 * sigma * sigma); //权重的指数部分
//使用AutoBuffer分配一段内存，这里多出4个空间的目的是为了方便后面平滑直方图的需要
AutoBuffer<float> buffer((4 * len) + n + 4);
//X保存水平差分，Y保存数值差分，Mag保存梯度幅度，Ori保存梯度角度，W保存高斯权重
float* X = buffer, * Y = buffer + len, * Mag = Y, * Ori = Y + len, * W = Ori + len;
float* temp_hist = W + len + 2; //临时保存直方图数据
for (int i = 0; i < n; ++i)
temp_hist[i] = 0.f; //数据清零
//计算邻域像素的水平差分和竖直差分
int k = 0;
for (int i = -radius; i < radius; ++i)
int y = pt.y + i; //邻域点的纵坐标，行
if (y <= 0 || y >= output_image.rows - 1)
continue;
for (int j = -radius; j < radius; ++j)
int x = pt.x + j; //邻域点的纵坐标，列
if (x <= 0 || x >= output_image.cols - 1)
continue;
//float dx = image.at<float>(y, x + 1) - image.at<float>(y, x - 1); //水平差分
float dx = output_image.at<float>(y - 1, x + 1) - output_image.at<float>(y - 1, x - 1)
+ 2 * output_image.at<float>(y, x + 1) - 2 * output_image.at<float>(y, x - 1)
+ output_image.at<float>(y + 1, x + 1) - output_image.at<float>(y + 1, x - 1);
float dy = output_image.at<float>(y + 1, x - 1) - output_image.at<float>(y - 1, x - 1)
+ 2 * output_image.at<float>(y + 1, x) - 2 * output_image.at<float>(y - 1, x) +
output_image.at<float>(y + 1, x + 1) - output_image.at<float>(y - 1, x + 1);
/*float dx = image.at<float>(y - 1, x + 1) - image.at<float>(y - 1, x - 1)
+ 2 * image.at<float>(y, x + 1) - 2 * image.at<float>(y, x - 1)
+ image.at<float>(y + 1, x + 1) - image.at<float>(y + 1, x - 1);
float dy = image.at<float>(y + 1, x - 1) - image.at<float>(y - 1, x - 1)
+ 2 * image.at<float>(y + 1, x) - 2 * image.at<float>(y - 1, x) +
image.at<float>(y + 1, x + 1) - image.at<float>(y - 1, x + 1);*/
//保存水平差分和竖直差分及其对应的权重
X[k] = dx;
Y[k] = dy;
W[k] = (i * i + j * j) * exp_scale;
++k;
len = k; //邻域内特征点的个数
cv::hal::exp(W, W, len); //计算邻域内所有像素的高斯权重
cv::hal::fastAtan2(Y, X, Ori, len, true); //计算邻域内所有像素的梯度方向，角度范围0-360度
cv::hal::magnitude32f(X, Y, Mag, len); //计算邻域内所有像素的梯度幅度，计算的是数学意义上的梯度
//遍历邻域的像素
for (int i = 0; i < len; ++i)
int bin = cvRound((n / 360.f) * Ori[i]); //利用像素的梯度方向，约束bin的范围在[0,(n-1)]
//像素点梯度方向为360度时，和0°一样
if (bin >= n)
bin = bin - n;
if (bin < 0)
bin = bin + n;
temp_hist[bin] = temp_hist[bin] + Mag[i] * W[i]; //统计邻域内像素各个方向在梯度直方图的幅值(加权后的幅值)
//平滑直方图
temp_hist[-1] = temp_hist[n - 1];
temp_hist[-2] = temp_hist[n - 2];
temp_hist[n] = temp_hist[0];
temp_hist[n + 1] = temp_hist[1];
for (int i = 0; i < n; ++i)
hist[i] = (temp_hist[i - 2] + temp_hist[i + 2]) * (1.f / 16.f) +
(temp_hist[i - 1] + temp_hist[i + 1]) * (4.f / 16.f) +
temp_hist[i] * (6.f / 16.f);
//获得直方图中最大值
float max_value = hist[0];
for (int i = 1; i < n; ++i)
if (hist[i] > max_value)
max_value = hist[i];
return max_value;
/******************该函数计算特征点主方向,用于LOGH版本*********************/
/*amplit表示特征点所在层的梯度幅度，即输入图像对应像素点的梯度存在了对应位置
orient表示特征点所在层的梯度方向，0-360度
point表示特征点坐标
scale表示特征点的所在层的尺度
hist表示生成的直方图
n表示主方向直方图bin个数
*/
static float calc_orient_hist(const Mat& amplit, const Mat& orient, Point pt, float scale, float* hist, int n)
//暂且认为是只进行了下采样，没有进行高斯模糊
int num_row = amplit.rows;
int num_col = amplit.cols;
Point point(cvRound(pt.x), cvRound(pt.y));
//int radius = cvRound(SAR_SIFT_FACT_RADIUS_ORI * scale);
int radius = cvRound(6 * scale);
radius = min(radius, min(num_row / 2, num_col / 2));
float gauss_sig = 2 * scale; //高斯加权标准差
float exp_temp = -1.f / (2 * gauss_sig * gauss_sig); //权重指数部分
//邻域区域
int radius_x_left = point.x - radius;
int radius_x_right = point.x + radius;
int radius_y_up = point.y - radius;
int radius_y_down = point.y + radius;
//防止越界
if (radius_x_left < 0)
radius_x_left = 0;
if (radius_x_right > num_col - 1)
radius_x_right = num_col - 1;
if (radius_y_up < 0)
radius_y_up = 0;
if (radius_y_down > num_row - 1)
radius_y_down = num_row - 1;
//此时特征点周围矩形区域相对于本矩形区域的中心
int center_x = point.x - radius_x_left;
int center_y = point.y - radius_y_up;
//矩形区域的边界，计算高斯权值
Range x_rng(-(point.x - radius_x_left), radius_x_right - point.x);
Range y_rng(-(point.y - radius_y_up), radius_y_down - point.y);
Mat gauss_weight;
gauss_weight.create(y_rng.end - y_rng.start + 1, x_rng.end - x_rng.start + 1, CV_32FC1);
//求各个像素点的高斯权重
for (int i = y_rng.start; i <= y_rng.end; ++i)
float* ptr_gauss = gauss_weight.ptr<float>(i - y_rng.start);
for (int j = x_rng.start; j <= x_rng.end; ++j)
ptr_gauss[j - x_rng.start] = exp((i * i + j * j) * exp_temp);
//索引特征点周围的像素梯度幅度，梯度方向
Mat sub_amplit = amplit(Range(radius_y_up, radius_y_down + 1), Range(radius_x_left, radius_x_right + 1));
Mat sub_orient = orient(Range(radius_y_up, radius_y_down + 1), Range(radius_x_left, radius_x_right + 1));
//Mat W = sub_amplit.mul(gauss_weight); //加入高斯权重，计算高斯权重时，正确匹配点对反而变少了
Mat W = sub_amplit; //没加高斯权重，梯度幅值
//计算直方图
AutoBuffer<float> buffer(n + 4);
float* temp_hist = buffer + 2;
for (int i = 0; i < n; ++i)
temp_hist[i] = 0.f;
for (int i = 0; i < sub_orient.rows; i++)
float* ptr_1 = W.ptr<float>(i);
float* ptr_2 = sub_orient.ptr<float>(i);
for (int j = 0; j < sub_orient.cols; j++)
if (((i - center_y) * (i - center_y) + (j - center_x) * (j - center_x)) < radius * radius)
int bin = cvRound(ptr_2[j] * n / 360.f);
if (bin > n)
bin = bin - n;
if (bin < 0)
bin = bin + n;
temp_hist[bin] += ptr_1[j];
//平滑直方图，可以防止突变
temp_hist[-1] = temp_hist[n - 1];
temp_hist[-2] = temp_hist[n - 2];
temp_hist[n] = temp_hist[0];
temp_hist[n + 1] = temp_hist[1];
for (int i = 0; i < n; ++i)
hist[i] = (temp_hist[i - 2] + temp_hist[i + 2]) * (1.f / 16.f) +
(temp_hist[i - 1] + temp_hist[i + 1]) * (4.f / 16.f) +
temp_hist[i] * (6.f / 16.f);
//获得直方图中最大值
float max_value = hist[0];
for (int i = 1; i < n; ++i)
if (hist[i] > max_value)
max_value = hist[i];
return max_value;
/****************************该函数精确定位特征点位置(x,y,scale)，用于后面特征点的检测*************************/
/*功能：确定特征点的位置，并通过主曲率消除边缘相应点,该版本是简化版
dog_pry表示DOG金字塔
kpt表示精确定位后该特征点的信息
octave表示初始特征点所在的组
layer表示初始特征点所在的层
row表示初始特征点在图像中的行坐标
col表示初始特征点在图像中的列坐标
nOctaveLayers表示DOG金字塔每组中间层数，默认是3
contrastThreshold表示对比度阈值，默认是0.04
edgeThreshold表示边缘阈值，默认是10
sigma表示高斯尺度空间最底层图像尺度，默认是1.6*/
static bool adjust_local_extrema_1(const vector<vector<Mat>>& dog_pyr, KeyPoint& kpt, int octave, int& layer, int& row,
int& col, int nOctaveLayers, float contrastThreshold, float edgeThreshold, float sigma)
float xi = 0, xr = 0, xc = 0;
int i = 0;
for (; i < MAX_INTERP_STEPS; ++i) //最大迭代次数
const Mat& img = dog_pyr[octave][layer]; //当前层图像索引
const Mat& prev = dog_pyr[octave][layer - 1]; //之前层图像索引
const Mat& next = dog_pyr[octave][layer + 1]; //下一层图像索引
//特征点位置x方向，y方向,尺度方向的一阶偏导数
float dx = (img.at<float>(row, col + 1) - img.at<float>(row, col - 1)) * (1.f / 2.f);
float dy = (img.at<float>(row + 1, col) - img.at<float>(row - 1, col)) * (1.f / 2.f);
float dz = (next.at<float>(row, col) - prev.at<float>(row, col)) * (1.f / 2.f);
//计算特征点位置二阶偏导数
float v2 = img.at<float>(row, col);
float dxx = img.at<float>(row, col + 1) + img.at<float>(row, col - 1) - 2 * v2;
float dyy = img.at<float>(row + 1, col) + img.at<float>(row - 1, col) - 2 * v2;
float dzz = prev.at<float>(row, col) + next.at<float>(row, col) - 2 * v2;
//计算特征点周围混合二阶偏导数
float dxy = (img.at<float>(row + 1, col + 1) + img.at<float>(row - 1, col - 1) -
img.at<float>(row + 1, col - 1) - img.at<float>(row - 1, col + 1)) * (1.f / 4.f);
float dxz = (next.at<float>(row, col + 1) + prev.at<float>(row, col - 1) -
next.at<float>(row, col - 1) - prev.at<float>(row, col + 1)) * (1.f / 4.f);
float dyz = (next.at<float>(row + 1, col) + prev.at<float>(row - 1, col) -
next.at<float>(row - 1, col) - prev.at<float>(row + 1, col)) * (1.f / 4.f);
Matx33f H(dxx, dxy, dxz, dxy, dyy, dyz, dxz, dyz, dzz);
Vec3f dD(dx, dy, dz);
Vec3f X = H.solve(dD, DECOMP_SVD);
xc = -X[0]; //x方向偏移量
xr = -X[1]; //y方向偏移量
xi = -X[2]; //尺度方向偏移量
//如果三个方向偏移量都小于0.5，说明已经找到特征点准确位置
if (abs(xc) < 0.5f && abs(xr) < 0.5f && abs(xi) < 0.5f)
break;
//如果其中一个方向的偏移量过大，则删除该点
if (abs(xc) > (float)(INT_MAX / 3) ||
abs(xr) > (float)(INT_MAX / 3) ||
abs(xi) > (float)(INT_MAX / 3))
return false;
col = col + cvRound(xc);
row = row + cvRound(xr);
layer = layer + cvRound(xi);
//如果特征点定位在边界区域，同样也需要删除
if (layer<1 || layer>nOctaveLayers ||
col<IMG_BORDER || col>img.cols - IMG_BORDER ||
row<IMG_BORDER || row>img.rows - IMG_BORDER)
return false;
//如果i=MAX_INTERP_STEPS，说明循环结束也没有满足条件，即该特征点需要被删除
if (i >= MAX_INTERP_STEPS)
return false;
/**************************再次删除低响应点(对比度较低的点)********************************/
//再次计算特征点位置x方向，y方向,尺度方向的一阶偏导数
//高对比度的特征对图像的变形是稳定的
const Mat& img = dog_pyr[octave][layer];
const Mat& prev = dog_pyr[octave][layer - 1];
const Mat& next = dog_pyr[octave][layer + 1];
float dx = (img.at<float>(row, col + 1) - img.at<float>(row, col - 1)) * (1.f / 2.f);
float dy = (img.at<float>(row + 1, col) - img.at<float>(row - 1, col)) * (1.f / 2.f);
float dz = (next.at<float>(row, col) - prev.at<float>(row, col)) * (1.f / 2.f);
Matx31f dD(dx, dy, dz);
float t = dD.dot(Matx31f(xc, xr, xi));
float contr = img.at<float>(row, col) + t * 0.5f; //特征点响应 |D(x~)| 即对比度
//Low建议contr阈值是0.03，但是RobHess等建议阈值为0.04/nOctaveLayers
if (abs(contr) < contrastThreshold / nOctaveLayers) //阈值设为0.03时特征点数量过多
return false;
/******************************删除边缘响应比较强的点************************************/
//再次计算特征点位置二阶偏导数，获取特征点出的 Hessian 矩阵，主曲率通过 2X2 的 Hessian 矩阵求出
//一个定义不好的高斯差分算子的极值在横跨边缘的地方有较大的主曲率而在垂直边缘的方向有较小的主曲率
float v2 = img.at<float>(row, col);
float dxx = img.at<float>(row, col + 1) + img.at<float>(row, col - 1) - 2 * v2;
float dyy = img.at<float>(row + 1, col) + img.at<float>(row - 1, col) - 2 * v2;
float dxy = (img.at<float>(row + 1, col + 1) + img.at<float>(row - 1, col - 1) -
img.at<float>(row + 1, col - 1) - img.at<float>(row - 1, col + 1)) * (1.f / 4.f);
float det = dxx * dyy - dxy * dxy;
float trace = dxx + dyy;
//主曲率和阈值的对比判定
if (det < 0 || (trace * trace * edgeThreshold >= det * (edgeThreshold + 1) * (edgeThreshold + 1)))
return false;
/*********到目前为止该特征的满足上面所有要求，保存该特征点信息***********/
kpt.pt.x = ((float)col + xc) * (1 << octave); //相对于最底层的图像的x坐标
kpt.pt.y = ((float)row + xr) * (1 << octave); //相对于最底层图像的y坐标
kpt.octave = octave + (layer << 8); //组号保存在低字节，层号保存在高字节
//相对于最底层图像的尺度
kpt.size = sigma * powf(2.f, (layer + xi) / nOctaveLayers) * (1 << octave);
kpt.response = abs(contr); //特征点响应值(对比度)
return true;
/****************************该函数精确定位特征点位置(x,y,scale)，用于后面特征点的检测*************************/
//该版本是 SIFT 原版，检测得到的特征点数量更多
static bool adjust_local_extrema_2(const vector<vector<Mat>>& dog_pyr, KeyPoint& kpt, int octave, int& layer, int& row,
int& col, int nOctaveLayers, float contrastThreshold, float edgeThreshold, float sigma)
const float img_scale = 1.f / (255 * SIFT_FIXPT_SCALE); //SIFT_FIXPT_SCALE=48
const float deriv_scale = img_scale * 0.5f;
const float second_deriv_scale = img_scale;
const float cross_deriv_scale = img_scale * 0.25f;
float xi = 0, xr = 0, xc = 0;
int i = 0;
for (; i < MAX_INTERP_STEPS; ++i) //最大迭代次数
const Mat& img = dog_pyr[octave][layer]; //当前层图像索引
const Mat& prev = dog_pyr[octave][layer - 1]; //之前层图像索引
const Mat& next = dog_pyr[octave][layer + 1]; //下一层图像索引
//计算一阶偏导数，通过临近点差分求得
float dx = (img.at<float>(row, col + 1) - img.at<float>(row, col - 1)) * deriv_scale;
float dy = (img.at<float>(row + 1, col) - img.at<float>(row - 1, col)) * deriv_scale;
float dz = (next.at<float>(row, col) - prev.at<float>(row, col)) * deriv_scale;
//计算特征点位置二阶偏导数
//float v2 = img.at<float>(row, col);
float v2 = (float)img.at<float>(row, col) * 2.f;
float dxx = (img.at<float>(row, col + 1) + img.at<float>(row, col - 1) - v2) * second_deriv_scale;
float dyy = (img.at<float>(row + 1, col) + img.at<float>(row - 1, col) - v2) * second_deriv_scale;
float dzz = (prev.at<float>(row, col) + next.at<float>(row, col) - v2) * second_deriv_scale;
//计算特征点周围混合二阶偏导数
float dxy = (img.at<float>(row + 1, col + 1) + img.at<float>(row - 1, col - 1) -
img.at<float>(row + 1, col - 1) - img.at<float>(row - 1, col + 1)) * cross_deriv_scale;
float dxz = (next.at<float>(row, col + 1) + prev.at<float>(row, col - 1) -
next.at<float>(row, col - 1) - prev.at<float>(row, col + 1)) * cross_deriv_scale;
float dyz = (next.at<float>(row + 1, col) + prev.at<float>(row - 1, col) -
next.at<float>(row - 1, col) - prev.at<float>(row + 1, col)) * cross_deriv_scale;
Matx33f H(dxx, dxy, dxz, dxy, dyy, dyz, dxz, dyz, dzz);
Vec3f dD(dx, dy, dz);
Vec3f X = H.solve(dD, DECOMP_SVD);
xc = -X[0]; //x方向偏移量
xr = -X[1]; //y方向偏移量
xi = -X[2]; //尺度方向偏移量
//如果三个方向偏移量都小于0.5，说明已经找到特征点准确位置
if (abs(xc) < 0.5f && abs(xr) < 0.5f && abs(xi) < 0.5f)
break;
//如果其中一个方向的偏移量过大，则删除该点
if (abs(xc) > (float)(INT_MAX / 3) ||
abs(xr) > (float)(INT_MAX / 3) ||
abs(xi) > (float)(INT_MAX / 3))
return false;
col = col + cvRound(xc);
row = row + cvRound(xr);
layer = layer + cvRound(xi);
//如果特征点定位在边界区域，同样也需要删除
if (layer<1 || layer>nOctaveLayers ||
col < IMG_BORDER || col >= img.cols - IMG_BORDER ||
row < IMG_BORDER || row >= img.rows - IMG_BORDER)
return false;
//如果i=MAX_INTERP_STEPS，说明循环结束也没有满足条件，即该特征点需要被删除
if (i >= MAX_INTERP_STEPS)
return false;
/**************************再次删除低响应点(对比度较低的点)********************************/
//再次计算特征点位置x方向，y方向,尺度方向的一阶偏导数
//高对比度的特征对图像的变形是稳定的
const Mat& img = dog_pyr[octave][layer];
const Mat& prev = dog_pyr[octave][layer - 1];
const Mat& next = dog_pyr[octave][layer + 1];
float dx = (img.at<float>(row, col + 1) - img.at<float>(row, col - 1)) * deriv_scale;
float dy = (img.at<float>(row + 1, col) - img.at<float>(row - 1, col)) * deriv_scale;
float dz = (next.at<float>(row, col) - prev.at<float>(row, col)) * deriv_scale;
Matx31f dD(dx, dy, dz);
float t = dD.dot(Matx31f(xc, xr, xi));
float contr = img.at<float>(row, col) + t * 0.5f; //特征点响应 |D(x~)| 即对比度
//Low建议contr阈值是0.03，但是RobHess等建议阈值为0.04/nOctaveLayers
if (abs(contr) < contrastThreshold / nOctaveLayers) //阈值设为0.03时特征点数量过多
return false;
/******************************删除边缘响应比较强的点************************************/
//再次计算特征点位置二阶偏导数，获取特征点出的 Hessian 矩阵，主曲率通过 2X2 的 Hessian 矩阵求出
//一个定义不好的高斯差分算子的极值在横跨边缘的地方有较大的主曲率而在垂直边缘的方向有较小的主曲率
float v2 = (float)img.at<float>(row, col) * 2.f;
float dxx = (img.at<float>(row, col + 1) + img.at<float>(row, col - 1) - v2) * second_deriv_scale;
float dyy = (img.at<float>(row + 1, col) + img.at<float>(row - 1, col) - v2) * second_deriv_scale;
float dxy = (img.at<float>(row + 1, col + 1) + img.at<float>(row - 1, col - 1) -
img.at<float>(row + 1, col - 1) - img.at<float>(row - 1, col + 1)) * cross_deriv_scale;
float det = dxx * dyy - dxy * dxy;
float trace = dxx + dyy;
//主曲率和阈值的对比判定
if (det <= 0 || (trace * trace * edgeThreshold >= det * (edgeThreshold + 1) * (edgeThreshold + 1)))
return false;
/*********保存该特征点信息***********/
kpt.pt.x = ((float)col + xc) * (1 << octave); //高斯金字塔坐标根据组数扩大相应的倍数
kpt.pt.y = ((float)row + xr) * (1 << octave);
// SIFT 描述子
kpt.octave = octave + (layer << 8) + (cvRound((xi + 0.5) * 255) << 16); //特征点被检测出时所在的金字塔组
kpt.size = sigma * powf(2.f, (layer + xi) / nOctaveLayers) * (1 << octave) * 2; //关键点邻域直径
kpt.response = abs(contr); //特征点响应值(对比度)
return true;
/************该函数在DOG金字塔上进行特征点检测，特征点精确定位，删除低对比度点，删除边缘响应较大点**********/
/*dog_pyr表示高斯金字塔 原始 SIFT 算子
gauss_pyr表示高斯金字塔
keypoints表示检测到的特征点*/
void mySift::find_scale_space_extrema(const vector<vector<Mat>>& dog_pyr, const vector<vector<Mat>>& gauss_pyr,
vector<KeyPoint>& keypoints)
int nOctaves = (int)dog_pyr.size(); //子八度个数
//Low文章建议thresholdThreshold是0.03，Rob Hess等人使用0.04/nOctaveLayers作为阈值
float threshold = (float)(contrastThreshold / nOctaveLayers);
const int n = ORI_HIST_BINS; //n=36
float hist[n];
KeyPoint kpt;
keypoints.clear(); //先清空keypoints
for (int i = 0; i < nOctaves; ++i) //对于每一组
for (int j = 1; j <= nOctaveLayers; ++j) //对于组内每一层
const Mat& curr_img = dog_pyr[i][j]; //当前层
const Mat& prev_img = dog_pyr[i][j - 1]; //上一层
const Mat& next_img = dog_pyr[i][j + 1]; //下一层
int num_row = curr_img.rows;
int num_col = curr_img.cols; //获得当前组图像的大小
size_t step = curr_img.step1(); //一行元素所占字节数
//遍历每一个尺度层中的有效像素，像素值
for (int r = IMG_BORDER; r < num_row - IMG_BORDER; ++r)
const float* curr_ptr = curr_img.ptr<float>(r); //指向的是第 r 行的起点，返回的是 float 类型的像素值
const float* prev_ptr = prev_img.ptr<float>(r - 1);
const float* next_ptr = next_img.ptr<float>(r + 1);
for (int c = IMG_BORDER; c < num_col - IMG_BORDER; ++c)
float val = curr_ptr[c]; //当前中心点响应值
//开始检测特征点
if (abs(val) > threshold &&
((val > 0 && val >= curr_ptr[c - 1] && val >= curr_ptr[c + 1] &&
val >= curr_ptr[c - step - 1] && val >= curr_ptr[c - step] && val >= curr_ptr[c - step + 1] &&
val >= curr_ptr[c + step - 1] && val >= curr_ptr[c + step] && val >= curr_ptr[c + step + 1] &&
val >= prev_ptr[c] && val >= prev_ptr[c - 1] && val >= prev_ptr[c + 1] &&
val >= prev_ptr[c - step - 1] && val >= prev_ptr[c - step] && val >= prev_ptr[c - step + 1] &&
val >= prev_ptr[c + step - 1] && val >= prev_ptr[c + step] && val >= prev_ptr[c + step + 1] &&
val >= next_ptr[c] && val >= next_ptr[c - 1] && val >= next_ptr[c + 1] &&
val >= next_ptr[c - step - 1] && val >= next_ptr[c - step] && val >= next_ptr[c - step + 1] &&
val >= next_ptr[c + step - 1] && val >= next_ptr[c + step] && val >= next_ptr[c + step + 1]) ||
(val < 0 && val <= curr_ptr[c - 1] && val <= curr_ptr[c + 1] &&
val <= curr_ptr[c - step - 1] && val <= curr_ptr[c - step] && val <= curr_ptr[c - step + 1] &&
val <= curr_ptr[c + step - 1] && val <= curr_ptr[c + step] && val <= curr_ptr[c + step + 1] &&
val <= prev_ptr[c] && val <= prev_ptr[c - 1] && val <= prev_ptr[c + 1] &&
val <= prev_ptr[c - step - 1] && val <= prev_ptr[c - step] && val <= prev_ptr[c - step + 1] &&
val <= prev_ptr[c + step - 1] && val <= prev_ptr[c + step] && val <= prev_ptr[c + step + 1] &&
val <= next_ptr[c] && val <= next_ptr[c - 1] && val <= next_ptr[c + 1] &&
val <= next_ptr[c - step - 1] && val <= next_ptr[c - step] && val <= next_ptr[c - step + 1] &&
val <= next_ptr[c + step - 1] && val <= next_ptr[c + step] && val <= next_ptr[c + step + 1])))
//++numKeys;
//获得特征点初始行号，列号，组号，组内层号
int octave = i, layer = j, r1 = r, c1 = c;
if (!adjust_local_extrema_1(dog_pyr, kpt, octave, layer, r1, c1,
nOctaveLayers, (float)contrastThreshold,
(float)edgeThreshold, (float)sigma))
continue; //如果该初始点不满足条件，则不保存改点
float scale = kpt.size / float(1 << octave); //该特征点相对于本组的尺度
//max_hist值对应的方向为主方向
float max_hist = clac_orientation_hist(gauss_pyr[octave][layer], Point(c1, r1), scale, n, hist);
//大于mag_thr值对应的方向为辅助方向
float mag_thr = max_hist * ORI_PEAK_RATIO; //主峰值 80% 的方向作为辅助方向
//遍历直方图中的 36 个bin
for (int i = 0; i < n; ++i)
int left = i > 0 ? i - 1 : n - 1;
int right = i < n - 1 ? i + 1 : 0;
//创建新的特征点，大于主峰值 80% 的方向，赋值给该特征点，作为一个新的特征点；即有多个特征点，位置、尺度相同，方向不同
if (hist[i] > hist[left] && hist[i] > hist[right] && hist[i] >= mag_thr)
float bin = i + 0.5f * (hist[left] - hist[right]) / (hist[left] + hist[right] - 2 * hist[i]);
bin = bin < 0 ? n + bin : bin >= n ? bin - n : bin;
kpt.angle = (360.f / n) * bin; //原始 SIFT 算子使用的特征点的主方向0-360度
keypoints.push_back(kpt); //保存该特征点
//cout << "初始满足要求特征点个数是: " << numKeys << endl;
/************该函数在DOG金字塔上进行特征点检测，特征点精确定位，删除低对比度点，删除边缘响应较大点**********/
//对特征点进行方向的细化 + 增加更多的主方向版本 ——此时细化是对最后要给关键点进行赋值时的细化
//还可以考虑直接对方向直方图进行细化
void mySift::find_scale_space_extrema1(const vector<vector<Mat>>& dog_pyr, vector<vector<Mat>>& gauss_pyr,
vector<KeyPoint>& keypoints)
int nOctaves = (int)dog_pyr.size(); //子八度个数
//Low文章建议thresholdThreshold是0.03，Rob Hess等人使用0.04/nOctaveLayers作为阈值
float threshold = (float)(contrastThreshold / nOctaveLayers);
const int n = ORI_HIST_BINS; //n=36
float hist[n];
KeyPoint kpt;
vector<Mat> amplit; //存放高斯差分金字塔每一层的梯度幅度图像
vector<Mat> orient; //存放高斯差分金字塔每一层的梯度方向图像
keypoints.clear(); //先清空keypoints
for (int i = 0; i < nOctaves; ++i) //对于每一组
for (int j = 1; j <= nOctaveLayers; ++j) //对于组内每一层
const Mat& curr_img = dog_pyr[i][j]; //当前层
const Mat& prev_img = dog_pyr[i][j - 1]; //上一层
const Mat& next_img = dog_pyr[i][j + 1]; //下一层
int num_row = curr_img.rows;
int num_col = curr_img.cols; //获得当前组图像的大小
size_t step = curr_img.step1(); //一行元素所占字节数
//遍历每一个尺度层中的有效像素，像素值
for (int r = IMG_BORDER; r < num_row - IMG_BORDER; ++r)
const float* curr_ptr = curr_img.ptr<float>(r); //指向的是第 r 行的起点，返回的是 float 类型的像素值
const float* prev_ptr = prev_img.ptr<float>(r);
const float* next_ptr = next_img.ptr<float>(r);
for (int c = IMG_BORDER; c < num_col - IMG_BORDER; ++c)
float val = curr_ptr[c]; //当前中心点响应值
//开始检测特征点
if (abs(val) > threshold &&
((val > 0 && val >= curr_ptr[c - 1] && val >= curr_ptr[c + 1] &&
val >= curr_ptr[c - step - 1] && val >= curr_ptr[c - step] && val >= curr_ptr[c - step + 1] &&
val >= curr_ptr[c + step - 1] && val >= curr_ptr[c + step] && val >= curr_ptr[c + step + 1] &&
val >= prev_ptr[c] && val >= prev_ptr[c - 1] && val >= prev_ptr[c + 1] &&
val >= prev_ptr[c - step - 1] && val >= prev_ptr[c - step] && val >= prev_ptr[c - step + 1] &&
val >= prev_ptr[c + step - 1] && val >= prev_ptr[c + step] && val >= prev_ptr[c + step + 1] &&
val >= next_ptr[c] && val >= next_ptr[c - 1] && val >= next_ptr[c + 1] &&
val >= next_ptr[c - step - 1] && val >= next_ptr[c - step] && val >= next_ptr[c - step + 1] &&
val >= next_ptr[c + step - 1] && val >= next_ptr[c + step] && val >= next_ptr[c + step + 1]) ||
(val < 0 && val <= curr_ptr[c - 1] && val <= curr_ptr[c + 1] &&
val <= curr_ptr[c - step - 1] && val <= curr_ptr[c - step] && val <= curr_ptr[c - step + 1] &&
val <= curr_ptr[c + step - 1] && val <= curr_ptr[c + step] && val <= curr_ptr[c + step + 1] &&
val <= prev_ptr[c] && val <= prev_ptr[c - 1] && val <= prev_ptr[c + 1] &&
val <= prev_ptr[c - step - 1] && val <= prev_ptr[c - step] && val <= prev_ptr[c - step + 1] &&
val <= prev_ptr[c + step - 1] && val <= prev_ptr[c + step] && val <= prev_ptr[c + step + 1] &&
val <= next_ptr[c] && val <= next_ptr[c - 1] && val <= next_ptr[c + 1] &&
val <= next_ptr[c - step - 1] && val <= next_ptr[c - step] && val <= next_ptr[c - step + 1] &&
val <= next_ptr[c + step - 1] && val <= next_ptr[c + step] && val <= next_ptr[c + step + 1])))
//++numKeys;
//获得特征点初始行号，列号，组号，组内层号
int octave = i, layer = j, r1 = r, c1 = c, nums = i * nOctaves + j;
if (!adjust_local_extrema_2(dog_pyr, kpt, octave, layer, r1, c1,
nOctaveLayers, (float)contrastThreshold,
(float)edgeThreshold, (float)sigma))
continue; //如果该初始点不满足条件，则不保存改点
float scale = kpt.size / float(1 << octave); //该特征点相对于本组的尺度
//计算梯度幅度和梯度方向
//amplit_orient(curr_img, amplit, orient, scale, nums);
//max_hist值对应的方向为主方向
float max_hist = clac_orientation_hist(gauss_pyr[octave][layer], Point(c1, r1), scale, n, hist);
//float max_hist = calc_orient_hist(amplit[nums], orient[nums], Point2f(c1, r1), scale, hist, n);
大于mag_thr值对应的方向为辅助方向
//float mag_thr = max_hist * ORI_PEAK_RATIO; //主峰值 80% 的方向作为辅助方向
//增加更多的主方向，以增加特征点对梯度差异的鲁棒性
float sum = 0.0; //直方图对应的幅值之和
float mag_thr = 0.0; //判断是否为主方向的阈值
for (int i = 0; i < n; ++i)
sum += hist[i];
mag_thr = 0.5 * (1.0 / 36) * sum;
//遍历直方图中的 36 个bin
for (int i = 0; i < n; ++i)
int left = i > 0 ? i - 1 : n - 1;
int rig