sift代码实现详解

1、创建高斯金字塔第-1组

    1.1、将源图片转成灰度图

void ConvertToGray(const Mat& src, Mat& dst)
{
    cv::Size size = src.size();
    if (dst.empty())
        dst.create(size, CV_64F);                              //[1]利用Mat类的成员函数创建Mat容器

    uchar*    srcData = src.data;                             //[2]指向存储所有像素值的矩阵的数据区域
    pixel_t*  dstData = (pixel_t*)dst.data;                   //[3]指向dst的矩阵体数据区域

    int       dstStep = dst.step / sizeof(dstData[0]);         //[4]一行图像含有多少字节数

    for (int j = 0; j<src.cols; j++)
    {
        for (int i = 0; i<src.rows; i++)
        {
            double b = (srcData + src.step*i + src.channels()*j)[0] / 255.0;
            double g = (srcData + src.step*i + src.channels()*j)[1] / 255.0;
            double r = (srcData + src.step*i + src.channels()*j)[2] / 255.0;
            *(dstData + dstStep*i + dst.channels()*j) = (r + g + b) / 3.0;
        }
    }
}

    1.2、进行上采样

/*************************************************************************************************************************
*模块说明：
*        线性插值放大，将图像放大2倍
**************************************************************************************************************************/
void UpSample(const Mat &src, Mat &dst)
{
    if (src.channels() != 1)
        return;
    dst.create(src.rows * 2, src.cols * 2, src.type());    //dst的行列是原图2倍大小的行列

    pixel_t* srcData = (pixel_t*)src.data;   
    pixel_t* dstData = (pixel_t*)dst.data;

    int srcStep = src.step / sizeof(srcData[0]);  //一行图像含有多少字节数
    int dstStep = dst.step / sizeof(dstData[0]);

    int m = 0, n = 0;
    for (int j = 0; j < src.cols - 1; j++, n += 2)
    {
        m = 0;
        for (int i = 0; i < src.rows - 1; i++, m += 2)
        {
/*第1步：n=0 m=0，n=2,m=0;  第2步：n=0 m=1，n=2,m=1; 第3步：n=1 m=0，n=3,m=0;  第4步：n=1 m=1，n=3,m=1;    第5步：最后2行2列  */
/*第1步：n=0 m=2，n=2,m=2;  第2步：n=0 m=3，n=2,m=3; 第3步：n=1 m=2，n=3,m=2;  第4步：n=1 m=3，n=3,m=3;    假设5*5的矩阵      */
            double sample = *(srcData + srcStep * i + src.channels() * j);
            *(dstData + dstStep * m + dst.channels() * n) = sample;

            double rs = *(srcData + srcStep * (i)+src.channels()*j) + (*(srcData + srcStep * (i + 1) + src.channels()*j));
            *(dstData + dstStep * (m + 1) + dst.channels() * n) = rs / 2;
            double cs = *(srcData + srcStep * i + src.channels()*(j)) + (*(srcData + srcStep * i + src.channels()*(j + 1)));
            *(dstData + dstStep * m + dst.channels() * (n + 1)) = cs / 2;

            double center = (*(srcData + srcStep * (i + 1) + src.channels() * j))
                + (*(srcData + srcStep * i + src.channels() * j))
                + (*(srcData + srcStep * (i + 1) + src.channels() * (j + 1)))
                + (*(srcData + srcStep * i + src.channels() * (j + 1)));

            *(dstData + dstStep * (m + 1) + dst.channels() * (n + 1)) = center / 4;

        }

    }



    if (dst.rows < 3 || dst.cols < 3)
        return;

    //最后两行两列
    for (int k = dst.rows - 1; k >= 0; k--)
    {
        *(dstData + dstStep *(k)+dst.channels()*(dst.cols - 2)) = *(dstData + dstStep *(k)+dst.channels()*(dst.cols - 3));
        *(dstData + dstStep *(k)+dst.channels()*(dst.cols - 1)) = *(dstData + dstStep *(k)+dst.channels()*(dst.cols - 3));
    }
    for (int k = dst.cols - 1; k >= 0; k--)
    {
        *(dstData + dstStep *(dst.rows - 2) + dst.channels()*(k)) = *(dstData + dstStep *(dst.rows - 3) + dst.channels()*(k));
        *(dstData + dstStep *(dst.rows - 1) + dst.channels()*(k)) = *(dstData + dstStep *(dst.rows - 3) + dst.channels()*(k));
    }
}

    1.3、求该层σ（sigma）

     1.4、进行高斯模糊

 计算金字塔层数（根据图像长、宽中最小的那个）

这里呢？low论文中高斯金字塔负一层创建,具体的可参考：https://www.cnblogs.com/xujianqing/p/5988627.html

SIFT算法原理详解

/*************************************************************************************************************************
*模块说明：
*        OpenCv中的高斯平滑函数
**************************************************************************************************************************/
void GaussianSmooth(const Mat &src, Mat &dst, double sigma)
{
	GaussianBlur(src, dst, Size(0, 0), sigma, sigma);
}
/*************************************************************************************************************************
*模块说明：
*        模块1--创建初始灰度图像------初始图像先将原图像灰度化，再扩大一倍后，使用高斯模糊平滑
*功能说明：
*        在最开始建立高斯金字塔的时候，要预先模糊输入的图像来作为第0个组的第0层图像，这时相当于丢弃了最高的空域的采样率。因
*        此，通常的做法是先将图像的尺度扩大一倍来生成第-1组。我们假定初始的输入图像为了抗击混淆现象，已经对其进行了sigma=0.5
*        的高斯模糊，如果输入图像的尺寸用双线性插值扩大一倍，那么相当于sigma=1.0
*这样做的原因:
*        在检测极值点前，对原始图像的高斯平滑以致图像丢失高频信息，所以Lowe建议在建立尺度空间前，首先对原始图像长宽扩展一倍，
*        以便保留原始图像的信息(特别是图像的高频信息，比如边缘，角点)，增加特征点的数量。
*函数功能:
*        这个函数的作用是用于创建高斯金字塔的-1层图像，在主函数中定义 sigma =1.6，INIT_SIGMA=O.5;
在lowe的论文中，他定义图片的尺度是1.6，被摄像头模糊的尺度是0.5，但是这个1.52对我们来说并没有啥意思，
lowe为了获得更多的特征点（姑且这么解释吧）他先对原始图像进行扩大一倍。然后呢再对它进行高斯平滑作为第一阶第一层，
高斯模糊尺度为=(1.6* 1.6 - (0.5 * 2) * (0.5 * 2))
**************************************************************************************************************************/
void CreateInitSmoothGray(const Mat &src, Mat &dst, double sigma = SIGMA)
{
	cv::Mat gray;                                  //[1]保存原始图像的灰度图像        
	cv::Mat up;                                    //[2]保存原始图像的上采样图像

	//std::cout << "CreateInitSmoothGray sigma     = " << sigma << std::endl;  //查看sigma初始值

	ConvertToGray(src, gray);                     //灰度图像
	UpSample(gray, up);                            //[3]图像的上采样
												   //[4]高斯金字塔的-1组的sigma的计算公式
	double  sigma_init = sqrt(sigma * sigma - (INIT_SIGMA * 2) * (INIT_SIGMA * 2));// INIT_SIGM初始sigma=0.5;  -1层的sigma

	GaussianSmooth(up, dst, sigma_init);           //[5]高斯平滑
}　

 2、创建高斯金字塔

    2.1、高斯模糊卷积

    2.2、降采样

隔点采用：

/*************************************************************************************************************************
*模块说明：
*        隔点采样
**************************************************************************************************************************/
void DownSample(const Mat& src, Mat& dst)
{
    if (src.channels() != 1)
        return;

    if (src.cols <= 1 || src.rows <= 1)
    {
        src.copyTo(dst);
        return;
    }

    dst.create((int)(src.rows / 2), (int)(src.cols / 2), src.type());   //dst的行列=src/2

    pixel_t* srcData = (pixel_t*)src.data;
    pixel_t* dstData = (pixel_t*)dst.data;

    int srcStep = src.step / sizeof(srcData[0]);       //src一行图像含有多少字节数
    int dstStep = dst.step / sizeof(dstData[0]);

    int m = 0, n = 0;
    for (int j = 0; j < src.cols; j += 2, n++)
    {
        m = 0;
        for (int i = 0; i < src.rows; i += 2, m++)
        {
            pixel_t sample = *(srcData + srcStep * i + src.channels() * j);
            if (m < dst.rows && n < dst.cols)
            {
                *(dstData + dstStep * m + dst.channels() * n) = sample;    //将src隔点采样到dst中去。刚好src行列/2
            }
        }
    }

}
/************************************

上述2步原理参考我的博文：SIFT解析（一）建立高斯金字塔

;

 

也就是说o代表哪一组，S表示哪一层，s表示那一层中第几幅图片，也就是第几层。

S=3,第0（o=0）组的（S+3）6幅图像（及第0组内有6层图像）分别表示为（s第几层的图片）

又

大家如果觉得还是不清楚，可以读一下百度文库中这篇文章中的高斯金字塔构建原理，讲解的很清楚：

 https://wenku.baidu.com/view/d7edd2464b73f242336c5ffa.html

这部分代码可以提现为

/*************************************************************************************************************************
*模块说明：
*        OpenCv中的高斯平滑函数
**************************************************************************************************************************/
void GaussianSmooth(const Mat &src, Mat &dst, double sigma)
{
    GaussianBlur(src, dst, Size(0, 0), sigma, sigma);
}
/*************************************************************************************************************************
*模块说明：
*        模块1--创建初始灰度图像------初始图像先将原图像灰度化，再扩大一倍后，使用高斯模糊平滑
*功能说明：
*        在最开始建立高斯金字塔的时候，要预先模糊输入的图像来作为第0个组的第0层图像，这时相当于丢弃了最高的空域的采样率。因
*        此，通常的做法是先将图像的尺度扩大一倍来生成第-1组。我们假定初始的输入图像为了抗击混淆现象，已经对其进行了sigma=0.5
*        的高斯模糊，如果输入图像的尺寸用双线性插值扩大一倍，那么相当于sigma=1.0
*这样做的原因:
*        在检测极值点前，对原始图像的高斯平滑以致图像丢失高频信息，所以Lowe建议在建立尺度空间前，首先对原始图像长宽扩展一倍，
*        以便保留原始图像的信息(特别是图像的高频信息，比如边缘，角点)，增加特征点的数量。
*函数功能:
*        这个函数的作用是用于创建高斯金字塔的-1层图像，在主函数中定义 sigma =1.6
在lowe的论文中，他定义图片的尺度是1.6，被摄像头模糊的尺度是0.5，但是这个1.52对我们来说并没有啥意思，
lowe为了获得更多的特征点（姑且这么解释吧）他先对原始图像进行扩大一倍。然后呢再对它进行高斯平滑作为第一阶第一层，
高斯模糊尺度为=(1.6* 1.6 - (0.5 * 2) * (0.5 * 2))
**************************************************************************************************************************/
void CreateInitSmoothGray(const Mat &src, Mat &dst, double sigma = SIGMA)
{
    cv::Mat gray;                                  //[1]保存原始图像的灰度图像        
    cv::Mat up;                                    //[2]保存原始图像的上采样图像

    //std::cout << "CreateInitSmoothGray sigma     = " << sigma << std::endl;  //查看sigma初始值

    ConvertToGray(src, gray);                     //灰度图像
    UpSample(gray, up);                            //[3]图像的上采样
                                                   //[4]高斯金字塔的-1组的sigma的计算公式
    double  sigma_init = sqrt(sigma * sigma - (INIT_SIGMA * 2) * (INIT_SIGMA * 2));// INIT_SIGM初始sigma=0.5;  -1层的sigma

    GaussianSmooth(up, dst, sigma_init);           //[5]高斯平滑
}
/*************************************************************************************************************************
*模块说明：
*        模块二：3.3 图像高斯金字塔的构建
在主函数中：int octaves = log((double)min(init_gray.rows, init_gray.cols)) / log(2.0) - 2;             //计算高斯金字塔的层数
高斯金字塔中每组图像中有三层/张图片 INTERVALS=3；SIGMA=1.6；
**************************************************************************************************************************/
void GaussianPyramid(const Mat &src, vector<Mat>&gauss_pyr, int octaves, int intervals = INTERVALS, double sigma = SIGMA)
{
    //std::cout << "sigma is     = " << sigma << std::endl;
    double *sigmas = new double[intervals + 3];
    double k = pow(2.0, 1.0 / intervals);      //pow（x,y）,x的y次方；k=2的（1/s)次方
//    std::cout << "k is     = " << k << std::endl;
    sigmas[0] = sigma;
//    std::cout << "sigma0 is     = " << sigmas[0] << std::endl;
    double sig_prev;
    double sig_total;
//用循环来表示第0组内各层图片的尺度,这个可以表示为0组的各层图像尺度；
    for (int i = 1; i < intervals + 3; i++)
    {
        sig_prev = pow(k, i - 1) * sigma;     //第i-1层尺度
        sig_total = sig_prev * k;             //第i层尺度
        sigmas[i] = sqrt(sig_total * sig_total - sig_prev * sig_prev);  //组内每层尺度坐标
    }

    for (int o = 0; o < octaves; o++)
    {
        //每组多三层
        for (int i = 0; i < intervals + 3; i++)
        {
            Mat mat;
            if (o == 0 && i == 0)
            {
                src.copyTo(mat);   //第0组，第0层，就是原图像
            }
            else if (i == 0)     //新的一组的首层图像是由上一组最后一层图像下采样得到
            {
                DownSample(gauss_pyr[(o - 1)*(intervals + 3) + intervals], mat);
            }
            else                 //对上一层图像gauss_pyr[i-1]进行参数为sig[i]的高斯平滑，得到当前层图像

            {
                GaussianSmooth(gauss_pyr[o * (intervals + 3) + i - 1], mat, sigmas[i]);
        
            }
            gauss_pyr.push_back(mat); //将mat中的像素添加到gauss_pyr数组里面；也就是说每组从第0层到第5层的图片像素装进vector数组中；
        }
    }

    delete[] sigmas;  //释放sigmas参数数组
}

但是在原理中，我发现一个很费解的问题

而在论文代码中，并没有看出在哪体现出来每一处的2^^o次方这个概念。后来查找文章才发现？

 在极值比较的过程中，每一组图像的首末两层是无法进行极值比较的，为了满足尺度变化的连续性；

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

这里有的童鞋不理解什么叫“为了满足尺度变化的连续性”，现在做仔细阐述：

假设s=3，也就是每个塔里有3层，则k=2^1/s=2^1/3，那么按照上图可得Gauss Space和DoG space 分别有3个（s个）和2个（s-1个）分量，

在DoG space中，1st-octave两项分别是σ,kσ; 2nd-octave两项分别是2σ,2kσ;由于无法比较极值，我们必须在高斯空间继续添加高斯模糊项，

使得形成σ,kσ,k2σ,k3σ,k4σ这样就可以选择DoG space中的中间三项kσ,k2σ,k3σ（只有左右都有才能有极值），

那么下一octave中（由上一层降采样获得）所得三项即为2kσ,2k2σ,2k3σ，其首项2kσ=2^4/3。刚好与上一octave末项k3σ=2^3/3尺度变化连续起来，

所以每次要在Gaussian space添加3项，每组（塔）共S+3层图像，相应的DoG金字塔有S+2层图像。

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

这个尺度因子都是在原图上进行的。而在算法实现过程中，采用高斯平滑是在上一层图像上再叠加高斯平滑。即我们在程序中看到的平滑因子为（4-5）

 我们在源代码上同一时候也没有看到组间的2倍的关系，实际在对每组的平滑因子都是一样的，2倍的关系是因为在降採样的过程中产生的，第二层的第一张图是由第一层的平滑因子为2σ的图像上（即倒数第三张）降採样得到，此时图像平滑因子为σ，所以继续採用以上的平滑因子。但在原图上看。形成了所有的空间尺度。

 这里一定要写下参考大神的博格：

