SIFT特征提取算法的C++与Matlab实现源码.zip

#include "sift.h" // keypoint void keypoint::coinvertToKeyPoint(vector<keypoint>& kpts, vector<KeyPoint>& KPts) { for (auto kpt : kpts) { KeyPoint KPt; KPt.pt.x = kpt.pt.y * pow(2.f, float(kpt.octave - 1)); KPt.pt.y = kpt.pt.x * pow(2.f, float(kpt.octave - 1)); KPt.size = kpt.scale * pow(2.f, float(kpt.octave - 1)) * 2; KPt.angle = 360 - kpt.angle; KPt.response = kpt.response; KPt.octave = kpt.oct_info; KPts.push_back(KPt); } } // pyramid void pyramid::appendTo(int oct, Mat& img) { Mat temp; img.copyTo(temp); pyr[oct].push_back(temp); } vector<Mat>& pyramid::operator[](int oct) { return pyr[oct]; } // sift // 检测开始 bool sift::detect(const Mat& image, vector<keypoint>& kpts, Mat& fvec) { image.copyTo(img_org); // 计算Octave 默认金字塔最上层 图像长宽较小值为8 Octaves = round(log( float(min(img_org.cols, img_org.rows)) ) / log(2.f) - 2.f); // 图像灰度化 并且 像素转化为float类型 Mat img_gray, img_gray_f; if (img_org.channels() == 3 || img_org.channels() == 4) cvtColor(img_org, img_gray, COLOR_BGR2GRAY); else img_org.copyTo(img_gray); img_gray.convertTo(img_gray_f, CV_32FC1); // INTER_LINEAR 双线性差值 resize(img_gray_f, img, Size(img_gray_f.cols * 2, img_gray_f.rows * 2), 0, 0, INTER_LINEAR); // 插值后滤波 double sigma_init = sqrt(max(Sigma * Sigma - 0.5 * 0.5 * 4, 0.01)); GaussianBlur(img, img, Size(2 * cvCeil(2 * sigma_init) + 1, 2 * cvCeil(2 * sigma_init) + 1), sigma_init, sigma_init); //cout << 1 << endl; // 预处理结束 构建金字塔 存于pyr_G 与 pyr_DoG buildPyramid(); //cout << 2 << endl; // 检测特征点 findFeaturePoints(kpts); //cout << 3 << endl; // 提取特征点处的128维特征向量 calcFeatureVector(kpts, fvec); //cout << 4 << endl; return true; } // 构建金字塔 void sift::buildPyramid() { // 先计算好sigma vector<double> sigma_i(Layers + 1); sigma_i[0] = Sigma; for (int lyr_i = 1; lyr_i < Layers + 1; lyr_i++) // 0 1 2 3 4 5 { double sigma_prev = pow(K, lyr_i - 1) * Sigma; double sigma_curr = K * sigma_prev; sigma_i[lyr_i] = sqrt(sigma_curr*sigma_curr - sigma_prev*sigma_prev); } // 金子塔初始化 pyr_G.clear(); pyr_DoG.clear(); // 生成金字塔 Mat img_i, img_DoG; img.copyTo(img_i); pyr_G.build(Octaves); // 确定金子塔层数 pyr_DoG.build(Octaves); for (int oct_i = 0; oct_i < Octaves; oct_i++) { pyr_G.appendTo(oct_i, img_i); // 每层第一张图像不需要高斯模糊 // 生成一个octave其他高斯模糊图像 同时生成这一个octaveDoG图像 for (int lyr_i = 1; lyr_i < Layers + 1; lyr_i++) // 0 1 2 3 4 5 { GaussianBlur(img_i, img_i, Size(2 * cvCeil(2 * sigma_i[lyr_i]) + 1, 2 * cvCeil(2 * sigma_i[lyr_i]) + 1), sigma_i[lyr_i], sigma_i[lyr_i]); pyr_G.appendTo(oct_i, img_i); subtract(img_i, pyr_G[oct_i][lyr_i - 1], img_DoG, noArray(), CV_32FC1); pyr_DoG.appendTo(oct_i, img_DoG); } // 降采样生成下一个octave的第一张图像 resize(pyr_G[oct_i][Layers - 2], img_i, Size(img_i.cols / 2, img_i.rows / 2), 0, 0, INTER_NEAREST); } } // 检测特征点 包括寻找极值点-> 筛选极值点-> 计算特征点主方向 void sift::findFeaturePoints(vector<keypoint>& kpts) { // 清空kpts kpts.clear(); // 迭代变量 int oct_i, lyr_i, r, c, pr, pc, kpt_i, ang_i, k; // 27个像素点容器 float pxs[27]; // 极值点暂存序列 vector<keypoint> kpts_temp; // 寻找极值点 int threshold = cvFloor(0.5 * 0.04 / S * 255); for (oct_i = 0; oct_i < Octaves; oct_i++) { for (lyr_i = 1; lyr_i < Layers - 1; lyr_i++) // 0 123 4 { const Mat& img_curr = pyr_DoG[oct_i][lyr_i]; // 遍历像素 排除最外层像素 for (r = 1; r < img_curr.rows - 1; r++) { for (c = 1; c < img_curr.cols - 1; c++) { // 取出比较像素块 const Mat& prev = pyr_DoG[oct_i][lyr_i - 1]; const Mat& curr = pyr_DoG[oct_i][lyr_i]; const Mat& next = pyr_DoG[oct_i][lyr_i + 1]; float px = curr.at<float>(r, c); if (abs(px) >= threshold) { // 取出比较像素 -1 0 1 for (pr = -1, k = 0; pr < 2; pr++) for (pc = -1; pc < 2; pc++) { pxs[k] = prev.at<float>(r + pr, c + pc); pxs[k + 1] = curr.at<float>(r + pr, c + pc); pxs[k + 2] = next.at<float>(r + pr, c + pc); k += 3; } if ((px >= pxs[0] && px >= pxs[1] && px >= pxs[2] && px >= pxs[3] && px >= pxs[4] && px >= pxs[5] && px >= pxs[6] && px >= pxs[7] && px >= pxs[8] && px >= pxs[9] && px >= pxs[10] && px >= pxs[11] && px >= pxs[12] && px >= pxs[14] && px >= pxs[15] && px >= pxs[16] && px >= pxs[17] && px >= pxs[18] && px >= pxs[19] && px >= pxs[20] && px >= pxs[21] && px >= pxs[22] && px >= pxs[23] && px >= pxs[24] && px >= pxs[25] && px >= pxs[26]) || (px <= pxs[0] && px <= pxs[1] && px <= pxs[2] && px <= pxs[3] && px <= pxs[4] && px <= pxs[5] && px <= pxs[6] && px <= pxs[7] && px <= pxs[8] && px <= pxs[9] && px <= pxs[10] && px <= pxs[11] && px <= pxs[12] && px <= pxs[14] && px <= pxs[15] && px <= pxs[16] && px <= pxs[17] && px <= pxs[18] && px <= pxs[19] && px <= pxs[20] && px <= pxs[21] && px <= pxs[22] && px <= pxs[23] && px <= pxs[24] && px <= pxs[25] && px <= pxs[26])) { keypoint kpt(oct_i, lyr_i, Point(r, c)); kpts_temp.push_back(kpt);// 如果是极值点先保存 } } } } } } for (kpt_i = 0; kpt_i < kpts_temp.size(); kpt_i++) { if (!filterExtrema(kpts_temp[kpt_i])) continue; vector<float> angs; calcMainOrientation(kpts_temp[kpt_i], angs); for (ang_i = 0; ang_i < angs.size(); ang_i++) { kpts_temp[kpt_i].angle = angs[ang_i]; kpts.push_back(kpts_temp[kpt_i]); } } } // 筛选极值点 若不合格返回false bool sift::filterExtrema(keypoint& kpt) { // 用到的阈值 float biasThreshold = 0.5f; float contrastThreshold = 0.04f; float edgeThreshold = 10.f; float normalscl = 1.f / 255; // 将图像从 0~255 归一化到 0~1 的因子 // 筛选步骤 1 去除与精确极值偏移较大的点 bool isdrop = true; // 是否删除该极值点标识 // 其他筛选步骤用得上的变量 预先声明 // 取极值点坐标 int kpt_r = kpt.pt.x; int kpt_c = kpt.pt.y; int kpt_lyr = kpt.layer; // 导数 Vec3f dD; float dxx, dyy, dss, dxy, dxs, dys; // 偏差值 Vec3f x_hat; // 5次插值逼近真实极值 for (int try_i = 0; try_i < 5; try_i++) { // 取DoG图像 每次都要重新取 因为layer会变 const Mat& DoG_prev = pyr_DoG[kpt.octave][kpt_lyr - 1]; const Mat& DoG_curr = pyr_DoG[kpt.octave][kpt_lyr]; const Mat& DoG_next = pyr_DoG[kpt.octave][kpt_lyr + 1]; // 计算一阶导 dD = Vec3f((DoG_curr.at<float>(kpt_r, kpt_c + 1) - DoG_curr.at<float>(kpt_r, kpt_c - 1)) * normalscl * 0.5f, (DoG_curr.at<float>(kpt_r + 1, kpt_c) - DoG_curr.at<float>(kpt_r - 1, kpt_c)) * normalscl * 0.5f, (DoG_prev.at<float>(kpt_r, kpt_c) - DoG_next.at<float>(kpt_r, kpt_c)) * normalscl * 0.5f); // 计算二阶导(Hessian矩阵) 注意分母为1 dxx = (DoG_curr.at<float>(kpt_r, kpt_c + 1) + DoG_curr.at<float>(kpt_r, kpt_c - 1) - 2 * DoG_curr.at<float>(kpt_r, kpt_c)) * normalscl; dyy = (DoG_curr.at<float>(kpt_r + 1, kpt_c) + DoG_curr.at<float>(kpt_r - 1, kpt_c) - 2 * DoG_curr.at<float>(kpt_r, kpt_c)) * normalscl; dss = (DoG_next.at<float>(kpt_r, kpt_c) + DoG_prev.at<float>(kpt_r, kpt_c) - 2 * DoG_curr.at<float>(kpt_r, kpt_c)) * normalscl; // 混合导 注意分母为4 dxy = (DoG_curr.at<float>(kpt_r + 1, kpt_c + 1) - DoG_curr.at<float>(kpt_r + 1, kpt_c - 1) - DoG_curr.at<float>(kpt_r - 1, kpt_c + 1) + DoG_curr.at<float>(kpt_r - 1, kpt_c - 1)) * normalscl * 0.25f; dxs = (DoG_next.at<float>(kpt_r, kpt_c + 1) - DoG_next.at<float>(kpt_r, kpt_c - 1) - DoG_prev.at<float>(kpt_r, kpt_c + 1) + DoG_prev.at<float>(kpt_r, kpt_c - 1)) * normalscl * 0.25f; dys = (DoG_next.at<float>(kpt_r + 1, kpt_c) - DoG_next.at<float>(kpt_r - 1, kpt_c) - DoG_prev.at<float>(kpt_r + 1, kpt_c) + DoG_prev.at<float>(kpt_r - 1, kpt_c)) * normalscl * 0.25f; // 由二阶导合成Hessian矩阵 Matx33f H(dxx, dxy, dxs, dxy, dyy, dys, dxs, dys, dss); // 求解偏差值 x_hat = Vec3f(H.solve(dD, DECOMP_LU)); for (int x = 0; x < 3; x++)// 注意这里有个 - �