SURF特征点检测算法实现代码源自于Opencv2.4.11版本sources\modules\nonfree\src\surf.cpp文件中
这篇文章只会对代码进行介绍（代码的介绍以注释的形式给出），SURF的具体原理可以参考这个——SURF算法

1 SURF特征检测器的构造函数

SURF特征检测器的构造函数有两个，一个是不带参数的，一个是带参数的：

// 默认构造函数(不带参数的构造函数)
SURF::SURF()
{hessianThreshold = 100;               // 默认为100，仅当hessian矩阵的行列式值大于该阈值时，才会被初步认定为特征点extended = false;                    // 默认为false，是否对特征描述符进行扩展，true代表128维，false代表64维upright = false;                    // 默认为false，是否计算特征点得方向，false代表计算，true代表不计算nOctaves = 4;                       // 默认为4，高斯金字塔组数nOctaveLayers = 3;                   // 默认为3，每组高斯金字塔包含的图像层数
}// 带参数的构造函数
SURF::SURF(double _threshold, int _nOctaves, int _nOctaveLayers, bool _extended, bool _upright)
{hessianThreshold = _threshold;        // 必须由用户传入，仅当hessian矩阵的行列式值大于该阈值时，才会被初步认定为特征点extended = _extended;               // 默认为true，是否对特征描述符进行扩展，true代表128维，false代表64维upright = _upright;                  // 默认为false，是否计算特征点得方向，false代表计算，true代表不计算nOctaves = _nOctaves;               // 默认为4，高斯金字塔组数nOctaveLayers = _nOctaveLayers;      // 默认为2，每组高斯金字塔包含的图像层数
}

2 特征点检测和描述

对原始图像进行SURF特征点检测，得到特征点和特征点描述符。

/** 函数：SURF::operator()* 功能：对图像_img进行SURF特征点检测，特征点检测结果存放在keypoints中，特征点描述符存放在_descriptors中* 参数：_img                   输入  待检测特征点的图像*      mask                    输入  掩码图像*      keypoints                输出  检测到的特征点*        _descriptors            输出  特征点描述符*         useProvidedKeypoints    输入  默认为false(没有找到接口可以改)* 返回值：无*/
void SURF::operator()(InputArray _img, InputArray _mask,CV_OUT vector<KeyPoint>& keypoints,OutputArray _descriptors,bool useProvidedKeypoints) const
{// mask1为二值化后的掩码图像，sum为原始图像的积分图像，msum为二值化后掩码图像的积分图像Mat img = _img.getMat(), mask = _mask.getMat(), mask1, sum, msum;// 判断_descriptors容器是否有效。若有效，则计算特征点的描述符，否则不计算bool doDescriptors = _descriptors.needed();// 判断传入图像是否为空，图像深度是否为8位CV_Assert(!img.empty() && img.depth() == CV_8U);// 判断图像是否为单通道(即灰度图或二值图)图像。若不是，则对图像进行转换if( img.channels() > 1 )cvtColor(img, img, COLOR_BGR2GRAY);// 掩码图像如果不为空，即传入了掩码图像，那么它必须为深度为8位的图像，且大小必须与待检测特征点的图像大小相同CV_Assert(mask.empty() || (mask.type() == CV_8U && mask.size() == img.size()));// hessian阈值必须大于0CV_Assert(hessianThreshold >= 0);CV_Assert(nOctaves > 0);CV_Assert(nOctaveLayers > 0);// 求积分图像integral(img, sum, CV_32S);if( !useProvidedKeypoints ){// 如果掩码图像不为空if( !mask.empty() ){// 将掩码图像的变为二值图像(图像数据由0/1组成)cv::min(mask, 1, mask1);// 求掩码图像的积分图像integral(mask1, msum, CV_32S);}// 利用积分图像进行SURF特征点检测，得到特征点检测结果，详细代码注释见2.1fastHessianDetector( sum, msum, keypoints, nOctaves, nOctaveLayers, (float)hessianThreshold );}// 获取检测出来的特征点个数int i, j, N = (int)keypoints.size();// 如果特征点个数大于0if( N > 0 ){Mat descriptors;bool _1d = false;// 是否扩展特征描述符的维数，64维/128维int dcols = extended ? 128 : 64;size_t dsize = dcols*sizeof(float);if( doDescriptors ){// 判断特征描述符向量是否为32位浮点型向量_1d = _descriptors.kind() == _InputArray::STD_VECTOR && _descriptors.type() == CV_32F;if( _1d ){// 分配N*dcols行，1列的32位浮点型矩阵空间_descriptors.create(N*dcols, 1, CV_32F);descriptors = _descriptors.getMat().reshape(1, N);}else{// 分配N行，dcols列的32位浮点型矩阵空间_descriptors.create(N, dcols, CV_32F);descriptors = _descriptors.getMat();}}// 不论用户是否需要描述符，都会调用SURFInvoker，因为它会计算每个特征点的方向，详细代码注释见2.2parallel_for_(Range(0, N), SURFInvoker(img, sum, keypoints, descriptors, extended, upright) );// 移除被标记待删除的特征点，即移除size属性小于等于0的特征点for( i = j = 0; i < N; i++ ){if( keypoints[i].size > 0 ){if( i > j ){keypoints[j] = keypoints[i];if( doDescriptors )memcpy( descriptors.ptr(j), descriptors.ptr(i), dsize);}j++;}}// 删除特征点后，重新调整特征点向量的大小，并将之前计算出来的描述符结果放入_descriptors向量中输出if( N > j ){N = j;keypoints.resize(N);if( doDescriptors ){Mat d = descriptors.rowRange(0, N);if( _1d )d = d.reshape(1, N*dcols);d.copyTo(_descriptors);}}}
}

2.1 特征点检测（fastHessianDetector函数注释）

