SURF论文实现与相关数学推导

SURF特征点识别论文实现

SURF特征点识别的主要思路

SURF的特征点检测方法脱胎于DoH特征点检测方法，DoH特征点检测方法计算图像点的Hessian矩阵的行列值来确定该像素点是否为特征点。

H ( L ) = [ L x x L x y L x y L y y ] ( 1 ) H(L)=\begin{bmatrix}Lxx & Lxy \\Lxy & Lyy\end{bmatrix}(1) H(L)=[LxxLxy​LxyLyy​](1)

其中 L x x Lxx Lxx代表高斯正态分布对x的二阶导数与该像素点的卷积值. L x y Lxy Lxy代表高斯正态分布的xy二阶导数与该像素点的卷积值，依次类推 L y y Lyy Lyy代表高斯正态分布对y的二阶导数与该像素点的卷积值。我们知道Hessian矩阵与多元函数的极值点密切相关。

1—如果Hessian矩阵是一个正定矩阵，那么该点是该多元函数的局部极小值点。

2—如果Hessina矩阵是一个负定矩阵，那么该点是该多元函数的局部极大值点。

3—如果Hessian矩阵是一个不定矩阵，那么该点不能确定是否为函数的极值点。

当然我们还要清楚一个事实一个函数与高斯正态分布的导数求卷积相当于对该函数直接求导。例如:

所以上述的Hessian矩阵 L x x , L y y , L y y Lxx,Lyy,Lyy Lxx,Lyy,Lyy由对应的高斯正态分布的导数与图像的卷积求得与直接对图像点进行求导的意义是一样的，但是SURF方法运用了积分图像与盒子滤波器可以简化运算，盒子滤波器是一种对于高斯正态分布导数的近似。

对于函数：

其对x的二阶偏导数：

用matlab画出二阶偏导数的图像：

接下来从Z轴方向观察上面的图形：

可以发现这就是盒子滤波器的精确版本，但是为了与积分图像结合，论文中对上述图像进行了简化.同样我们可以得 L y y , L x y Lyy,Lxy Lyy,Lxy的图像与其盒子滤波器的近似

L y y Lyy Lyy

L x y Lxy Lxy

从Z轴观察上述图像分别为：

其对应的盒子滤波器的近似图像我们称之为 D x x , D y y , D x y Dxx,Dyy,Dxy Dxx,Dyy,Dxy图像如下图所示：

对于公式(1)我们可以得到如下的数学推导：

接下来我们要了解一下积分图像：

可以看到每一点(i,j)的值等于该点左上角所以像素点灰度值的和，当然积分图像是在灰度图像上进行计算的,因为灰度图是单通道的。我们在计算图像卷积的过程中结合盒子滤波器可以快速的进行卷积运算。

这一部分的代码实现如下所示