https://www.cnblogs.com/Alliswell-WP/p/SIFT_code1.html     //这个是很多问题的解析,包括各种问题；
https://blog.csdn.net/xw20084898/article/details/16832755  //这个参考了主要是尺度变换的连续性

 3、求高斯差分金字塔（用DoG）

两个高斯金字塔层相减生成一个差分金字塔。G（x, y, σ）已由上面算出

G（x, y, σ1）- G（x, y, σ2）

/*************************************************************************************************************************
*模块说明：
*        图像的差分
**************************************************************************************************************************/
void Sub(const Mat& a, const Mat& b, Mat & c)
{
    if (a.rows != b.rows || a.cols != b.cols || a.type() != b.type())
        return;
    if (!c.empty())
        return;
    c.create(a.size(), a.type());

    pixel_t* ap = (pixel_t*)a.data;
    pixel_t* ap_end = (pixel_t*)a.dataend;
    pixel_t* bp = (pixel_t*)b.data;
    pixel_t* cp = (pixel_t*)c.data;
    int step = a.step / sizeof(pixel_t);

    while (ap != ap_end)
    {
        *cp++ = *ap++ - *bp++;
    }
}
/*************************************************************************************************************************
*模块说明：
*       模块三：3.4 高斯差分金字塔
*功能说明：
*       1--2002年，Mikolajczyk在详细的实验比较中发现尺度归一化的高斯拉普拉斯函数的极大值和极小值同其他的特征提取函数，例如：
*          梯度，Hessian或者Harris角点特征比较，能够产生稳定的图像特征。
*       2--而Lindberg早在1994年发现高斯差分函数（Difference of Gaussian,简称DOG算子）与尺度归一化的拉普拉斯函数非常相似，因此，
*          利用图像间的差分代替微分。
**************************************************************************************************************************/
void DogPyramid(const vector<Mat>& gauss_pyr, vector<Mat>& dog_pyr, int octaves, int intervals = INTERVALS)
{
    for (int o = 0; o < octaves; o++)
    {
        for (int i = 1; i < intervals + 3; i++)
        {
            Mat mat;
            Sub(gauss_pyr[o*(intervals + 3) + i], gauss_pyr[o*(intervals + 3) + i - 1], mat);
            dog_pyr.push_back(mat);
        }
    }
}　

 总结：

高斯金字塔共分o组(Octave), 每组又分S层(Layer)。 组内各层图像的分辨率是相同的，即长和宽相同，但尺度逐渐增加，即越往塔顶图像越模糊。而下-组的图像是由上一组图像按照隔点降采样得到的，即图像的长和宽分别减半。

 octaves = log((double)min(init_gray.rows, init_gray.cols)) / log(2.0) - 2;             //计算高斯金字塔的层数

本实验中o=6，S=s+3=6,小写的s是高斯差分金字塔的层数；

但是我们在代码中可以看到，每组中由6层，sub函数差分图像后每组还剩5层；

那为什么我们说小写的s（s=3）是高斯差分金字塔的层数；这就体现在后面的极值点检测中；

4、极值点检测

     4.1、排除阈值小的点

     4.2、判断是否是极值

    4.3、修正极值点，删除不稳定的点

上述2步原理参考我的博文：SIFT算法原理（2）-极值点的精确定位

首先，第一个问题：

问题描述：如何找寻关键点？在几层之间对比？思路是什么？

答：极值的检测是在DoG空间进行的，检测是以前点为中心，3pixel*3pixel*3pixel的立方体为邻域，判断当前点是否为局部最大或最小。如下图所示，橘黄色为当前检测点，绿色点为其邻域。因为要比较当前点的上下层图像，所以极值检测从DoG每层的第2幅图像开始，终止于每层的倒数第2幅图像（第1幅没有下层，最后1幅没有上层，无法比较）。所以实际上我们极值点只在中间3层进行检测

极值点初步检测：

/*************************************************************************************************************************
*模块说明：
*       模块四：3.5 空间极值点的检测(关键点的初步探查)
*功能说明：
*       1--关键点是由DOG空间的局部极值点组成的，关键点的初步探查是通过同一组内各DoG相邻两层图像之间的比较完成的。为了寻找DoG
*          函数的极值点，每一个像素点都要和它所有相邻的点比较，看其是否比它的图像域和尺度域相邻的点大还是小。
*       2--当然这样产生的极值点并不全都是稳定的特征点，因为某些极值点相应较弱，而且DOG算子会产生较强的边缘响应。
**************************************************************************************************************************/
void DetectionLocalExtrema(const vector<Mat>& dog_pyr, vector<Keypoint>& extrema, int octaves, int intervals = INTERVALS)
{

	double  thresh = 0.5 * DXTHRESHOLD / intervals;  //|D(x^)| < 0.03 ; DXTHRESHOLD=0.03极值点太多

	for (int o = 0; o < octaves; o++)
	{
		//第一层和最后一层极值忽略
		for (int i = 1; i < (intervals + 2) - 1; i++)             //Dog图像每一组5层（差分6-1），又第一层和最后一程忽略；
		{
			int index = o*(intervals + 2) + i;                              //[1]图片索引的定位
			pixel_t *data = (pixel_t *)dog_pyr[index].data;                //[2]获取图片的矩阵体的首地址
		//	std::cout << "首地址为      = " << sizeof(data[0]) << std::endl;

			int step = dog_pyr[index].step / sizeof(data[0]);           //[3]说明矩阵在存储空间中的存储是以线性空间的方式存放的


			for (int y = IMG_BORDER; y < dog_pyr[index].rows - IMG_BORDER; y++)   //IMG_BORDER=5 ，这里我查找资料说的是指检测边界以内的点，因为边界易受干扰
			{
				for (int x = IMG_BORDER; x < dog_pyr[index].cols - IMG_BORDER; x++)
				{
				//	std::cout << "x is     = " << x << std::endl;
					pixel_t val = *(data + y*step + x);                     //遍历像素点
					if (std::fabs(val) > thresh)                           //[4]排除阈值过小的点，fabs（）取绝对值
					{
						if (isExtremum(x, y, dog_pyr, index))                //[5]判断是否是极值
						{
							Keypoint *extrmum = InterploationExtremum(x, y, dog_pyr, index, o, i);//极值点x,y；当前组当前层；
							if (extrmum)
							{
								if (passEdgeResponse(extrmum->x, extrmum->y, dog_pyr, index))
								{
									extrmum->val = *(data + extrmum->y*step + extrmum->x);
									extrema.push_back(*extrmum);
								}

								delete extrmum;

							}//extrmum
						}//isExtremum
					}//std::fabs
				}//for x
			}//for y

		}
	}
}

isExtremum这个函数主要是判断极值点，主要原理参考图2；代码如下：

/*************************************************************************************************************************
*模块说明：
*       模块四的第一步:3.4-空间极值点的初步检测(关键点的初步探查)
*功能说明:
*       1--在高斯差分函数之后，得到一系列的关键点的疑似点，我们需要对这些关键点的疑似点初步进行检测和筛选
*      检测是以前点为中心，3pixel*3pixel*3pixel的立方体为邻域，判断当前点是否为局部最大或最小。i,j,k=-1,0,1表示当前点的邻域
**************************************************************************************************************************/
bool isExtremum(int x, int y, const vector<Mat>& dog_pyr, int index)
{
	pixel_t * data = (pixel_t *)dog_pyr[index].data;
	int      step = dog_pyr[index].step / sizeof(data[0]);
	pixel_t   val = *(data + y*step + x);

	if (val > 0)
	{
		for (int i = -1; i <= 1; i++)    //层
		{
			int stp = dog_pyr[index + i].step / sizeof(data[0]);
			for (int j = -1; j <= 1; j++)     //行
			{
				for (int k = -1; k <= 1; k++)  //列
				{
					if (val < *((pixel_t*)dog_pyr[index + i].data + stp*(y + j) + (x + k)))
					{
						return false;   //不是极值点退出
					}
				}
			}
		}
	}
	else
	{
		for (int i = -1; i <= 1; i++)
		{
			int stp = dog_pyr[index + i].step / sizeof(pixel_t);
			for (int j = -1; j <= 1; j++)
			{
				for (int k = -1; k <= 1; k++)
				{
					//检查最小极值
					if (val > *((pixel_t*)dog_pyr[index + i].data + stp*(y + j) + (x + k)))
					{
						return false;  //不是极值点退出
					}
				}
			}
		}
	}
	return true;  //是级值点，返回
}

一个定义不好的高斯差分算子的极值在横跨边缘的地方有较大的主曲率，而在垂直边缘的方向有较小的主曲率。DOG算子会产生较强的边缘响应，需要剔除不稳定的边缘响应点。获取特征点处的Hessian矩阵，主曲率通过一个2x2 的Hessian矩阵H求出（D的主曲率和H的特征值成正比）；

原理参考：SIFT算法原理（2）-极值点的精确定位

代码如下：

/*************************************************************************************************************************
*模块说明：
*       模块四的第三步:4.2--消除边缘响应点
*功能说明:
*       1）一个定义不好的高斯差分算子的极值在横跨边缘的地方有较大的住主曲率，在垂直边缘的方向有较小的主曲率。
*       2）DOG算子会产生较强的边缘响应，需要剔除不稳定的边缘响应点，获取特征点处的Hessian矩阵，主曲率通过一个2*2的Hessian矩
*          阵H求出
*       3）主曲率D和Hessian矩阵的特征值成正比，公式(r+1)*(r+1)/r的值在两个特征值相等时最小；这个值越大，说明两个特征值的比值
*          越大，即在某一个方向的梯度值越大，而在另一个方向的梯度值越小，而边缘恰恰就是这种情况。所以，为了剔除边缘响应点，
*          需要让该比值小于一定的阈值，因此，为了检测主曲率是否在某阈值r下，只需检测。CSDN论文中的公式(4-7成立)，成立，将关
*          键点保留，反之，剔除。
**************************************************************************************************************************/
#define DAt(x, y) (*(data+(y)*step+(x))) 
bool passEdgeResponse(int x, int y, const vector<Mat>& dog_pyr, int index, double r = RATIO)
{
	pixel_t *data = (pixel_t *)dog_pyr[index].data;  //极值点的图片索引定位
	int step = dog_pyr[index].step / sizeof(data[0]);
	pixel_t val = *(data + y*step + x);    //极值点

	double Dxx;
	double Dyy;
	double Dxy;
	double Tr_h;                                                         //[1]Hessian矩阵的迹
	double Det_h;                                                        //[2]Hessian矩阵所对应的行列式的值

	Dxx = DAt(x + 1, y) + DAt(x - 1, y) - 2 * val;
	Dyy = DAt(x, y + 1) + DAt(x, y - 1) - 2 * val;
	Dxy = (DAt(x + 1, y + 1) + DAt(x - 1, y - 1) - DAt(x - 1, y + 1) - DAt(x + 1, y - 1)) / 4.0;

	Tr_h = Dxx + Dyy;
	Det_h = Dxx * Dyy - Dxy * Dxy;

	if (Det_h <= 0)return false;

	if (Tr_h * Tr_h / Det_h < (r + 1) * (r + 1) / r) return true;

	return false;
}

5、计算极值点尺度，图像特征缩放

/*************************************************************************************************************************
*模块说明：
*       模块五：计算特征点对应的尺度

*
**************************************************************************************************************************/
void CalculateScale(vector<Keypoint>& features, double sigma = SIGMA, int intervals = INTERVALS)
{
	double intvl = 0;
	for (int i = 0; i < features.size(); i++)
	{
		intvl = features[i].interval + features[i].offset_interval;
		features[i].scale = sigma * pow(2.0, features[i].octave + intvl / intervals);
		features[i].octave_scale = sigma * pow(2.0, intvl / intervals);
	}

}

//对扩大的图像特征缩放
/*将特征点序列中每个特征点的坐标减半(图像金字塔设置了将图像放大为原图的2倍时，特征点检测完之后调用)
/features：特征点序列*/
void HalfFeatures(vector<Keypoint>& features)
{
	for (int i = 0; i < features.size(); i++)
	{
		features[i].dx /= 2;
		features[i].dy /= 2;

		features[i].scale /= 2;
	}
}　

关键点的x代表cols,列；关键点的y代表rows,行

features[i].scale = sigma * pow(2.0, features[i].octave + intvl / intervals);//尺度

   features[i].octave_scale = sigma * pow(2.0, intvl / intervals);

6、关键点方向分配

   6.1、计算梯度直方图

   6.2、对直方图做两次高斯平滑

   6.3、求直方图中的主方向

   6.4、使求出的主方向更加精确 

原理可参考：RobHess的SIFT代码解析步骤三

主代码：

/********************************************************************************************************************************
*模块说明:
*        模块六：5 关键点方向分配
*功能说明：
*        1）为了使描述符具有旋转不变性，需要利用图像的局部特征为每一个关键点分配一个基准方向。使用图像梯度的方法求取局部结构的稳定
*           方向。
*        2）对于在DOG金字塔中检测出来的关键点，采集其所在高斯金字塔图像3sigma邻域窗口内像素的梯度和方向梯度和方向特征。
*        3）梯度的模和方向如下所示:
*        4) 在完成关键点的梯度计算后，使用直方图统计邻域内像素的梯度和方向。梯度直方图将0~360度的方向范围分为36个柱，其中每柱10度，
*           如图5.1所示，直方图的峰值方向代表了关键点的主方向
         cvRound()：返回跟参数最接近的整数值，即四舍五入；
*********************************************************************************************************************************/
void OrientationAssignment(vector<Keypoint>& extrema, vector<Keypoint>& features, const vector<Mat>& gauss_pyr)
{
	int n = extrema.size();
	double *hist;

	for (int i = 0; i < n; i++)
	{
		//ORI_HIST_BINS=36;ORI_WINDOW_RADIUS=3.0*1.5,特征点方向赋值过程中，搜索邻域的半径为：3*1.5*σ
		//ORI_SIGMA_TIMES=1.5,
		hist = CalculateOrientationHistogram(gauss_pyr[extrema[i].octave*(INTERVALS + 3) + extrema[i].interval],
			extrema[i].x, extrema[i].y, ORI_HIST_BINS, cvRound(ORI_WINDOW_RADIUS*extrema[i].octave_scale),
			ORI_SIGMA_TIMES*extrema[i].octave_scale);                             //[1]计算梯度的方向直方图

		for (int j = 0; j < ORI_SMOOTH_TIMES; j++)
			GaussSmoothOriHist(hist, ORI_HIST_BINS);                              //[2]对方向直方图进行高斯平滑
		double highest_peak = DominantDirection(hist, ORI_HIST_BINS);            //[3]求取方向直方图中的峰值
																				 //[4]计算更加精确的关键点主方向
		CalcOriFeatures(extrema[i], features, hist, ORI_HIST_BINS, highest_peak*ORI_PEAK_RATIO); //ORI_PEAK_RATIO=0.8
		//精确关键点的主方向，抛物插值
		delete[] hist;

	}
}

调用代码相关代码如下：

问题描述：为什么要计算特征梯度直方图？怎么计算？

答：在完成特征点邻域的高斯图像的梯度计算后，使用直方图统计邻域内像素的梯度方向和幅值。梯度方向直方图的横轴是梯度方向角，纵轴是梯度方向角对应的（带高斯权重）梯度幅值累加值。梯度方向直方图将。0°~360°的范围，分为36个柱，每10°为一个柱。直方图的峰值代表了该特征点处邻域内图像梯度的主方向，也即该特征点的主方向，

 关键点的直方图计算