利用积分图像进行SURF特征点检测，得到特征点检测结果。

/** 函数：fastHessianDetector* 功能：利用积分图像进行SURF特征点检测，特征点检测结果存放在keypoints中* 参数：sum                    输入  待检测特征点图像的积分图像*      mask_sum                输入  掩码图像的积分图像*      keypoints               输出  检测到的特征点*        nOctaves                输入  图像组数*       nOctaveLayers           输入  每一组图像的层数*       hessianThreshold        输入  hessian阈值* 返回值：无*/
static void fastHessianDetector( const Mat& sum, const Mat& mask_sum, vector<KeyPoint>& keypoints,int nOctaves, int nOctaveLayers, float hessianThreshold )
{/* 按照采样步长，沿第一组(octave)图像的x轴和y轴进行采样。每增加一组，采样步长翻倍。注意：增加采样步长能够提高处理速度，但是特征点提取的可靠性也会下降*/const int SAMPLE_STEP0 = 1;// 图像的全部层数 = 图像组数 x (每组图像层数 + 2)int nTotalLayers = (nOctaveLayers+2)*nOctaves;// 图像的中间层数 = 每组图像层数 x 图像组数int nMiddleLayers = nOctaveLayers*nOctaves;vector<Mat> dets(nTotalLayers);vector<Mat> traces(nTotalLayers);vector<int> sizes(nTotalLayers);vector<int> sampleSteps(nTotalLayers);vector<int> middleIndices(nMiddleLayers);keypoints.clear();// 为每一层图像分配空间，并且计算每一层图像的相关属性int index = 0, middleIndex = 0, step = SAMPLE_STEP0;for( int octave = 0; octave < nOctaves; octave++ ){for( int layer = 0; layer < nOctaveLayers+2; layer++ ){// 积分图像要比对应的原始图像的宽和高都大一个像素，所以这里要减1dets[index].create( (sum.rows-1)/step, (sum.cols-1)/step, CV_32F );traces[index].create( (sum.rows-1)/step, (sum.cols-1)/step, CV_32F );// SURF_HAAR_SIZE0为第一组第一层图像的小波大小，默认为9，即9x9 (9x9是对于σ=1.2的高斯二阶微分滤波器而言)// SURF_HAAR_SIZE_INC为层与层之间小波大小增加的步长，这个值必须是偶数，保证每一组的图像小波大小// 要么全是奇数要么全是偶数，进一步确保层与层之间能够正确对齐，默认值为6sizes[index] = (SURF_HAAR_SIZE0 + SURF_HAAR_SIZE_INC*layer) << octave; // 左移1位代表x2// 记录步长sampleSteps[index] = step;// 记录每组中除了第0层和最后一层的图像层对应的索引号if( 0 < layer && layer <= nOctaveLayers )middleIndices[middleIndex++] = index;index++;}// 下一组与上一组图像之间的采样步长是两倍的关系，初始步长为1step *= 2;}// 计算每一层的hessian矩阵行列式值和迹，详细计算过程请见2.1.1 calcLayerDetAndTrace()函数注释parallel_for_( Range(0, nTotalLayers),SURFBuildInvoker(sum, sizes, sampleSteps, dets, traces) );// 寻找hessian矩阵所有行列式值的最大值，详细计算过程请见2.1.2 findMaximaInLayer()函数注释parallel_for_( Range(0, nMiddleLayers),SURFFindInvoker(sum, mask_sum, dets, traces, sizes,sampleSteps, middleIndices, keypoints,nOctaveLayers, hessianThreshold) );// 将特征点按照从“大”到“小”的顺序排序，具体排序规则请见2.1.3 KeypointGreater结构体注释std::sort(keypoints.begin(), keypoints.end(), KeypointGreater());
}

2.1.1 calcLayerDetAndTrace函数注释

计算每一层的hessian矩阵行列式值和迹。

放一个盒型滤波器的图在这里，方便阅读代码时进行参照（图源为文献：Speeded-Up Robust Features (SURF)）。
盒型滤波器是高斯二阶微分进行离散化和裁剪得到的，可用来近似代替高斯滤波板，与原图进行卷积计算。
注意：图中灰色部分值为0，白色部分值为正，黑色部分值为负。