struct SURFHaar {
	int p0, p1, p2, p3;
	float w;
	SURFHaar() :p0(0), p1(0), p2(0), p3(0), w(0) {};
};//定义盒子滤波器的四个顶点
float calHaarResponse(const int* origin, const SURFHaar* f, int n) {
	double d = 0;
	for (int k = 0; k < n; ++k)
		d += (origin[(f[k].p0)] - origin[f[k].p1] - origin[f[k].p2] + origin[f[k].p3])*f[k].w;
	return (float)d;
}//计算盒子卷积响应值
static void resizeHaarSizes(const int src[][5], SURFHaar* dst, int n, int oldSize, int newSize, int widthStep)
{
	float ratio = (float)newSize / oldSize;
	for (int k = 0; k < n; k++)
	{
		int dx1 = cvRound(ratio*src[k][1]);
		int dy1 = cvRound(ratio*src[k][0]);
		int dx2 = cvRound(ratio*src[k][2]);
		int dy2 = cvRound(ratio*src[k][3]);
		dst[k].p0 = dy1 * widthStep + dx1;
		dst[k].p1 = dy2 * widthStep + dx1;
		dst[k].p2 = dy1 * widthStep + dx2;
		dst[k].p3 = dy2 * widthStep + dx2;
		dst[k].w = src[k][4]/((float)(dx2-dx1)*(dy2-dy1)) ;
		//dst[k].w = src[k][5];
	}
}//对滤波器进行初始化
static void calcLayerDetAndTrace(const Mat& sum, int size, int sampleStep, Mat& det, Mat& trace) {
	const int NX = 3, NY = 3, NXY = 4;
	const int dx_s[NX][5] = { {0, 2, 3, 7, 1}, {3, 2, 6, 7, -2}, {6, 2, 9, 7, 1} };//(0,2)(3,7),1,(3,2)(6,7),-2,(6,2),(9,7),1
	const int dy_s[NY][5] = { {2, 0, 7, 3, 1}, {2, 3, 7, 6, -2}, {2, 6, 7, 9, 1} };//(2,0)(7,3),1,(2,3)(7,6),-2,(2,6),(7,9),1
	const int dxy_s[NXY][5] = { {1, 1, 4, 4, 1}, {5, 1, 8, 4, -1}, {1, 5, 4, 8, -1}, {5, 5, 8, 8, 1} };//(1,1),(4,4),1,(5,1),(8,4),-1,(5,5),(8,8),1
	SURFHaar Dx[NX], DY[NY], DZ[NXY];
	
	if (size > sum.cols - 1 || size > sum.rows - 1)
		return;
	resizeHaarSizes(dx_s, Dx, NX, 9, size, sum.cols);
	resizeHaarSizes(dy_s, DY, NY, 9, size, sum.cols);
	resizeHaarSizes(dxy_s, DZ, NXY, 9, size, sum.cols);
	int sample_r = 1 + (sum.rows - size - 1) / sampleStep;
	int sample_c = 1 + (sum.cols - size - 1) / sampleStep;
	//int sample_r = 1 + (sum.rows - 1) / sampleStep;
	//int sample_c = 1 + (sum.cols -  1) / sampleStep;
	int margin = (size / 2) / sampleStep;
	for (int r = 0; r < sample_r; ++r) {
		const int* sum_ptr = sum.ptr<int>(r*sampleStep);
		float * det_ptr = &det.at<float>(r + margin, margin);
		float * trace_ptr = &trace.at<float>(r + margin, margin);
		for (int c = 0; c < sample_c; ++c) {
			float dx = calHaarResponse(sum_ptr, Dx, 3);
			float dy = calHaarResponse(sum_ptr, DY, 3);
			float dxy = calHaarResponse(sum_ptr, DZ, 4);
			sum_ptr += sampleStep;
			det_ptr[c] = dx * dy - 0.81f*dxy*dxy;
			trace_ptr[c] = dx + dy;
		}
	}
}//计算Det(H)
           

通常情况下要想获取不同尺寸的斑点，必须建立图像的尺度空间金字塔。一般的方法是通过采用不同σ的高斯函数，对图像进行平滑滤波，然后重采样图像以获得更高一层的金字塔图像。Lowe在SIFT方法中就是通过相邻两层图像金字塔相减得到DoG图像，然后再在DoG图像上进行斑点和边缘性检测工作。

由于SURF方法采用了盒子滤波器和积分图像，所以，并不需要像SIFT算法那样去直接建立图像金字塔，而是采用不断增大盒子滤波器模板尺寸的间接方法。通过不同尺寸盒子滤波模板与积分图像求取Hessian矩阵行列式的响应值，然后在响应图像上采用3D非最大值抑制，求取各种不同层次的斑点。

在论文中，首先使用9*9的模板对图像进行滤波，其结果作为最初始的尺寸空间层（此时，尺度值σ=1.2的高斯微分），后续的层将通过逐步放大滤波模板尺寸，以及放大后的模板不断与图像进行滤波得到。由于采用盒子滤波和积分图像，滤波过程并不随着滤波模板尺寸的增加而使运算工作量增加。

与SIFT算法类似，我们需要将尺度空间划分成若干组(Octaves).一个组代表了逐步放大的滤波模板对同一输入图像进行滤波的一系列响应图。每个组又由若干固定的层组成。由于积分图像离散化的原因，两个相邻层之间的最小尺度变化量是由高斯二阶微分滤波器的在微分方向上对正负斑点响应长度Lo决定的。

他是盒子滤波模板尺寸的1/3.对于99的盒子滤波模板Lo为3.下一层的响应长度至少应该在Lo的基础上增加两个像元，以保证一边一个像元，即Lo=5，这样，模板的尺寸就为1515，依次类推，我们得到一个尺寸逐渐增大的模板序列，他们尺寸分别为：99，1515，2121，2727，黑色，白色的区域值加偶数个像元，以保证一个中心像元的存在。

