SURF特征匹原理及源代码
SURF是对SIFT的改进，相对于SIFT，主要优点是速度更快，更适合做实时特征检查。
一、SURF原理：
相对于SIFT，SUFT采用Hssian算法检测关键点，很大程度上提高了程序的运行速度，同时，在尺度空间的构建上，SIFT通过改变高斯卷积核的大小，构建不同的组，下文进行详细介绍。
二、SURF的实现：
与SIFT步骤相同，可分为四步：尺度空间构建和极值检测、特征点精确定位和过滤、特征方向赋值和特征点描述。
尺度空间构建和极值检测：
SIFT通过高斯滤波后的同组图像中的每层相减得到DOG空间来提取特征点，SURF通过Harris矩阵（极大提升速度的原因），求取X和Y的导数，进行角点的大致判断。下面简要介绍实现步骤：
Harris矩阵的表示： 
上式中的L表示在经过高斯滤波后的图像的像素值的大小，Lxx、Lxy、Lyy分别表示在x方向的二次偏导导、在x和y方向上的二次偏导和在y方向的二次偏导。在进行Harris矩阵求取时，需满足尺度不变特性。即需先进行高斯滤波，得到多尺度的金字塔图像（具体构造会在下面介绍）。根据上式可得到Harris矩阵表示为行列式的值的大小：
Det(H) = LxxLyy - LxyLxy
根据图像像素点的离散性可知，图像像素值的求导可通过相邻像素点像素值的相减得到：
Lx = L(x+1,y) - L(x,y) Lxx = L(x+1,y) - L(x,y) - L(x,y)+L(x-1,y) = L(x+1,y)+L(x-1,y)-2L(x,y)
同理可求得Lxx、Lxy和Lyy的表达式，上式中的L均为经过高斯滤波后的像素值大小。由上面表达式可知，得到Harris响应值，需要先经过高斯平滑和二次求导两步操作，SURF算法作者采用箱式滤波模板近似代替上面两步，加快了求导的速度。下面暂时介绍SURF中箱式滤波原理（这个箱式滤波，在SURF中的使用方法，，，头疼、自闭、脱发！！！！！！！）：
SURF中箱式滤波：
箱式滤波可简单理解为将像素点周围划分为多个区域，再将每个局部区域内的像素值相加再乘以一个权重值后，进行相加。
上图中，第一幅表示为经过高斯滤波后在y方向的二阶求导，第二幅图像为经过高斯滤波后，在xy方向的求导，第三幅图像表示为箱式滤波在y方向的求导，第四幅图表示箱式滤波为在xy方向上的求导。在采用箱式滤波的过程中，针对每个局部区域像素值的计算，可通过积分图像进行快速计算，加快了计算的速度。（积分图像在另一篇博客中已有介绍）
在采用箱式卷积核替换后，Harris响应值的计算可用下式表示：
Det(H) = DxxDyy -(wDxy)2 //其中D表示盒子模板与图像的近似卷积 其中w为更好模拟而设的变量（原算法中建议取0.9）
通过判断 Det(H)结果的符号，可简要判断是否为局部极值点，若为负值，则不是局部极值点，若为正值，则可能归类为局部极值点。
尺寸空间构建：
相对于SIFT，SURF构建的图像金字塔不同组图像的尺寸大小相同，不同的是卷积核的大小（在此，，，查了很多博客，，，发现都写成高斯卷积核大小，，，很想骂人，，，在开始写为高斯卷积核，在计算Harris相应值时，又改为箱式滤波核，纠结了很久，看了不少博客，最终应该表示为箱式卷积核的大小）。下图为一个简要概括：
非极值抑制：
通过上面得到极值点，将其响应值大小与上下两层响应值大小进行比较，进行判断是否为极值。（原理和SIFT相同）
主方向的计算：
为确保得到的描述子具有旋转不变性，需对每一个特征点分配一个主方向。具体方法为：以特征点为中心，以6s（s为特征点的尺度）为半径的圆形内，对图像进行Haar小波响应运算（相当于对图像进行梯度运算，但采用积分图像，更有效率。Haar小波响应运算暂不介绍，如有兴趣，可网上百度查询）。后设置一个以方向为中心，以pi/3为张角的扇形滑动窗口，通过一定的弧度（大概0.2）为步长，旋转这个滑动窗口，对窗口内的图像Haar小波的响应值进行累加。主方向最大的Harr响应累加值对应的方向即为主方向。
SURF特征点描述子计算：
在SURF算法中，也是在特征点周围取一个正方形框，框的边长为20s(s是所检测到该特征点所在的尺度)。该框带方向，方向当然就是上面一步检测出来的主方向。然后把该框分为16个子区域，每个子区域统计25个像素的水平方向和垂直方向的haar小波特征，这里的水平方向和垂直方向是相对于主方向而言的。该haar小波特征为水平方向值之和、水平方向绝对值之和和垂直方向之和、垂直方向绝对值之和。这样每个小区域就有4个值，所以每个特征点就是16*4维的向量，相比SIFT而言，少了一半，这在特征匹配的过程中会大大的加速匹配的速度。
在OpenSURF的实现源码中采用的方式，通过点旋转公式，把点旋转到主方向上并进行最近邻插值的对应点，公式如下：
OPENCV-contrib内接口：
cv::Ptr<cv::xfeatures2d::SURF> surf = cv::xfeatures2d::SURF::create();
完整实现：
#include <opencv2/opencv.hpp> #include <opencv2/xfeatures2d.hpp> int main() { cv::Mat imageL = cv::imread("Ｔ.png"); cv::Mat imageR = cv::imread("-DT.png"); //提取特征点方法 cv::Ptr<cv::xfeatures2d::SURF> surf = cv::xfeatures2d::SURF::create(); //特征点 std::vector<cv::KeyPoint> keyPointL, keyPointR; //单独提取特征点 surf->detect(imageL, keyPointL); surf->detect(imageR, keyPointR); //画特征点 cv::Mat keyPointImageL; cv::Mat keyPointImageR; drawKeypoints(imageL, keyPointL, keyPointImageL, cv::Scalar::all(-1), cv::DrawMatchesFlags::DRAW_RICH_KEYPOINTS); drawKeypoints(imageR, keyPointR, keyPointImageR, cv::Scalar::all(-1), cv::DrawMatchesFlags::DRAW_RICH_KEYPOINTS); //显示窗口 cv::namedWindow("KeyPoints of imageL"); cv::namedWindow("KeyPoints of imageR"); //显示特征点 cv::imshow("KeyPoints of imageL", keyPointImageL); cv::imshow("KeyPoints of imageR", keyPointImageR); //特征点匹配 cv::Mat despL, despR; //提取特征点并计算特征描述子 surf->detectAndCompute(imageL, cv::Mat(), keyPointL, despL); surf->detectAndCompute(imageR, cv::Mat(), keyPointR, despR); std::vector<cv::DMatch> matches; //如果采用flannBased方法 那么 desp通过orb的到的类型不同需要先转换类型 if (despL.type() != CV_32F || despR.type() != CV_32F) { despL.convertTo(despL, CV_32F); despR.convertTo(despR, CV_32F); } cv::Ptr<cv::DescriptorMatcher> matcher = cv::DescriptorMatcher::create("FlannBased"); matcher->match(despL, despR, matches); //计算特征点距离的最大值 double maxDist = 0; for (int i = 0; i < despL.rows; i++) { double dist = matches[i].distance; if (dist > maxDist) maxDist = dist; } //挑选好的匹配点 std::vector< cv::DMatch > good_matches; for (int i = 0; i < despL.rows; i++) { if (matches[i].distance < 0.03*maxDist) { good_matches.push_back(matches[i]); } } cv::Mat imageOutput; cv::drawMatches(imageL, keyPointL, imageR, keyPointR, good_matches, imageOutput); cv::namedWindow("picture of matching"); cv::imshow("picture of matching", imageOutput); cv::waitKey(0); return 0; }
源代码：
struct SURFInvoker { enum{ORI_RADIUS = 6, ORI_WIN = 60, PATCH_SZ = 20}; // Parameters const Mat* img; const Mat* sum; vector<KeyPoint>* keypoints; Mat* descriptors; bool extended; bool upright; // Pre-calculated values int nOriSamples; vector<Point> apt; // 特征点周围用于描述方向的邻域的点 vector<float> aptw; // 描述 方向时的 高斯 权 vector<float> DW; 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; // 用于描述特征点的 方向的 邻域大小： 12*sigma+1 (sigma =1.2) 因为高斯加权的核的参数为2sigma // nOriSampleBound为 矩形框内点的个数 const int nOriSampleBound = (2 * ORI_RADIUS + 1)*(2 * ORI_RADIUS + 1); // 这里把s近似为1 ORI_DADIUS = 6s // 分配大小 apt.resize(nOriSampleBound); aptw.resize(nOriSampleBound); DW.resize(PATCH_SZ*PATCH_SZ); // PATHC_SZ为特征描述子的 区域大小 20s(s 这里初始为1了) /* 计算特征点方向用的 高斯分布 权值与坐标 */ Mat G_ori = getGaussianKernel(2 * ORI_RADIUS + 1, SURF_ORI_SIGMA, CV_32F); // SURF_ORI_SIGMA = 1.2 *2 =2.5 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); // 下面这里有个坐标转换，因为i,j都是从-ORI_RADIUS开始的。 aptw[nOriSamples++] = G_ori.at<float>(i + ORI_RADIUS, 0) * G_ori.at<float>(j + ORI_RADIUS, 0); } } } CV_Assert(nOriSamples <= nOriSampleBound); // nOriSamples为圆形区域内的点，nOriSampleBound是正方形区域的点 /* 用于特征点描述子的高斯 权值 */ Mat G_desc = getGaussianKernel(PATCH_SZ, SURF_DESC_SIGMA, CV_32F); // 用于生成特征描述子的 高斯加权 sigma = 3.3s (s初取1) 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); } /* x与y方向上的 Harr小波,参数为4s */ const 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 } }; float X[nOriSampleBound], Y[nOriSampleBound], angle[nOriSampleBound]; // 用于计算特生点主方向 uchar PATCH[PATCH_SZ + 1][PATCH_SZ + 1]; float DX[PATCH_SZ][PATCH_SZ], DY[PATCH_SZ][PATCH_SZ]; // 20s * 20s区域的 梯度值 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); int dsize = extended ? 128 : 64; int k, k1 = 0, k2 = (int)(*keypoints).size();// k2为Harr小波的 模板尺寸 float maxSize = 0; for (k = k1; k < k2; k++) { maxSize = std::max(maxSize, (*keypoints)[k].size); } // maxSize*1.2/9 表示最大的尺度 s int imaxSize = std::max(cvCeil((PATCH_SZ + 1)*maxSize*1.2f / 9.0f), 1); 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]; KeyPoint& kp = (*keypoints)[k]; float size = kp.size; Point2f center = kp.pt; /* s是当前层的尺度参数 1.2是第一层的参数，9是第一层的模板大小*/ float s = size*1.2f / 9.0f; /* grad_wav_size是 harr梯度模板的大小 边长为 4s */ int grad_wav_size = 2 * cvRound(2 * s); if (sum->rows < grad_wav_size || sum->cols < grad_wav_size) { /* when grad_wav_size is too big, * the sampling of gradient will be meaningless * mark keypoint for deletion. */ kp.size = -1; continue; } float descriptor_dir = 360.f - 90.f; if (upright == 0) { // 这一步 是计算梯度值，先将harr模板放大，再根据积分图计算，与前面求D_x，D_y一致类似 resizeHaarPattern(dx_s, dx_t, NX, 4, grad_wav_size, sum->cols); resizeHaarPattern(dy_s, dy_t, NY, 4, grad_wav_size, sum->cols); 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; const int* ptr = &sum->at<int>(y, x); float vx = calcHaarPattern(ptr, dx_t, 2); float vy = calcHaarPattern(ptr, dy_t, 2); X[nangle] = vx*aptw[kk]; Y[nangle] = vy*aptw[kk]; nangle++; } if (nangle == 0) { // No gradient could be sampled because the keypoint is too // near too one or more of the sides of the image. As we // therefore cannot find a dominant direction, we skip this // keypoint and mark it for later deletion from the sequence. kp.size = -1; continue; } matX.cols = matY.cols = _angle.cols = nangle; // 计算邻域内每个点的 梯度角度 cvCartToPolar(&matX, &matY, 0, &_angle, 1); float bestx = 0, besty = 0, descriptor_mod = 0; for (i = 0; i < 360; i += SURF_ORI_SEARCH_INC) // SURF_ORI_SEARCH_INC 为扇形区域扫描的步长 { float sumx = 0, sumy = 0, temp_mod; for (j = 0; j < nangle; j++) { // d是 分析到的那个点与 现在主方向的偏度 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; // descriptor_mod 是最大峰值 if (temp_mod > descriptor_mod) { descriptor_mod = temp_mod; bestx = sumx; besty = sumy; } } descriptor_dir = fastAtan2(-besty, bestx); } kp.angle = descriptor_dir; if (!descriptors || !descriptors->data) continue; /* 用特征点周围20*s为边长的邻域 计算特征描述子 */ int 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 dir float cos_dir = std::cos(descriptor_dir); float win_offset = -(float)(win_size - 1) / 2; 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; 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); 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 { 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); } } } } else { float 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); WIN[i*win_size + j] = img->at<uchar>(y, x); } } } // Scale the window to size PATCH_SZ so each pixel's size is s. This // makes calculating the gradients with wavelets of size 2s easy resize(win, _patch, _patch.size(), 0, 0, INTER_AREA); // Calculate gradients in x and y with wavelets of size 2s for (i = 0; i < PATCH_SZ; i++) for (j = 0; j < PATCH_SZ; j++) { float dw = DW[i*PATCH_SZ + j]; // 高斯加权系数 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; DX[i][j] = vx; DY[i][j] = vy; } // Construct the descriptor vec = descriptors->ptr<float>(k); for (kk = 0; kk < dsize; kk++) vec[kk] = 0; double square_mag = 0; if (extended) { // 128维描述子，考虑dx与dy的正负号 for (i = 0; i < 4; i++) for (j = 0; j < 4; j++) { // 每个方块内是一个5s * 5s的区域，每个方法由8个特征描述 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) { vec[0] += tx; vec[1] += (float)fabs(tx); } else { vec[2] += tx; vec[3] += (float)fabs(tx); } if (tx >= 0) { vec[4] += ty; vec[5] += (float)fabs(ty); } else { vec[6] += ty; vec[7] += (float)fabs(ty); } } } for (kk = 0; kk < 8; kk++) square_mag += vec[kk] * vec[kk]; vec += 8; } } else { // 64位描述子 for (i = 0; i < 4; i++) for (j = 0; j < 4; j++) { 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]; vec[0] += tx; vec[1] += ty; vec[2] += (float)fabs(tx); vec[3] += (float)fabs(ty); } } for (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; } } };
参考：
https://blog.csdn.net/songzitea/article/details/16986423
https://blog.csdn.net/msq19895070/article/details/24290599
https://wenku.baidu.com/view/d71c9d2e905f804d2b160b4e767f5acfa1c783b4.html?rec_flag=default
https://blog.csdn.net/qq_37764129/article/details/80969515