/** 函数：calcLayerDetAndTrace* 功能：计算每一层的hessian矩阵行列式值和迹* 参数：sum               输入  待检测特征点图像的积分图像*      size                输入  每一层图像对应的小波大小/盒型滤波器大小(默认初始大小为9，即9x9，逐层递增)*      sampleStep         输入  每一层图像对应的采样步长(初始步长为1，下一组图像的采样步长是与上一组图像的两倍)*       det                 输出  hessian矩阵的行列式*      trace               输出  hessian矩阵的迹* 返回值：无*/
static void calcLayerDetAndTrace( const Mat& sum, int size, int sampleStep,Mat& det, Mat& trace )
{// 小波模板参数个数，即一个盒型滤波器中，除灰色外，颜色块的个数const int NX=3, NY=3, NXY=4;// 每一行的第一个参数为dx1，第二个参数为dy1，第三个参数为dx2，第四个参数为dy2，第五个参数为权重系数。// 每一行的含义为：权重系数在该模板中的区域位置。(dx1, dy1)表示区域的左上角点坐标，(dx2, dy2)表示区域的右下角点坐标+1。// 具体请参照上面给出的盒型滤波器模板图const int dx_s[NX][5] = { {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} };   //图左一的转置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} };    //图右一// SurfHF为结构体，包含4个整型(p0/p1/p2/p3)和1个浮点型(w)数据// 其中p0~p3分别存放左上、左下、右上和右下位置(便于使用积分图像)，w存放权重值SurfHF Dx[NX], Dy[NY], Dxy[NXY];if( size > sum.rows-1 || size > sum.cols-1 )return;// 重新调整Haar小波模板大小(重新调整盒型滤波器模板大小)，并将其从矩阵格式转换为SurfHF格式// 第一个参数为输入的原始小波模板，第二个参数为输出的调整后的小波模板，第三个参数为小波模板参数个数(除灰色外，颜色块的个数)// 第四个参数为原始小波大小，第五个参数为调整后的小波大小，第六个参数为积分图像宽度resizeHaarPattern( dx_s , Dx , NX , 9, size, sum.cols );resizeHaarPattern( dy_s , Dy , NY , 9, size, sum.cols );resizeHaarPattern( dxy_s, Dxy, NXY, 9, size, sum.cols );// 积分图像sum比原始图像大1个像素int samples_i = 1+(sum.rows-1-size)/sampleStep;int samples_j = 1+(sum.cols-1-size)/sampleStep;// 忽略卷积核超出图像的部分int margin = (size/2)/sampleStep;for( int i = 0; i < samples_i; i++ ){// 积分图像指针指向第i*sampleStep行，第0列，即行采样 (此处是“sampleStep为采样步长”的最好理解)const int* sum_ptr = sum.ptr<int>(i*sampleStep);// 行列式图像指针和迹图像指针指向第i+margin行，第margin列float* det_ptr = &det.at<float>(i+margin, margin);float* trace_ptr = &trace.at<float>(i+margin, margin);for( int j = 0; j < samples_j; j++ ){// 计算Dxx，计算方式为：// 首先，利用原始积分图像，将一个盒型滤波器中不同颜色块区域包含的像素的灰度值之和乘以该颜色块对应的权重系数，// 然后，将一个盒型滤波器中不同颜色块求得的对应的值相加，总和即为所求值。// 代码：d += (sum_ptr[Dx[k].p0] + sum_ptr[Dx[k].p3] - sum_ptr[Dx[k].p1] - sum_ptr[Dx[k].p2]) * f[k].w;float dx  = calcHaarPattern( sum_ptr, Dx , 3 );// 计算Dyyfloat dy  = calcHaarPattern( sum_ptr, Dy , 3 );// 计算Dxyfloat dxy = calcHaarPattern( sum_ptr, Dxy, 4 );// 列采样sum_ptr += sampleStep;// 计算行列式值：Det(H) = Dxx*Dyy - (w*Dxy)^2，其中w约为0.9，目的是为了平衡因使用盒式滤波器近似所带来的误差det_ptr[j] = dx*dy - 0.81f*dxy*dxy;// 计算迹：Tr(H) = Dxx + Dyytrace_ptr[j] = dx + dy;}}
}

2.1.2 findMaximaInLayer函数注释

寻找hessian矩阵所有行列式值的最大值，并进行非极大值抑制，将得到的特征点放入特称点向量输出。

/** 函数：findMaximaInLayer* 功能：寻找hessian矩阵所有行列式值的最大值，并进行非极大值抑制，将得到的特征点放入特称点向量输出* 参数：sum                 输入  待检测特征点图像的积分图像*      mask_sum                输入  掩码图像的积分图像*      dets                    输入  行列式图像*      traces                  输入  迹图像*        sizes                   输入  图像对应的小波大小/盒型滤波器大小(默认初始大小为9，即9x9，逐层递增)*      keypoints             输出  检测到的特征点*        octave                  输入  图像组号*       layer                   输入  中间层(除第0层和最后一层外)的层号*         hessianThreshold        输入  hessian阈值*      sampleStep              输入  图像对应的采样步长(初始步长为1，下一组图像的采样步长是与上一组图像的两倍)* 返回值：无*/
void SURFFindInvoker::findMaximaInLayer( const Mat& sum, const Mat& mask_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} };SurfHF Dm;// 找到层号对应的小波/盒型滤波器大小int size = sizes[layer];// 积分图像sum比原始图像大1个像素int layer_rows = (sum.rows-1)/sampleStep;int layer_cols = (sum.cols-1)/sampleStep;// 忽略没有3x3x3大小邻域的像素int margin = (sizes[layer+1]/2)/sampleStep+1;// 如果掩码图像的积分图像不为空，就重新调整Haar小波模板大小(重新调整盒型滤波器模板大小)，并将其转换为SurfHF格式if( !mask_sum.empty() )resizeHaarPattern( dm, &Dm, NM, 9, size, mask_sum.cols );// Mat中step为图像每一行中所有元素的字节总量，// elemSize()返回的是以8位(一个字节)为一个单位，乘以通道数和8位的整数倍，表示矩阵中每一个元素的数据大小，// 此处的step表示一幅图像一行有多少个元素。int step = (int)(dets[layer].step/dets[layer].elemSize());for( int i = margin; i < layer_rows - margin; i++ ){// 行列式图像和迹图像指针均指向第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++ ){// 获取行列式图像第i行第j列的值float val0 = det_ptr[j];// 如果该值大于用户设定的或默认的hessian阈值if( val0 > hessianThreshold ){// 获取小波/盒型滤波器在原始积分图像中的起始坐标(此处包含一些整型除法，若想简化以下步骤，请务必检查计算结果!!)int sum_i = sampleStep*(i-(size/2)/sampleStep);int sum_j = sampleStep*(j-(size/2)/sampleStep);// 获取最大值附近的3x3x3邻域，存放在矩阵N9[3][9]中，其中，极大值存放在N9[1][4]中。const float *det1 = &dets[layer-1].at<float>(i, j);   // (i, j)点上一层的行列式值const float *det2 = &dets[layer].at<float>(i, j);  // (i, j)点当前层的行列式值const float *det3 = &dets[layer+1].at<float>(i, j);   // (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( !mask_sum.empty() ){const int* mask_ptr = &mask_sum.at<int>(sum_i, sum_j);float mval = calcHaarPattern( mask_ptr, &Dm, 1 );// 如果计算出来的值小于0.5，则跳过后续步骤，继续循环if( mval < 0.5 )continue;}// 非极大值抑制，即判断val0在其3x3x3邻域中是否为最大值。val0在N9[1][4]的位置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(float x[x坐标], float y[y坐标], float _size[小波模板大小/盒型滤波器大小], // float _angle=-1[角度], float _response=0[响应值], int _octave=0[图像组号], // int _class_id=-1[分类号，迹值大于0则为1，小于0则为-1])KeyPoint kpt( center_j, center_i, (float)sizes[layer],-1, val0, octave, CV_SIGN(trace_ptr[j]) );// 在3x3x3邻域内插值最大值位置// 插值方式：// 第一步：求x/y/s方向的负一阶导数，得到向量b(-dx, -dy, -ds)；// 第二步：求二阶导dxx，dxy，dxs，dxy，dyy，dys，dxs，dys，dss，组成3x3的矩阵A；// 第三步：采用LU分解的方法解Ax=b，得到x (此处，x是插值的极值坐标相对于初始估计的偏移量)；// 第四步：判断x向量中是否每个值均不为0，且绝对值小于1，// 如果满足条件，则对传入的特征点kpt的坐标值和size属性进行更改，并返回1，否则，返回0。int ds = size - sizes[layer-1];int interp_ok = interpolateKeypoint( N9, sampleStep, sampleStep, ds, kpt );// 如果插值成功，则将该点作为特征点放入特征点向量(注意：有时特征点插值函数也会返回负值)if( interp_ok  ){/*printf( "KeyPoint %f %f %d\n", point.pt.x, point.pt.y, point.size );*/cv::AutoLock lock(findMaximaInLayer_m);keypoints.push_back(kpt);}}}}}
}