/********************************************************************************************************************************
*模块说明:
*        模块六---步骤1：计算梯度的方向直方图
*功能说明：
*        1）直方图以每10度为一个柱，共36个柱，柱代表的方向为为像素点的梯度方向，柱的长短代表了梯度幅值。
*        2）根据Lowe的建议，直方图统计采用3*1.5*sigma
*        3）在直方图统计时，每相邻三个像素点采用高斯加权，根据Lowe的建议，模板采用[0.25,0.5,0.25],并且连续加权两次
*结    论：
*        图像的关键点检测完毕后，每个关键点就拥有三个信息：位置、尺度、方向；同时也就使关键点具备平移、缩放和旋转不变性
bin=36.
*********************************************************************************************************************************/
double* CalculateOrientationHistogram(const Mat& gauss, int x, int y, int bins, int radius, double sigma)
{
	//std::cout << "首地址1为      = " << radius << std::endl;  //radius=3*1.5?*?σ(extrema[i].octave_scale)
    //std::cout << "首地址2为      = " << bins << std::endl;   // sigma=1.5*σ(extrema[i].octave_scale)
	double* hist = new double[bins];                           //[1]动态分配一个double类型的数组，用来存放梯度直方图
	for (int i = 0; i < bins; i++)                               //[2]给这个数组初始化
		*(hist + i) = 0.0;

	double  mag;                                                //[3]关键点的梯度幅值                                          
	double  ori;                                                //[4]关键点的梯度方向
	double  weight;

	int           bin;
	const double PI2 = 2.0*CV_PI;  //2Π
	double       econs = -1.0 / (2.0*sigma*sigma);

	for (int i = -radius; i <= radius; i++)
	{
		for (int j = -radius; j <= radius; j++)
		{
			if (CalcGradMagOri(gauss, x + i, y + j, mag, ori))       //[5]计算该关键点的梯度幅值和方向
			{
			//	std::cout << "ori为      = " << ori << std::endl;
				weight = exp((i*i + j*j)*econs);              // 该点的梯度幅值权重
				bin = cvRound(bins * (CV_PI - ori) / PI2);     //[6]对一个double行的数进行四舍五入，返回一个整形的数;
			//	bin = cvRound(bins * (CV_PI + ori) / PI2);     //[6]对一个double行的数进行四舍五入，返回一个整形的数;
				//计算梯度的方向对应的直方图中的bin下标,ori范围(-Π，Π]，+Π后范围(0，2Π]
			//	std::cout << "bin1为      = " << bin << std::endl;
				bin = bin < bins ? bin : 0;            //把bin的下标归到0~36；bin<bins，bin=bin;否则bin=0；
				//std::cout << "bin2为      = " << bin << std::endl;

				hist[bin] += mag * weight;                      //[7]统计梯度的方向直方图
			}
		}
	}

	return hist;
}

 这个公式主要实在计算特征值的方向和幅值，体现在CalcGradMagOri（）函数上，代码如下：

/********************************************************************************************************************************
*模块说明:
*        模块六---步骤2：计算关键点的梯度和梯度方向
*功能说明：
*        1）计算关键点(x,y)处的梯度幅值和梯度方向
*        2）将所计算出来的梯度幅值和梯度方向保存在变量mag和ori中
*********************************************************************************************************************************/
bool CalcGradMagOri(const Mat& gauss, int x, int y, double& mag, double& ori)
{
	if (x > 0 && x < gauss.cols - 1 && y > 0 && y < gauss.rows - 1)
	{
		pixel_t *data = (pixel_t*)gauss.data;
		int step = gauss.step / sizeof(*data);

		double dx = *(data + step*y + (x + 1)) - (*(data + step*y + (x - 1)));           //[1]利用X方向上的差分代替微分dx
		double dy = *(data + step*(y + 1) + x) - (*(data + step*(y - 1) + x));           //[2]利用Y方向上的差分代替微分dy

		mag = sqrt(dx*dx + dy*dy);                                          //[3]计算该关键点的梯度幅值
		ori = atan2(dy, dx);                                                //[4]计算该关键点的梯度方向
		return true;
	}
	else
		return false;
}

atan2相关知识汇总——https://blog.csdn.net/weixin_42142612/article/details/80972768

对直方图做两次高斯平滑

LOWE建议对直方图进行平滑，降低突变的影响，弥补因没有仿射不变性而产生的特征点不稳定的问题

hist[0]=0.25*hist[35]+0.5*hist[0]+0.25*hist[1];

/********************************************************************************************************************************
*模块说明:
*        模块六---步骤3：对梯度方向直方图进行连续两次的高斯平滑
*功能说明：
*        1）在直方图统计时，每相邻三个像素点采用高斯加权，根据Lowe的建议，模板采用[0.25,0.5,0.25],并且连续加权两次
*        2）对直方图进行两次平滑
*********************************************************************************************************************************/
void GaussSmoothOriHist(double *hist, int n)
{
	double prev = hist[n - 1];
	double temp;
	double h0 = hist[0];

	for (int i = 0; i < n; i++)
	{
		temp = hist[i];
		hist[i] = 0.25 * prev + 0.5 * hist[i] + 0.25 * (i + 1 >= n ? h0 : hist[i + 1]);//对方向直方图进行高斯平滑
		prev = temp;
	}
}

   6.3、求直方图中的主方向

   6.4、使求出的主方向更加精确 

方向直方图的峰值则代表了该特征点处邻域梯度的方向，以直方图中最大值作为该关键点的主方向。为了增强匹配的鲁棒性，只保留峰值大于主方向峰值80％的方向（另一个方向的直方图的幅值大于主方向峰值80％的方向）作为该关键点的辅方向。

 

/********************************************************************************************************************************
*模块说明:
*        模块六---步骤5：计算更加精确的关键点主方向----抛物插值
*功能说明：
*        1）方向直方图的峰值则代表了该特征点的方向，以直方图中的最大值作为该关键点的主方向。为了增强匹配的鲁棒性，只保留峰值大于主
*           方向峰值80%的方向作为改关键点的辅方向。因此，对于同一梯度值得多个峰值的关键点位置，在相同位置和尺度将会有多个关键点被
*           创建但方向不同。仅有15%的关键点被赋予多个方向，但是可以明显的提高关键点的稳定性。
*        2）在实际编程中，就是把该关键点复制成多份关键点，并将方向值分别赋给这些复制后的关键点
*        3）并且，离散的梯度直方图要进行【插值拟合处理】，来求得更加精确的方向角度值
*********************************************************************************************************************************/
#define Parabola_Interpolate(l, c, r) (0.5*((l)-(r))/((l)-2.0*(c)+(r))) //求小柱子在直方图中的索引号
void CalcOriFeatures(const Keypoint& keypoint, vector<Keypoint>& features, const double *hist, int n, double mag_thr)
{
	double  bin;
	double  PI2 = CV_PI * 2.0;
	int l;
	int r;
	//遍历直方图
	for (int i = 0; i < n; i++)
	{
		l = (i == 0) ? n - 1 : i - 1; //前一个(左边的)bin的下标，
		r = (i + 1) % n;              //后一个(右边的)bin的下标
		//i=0,l=35,r=1;

		//hist[i]是极值
		//若当前的bin是局部极值(比前一个和后一个bin都大)，并且值大于给定的幅值阈值，则新生成一个特征点并添加到特征点序列末尾
		if (hist[i] > hist[l] && hist[i] > hist[r] && hist[i] >= mag_thr)  //mag_thr=0.8*关键点最大值
		{
			//根据左、中、右三个bin的值对当前bin进行直方图插值
			bin = i + Parabola_Interpolate(hist[l], hist[i], hist[r]);
			bin = (bin < 0) ? (bin + n) : (bin >= n ? (bin - n) : bin); //将插值结果规范到[0,n]内

			Keypoint new_key;

			CopyKeypoint(keypoint, new_key); //复制特征点

			new_key.ori = ((PI2 * bin) / n) - CV_PI; //方向还原为（-Π，Π]，bin最大值36，2Π*36/36-Π=Π，bin最小值0，2Π*0/36-Π=-Π
			features.push_back(new_key); //插入到特征点序列末尾
		}
	}
}

//复制关键点
void CopyKeypoint(const Keypoint& src, Keypoint& dst)
{
	dst.dx = src.dx;
	dst.dy = src.dy;

	dst.interval = src.interval;
	dst.octave = src.octave;
	dst.octave_scale = src.octave_scale;
	dst.offset_interval = src.offset_interval;

	dst.offset_x = src.offset_x;
	dst.offset_y = src.offset_y;

	dst.ori = src.ori;
	dst.scale = src.scale;
	dst.val = src.val;
	dst.x = src.x;
	dst.y = src.y;
}

再次讨论一下features

关键点的x代表cols,列；关键点的y代表rows,行


features[i].scale = sigma * pow(2.0, features[i].octave + intvl / intervals);//尺度


features[i].octave_scale = sigma * pow(2.0, intvl / intervals);

梯度幅值和梯度方向保存在变量mag和ori中

7、关键点描述

   7.1、确定描述子所需的邻域区域

问题1：所需的图像区域半径radius怎么计算？

答：将关键点附近的区域划分为d*d(Lowe建议d=4)个子区域，每个子区域作为一个种子点，每个种子点有8个方向。考虑到实际计算时，需要采用三线性插值，所需图像窗口边长为3x3xσ_oct x(d+1)  。在考虑到旋转因素(方便下一步将坐标轴旋转到关键点的方向)，如下图所示，实际计算所需的图像区域半径为：

部分代码：

//为了保证特征描述子具有旋转不变性，要以特征点为中心，在附近邻域内旋转θ角，即旋转为特征点的方向
	double cos_ori = cos(ori);
	double sin_ori = sin(ori);

	//6.1高斯权值，sigma等于描述字窗口宽度的一半
	double sigma = 0.5 * width;
	double conste = -1.0 / (2 * sigma*sigma);

	double PI2 = CV_PI * 2;

	//计算特征描述子过程中，特征点周围的width*width个区域中，每个区域的宽度为3*σ个像素，
	double sub_hist_width = DESCR_SCALE_ADJUST * octave_scale; //DESCR_SCALE_ADJUST=3

	//【1】计算描述子所需的图像领域区域的半径
	//考虑到要进行双线性插值，每个区域的宽度应为:DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )
	//在考虑到旋转因素，每个区域的宽度应为：DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )* sqrt(2)
	//所以搜索的半径是：DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )* sqrt(2) / 2
	int    radius = (sub_hist_width*sqrt(2.0)*(width + 1)) / 2.0 + 0.5;    //[1]0.5取四舍五入

7.2、（依据6求出的方向）旋转坐标轴，使具有旋转不变性

将坐标轴旋转为关键点的方向，以确保旋转不变性。

/遍历每个区域的像素
	for (int i = -radius; i <= radius; i++)
	{
		for (int j = -radius; j <= radius; j++)
		{
			//坐标旋转为主方向,以确保旋转不变性。
			double rot_x = (cos_ori * j - sin_ori * i) / sub_hist_width;
			double rot_y = (sin_ori * j + cos_ori * i) / sub_hist_width;

			//将邻域内的采样点分配到对应的子区域内，将子区域内的梯度值分配到8个方向上，计算其权值
			double xbin = rot_x + width / 2 - 0.5;                         //[2]xbin,ybin为落在4*4窗口中的下标值
			double ybin = rot_y + width / 2 - 0.5;

7.3、将邻域内采样点分配到旋转后的对应子区域，确认种子点

if (xbin > -1.0 && xbin < width && ybin > -1.0 && ybin < width)
			{
				if (CalcGradMagOri(gauss, x + j, y + i, grad_mag, grad_ori)) //[3]计算关键点的梯度方向
				{
					grad_ori = (CV_PI - grad_ori) - ori;
					while (grad_ori < 0.0)
						grad_ori += PI2;
					while (grad_ori >= PI2)
						grad_ori -= PI2;              //把梯度方向归入[0，2Π]

					double obin = grad_ori * (bins / PI2);   
					//梯度方向*(8/2Π)归到[0，8]，8为直方图有8个bin，2Π为角度取值范围的长度，把方向的索引归于0~8.
					//conste = -1.0 / (2 * sigma*sigma);
					//计算权重
					double weight = exp(conste*(rot_x*rot_x + rot_y * rot_y));　

 7.4、计算种子点的8个方向的梯度信息

三线性差值：

/*
功能说明：
将一个条目插入到形成特征描述符的方向直方图数组中，该条目分配最多8个方向，对于每个维度每个方向乘以（1-d） 的权重，
其中d是以bin测量的bin中心值距离。
//hist为方向直方图的三维数组
//xbin子行
//ybin子列
//obin子方向
//mag权重
//d方向直方图的宽度
//bins每个直方图中bin的个数
三线性插值
*/
void InterpHistEntry(double ***hist, double xbin, double ybin, double obin, double mag, int bins, int d)
{
	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(ybin);  //向下取整
	c0 = cvFloor(xbin);
	o0 = cvFloor(obin);
	d_r = ybin - r0;   //小数余项
	d_c = xbin - c0;
	d_o = obin - o0;

//这里也可以看出前面rbin，cbin为何减0.5，这样原点上这点d_r,d_c均为0.5，即原点上这点的方向将被平均分配叠加在它周围4个直方图上面。

	/*
	做插值：
	xbin,ybin,obin:种子点所在子窗口的位置和方向
	所有种子点都将落在4*4的窗口中
	r0,c0取不大于xbin，ybin的正整数
	r0,c0只能取到0,1,2
	xbin,ybin在(-1, 2)
	r0取不大于xbin的正整数时。
	r0+0 <= xbin <= r0+1
	mag在区间[r0,r1]上做插值
	obin同理
	*/

	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);  //对近邻两行的贡献因子
			row = hist[rb];    //第rb行上的hist
			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];   //h表示第rb行cb列上的hist
					for (o = 0; o <= 1; o++)   //沿着直方图方向不超过1个单位长度
					{
						ob = (o0 + o) % bins;        //bins=8，8个柱子
						v_o = v_c * ((o == 0) ? 1.0 - d_o : d_o);  //对近邻两个方向的贡献因子
						h[ob] += v_o;   
					}
				}
			}
		}
	}


}

这里三线性差值还不是很懂，后续再继续学习

前面几个步骤的总代码：

