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描述子的生成方法