2.1.3 KeypointGreater结构体注释

关于排序特征点“大小”排序。

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;// 第四顺位，按特征点的Y坐标排序，Y坐标越大，该特征点越“大”if(kp1.pt.y < kp2.pt.y) return false;if(kp1.pt.y > kp2.pt.y) return true;// 第五顺位，按特征点的X坐标排序，X坐标越小，该特征点越“大”return kp1.pt.x < kp2.pt.x;}
};

2.2 特征点描述

特征点描述代码分为两部分，第一部分是构造函数，用于参数的赋值、空间的分配以及一些预处理；第二部分是真正的执行函数，进行特征描述符的计算，返回输出描述子。

2.2.1 构造函数

构造函数中除了赋值、空间分配外，最主要的就是调用getGaussianKernel()函数得到高斯模板系数。

放一个特征点邻域模板图，如下图所示。其中，(0, 0)为特征点的坐标，其它点为该特征点的邻域点。

/** 函数：SURFInvoker构造函数* 功能：为特征点描述工作做准备* 参数：_img             输入  原始图像*       _sum                输入  原始图像对应的积分图像*        _keypoints          输入  用于存放特征点的向量*      _descriptors       输入  用于存放特征点描述符的向量*      _extended           输入  是否扩展描述符的维度(true为128维，false为64维)*         _upright            输入  是否为每一个特征点计算方向，默认为false(false为计算，true为不计算)* 返回值：无*/
SURFInvoker( const Mat& _img, const Mat& _sum,vector<KeyPoint>& _keypoints, Mat& _descriptors,bool _extended, bool _upright )
{// 赋值keypoints = &_keypoints;descriptors = &_descriptors;img = &_img;sum = &_sum;extended = _extended;upright = _upright;// 定义在以ORI_RADIUS为半径的圆中的网格点数量(ORI_RADIUS = 6)，或者可以理解为以ORI_RADIUS为半径的圆的最小外接矩形的大小const int nOriSampleBound = (2*ORI_RADIUS+1)*(2*ORI_RADIUS+1);// 分配空间apt.resize(nOriSampleBound); // 存储特征点的每一个邻域点的相对坐标，邻域点具体分布请参考上面给出的特征点邻域模板图aptw.resize(nOriSampleBound);    // 存储特征点的每一个邻域点对应的二维高斯权重值DW.resize(PATCH_SZ*PATCH_SZ);  // 存储用来跟特征点描述符相乘的二维高斯滤波系数// 计算高斯模板尺寸为2*ORI_RADIUS+1 = 2*6+1 = 13、高斯标准差为SURF_ORI_SIGMA = 2.5的高斯滤波系数// 详细计算过程见getGaussianKernel函数注释Mat G_ori = getGaussianKernel( 2*ORI_RADIUS+1, SURF_ORI_SIGMA, CV_32F );nOriSamples = 0;// 计算二维高斯滤波系数(二维高斯滤波系数 = 两个一维高斯滤波系数相乘)，生成后面将用于计算方向的坐标和权重for( int i = -ORI_RADIUS; i <= ORI_RADIUS; i++ ){for( int j = -ORI_RADIUS; j <= ORI_RADIUS; j++ ){if( i*i + j*j <= ORI_RADIUS*ORI_RADIUS ){// 存储特征点的每一个邻域点的相对坐标，邻域点具体分布请参考上面给出的特征点邻域模板图apt[nOriSamples] = cvPoint(i,j);// 存储特征点的每一个邻域点对应的二维高斯权重值aptw[nOriSamples++] = G_ori.at<float>(i+ORI_RADIUS,0) * G_ori.at<float>(j+ORI_RADIUS,0);}}}// nOriSamples为邻域点的个数，nOriSampleBound以ORI_RADIUS为半径的圆的最小外接矩形的大小/矩形区域的网格数// 判断向量大小，确保不会发生内存溢出问题CV_Assert( nOriSamples <= nOriSampleBound );// 计算高斯模板尺寸为PATCH_SZ = 20(PATHC_SZ为特征描述符的大小)、高斯标准差为SURF_DESC_SIGMA = 3.3的高斯滤波系数Mat G_desc = getGaussianKernel( PATCH_SZ, SURF_DESC_SIGMA, CV_32F );// 计算二维高斯滤波系数，生成特征点描述符的权重，在执行函数中会用到for( int i = 0; i < PATCH_SZ; i++ ){for( int j = 0; j < PATCH_SZ; j++ )// 存储用来与特征点描述符相乘的二维高斯滤波系数DW[i*PATCH_SZ+j] = G_desc.at<float>(i,0) * G_desc.at<float>(j,0);}
}/** 函数：getGaussianKernel* 功能：计算并返回nx1大小的高斯滤波系数，计算公式：ΣGi = α * e^(-((i-(n-1)/2)^2)/(2*sigma^2))* 参数：n               输入  高斯滤波系数的尺寸，这个参数必须是正奇数*        sigma           输入  高斯标准差，如果该值为非正数，则sigma = 0.3*((n-1)*0.5 - 1) + 0.8*      ktype           输入  高斯滤波系数的类型，可以为CV_32f或CV_64F* 返回值：返回归一化的大小为nx1的高斯滤波系数*/
cv::Mat cv::getGaussianKernel( int n, double sigma, int ktype )
{const int SMALL_GAUSSIAN_SIZE = 7;// 构建一个简单的高斯滤波系数矩阵，其每一行对称且和为1static const float small_gaussian_tab[][SMALL_GAUSSIAN_SIZE] ={{1.f},{0.25f, 0.5f, 0.25f},{0.0625f, 0.25f, 0.375f, 0.25f, 0.0625f},{0.03125f, 0.109375f, 0.21875f, 0.28125f, 0.21875f, 0.109375f, 0.03125f}};// 判断高斯滤波系数的尺寸n是否为奇数且小于7，高斯标准差是否为非正值// 如果满足条件，则使用事先构建好的高斯滤波系数矩阵；否则，重新计算系数const float* fixed_kernel = n % 2 == 1 && n <= SMALL_GAUSSIAN_SIZE && sigma <= 0 ?small_gaussian_tab[n>>1] : 0;CV_Assert( ktype == CV_32F || ktype == CV_64F );Mat kernel(n, 1, ktype);float* cf = (float*)kernel.data;    // CV_32F类型double* cd = (double*)kernel.data;  // CV_64F类型// 如果sigma为非正值，那么就使用高斯滤波系数的尺寸计算得到新的sigma值double sigmaX = sigma > 0 ? sigma : ((n-1)*0.5 - 1)*0.3 + 0.8;// 事先计算出计算公式中指数的分母部分：scale2X = -(2*sigma^2)double scale2X = -0.5/(sigmaX*sigmaX);double sum = 0;// 计算高斯滤波系数int i;for( i = 0; i < n; i++ ){// 计算出计算公式中指数的分子部分(不带平方)：x = i-(n-1)/2double x = i - (n-1)*0.5;// 若在之前进行判断时满足条件，则使用事先构建好的高斯滤波系数矩阵中的值；否则，按照计算公式重新计算// t = e^(scale2X*x*x) = e^(-((i-(n-1)/2)^2)/(2*sigma^2))double t = fixed_kernel ? (double)fixed_kernel[i] : std::exp(scale2X*x*x);// 根据用户所需要的高斯滤波系数的类型赋值，并求和Σ，sum = Σt = Σ e^(-((i-(n-1)/2)^2)/(2*sigma^2))if( ktype == CV_32F ){cf[i] = (float)t;sum += cf[i];}else{cd[i] = t;sum += cd[i];}}// 将高斯滤波系数进行归一化计算，保证所有高斯滤波系数加起来的和值为1sum = 1./sum;for( i = 0; i < n; i++ ){if( ktype == CV_32F )cf[i] = (float)(cf[i]*sum);elsecd[i] *= sum;}// 返回归一化的大小为nx1的高斯滤波系数return kernel;
}