/********************************************************************************************************************************
*模块说明:
*        模块七--步骤1:计算描述子的直方图
*功能说明：
gauss：图像指针

x：特征点所在的行

y：特征点所在的列

ori：特征点的主方向

octave_scale：特征点的尺度//特征点所在组的尺度

width：计算方向直方图时，将特征点附近划分为width*width个区域，每个区域生成一个直方图，默认width为4

bins：每个直方图中bins的个数

返回值：double类型的三维数组，即一个d*d的二维数组，数组中每个元素是一个有n个bins的直方图数组
*********************************************************************************************************************************/
double*** CalculateDescrHist(const Mat& gauss, int x, int y, double octave_scale, double ori, int bins, int width)
{
	double ***hist = new double**[width];   //width*width*bins的三维直方图数组 //为第一维分配空间

	for (int i = 0; i < width; i++)
	{
		hist[i] = new double*[width];  //为第二维分配空间
		for (int j = 0; j < width; j++)
		{
			hist[i][j] = new double[bins]; //为第三维分配空间
		}
	}
	
	for (int r = 0; r < width; r++)
		for (int c = 0; c < width; c++)
			for (int o = 0; o < bins; o++)
				hist[r][c][o] = 0.0;

	//为了保证特征描述子具有旋转不变性，要以特征点为中心，在附近邻域内旋转θ角，即旋转为特征点的方向
	double cos_ori = cos(ori);
	double sin_ori = sin(ori);

	//6.1高斯权值，sigma等于描述字窗口宽度的一半
	double sigma = 0.5 * width;
	double conste = -1.0 / (2 * sigma*sigma);

	double PI2 = CV_PI * 2;

	//计算特征描述子过程中，特征点周围的width*width个区域中，每个区域的宽度为3*σ个像素，
	double sub_hist_width = DESCR_SCALE_ADJUST * octave_scale; //DESCR_SCALE_ADJUST=3

	//【1】计算描述子所需的图像领域区域的半径
	//考虑到要进行双线性插值，每个区域的宽度应为:DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )
	//在考虑到旋转因素，每个区域的宽度应为：DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )* sqrt(2)
	//所以搜索的半径是：DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )* sqrt(2) / 2
	int    radius = (sub_hist_width*sqrt(2.0)*(width + 1)) / 2.0 + 0.5;    //[1]0.5取四舍五入
	double grad_ori;
	double grad_mag;

	//遍历每个区域的像素
	for (int i = -radius; i <= radius; i++)
	{
		for (int j = -radius; j <= radius; j++)
		{
			//坐标旋转为主方向,以确保旋转不变性。
			double rot_x = (cos_ori * j - sin_ori * i) / sub_hist_width;
			double rot_y = (sin_ori * j + cos_ori * i) / sub_hist_width;

			//将邻域内的采样点分配到对应的子区域内，将子区域内的梯度值分配到8个方向上，计算其权值
			double xbin = rot_x + width / 2 - 0.5;                         //[2]xbin,ybin为落在4*4窗口中的下标值
			double ybin = rot_y + width / 2 - 0.5;

			if (xbin > -1.0 && xbin < width && ybin > -1.0 && ybin < width)
			{
				if (CalcGradMagOri(gauss, x + j, y + i, grad_mag, grad_ori)) //[3]计算关键点的梯度方向
				{
					grad_ori = (CV_PI - grad_ori) - ori;
					while (grad_ori < 0.0)
						grad_ori += PI2;
					while (grad_ori >= PI2)
						grad_ori -= PI2;              //把梯度方向归入[0，2Π]

					double obin = grad_ori * (bins / PI2);   
					//梯度方向*(8/2Π)归到[0，8]，8为直方图有8个bin，2Π为角度取值范围的长度，把方向的索引归于0~8.
					//conste = -1.0 / (2 * sigma*sigma);
					//计算权重
					double weight = exp(conste*(rot_x*rot_x + rot_y * rot_y));

					InterpHistEntry(hist, xbin, ybin, obin, grad_mag*weight, bins, width); //线性差值求每个种子点八个方向的梯度。

				}
			}
		}
	}

	return hist;
}

7.5、将8个方向的梯度变为特征向量

 7.6、归一化消除光照影响

 7.7、描述子门限，再次归一化

特征描述子

如上统计的4*4*8=128个梯度信息即为该关键点的特征向量。
      特征向量形成后，为了去除光照变化的影响，需要对它们进行归一化处理，对于图像灰度值整体漂移，图像各点的梯度是邻域像素相减得到，所以也能去除。得到的描述子向量为H=(h1,h2,.......,h128)，归一化后的特征向量为L=(L1,L2,......,L128)，则

描述子的门限化

非线性光照，相机饱和度变化对造成某些方向的梯度值过大，而对方向的影响微弱。因此设置门限值(向量归一化后，一般取0.2)截断较大的梯度值(大于0.2的则就令它等于0.2，小于0.2的则保持不变)。然后再进行一次归一化处理，提高特征的鉴别性。

 

//归一化特征点的特征描述子，即将特征描述子数组中每个元素除以特征描述子的模
void NormalizeDescr(Keypoint& feat)
{
	double cur, len_inv, len_sq = 0.0;
	int i, d = feat.descr_length;   //特征描述子的维数
 
//求特征描述子的模
	for (i = 0; i < d; i++)
	{
		cur = feat.descriptor[i];
		len_sq += cur*cur;
	}
	len_inv = 1.0 / sqrt(len_sq);  //sqrt( len_sq )为特征描述子的模
	for (i = 0; i < d; i++)
		feat.descriptor[i] *= len_inv; //特征描述子中每个元素除以特征描述子的模，完成归一化
}
/********************************************************************************************************************************
*模块说明:
*        模块七--步骤2:直方图到描述子的转换
*功能说明：
将某特征点的方向直方图转换为特征描述子向量，对特征描述子归一化并将所有元素转化为整型，存入指定特征点中
参数：
hist：width*width*bins的三维直方图数组
width：计算方向直方图时，将特征点附近划分为width*width个区域，每个区域生成一个直方图
bins：每个直方图的bin个数
feature：特征点指针，将计算好的特征描述子存入其中
*********************************************************************************************************************************/
void HistToDescriptor(double ***hist, int width, int bins, Keypoint& feature)
{
	int int_val, i, r, c, o, k = 0;
//遍历width*width*bins的三维直方图数组，将其中的所有数据(一般是128个)都存入feature结构的descriptor成员中
	for (r = 0; r < width; r++)
		for (c = 0; c < width; c++)
			for (o = 0; o < bins; o++)
			{
				feature.descriptor[k++] = hist[r][c][o];
			}
	
	feature.descr_length = k;  //特征描述子的维数，一般是128
	NormalizeDescr(feature);                           //[1]描述子特征向量归一化
  //归一化特征点的特征描述子，即将特征描述子数组中每个元素除以特征描述子的模

  //	DESCR_MAG_THR 0.2
  //遍历特征描述子向量，将超过阈值DESCR_MAG_THR的元素强行赋值DESCR_MAG_THR
	for (i = 0; i < k; i++)                           //[2]描述子向量门限
		if (feature.descriptor[i] > DESCR_MAG_THR)
			feature.descriptor[i] = DESCR_MAG_THR;

	NormalizeDescr(feature);                           //[3]描述子进行最后一次的归一化操作

//遍历特征描述子向量，每个元素乘以系数INT_DESCR_FCTR( 512.0)来变为整型，并且最大值不能超过255
	for (i = 0; i < k; i++)                           //[4]将单精度浮点型的描述子转换为整形的描述子
	{
		int_val = INT_DESCR_FCTR * feature.descriptor[i];
		feature.descriptor[i] = min(255, int_val);
	}
}　

上述3步原理参考我的博文：SIFT算法原理（3）-确定关键点的主方位，构建关键点描述符

如何画出已找到的关键点：

一些画出点的代码信息：

opencv数据结构CvScalar，【杂谈opencv】OpenCV中的cvRound()、cvFloor()、 cvCeil()函数讲解，

opencv学习中——CvPoint、CvSize、CvRect、CV_RGB、cvRectangle

代码实现：

/********************************************************************************************************************************
*模块说明:
*        模块七:6 关键点描述
*功能说明：
*        1）通过以上步骤，对于一个关键点，拥有三个信息：位置、尺度、方向
*        2）接下来就是为每个关键点建立一个描述符，用一组向量来将这个关键点描述出来，使其不随各种变化而变化，比如光照、视角变化等等
*        3）这个描述子不但包括关键点，也包含关键点周围对其贡献的像素点，并且描述符应该有较高的独特性，以便于特征点正确的匹配概率
*        1）SIFT描述子----是关键点邻域高斯图像梯度统计结果的一种表示。
*        2）通过对关键点周围图像区域分块，计算块内梯度直方图，生成具有独特性的向量
*        3）这个向量是该区域图像信息的一种表述和抽象，具有唯一性。
*Lowe论文：
*    Lowe建议描述子使用在关键点尺度空间内4*4的窗口中计算的8个方向的梯度信息，共4*4*8=128维向量来表征。具体的步骤如下所示:
*        1)确定计算描述子所需的图像区域
*        2）将坐标轴旋转为关键点的方向，以确保旋转不变性，如CSDN博文中的图6.2所示；旋转后邻域采样点的新坐标可以通过公式(6-2)计算
*        3）将邻域内的采样点分配到对应的子区域，将子区域内的梯度值分配到8个方向上，计算其权值
*        4）插值计算每个种子点八个方向的梯度
*        5）如上统计的4*4*8=128个梯度信息即为该关键点的特征向量。特征向量形成后，为了去除光照变化的影响，需要对它们进行归一化处理，
*           对于图像灰度值整体漂移，图像各点的梯度是邻域像素相减得到的，所以也能去除。得到的描述子向量为H，归一化之后的向量为L
*        6）描述子向量门限。非线性光照，相机饱和度变化对造成某些方向的梯度值过大，而对方向的影响微弱。因此，设置门限值（向量归一化
*           后，一般取0.2）截断较大的梯度值。然后，在进行一次归一化处理，提高特征的鉴别性。
*        7）按特征点的尺度对特征描述向量进行排序
*        8）至此，SIFT特征描述向量生成。
将关键点附近的区域划分为d*d(Lowe建议d=4)个子区域，每个子区域作为一个种子点，每个种子点有8个方向。
*********************************************************************************************************************************/
void DescriptorRepresentation(vector<Keypoint>& features, const vector<Mat>& gauss_pyr, int bins, int width)   //bins=8,width=4
{
	double ***hist;

	for (int i = 0; i < features.size(); i++)
	{                                                                       //[1]计算描述子的直方图
		//gauss_pyr[features[i].octave*(INTERVALS + 3) + features[i].interval  计算特征点所在矩阵中的位置；
		hist = CalculateDescrHist(gauss_pyr[features[i].octave*(INTERVALS + 3) + features[i].interval],
			features[i].x, features[i].y, features[i].octave_scale, features[i].ori, bins, width);

		HistToDescriptor(hist, width, bins, features[i]);                   //[2]直方图到描述子的转换

/*释放计算特征描述子过程中用到的方向直方图的内存空间
hist：方向直方图的指针，是一个width*width*bins的三维直方图数组
*/
		for (int j = 0; j < width; j++)
		{
			for (int k = 0; k < width; k++)
			{
				delete[] hist[j][k]; //释放第三维的内存
			}
			delete[] hist[j];   //释放第二维的内存
		}
		delete[] hist;  //释放第一维的内存
	}
}

上面多维数组可参考：图解c/c＋＋多级指针与“多维”数组

排序函数用法：
sort函数的用法(C++排序库函数的调用)

sift代码调用总：

/*******************************************************************************************************************
*函数说明:
*        最大的模块1：SIFT算法模块
*函数参数说明:
*        1---const Mat &src---------------准备进行特征点检测的原始图片
*        2---Vector<Keypoint>& features---用来存储检测出来的关键点
*        3---double sigma-----------------sigma值
*        4---int intervals----------------关键点所在的层数
********************************************************************************************************************/
void Sift(const Mat &src, vector<Keypoint>& features, double sigma, int intervals)
{
	std::cout << "【Step_one】Create -1 octave gaussian pyramid image" << std::endl;

	std::cout << "sigma      = " << sigma << std::endl;
	std::cout << "intervals     = " << intervals << std::endl;

	cv::Mat          init_gray;
	CreateInitSmoothGray(src, init_gray, sigma);   //灰度图
	int octaves = log((double)min(init_gray.rows, init_gray.cols)) / log(2.0) - 2;             //计算高斯金字塔的层数
	std::cout << "【1】The height and width of init_gray_img = " << init_gray.rows << "*" << init_gray.cols << std::endl;
	std::cout << "【2】The octaves of the gauss pyramid      = " << octaves << std::endl;


	std::cout << "【Step_two】Building gaussian pyramid ..." << std::endl;
	std::vector<Mat> gauss_pyr;
	GaussianPyramid(init_gray, gauss_pyr, octaves, intervals, sigma);    //高斯金字塔
	write_pyr(gauss_pyr, "gausspyramid");


	std::cout << "【Step_three】Building difference of gaussian pyramid..." << std::endl;
	vector<Mat> dog_pyr;
	DogPyramid(gauss_pyr, dog_pyr, octaves, intervals);     //差分金字塔
	write_pyr(dog_pyr, "dogpyramid");



	std::cout << "【Step_four】Deatecting local extrema..." << std::endl;
	vector<Keypoint> extrema;
	DetectionLocalExtrema(dog_pyr, extrema, octaves, intervals);   //极值点初探（排除较小点，找到极值点，消除边缘响应）
	std::cout << "【3】keypoints cout: " << extrema.size() << " --" << std::endl;
	std::cout << "【4】extrema detection finished." << std::endl;
	std::cout << "【5】please look dir gausspyramid, dogpyramid and extrema.txt.--" << endl;



	std::cout << "【Step_five】CalculateScale..." << std::endl;
	CalculateScale(extrema, sigma, intervals);   //计算特征点的尺度
	HalfFeatures(extrema);            //图像缩放



	std::cout << "【Step_six】Orientation assignment..." << std::endl;
	OrientationAssignment(extrema, features, gauss_pyr);    //关键点直方图计算，使用的extrema是缩放后的图像
	std::cout << "【6】features count: " << features.size() << std::endl;



	std::cout << "【Step_seven】DescriptorRepresentation..." << std::endl;
	//DESCR_HIST_BINS=8;DESCR_WINDOW_WIDTH=4;
	DescriptorRepresentation(features, gauss_pyr, DESCR_HIST_BINS, DESCR_WINDOW_WIDTH);  //关键点主方向的描述
	sort(features.begin(), features.end(), FeatureCmp);  //排序， FeatureCmp是排序方法
	cout << "finished." << endl;
}

　　

最后放一下总的代码

vs2015+0pencv3.2.0实现

sift.cpp

#include"sift.h"
#include<fstream>
#include<iostream>
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <opencv2/core/core.hpp>
#include <stdio.h>
#include <stdlib.h>

using namespace std;
using namespace cv;

/*************************************************************************************************************************
*模块说明：
*        主函数；
**************************************************************************************************************************/
int main(int argc, char **argv)
{
	cv::Mat src = imread("E:\VS2015Opencv\vs2015\project\picture\05.jpg");

	if (src.empty())
	{
		cout << "jpg open error! " << endl;
		getchar();
		return -1;
	}

	if (src.channels() != 3) return -2;

	vector<Keypoint> features;

	Sift(src, features, 1.6);                           //【1】SIFT特征点检测和特征点描述
	DrawKeyPoints(src, features);                       //【2】画出关键点(特征点)
	DrawSiftFeatures(src, features);                    //【3】画出SIFT特征点
	write_features(features, "descriptor.txt");         //【4】将特征点写入文本

	cv::imshow("src", src);
	cv::waitKey();

	return 0;
}