当然这个过程在我们上述的resizeHaarSizes函数中有其实现过程。利用这个模板序列，就可以构建尺度空间了，首先使用99的模板，然后依次使用1515，2121，2727的模板。这样便产生了一组由Hessian矩阵行列式值组成的响应图像。利用这组响应图像就可以应用3D非最大值抑制搜素技术，在位置和尺度上进行斑点搜索，产生的一系列盒子滤波模板尺寸如下图所示：

关于图像的尺度问题，我有几点自己的看法。对于同一个景物，虽然由于成像传感器观察距离远近的不同，焦距的不同，使得成像分辨率不同，因此，特征尺度也不同。但是特征尺度值所对应的区域则是相同的，如下图所示：

两个蓝色圆圈代表了同一特征的不同尺度，也就是说两个圆圈代表同一个物体(^^格瓦斯饮料瓶)，但由于一个距离手机相机距离近，一个距离手机镜头远，所以特征大小不一，我们需要对同一特征的不同大小进行处理，具体就是选取一系列不同的尺度的高斯卷积函数对图像进行滤波处理得到滤波图像，然后进行3D非极大值抑制得到相对准确的响应值。

static void findMaxLayer(const Mat& sum,
	const vector<Mat>& dets, const vector<Mat>& traces,
	const vector<int>& sizes, vector<KeyPoint>& keypoints,
	int octave, int layer, float hessianThreshold, int sampleStep) {
	const int nm = 1;
	//const int DM[nm][5] = { {0,0,9,9,1} };
	//SURFHaar DM;
	int size = sizes[layer];
	int layer_rows = (sum.rows - -1) / sampleStep;
	int layer_cols = (sum.cols - 1) / sampleStep;
	int margin = (sizes[layer + 1] / 2) / sampleStep;
	int step = (int)(dets[layer].step / dets[layer].elemSize());//每一行的有多少元素
	for (int i = margin; i < layer_rows - margin; ++i) {
		const float* det_ptr = dets[layer].ptr<float>(i);
		const float* trace_ptr = traces[layer].ptr<float>(i);
		for (int j = margin; j < layer_cols - margin; ++j) {
			float val0 = det_ptr[j];
			if (val0 > hessianThreshold) {
				int sum_i = (i - (size / 2) / sampleStep)*sampleStep;
				int sum_j = (j - (size / 2 / sampleStep))*sampleStep;
				const float* det1 = &dets[layer - 1].at<float>(i, j);
				const float* det2 = &dets[layer ].at<float>(i, j);
				const float* det3 = &dets[layer + 1].at<float>(i, j);
				float N9[3][9] = { {det1[-step - 1],det1[-step],det1[-step + 1],det1[-1],det1[0],det1[1],det1[step - 1],det1[step],det1[step + 1]},
				{det2[-step - 1],det2[-step],det2[-step + 1],det2[-1],det2[0],det2[1],det2[step - 1],det2[step],det2[step + 1]},
				{det3[-step - 1],det3[-step],det3[-step + 1],det3[-1],det3[0],det3[1],det3[step - 1],det3[step],det3[step + 1]} };
				if (val0 > N9[0][0] && val0 > N9[0][1] && val0 > N9[0][2] &&
					val0 > N9[0][3] && val0 > N9[0][4] && val0 > N9[0][5] &&
					val0 > N9[0][6] && val0 > N9[0][7] && val0 > N9[0][8] &&
					val0 > N9[1][0] && val0 > N9[1][1] && val0 > N9[1][2] &&
					val0 > N9[1][3] && val0 > N9[1][5] &&
					val0 > N9[1][6] && val0 > N9[1][7] && val0 > N9[1][8] &&
					val0 > N9[2][0] && val0 > N9[2][1] && val0 > N9[2][2] &&
					val0 > N9[2][3] && val0 > N9[2][4] && val0 > N9[2][5] &&
					val0 > N9[2][6] && val0 > N9[2][7] && val0 > N9[2][8])
				{
					float center_i = sum_i + (size - 1)*0.5f;
					float center_j = sum_j + (size - 1)*0.5f;
					KeyPoint kpt(center_j, center_i,(float)sizes[layer], -1, val0, octave, CV_SIGN(trace_ptr[j]));//#define  CV_CMP(a,b) (((a) > (b)) - ((a) < (b)))
					int ds = size - sizes[layer - 1];//#define  CV_SIGN(a)  CV_CMP((a),0)
					int inter = interpolateKeypoint(N9, sampleStep, sampleStep, ds, kpt);
					if (inter) {
						keypoints.push_back(kpt);
					}
				}

			}
		}
	}
}//3D非极大值抑制
           