2.2.2 执行函数

特征点描述真正的执行函数，进行特征描述符的计算。

Haar梯度小波模板如下图所示，用于计算X方向与Y方向的响应（图源为文献：Speeded-Up Robust Features (SURF)），放在此处以供参考，便于代码阅读。
图中，黑色部分权重值为-1，白色部分权重值为1。

主方向搜索图示如下，（图源为文献：Speeded-Up Robust Features (SURF)）。
主方向计算方法：以60度为滑动窗口，5度为步进步长，将特征点邻域（圆形）遍历搜索一周，查找主方向，主方向是加权梯度和的幅值最大的方向。

void operator()(const Range& range) const
{// 构建X方向与Y方向梯度的Haar小波模板，请参考上面给出的Haar梯度小波模板图，该模板原始大小为4const int NX=2, NY=2;const int dx_s[NX][5] = {{0, 0, 2, 4, -1}, {2, 0, 4, 4, 1}};const int dy_s[NY][5] = {{0, 0, 4, 2, 1}, {0, 2, 4, 4, -1}};// 在数组长度中用nOriSampleBound比用nOriSamples优化效果更好，可能是因为在编译时nOriSampleBound已知为常数const int nOriSampleBound =(2*ORI_RADIUS+1)*(2*ORI_RADIUS+1);// 分配空间// X/Y数组分别存储邻域点X/Y方向的梯度与二维高斯权重值相乘后的结果，即加权梯度// angle数组存放每一个邻域点的加权梯度对应的角度float X[nOriSampleBound], Y[nOriSampleBound], angle[nOriSampleBound];// PATCH数组用于存放20*20大小的原始图像灰度值(先获取20s*20s大小的原始图像数据，再缩放至20*20大小)uchar PATCH[PATCH_SZ+1][PATCH_SZ+1];// DX用于存放X方向的梯度加权值，DY用于存放Y方向的梯度加权值float DX[PATCH_SZ][PATCH_SZ], DY[PATCH_SZ][PATCH_SZ];CvMat matX = cvMat(1, nOriSampleBound, CV_32F, X);CvMat matY = cvMat(1, nOriSampleBound, CV_32F, Y);CvMat _angle = cvMat(1, nOriSampleBound, CV_32F, angle);Mat _patch(PATCH_SZ+1, PATCH_SZ+1, CV_8U, PATCH);// 确定特征描述符的维度大小 128或64int dsize = extended ? 128 : 64;// 获取特征点中size(Haar小波模板大小/盒型滤波器模板大小)属性最大值int k, k1 = range.start, k2 = range.end;float maxSize = 0;for( k = k1; k < k2; k++ ){maxSize = std::max(maxSize, (*keypoints)[k].size);}// s = 1.2∗L/9，其中L = maxSize，s表示特征点的尺度// imaxSize是特征描述符的区域大小PATCH_SZ相对最大尺度进行缩放后的结果int imaxSize = std::max(cvCeil((PATCH_SZ+1)*maxSize*1.2f/9.0f), 1);// imaxSize*imaxSize为最大特征描述符的窗口大小Ptr<CvMat> winbuf = cvCreateMat( 1, imaxSize*imaxSize, CV_8U );// 开始为每一个特征点计算主方向和特征描述符for( k = k1; k < k2; k++ ){int i, j, kk, nangle;float* vec;SurfHF dx_t[NX], dy_t[NY];// 获取当前特征点对应的Haar小波模板大小/盒型滤波器模板大小以及特征点坐标KeyPoint& kp = (*keypoints)[k];float size = kp.size;Point2f center = kp.pt;// 计算尺度s，用来计算主方向和构建特征点描述符的采样间隔和Haar小波大小是相对于尺度s来定义的float s = size*1.2f/9.0f;// 为了定义主方向，X方向的梯度和Y方向的梯度使用大小为4s的小波在半径为6s的圆中进行采样// 保证梯度小波大小为偶数，使小波模板以其中心点对称，此处梯度小波模板大小为4s(cvRound函数将传入值进行四舍五入取整)int grad_wav_size = 2*cvRound( 2*s );if( sum->rows < grad_wav_size || sum->cols < grad_wav_size ){// 当梯度小波模板大小过大时，采样梯度无意义，标记该特征点，后面删除kp.size = -1;continue;}float descriptor_dir = 360.f - 90.f;// 开始计算特征点的主方向(upright默认为false)if (upright == 0){// 重新调整梯度小波模板大小，并将其转换为SurfHF格式resizeHaarPattern( dx_s, dx_t, NX, 4, grad_wav_size, sum->cols );resizeHaarPattern( dy_s, dy_t, NY, 4, grad_wav_size, sum->cols );// nOriSamples为邻域点的个数for( kk = 0, nangle = 0; kk < nOriSamples; kk++ ){// 首先，特征点坐标 + 经过尺度缩放后的邻域点相对坐标 -> 得到邻域点绝对坐标(此时得到的坐标在梯度小波模板中心)；// 然后，邻域点绝对坐标 - 梯度小波模板大小的一半 -> 将邻域绝对坐标移动到梯度小波模板的左上角int x = cvRound( center.x + apt[kk].x*s - (float)(grad_wav_size-1)/2 );int y = cvRound( center.y + apt[kk].y*s - (float)(grad_wav_size-1)/2 );// 若坐标超过积分图像边界，则不做处理if( y < 0 || y >= sum->rows - grad_wav_size ||x < 0 || x >= sum->cols - grad_wav_size )continue;// 指针指向(移至梯度小波模板左上角后的)邻域点在积分图像中的坐标(x, y)const int* ptr = &sum->at<int>(y, x);// 分别计算该邻域点在X方向和Y方向的梯度float vx = calcHaarPattern( ptr, dx_t, 2 );float vy = calcHaarPattern( ptr, dy_t, 2 );// X方向和Y方向的梯度分别与其对应的权值相乘，存入数组X[nangle] = vx*aptw[kk];Y[nangle] = vy*aptw[kk];nangle++;}// 若特征点过于接近图像边缘，则无法得出梯度结果，主方向也无法计算得出，故跳过并标记该特征点，后面删除if( nangle == 0 ){kp.size = -1;continue;}matX.cols = matY.cols = _angle.cols = nangle;// 笛卡尔坐标系 -> 极坐标系 (cvCartToPolar函数输入X坐标和Y坐标，输出幅值和角度，此处没有要幅值，只要了角度)// 角度计算公式：angle(I) = atan2(y(I); x(I))*[180/pi]cvCartToPolar( &matX, &matY, 0, &_angle, 1 );float bestx = 0, besty = 0, descriptor_mod = 0;// SURF_ORI_SEARCH_INC = 5，以60度为滑动窗口，5度为步进步长，将特征点邻域遍历搜索一周，查找主方向// 请参考上面给出的主方向搜索图示for( i = 0; i < 360; i += SURF_ORI_SEARCH_INC ){float sumx = 0, sumy = 0, temp_mod;for( j = 0; j < nangle; j++ ){int d = std::abs(cvRound(angle[j]) - i);if( d < ORI_WIN/2 || d > 360-ORI_WIN/2 ){sumx += X[j];sumy += Y[j];}}// 求当前方向的幅值大小temp_mod = sumx*sumx + sumy*sumy;// 若当前方向的幅值大小比历史记录的幅值大小更大，则更新历史记录的幅值大小if( temp_mod > descriptor_mod ){descriptor_mod = temp_mod;bestx = sumx;besty = sumy;}}// 该特征点的主方向就是历史记录的幅值最大的方向，为什么是-besty???descriptor_dir = fastAtan2( -besty, bestx );}// 至此，特征点的主方向计算完毕，记录该特征点的主方向kp.angle = descriptor_dir;if( !descriptors || !descriptors->data )continue;// 开始计算特征描述符// 准备工作一：以特征点为中心，选取一个正方形窗口，方向为特征点的主方向，边长为20s，并向其中填充原始图像数据// 特征描述符窗口边长为20sint win_size = (int)((PATCH_SZ+1)*s);// 确保当前特征描述符窗口大小不超过最大特征描述符窗口大小，防止内存溢出CV_Assert( winbuf->cols >= win_size*win_size );Mat win(win_size, win_size, CV_8U, winbuf->data.ptr);// 如果前面计算了特征点的方向if( !upright ){// 度 -> 弧度descriptor_dir *= (float)(CV_PI/180);float sin_dir = -std::sin(descriptor_dir);    // 注意，这里sin取了负号float cos_dir =  std::cos(descriptor_dir);float win_offset = -(float)(win_size-1)/2;// 计算特征描述符矩形框的起始点在原始图像中的位置// 进行仿射变换 [start_x] = [cos_dir -sin_dir][win_offset] + [center.x]//             [start_y]   [sin_dir  cos_dir][win_offset]   [center.y]float start_x = center.x + win_offset*cos_dir + win_offset*sin_dir;float start_y = center.y - win_offset*sin_dir + win_offset*cos_dir;uchar* WIN = win.data;
#if 0// 最近邻版本(更快)for( i = 0; i < win_size; i++, start_x += sin_dir, start_y += cos_dir ){float pixel_x = start_x;float pixel_y = start_y;for( j = 0; j < win_size; j++, pixel_x += cos_dir, pixel_y -= sin_dir ){// 保证x和y不超过图像边界int x = std::min(std::max(cvRound(pixel_x), 0), img->cols-1);int y = std::min(std::max(cvRound(pixel_y), 0), img->rows-1);// 存储一个特征描述符窗口大小的原始图像像素值WIN[i*win_size + j] = img->at<uchar>(y, x);}}
#else// 双线性插值版本int ncols1 = img->cols-1, nrows1 = img->rows-1;// 获取一行图像的字节数size_t imgstep = img->step;for( i = 0; i < win_size; i++, start_x += sin_dir, start_y += cos_dir ){double pixel_x = start_x;double pixel_y = start_y;for( j = 0; j < win_size; j++, pixel_x += cos_dir, pixel_y -= sin_dir ){int ix = cvFloor(pixel_x), iy = cvFloor(pixel_y);// 如果ix和iy没有超过图像最大边界，就使用双线性插值的计算结果充特征描述符矩形框if( (unsigned)ix < (unsigned)ncols1 &&(unsigned)iy < (unsigned)nrows1 ){float a = (float)(pixel_x - ix), b = (float)(pixel_y - iy);const uchar* imgptr = &img->at<uchar>(iy, ix);// 双线性插值得到图像数据WIN[i*win_size + j] = (uchar)cvRound(imgptr[0]*(1.f - a)*(1.f - b) +imgptr[1]*a*(1.f - b) +imgptr[imgstep]*(1.f - a)*b +imgptr[imgstep+1]*a*b);}else{// 如果ix和iy超过了图像最大边界，就使用最近邻像素填充特征描述符矩形框// 保证x和y不超过图像边界int x = std::min(std::max(cvRound(pixel_x), 0), ncols1);int y = std::min(std::max(cvRound(pixel_y), 0), nrows1);// 直接存储原始图像在该坐标点的灰度值数据WIN[i*win_size + j] = img->at<uchar>(y, x);}}}
#endif}// 如果前面没有计算特征点的方向else{// 构建特征描述符矩形框。此时descriptor_dir = 90 grad，sin_dir = 1，cos_dir = 0float win_offset = -(float)(win_size-1)/2;// 起始点为左下角int start_x = cvRound(center.x + win_offset);int start_y = cvRound(center.y - win_offset);uchar* WIN = win.data;for( i = 0; i < win_size; i++, start_x++ ){int pixel_x = start_x;int pixel_y = start_y;for( j = 0; j < win_size; j++, pixel_y-- ){// 防止坐标越界int x = MAX( pixel_x, 0 );int y = MAX( pixel_y, 0 );x = MIN( x, img->cols-1 );y = MIN( y, img->rows-1 );// 提取对应坐标的图像数据，填充矩形框(y是行坐标，x是列坐标)WIN[i*win_size + j] = img->at<uchar>(y, x);}}}// 将特征描述符窗口大小由20s*20s缩放为20*20，那么每个像素的大小就是s，便于使用2s大小的小波进行梯度计算// 此处采用的缩放方式是：区域插值resize(win, _patch, _patch.size(), 0, 0, INTER_AREA);// 准备工作二：使用2s大小的梯度小波计算X和Y方向的梯度for( i = 0; i < PATCH_SZ; i++ )for( j = 0; j < PATCH_SZ; j++ ){// 获取权值float dw = DW[i*PATCH_SZ + j];// 计算计算X和Y方向的梯度，并乘上权值float vx = (PATCH[i][j+1] - PATCH[i][j] + PATCH[i+1][j+1] - PATCH[i+1][j])*dw;float vy = (PATCH[i+1][j] - PATCH[i][j] + PATCH[i+1][j+1] - PATCH[i][j+1])*dw;// 分别存储X方向和Y方向的加权梯度值DX[i][j] = vx;DY[i][j] = vy;}// 正式计算特征描述符vec = descriptors->ptr<float>(k);   // 指针指向特征描述符矩阵的第k行，代表是第k个特征点的描述符for( kk = 0; kk < dsize; kk++ )vec[kk] = 0;double square_mag = 0;if( extended ){// 128维描述符// 20*20的特征描述符窗口被划分为4*4=16个子块for( i = 0; i < 4; i++ )for( j = 0; j < 4; j++ ){// 每个子块中有5*5个大小为s的像素点for(int y = i*5; y < i*5+5; y++ ){for(int x = j*5; x < j*5+5; x++ ){float tx = DX[y][x], ty = DY[y][x];if( ty >= 0 ){// 统计Y方向的加权梯度值大于等于0时的Σx和Σ|x|vec[0] += tx;vec[1] += (float)fabs(tx);} else {// 统计Y方向的加权梯度值小于0时的Σx和Σ|x|vec[2] += tx;vec[3] += (float)fabs(tx);}if ( tx >= 0 ){// 统计X方向的加权梯度值大于等于0时的Σy和Σ|y|vec[4] += ty;vec[5] += (float)fabs(ty);} else {// 统计X方向的加权梯度值小于0时的Σy和Σ|y|vec[6] += ty;vec[7] += (float)fabs(ty);}}}// 计算幅值// (vec[0])^2+(vec[1])^2+(vec[2])^2+(vec[3])^2+(vec[4])^2+(vec[5])^2+(vec[6])^2+(vec[7])^2for( kk = 0; kk < 8; kk++ )square_mag += vec[kk]*vec[kk];// 指针偏移，计算下一个子块vec += 8;}}else{// 64维描述符// 20*20的特征描述符窗口被划分为4*4=16个子块for( i = 0; i < 4; i++ )for( j = 0; j < 4; j++ ){// 每个子块中有5*5个大小为s的像素点for(int y = i*5; y < i*5+5; y++ ){for(int x = j*5; x < j*5+5; x++ ){// 统计Σx、Σy、Σ|x|和Σ|y|float tx = DX[y][x], ty = DY[y][x];vec[0] += tx; vec[1] += ty;vec[2] += (float)fabs(tx); vec[3] += (float)fabs(ty);}}// 计算幅值 (Σx)^2 + (Σy)^2 + (Σ|x|)^2 + (Σ|y|)^2for( kk = 0; kk < 4; kk++ )square_mag += vec[kk]*vec[kk];// 指针偏移，计算下一个子块vec+=4;}}// 将特征描述符单位化，使其具有对比度不变性vec = descriptors->ptr<float>(k);float scale = (float)(1./(sqrt(square_mag) + DBL_EPSILON));for( kk = 0; kk < dsize; kk++ )vec[kk] *= scale;}
}