/*************************************************************************************************************************
*模块说明：
*        图像灰度化函数----将彩色图像转为灰度图像
**************************************************************************************************************************/
void ConvertToGray(const Mat& src, Mat& dst)
{
	cv::Size size = src.size();
	if (dst.empty())
		dst.create(size, CV_64F);                              //[1]利用Mat类的成员函数创建Mat容器

	uchar*    srcData = src.data;                             //[2]指向存储所有像素值的矩阵的数据区域
	pixel_t*  dstData = (pixel_t*)dst.data;                   //[3]指向dst的矩阵体数据区域

	int       dstStep = dst.step / sizeof(dstData[0]);         //[4]一行图像含有多少字节数

	for (int j = 0; j<src.cols; j++)
	{
		for (int i = 0; i<src.rows; i++)
		{
			double b = (srcData + src.step*i + src.channels()*j)[0] / 255.0;
			double g = (srcData + src.step*i + src.channels()*j)[1] / 255.0;
			double r = (srcData + src.step*i + src.channels()*j)[2] / 255.0;
			*(dstData + dstStep*i + dst.channels()*j) = (r + g + b) / 3.0;
		}
	}
}
/*************************************************************************************************************************
*模块说明：
*        隔点采样
**************************************************************************************************************************/
void DownSample(const Mat& src, Mat& dst)
{
	if (src.channels() != 1)
		return;

	if (src.cols <= 1 || src.rows <= 1)
	{
		src.copyTo(dst);
		return;
	}

	dst.create((int)(src.rows / 2), (int)(src.cols / 2), src.type());   //dst的行列=src/2

	pixel_t* srcData = (pixel_t*)src.data;
	pixel_t* dstData = (pixel_t*)dst.data;

	int srcStep = src.step / sizeof(srcData[0]);       //src一行图像含有多少字节数
	int dstStep = dst.step / sizeof(dstData[0]);

	int m = 0, n = 0;
	for (int j = 0; j < src.cols; j += 2, n++)
	{
		m = 0;
		for (int i = 0; i < src.rows; i += 2, m++)
		{
			pixel_t sample = *(srcData + srcStep * i + src.channels() * j);
			if (m < dst.rows && n < dst.cols)
			{
				*(dstData + dstStep * m + dst.channels() * n) = sample;    //将src隔点采样到dst中去。刚好src行列/2
			}
		}
	}

}
/*************************************************************************************************************************
*模块说明：
*        线性插值放大，将图像放大2倍
**************************************************************************************************************************/
void UpSample(const Mat &src, Mat &dst)
{
	if (src.channels() != 1)
		return;
	dst.create(src.rows * 2, src.cols * 2, src.type());    //dst的行列是原图2倍大小的行列

	pixel_t* srcData = (pixel_t*)src.data;   
	pixel_t* dstData = (pixel_t*)dst.data;

	int srcStep = src.step / sizeof(srcData[0]);  //一行图像含有多少字节数
	int dstStep = dst.step / sizeof(dstData[0]);

	int m = 0, n = 0;
	for (int j = 0; j < src.cols - 1; j++, n += 2)
	{
		m = 0;
		for (int i = 0; i < src.rows - 1; i++, m += 2)
		{
/*第1步：n=0 m=0，n=2,m=0;  第2步：n=0 m=1，n=2,m=1; 第3步：n=1 m=0，n=3,m=0;  第4步：n=1 m=1，n=3,m=1;    第5步：最后2行2列  */
/*第1步：n=0 m=2，n=2,m=2;  第2步：n=0 m=3，n=2,m=3; 第3步：n=1 m=2，n=3,m=2;  第4步：n=1 m=3，n=3,m=3;    假设5*5的矩阵      */
			double sample = *(srcData + srcStep * i + src.channels() * j);
			*(dstData + dstStep * m + dst.channels() * n) = sample;

			double rs = *(srcData + srcStep * (i)+src.channels()*j) + (*(srcData + srcStep * (i + 1) + src.channels()*j));
			*(dstData + dstStep * (m + 1) + dst.channels() * n) = rs / 2;
			double cs = *(srcData + srcStep * i + src.channels()*(j)) + (*(srcData + srcStep * i + src.channels()*(j + 1)));
			*(dstData + dstStep * m + dst.channels() * (n + 1)) = cs / 2;

			double center = (*(srcData + srcStep * (i + 1) + src.channels() * j))
				+ (*(srcData + srcStep * i + src.channels() * j))
				+ (*(srcData + srcStep * (i + 1) + src.channels() * (j + 1)))
				+ (*(srcData + srcStep * i + src.channels() * (j + 1)));

			*(dstData + dstStep * (m + 1) + dst.channels() * (n + 1)) = center / 4;

		}

	}



	if (dst.rows < 3 || dst.cols < 3)
		return;

	//最后两行两列
	for (int k = dst.rows - 1; k >= 0; k--)
	{
		*(dstData + dstStep *(k)+dst.channels()*(dst.cols - 2)) = *(dstData + dstStep *(k)+dst.channels()*(dst.cols - 3));
		*(dstData + dstStep *(k)+dst.channels()*(dst.cols - 1)) = *(dstData + dstStep *(k)+dst.channels()*(dst.cols - 3));
	}
	for (int k = dst.cols - 1; k >= 0; k--)
	{
		*(dstData + dstStep *(dst.rows - 2) + dst.channels()*(k)) = *(dstData + dstStep *(dst.rows - 3) + dst.channels()*(k));
		*(dstData + dstStep *(dst.rows - 1) + dst.channels()*(k)) = *(dstData + dstStep *(dst.rows - 3) + dst.channels()*(k));
	}

}

/*************************************************************************************************************************
*模块说明：
*        OpenCv中的高斯平滑函数
**************************************************************************************************************************/
void GaussianSmooth(const Mat &src, Mat &dst, double sigma)
{
	GaussianBlur(src, dst, Size(0, 0), sigma, sigma);
}
/*************************************************************************************************************************
*模块说明：
*        模块1--创建初始灰度图像------初始图像先将原图像灰度化，再扩大一倍后，使用高斯模糊平滑
*功能说明：
*        在最开始建立高斯金字塔的时候，要预先模糊输入的图像来作为第0个组的第0层图像，这时相当于丢弃了最高的空域的采样率。因
*        此，通常的做法是先将图像的尺度扩大一倍来生成第-1组。我们假定初始的输入图像为了抗击混淆现象，已经对其进行了sigma=0.5
*        的高斯模糊，如果输入图像的尺寸用双线性插值扩大一倍，那么相当于sigma=1.0
*这样做的原因:
*        在检测极值点前，对原始图像的高斯平滑以致图像丢失高频信息，所以Lowe建议在建立尺度空间前，首先对原始图像长宽扩展一倍，
*        以便保留原始图像的信息(特别是图像的高频信息，比如边缘，角点)，增加特征点的数量。
*函数功能:
*        这个函数的作用是用于创建高斯金字塔的-1层图像，在主函数中定义 sigma =1.6
在lowe的论文中，他定义图片的尺度是1.6，被摄像头模糊的尺度是0.5，但是这个1.52对我们来说并没有啥意思，
lowe为了获得更多的特征点（姑且这么解释吧）他先对原始图像进行扩大一倍。然后呢再对它进行高斯平滑作为第一阶第一层，
高斯模糊尺度为=(1.6* 1.6 - (0.5 * 2) * (0.5 * 2))
**************************************************************************************************************************/
void CreateInitSmoothGray(const Mat &src, Mat &dst, double sigma = SIGMA)
{
	cv::Mat gray;                                  //[1]保存原始图像的灰度图像        
	cv::Mat up;                                    //[2]保存原始图像的上采样图像

	//std::cout << "CreateInitSmoothGray sigma     = " << sigma << std::endl;  //查看sigma初始值

	ConvertToGray(src, gray);                     //灰度图像
	UpSample(gray, up);                            //[3]图像的上采样
												   //[4]高斯金字塔的-1组的sigma的计算公式
	double  sigma_init = sqrt(sigma * sigma - (INIT_SIGMA * 2) * (INIT_SIGMA * 2));// INIT_SIGM初始sigma=0.5;  -1层的sigma

	GaussianSmooth(up, dst, sigma_init);           //[5]高斯平滑
}
/*************************************************************************************************************************
*模块说明：
*        模块二：3.3 图像高斯金字塔的构建
在主函数中：int octaves = log((double)min(init_gray.rows, init_gray.cols)) / log(2.0) - 2;             //计算高斯金字塔的层数
高斯金字塔中每组图像中有三层/张图片 INTERVALS=3；SIGMA=1.6；
**************************************************************************************************************************/
void GaussianPyramid(const Mat &src, vector<Mat>&gauss_pyr, int octaves, int intervals = INTERVALS, double sigma = SIGMA)
{
	//std::cout << "sigma is     = " << sigma << std::endl;
	double *sigmas = new double[intervals + 3];
	double k = pow(2.0, 1.0 / intervals);      //pow（x,y）,x的y次方；k=2的（1/s)次方
//	std::cout << "k is     = " << k << std::endl;
	sigmas[0] = sigma;
//	std::cout << "sigma0 is     = " << sigmas[0] << std::endl;
	double sig_prev;
	double sig_total;
//用循环来表示第0组内各层图片的尺度，这个可以表示为0组的各层图像尺度；
	for (int i = 1; i < intervals + 3; i++)
	{
		sig_prev = pow(k, i - 1) * sigma;     //第i-1层尺度
		sig_total = sig_prev * k;             //第i层尺度
		sigmas[i] = sqrt(sig_total * sig_total - sig_prev * sig_prev);  //组内每层尺度坐标
	}

	for (int o = 0; o < octaves; o++)
	{
		//每组多三层
		for (int i = 0; i < intervals + 3; i++)
		{
			Mat mat;
			if (o == 0 && i == 0)    //第0组，第0层，就是原图像
			{
				src.copyTo(mat);
			}
			else if (i == 0)           //新的一组的首层图像是由上一组最后一层图像下采样得到
			{
				DownSample(gauss_pyr[(o - 1)*(intervals + 3) + intervals], mat); 
			}
			else        //对上一层图像gauss_pyr[i-1]进行参数为sig[i]的高斯平滑，得到当前层图像
			{ 
				GaussianSmooth(gauss_pyr[o * (intervals + 3) + i - 1], mat, sigmas[i]);
			//	std::cout << "sigmas[i] is     = " << sigmas[i] << std::endl;
			}
			gauss_pyr.push_back(mat);  //将mat中的像素添加到gauss_pyr数组里面；也就是说每组从第0层到第5层的图片像素装进vector数组中；
	
		}
		
	}

	delete[] sigmas;
}
/*************************************************************************************************************************
*模块说明：
*        图像的差分
**************************************************************************************************************************/
void Sub(const Mat& a, const Mat& b, Mat & c)
{
	if (a.rows != b.rows || a.cols != b.cols || a.type() != b.type())
		return;
	if (!c.empty())
		return;
	c.create(a.size(), a.type());

	pixel_t* ap = (pixel_t*)a.data;
	pixel_t* ap_end = (pixel_t*)a.dataend;
	pixel_t* bp = (pixel_t*)b.data;
	pixel_t* cp = (pixel_t*)c.data;
	int step = a.step / sizeof(pixel_t);

	while (ap != ap_end)
	{
		*cp++ = *ap++ - *bp++;
	}
}
/*************************************************************************************************************************
*模块说明：
*       模块三：3.4 高斯差分金字塔
*功能说明：
*       1--2002年，Mikolajczyk在详细的实验比较中发现尺度归一化的高斯拉普拉斯函数的极大值和极小值同其他的特征提取函数，例如：
*          梯度，Hessian或者Harris角点特征比较，能够产生稳定的图像特征。
*       2--而Lindberg早在1994年发现高斯差分函数（Difference of Gaussian,简称DOG算子）与尺度归一化的拉普拉斯函数非常相似，因此，
*          利用图像间的差分代替微分。
**************************************************************************************************************************/
void DogPyramid(const vector<Mat>& gauss_pyr, vector<Mat>& dog_pyr, int octaves, int intervals = INTERVALS)
{
	for (int o = 0; o < octaves; o++)
	{
		for (int i = 1; i < intervals + 3; i++)
		{
			Mat mat;
			Sub(gauss_pyr[o*(intervals + 3) + i], gauss_pyr[o*(intervals + 3) + i - 1], mat);
			dog_pyr.push_back(mat);
		}
	}
}
/*************************************************************************************************************************
*模块说明：
*       模块四的第一步:3.4-空间极值点的初步检测(关键点的初步探查)
*功能说明:
*       1--在高斯差分函数之后，得到一系列的关键点的疑似点，我们需要对这些关键点的疑似点初步进行检测和筛选
*      检测是以前点为中心，3pixel*3pixel*3pixel的立方体为邻域，判断当前点是否为局部最大或最小。i,j,k=-1,0,1表示当前点的邻域
**************************************************************************************************************************/
bool isExtremum(int x, int y, const vector<Mat>& dog_pyr, int index)
{
	pixel_t * data = (pixel_t *)dog_pyr[index].data;
	int      step = dog_pyr[index].step / sizeof(data[0]);
	pixel_t   val = *(data + y*step + x);

	if (val > 0)
	{
		for (int i = -1; i <= 1; i++)    //层
		{
			int stp = dog_pyr[index + i].step / sizeof(data[0]);
			for (int j = -1; j <= 1; j++)     //行
			{
				for (int k = -1; k <= 1; k++)  //列
				{
					if (val < *((pixel_t*)dog_pyr[index + i].data + stp*(y + j) + (x + k)))
					{
						return false;   //不是极值点退出
					}
				}
			}
		}
	}
	else
	{
		for (int i = -1; i <= 1; i++)
		{
			int stp = dog_pyr[index + i].step / sizeof(pixel_t);
			for (int j = -1; j <= 1; j++)
			{
				for (int k = -1; k <= 1; k++)
				{
					//检查最小极值
					if (val > *((pixel_t*)dog_pyr[index + i].data + stp*(y + j) + (x + k)))
					{
						return false;  //不是极值点退出
					}
				}
			}
		}
	}
	return true;  //是级值点，返回
}
/*************************************************************************************************************************
*模块说明：
*       模块四的第三步:4.2--消除边缘响应点
*功能说明:
*       1）一个定义不好的高斯差分算子的极值在横跨边缘的地方有较大的住主曲率，在垂直边缘的方向有较小的主曲率。
*       2）DOG算子会产生较强的边缘响应，需要剔除不稳定的边缘响应点，获取特征点处的Hessian矩阵，主曲率通过一个2*2的Hessian矩
*          阵H求出
*       3）主曲率D和Hessian矩阵的特征值成正比，公式(r+1)*(r+1)/r的值在两个特征值相等时最小；这个值越大，说明两个特征值的比值
*          越大，即在某一个方向的梯度值越大，而在另一个方向的梯度值越小，而边缘恰恰就是这种情况。所以，为了剔除边缘响应点，
*          需要让该比值小于一定的阈值，因此，为了检测主曲率是否在某阈值r下，只需检测。CSDN论文中的公式(4-7成立)，成立，将关
*          键点保留，反之，剔除。
**************************************************************************************************************************/
#define DAt(x, y) (*(data+(y)*step+(x))) 
bool passEdgeResponse(int x, int y, const vector<Mat>& dog_pyr, int index, double r = RATIO)
{
	pixel_t *data = (pixel_t *)dog_pyr[index].data;  //极值点的图片索引定位
	int step = dog_pyr[index].step / sizeof(data[0]);
	pixel_t val = *(data + y*step + x);    //极值点

	double Dxx;
	double Dyy;
	double Dxy;
	double Tr_h;                                                         //[1]Hessian矩阵的迹
	double Det_h;                                                        //[2]Hessian矩阵所对应的行列式的值

	Dxx = DAt(x + 1, y) + DAt(x - 1, y) - 2 * val;
	Dyy = DAt(x, y + 1) + DAt(x, y - 1) - 2 * val;
	Dxy = (DAt(x + 1, y + 1) + DAt(x - 1, y - 1) - DAt(x - 1, y + 1) - DAt(x + 1, y - 1)) / 4.0;

	Tr_h = Dxx + Dyy;
	Det_h = Dxx * Dyy - Dxy * Dxy;

	if (Det_h <= 0)return false;

	if (Tr_h * Tr_h / Det_h < (r + 1) * (r + 1) / r) return true;

	return false;
}
/*************************************************************************************************************************
*模块说明：
*       有限差分求导？
**************************************************************************************************************************/
#define Hat(i, j) (*(H+(i)*3 + (j)))

double PyrAt(const vector<Mat>& pyr, int index, int x, int y)
{
	pixel_t *data = (pixel_t*)pyr[index].data;
	int      step = pyr[index].step / sizeof(data[0]);
	pixel_t   val = *(data + y*step + x);

	return val;
}
/*************************************************************************************************************************
*模块说明：
*       有限差分求导？
**************************************************************************************************************************/
#define At(index, x, y) (PyrAt(dog_pyr, (index), (x), (y)))