但是我们应该清楚我们的上述的操作最小单元是一个像素点，如果这个响应最大值实际上在两个像素点之间，那么我们的结果很显然就是不准确的，事实上在大多数的情况下，最大响应值对应在图像坐标上的位置都不是整数，那么就需要相应的插值运算。具体计算步骤如下所示：

具体的实现代码如下所示：

static int interpolateKeypoint(float N9[3][9], int dx, int dy, int ds, KeyPoint& kpt) {
	Vec3f b(-(N9[1][5] - N9[1][3]) / 2, -(N9[1][7] - N9[1][1]) / 2, -(N9[2][4] - N9[0][4]) / 2);
	Matx33f A((N9[1][5] + N9[1][3] - 2 * N9[1][4]), (N9[1][8] - N9[1][6] - N9[1][2] + N9[1][0]) / 4, (N9[2][5] - N9[2][3] - N9[0][5] + N9[0][3]) / 4,
		(N9[1][8] - N9[1][2] - N9[1][6] + N9[1][0]) / 4, (N9[1][7] + N9[1][1] - 2 * N9[1][4]), (N9[2][7] - N9[2][1] - N9[0][7] + N9[0][1]) / 4,
		(N9[2][5] - N9[2][3] - N9[0][5] + N9[0][3]) / 4, (N9[2][7] - N9[2][1] - N9[0][7] + N9[0][1]) / 4, (N9[0][4] + N9[2][4] - 2 * N9[1][4]));
	Vec3f x = A.solve(b, DECOMP_LU);
	bool ok = (x[0] != 0 || x[1] != 0 || x[2] != 0) && std::abs(x[0]) <= 1 && std::abs(x[1]) <= 1 && std::abs(x[2]) <=1;
	if (ok) {
		kpt.pt.x += 0.5*dx * x[0];
		kpt.pt.y += 0.5*dy * x[1];
		kpt.size = (float)cvRound(kpt.size + 0.5*ds * x[2]);
	}
	return ok;
}
           

到这里我们特征点识别的主要步骤就很清晰了，首先根据基于积分图像与盒子滤波器的方法求出Hessian矩阵的近似值，为了构建图像的尺度空间，我们选取一系列不同的尺度值对图像进行滤波处理，尺度值不同则盒子滤波器的大小也不同，将这些不同大小的盒子滤波器处理后的图像分成不同的组，每一组有不小于三层(只有大于三层才能进行3D非极大值抑制)。但是这样寻找的极值点不是最优点，我们需要在亚像素级别对求得的极大值进行插值运算，将以上代码整合一起的函数如下所示：

static void caltotalHessian(const Mat& sum, vector<KeyPoint>&keypoints, int nOctave, int octaveLayer, float hessianThreshold) {
	const int sample_step = 1;
	int nTotalayer = nOctave * (octaveLayer + 2);
	int nMiddleLayer = nOctave * octaveLayer;
	vector<Mat>ndets(nTotalayer);
	vector<Mat>ntraces(nTotalayer);
	vector<int>nsizes(nTotalayer);
	vector<int>nsamplesteps(nTotalayer);
	vector<int>nmiddleindices(nMiddleLayer);
	keypoints.clear();
	int index = 0, middleIndex = 0, step = sample_step;
	for (int o = 0; o < nOctave; ++o) {
		for (int l = 0; l < octaveLayer + 2; ++l) {
			ndets[index].create((sum.rows - 1) / step, (sum.cols - 1) / step, CV_32FC1);
			ntraces[index].create((sum.rows - 1) / step, (sum.cols - 1) / step, CV_32FC1);
			nsizes[index] = (SURF_HAAR_SIZE0 + SURF_HAAR_SIZE_INC * l) << o;
			nsamplesteps[index] = step;
			if (0 < l && l <= octaveLayer)
				nmiddleindices[middleIndex++] = index;
			index++;
			
		}
		step *= 2;
	}
	for (int i = 0; i < nTotalayer; ++i) {
	
		calcLayerDetAndTrace(sum, nsizes[i], nsamplesteps[i], ndets[i], ntraces[i]);
	}
	for(int j = 0; j < middleIndex; ++j ){
		int layer = nmiddleindices[j];
	
		int octave = j / octaveLayer;
		int s_step = nsamplesteps[layer];
		
		findMaxLayer(sum, ndets, ntraces,nsizes, keypoints, octave, layer, hessianThreshold, s_step);
	}
	
	
}
           