SURF算法之Opencv代码详解相关推荐

kmeans python interation flag_机器学习经典算法-logistic回归代码详解
一.算法简要我们希望有这么一种函数:接受输入然后预测出类别,这样用于分类.这里,用到了数学中的sigmoid函数,sigmoid函数的具体表达式和函数图象如下: 可以较为清楚的看到,当输入的x小于0 ...
天津理工大学《操作系统》实验二，存储器的分配与回收算法实现，代码详解，保姆式注释讲解
天津理工大学<操作系统>实验二,存储器的分配与回收算法实现,代码详解,保姆式注释讲解实验内容 1．本实验是模拟操作系统的主存分配,运用可变分区的存储管理算法设计主存分配和回收程序,并不 ...
MeanTeacher文章解读+算法流程+核心代码详解
MeanTeacher 本博客仅做算法流程疏导,具体细节请参见原文原文原文链接点这里 Github 代码 Github代码点这里解读论文解读点这里算法流程代码详解 train_transf ...
python自然语言处理实战核心技术与算法——HMM模型代码详解
本人初学NLP,当我看着<python自然语言处理实战核心技术与算法>书上这接近200行的代码看着有点头皮发麻,于是我读了接近一天基本把每行代码的含义给读的个七七八八,考虑到可能会有人和我 ...
opencv实现camshift算法，以及代码详解
大家好哦,小编来讲程序啦,好好看哦,慢慢体会...一.说明实验介绍本次实验将使用利用 OpenCV 来实现对视频中感兴趣的动态物体的追踪. 首先我贴出我的水下追踪效果图吧在图中我们可以看到鱼的追踪轨 ...
【MATLAB】Parzen窗与K近邻算法原理与代码详解
文章目录 1.非参数估计原理 2.Parzen窗 2.1.算法原理 2.2.Matlab实现与参数探究 3.K近邻 3.1.算法原理 3.2.Matlab实现与参数探究 1.非参数估计原理 \qqua ...
FixMatch文章解读+算法流程+核心代码详解
FixMatch 本博客仅做算法流程疏导,具体细节请参见原文原文查看原文点这里 Github代码 Github代码点这里解读 FixMatch算法抓住了半监督算法的两个重要观点,第一个是一致性正 ...
K均值(K-means)聚类算法原理与代码详解
0. 算法原理: 上述过程简单描述: a: 初始数据 b: 选择质点 c: 根据质点划分 d: 求均值,更新质心点 e: 划分 f: 更新质心点 1. 代码实现: # K means 教程# 0. 引 ...
KNN算法原理和代码详解
原理有这样一条河流like that,河流的左边是rich 人家,河流的右边是poor 人家,这时新搬来一家小甲,这个算法是看小甲是有钱人家还是没钱人家. 要解决这个问题,那么就可以说立着他最近的几 ...
基于Matlab的遗传算法优化BP神经网络的算法实现（附算法介绍与代码详解）
目录一.内容提要二.算法简介 2.1 遗传算法(Genetic Algorithm,GA) 2.2 BP(Back Propagation)神经网络三.实例计算四.代码解读代码运行代码获取 ...

SURF算法之Opencv代码详解