//3维D(x)一阶偏导,dx列向量
void DerivativeOf3D(int x, int y, const vector<Mat>& dog_pyr, int index, double *dx)
{
	double Dx = (At(index, x + 1, y) - At(index, x - 1, y)) / 2.0;
	double Dy = (At(index, x, y + 1) - At(index, x, y - 1)) / 2.0;
	double Ds = (At(index + 1, x, y) - At(index - 1, x, y)) / 2.0;

	dx[0] = Dx;
	dx[1] = Dy;
	dx[2] = Ds;
}

//3维D(x)二阶偏导，即Hessian矩阵
void Hessian3D(int x, int y, const vector<Mat>& dog_pyr, int index, double *H)
{
	double val, Dxx, Dyy, Dss, Dxy, Dxs, Dys;

	val = At(index, x, y);

	Dxx = At(index, x + 1, y) + At(index, x - 1, y) - 2 * val;
	Dyy = At(index, x, y + 1) + At(index, x, y - 1) - 2 * val;
	Dss = At(index + 1, x, y) + At(index - 1, x, y) - 2 * val;

	Dxy = (At(index, x + 1, y + 1) + At(index, x - 1, y - 1)
		- At(index, x + 1, y - 1) - At(index, x - 1, y + 1)) / 4.0;

	Dxs = (At(index + 1, x + 1, y) + At(index - 1, x - 1, y)
		- At(index - 1, x + 1, y) - At(index + 1, x - 1, y)) / 4.0;

	Dys = (At(index + 1, x, y + 1) + At(index - 1, x, y - 1)
		- At(index + 1, x, y - 1) - At(index - 1, x, y + 1)) / 4.0;

	Hat(0, 0) = Dxx;
	Hat(1, 1) = Dyy;
	Hat(2, 2) = Dss;

	Hat(1, 0) = Hat(0, 1) = Dxy;
	Hat(2, 0) = Hat(0, 2) = Dxs;
	Hat(2, 1) = Hat(1, 2) = Dys;
}
/*************************************************************************************************************************
*模块说明：
*       4.4 三阶矩阵求逆
**************************************************************************************************************************/
#define HIat(i, j) (*(H_inve+(i)*3 + (j)))
//3*3阶矩阵求逆
bool Inverse3D(const double *H, double *H_inve)
{

	double A = Hat(0, 0)*Hat(1, 1)*Hat(2, 2)
		+ Hat(0, 1)*Hat(1, 2)*Hat(2, 0)
		+ Hat(0, 2)*Hat(1, 0)*Hat(2, 1)
		- Hat(0, 0)*Hat(1, 2)*Hat(2, 1)
		- Hat(0, 1)*Hat(1, 0)*Hat(2, 2)
		- Hat(0, 2)*Hat(1, 1)*Hat(2, 0);

	if (fabs(A) < 1e-10) return false;

	HIat(0, 0) = Hat(1, 1) * Hat(2, 2) - Hat(2, 1)*Hat(1, 2);
	HIat(0, 1) = -(Hat(0, 1) * Hat(2, 2) - Hat(2, 1) * Hat(0, 2));
	HIat(0, 2) = Hat(0, 1) * Hat(1, 2) - Hat(0, 2)*Hat(1, 1);

	HIat(1, 0) = Hat(1, 2) * Hat(2, 0) - Hat(2, 2)*Hat(1, 0);
	HIat(1, 1) = -(Hat(0, 2) * Hat(2, 0) - Hat(0, 0) * Hat(2, 2));
	HIat(1, 2) = Hat(0, 2) * Hat(1, 0) - Hat(0, 0)*Hat(1, 2);

	HIat(2, 0) = Hat(1, 0) * Hat(2, 1) - Hat(1, 1)*Hat(2, 0);
	HIat(2, 1) = -(Hat(0, 0) * Hat(2, 1) - Hat(0, 1) * Hat(2, 0));
	HIat(2, 2) = Hat(0, 0) * Hat(1, 1) - Hat(0, 1)*Hat(1, 0);

	for (int i = 0; i < 9; i++)
	{
		*(H_inve + i) /= A;
	}
	return true;
}
/*************************************************************************************************************************
*模块说明：
*
**************************************************************************************************************************/
//计算x^
void GetOffsetX(int x, int y, const vector<Mat>& dog_pyr, int index, double *offset_x)
{
	//x^ = -H^(-1) * dx; dx = (Dx, Dy, Ds)^T
	double H[9], H_inve[9] = { 0 };
	Hessian3D(x, y, dog_pyr, index, H);
	Inverse3D(H, H_inve);
	double dx[3];
	DerivativeOf3D(x, y, dog_pyr, index, dx);

	for (int i = 0; i < 3; i++)
	{
		offset_x[i] = 0.0;
		for (int j = 0; j < 3; j++)
		{
			offset_x[i] += H_inve[i * 3 + j] * dx[j];
		}
		offset_x[i] = -offset_x[i];
	}
}

//计算|D(x^)|
double GetFabsDx(int x, int y, const vector<Mat>& dog_pyr, int index, const double* offset_x)
{
	//|D(x^)|=D + 0.5 * dx * offset_x; dx=(Dx, Dy, Ds)^T
	double dx[3];
	DerivativeOf3D(x, y, dog_pyr, index, dx);

	double term = 0.0;
	for (int i = 0; i < 3; i++)
		term += dx[i] * offset_x[i];

	pixel_t *data = (pixel_t *)dog_pyr[index].data;
	int step = dog_pyr[index].step / sizeof(data[0]);
	pixel_t val = *(data + y*step + x);

	return fabs(val + 0.5 * term);
}
/*************************************************************************************************************************
*模块说明：
*       模块四的第二步:修正极值点，删除不稳定的点
*功能说明:
*       1--根据高斯差分函数产生的极值点并不全都是稳定的特征点，因为某些极值点的响应较弱，而且DOG算子会产生较强的边缘响应
*       2--以上方法检测到的极值点是离散空间的极值点，下面通过拟合三维二次函数来精确定位关键点的位置和尺度，同时去除对比度
*          低和不稳定的边缘响应点(因为DOG算子会产生较强的边缘响应)，以增强匹配的稳定性、提高抗噪声的能力。
*       3--修正极值点，删除不稳定点，|D(x)| < 0.03 Lowe 2004
**************************************************************************************************************************/
Keypoint* InterploationExtremum(int x, int y, const vector<Mat>& dog_pyr, int index, int octave, int interval, double dxthreshold = DXTHRESHOLD)
{
	//计算x=(x,y,sigma)^T
	//x^ = -H^(-1) * dx; dx = (Dx, Dy, Ds)^T
	double offset_x[3] = { 0 };

	const Mat &mat = dog_pyr[index];

	int idx = index;
	int intvl = interval;
	int i = 0;

	while (i < MAX_INTERPOLATION_STEPS)
	{
		GetOffsetX(x, y, dog_pyr, idx, offset_x);
		//4. Accurate keypoint localization.  Lowe
		//如果offset_x 的任一维度大于0.5，it means that the extremum lies closer to a different sample point.
		if (fabs(offset_x[0]) < 0.5 && fabs(offset_x[1]) < 0.5 && fabs(offset_x[2]) < 0.5)
			break;

		//用周围的点代替
		x += cvRound(offset_x[0]);
		y += cvRound(offset_x[1]);
		interval += cvRound(offset_x[2]);

		idx = index - intvl + interval;
		//此处保证检测边时 x+1,y+1和x-1, y-1有效
		if (interval < 1 || interval > INTERVALS || x >= mat.cols - 1 || x < 2 || y >= mat.rows - 1 || y < 2)
		{
			return NULL;
		}

		i++;
	}

	//窜改失败
	if (i >= MAX_INTERPOLATION_STEPS)
		return NULL;

	//rejecting unstable extrema
	//|D(x^)| < 0.03取经验值
	if (GetFabsDx(x, y, dog_pyr, idx, offset_x) < dxthreshold / INTERVALS)
	{
		return NULL;
	}

	Keypoint *keypoint = new Keypoint;


	keypoint->x = x;
	keypoint->y = y;

	keypoint->offset_x = offset_x[0];
	keypoint->offset_y = offset_x[1];

	keypoint->interval = interval;
	keypoint->offset_interval = offset_x[2];

	keypoint->octave = octave;

	keypoint->dx = (x + offset_x[0])*pow(2.0, octave);
	keypoint->dy = (y + offset_x[1])*pow(2.0, octave);

	return keypoint;
}
/*************************************************************************************************************************
*模块说明：
*       模块四：3.5 空间极值点的检测(关键点的初步探查)
*功能说明：
*       1--关键点是由DOG空间的局部极值点组成的，关键点的初步探查是通过同一组内各DoG相邻两层图像之间的比较完成的。为了寻找DoG
*          函数的极值点，每一个像素点都要和它所有相邻的点比较，看其是否比它的图像域和尺度域相邻的点大还是小。
*       2--当然这样产生的极值点并不全都是稳定的特征点，因为某些极值点相应较弱，而且DOG算子会产生较强的边缘响应。
**************************************************************************************************************************/
void DetectionLocalExtrema(const vector<Mat>& dog_pyr, vector<Keypoint>& extrema, int octaves, int intervals = INTERVALS)
{

	double  thresh = 0.5 * DXTHRESHOLD / intervals;  //|D(x^)| < 0.03 ; DXTHRESHOLD=0.03极值点太多

	for (int o = 0; o < octaves; o++)
	{
		//第一层和最后一层极值忽略
		for (int i = 1; i < (intervals + 2) - 1; i++)             //Dog图像每一组5层（差分6-1），又第一层和最后一程忽略；
		{
			int index = o*(intervals + 2) + i;                              //[1]图片索引的定位
			pixel_t *data = (pixel_t *)dog_pyr[index].data;                //[2]获取图片的矩阵体的首地址
		//	std::cout << "首地址为      = " << sizeof(data[0]) << std::endl;

			int step = dog_pyr[index].step / sizeof(data[0]);           //[3]说明矩阵在存储空间中的存储是以线性空间的方式存放的


			for (int y = IMG_BORDER; y < dog_pyr[index].rows - IMG_BORDER; y++)   //IMG_BORDER=5 ，这里我查找资料说的是指检测边界以内的点，因为边界易受干扰
			{
				for (int x = IMG_BORDER; x < dog_pyr[index].cols - IMG_BORDER; x++)
				{
				//	std::cout << "x is     = " << x << std::endl;
					pixel_t val = *(data + y*step + x);                     //遍历像素点
					if (std::fabs(val) > thresh)                           //[4]排除阈值过小的点，fabs（）取绝对值
					{
						if (isExtremum(x, y, dog_pyr, index))                //[5]判断是否是极值
						{
							Keypoint *extrmum = InterploationExtremum(x, y, dog_pyr, index, o, i);//极值点x,y；当前组当前层；
							if (extrmum)
							{
								if (passEdgeResponse(extrmum->x, extrmum->y, dog_pyr, index))
								{
									extrmum->val = *(data + extrmum->y*step + extrmum->x);
									extrema.push_back(*extrmum);
								}

								delete extrmum;

							}//extrmum
						}//isExtremum
					}//std::fabs
				}//for x
			}//for y

		}
	}
}
/*************************************************************************************************************************
*模块说明：
*       模块五：计算特征点对应的尺度

*
**************************************************************************************************************************/
void CalculateScale(vector<Keypoint>& features, double sigma = SIGMA, int intervals = INTERVALS)
{
	double intvl = 0;
	for (int i = 0; i < features.size(); i++)
	{
		intvl = features[i].interval + features[i].offset_interval;
		features[i].scale = sigma * pow(2.0, features[i].octave + intvl / intervals);
		features[i].octave_scale = sigma * pow(2.0, intvl / intervals);
	}

}