接下来我们测试一下代码：

#include "pch.h"
#include<stdio.h>
#include<stdlib.h>
#include <iostream>
#include<opencv2/opencv.hpp>
#include<opencv2/highgui/highgui.hpp>

using namespace std;
using namespace cv;
static const int SURF_HAAR_SIZE0 = 9;
static const int SURF_HAAR_SIZE_INC = 6;

struct SURFHaar {
	int p0, p1, p2, p3;
	float w;
	SURFHaar() :p0(0), p1(0), p2(0), p3(0), w(0) {};
};
float calHaarResponse(const int* origin, const SURFHaar* f, int n) {
	double d = 0;
	for (int k = 0; k < n; ++k)
		d += (origin[(f[k].p0)] - origin[f[k].p1] - origin[f[k].p2] + origin[f[k].p3])*f[k].w;
	return (float)d;
}

static void resizeHaarSizes(const int src[][5], SURFHaar* dst, int n, int oldSize, int newSize, int widthStep)
{
	float ratio = (float)newSize / oldSize;
	for (int k = 0; k < n; k++)
	{
		int dx1 = cvRound(ratio*src[k][1]);
		int dy1 = cvRound(ratio*src[k][0]);
		int dx2 = cvRound(ratio*src[k][3]);
		int dy2 = cvRound(ratio*src[k][2]);
		dst[k].p0 = dy1 * widthStep + dx1;
		dst[k].p1 = dy2 * widthStep + dx1;
		dst[k].p2 = dy1 * widthStep + dx2;
		dst[k].p3 = dy2 * widthStep + dx2;
		dst[k].w = src[k][4]/((float)(dx2-dx1)*(dy2-dy1)) ;
		//dst[k].w = src[k][4];
	}
}
static void calcLayerDetAndTrace(const Mat& sum, int size, int sampleStep, Mat& det, Mat& trace) {
	const int NX = 3, NY = 3, NXY = 4;
	const int dx_s[NX][5] = { {0, 2, 3, 7, 1}, {3, 2, 6, 7, -2}, {6, 2, 9, 7, 1} };//(0,2)(3,7),1,(3,2)(6,7),-2,(6,2),(9,7),1
	const int dy_s[NY][5] = { {2, 0, 7, 3, 1}, {2, 3, 7, 6, -2}, {2, 6, 7, 9, 1} };//(2,0)(7,3),1,(2,3)(7,6),-2,(2,6),(7,9),1
	const int dxy_s[NXY][5] = { {1, 1, 4, 4, 1}, {5, 1, 8, 4, -1}, {1, 5, 4, 8, -1}, {5, 5, 8, 8, 1} };//(1,1),(4,4),1,(5,1),(8,4),-1,(5,5),(8,8),1
	SURFHaar Dx[NX], DY[NY], DZ[NXY];
	
	if (size > sum.cols - 1 || size > sum.rows - 1)
		return;
	resizeHaarSizes(dx_s, Dx, NX, 9, size, sum.cols);
	resizeHaarSizes(dy_s, DY, NY, 9, size, sum.cols);
	resizeHaarSizes(dxy_s, DZ, NXY, 9, size, sum.cols);
	int sample_r = 1 + (sum.rows - size - 1) / sampleStep;
	int sample_c = 1 + (sum.cols - size - 1) / sampleStep;
	//int sample_r = 1 + (sum.rows - 1) / sampleStep;
	//int sample_c = 1 + (sum.cols -  1) / sampleStep;
	int margin = (size / 2) / sampleStep;
	for (int r = 0; r < sample_r; ++r) {
		const int* sum_ptr = sum.ptr<int>(r*sampleStep);
		float * det_ptr = &det.at<float>(r + margin, margin);
		float * trace_ptr = &trace.at<float>(r + margin, margin);
		for (int c = 0; c < sample_c; ++c) {
			float dx = calHaarResponse(sum_ptr, Dx, 3);
			float dy = calHaarResponse(sum_ptr, DY, 3);
			float dxy = calHaarResponse(sum_ptr, DZ, 4);
			sum_ptr += sampleStep;
			det_ptr[c] = dx * dy - 0.81f*dxy*dxy;
			trace_ptr[c] = dx + dy;
		}
	}
}
static int interpolateKeypoint(float N9[3][9], int dx, int dy, int ds, KeyPoint& kpt) {
	Vec3f b(-(N9[1][5] - N9[1][3]) / 2, -(N9[1][7] - N9[1][1]) / 2, -(N9[2][4] - N9[0][4]) / 2);
	Matx33f A((N9[1][5] + N9[1][3] - 2 * N9[1][4]), (N9[1][8] - N9[1][6] - N9[1][2] + N9[1][0]) / 4, (N9[2][5] - N9[2][3] - N9[0][5] + N9[0][3]) / 4,
		(N9[1][8] - N9[1][2] - N9[1][6] + N9[1][0]) / 4, (N9[1][7] + N9[1][1] - 2 * N9[1][4]), (N9[2][7] - N9[2][1] - N9[0][7] + N9[0][1]) / 4,
		(N9[2][5] - N9[2][3] - N9[0][5] + N9[0][3]) / 4, (N9[2][7] - N9[2][1] - N9[0][7] + N9[0][1]) / 4, (N9[0][4] + N9[2][4] - 2 * N9[1][4]));
	Vec3f x = A.solve(b, DECOMP_LU);
	bool ok = (x[0] != 0 || x[1] != 0 || x[2] != 0) && std::abs(x[0]) <= 1 && std::abs(x[1]) <= 1 && std::abs(x[2]) <=1;
	if (ok) {
		kpt.pt.x += 0.5*dx * x[0];
		kpt.pt.y += 0.5*dy * x[1];
		kpt.size = (float)cvRound(kpt.size + 0.5*ds * x[2]);
	}
	return ok;
}
static void findMaxLayer(const Mat& sum,
	const vector<Mat>& dets, const vector<Mat>& traces,
	const vector<int>& sizes, vector<KeyPoint>& keypoints,
	int octave, int layer, float hessianThreshold, int sampleStep) {
	const int nm = 1;
	//const int DM[nm][5] = { {0,0,9,9,1} };
	//SURFHaar DM;
	int size = sizes[layer];
	int layer_rows = (sum.rows - -1) / sampleStep;
	int layer_cols = (sum.cols - 1) / sampleStep;
	int margin = (sizes[layer + 1] / 2) / sampleStep;
	int step = (int)(dets[layer].step / dets[layer].elemSize());//每一行的有多少元素
	for (int i = margin; i < layer_rows - margin; ++i) {
		const float* det_ptr = dets[layer].ptr<float>(i);
		const float* trace_ptr = traces[layer].ptr<float>(i);
		for (int j = margin; j < layer_cols - margin; ++j) {
			float val0 = det_ptr[j];
			if (val0 > hessianThreshold) {
				int sum_i = (i - (size / 2) / sampleStep)*sampleStep;
				int sum_j = (j - (size / 2 / sampleStep))*sampleStep;
				const float* det1 = &dets[layer - 1].at<float>(i, j);
				const float* det2 = &dets[layer ].at<float>(i, j);
				const float* det3 = &dets[layer + 1].at<float>(i, j);
				float N9[3][9] = { {det1[-step - 1],det1[-step],det1[-step + 1],det1[-1],det1[0],det1[1],det1[step - 1],det1[step],det1[step + 1]},
				{det2[-step - 1],det2[-step],det2[-step + 1],det2[-1],det2[0],det2[1],det2[step - 1],det2[step],det2[step + 1]},
				{det3[-step - 1],det3[-step],det3[-step + 1],det3[-1],det3[0],det3[1],det3[step - 1],det3[step],det3[step + 1]} };
				if (val0 > N9[0][0] && val0 > N9[0][1] && val0 > N9[0][2] &&
					val0 > N9[0][3] && val0 > N9[0][4] && val0 > N9[0][5] &&
					val0 > N9[0][6] && val0 > N9[0][7] && val0 > N9[0][8] &&
					val0 > N9[1][0] && val0 > N9[1][1] && val0 > N9[1][2] &&
					val0 > N9[1][3] && val0 > N9[1][5] &&
					val0 > N9[1][6] && val0 > N9[1][7] && val0 > N9[1][8] &&
					val0 > N9[2][0] && val0 > N9[2][1] && val0 > N9[2][2] &&
					val0 > N9[2][3] && val0 > N9[2][4] && val0 > N9[2][5] &&
					val0 > N9[2][6] && val0 > N9[2][7] && val0 > N9[2][8])
				{
					float center_i = sum_i + (size - 1)*0.5f;
					float center_j = sum_j + (size - 1)*0.5f;
					KeyPoint kpt(center_j, center_i,(float)sizes[layer], -1, val0, octave, CV_SIGN(trace_ptr[j]));//#define  CV_CMP(a,b) (((a) > (b)) - ((a) < (b)))
					int ds = size - sizes[layer - 1];//#define  CV_SIGN(a)  CV_CMP((a),0)
					int inter = interpolateKeypoint(N9, sampleStep, sampleStep, ds, kpt);
					if (inter) {
						keypoints.push_back(kpt);
					}
				}

			}
		}
	}
}