//对扩大的图像特征缩放
/*将特征点序列中每个特征点的坐标减半(图像金字塔设置了将图像放大为原图的2倍时，特征点检测完之后调用)
/features：特征点序列*/
void HalfFeatures(vector<Keypoint>& features)
{
	for (int i = 0; i < features.size(); i++)
	{
		features[i].dx /= 2;
		features[i].dy /= 2;

		features[i].scale /= 2;
	}
}
/********************************************************************************************************************************
*模块说明:
*        模块六---步骤2：计算关键点的梯度和梯度方向
*功能说明：
*        1）计算关键点(x,y)处的梯度幅值和梯度方向
*        2）将所计算出来的梯度幅值和梯度方向保存在变量mag和ori中
*********************************************************************************************************************************/
bool CalcGradMagOri(const Mat& gauss, int x, int y, double& mag, double& ori)
{
	if (x > 0 && x < gauss.cols - 1 && y > 0 && y < gauss.rows - 1)
	{
		pixel_t *data = (pixel_t*)gauss.data;
		int step = gauss.step / sizeof(*data);

		double dx = *(data + step*y + (x + 1)) - (*(data + step*y + (x - 1)));           //[1]利用X方向上的差分代替微分dx
		double dy = *(data + step*(y + 1) + x) - (*(data + step*(y - 1) + x));           //[2]利用Y方向上的差分代替微分dy

		mag = sqrt(dx*dx + dy*dy);                                          //[3]计算该关键点的梯度幅值
		ori = atan2(dy, dx);                                                //[4]计算该关键点的梯度方向
		return true;
	}
	else
		return false;
}
/********************************************************************************************************************************
*模块说明:
*        模块六---步骤1：计算梯度的方向直方图
*功能说明：
*        1）直方图以每10度为一个柱，共36个柱，柱代表的方向为为像素点的梯度方向，柱的长短代表了梯度幅值。
*        2）根据Lowe的建议，直方图统计采用3*1.5*sigma
*        3）在直方图统计时，每相邻三个像素点采用高斯加权，根据Lowe的建议，模板采用[0.25,0.5,0.25],并且连续加权两次
*结    论：
*        图像的关键点检测完毕后，每个关键点就拥有三个信息：位置、尺度、方向；同时也就使关键点具备平移、缩放和旋转不变性
bin=36.
*********************************************************************************************************************************/
double* CalculateOrientationHistogram(const Mat& gauss, int x, int y, int bins, int radius, double sigma)
{
	//std::cout << "首地址1为      = " << radius << std::endl;  //radius=3*1.5?*?σ(extrema[i].octave_scale)
    //std::cout << "首地址2为      = " << bins << std::endl;   // sigma=1.5*σ(extrema[i].octave_scale)
	double* hist = new double[bins];                           //[1]动态分配一个double类型的数组，用来存放梯度直方图
	for (int i = 0; i < bins; i++)                               //[2]给这个数组初始化
		*(hist + i) = 0.0;

	double  mag;                                                //[3]关键点的梯度幅值                                          
	double  ori;                                                //[4]关键点的梯度方向
	double  weight;

	int           bin;
	const double PI2 = 2.0*CV_PI;  //2Π
	double       econs = -1.0 / (2.0*sigma*sigma);

	for (int i = -radius; i <= radius; i++)
	{
		for (int j = -radius; j <= radius; j++)
		{
			if (CalcGradMagOri(gauss, x + i, y + j, mag, ori))       //[5]计算该关键点的梯度幅值和方向
			{
			//	std::cout << "ori为      = " << ori << std::endl;
				weight = exp((i*i + j*j)*econs);              // 该点的梯度幅值权重
				bin = cvRound(bins * (CV_PI - ori) / PI2);     //[6]对一个double行的数进行四舍五入，返回一个整形的数;
			//	bin = cvRound(bins * (CV_PI + ori) / PI2);     //[6]对一个double行的数进行四舍五入，返回一个整形的数;
				//计算梯度的方向对应的直方图中的bin下标,ori范围(-Π，Π]，+Π后范围(0，2Π]
			//	std::cout << "bin1为      = " << bin << std::endl;
				bin = bin < bins ? bin : 0;            //把bin的下标归到0~36；bin<bins，bin=bin;否则bin=0；
				//std::cout << "bin2为      = " << bin << std::endl;

				hist[bin] += mag * weight;                      //[7]统计梯度的方向直方图
			}
		}
	}

	return hist;
}
/********************************************************************************************************************************
*模块说明:
*        模块六---步骤3：对梯度方向直方图进行连续两次的高斯平滑
*功能说明：
*        1）在直方图统计时，每相邻三个像素点采用高斯加权，根据Lowe的建议，模板采用[0.25,0.5,0.25],并且连续加权两次
*        2）对直方图进行两次平滑
*********************************************************************************************************************************/
void GaussSmoothOriHist(double *hist, int n)
{
	double prev = hist[n - 1];
	double temp;
	double h0 = hist[0];

	for (int i = 0; i < n; i++)
	{
		temp = hist[i];
		hist[i] = 0.25 * prev + 0.5 * hist[i] + 0.25 * (i + 1 >= n ? h0 : hist[i + 1]);//对方向直方图进行高斯平滑
		prev = temp;
	}
}
/********************************************************************************************************************************
*模块说明:
*        模块六---步骤4：计算方向直方图中的主方向
*********************************************************************************************************************************/
double DominantDirection(double *hist, int n)
{
	double maxd = hist[0];
	for (int i = 1; i < n; i++)
	{
		if (hist[i] > maxd)                            //求取36个柱中的最大峰值
			maxd = hist[i];
	}
	return maxd;
}
//复制关键点
void CopyKeypoint(const Keypoint& src, Keypoint& dst)
{
	dst.dx = src.dx;
	dst.dy = src.dy;

	dst.interval = src.interval;
	dst.octave = src.octave;
	dst.octave_scale = src.octave_scale;
	dst.offset_interval = src.offset_interval;

	dst.offset_x = src.offset_x;
	dst.offset_y = src.offset_y;

	dst.ori = src.ori;
	dst.scale = src.scale;
	dst.val = src.val;
	dst.x = src.x;
	dst.y = src.y;
}
/********************************************************************************************************************************
*模块说明:
*        模块六---步骤5：计算更加精确的关键点主方向----抛物插值
*功能说明：
*        1）方向直方图的峰值则代表了该特征点的方向，以直方图中的最大值作为该关键点的主方向。为了增强匹配的鲁棒性，只保留峰值大于主
*           方向峰值80%的方向作为改关键点的辅方向。因此，对于同一梯度值得多个峰值的关键点位置，在相同位置和尺度将会有多个关键点被
*           创建但方向不同。仅有15%的关键点被赋予多个方向，但是可以明显的提高关键点的稳定性。
*        2）在实际编程中，就是把该关键点复制成多份关键点，并将方向值分别赋给这些复制后的关键点
*        3）并且，离散的梯度直方图要进行【插值拟合处理】，来求得更加精确的方向角度值
*********************************************************************************************************************************/
#define Parabola_Interpolate(l, c, r) (0.5*((l)-(r))/((l)-2.0*(c)+(r))) //求小柱子在直方图中的索引号
void CalcOriFeatures(const Keypoint& keypoint, vector<Keypoint>& features, const double *hist, int n, double mag_thr)
{
	double  bin;
	double  PI2 = CV_PI * 2.0;
	int l;
	int r;
	//遍历直方图
	for (int i = 0; i < n; i++)
	{
		l = (i == 0) ? n - 1 : i - 1; //前一个(左边的)bin的下标，
		r = (i + 1) % n;              //后一个(右边的)bin的下标
		//i=0,l=35,r=1;

		//hist[i]是极值
		//若当前的bin是局部极值(比前一个和后一个bin都大)，并且值大于给定的幅值阈值，则新生成一个特征点并添加到特征点序列末尾
		if (hist[i] > hist[l] && hist[i] > hist[r] && hist[i] >= mag_thr)  //mag_thr=0.8*关键点最大值
		{
			//根据左、中、右三个bin的值对当前bin进行直方图插值
			bin = i + Parabola_Interpolate(hist[l], hist[i], hist[r]);
			bin = (bin < 0) ? (bin + n) : (bin >= n ? (bin - n) : bin); //将插值结果规范到[0,n]内

			Keypoint new_key;

			CopyKeypoint(keypoint, new_key); //复制特征点

			new_key.ori = ((PI2 * bin) / n) - CV_PI; //方向还原为（-Π，Π]，bin最大值36，2Π*36/36-Π=Π，bin最小值0，2Π*0/36-Π=-Π
			features.push_back(new_key); //插入到特征点序列末尾
		}
	}
}
/********************************************************************************************************************************
*模块说明:
*        模块六：5 关键点方向分配
*功能说明：
*        1）为了使描述符具有旋转不变性，需要利用图像的局部特征为每一个关键点分配一个基准方向。使用图像梯度的方法求取局部结构的稳定
*           方向。
*        2）对于在DOG金字塔中检测出来的关键点，采集其所在高斯金字塔图像3sigma邻域窗口内像素的梯度和方向梯度和方向特征。
*        3）梯度的模和方向如下所示:
*        4) 在完成关键点的梯度计算后，使用直方图统计邻域内像素的梯度和方向。梯度直方图将0~360度的方向范围分为36个柱，其中每柱10度，
*           如图5.1所示，直方图的峰值方向代表了关键点的主方向
         cvRound()：返回跟参数最接近的整数值，即四舍五入；
*********************************************************************************************************************************/
void OrientationAssignment(vector<Keypoint>& extrema, vector<Keypoint>& features, const vector<Mat>& gauss_pyr)
{
	int n = extrema.size();
	double *hist;

	for (int i = 0; i < n; i++)
	{
		//ORI_HIST_BINS=36;ORI_WINDOW_RADIUS=3.0*1.5,特征点方向赋值过程中，搜索邻域的半径为：3*1.5*σ
		//ORI_SIGMA_TIMES=1.5,
		//gauss_pyr[extrema[i].octave*(INTERVALS + 3) + extrema[i].interval是在计算关键点在数组中的位置；
		hist = CalculateOrientationHistogram(gauss_pyr[extrema[i].octave*(INTERVALS + 3) + extrema[i].interval],
			extrema[i].x, extrema[i].y, ORI_HIST_BINS, cvRound(ORI_WINDOW_RADIUS*extrema[i].octave_scale),
			ORI_SIGMA_TIMES*extrema[i].octave_scale);                             //[1]计算梯度的方向直方图

		for (int j = 0; j < ORI_SMOOTH_TIMES; j++)
			GaussSmoothOriHist(hist, ORI_HIST_BINS);                              //[2]对方向直方图进行高斯平滑
		double highest_peak = DominantDirection(hist, ORI_HIST_BINS);            //[3]求取方向直方图中的峰值
																				 //[4]计算更加精确的关键点主方向
		CalcOriFeatures(extrema[i], features, hist, ORI_HIST_BINS, highest_peak*ORI_PEAK_RATIO); //ORI_PEAK_RATIO=0.8
		//精确关键点的主方向，抛物插值
		delete[] hist;

	}
}


/*
功能说明：
将一个条目插入到形成特征描述符的方向直方图数组中，该条目分配最多8个方向，对于每个维度每个方向乘以（1-d） 的权重，
其中d是以bin测量的bin中心值距离。
//hist为方向直方图的三维数组
//xbin子行
//ybin子列
//obin子方向
//mag权重
//d方向直方图的宽度
//bins每个直方图中bin的个数
三线性插值
*/
void InterpHistEntry(double ***hist, double xbin, double ybin, double obin, double mag, int bins, int d)
{
	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(ybin);  //向下取整
	c0 = cvFloor(xbin);
	o0 = cvFloor(obin);
	d_r = ybin - r0;   //小数余项
	d_c = xbin - c0;
	d_o = obin - o0;

//这里也可以看出前面rbin，cbin为何减0.5，这样原点上这点d_r,d_c均为0.5，即原点上这点的方向将被平均分配叠加在它周围4个直方图上面。

	/*
	做插值：
	xbin,ybin,obin:种子点所在子窗口的位置和方向
	所有种子点都将落在4*4的窗口中
	r0,c0取不大于xbin，ybin的正整数
	r0,c0只能取到0,1,2
	xbin,ybin在(-1, 2)
	r0取不大于xbin的正整数时。
	r0+0 <= xbin <= r0+1
	mag在区间[r0,r1]上做插值
	obin同理
	*/

	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);  //对近邻两行的贡献因子
			row = hist[rb];    //第rb行上的hist
			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];   //h表示第rb行cb列上的hist
					for (o = 0; o <= 1; o++)   //沿着直方图方向不超过1个单位长度
					{
						ob = (o0 + o) % bins;        //bins=8，8个柱子
						v_o = v_c * ((o == 0) ? 1.0 - d_o : d_o);  //对近邻两个方向的贡献因子
						h[ob] += v_o;   
					}
				}
			}
		}
	}


}
/********************************************************************************************************************************
*模块说明:
*        模块七--步骤1:计算描述子的直方图
*功能说明：
gauss：图像指针

x：特征点所在的行

y：特征点所在的列

ori：特征点的主方向

octave_scale：特征点的尺度//特征点所在组的尺度

width：计算方向直方图时，将特征点附近划分为width*width个区域，每个区域生成一个直方图，默认width为4

bins：每个直方图中bins的个数

返回值：double类型的三维数组，即一个d*d的二维数组，数组中每个元素是一个有n个bins的直方图数组
*********************************************************************************************************************************/
double*** CalculateDescrHist(const Mat& gauss, int x, int y, double octave_scale, double ori, int bins, int width)
{
	double ***hist = new double**[width];   //width*width*bins的三维直方图数组 //为第一维分配空间

	for (int i = 0; i < width; i++)
	{
		hist[i] = new double*[width];  //为第二维分配空间
		for (int j = 0; j < width; j++)
		{
			hist[i][j] = new double[bins]; //为第三维分配空间
		}
	}
	
	for (int r = 0; r < width; r++)
		for (int c = 0; c < width; c++)
			for (int o = 0; o < bins; o++)
				hist[r][c][o] = 0.0;

	//为了保证特征描述子具有旋转不变性，要以特征点为中心，在附近邻域内旋转θ角，即旋转为特征点的方向
	double cos_ori = cos(ori);
	double sin_ori = sin(ori);

	//6.1高斯权值，sigma等于描述字窗口宽度的一半
	double sigma = 0.5 * width;
	double conste = -1.0 / (2 * sigma*sigma);

	double PI2 = CV_PI * 2;

	//计算特征描述子过程中，特征点周围的width*width个区域中，每个区域的宽度为3*σ个像素，
	double sub_hist_width = DESCR_SCALE_ADJUST * octave_scale; //DESCR_SCALE_ADJUST=3

	//【1】计算描述子所需的图像领域区域的半径
	//考虑到要进行双线性插值，每个区域的宽度应为:DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )
	//在考虑到旋转因素，每个区域的宽度应为：DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )* sqrt(2)
	//所以搜索的半径是：DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )* sqrt(2) / 2
	int    radius = (sub_hist_width*sqrt(2.0)*(width + 1)) / 2.0 + 0.5;    //[1]0.5取四舍五入
	double grad_ori;
	double grad_mag;

	//遍历每个区域的像素
	for (int i = -radius; i <= radius; i++)
	{
		for (int j = -radius; j <= radius; j++)
		{
			//坐标旋转为主方向,以确保旋转不变性。
			double rot_x = (cos_ori * j - sin_ori * i) / sub_hist_width;
			double rot_y = (sin_ori * j + cos_ori * i) / sub_hist_width;

			//将邻域内的采样点分配到对应的子区域内，将子区域内的梯度值分配到8个方向上，计算其权值
			double xbin = rot_x + width / 2 - 0.5;                         //[2]xbin,ybin为落在4*4窗口中的下标值
			double ybin = rot_y + width / 2 - 0.5;

			if (xbin > -1.0 && xbin < width && ybin > -1.0 && ybin < width)
			{
				if (CalcGradMagOri(gauss, x + j, y + i, grad_mag, grad_ori)) //[3]计算关键点的梯度方向
				{
					grad_ori = (CV_PI - grad_ori) - ori;
					while (grad_ori < 0.0)
						grad_ori += PI2;
					while (grad_ori >= PI2)
						grad_ori -= PI2;              //把梯度方向归入[0，2Π]

					double obin = grad_ori * (bins / PI2);   
					//梯度方向*(8/2Π)归到[0，8]，8为直方图有8个bin，2Π为角度取值范围的长度，把方向的索引归于0~8.
					//conste = -1.0 / (2 * sigma*sigma);
					//计算权重
					double weight = exp(conste*(rot_x*rot_x + rot_y * rot_y));

					InterpHistEntry(hist, xbin, ybin, obin, grad_mag*weight, bins, width); //线性差值求每个种子点八个方向的梯度。

				}
			}
		}
	}

	return hist;
}

//归一化特征点的特征描述子，即将特征描述子数组中每个元素除以特征描述子的模
void NormalizeDescr(Keypoint& feat)
{
	double cur, len_inv, len_sq = 0.0;
	int i, d = feat.descr_length;   //特征描述子的维数
 
//求特征描述子的模
	for (i = 0; i < d; i++)
	{
		cur = feat.descriptor[i];
		len_sq += cur*cur;
	}
	len_inv = 1.0 / sqrt(len_sq);  //sqrt( len_sq )为特征描述子的模
	for (i = 0; i < d; i++)
		feat.descriptor[i] *= len_inv; //特征描述子中每个元素除以特征描述子的模，完成归一化
}
/********************************************************************************************************************************
*模块说明:
*        模块七--步骤2:直方图到描述子的转换
*功能说明：
将某特征点的方向直方图转换为特征描述子向量，对特征描述子归一化并将所有元素转化为整型，存入指定特征点中
参数：
hist：width*width*bins的三维直方图数组
width：计算方向直方图时，将特征点附近划分为width*width个区域，每个区域生成一个直方图
bins：每个直方图的bin个数
feature：特征点指针，将计算好的特征描述子存入其中
*********************************************************************************************************************************/
void HistToDescriptor(double ***hist, int width, int bins, Keypoint& feature)
{
	int int_val, i, r, c, o, k = 0;
//遍历width*width*bins的三维直方图数组，将其中的所有数据(一般是128个)都存入feature结构的descriptor成员中
	for (r = 0; r < width; r++)
		for (c = 0; c < width; c++)
			for (o = 0; o < bins; o++)
			{
				feature.descriptor[k++] = hist[r][c][o];
			}
	
	feature.descr_length = k;  //特征描述子的维数，一般是128
	NormalizeDescr(feature);                           //[1]描述子特征向量归一化
  //归一化特征点的特征描述子，即将特征描述子数组中每个元素除以特征描述子的模

  //	DESCR_MAG_THR 0.2
  //遍历特征描述子向量，将超过阈值DESCR_MAG_THR的元素强行赋值DESCR_MAG_THR
	for (i = 0; i < k; i++)                           //[2]描述子向量门限
		if (feature.descriptor[i] > DESCR_MAG_THR)
			feature.descriptor[i] = DESCR_MAG_THR;

	NormalizeDescr(feature);                           //[3]描述子进行最后一次的归一化操作

//遍历特征描述子向量，每个元素乘以系数INT_DESCR_FCTR( 512.0)来变为整型，并且最大值不能超过255
	for (i = 0; i < k; i++)                           //[4]将单精度浮点型的描述子转换为整形的描述子
	{
		int_val = INT_DESCR_FCTR * feature.descriptor[i];
		feature.descriptor[i] = min(255, int_val);
	}
}
/********************************************************************************************************************************
*模块说明:
*        模块七:6 关键点描述
*功能说明：
*        1）通过以上步骤，对于一个关键点，拥有三个信息：位置、尺度、方向
*        2）接下来就是为每个关键点建立一个描述符，用一组向量来将这个关键点描述出来，使其不随各种变化而变化，比如光照、视角变化等等
*        3）这个描述子不但包括关键点，也包含关键点周围对其贡献的像素点，并且描述符应该有较高的独特性，以便于特征点正确的匹配概率
*        1）SIFT描述子----是关键点邻域高斯图像梯度统计结果的一种表示。
*        2）通过对关键点周围图像区域分块，计算块内梯度直方图，生成具有独特性的向量
*        3）这个向量是该区域图像信息的一种表述和抽象，具有唯一性。
*Lowe论文：
*    Lowe建议描述子使用在关键点尺度空间内4*4的窗口中计算的8个方向的梯度信息，共4*4*8=128维向量来表征。具体的步骤如下所示:
*        1)确定计算描述子所需的图像区域
*        2）将坐标轴旋转为关键点的方向，以确保旋转不变性，如CSDN博文中的图6.2所示；旋转后邻域采样点的新坐标可以通过公式(6-2)计算
*        3）将邻域内的采样点分配到对应的子区域，将子区域内的梯度值分配到8个方向上，计算其权值
*        4）插值计算每个种子点八个方向的梯度
*        5）如上统计的4*4*8=128个梯度信息即为该关键点的特征向量。特征向量形成后，为了去除光照变化的影响，需要对它们进行归一化处理，
*           对于图像灰度值整体漂移，图像各点的梯度是邻域像素相减得到的，所以也能去除。得到的描述子向量为H，归一化之后的向量为L
*        6）描述子向量门限。非线性光照，相机饱和度变化对造成某些方向的梯度值过大，而对方向的影响微弱。因此，设置门限值（向量归一化
*           后，一般取0.2）截断较大的梯度值。然后，在进行一次归一化处理，提高特征的鉴别性。
*        7）按特征点的尺度对特征描述向量进行排序
*        8）至此，SIFT特征描述向量生成。
将关键点附近的区域划分为d*d(Lowe建议d=4)个子区域，每个子区域作为一个种子点，每个种子点有8个方向。
*********************************************************************************************************************************/
void DescriptorRepresentation(vector<Keypoint>& features, const vector<Mat>& gauss_pyr, int bins, int width)   //bins=8,width=4
{
	double ***hist;

	for (int i = 0; i < features.size(); i++)
	{                                                                       //[1]计算描述子的直方图
		//gauss_pyr[features[i].octave*(INTERVALS + 3) + features[i].interval  计算特征点所在矩阵中的位置；
		hist = CalculateDescrHist(gauss_pyr[features[i].octave*(INTERVALS + 3) + features[i].interval],
			features[i].x, features[i].y, features[i].octave_scale, features[i].ori, bins, width);

		HistToDescriptor(hist, width, bins, features[i]);                   //[2]直方图到描述子的转换

/*释放计算特征描述子过程中用到的方向直方图的内存空间
hist：方向直方图的指针，是一个width*width*bins的三维直方图数组
*/
		for (int j = 0; j < width; j++)
		{
			for (int k = 0; k < width; k++)
			{
				delete[] hist[j][k]; //释放第三维的内存
			}
			delete[] hist[j];   //释放第二维的内存
		}
		delete[] hist;  //释放第一维的内存
	}
}


/*比较函数，将特征点按尺度的降序排列
参数：
f1：第一个特征点的指针
f2：第二个特征点的指针
返回值：如果f1的尺度小于f2的尺度，返回1；否则返回-1；若相等返回0
*/
bool FeatureCmp(Keypoint& f1, Keypoint& f2)
{
	return f1.scale < f2.scale;
}
/*******************************************************************************************************************
*函数说明:
*        最大的模块1：SIFT算法模块
*函数参数说明:
*        1---const Mat &src---------------准备进行特征点检测的原始图片
*        2---Vector<Keypoint>& features---用来存储检测出来的关键点
*        3---double sigma-----------------sigma值
*        4---int intervals----------------关键点所在的层数
********************************************************************************************************************/
void Sift(const Mat &src, vector<Keypoint>& features, double sigma, int intervals)
{
	std::cout << "【Step_one】Create -1 octave gaussian pyramid image" << std::endl;

	std::cout << "sigma      = " << sigma << std::endl;
	std::cout << "intervals     = " << intervals << std::endl;

	cv::Mat          init_gray;
	CreateInitSmoothGray(src, init_gray, sigma);   //灰度图
	int octaves = log((double)min(init_gray.rows, init_gray.cols)) / log(2.0) - 2;             //计算高斯金字塔的层数
	std::cout << "【1】The height and width of init_gray_img = " << init_gray.rows << "*" << init_gray.cols << std::endl;
	std::cout << "【2】The octaves of the gauss pyramid      = " << octaves << std::endl;


	std::cout << "【Step_two】Building gaussian pyramid ..." << std::endl;
	std::vector<Mat> gauss_pyr;
	GaussianPyramid(init_gray, gauss_pyr, octaves, intervals, sigma);    //高斯金字塔
	write_pyr(gauss_pyr, "gausspyramid");


	std::cout << "【Step_three】Building difference of gaussian pyramid..." << std::endl;
	vector<Mat> dog_pyr;
	DogPyramid(gauss_pyr, dog_pyr, octaves, intervals);     //差分金字塔
	write_pyr(dog_pyr, "dogpyramid");



	std::cout << "【Step_four】Deatecting local extrema..." << std::endl;
	vector<Keypoint> extrema;
	DetectionLocalExtrema(dog_pyr, extrema, octaves, intervals);   //极值点初探（排除较小点，找到极值点，消除边缘响应）
	std::cout << "【3】keypoints cout: " << extrema.size() << " --" << std::endl;
	std::cout << "【4】extrema detection finished." << std::endl;
	std::cout << "【5】please look dir gausspyramid, dogpyramid and extrema.txt.--" << endl;



	std::cout << "【Step_five】CalculateScale..." << std::endl;
	CalculateScale(extrema, sigma, intervals);   //计算特征点的尺度
	HalfFeatures(extrema);            //图像缩放



	std::cout << "【Step_six】Orientation assignment..." << std::endl;
	OrientationAssignment(extrema, features, gauss_pyr);    //关键点直方图计算，使用的extrema是缩放后的图像
	std::cout << "【6】features count: " << features.size() << std::endl;



	std::cout << "【Step_seven】DescriptorRepresentation..." << std::endl;
	//DESCR_HIST_BINS=8;DESCR_WINDOW_WIDTH=4;
	DescriptorRepresentation(features, gauss_pyr, DESCR_HIST_BINS, DESCR_WINDOW_WIDTH);  //关键点主方向的描述
	sort(features.begin(), features.end(), FeatureCmp);  //排序， FeatureCmp是排序方法
	cout << "finished." << endl;
}
/*******************************************************************************************************************
*函数说明:
*        画出SIFT特征点的具体函数
********************************************************************************************************************/
void DrawSiftFeature(Mat& src, Keypoint& feat, CvScalar color)
{
	int len, hlen, blen, start_x, start_y, end_x, end_y, h1_x, h1_y, h2_x, h2_y;
	double scl, ori;
	double scale = 5.0;
	double hscale = 0.75;
	CvPoint start, end, h1, h2;

	/* compute points for an arrow scaled and rotated by feat's scl and ori */
	start_x = cvRound(feat.dx);
	start_y = cvRound(feat.dy);
	scl = feat.scale;
	ori = feat.ori;
	len = cvRound(scl * scale);
	hlen = cvRound(scl * hscale);
	blen = len - hlen;
	end_x = cvRound(len *  cos(ori)) + start_x;
	end_y = cvRound(len * -sin(ori)) + start_y;
	h1_x = cvRound(blen *  cos(ori + CV_PI / 18.0)) + start_x;
	h1_y = cvRound(blen * -sin(ori + CV_PI / 18.0)) + start_y;
	h2_x = cvRound(blen *  cos(ori - CV_PI / 18.0)) + start_x;
	h2_y = cvRound(blen * -sin(ori - CV_PI / 18.0)) + start_y;
	start = cvPoint(start_x, start_y);
	end = cvPoint(end_x, end_y);
	h1 = cvPoint(h1_x, h1_y);
	h2 = cvPoint(h2_x, h2_y);

	line(src, start, end, color, 1, 8, 0);
	line(src, end, h1, color, 1, 8, 0);
	line(src, end, h2, color, 1, 8, 0);
}
/*******************************************************************************************************************
*函数说明:
*         最大的模块3：画出SIFT特征点
********************************************************************************************************************/
void DrawSiftFeatures(Mat& src, vector<Keypoint>& features)
{
	CvScalar color = CV_RGB(0, 255, 0);
	for (int i = 0; i < features.size(); i++)
	{
		DrawSiftFeature(src, features[i], color);
	}
}
/*******************************************************************************************************************
*函数说明:
*         最大的模块2：画出关键点KeyPoints
********************************************************************************************************************/
void DrawKeyPoints(Mat &src, vector<Keypoint>& keypoints)
{
	int j = 0;
	for (int i = 0; i < keypoints.size(); i++)
	{

		CvScalar color = { 255, 0 ,0 };
		circle(src, Point(keypoints[i].dx, keypoints[i].dy), 3, color);
		j++;
	}
}

const char* GetFileName(const char* dir, int i)
{
	char *name = new char[50];
	sprintf(name, "%s\%d.jpg", dir, i);
	return name;
}

void cv64f_to_cv8U(const Mat& src, Mat& dst)
{
	double* data = (double*)src.data;
	int step = src.step / sizeof(*data);

	if (!dst.empty())
		return;
	dst.create(src.size(), CV_8U);

	uchar* dst_data = dst.data;

	for (int i = 0, m = 0; i < src.cols; i++, m++)
	{
		for (int j = 0, n = 0; j < src.rows; j++, n++)
		{
			*(dst_data + dst.step*j + i) = (uchar)(*(data + step*j + i) * 255);
		}
	}
}


//通过转换后保存的图像，会失真,和imshow显示出的图像相差很大
void writecv64f(const char* filename, const Mat& mat)
{
	Mat dst;
	cv64f_to_cv8U(mat, dst);
	imwrite(filename, dst);
}

void write_pyr(const vector<Mat>& pyr, const char* dir)
{
	for (int i = 0; i < pyr.size(); i++)
	{
		writecv64f(GetFileName(dir, i), pyr[i]);
	}
}

void read_features(vector<Keypoint>&features, const char*file)
{
	ifstream in(file);
	int n = 0, dims = 0;
	in >> n >> dims;
	cout << n << " " << dims << endl;
	for (int i = 0; i < n; i++)
	{
		Keypoint key;
		in >> key.dy >> key.dx >> key.scale >> key.ori;
		for (int j = 0; j < dims; j++)
		{
			in >> key.descriptor[j];
		}
		features.push_back(key);
	}
	in.close();
}
/*******************************************************************************************************************
*函数说明:
*         最大的模块4：将特征点写入文本文件
********************************************************************************************************************/
void write_features(const vector<Keypoint>&features, const char*file)
{
	ofstream dout(file);
	dout << features.size() << " " << FEATURE_ELEMENT_LENGTH << endl;
	for (int i = 0; i < features.size(); i++)
	{
		dout << features[i].dy << " " << features[i].dx << " " << features[i].scale << " " << features[i].ori << endl;
		for (int j = 0; j < FEATURE_ELEMENT_LENGTH; j++)
		{
			if (j % 20 == 0)
				dout << endl;
			dout << features[i].descriptor[j] << " ";
		}
		dout << endl;
	}
	dout.close();
}

　　

sift.h

#ifndef SIFT_H
#define SIFT_H


#include <opencv2/core/core.hpp>
#include "opencv2/highgui/highgui.hpp"



using namespace cv;
using namespace std;

typedef double pixel_t;                             //【1】像素类型

#define INIT_SIGMA 0.5                               //【2】初始sigma
#define SIGMA 1.6
#define INTERVALS 3                                  //【3】高斯金字塔中每组图像中有三层/张图片

#define RATIO 10                                     //【4】半径r
#define MAX_INTERPOLATION_STEPS 5                    //【5】最大空间间隔
#define DXTHRESHOLD 0.03                             //【6】|D(x^)| < 0.03   0.03极值点太多

#define ORI_HIST_BINS 36                             //【7】bings=36
#define ORI_SIGMA_TIMES 1.5
#define ORI_WINDOW_RADIUS 3.0 * ORI_SIGMA_TIMES 
#define ORI_SMOOTH_TIMES 2
#define ORI_PEAK_RATIO 0.8
#define FEATURE_ELEMENT_LENGTH 128
#define DESCR_HIST_BINS 8
#define IMG_BORDER 5 
#define DESCR_WINDOW_WIDTH 4
#define DESCR_SCALE_ADJUST 3
#define DESCR_MAG_THR 0.2
#define INT_DESCR_FCTR 512.0
													/*********************************************************************************************
													*模块说明：
													*        关键点/特征点的结构体声明
													*注意点1：
													*        在高斯金字塔构建的过程中，一幅图像可以产生好几组(octave)图像，一组图像包含几层(inteval)
													*        图像
													*********************************************************************************************/
struct Keypoint
{
	int    octave;                                        //【1】关键点所在组
	int    interval;                                      //【2】关键点所在层
	double offset_interval;                               //【3】调整后的层的增量

	int    x;                                             //【4】x,y坐标,根据octave和interval可取的层内图像
	int    y;
	double scale;                                         //【5】空间尺度坐标scale = sigma0*pow(2.0, o+s/S)

	double dx;                                            //【6】特征点坐标，该坐标被缩放成原图像大小 
	double dy;

	double offset_x;
	double offset_y;

	//============================================================
	//1---高斯金字塔组内各层尺度坐标，不同组的相同层的sigma值相同
	//2---关键点所在组的组内尺度
	//============================================================
	double octave_scale;                                  //【1】offset_i;
	double ori;                                           //【2】方向
	int    descr_length;
	double descriptor[FEATURE_ELEMENT_LENGTH];            //【3】特征点描述符            
	double val;                                           //【4】极值
};
/*********************************************************************************
*模块说明：
*        SIFT算法内,所有成员函数的声明
*********************************************************************************/
void read_features(vector<Keypoint>&features, const char*file);
void write_features(const vector<Keypoint>&features, const char*file);

void testInverse3D();

void write_pyr(const vector<Mat>& pyr, const char* dir);
void DrawKeyPoints(Mat &src, vector<Keypoint>& keypoints);

const char* GetFileName(const char* dir, int i);

void ConvertToGray(const Mat& src, Mat& dst);
void DownSample(const Mat& src, Mat& dst);
void UpSample(const Mat& src, Mat& dst);

void GaussianTemplateSmooth(const Mat &src, Mat &dst, double);
void GaussianSmooth2D(const Mat &src, Mat &dst, double sigma);
void GaussianSmooth(const Mat &src, Mat &dst, double sigma);

void Sift(const Mat &src, vector<Keypoint>& features, double sigma = SIGMA, int intervals = INTERVALS);

void CreateInitSmoothGray(const Mat &src, Mat &dst, double);
void GaussianPyramid(const Mat &src, vector<Mat>&gauss_pyr, int octaves, int intervals, double sigma);

void Sub(const Mat& a, const Mat& b, Mat & c);

void DogPyramid(const vector<Mat>& gauss_pyr, vector<Mat>& dog_pyr, int octaves, int intervals);
void DetectionLocalExtrema(const vector<Mat>& dog_pyr, vector<Keypoint>& extrema, int octaves, int intervals);
void DrawSiftFeatures(Mat& src, vector<Keypoint>& features);

#endif

　　

参考：

链接：https://blog.csdn.net/jkldlnx/article/details/72388669

OpenCV成长之路(9)：特征点检测与图像匹配

最后c代码解释：

https://wenku.baidu.com/view/d7edd2464b73f242336c5ffa.html

https://www.cnblogs.com/Alliswell-WP/tag/SIFT/

重要要求