static void caltotalHessian(const Mat& sum, vector<KeyPoint>&keypoints, int nOctave, int octaveLayer, float hessianThreshold) {
	const int sample_step = 1;
	int nTotalayer = nOctave * (octaveLayer + 2);
	int nMiddleLayer = nOctave * octaveLayer;
	vector<Mat>ndets(nTotalayer);
	vector<Mat>ntraces(nTotalayer);
	vector<int>nsizes(nTotalayer);
	vector<int>nsamplesteps(nTotalayer);
	vector<int>nmiddleindices(nMiddleLayer);
	keypoints.clear();
	int index = 0, middleIndex = 0, step = sample_step;
	for (int o = 0; o < nOctave; ++o) {
		for (int l = 0; l < octaveLayer + 2; ++l) {
			ndets[index].create((sum.rows - 1) / step, (sum.cols - 1) / step, CV_32FC1);
			ntraces[index].create((sum.rows - 1) / step, (sum.cols - 1) / step, CV_32FC1);
			nsizes[index] = (SURF_HAAR_SIZE0 + SURF_HAAR_SIZE_INC * l) << o;
			nsamplesteps[index] = step;
			if (0 < l && l <= octaveLayer)
				nmiddleindices[middleIndex++] = index;
			index++;
			
		}
		step *= 2;
	}
	for (int i = 0; i < nTotalayer; ++i) {
	
		calcLayerDetAndTrace(sum, nsizes[i], nsamplesteps[i], ndets[i], ntraces[i]);
	}
	for(int j = 0; j < middleIndex; ++j ){
		int layer = nmiddleindices[j];
	
		int octave = j / octaveLayer;
		int s_step = nsamplesteps[layer];
		
		findMaxLayer(sum, ndets, ntraces,nsizes, keypoints, octave, layer, hessianThreshold, s_step);
	}
	
	
}
struct KeypointGreater
{
	inline bool operator()(const KeyPoint& kp1, const KeyPoint& kp2) const
	{
		if (kp1.response > kp2.response) return true;
		if (kp1.response < kp2.response) return false;
		if (kp1.size > kp2.size) return true;
		if (kp1.size < kp2.size) return false;
		if (kp1.octave > kp2.octave) return true;
		if (kp1.octave < kp2.octave) return false;
		if (kp1.pt.y < kp2.pt.y) return false;
		if (kp1.pt.y > kp2.pt.y) return true;
		return kp1.pt.x < kp2.pt.x;
	}
};
int main()
{
	
	Mat src1 = imread("C:\\Users\\马良\\Desktop\\7.jpg");

	Mat sum1,src1Gray；
	vector<KeyPoint>keyPoints1；
	cvtColor(src1, src1Gray, CV_RGB2GRAY);

	integral(src1Gray,sum1,CV_32S);

	int nOctave = 3;
	int octaveLayer = 3;
	float hessianThreshold =5000;
	caltotalHessian(sum1,keyPoints1, nOctave, octaveLayer, hessianThreshold);
	drawKeypoints(src1,keyPoints1, dst);
	imshow("src1", src1);
	imshow("dst", dst);
	imwrite("C:\\Users\\马良\\Desktop\\dst.jpg", dst);
	waitKey(0);
	return 0；
}
           

结果如下所示：

可以看到，识别到的特征点多是一些边缘点和亮度较高的点，所以SURF方法效果还是非常不错。的在我下一篇博客中将介绍SURF描述子的生成方法