說明：本文旨在給出?SIFT 算法的具體實現，而在 SIFT 詳解上只是做出簡單介紹，在這里可以給大家推薦一篇好文：https://blog.csdn.net/zddblog/article/details/7521424；結合這篇文章和下文的具體代碼實現，我相信你能很快掌握并運用 SIFT 算法，加油哦！！！如有疑問大家可以一起討論學習！！！

一、SIFT算法簡介

????????SIFT，即尺度不變特征變換（Scale-invariant feature transform，SIFT），是用于圖像處理領域的一種描述。這種描述具有尺度不變性，可在圖像中檢測出關鍵點，是一種局部特征描述子。該方法于1999年由David Lowe?首先發表于計算機視覺國際會議（International Conference on Computer Vision，ICCV），2004年再次經David Lowe整理完善后發表于International journal of computer vision（IJCV）。截止2014年8月，該論文單篇被引次數達25000余次。

1、SIFT算法的特點

（1） SIFT特征是圖像的局部特征，其對旋轉、尺度縮放、亮度變化保持不變性，對視角變化、仿射變換、噪聲也保持一定程度的穩定性；

（2）獨特性（Distinctiveness）好，信息量豐富，適用于在海量特征數據庫中進行快速、準確的匹配；

（3）多量性，即使少數的幾個物體也可以產生大量的SIFT特征向量；

（4）高速性，經優化的SIFT匹配算法甚至可以達到實時的要求；

（5）可擴展性，可以很方便的與其他形式的特征向量進行聯合。

2、SIFT算法可以解決的問題

????????目標的自身狀態、場景所處的環境和成像器材的成像特性等因素影響圖像配準/目標識別跟蹤的性能。而SIFT算法在一定程度上可解決：

（1）目標的旋轉、縮放、平移（RST）????????

（2）圖像仿射/投影變換（視點viewpoint）

（3）光照影響（illumination）

（4）目標遮擋（occlusion）

（5）雜物場景（clutter）

（6）噪聲

????????SIFT算法的實質是在不同的尺度空間上查找關鍵點(特征點)，并計算出關鍵點的方向。SIFT所查找到的關鍵點是一些十分突出，不會因光照，仿射變換和噪音等因素而變化的點，如角點、邊緣點、暗區的亮點及亮區的暗點等。?

二、SIFT算法分為如下四步

1. 尺度空間極值檢測

????????搜索所有尺度上的圖像位置。通過高斯微分函數來識別潛在的對于尺度和旋轉不變的興趣點。

2. 關鍵點定位

????????在每個候選的位置上，通過一個擬合精細的模型來確定位置和尺度。關鍵點的選擇依據于它們的穩定程度。

3. 方向確定

????????基于圖像局部的梯度方向，分配給每個關鍵點位置一個或多個方向。所有后面的對圖像數據的操作都相對于關鍵點的方向、尺度和位置進行變換，從而提供對于這些變換的不變性。

4. 關鍵點描述

????????在每個關鍵點周圍的鄰域內，在選定的尺度上測量圖像局部的梯度。這些梯度被變換成一種表示，這種表示允許比較大的局部形狀的變形和光照變化。

三、SIFT具體實現

說明：代碼實現主要有七個源文件，我和大家一樣有時候下載了好幾個資源然而沒有一個能用的，導致最后都不知道要下載哪一個、是否要下載；所以這里直接把代碼貼出來了！

1、mySift.h

#pragma once //防止頭文件重復包含和下面作用一樣#include<iostream> #include<vector> #include<opencv2\core\core.hpp> #include<opencv2\features2d\features2d.hpp>using namespace std; using namespace cv;/*************************定義常量*****************************///高斯核大小和標準差關系，size=2*(GAUSS_KERNEL_RATIO*sigma)+1,經常設置GAUSS_KERNEL_RATIO=2-3之間 const double GAUSS_KERNEL_RATIO = 3;const int MAX_OCTAVES = 8; //金字塔最大組數const float CONTR_THR = 0.04f; //默認是的對比度閾值(D(x))const float CURV_THR = 10.0f; //關鍵點主曲率閾值const float INIT_SIGMA = 0.5f; //輸入圖像的初始尺度const int IMG_BORDER = 2; //圖像邊界忽略的寬度，也可以改為 1 const int MAX_INTERP_STEPS = 5; //關鍵點精確插值次數const int ORI_HIST_BINS = 36; //計算特征點方向直方圖的BINS個數const float ORI_SIG_FCTR = 1.5f; //計算特征點主方向時候，高斯窗口的標準差因子const float ORI_RADIUS = 3 * ORI_SIG_FCTR; //計算特征點主方向時，窗口半徑因子const float ORI_PEAK_RATIO = 0.8f; //計算特征點主方向時，直方圖的峰值比const int DESCR_WIDTH = 4; //描述子直方圖的網格大小(4x4)const int DESCR_HIST_BINS = 8; //每個網格中直方圖角度方向的維度const float DESCR_MAG_THR = 0.2f; //描述子幅度閾值const float DESCR_SCL_FCTR = 3.0f; //計算描述子時，每個網格的大小因子const int SAR_SIFT_GLOH_ANG_GRID = 8; //GLOH網格沿角度方向等分區間個數const int SAR_SIFT_DES_ANG_BINS = 8; //像素梯度方向在0-360度內等分區間個數const float SAR_SIFT_RADIUS_DES = 12.0f; //描述子鄰域半徑 const int Mmax = 8; //像素梯度方向在0-360度內等分區間個數const double T = 100.0; //sobel算子去除冗余特征點的閾值const float SAR_SIFT_GLOH_RATIO_R1_R2 = 0.73f;//GLOH網格中間圓半徑和外圓半徑之比const float SAR_SIFT_GLOH_RATIO_R1_R3 = 0.25f;//GLOH網格最內層圓半徑和外圓半徑之比#define Feature_Point_Minimum 1500 //輸入圖像特征點最小個數#define We 0.2#define Wn 0.5#define Row_num 3#define Col_num 3#define SIFT_FIXPT_SCALE 48 //不理解，后面可查原論文/************************sift類*******************************/ class mySift {public://默認構造函數mySift(int nfeatures = 0, int nOctaveLayers = 3, double contrastThreshold = 0.03,double edgeThreshold = 10, double sigma = 1.6, bool double_size = true) :nfeatures(nfeatures),nOctaveLayers(nOctaveLayers), contrastThreshold(contrastThreshold),edgeThreshold(edgeThreshold), sigma(sigma), double_size(double_size) {}//獲得尺度空間每組中間層數int get_nOctave_layers() { return nOctaveLayers; }//獲得圖像尺度是否擴大一倍bool get_double_size() { return double_size; }//計算金字塔組數int num_octaves(const Mat& image);//生成高斯金字塔第一組，第一層圖像void create_initial_image(const Mat& image, Mat& init_image);//使用 sobel 算子創建高斯金字塔第一層圖像，以減少冗余特征點void sobel_create_initial_image(const Mat& image, Mat& init_image);//創建高斯金字塔void build_gaussian_pyramid(const Mat& init_image, vector<vector<Mat>>& gauss_pyramid, int nOctaves);//創建高斯差分金字塔void build_dog_pyramid(vector<vector<Mat>>& dog_pyramid, const vector<vector<Mat>>& gauss_pyramid);//該函數生成高斯差分金字塔void amplit_orient(const Mat& image, vector<Mat>& amplit, vector<Mat>& orient, float scale, int nums);//DOG金字塔特征點檢測void find_scale_space_extrema(const vector<vector<Mat>>& dog_pyr, const vector<vector<Mat>>& gauss_pyr,vector<KeyPoint>& keypoints);//DOG金字塔特征點檢測，特征點方向細化版void find_scale_space_extrema1(const vector<vector<Mat>>& dog_pyr, vector<vector<Mat>>& gauss_pyr,vector<KeyPoint>& keypoints);//計算特征點的描述子void calc_sift_descriptors(const vector<vector<Mat>>& gauss_pyr, vector<KeyPoint>& keypoints,Mat& descriptors, const vector<Mat>& amplit, const vector<Mat>& orient);//構建空間尺度—主要是為了獲取 amplit 和 orient 在使用 GLOH 描述子的時候使用void build_sar_sift_space(const Mat& image, vector<Mat>& sar_harris_fun, vector<Mat>& amplit, vector<Mat>& orient);//GLOH 計算一個特征描述子void calc_gloh_descriptors(const vector<Mat>& amplit, const vector<Mat>& orient, const vector<KeyPoint>& keys, Mat& descriptors);//特征點檢測void detect(const Mat& image, vector<vector<Mat>>& gauss_pyr, vector<vector<Mat>>& dog_pyr, vector<KeyPoint>& keypoints,vector<vector<vector<float>>>& all_cell_contrasts,vector<vector<float>>& average_contrast, vector<vector<int>>& n_cells, vector<int>& num_cell, vector<vector<int>>& available_n,vector<int>& available_num, vector<KeyPoint>& final_keypoints,vector<KeyPoint>& Final_keypoints, vector<KeyPoint>& Final_Keypoints);//特征點描述void comput_des(const vector<vector<Mat>>& gauss_pyr, vector<KeyPoint>& final_keypoints, const vector<Mat>& amplit,const vector<Mat>& orient, Mat& descriptors);private:int nfeatures; //設定檢測的特征點的個數值,如果此值設置為0，則不影響結果int nOctaveLayers; //每組金字塔中間層數double contrastThreshold; //對比度閾值（D(x)）double edgeThreshold; //特征點邊緣曲率閾值double sigma; //高斯尺度空間初始層的尺度bool double_size; //是否上采樣原始圖像};//注意類結束的分號

2、mySift.cpp

#include"mySift.h"#include<string> #include<cmath> #include<iostream> //輸入輸出 #include<vector> //vector #include<algorithm> #include<numeric> //用于容器元素求和#include<opencv2\core\hal\hal.hpp> #include<opencv2\core\hal\intrin.hpp>#include<opencv2\opencv.hpp> #include<opencv2\core\core.hpp> //opencv基本數據結構 #include<opencv2\highgui\highgui.hpp> //圖像界面 #include<opencv2\imgproc\imgproc.hpp> //基本圖像處理函數 #include<opencv2\features2d\features2d.hpp> //特征提取 #include<opencv2\imgproc\types_c.h> #include<opencv2\videoio.hpp> #include <highgui\highgui.hpp> #include <imgproc\imgproc.hpp>/******************根據輸入圖像大小計算高斯金字塔的組數****************************/ /*image表示原始輸入灰度圖像,inline函數必須在聲明處定義 double_size_image表示是否在構建金字塔之前上采樣原始圖像 */ int mySift::num_octaves(const Mat& image) {int temp;float size_temp = (float)min(image.rows, image.cols);temp = cvRound(log(size_temp) / log((float)2) - 2);if (double_size)temp += 1;if (temp > MAX_OCTAVES) //尺度空間最大組數設置為MAX_OCTAVEStemp = MAX_OCTAVES;return temp; }/************************sobel濾波器計算高斯尺度空間圖像梯度大小*****************************/ void sobelfilter(Mat& image, Mat& G) {// image 是經過歸一化后的數據 (0,1)int rows = image.rows;int cols = image.cols;float dx = 0.0, dy = 0.0;//cv::Mat Gx = cv::Mat::zeros(rows, cols, CV_32FC1); //包含了圖像像素在水平方向上的導數的近似值的圖像//cv::Mat Gy = cv::Mat::zeros(rows, cols, CV_32FC1); //包含了圖像像素在垂直方向上的導數的近似值的圖像G = Mat::zeros(rows, cols, CV_32FC1); //在每個像素點處的灰度大小由Gx和Gy共同決定double v = 0.0, vx, vy;//利用 sobel 算子求梯度幅度圖像for (int i = 0; i < rows; i++){for (int j = 0; j < cols; j++){v = 0.0;if (i == 0 || j == 0 || i == rows - 1 || j == cols - 1){G.at<float>(i, j) = 0.0;}else{float dx = image.at<float>(i - 1, j + 1) - image.at<float>(i - 1, j - 1)+ 2 * image.at<float>(i, j + 1) - 2 * image.at<float>(i, j - 1)+ image.at<float>(i + 1, j + 1) - image.at<float>(i + 1, j - 1);float dy = image.at<float>(i + 1, j - 1) - image.at<float>(i - 1, j - 1)+ 2 * image.at<float>(i + 1, j) - 2 * image.at<float>(i - 1, j) +image.at<float>(i + 1, j + 1) - image.at<float>(i - 1, j + 1);v = abs(dx) + abs(dy); //簡化后 G = |Gx| + |Gy|//保證像素值在有效范圍內v = fmax(v, 0); //返回浮點數中較大的一個v = fmin(v, 255); //返回浮點數中較小的一個if (v > T) //T為閾值等于50G.at<float>(i, j) = (float)v;elseG.at<float>(i, j) = 0.0;}}}//水平方向上的導數的近似值的圖像/*for (int i = 0; i < rows; i++){for (int j = 0; j < cols; j++){vx = 0;if (i == 0 || j == 0 || i == rows - 1 || j == cols - 1)Gx.at<float>(i, j) = 0;else{dx = image.at<float>(i - 1, j + 1) - image.at<float>(i - 1, j - 1)+ 2 * image.at<float>(i, j + 1) - 2 * image.at<float>(i, j - 1)+ image.at<float>(i + 1, j + 1) - image.at<float>(i + 1, j - 1);vx = abs(dx);vx = fmax(vx, 0); vx = fmin(vx, 255);Gx.at<float>(i, j) = (float)vx;}}}*///垂直方向上的導數的近似值的圖像/*for (int i = 0; i < rows; i++){for (int j = 0; j < cols; j++){vy = 0;if (i == 0 || j == 0 || i == rows - 1 || j == cols - 1)Gy.at<float>(i, j) = 0;else{dy = image.at<float>(i + 1, j - 1) - image.at<float>(i - 1, j - 1)+ 2 * image.at<float>(i + 1, j) - 2 * image.at<float>(i - 1, j) +image.at<float>(i + 1, j + 1) - image.at<float>(i - 1, j + 1);vy = abs(dy);vy = fmax(vy, 0); vx = fmin(vy, 255);Gy.at<float>(i, j) = (float)vy;}}}*///cv::imshow("Gx", Gx); // horizontal//cv::imshow("Gy", Gy); // vertical//cv::imshow("G", G); // gradient }/*********該函數根據尺度和窗口半徑生成ROEWA濾波模板************/ /*size表示核半徑，因此核寬度是2*size+1scale表示指數權重參數kernel表示生成的濾波核*/ static void roewa_kernel(int size, float scale, Mat& kernel) {kernel.create(2 * size + 1, 2 * size + 1, CV_32FC1);for (int i = -size; i <= size; ++i){float* ptr_k = kernel.ptr<float>(i + size);for (int j = -size; j <= size; ++j){ptr_k[j + size] = exp(-1.f * (abs(i) + abs(j)) / scale);}} }/************************創建高斯金字塔第一組，第一層圖像************************************/ /*image表示輸入原始圖像init_image表示生成的高斯尺度空間的第一層圖像*/ void mySift::create_initial_image(const Mat& image, Mat& init_image) {Mat gray_image;if (image.channels() != 1)cvtColor(image, gray_image, CV_RGB2GRAY); //轉換為灰度圖像elsegray_image = image.clone();Mat floatImage; //轉換到0-1之間的浮點類型數據歸一化，方便接下來的處理//float_image=(float)gray_image*(1.0/255.0)gray_image.convertTo(floatImage, CV_32FC1, 1.0 / 255.0, 0);double sig_diff = 0;if (double_size){Mat temp_image;//通過插值的方法改變圖像尺寸的大小resize(floatImage, temp_image, Size(2 * floatImage.cols, 2 * floatImage.rows), 0.0, 0.0, INTER_LINEAR);//高斯平滑的標準差，值較大時平滑效果比較明顯sig_diff = sqrt(sigma * sigma - 4.0 * INIT_SIGMA * INIT_SIGMA);//高斯濾波窗口大小選擇很重要，這里選擇(4*sig_diff_1+1)-(6*sig_diff+1)之間，且四舍五入int kernel_width = 2 * cvRound(GAUSS_KERNEL_RATIO * sig_diff) + 1;Size kernel_size(kernel_width, kernel_width);//對圖像進行平滑處理(高斯模糊)，即降低圖像的分辨率,高斯模糊是實現尺度變換的唯一變換核，并其實唯一的線性核GaussianBlur(temp_image, init_image, kernel_size, sig_diff, sig_diff);}else{sig_diff = sqrt(sigma * sigma - 1.0 * INIT_SIGMA * INIT_SIGMA);//高斯濾波窗口大小選擇很重要，這里選擇(4*sig_diff_1+1)-(6*sig_diff+1)之間int kernel_width = 2 * cvRound(GAUSS_KERNEL_RATIO * sig_diff) + 1;Size kernel_size(kernel_width, kernel_width);GaussianBlur(floatImage, init_image, kernel_size, sig_diff, sig_diff);} }/************************使用 sobel 算子創建高斯金字塔第一組，第一層圖像****************************/ //目的是為了減少冗余特征點 void mySift::sobel_create_initial_image(const Mat& image, Mat& init_image) {Mat gray_image, gray_images; //gray_images用于存放經過sobel算子操作后的圖像if (image.channels() != 1)cvtColor(image, gray_image, CV_RGB2GRAY); //轉換為灰度圖像elsegray_image = image.clone();sobelfilter(gray_image, gray_images);Mat floatImage; //轉換到0-1之間的浮點類型數據歸一化，方便接下來的處理//float_image=(float)gray_image*(1.0/255.0)gray_images.convertTo(floatImage, CV_32FC1, 1.0 / 255.0, 0);double sig_diff = 0;if (double_size){Mat temp_image;//通過插值的方法改變圖像尺寸的大小resize(floatImage, temp_image, Size(2 * floatImage.cols, 2 * floatImage.rows), 0.0, 0.0, INTER_LINEAR);//高斯平滑的標準差，值較大時平滑效果比較明顯sig_diff = sqrt(sigma * sigma - 4.0 * INIT_SIGMA * INIT_SIGMA);//高斯濾波窗口大小選擇很重要，這里選擇(4*sig_diff_1+1)-(6*sig_diff+1)之間，且四舍五入int kernel_width = 2 * cvRound(GAUSS_KERNEL_RATIO * sig_diff) + 1;Size kernel_size(kernel_width, kernel_width);//對圖像進行平滑處理(高斯模糊)，即降低圖像的分辨率,高斯模糊是實現尺度變換的唯一變換核，并其實唯一的線性核GaussianBlur(temp_image, init_image, kernel_size, sig_diff, sig_diff);}else{sig_diff = sqrt(sigma * sigma - 1.0 * INIT_SIGMA * INIT_SIGMA);//高斯濾波窗口大小選擇很重要，這里選擇(4*sig_diff_1+1)-(6*sig_diff+1)之間int kernel_width = 2 * cvRound(GAUSS_KERNEL_RATIO * sig_diff) + 1;Size kernel_size(kernel_width, kernel_width);GaussianBlur(floatImage, init_image, kernel_size, sig_diff, sig_diff);} }/**************************生成高斯金字塔*****************************************//*init_image表示已經生成的高斯金字塔第一層圖像gauss_pyramid表示生成的高斯金字塔nOctaves表示高斯金字塔的組數*/ void mySift::build_gaussian_pyramid(const Mat& init_image, vector<vector<Mat>>& gauss_pyramid, int nOctaves) {vector<double> sig;sig.push_back(sigma);double k = pow(2.0, 1.0 / nOctaveLayers); //高斯金字塔每一層的系數 kfor (int i = 1; i < nOctaveLayers + 3; ++i){double prev_sig = pow(k, i - 1) * sigma; //每一個尺度層的尺度double curr_sig = k * prev_sig;//組內每層的尺度坐標計算公式sig.push_back(sqrt(curr_sig * curr_sig - prev_sig * prev_sig));}gauss_pyramid.resize(nOctaves);for (int i = 0; i < nOctaves; ++i){gauss_pyramid[i].resize(nOctaveLayers + 3);}for (int i = 0; i < nOctaves; ++i) //對于每一組{for (int j = 0; j < nOctaveLayers + 3; ++j) //對于組內的每一層{if (i == 0 && j == 0) //第一組，第一層gauss_pyramid[0][0] = init_image;else if (j == 0){resize(gauss_pyramid[i - 1][3], gauss_pyramid[i][0],Size(gauss_pyramid[i - 1][3].cols / 2,gauss_pyramid[i - 1][3].rows / 2), 0, 0, INTER_LINEAR);}else{//高斯濾波窗口大小選擇很重要，這里選擇(4*sig_diff_1+1)-(6*sig_diff+1)之間int kernel_width = 2 * cvRound(GAUSS_KERNEL_RATIO * sig[j]) + 1;Size kernel_size(kernel_width, kernel_width);GaussianBlur(gauss_pyramid[i][j - 1], gauss_pyramid[i][j], kernel_size, sig[j], sig[j]);}}} }/*******************生成高斯差分金字塔，即LOG金字塔*************************/ /*dog_pyramid表示DOG金字塔gauss_pyramin表示高斯金字塔*/ void mySift::build_dog_pyramid(vector<vector<Mat>>& dog_pyramid, const vector<vector<Mat>>& gauss_pyramid) {vector<vector<Mat>>::size_type nOctaves = gauss_pyramid.size();for (vector<vector<Mat>>::size_type i = 0; i < nOctaves; ++i){//用于存放每一個梯度中的所有尺度層vector<Mat> temp_vec;for (auto j = 0; j < nOctaveLayers + 2; ++j){Mat temp_img = gauss_pyramid[i][j + 1] - gauss_pyramid[i][j];temp_vec.push_back(temp_img);}dog_pyramid.push_back(temp_vec);temp_vec.clear();} }/***********生成高斯差分金字塔當前層對應的梯度幅度圖像和梯度方向圖像***********/ /*image為高斯差分金字塔當前層圖像*amplit為當前層梯度幅度圖像*orient為當前層梯度方向圖像*scale當前層尺度*nums為相對底層的層數*/ void mySift::amplit_orient(const Mat& image, vector<Mat>& amplit, vector<Mat>& orient, float scale, int nums) {//分配內存amplit.resize(Mmax * nOctaveLayers);orient.resize(Mmax * nOctaveLayers);int radius = cvRound(2 * scale);Mat kernel; //kernel(2 * radius + 1, 2 * radius + 1, CV_32FC1);roewa_kernel(radius, scale, kernel); //返回濾波核，也即指數部分,存放在矩陣的右下角//四個濾波模板生成Mat W34 = Mat::zeros(2 * radius + 1, 2 * radius + 1, CV_32FC1); //把kernel矩陣下半部分復制到對應部分Mat W12 = Mat::zeros(2 * radius + 1, 2 * radius + 1, CV_32FC1); //把kernel矩陣上半部分復制到對應部分Mat W14 = Mat::zeros(2 * radius + 1, 2 * radius + 1, CV_32FC1); //把kernel矩陣右半部分復制到對應部分Mat W23 = Mat::zeros(2 * radius + 1, 2 * radius + 1, CV_32FC1); //把kernel矩陣左半部分復制到對應部分kernel(Range(radius + 1, 2 * radius + 1), Range::all()).copyTo(W34(Range(radius + 1, 2 * radius + 1), Range::all()));kernel(Range(0, radius), Range::all()).copyTo(W12(Range(0, radius), Range::all()));kernel(Range::all(), Range(radius + 1, 2 * radius + 1)).copyTo(W14(Range::all(), Range(radius + 1, 2 * radius + 1)));kernel(Range::all(), Range(0, radius)).copyTo(W23(Range::all(), Range(0, radius)));//濾波Mat M34, M12, M14, M23;double eps = 0.00001;//float_image為圖像歸一化后的圖像數據，做卷積運算filter2D(image, M34, CV_32FC1, W34, Point(-1, -1), eps);filter2D(image, M12, CV_32FC1, W12, Point(-1, -1), eps);filter2D(image, M14, CV_32FC1, W14, Point(-1, -1), eps);filter2D(image, M23, CV_32FC1, W23, Point(-1, -1), eps);//計算水平梯度和豎直梯度Mat Gx, Gy;log((M14) / (M23), Gx);log((M34) / (M12), Gy);//計算梯度幅度和梯度方向magnitude(Gx, Gy, amplit[nums]); //梯度幅度圖像，平方和開平方phase(Gx, Gy, orient[nums], true); //梯度方向圖像}/***********************該函數計算尺度空間特征點的主方向，用于后面特征點的檢測***************************/ /*image表示特征點所在位置的高斯圖像，后面可對著源碼進行修改pt表示特征點的位置坐標(x,y)scale特征點的尺度n表示直方圖bin個數hist表示計算得到的直方圖函數返回值是直方圖hist中的最大數值*/ static float clac_orientation_hist(const Mat& image, Point pt, float scale, int n, float* hist) {int radius = cvRound(ORI_RADIUS * scale); //特征點鄰域半徑(3*1.5*scale)int len = (2 * radius + 1) * (2 * radius + 1); //特征點鄰域像素總個數（最大值）float sigma = ORI_SIG_FCTR * scale; //特征點鄰域高斯權重標準差(1.5*scale)float exp_scale = -1.f / (2 * sigma * sigma); //權重的指數部分//使用AutoBuffer分配一段內存，這里多出4個空間的目的是為了方便后面平滑直方圖的需要AutoBuffer<float> buffer((4 * len) + n + 4);//X保存水平差分，Y保存數值差分，Mag保存梯度幅度，Ori保存梯度角度，W保存高斯權重float* X = buffer, * Y = buffer + len, * Mag = Y, * Ori = Y + len, * W = Ori + len;float* temp_hist = W + len + 2; //臨時保存直方圖數據for (int i = 0; i < n; ++i)temp_hist[i] = 0.f; //數據清零//計算鄰域像素的水平差分和豎直差分int k = 0;for (int i = -radius; i < radius; ++i){int y = pt.y + i; //鄰域點的縱坐標if (y <= 0 || y >= image.rows - 1)continue;for (int j = -radius; j < radius; ++j){int x = pt.x + j;if (x <= 0 || x >= image.cols - 1)continue;float dx = image.at<float>(y, x + 1) - image.at<float>(y, x - 1); //水平差分float dy = image.at<float>(y + 1, x) - image.at<float>(y - 1, x); //豎直差分//保存水平差分和豎直差分及其對應的權重X[k] = dx;Y[k] = dy;W[k] = (i * i + j * j) * exp_scale;++k;}}len = k; //鄰域內特征點的個數cv::hal::exp(W, W, len); //計算鄰域內所有像素的高斯權重cv::hal::fastAtan2(Y, X, Ori, len, true); //計算鄰域內所有像素的梯度方向，角度范圍0-360度cv::hal::magnitude32f(X, Y, Mag, len); //計算鄰域內所有像素的梯度幅度，計算的是數學意義上的梯度//遍歷鄰域的像素for (int i = 0; i < len; ++i){int bin = cvRound((n / 360.f) * Ori[i]); //利用像素的梯度方向，約束bin的范圍在[0,(n-1)]//像素點梯度方向為360度時，和0°一樣if (bin >= n)bin = bin - n;if (bin < 0)bin = bin + n;temp_hist[bin] = temp_hist[bin] + Mag[i] * W[i]; //統計鄰域內像素各個方向在梯度直方圖的幅值(加權后的幅值)}//平滑直方圖temp_hist[-1] = temp_hist[n - 1];temp_hist[-2] = temp_hist[n - 2];temp_hist[n] = temp_hist[0];temp_hist[n + 1] = temp_hist[1];for (int i = 0; i < n; ++i){hist[i] = (temp_hist[i - 2] + temp_hist[i + 2]) * (1.f / 16.f) +(temp_hist[i - 1] + temp_hist[i + 1]) * (4.f / 16.f) +temp_hist[i] * (6.f / 16.f);}//獲得直方圖中最大值float max_value = hist[0];for (int i = 1; i < n; ++i){if (hist[i] > max_value)max_value = hist[i];}return max_value; }/***********************使用 sobel 濾波器定義的新梯度計算尺度空間特征點的主方向**************************/ static float clac_orientation_hist_2(Mat& image, Point pt, float scale, int n, float* hist) {Mat output_image; //使用 sobel 濾波器計算的圖像的梯度幅度圖像sobelfilter(image, output_image); //使用 sobel 濾波器求高斯差分圖像的梯度幅度圖像int radius = cvRound(ORI_RADIUS * scale); //特征點鄰域半徑(3*1.5*scale)int len = (2 * radius + 1) * (2 * radius + 1); //特征點鄰域像素總個數（最大值）float sigma = ORI_SIG_FCTR * scale; //特征點鄰域高斯權重標準差(1.5*scale)float exp_scale = -1.f / (2 * sigma * sigma); //權重的指數部分//使用AutoBuffer分配一段內存，這里多出4個空間的目的是為了方便后面平滑直方圖的需要AutoBuffer<float> buffer((4 * len) + n + 4);//X保存水平差分，Y保存數值差分，Mag保存梯度幅度，Ori保存梯度角度，W保存高斯權重float* X = buffer, * Y = buffer + len, * Mag = Y, * Ori = Y + len, * W = Ori + len;float* temp_hist = W + len + 2; //臨時保存直方圖數據for (int i = 0; i < n; ++i)temp_hist[i] = 0.f; //數據清零//計算鄰域像素的水平差分和豎直差分int k = 0;for (int i = -radius; i < radius; ++i){int y = pt.y + i; //鄰域點的縱坐標，行if (y <= 0 || y >= output_image.rows - 1)continue;for (int j = -radius; j < radius; ++j){int x = pt.x + j; //鄰域點的縱坐標，列if (x <= 0 || x >= output_image.cols - 1)continue;//float dx = image.at<float>(y, x + 1) - image.at<float>(y, x - 1); //水平差分float dx = output_image.at<float>(y - 1, x + 1) - output_image.at<float>(y - 1, x - 1)+ 2 * output_image.at<float>(y, x + 1) - 2 * output_image.at<float>(y, x - 1)+ output_image.at<float>(y + 1, x + 1) - output_image.at<float>(y + 1, x - 1);float dy = output_image.at<float>(y + 1, x - 1) - output_image.at<float>(y - 1, x - 1)+ 2 * output_image.at<float>(y + 1, x) - 2 * output_image.at<float>(y - 1, x) +output_image.at<float>(y + 1, x + 1) - output_image.at<float>(y - 1, x + 1);/*float dx = image.at<float>(y - 1, x + 1) - image.at<float>(y - 1, x - 1)+ 2 * image.at<float>(y, x + 1) - 2 * image.at<float>(y, x - 1)+ image.at<float>(y + 1, x + 1) - image.at<float>(y + 1, x - 1);float dy = image.at<float>(y + 1, x - 1) - image.at<float>(y - 1, x - 1)+ 2 * image.at<float>(y + 1, x) - 2 * image.at<float>(y - 1, x) +image.at<float>(y + 1, x + 1) - image.at<float>(y - 1, x + 1);*///保存水平差分和豎直差分及其對應的權重X[k] = dx;Y[k] = dy;W[k] = (i * i + j * j) * exp_scale;++k;}}len = k; //鄰域內特征點的個數cv::hal::exp(W, W, len); //計算鄰域內所有像素的高斯權重cv::hal::fastAtan2(Y, X, Ori, len, true); //計算鄰域內所有像素的梯度方向，角度范圍0-360度cv::hal::magnitude32f(X, Y, Mag, len); //計算鄰域內所有像素的梯度幅度，計算的是數學意義上的梯度//遍歷鄰域的像素for (int i = 0; i < len; ++i){int bin = cvRound((n / 360.f) * Ori[i]); //利用像素的梯度方向，約束bin的范圍在[0,(n-1)]//像素點梯度方向為360度時，和0°一樣if (bin >= n)bin = bin - n;if (bin < 0)bin = bin + n;temp_hist[bin] = temp_hist[bin] + Mag[i] * W[i]; //統計鄰域內像素各個方向在梯度直方圖的幅值(加權后的幅值)}//平滑直方圖temp_hist[-1] = temp_hist[n - 1];temp_hist[-2] = temp_hist[n - 2];temp_hist[n] = temp_hist[0];temp_hist[n + 1] = temp_hist[1];for (int i = 0; i < n; ++i){hist[i] = (temp_hist[i - 2] + temp_hist[i + 2]) * (1.f / 16.f) +(temp_hist[i - 1] + temp_hist[i + 1]) * (4.f / 16.f) +temp_hist[i] * (6.f / 16.f);}//獲得直方圖中最大值float max_value = hist[0];for (int i = 1; i < n; ++i){if (hist[i] > max_value)max_value = hist[i];}return max_value; }/******************該函數計算特征點主方向,用于LOGH版本*********************/ /*amplit表示特征點所在層的梯度幅度，即輸入圖像對應像素點的梯度存在了對應位置orient表示特征點所在層的梯度方向，0-360度point表示特征點坐標scale表示特征點的所在層的尺度hist表示生成的直方圖n表示主方向直方圖bin個數*/ static float calc_orient_hist(const Mat& amplit, const Mat& orient, Point pt, float scale, float* hist, int n) {//暫且認為是只進行了下采樣，沒有進行高斯模糊int num_row = amplit.rows;int num_col = amplit.cols;Point point(cvRound(pt.x), cvRound(pt.y));//int radius = cvRound(SAR_SIFT_FACT_RADIUS_ORI * scale);int radius = cvRound(6 * scale);radius = min(radius, min(num_row / 2, num_col / 2));float gauss_sig = 2 * scale; //高斯加權標準差float exp_temp = -1.f / (2 * gauss_sig * gauss_sig); //權重指數部分//鄰域區域int radius_x_left = point.x - radius;int radius_x_right = point.x + radius;int radius_y_up = point.y - radius;int radius_y_down = point.y + radius;//防止越界if (radius_x_left < 0)radius_x_left = 0;if (radius_x_right > num_col - 1)radius_x_right = num_col - 1;if (radius_y_up < 0)radius_y_up = 0;if (radius_y_down > num_row - 1)radius_y_down = num_row - 1;//此時特征點周圍矩形區域相對于本矩形區域的中心int center_x = point.x - radius_x_left;int center_y = point.y - radius_y_up;//矩形區域的邊界，計算高斯權值Range x_rng(-(point.x - radius_x_left), radius_x_right - point.x);Range y_rng(-(point.y - radius_y_up), radius_y_down - point.y);Mat gauss_weight;gauss_weight.create(y_rng.end - y_rng.start + 1, x_rng.end - x_rng.start + 1, CV_32FC1);//求各個像素點的高斯權重for (int i = y_rng.start; i <= y_rng.end; ++i){float* ptr_gauss = gauss_weight.ptr<float>(i - y_rng.start);for (int j = x_rng.start; j <= x_rng.end; ++j)ptr_gauss[j - x_rng.start] = exp((i * i + j * j) * exp_temp);}//索引特征點周圍的像素梯度幅度，梯度方向Mat sub_amplit = amplit(Range(radius_y_up, radius_y_down + 1), Range(radius_x_left, radius_x_right + 1));Mat sub_orient = orient(Range(radius_y_up, radius_y_down + 1), Range(radius_x_left, radius_x_right + 1));//Mat W = sub_amplit.mul(gauss_weight); //加入高斯權重，計算高斯權重時，正確匹配點對反而變少了Mat W = sub_amplit; //沒加高斯權重，梯度幅值//計算直方圖AutoBuffer<float> buffer(n + 4);float* temp_hist = buffer + 2;for (int i = 0; i < n; ++i)temp_hist[i] = 0.f;for (int i = 0; i < sub_orient.rows; i++){float* ptr_1 = W.ptr<float>(i);float* ptr_2 = sub_orient.ptr<float>(i);for (int j = 0; j < sub_orient.cols; j++){if (((i - center_y) * (i - center_y) + (j - center_x) * (j - center_x)) < radius * radius){int bin = cvRound(ptr_2[j] * n / 360.f);if (bin > n)bin = bin - n;if (bin < 0)bin = bin + n;temp_hist[bin] += ptr_1[j];}}}//平滑直方圖，可以防止突變temp_hist[-1] = temp_hist[n - 1];temp_hist[-2] = temp_hist[n - 2];temp_hist[n] = temp_hist[0];temp_hist[n + 1] = temp_hist[1];for (int i = 0; i < n; ++i){hist[i] = (temp_hist[i - 2] + temp_hist[i + 2]) * (1.f / 16.f) +(temp_hist[i - 1] + temp_hist[i + 1]) * (4.f / 16.f) +temp_hist[i] * (6.f / 16.f);}//獲得直方圖中最大值float max_value = hist[0];for (int i = 1; i < n; ++i){if (hist[i] > max_value)max_value = hist[i];}return max_value; }/****************************該函數精確定位特征點位置(x,y,scale)，用于后面特征點的檢測*************************/ /*功能：確定特征點的位置，并通過主曲率消除邊緣相應點,該版本是簡化版dog_pry表示DOG金字塔kpt表示精確定位后該特征點的信息octave表示初始特征點所在的組layer表示初始特征點所在的層row表示初始特征點在圖像中的行坐標col表示初始特征點在圖像中的列坐標nOctaveLayers表示DOG金字塔每組中間層數，默認是3contrastThreshold表示對比度閾值，默認是0.04edgeThreshold表示邊緣閾值，默認是10sigma表示高斯尺度空間最底層圖像尺度，默認是1.6*/ static bool adjust_local_extrema_1(const vector<vector<Mat>>& dog_pyr, KeyPoint& kpt, int octave, int& layer, int& row,int& col, int nOctaveLayers, float contrastThreshold, float edgeThreshold, float sigma) {float xi = 0, xr = 0, xc = 0;int i = 0;for (; i < MAX_INTERP_STEPS; ++i) //最大迭代次數{const Mat& img = dog_pyr[octave][layer]; //當前層圖像索引const Mat& prev = dog_pyr[octave][layer - 1]; //之前層圖像索引const Mat& next = dog_pyr[octave][layer + 1]; //下一層圖像索引//特征點位置x方向，y方向,尺度方向的一階偏導數float dx = (img.at<float>(row, col + 1) - img.at<float>(row, col - 1)) * (1.f / 2.f);float dy = (img.at<float>(row + 1, col) - img.at<float>(row - 1, col)) * (1.f / 2.f);float dz = (next.at<float>(row, col) - prev.at<float>(row, col)) * (1.f / 2.f);//計算特征點位置二階偏導數float v2 = img.at<float>(row, col);float dxx = img.at<float>(row, col + 1) + img.at<float>(row, col - 1) - 2 * v2;float dyy = img.at<float>(row + 1, col) + img.at<float>(row - 1, col) - 2 * v2;float dzz = prev.at<float>(row, col) + next.at<float>(row, col) - 2 * v2;//計算特征點周圍混合二階偏導數float dxy = (img.at<float>(row + 1, col + 1) + img.at<float>(row - 1, col - 1) -img.at<float>(row + 1, col - 1) - img.at<float>(row - 1, col + 1)) * (1.f / 4.f);float dxz = (next.at<float>(row, col + 1) + prev.at<float>(row, col - 1) -next.at<float>(row, col - 1) - prev.at<float>(row, col + 1)) * (1.f / 4.f);float dyz = (next.at<float>(row + 1, col) + prev.at<float>(row - 1, col) -next.at<float>(row - 1, col) - prev.at<float>(row + 1, col)) * (1.f / 4.f);Matx33f H(dxx, dxy, dxz, dxy, dyy, dyz, dxz, dyz, dzz);Vec3f dD(dx, dy, dz);Vec3f X = H.solve(dD, DECOMP_SVD);xc = -X[0]; //x方向偏移量xr = -X[1]; //y方向偏移量xi = -X[2]; //尺度方向偏移量//如果三個方向偏移量都小于0.5，說明已經找到特征點準確位置if (abs(xc) < 0.5f && abs(xr) < 0.5f && abs(xi) < 0.5f)break;//如果其中一個方向的偏移量過大，則刪除該點if (abs(xc) > (float)(INT_MAX / 3) ||abs(xr) > (float)(INT_MAX / 3) ||abs(xi) > (float)(INT_MAX / 3))return false;col = col + cvRound(xc);row = row + cvRound(xr);layer = layer + cvRound(xi);//如果特征點定位在邊界區域，同樣也需要刪除if (layer<1 || layer>nOctaveLayers ||col<IMG_BORDER || col>img.cols - IMG_BORDER ||row<IMG_BORDER || row>img.rows - IMG_BORDER)return false;}//如果i=MAX_INTERP_STEPS，說明循環結束也沒有滿足條件，即該特征點需要被刪除if (i >= MAX_INTERP_STEPS)return false;/**************************再次刪除低響應點(對比度較低的點)********************************///再次計算特征點位置x方向，y方向,尺度方向的一階偏導數//高對比度的特征對圖像的變形是穩定的{const Mat& img = dog_pyr[octave][layer];const Mat& prev = dog_pyr[octave][layer - 1];const Mat& next = dog_pyr[octave][layer + 1];float dx = (img.at<float>(row, col + 1) - img.at<float>(row, col - 1)) * (1.f / 2.f);float dy = (img.at<float>(row + 1, col) - img.at<float>(row - 1, col)) * (1.f / 2.f);float dz = (next.at<float>(row, col) - prev.at<float>(row, col)) * (1.f / 2.f);Matx31f dD(dx, dy, dz);float t = dD.dot(Matx31f(xc, xr, xi));float contr = img.at<float>(row, col) + t * 0.5f; //特征點響應 |D(x~)| 即對比度//Low建議contr閾值是0.03，但是RobHess等建議閾值為0.04/nOctaveLayersif (abs(contr) < contrastThreshold / nOctaveLayers) //閾值設為0.03時特征點數量過多return false;/******************************刪除邊緣響應比較強的點************************************///再次計算特征點位置二階偏導數，獲取特征點出的 Hessian 矩陣，主曲率通過 2X2 的 Hessian 矩陣求出//一個定義不好的高斯差分算子的極值在橫跨邊緣的地方有較大的主曲率而在垂直邊緣的方向有較小的主曲率float v2 = img.at<float>(row, col);float dxx = img.at<float>(row, col + 1) + img.at<float>(row, col - 1) - 2 * v2;float dyy = img.at<float>(row + 1, col) + img.at<float>(row - 1, col) - 2 * v2;float dxy = (img.at<float>(row + 1, col + 1) + img.at<float>(row - 1, col - 1) -img.at<float>(row + 1, col - 1) - img.at<float>(row - 1, col + 1)) * (1.f / 4.f);float det = dxx * dyy - dxy * dxy;float trace = dxx + dyy;//主曲率和閾值的對比判定if (det < 0 || (trace * trace * edgeThreshold >= det * (edgeThreshold + 1) * (edgeThreshold + 1)))return false;/*********到目前為止該特征的滿足上面所有要求，保存該特征點信息***********/kpt.pt.x = ((float)col + xc) * (1 << octave); //相對于最底層的圖像的x坐標kpt.pt.y = ((float)row + xr) * (1 << octave); //相對于最底層圖像的y坐標kpt.octave = octave + (layer << 8); //組號保存在低字節，層號保存在高字節//相對于最底層圖像的尺度kpt.size = sigma * powf(2.f, (layer + xi) / nOctaveLayers) * (1 << octave);kpt.response = abs(contr); //特征點響應值(對比度)return true;}}/****************************該函數精確定位特征點位置(x,y,scale)，用于后面特征點的檢測*************************/ //該版本是 SIFT 原版，檢測得到的特征點數量更多 static bool adjust_local_extrema_2(const vector<vector<Mat>>& dog_pyr, KeyPoint& kpt, int octave, int& layer, int& row,int& col, int nOctaveLayers, float contrastThreshold, float edgeThreshold, float sigma) {const float img_scale = 1.f / (255 * SIFT_FIXPT_SCALE); //SIFT_FIXPT_SCALE=48const float deriv_scale = img_scale * 0.5f;const float second_deriv_scale = img_scale;const float cross_deriv_scale = img_scale * 0.25f;float xi = 0, xr = 0, xc = 0;int i = 0;for (; i < MAX_INTERP_STEPS; ++i) //最大迭代次數{const Mat& img = dog_pyr[octave][layer]; //當前層圖像索引const Mat& prev = dog_pyr[octave][layer - 1]; //之前層圖像索引const Mat& next = dog_pyr[octave][layer + 1]; //下一層圖像索引//計算一階偏導數，通過臨近點差分求得float dx = (img.at<float>(row, col + 1) - img.at<float>(row, col - 1)) * deriv_scale;float dy = (img.at<float>(row + 1, col) - img.at<float>(row - 1, col)) * deriv_scale;float dz = (next.at<float>(row, col) - prev.at<float>(row, col)) * deriv_scale;//計算特征點位置二階偏導數//float v2 = img.at<float>(row, col);float v2 = (float)img.at<float>(row, col) * 2.f;float dxx = (img.at<float>(row, col + 1) + img.at<float>(row, col - 1) - v2) * second_deriv_scale;float dyy = (img.at<float>(row + 1, col) + img.at<float>(row - 1, col) - v2) * second_deriv_scale;float dzz = (prev.at<float>(row, col) + next.at<float>(row, col) - v2) * second_deriv_scale;//計算特征點周圍混合二階偏導數float dxy = (img.at<float>(row + 1, col + 1) + img.at<float>(row - 1, col - 1) -img.at<float>(row + 1, col - 1) - img.at<float>(row - 1, col + 1)) * cross_deriv_scale;float dxz = (next.at<float>(row, col + 1) + prev.at<float>(row, col - 1) -next.at<float>(row, col - 1) - prev.at<float>(row, col + 1)) * cross_deriv_scale;float dyz = (next.at<float>(row + 1, col) + prev.at<float>(row - 1, col) -next.at<float>(row - 1, col) - prev.at<float>(row + 1, col)) * cross_deriv_scale;Matx33f H(dxx, dxy, dxz, dxy, dyy, dyz, dxz, dyz, dzz);Vec3f dD(dx, dy, dz);Vec3f X = H.solve(dD, DECOMP_SVD);xc = -X[0]; //x方向偏移量xr = -X[1]; //y方向偏移量xi = -X[2]; //尺度方向偏移量//如果三個方向偏移量都小于0.5，說明已經找到特征點準確位置if (abs(xc) < 0.5f && abs(xr) < 0.5f && abs(xi) < 0.5f)break;//如果其中一個方向的偏移量過大，則刪除該點if (abs(xc) > (float)(INT_MAX / 3) ||abs(xr) > (float)(INT_MAX / 3) ||abs(xi) > (float)(INT_MAX / 3))return false;col = col + cvRound(xc);row = row + cvRound(xr);layer = layer + cvRound(xi);//如果特征點定位在邊界區域，同樣也需要刪除if (layer<1 || layer>nOctaveLayers ||col < IMG_BORDER || col >= img.cols - IMG_BORDER ||row < IMG_BORDER || row >= img.rows - IMG_BORDER)return false;}//如果i=MAX_INTERP_STEPS，說明循環結束也沒有滿足條件，即該特征點需要被刪除if (i >= MAX_INTERP_STEPS)return false;/**************************再次刪除低響應點(對比度較低的點)********************************///再次計算特征點位置x方向，y方向,尺度方向的一階偏導數//高對比度的特征對圖像的變形是穩定的const Mat& img = dog_pyr[octave][layer];const Mat& prev = dog_pyr[octave][layer - 1];const Mat& next = dog_pyr[octave][layer + 1];float dx = (img.at<float>(row, col + 1) - img.at<float>(row, col - 1)) * deriv_scale;float dy = (img.at<float>(row + 1, col) - img.at<float>(row - 1, col)) * deriv_scale;float dz = (next.at<float>(row, col) - prev.at<float>(row, col)) * deriv_scale;Matx31f dD(dx, dy, dz);float t = dD.dot(Matx31f(xc, xr, xi));float contr = img.at<float>(row, col) + t * 0.5f; //特征點響應 |D(x~)| 即對比度//Low建議contr閾值是0.03，但是RobHess等建議閾值為0.04/nOctaveLayersif (abs(contr) < contrastThreshold / nOctaveLayers) //閾值設為0.03時特征點數量過多return false;/******************************刪除邊緣響應比較強的點************************************///再次計算特征點位置二階偏導數，獲取特征點出的 Hessian 矩陣，主曲率通過 2X2 的 Hessian 矩陣求出//一個定義不好的高斯差分算子的極值在橫跨邊緣的地方有較大的主曲率而在垂直邊緣的方向有較小的主曲率float v2 = (float)img.at<float>(row, col) * 2.f;float dxx = (img.at<float>(row, col + 1) + img.at<float>(row, col - 1) - v2) * second_deriv_scale;float dyy = (img.at<float>(row + 1, col) + img.at<float>(row - 1, col) - v2) * second_deriv_scale;float dxy = (img.at<float>(row + 1, col + 1) + img.at<float>(row - 1, col - 1) -img.at<float>(row + 1, col - 1) - img.at<float>(row - 1, col + 1)) * cross_deriv_scale;float det = dxx * dyy - dxy * dxy;float trace = dxx + dyy;//主曲率和閾值的對比判定if (det <= 0 || (trace * trace * edgeThreshold >= det * (edgeThreshold + 1) * (edgeThreshold + 1)))return false;/*********保存該特征點信息***********/kpt.pt.x = ((float)col + xc) * (1 << octave); //高斯金字塔坐標根據組數擴大相應的倍數kpt.pt.y = ((float)row + xr) * (1 << octave);// SIFT 描述子kpt.octave = octave + (layer << 8) + (cvRound((xi + 0.5) * 255) << 16); //特征點被檢測出時所在的金字塔組kpt.size = sigma * powf(2.f, (layer + xi) / nOctaveLayers) * (1 << octave) * 2; //關鍵點鄰域直徑kpt.response = abs(contr); //特征點響應值(對比度)return true;}/************該函數在DOG金字塔上進行特征點檢測，特征點精確定位，刪除低對比度點，刪除邊緣響應較大點**********/ /*dog_pyr表示高斯金字塔 原始 SIFT 算子gauss_pyr表示高斯金字塔keypoints表示檢測到的特征點*/ void mySift::find_scale_space_extrema(const vector<vector<Mat>>& dog_pyr, const vector<vector<Mat>>& gauss_pyr,vector<KeyPoint>& keypoints) {int nOctaves = (int)dog_pyr.size(); //子八度個數//Low文章建議thresholdThreshold是0.03，Rob Hess等人使用0.04/nOctaveLayers作為閾值float threshold = (float)(contrastThreshold / nOctaveLayers);const int n = ORI_HIST_BINS; //n=36float hist[n];KeyPoint kpt;keypoints.clear(); //先清空keypointsfor (int i = 0; i < nOctaves; ++i) //對于每一組{for (int j = 1; j <= nOctaveLayers; ++j) //對于組內每一層{const Mat& curr_img = dog_pyr[i][j]; //當前層const Mat& prev_img = dog_pyr[i][j - 1]; //上一層const Mat& next_img = dog_pyr[i][j + 1]; //下一層int num_row = curr_img.rows;int num_col = curr_img.cols; //獲得當前組圖像的大小size_t step = curr_img.step1(); //一行元素所占字節數//遍歷每一個尺度層中的有效像素，像素值for (int r = IMG_BORDER; r < num_row - IMG_BORDER; ++r){const float* curr_ptr = curr_img.ptr<float>(r); //指向的是第 r 行的起點，返回的是 float 類型的像素值const float* prev_ptr = prev_img.ptr<float>(r - 1);const float* next_ptr = next_img.ptr<float>(r + 1);for (int c = IMG_BORDER; c < num_col - IMG_BORDER; ++c){float val = curr_ptr[c]; //當前中心點響應值//開始檢測特征點if (abs(val) > threshold &&((val > 0 && val >= curr_ptr[c - 1] && val >= curr_ptr[c + 1] &&val >= curr_ptr[c - step - 1] && val >= curr_ptr[c - step] && val >= curr_ptr[c - step + 1] &&val >= curr_ptr[c + step - 1] && val >= curr_ptr[c + step] && val >= curr_ptr[c + step + 1] &&val >= prev_ptr[c] && val >= prev_ptr[c - 1] && val >= prev_ptr[c + 1] &&val >= prev_ptr[c - step - 1] && val >= prev_ptr[c - step] && val >= prev_ptr[c - step + 1] &&val >= prev_ptr[c + step - 1] && val >= prev_ptr[c + step] && val >= prev_ptr[c + step + 1] &&val >= next_ptr[c] && val >= next_ptr[c - 1] && val >= next_ptr[c + 1] &&val >= next_ptr[c - step - 1] && val >= next_ptr[c - step] && val >= next_ptr[c - step + 1] &&val >= next_ptr[c + step - 1] && val >= next_ptr[c + step] && val >= next_ptr[c + step + 1]) ||(val < 0 && val <= curr_ptr[c - 1] && val <= curr_ptr[c + 1] &&val <= curr_ptr[c - step - 1] && val <= curr_ptr[c - step] && val <= curr_ptr[c - step + 1] &&val <= curr_ptr[c + step - 1] && val <= curr_ptr[c + step] && val <= curr_ptr[c + step + 1] &&val <= prev_ptr[c] && val <= prev_ptr[c - 1] && val <= prev_ptr[c + 1] &&val <= prev_ptr[c - step - 1] && val <= prev_ptr[c - step] && val <= prev_ptr[c - step + 1] &&val <= prev_ptr[c + step - 1] && val <= prev_ptr[c + step] && val <= prev_ptr[c + step + 1] &&val <= next_ptr[c] && val <= next_ptr[c - 1] && val <= next_ptr[c + 1] &&val <= next_ptr[c - step - 1] && val <= next_ptr[c - step] && val <= next_ptr[c - step + 1] &&val <= next_ptr[c + step - 1] && val <= next_ptr[c + step] && val <= next_ptr[c + step + 1]))){//++numKeys;//獲得特征點初始行號，列號，組號，組內層號int octave = i, layer = j, r1 = r, c1 = c;if (!adjust_local_extrema_1(dog_pyr, kpt, octave, layer, r1, c1,nOctaveLayers, (float)contrastThreshold,(float)edgeThreshold, (float)sigma)){continue; //如果該初始點不滿足條件，則不保存改點}float scale = kpt.size / float(1 << octave); //該特征點相對于本組的尺度//max_hist值對應的方向為主方向float max_hist = clac_orientation_hist(gauss_pyr[octave][layer], Point(c1, r1), scale, n, hist);//大于mag_thr值對應的方向為輔助方向float mag_thr = max_hist * ORI_PEAK_RATIO; //主峰值 80% 的方向作為輔助方向//遍歷直方圖中的 36 個binfor (int i = 0; i < n; ++i){int left = i > 0 ? i - 1 : n - 1;int right = i < n - 1 ? i + 1 : 0;//創建新的特征點，大于主峰值 80% 的方向，賦值給該特征點，作為一個新的特征點；即有多個特征點，位置、尺度相同，方向不同if (hist[i] > hist[left] && hist[i] > hist[right] && hist[i] >= mag_thr){float bin = i + 0.5f * (hist[left] - hist[right]) / (hist[left] + hist[right] - 2 * hist[i]);bin = bin < 0 ? n + bin : bin >= n ? bin - n : bin;kpt.angle = (360.f / n) * bin; //原始 SIFT 算子使用的特征點的主方向0-360度keypoints.push_back(kpt); //保存該特征點}}}}}}}//cout << "初始滿足要求特征點個數是: " << numKeys << endl; }/************該函數在DOG金字塔上進行特征點檢測，特征點精確定位，刪除低對比度點，刪除邊緣響應較大點**********/ //對特征點進行方向的細化 + 增加更多的主方向版本 ——此時細化是對最后要給關鍵點進行賦值時的細化 //還可以考慮直接對方向直方圖進行細化 void mySift::find_scale_space_extrema1(const vector<vector<Mat>>& dog_pyr, vector<vector<Mat>>& gauss_pyr,vector<KeyPoint>& keypoints) {int nOctaves = (int)dog_pyr.size(); //子八度個數//Low文章建議thresholdThreshold是0.03，Rob Hess等人使用0.04/nOctaveLayers作為閾值float threshold = (float)(contrastThreshold / nOctaveLayers);const int n = ORI_HIST_BINS; //n=36float hist[n];KeyPoint kpt;vector<Mat> amplit; //存放高斯差分金字塔每一層的梯度幅度圖像vector<Mat> orient; //存放高斯差分金字塔每一層的梯度方向圖像keypoints.clear(); //先清空keypointsfor (int i = 0; i < nOctaves; ++i) //對于每一組{for (int j = 1; j <= nOctaveLayers; ++j) //對于組內每一層{const Mat& curr_img = dog_pyr[i][j]; //當前層const Mat& prev_img = dog_pyr[i][j - 1]; //上一層const Mat& next_img = dog_pyr[i][j + 1]; //下一層int num_row = curr_img.rows;int num_col = curr_img.cols; //獲得當前組圖像的大小size_t step = curr_img.step1(); //一行元素所占字節數//遍歷每一個尺度層中的有效像素，像素值for (int r = IMG_BORDER; r < num_row - IMG_BORDER; ++r){const float* curr_ptr = curr_img.ptr<float>(r); //指向的是第 r 行的起點，返回的是 float 類型的像素值const float* prev_ptr = prev_img.ptr<float>(r);const float* next_ptr = next_img.ptr<float>(r);for (int c = IMG_BORDER; c < num_col - IMG_BORDER; ++c){float val = curr_ptr[c]; //當前中心點響應值//開始檢測特征點if (abs(val) > threshold &&((val > 0 && val >= curr_ptr[c - 1] && val >= curr_ptr[c + 1] &&val >= curr_ptr[c - step - 1] && val >= curr_ptr[c - step] && val >= curr_ptr[c - step + 1] &&val >= curr_ptr[c + step - 1] && val >= curr_ptr[c + step] && val >= curr_ptr[c + step + 1] &&val >= prev_ptr[c] && val >= prev_ptr[c - 1] && val >= prev_ptr[c + 1] &&val >= prev_ptr[c - step - 1] && val >= prev_ptr[c - step] && val >= prev_ptr[c - step + 1] &&val >= prev_ptr[c + step - 1] && val >= prev_ptr[c + step] && val >= prev_ptr[c + step + 1] &&val >= next_ptr[c] && val >= next_ptr[c - 1] && val >= next_ptr[c + 1] &&val >= next_ptr[c - step - 1] && val >= next_ptr[c - step] && val >= next_ptr[c - step + 1] &&val >= next_ptr[c + step - 1] && val >= next_ptr[c + step] && val >= next_ptr[c + step + 1]) ||(val < 0 && val <= curr_ptr[c - 1] && val <= curr_ptr[c + 1] &&val <= curr_ptr[c - step - 1] && val <= curr_ptr[c - step] && val <= curr_ptr[c - step + 1] &&val <= curr_ptr[c + step - 1] && val <= curr_ptr[c + step] && val <= curr_ptr[c + step + 1] &&val <= prev_ptr[c] && val <= prev_ptr[c - 1] && val <= prev_ptr[c + 1] &&val <= prev_ptr[c - step - 1] && val <= prev_ptr[c - step] && val <= prev_ptr[c - step + 1] &&val <= prev_ptr[c + step - 1] && val <= prev_ptr[c + step] && val <= prev_ptr[c + step + 1] &&val <= next_ptr[c] && val <= next_ptr[c - 1] && val <= next_ptr[c + 1] &&val <= next_ptr[c - step - 1] && val <= next_ptr[c - step] && val <= next_ptr[c - step + 1] &&val <= next_ptr[c + step - 1] && val <= next_ptr[c + step] && val <= next_ptr[c + step + 1]))){//++numKeys;//獲得特征點初始行號，列號，組號，組內層號int octave = i, layer = j, r1 = r, c1 = c, nums = i * nOctaves + j;if (!adjust_local_extrema_2(dog_pyr, kpt, octave, layer, r1, c1,nOctaveLayers, (float)contrastThreshold,(float)edgeThreshold, (float)sigma)){continue; //如果該初始點不滿足條件，則不保存改點}float scale = kpt.size / float(1 << octave); //該特征點相對于本組的尺度//計算梯度幅度和梯度方向//amplit_orient(curr_img, amplit, orient, scale, nums);//max_hist值對應的方向為主方向float max_hist = clac_orientation_hist(gauss_pyr[octave][layer], Point(c1, r1), scale, n, hist);//float max_hist = calc_orient_hist(amplit[nums], orient[nums], Point2f(c1, r1), scale, hist, n);大于mag_thr值對應的方向為輔助方向//float mag_thr = max_hist * ORI_PEAK_RATIO; //主峰值 80% 的方向作為輔助方向//增加更多的主方向，以增加特征點對梯度差異的魯棒性float sum = 0.0; //直方圖對應的幅值之和float mag_thr = 0.0; //判斷是否為主方向的閾值for (int i = 0; i < n; ++i){sum += hist[i];}mag_thr = 0.5 * (1.0 / 36) * sum;//遍歷直方圖中的 36 個binfor (int i = 0; i < n; ++i){int left = i > 0 ? i - 1 : n - 1;int right = i < n - 1 ? i + 1 : 0;//創建新的特征點，大于主峰值 80% 的方向，賦值給該特征點，作為一個新的特征點；即有多個特征點，位置、尺度相同，方向不同if (hist[i] > hist[left] && hist[i] > hist[right] && hist[i] >= mag_thr){float bin = i + 0.5f * (hist[left] - hist[right]) / (hist[left] + hist[right] - 2 * hist[i]);bin = bin < 0 ? n + bin : bin >= n ? bin - n : bin;//修改的地方，特征點的主方向修改為了0-180度，相當于對方向做了一個細化float angle = (360.f / n) * bin;if (angle >= 1 && angle <= 180){kpt.angle = angle;}else if (angle > 180 && angle < 360){kpt.angle = 360 - angle;}//kpt.angle = (360.f / n) * bin; //原始 SIFT 算子使用的特征點的主方向0-360度keypoints.push_back(kpt); //保存該特征點}}}}}}}//cout << "初始滿足要求特征點個數是: " << numKeys << endl; }/*該函數生成matlab中的meshgrid函數*/ /*x_range表示x方向的范圍 y_range表示y方向的范圍 X表示生成的沿x軸變化的網格 Y表示生成沿y軸變換的網格 */ static void meshgrid(const Range& x_range, const Range& y_range, Mat& X, Mat& Y) {int x_start = x_range.start, x_end = x_range.end;int y_start = y_range.start, y_end = y_range.end;int width = x_end - x_start + 1;int height = y_end - y_start + 1;X.create(height, width, CV_32FC1);Y.create(height, width, CV_32FC1);for (int i = y_start; i <= y_end; i++){float* ptr_1 = X.ptr<float>(i - y_start);for (int j = x_start; j <= x_end; ++j)ptr_1[j - x_start] = j * 1.0f;}for (int i = y_start; i <= y_end; i++){float* ptr_2 = Y.ptr<float>(i - y_start);for (int j = x_start; j <= x_end; ++j)ptr_2[j - x_start] = i * 1.0f;} }/******************************計算一個特征點的描述子***********************************/ /*gauss_image表示特征點所在的高斯層main_angle表示特征點的主方向，角度范圍是0-360度pt表示特征點在高斯圖像上的坐標，相對與本組，不是相對于最底層scale表示特征點所在層的尺度，相對于本組，不是相對于最底層d表示特征點鄰域網格寬度n表示每個網格內像素梯度角度等分個數descriptor表示生成的特征點的描述子*/ static void calc_sift_descriptor(const Mat& gauss_image, float main_angle, Point2f pt,float scale, int d, int n, float* descriptor) {Point ptxy(cvRound(pt.x), cvRound(pt.y)); //坐標取整float cos_t = cosf(-main_angle * (float)(CV_PI / 180)); //把角度轉化為弧度，計算主方向的余弦float sin_t = sinf(-main_angle * (float)(CV_PI / 180)); //把角度轉化為弧度，計算主方向的正弦float bins_per_rad = n / 360.f; // n = 8 ，梯度直方圖分為 8 個方向float exp_scale = -1.f / (d * d * 0.5f); //權重指數部分float hist_width = DESCR_SCL_FCTR * scale; //特征點鄰域內子區域邊長，子區域的邊長int radius = cvRound(hist_width * (d + 1) * sqrt(2) * 0.5f);//特征點鄰域半徑(d+1)*(d+1)，四舍五入int rows = gauss_image.rows, cols = gauss_image.cols; //當前高斯層行、列信息//特征點鄰域半徑radius = min(radius, (int)sqrt((double)rows * rows + cols * cols));cos_t = cos_t / hist_width;sin_t = sin_t / hist_width;int len = (2 * radius + 1) * (2 * radius + 1); //鄰域內總像素數，為后面動態分配內存使用int histlen = (d + 2) * (d + 2) * (n + 2); //值為 360 AutoBuffer<float> buf(6 * len + histlen);//X保存水平差分，Y保存豎直差分，Mag保存梯度幅度，Angle保存特征點方向, W保存高斯權重float* X = buf, * Y = buf + len, * Mag = Y, * Angle = Y + len, * W = Angle + len;float* RBin = W + len, * CBin = RBin + len, * hist = CBin + len;//首先清空直方圖數據for (int i = 0; i < d + 2; ++i) // i 對應 row{for (int j = 0; j < d + 2; ++j) // j 對應 col{for (int k = 0; k < n + 2; ++k)hist[(i * (d + 2) + j) * (n + 2) + k] = 0.f;}}//把鄰域內的像素分配到相應子區域內，計算子區域內每個像素點的權重(子區域即 d*d 中每一個小方格)int k = 0;//實際上是在 4 x 4 的網格中找 16 個種子點，每個種子點都在子網格的正中心，//通過三線性插值對不同種子點間的像素點進行加權作用到不同的種子點上繪制方向直方圖for (int i = -radius; i < radius; ++i) // i 對應 y 行坐標{for (int j = -radius; j < radius; ++j) // j 對應 x 列坐標{float c_rot = j * cos_t - i * sin_t; //旋轉后鄰域內采樣點的 x 坐標float r_rot = j * sin_t + i * cos_t; //旋轉后鄰域內采樣點的 y 坐標//旋轉后 5 x 5 的網格中的所有像素點被分配到 4 x 4 的網格中float cbin = c_rot + d / 2 - 0.5f; //旋轉后的采樣點落在子區域的 x 坐標float rbin = r_rot + d / 2 - 0.5f; //旋轉后的采樣點落在子區域的 y 坐標int r = ptxy.y + i, c = ptxy.x + j; //ptxy是高斯金字塔中的坐標//這里rbin，cbin范圍是(-1,d)if (rbin > -1 && rbin < d && cbin>-1 && cbin < d && r>0 && r < rows - 1 && c>0 && c < cols - 1){float dx = gauss_image.at<float>(r, c + 1) - gauss_image.at<float>(r, c - 1);float dy = gauss_image.at<float>(r + 1, c) - gauss_image.at<float>(r - 1, c);X[k] = dx; //鄰域內所有像素點的水平差分Y[k] = dy; //鄰域內所有像素點的豎直差分CBin[k] = cbin; //鄰域內所有采樣點落在子區域的 x 坐標RBin[k] = rbin; //鄰域內所有采樣點落在子區域的 y 坐標W[k] = (c_rot * c_rot + r_rot * r_rot) * exp_scale; //高斯權值的指數部分++k;}}}//計算采樣點落在子區域的像素梯度幅度，梯度角度，和高斯權值len = k;cv::hal::exp(W, W, len); //鄰域內所有采樣點落在子區域的像素的高斯權值cv::hal::fastAtan2(Y, X, Angle, len, true); //鄰域內所有采樣點落在子區域的像素的梯度方向，角度范圍是0-360度cv::hal::magnitude(X, Y, Mag, len); //鄰域內所有采樣點落在子區域的像素的梯度幅度//實際上是在 4 x 4 的網格中找 16 個種子點，每個種子點都在子網格的正中心，//通過三線性插值對不同種子點間的像素點進行加權作用到不同的種子點上繪制方向直方圖//計算每個特征點的特征描述子for (k = 0; k < len; ++k){float rbin = RBin[k], cbin = CBin[k]; //子區域內像素點坐標，rbin,cbin范圍是(-1,d)改進的地方，對方向進行了一個細化，也是為了增加對梯度差異的魯棒性//if (Angle[k] > 180 && Angle[k] < 360)// Angle[k] = 360 - Angle[k];//子區域內像素點處理后的方向float temp = Angle[k] - main_angle;/*if (temp > 180 && temp < 360)temp = 360 - temp;*/float obin = temp * bins_per_rad; //指定方向的數量后，鄰域內像素點對應的方向float mag = Mag[k] * W[k]; //子區域內像素點乘以權值后的梯度幅值int r0 = cvFloor(rbin); //ro取值集合是{-1，0，1，2，3}，沒太懂為什么？int c0 = cvFloor(cbin); //c0取值集合是{-1，0，1，2，3}int o0 = cvFloor(obin);rbin = rbin - r0; //子區域內像素點坐標的小數部分，用于線性插值，分配像素點的作用cbin = cbin - c0;obin = obin - o0; //子區域方向的小數部分//限制范圍為梯度直方圖橫坐標[0,n)，8 個方向直方圖if (o0 < 0)o0 = o0 + n;if (o0 >= n)o0 = o0 - n;//三線性插值用于計算落在兩個子區域之間的像素對兩個子區域的作用，并把其分配到對應子區域的8個方向上//像素對應的信息通過加權分配給其周圍的種子點，并把相應方向的梯度值進行累加 float v_r1 = mag * rbin; //第二行分配的值float v_r0 = mag - v_r1; //第一行分配的值float v_rc11 = v_r1 * cbin; //第二行第二列分配的值，右下角種子點float v_rc10 = v_r1 - v_rc11; //第二行第一列分配的值，左下角種子點float v_rc01 = v_r0 * cbin; //第一行第二列分配的值，右上角種子點float v_rc00 = v_r0 - v_rc01; //第一行第一列分配的值，左上角種子點//一個像素點的方向為每個種子點的兩個方向做出貢獻float v_rco111 = v_rc11 * obin; //右下角種子點第二個方向上分配的值float v_rco110 = v_rc11 - v_rco111; //右下角種子點第一個方向上分配的值float v_rco101 = v_rc10 * obin;float v_rco100 = v_rc10 - v_rco101;float v_rco011 = v_rc01 * obin;float v_rco010 = v_rc01 - v_rco011;float v_rco001 = v_rc00 * obin;float v_rco000 = v_rc00 - v_rco001;//該像素所在網格的索引int idx = ((r0 + 1) * (d + 2) + c0 + 1) * (n + 2) + o0;hist[idx] += v_rco000;hist[idx + 1] += v_rco001;hist[idx + n + 2] += v_rco010;hist[idx + n + 3] += v_rco011;hist[idx + (d + 2) * (n + 2)] += v_rco100;hist[idx + (d + 2) * (n + 2) + 1] += v_rco101;hist[idx + (d + 3) * (n + 2)] += v_rco110;hist[idx + (d + 3) * (n + 2) + 1] += v_rco111;}//由于圓周循環的特性，對計算以后幅角小于 0 度或大于 360 度的值重新進行調整，使//其在 0～360 度之間for (int i = 0; i < d; ++i){for (int j = 0; j < d; ++j){//類似于 hist[0][2][3] 第 0 行，第 2 列，種子點直方圖中的第 3 個 binint idx = ((i + 1) * (d + 2) + (j + 1)) * (n + 2);hist[idx] += hist[idx + n];//hist[idx + 1] += hist[idx + n + 1];//opencv源碼中這句話是多余的,hist[idx + n + 1]永遠是0.0for (k = 0; k < n; ++k)descriptor[(i * d + j) * n + k] = hist[idx + k];}}//對描述子進行歸一化int lenght = d * d * n;float norm = 0;//計算特征描述向量的模值的平方for (int i = 0; i < lenght; ++i){norm = norm + descriptor[i] * descriptor[i];}norm = sqrt(norm); //特征描述向量的模值//此次歸一化能去除光照的影響for (int i = 0; i < lenght; ++i){descriptor[i] = descriptor[i] / norm;}//閾值截斷,去除特征描述向量中大于 0.2 的值，能消除非線性光照的影響(相機飽和度對某些放的梯度影響較大，對方向的影響較小)for (int i = 0; i < lenght; ++i){descriptor[i] = min(descriptor[i], DESCR_MAG_THR);}//再次歸一化，能夠提高特征的獨特性norm = 0;for (int i = 0; i < lenght; ++i){norm = norm + descriptor[i] * descriptor[i];}norm = sqrt(norm);for (int i = 0; i < lenght; ++i){descriptor[i] = descriptor[i] / norm;} }/*************************該函數計算每個特征點的特征描述子*****************************/ /*amplit表示特征點所在層的梯度幅度圖像orient表示特征點所在層的梯度角度圖像pt表示特征點的位置scale表示特征點所在層的尺度main_ori表示特征點的主方向，0-360度d表示GLOH角度方向區間個數，默認是8，n表示每個網格內角度在0-360度之間等分個數，n默認是8*/ static void calc_gloh_descriptor(const Mat& amplit, const Mat& orient, Point2f pt, float scale, float main_ori, int d, int n, float* ptr_des) {Point point(cvRound(pt.x), cvRound(pt.y));//特征點旋轉方向余弦和正弦float cos_t = cosf(-main_ori / 180.f * (float)CV_PI);float sin_t = sinf(-main_ori / 180.f * (float)CV_PI);int num_rows = amplit.rows;int num_cols = amplit.cols;int radius = cvRound(SAR_SIFT_RADIUS_DES * scale);radius = min(radius, min(num_rows / 2, num_cols / 2));//特征點鄰域半徑int radius_x_left = point.x - radius;int radius_x_right = point.x + radius;int radius_y_up = point.y - radius;int radius_y_down = point.y + radius;//防止越界if (radius_x_left < 0)radius_x_left = 0;if (radius_x_right > num_cols - 1)radius_x_right = num_cols - 1;if (radius_y_up < 0)radius_y_up = 0;if (radius_y_down > num_rows - 1)radius_y_down = num_rows - 1;//此時特征點周圍本矩形區域的中心，相對于該矩形int center_x = point.x - radius_x_left;int center_y = point.y - radius_y_up;//特征點周圍區域內像素梯度幅度，梯度角度Mat sub_amplit = amplit(Range(radius_y_up, radius_y_down + 1), Range(radius_x_left, radius_x_right + 1));Mat sub_orient = orient(Range(radius_y_up, radius_y_down + 1), Range(radius_x_left, radius_x_right + 1));//以center_x和center_y位中心，對下面矩形區域進行旋轉Range x_rng(-(point.x - radius_x_left), radius_x_right - point.x);Range y_rng(-(point.y - radius_y_up), radius_y_down - point.y);Mat X, Y;meshgrid(x_rng, y_rng, X, Y);Mat c_rot = X * cos_t - Y * sin_t;Mat r_rot = X * sin_t + Y * cos_t;Mat GLOH_angle, GLOH_amplit;phase(c_rot, r_rot, GLOH_angle, true);//角度在0-360度之間GLOH_amplit = c_rot.mul(c_rot) + r_rot.mul(r_rot);//為了加快速度，沒有計算開方//三個圓半徑平方float R1_pow = (float)radius * radius;//外圓半徑平方float R2_pow = pow(radius * SAR_SIFT_GLOH_RATIO_R1_R2, 2.f);//中間圓半徑平方float R3_pow = pow(radius * SAR_SIFT_GLOH_RATIO_R1_R3, 2.f);//內圓半徑平方int sub_rows = sub_amplit.rows;int sub_cols = sub_amplit.cols;//開始構建描述子,在角度方向對描述子進行插值int len = (d * 2 + 1) * n;AutoBuffer<float> hist(len);for (int i = 0; i < len; ++i)//清零hist[i] = 0;for (int i = 0; i < sub_rows; ++i){float* ptr_amplit = sub_amplit.ptr<float>(i);float* ptr_orient = sub_orient.ptr<float>(i);float* ptr_GLOH_amp = GLOH_amplit.ptr<float>(i);float* ptr_GLOH_ang = GLOH_angle.ptr<float>(i);for (int j = 0; j < sub_cols; ++j){if (((i - center_y) * (i - center_y) + (j - center_x) * (j - center_x)) < radius * radius){float pix_amplit = ptr_amplit[j];//該像素的梯度幅度float pix_orient = ptr_orient[j];//該像素的梯度方向float pix_GLOH_amp = ptr_GLOH_amp[j];//該像素在GLOH網格中的半徑位置float pix_GLOH_ang = ptr_GLOH_ang[j];//該像素在GLOH網格中的位置方向int rbin, cbin, obin;rbin = pix_GLOH_amp < R3_pow ? 0 : (pix_GLOH_amp > R2_pow ? 2 : 1);//rbin={0,1,2}cbin = cvRound(pix_GLOH_ang * d / 360.f);cbin = cbin > d ? cbin - d : (cbin <= 0 ? cbin + d : cbin);//cbin=[1,d]obin = cvRound(pix_orient * n / 360.f);obin = obin > n ? obin - n : (obin <= 0 ? obin + n : obin);//obin=[1,n]if (rbin == 0)//內圓hist[obin - 1] += pix_amplit;else{int idx = ((rbin - 1) * d + cbin - 1) * n + n + obin - 1;hist[idx] += pix_amplit;}}}}//對描述子進行歸一化float norm = 0;for (int i = 0; i < len; ++i){norm = norm + hist[i] * hist[i];}norm = sqrt(norm);for (int i = 0; i < len; ++i){hist[i] = hist[i] / norm;}//閾值截斷for (int i = 0; i < len; ++i){hist[i] = min(hist[i], DESCR_MAG_THR);}//再次歸一化norm = 0;for (int i = 0; i < len; ++i){norm = norm + hist[i] * hist[i];}norm = sqrt(norm);for (int i = 0; i < len; ++i){ptr_des[i] = hist[i] / norm;}}/******************************計算一個特征點的描述子—改進版***********************************/ static void improve_calc_sift_descriptor(const Mat& gauss_image, float main_angle, Point2f pt,float scale, int d, int n, float* descriptor) {int n1 = 16, n2 = 6, n3 = 4;Point ptxy(cvRound(pt.x), cvRound(pt.y)); //坐標取整float cos_t = cosf(-main_angle * (float)(CV_PI / 180)); //計算主方向的余弦float sin_t = sinf(-main_angle * (float)(CV_PI / 180)); //計算主方向的正弦float bin_per_rad_1 = n1 / 360.f; //n=8float bin_per_rad_2 = n2 / 360.f; //原理特征點部分閾值float bin_per_rad_3 = n3 / 360.f; //原理特征點部分閾值float exp_scale = -1.f / (d * d * 0.5f); //權重指數部分float hist_width = DESCR_SCL_FCTR * scale; //子區域邊長，子區域的面積也即采樣像素點個數int radius = cvRound(hist_width * (d + 1) * sqrt(2) * 0.5f);//特征點鄰域半徑(d+1)*(d+1)int rows = gauss_image.rows, cols = gauss_image.cols;//特征點鄰域半徑radius = min(radius, (int)sqrt((double)rows * rows + cols * cols));cos_t = cos_t / hist_width;sin_t = sin_t / hist_width;int len = (2 * radius + 1) * (2 * radius + 1); //鄰域內總像素數int histlen = (d + 2) * (d + 2) * (n1 + 2);AutoBuffer<float> buf(6 * len + histlen);//X保存水平差分，Y保存豎直差分，Mag保存梯度幅度，Angle保存特征點方向, W保存高斯權重float* X = buf, * Y = buf + len, * Mag = Y, * Angle = Y + len, * W = Angle + len;float* X2 = buf, * Y2 = buf + len, * Mag2 = Y, * Angle2 = Y + len, * W2 = Angle + len;float* X3 = buf, * Y3 = buf + len, * Mag3 = Y, * Angle3 = Y + len, * W3 = Angle + len;float* RBin = W + len, * CBin = RBin + len, * hist = CBin + len;float* RBin2 = W + len, * CBin2 = RBin + len;float* RBin3 = W + len, * CBin3 = RBin + len;//首先清空直方圖數據for (int i = 0; i < d + 2; ++i){for (int j = 0; j < d + 2; ++j){for (int k = 0; k < n + 2; ++k)hist[(i * (d + 2) + j) * (n + 2) + k] = 0.f;}}//把鄰域內的像素分配到相應子區域內，計算子區域內每個像素點的權重(子區域即 d*d 中每一個小方格)int k1 = 0, k2 = 0, k3 = 0;vector<int> v; //存放外鄰域像素點對應的序號for (int i = -radius; i < radius; ++i){for (int j = -radius; j < radius; ++j){float c_rot = j * cos_t - i * sin_t; //旋轉后鄰域內采樣點的 x 坐標float r_rot = j * sin_t + i * cos_t; //旋轉后鄰域內采樣點的 y 坐標float rbin = r_rot + d / 2 - 0.5f; //旋轉后的采樣點落在子區域的 y 坐標float cbin = c_rot + d / 2 - 0.5f; //旋轉后的采樣點落在子區域的 x 坐標int r = ptxy.y + i, c = ptxy.x + j; //ptxy是高斯金字塔中的坐標//對離中心點近的部分進行操作if (abs(i) < (radius / 3) && abs(j) < (radius / 3)){//這里rbin,cbin范圍是(-1,d)if (rbin > -1 && rbin < d && cbin>-1 && cbin < d &&r>0 && r < rows - 1 && c>0 && c < cols - 1){float dx = gauss_image.at<float>(r, c + 1) - gauss_image.at<float>(r, c - 1);float dy = gauss_image.at<float>(r + 1, c) - gauss_image.at<float>(r - 1, c);X[k1] = dx; //鄰域內所有像素點的水平差分Y[k1] = dy; //鄰域內所有像素點的豎直差分RBin[k1] = rbin; //鄰域內所有采樣點落在子區域的 y 坐標CBin[k1] = cbin; //鄰域內所有采樣點落在子區域的 x 坐標//高斯權值的指數部分W[k1] = (c_rot * c_rot + r_rot * r_rot) * exp_scale;++k1;}}//對離中心點遠的部分進行操作else if (abs(i) < (2 * radius / 3) && abs(i) > (radius / 3) && abs(j) < (2 * radius / 3) && abs(j) > (radius / 3)){//這里rbin,cbin范圍是(-1,d)if (rbin > -1 && rbin < d && cbin>-1 && cbin < d &&r>0 && r < rows - 1 && c>0 && c < cols - 1){float dx = gauss_image.at<float>(r, c + 1) - gauss_image.at<float>(r, c - 1);float dy = gauss_image.at<float>(r + 1, c) - gauss_image.at<float>(r - 1, c);X2[k2] = dx; //鄰域內所有像素點的水平差分Y2[k2] = dy; //鄰域內所有像素點的豎直差分RBin2[k2] = rbin; //鄰域內所有采樣點落在子區域的 y 坐標CBin2[k2] = cbin; //鄰域內所有采樣點落在子區域的 x 坐標//高斯權值的指數部分W2[k2] = (c_rot * c_rot + r_rot * r_rot) * exp_scale;++k2;}}else{//這里rbin,cbin范圍是(-1,d)if (rbin > -1 && rbin < d && cbin>-1 && cbin < d &&r>0 && r < rows - 1 && c>0 && c < cols - 1){float dx = gauss_image.at<float>(r, c + 1) - gauss_image.at<float>(r, c - 1);float dy = gauss_image.at<float>(r + 1, c) - gauss_image.at<float>(r - 1, c);X3[k3] = dx; //鄰域內所有像素點的水平差分Y3[k3] = dy; //鄰域內所有像素點的豎直差分RBin3[k3] = rbin; //鄰域內所有采樣點落在子區域的 y 坐標CBin3[k3] = cbin; //鄰域內所有采樣點落在子區域的 x 坐標//高斯權值的指數部分W3[k3] = (c_rot * c_rot + r_rot * r_rot) * exp_scale;++k3;}}}}//兩個區域內數組的合并拼接for (int k = 0; k < k2; k++){X[k1 + k] = X2[k];Y[k1 + k] = Y2[k];RBin[k1 + k] = RBin2[k];CBin[k1 + k] = CBin2[k];W[k1 + k] = W2[k];}for (int k = 0; k < k3; k++){X[k1 + k2 + k] = X3[k];Y[k1 + k2 + k] = Y3[k];RBin[k1 + k2 + k] = RBin3[k];CBin[k1 + k2 + k] = CBin3[k];W[k1 + k2 + k] = W3[k];}//計算采樣點落在子區域的像素梯度幅度，梯度角度，和高斯權值len = k1 + k2 + k3;cv::hal::exp(W, W, len); //鄰域內所有采樣點落在子區域的像素的高斯權值cv::hal::fastAtan2(Y, X, Angle, len, true); //鄰域內所有采樣點落在子區域的像素的梯度方向，角度范圍是0-360度cv::hal::magnitude(X, Y, Mag, len); //鄰域內所有采樣點落在子區域的像素的梯度幅度//計算每個特征點的特征描述子for (int k = 0; k < len; ++k){float rbin = RBin[k], cbin = CBin[k]; //子區域內像素點坐標，rbin,cbin范圍是(-1,d)//離特征點進的鄰域if (k < k1){//子區域內像素點處理后的方向float obin = (Angle[k] - main_angle) * bin_per_rad_1;float mag = Mag[k] * W[k]; //子區域內像素點乘以權值后的梯度幅值int r0 = cvFloor(rbin); //ro取值集合是{-1,0,1,2，3}，向下取整int c0 = cvFloor(cbin); //c0取值集合是{-1，0，1，2，3}int o0 = cvFloor(obin);rbin = rbin - r0; //子區域內像素點坐標的小數部分，用于線性插值cbin = cbin - c0;obin = obin - o0;//限制范圍為梯度直方圖橫坐標[0,n)if (o0 < 0)o0 = o0 + n1;if (o0 >= n1)o0 = o0 - n1;//三線性插值用于計算落在兩個子區域之間的像素對兩個子區域的作用，并把其分配到對應子區域的8個方向上//使用三線性插值(即三維)方法，計算直方圖float v_r1 = mag * rbin; //第二行分配的值float v_r0 = mag - v_r1; //第一行分配的值float v_rc11 = v_r1 * cbin; //第二行第二列分配的值float v_rc10 = v_r1 - v_rc11; //第二行第一列分配的值float v_rc01 = v_r0 * cbin; //第一行第二列分配的值float v_rc00 = v_r0 - v_rc01;float v_rco111 = v_rc11 * obin; //第二行第二列第二個方向上分配的值，每個采樣點去鄰近兩個方向float v_rco110 = v_rc11 - v_rco111; //第二行第二列第一個方向上分配的值float v_rco101 = v_rc10 * obin;float v_rco100 = v_rc10 - v_rco101;float v_rco011 = v_rc01 * obin;float v_rco010 = v_rc01 - v_rco011;float v_rco001 = v_rc00 * obin;float v_rco000 = v_rc00 - v_rco001;//該像素所在網格的索引int idx = ((r0 + 1) * (d + 2) + c0 + 1) * (n1 + 2) + o0;hist[idx] += v_rco000;hist[idx + 1] += v_rco001;hist[idx + n1 + 2] += v_rco010;hist[idx + n1 + 3] += v_rco011;hist[idx + (d + 2) * (n1 + 2)] += v_rco100;hist[idx + (d + 2) * (n1 + 2) + 1] += v_rco101;hist[idx + (d + 3) * (n1 + 2)] += v_rco110;hist[idx + (d + 3) * (n1 + 2) + 1] += v_rco111;}//離特征點遠的鄰域else if (k >= k1 && k < k2){//子區域內像素點處理后的方向float obin = (Angle[k] - main_angle) * bin_per_rad_2;float mag = Mag[k] * W[k]; //子區域內像素點乘以權值后的梯度幅值int r0 = cvFloor(rbin); //ro取值集合是{-1,0,1,2，3}，向下取整int c0 = cvFloor(cbin); //c0取值集合是{-1，0，1，2，3}int o0 = cvFloor(obin);rbin = rbin - r0; //子區域內像素點坐標的小數部分，用于線性插值cbin = cbin - c0;obin = obin - o0;//限制范圍為梯度直方圖橫坐標[0,n)if (o0 < 0)o0 = o0 + n2;if (o0 >= n1)o0 = o0 - n2;//三線性插值用于計算落在兩個子區域之間的像素對兩個子區域的作用，并把其分配到對應子區域的8個方向上//使用三線性插值(即三維)方法，計算直方圖float v_r1 = mag * rbin; //第二行分配的值float v_r0 = mag - v_r1; //第一行分配的值float v_rc11 = v_r1 * cbin; //第二行第二列分配的值float v_rc10 = v_r1 - v_rc11; //第二行第一列分配的值float v_rc01 = v_r0 * cbin; //第一行第二列分配的值float v_rc00 = v_r0 - v_rc01;float v_rco111 = v_rc11 * obin; //第二行第二列第二個方向上分配的值，每個采樣點去鄰近兩個方向float v_rco110 = v_rc11 - v_rco111; //第二行第二列第一個方向上分配的值float v_rco101 = v_rc10 * obin;float v_rco100 = v_rc10 - v_rco101;float v_rco011 = v_rc01 * obin;float v_rco010 = v_rc01 - v_rco011;float v_rco001 = v_rc00 * obin;float v_rco000 = v_rc00 - v_rco001;//該像素所在網格的索引int idx = ((r0 + 1) * (d + 2) + c0 + 1) * (n2 + 2) + o0;hist[idx] += v_rco000;hist[idx + 1] += v_rco001;hist[idx + n2 + 2] += v_rco010;hist[idx + n2 + 3] += v_rco011;hist[idx + (d + 2) * (n2 + 2)] += v_rco100;hist[idx + (d + 2) * (n2 + 2) + 1] += v_rco101;hist[idx + (d + 3) * (n2 + 2)] += v_rco110;hist[idx + (d + 3) * (n2 + 2) + 1] += v_rco111;}else{//子區域內像素點處理后的方向float obin = (Angle[k] - main_angle) * bin_per_rad_3;float mag = Mag[k] * W[k]; //子區域內像素點乘以權值后的梯度幅值int r0 = cvFloor(rbin); //ro取值集合是{-1,0,1,2，3}，向下取整int c0 = cvFloor(cbin); //c0取值集合是{-1，0，1，2，3}int o0 = cvFloor(obin);rbin = rbin - r0; //子區域內像素點坐標的小數部分，用于線性插值cbin = cbin - c0;obin = obin - o0;//限制范圍為梯度直方圖橫坐標[0,n)if (o0 < 0)o0 = o0 + n3;if (o0 >= n1)o0 = o0 - n3;//三線性插值用于計算落在兩個子區域之間的像素對兩個子區域的作用，并把其分配到對應子區域的8個方向上//使用三線性插值(即三維)方法，計算直方圖float v_r1 = mag * rbin; //第二行分配的值float v_r0 = mag - v_r1; //第一行分配的值float v_rc11 = v_r1 * cbin; //第二行第二列分配的值float v_rc10 = v_r1 - v_rc11; //第二行第一列分配的值float v_rc01 = v_r0 * cbin; //第一行第二列分配的值float v_rc00 = v_r0 - v_rc01;float v_rco111 = v_rc11 * obin; //第二行第二列第二個方向上分配的值，每個采樣點去鄰近兩個方向float v_rco110 = v_rc11 - v_rco111; //第二行第二列第一個方向上分配的值float v_rco101 = v_rc10 * obin;float v_rco100 = v_rc10 - v_rco101;float v_rco011 = v_rc01 * obin;float v_rco010 = v_rc01 - v_rco011;float v_rco001 = v_rc00 * obin;float v_rco000 = v_rc00 - v_rco001;//該像素所在網格的索引int idx = ((r0 + 1) * (d + 2) + c0 + 1) * (n3 + 2) + o0;hist[idx] += v_rco000;hist[idx + 1] += v_rco001;hist[idx + n3 + 2] += v_rco010;hist[idx + n3 + 3] += v_rco011;hist[idx + (d + 2) * (n3 + 2)] += v_rco100;hist[idx + (d + 2) * (n3 + 2) + 1] += v_rco101;hist[idx + (d + 3) * (n3 + 2)] += v_rco110;hist[idx + (d + 3) * (n3 + 2) + 1] += v_rco111;}}//由于圓周循環的特性，對計算以后幅角小于 0 度或大于 360 度的值重新進行調整，使//其在 0～360 度之間for (int i = 0; i < d; ++i){for (int j = 0; j < d; ++j){int idx = ((i + 1) * (d + 2) + (j + 1)) * (n + 2);hist[idx] += hist[idx + n];//hist[idx + 1] += hist[idx + n + 1];//opencv源碼中這句話是多余的,hist[idx + n + 1]永遠是0.0for (int k = 0; k < n; ++k)descriptor[(i * d + j) * n + k] = hist[idx + k];}}//對描述子進行歸一化int lenght = d * d * n;float norm = 0;//計算特征描述向量的模值的平方for (int i = 0; i < lenght; ++i){norm = norm + descriptor[i] * descriptor[i];}norm = sqrt(norm); //特征描述向量的模值//此次歸一化能去除光照的影響for (int i = 0; i < lenght; ++i){descriptor[i] = descriptor[i] / norm;}//閾值截斷,去除特征描述向量中大于 0.2 的值，能消除非線性光照的影響(相機飽和度對某些放的梯度影響較大，對方向的影響較小)for (int i = 0; i < lenght; ++i){descriptor[i] = min(descriptor[i], DESCR_MAG_THR);}//再次歸一化，能夠提高特征的獨特性norm = 0;for (int i = 0; i < lenght; ++i){norm = norm + descriptor[i] * descriptor[i];}norm = sqrt(norm);for (int i = 0; i < lenght; ++i){descriptor[i] = descriptor[i] / norm;} }/********************************該函數計算所有特征點特征描述子***************************/ /*gauss_pyr表示高斯金字塔keypoints表示特征點、descriptors表示生成的特征點的描述子*/ void mySift::calc_sift_descriptors(const vector<vector<Mat>>& gauss_pyr, vector<KeyPoint>& keypoints,Mat& descriptors, const vector<Mat>& amplit, const vector<Mat>& orient) {int d = DESCR_WIDTH; //d=4,特征點鄰域網格個數是d x dint n = DESCR_HIST_BINS; //n=8,每個網格特征點梯度角度等分為8個方向descriptors.create(keypoints.size(), d * d * n, CV_32FC1); //分配空間for (size_t i = 0; i < keypoints.size(); ++i) //對于每一個特征點{int octaves, layer;//得到特征點所在的組號，層號octaves = keypoints[i].octave & 255;layer = (keypoints[i].octave >> 8) & 255;//得到特征點相對于本組的坐標，不是最底層Point2f pt(keypoints[i].pt.x / (1 << octaves), keypoints[i].pt.y / (1 << octaves));float scale = keypoints[i].size / (1 << octaves); //得到特征點相對于本組的尺度float main_angle = keypoints[i].angle; //特征點主方向//計算該點的描述子calc_sift_descriptor(gauss_pyr[octaves][layer], main_angle, pt, scale, d, n, descriptors.ptr<float>((int)i));//improve_calc_sift_descriptor(gauss_pyr[octaves][layer], main_angle, pt, scale, d, n, descriptors.ptr<float>((int)i));//calc_gloh_descriptor(amplit[octaves], orient[octaves], pt, scale, main_angle, d, n, descriptors.ptr<float>((int)i));if (double_size)//如果圖像尺寸擴大一倍{keypoints[i].pt.x = keypoints[i].pt.x / 2.f;keypoints[i].pt.y = keypoints[i].pt.y / 2.f;}}}/*************該函數構建SAR_SIFT尺度空間*****************/ /*image表示輸入的原始圖像sar_harris_fun表示尺度空間的Sar_harris函數amplit表示尺度空間像素的梯度幅度orient表示尺度空間像素的梯度方向*/ void mySift::build_sar_sift_space(const Mat& image, vector<Mat>& sar_harris_fun, vector<Mat>& amplit, vector<Mat>& orient) {//轉換輸入圖像格式Mat gray_image;if (image.channels() != 1)cvtColor(image, gray_image, CV_RGB2GRAY);elsegray_image = image;//把圖像轉換為0-1之間的浮點數據Mat float_image;double ratio = pow(2, 1.0 / 3.0); //相鄰兩層的尺度比,默認是2^(1/3)//在這里轉換為0-1之間的浮點數據和轉換為0-255之間的浮點數據，效果是一樣的//gray_image.convertTo(float_image, CV_32FC1, 1.f / 255.f, 0.f);//轉換為0-1之間gray_image.convertTo(float_image, CV_32FC1, 1, 0.f); //轉換為0-255之間的浮點數//分配內存sar_harris_fun.resize(Mmax);amplit.resize(Mmax);orient.resize(Mmax);for (int i = 0; i < Mmax; ++i){float scale = (float)sigma * (float)pow(ratio, i); //獲得當前層的尺度int radius = cvRound(2 * scale);Mat kernel;roewa_kernel(radius, scale, kernel);//四個濾波模板生成Mat W34 = Mat::zeros(2 * radius + 1, 2 * radius + 1, CV_32FC1);Mat W12 = Mat::zeros(2 * radius + 1, 2 * radius + 1, CV_32FC1);Mat W14 = Mat::zeros(2 * radius + 1, 2 * radius + 1, CV_32FC1);Mat W23 = Mat::zeros(2 * radius + 1, 2 * radius + 1, CV_32FC1);kernel(Range(radius + 1, 2 * radius + 1), Range::all()).copyTo(W34(Range(radius + 1, 2 * radius + 1), Range::all()));kernel(Range(0, radius), Range::all()).copyTo(W12(Range(0, radius), Range::all()));kernel(Range::all(), Range(radius + 1, 2 * radius + 1)).copyTo(W14(Range::all(), Range(radius + 1, 2 * radius + 1)));kernel(Range::all(), Range(0, radius)).copyTo(W23(Range::all(), Range(0, radius)));//濾波Mat M34, M12, M14, M23;double eps = 0.00001;filter2D(float_image, M34, CV_32FC1, W34, Point(-1, -1), eps);filter2D(float_image, M12, CV_32FC1, W12, Point(-1, -1), eps);filter2D(float_image, M14, CV_32FC1, W14, Point(-1, -1), eps);filter2D(float_image, M23, CV_32FC1, W23, Point(-1, -1), eps);//計算水平梯度和豎直梯度Mat Gx, Gy;log((M14) / (M23), Gx);log((M34) / (M12), Gy);//計算梯度幅度和梯度方向magnitude(Gx, Gy, amplit[i]);phase(Gx, Gy, orient[i], true);//構建sar-Harris矩陣//Mat Csh_11 = log(scale)*log(scale)*Gx.mul(Gx);//Mat Csh_12 = log(scale)*log(scale)*Gx.mul(Gy);//Mat Csh_22 = log(scale)*log(scale)*Gy.mul(Gy);Mat Csh_11 = scale * scale * Gx.mul(Gx);Mat Csh_12 = scale * scale * Gx.mul(Gy);Mat Csh_22 = scale * scale * Gy.mul(Gy);//此時閾值為0.8//Mat Csh_11 = Gx.mul(Gx);//Mat Csh_12 = Gx.mul(Gy);//Mat Csh_22 = Gy.mul(Gy);//此時閾值為0.8/100//高斯權重float gauss_sigma = sqrt(2.f) * scale;int size = cvRound(3 * gauss_sigma);Size kern_size(2 * size + 1, 2 * size + 1);GaussianBlur(Csh_11, Csh_11, kern_size, gauss_sigma, gauss_sigma);GaussianBlur(Csh_12, Csh_12, kern_size, gauss_sigma, gauss_sigma);GaussianBlur(Csh_22, Csh_22, kern_size, gauss_sigma, gauss_sigma);/*Mat gauss_kernel;//自定義圓形高斯核gauss_circle(size, gauss_sigma, gauss_kernel);filter2D(Csh_11, Csh_11, CV_32FC1, gauss_kernel);filter2D(Csh_12, Csh_12, CV_32FC1, gauss_kernel);filter2D(Csh_22, Csh_22, CV_32FC1, gauss_kernel);*/Mat Csh_21 = Csh_12;//構建sar_harris函數Mat temp_add = Csh_11 + Csh_22;double d = 0.04; //sar_haiirs函數表達式中的任意參數，默認是0.04sar_harris_fun[i] = Csh_11.mul(Csh_22) - Csh_21.mul(Csh_12) - (float)d * temp_add.mul(temp_add);} }/***************該函數計算所有特征點的特征向量*************/ /*amplit表示尺度空間像素幅度orient表示尺度空間像素梯度角度keys表示檢測到的特征點descriptors表示特征點描述子向量，【M x N】,M表示描述子個數，N表示描述子維度*/ void mySift::calc_gloh_descriptors(const vector<Mat>& amplit, const vector<Mat>& orient, const vector<KeyPoint>& keys, Mat& descriptors) {int d = SAR_SIFT_GLOH_ANG_GRID; //d=4或者d=8int n = SAR_SIFT_DES_ANG_BINS; //n=8默認int num_keys = (int)keys.size();int grids = 2 * d + 1;//descriptors.create(num_keys, grids * n, CV_32FC1);descriptors.create(num_keys, grids * n, CV_32FC1);for (int i = 0; i < num_keys; ++i){int octaves = keys[i].octave & 255; //特征點所在層float* ptr_des = descriptors.ptr<float>(i);float scale = keys[i].size / (1 << octaves); //得到特征點相對于本組的尺度; //特征點所在層的尺度float main_ori = keys[i].angle; //特征點主方向//得到特征點相對于本組的坐標，不是最底層Point2f point(keys[i].pt.x / (1 << octaves), keys[i].pt.y / (1 << octaves));cout << "layer=" << octaves << endl;cout << "scale=" << scale << endl;//計算該特征點的特征描述子calc_gloh_descriptor(amplit[octaves], orient[octaves], point, scale, main_ori, d, n, ptr_des);} }//特征點檢測和特征點描述把整個 SIFT 算子都涵蓋在內了///******************************特征點檢測*********************************/ /*image表示輸入的圖像gauss_pyr表示生成的高斯金字塔dog_pyr表示生成的高斯差分DOG金字塔keypoints表示檢測到的特征點vector<float>& cell_contrast 用于存放一個單元格中所有特征點的對比度vector<float>& cell_contrasts用于存放一個尺度層中所有單元格中特征點的對比度vector<vector<vector<float>>>& all_cell_contrasts用于存放所有尺度層中所有單元格的對比度vector<vector<float>>& average_contrast用于存放所有尺度層中多有單元格的平均對比度*/ void mySift::detect(const Mat& image, vector<vector<Mat>>& gauss_pyr, vector<vector<Mat>>& dog_pyr, vector<KeyPoint>& keypoints,vector<vector<vector<float>>>& all_cell_contrasts,vector<vector<float>>& average_contrast, vector<vector<int>>& n_cells, vector<int>& num_cell, vector<vector<int>>& available_n,vector<int>& available_num, vector<KeyPoint>& final_keypoints,vector<KeyPoint>& Final_keypoints, vector<KeyPoint>& Final_Keypoints) {if (image.empty() || image.depth() != CV_8U)CV_Error(CV_StsBadArg, "輸入圖像為空，或者圖像深度不是CV_8U");//計算高斯金字塔組數int nOctaves;nOctaves = num_octaves(image);//生成高斯金字塔第一層圖像Mat init_gauss;create_initial_image(image, init_gauss);//生成高斯尺度空間圖像build_gaussian_pyramid(init_gauss, gauss_pyr, nOctaves);//生成高斯差分金字塔(DOG金字塔，or LOG金字塔)build_dog_pyramid(dog_pyr, gauss_pyr);//在DOG金字塔上檢測特征點find_scale_space_extrema1(dog_pyr, gauss_pyr, keypoints); //原始 SIFT 算子// UR_SIFT,僅僅是檢測特征點，并未對不理想的點進行篩選/*UR_SIFT_find_scale_space_extrema(dog_pyr, gauss_pyr, keypoints, all_cell_contrasts,average_contrast, n_cells, num_cell, available_n, available_num, final_keypoints, Final_keypoints, Final_Keypoints, N);*///如果指定了特征點個數,并且設定的設置小于默認檢測的特征點個數if (nfeatures != 0 && nfeatures < (int)keypoints.size()){//特征點響應值(即對比度，對比度越小越不穩定)從大到小排序sort(keypoints.begin(), keypoints.end(), [](const KeyPoint& a, const KeyPoint& b) {return a.response > b.response; });//刪除點多余的特征點keypoints.erase(keypoints.begin() + nfeatures, keypoints.end());} }/**********************特征點描述*******************/ /*gauss_pyr表示高斯金字塔keypoints表示特征點集合descriptors表示特征點的描述子*/ void mySift::comput_des(const vector<vector<Mat>>& gauss_pyr, vector<KeyPoint>& final_keypoints, const vector<Mat>& amplit,const vector<Mat>& orient, Mat& descriptors) {calc_sift_descriptors(gauss_pyr, final_keypoints, descriptors, amplit, orient); }

3、 myMatch.h

#pragma once#include<opencv2\core\core.hpp> #include<opencv2\features2d\features2d.hpp>#include<vector> #include<string> #include<iostream>const double dis_ratio1 = 0.75; //最近鄰和次近鄰距離比閾值，就目前測試來看 dis_ratio = 0.75 時正確匹配的數量相對較多 const double dis_ratio2 = 0.8; const double dis_ratio3 = 0.9;const float ransac_error = 1.5; //ransac算法誤差閾值const double FSC_ratio_low = 0.8;const double FSC_ratio_up = 1;const int pointsCount = 9; // k均值聚類數據點個數 const int clusterCount = 3; // k均值聚類質心的個數enum DIS_CRIT { Euclidean = 0, COS }; //距離度量準則using namespace std; using namespace cv;class myMatch { public:myMatch();~myMatch();//該函數根據正確的匹配點對，計算出圖像之間的變換關系Mat LMS(const Mat& match1_xy, const Mat& match2_xy, string model, float& rmse);//改進版LMS超定方程Mat improve_LMS(const Mat& match1_xy, const Mat& match2_xy, string model, float& rmse);//該函數刪除錯誤的匹配點對Mat ransac(const vector<Point2f>& points_1, const vector<Point2f>& points_2, string model, float threshold, vector<bool>& inliers, float& rmse);//繪制棋盤圖像void mosaic_map(const Mat& image_1, const Mat& image_2, Mat& chessboard_1, Mat& chessboard_2, Mat& mosaic_image, int width);//該函數把兩幅配準后的圖像進行融合鑲嵌void image_fusion(const Mat& image_1, const Mat& image_2, const Mat T, Mat& fusion_image, Mat& matched_image);//該函數進行描述子的最近鄰和次近鄰匹配void match_des(const Mat& des_1, const Mat& des_2, vector<vector<DMatch>>& dmatchs, DIS_CRIT dis_crite);//建立尺度直方圖、ROM 直方圖void scale_ROM_Histogram(const vector<DMatch>& matches, float* scale_hist, float* ROM_hist, int n);//該函數刪除錯誤匹配點對，并進行配準Mat match(const Mat& image_1, const Mat& image_2, const vector<vector<DMatch>>& dmatchs, vector<KeyPoint> keys_1,vector<KeyPoint> keys_2, string model, vector<DMatch>& right_matchs, Mat& matched_line, vector<DMatch>& init_matchs); };

4、 myMatch.cpp

#include"myMatch.h"#include<opencv2\imgproc\types_c.h> #include<opencv2\imgproc\imgproc.hpp> #include<opencv2\highgui\highgui.hpp> #include<opencv2\features2d\features2d.hpp> //特征提取#include<algorithm> #include<vector> #include<cmath> #include<opencv.hpp>#include <numeric> //用于容器元素求和using namespace std; using namespace cv; //using namespace gpu;RNG rng(100);myMatch::myMatch() { }/********該函數根據正確的匹配點對，計算出圖像之間的變換關系********/ /*注意：輸入幾個點都能計算對應的 x 矩陣，(2N,8)*(8,1)=(2N,1)match1_xy表示參考圖像特征點坐標集合,[M x 2]矩陣，M表示特征的個數match2_xy表示待配準圖像特征點集合，[M x 2]矩陣，M表示特征點集合model表示變換類型，“相似變換”,"仿射變換","透視變換"rmse表示均方根誤差返回值為計算得到的3 x 3矩陣參數*/ Mat myMatch::LMS(const Mat& match1_xy, const Mat& match2_xy, string model, float& rmse) {if (match1_xy.rows != match2_xy.rows)CV_Error(CV_StsBadArg, "LMS模塊輸入特征點對個數不一致！");if (!(model == string("affine") || model == string("similarity") ||model == string("perspective") || model == string("projective")))CV_Error(CV_StsBadArg, "LMS模塊圖像變換類型輸入錯誤！");const int N = match1_xy.rows; //特征點個數Mat match2_xy_trans, match1_xy_trans; //特征點坐標轉置transpose(match1_xy, match1_xy_trans); //矩陣轉置(2,M)transpose(match2_xy, match2_xy_trans);Mat change = Mat::zeros(3, 3, CV_32FC1); //變換矩陣//A*X=B,接下來部分仿射變換和透視變換一樣,如果特征點個數是M，則A=[2*M,6]矩陣//A=[x1,y1,0,0,1,0;0,0,x1,y1,0,1;.....xn,yn,0,0,1,0;0,0,xn,yn,0,1]，應該是改版過的Mat A = Mat::zeros(2 * N, 6, CV_32FC1);for (int i = 0; i < N; ++i){A.at<float>(2 * i, 0) = match2_xy.at<float>(i, 0);//xA.at<float>(2 * i, 1) = match2_xy.at<float>(i, 1);//yA.at<float>(2 * i, 4) = 1.f;A.at<float>(2 * i + 1, 2) = match2_xy.at<float>(i, 0);A.at<float>(2 * i + 1, 3) = match2_xy.at<float>(i, 1);A.at<float>(2 * i + 1, 5) = 1.f;}//如果特征點個數是M,那個B=[2*M,1]矩陣//B=[u1,v1,u2,v2,.....,un,vn]Mat B;B.create(2 * N, 1, CV_32FC1);for (int i = 0; i < N; ++i){B.at<float>(2 * i, 0) = match1_xy.at<float>(i, 0); //xB.at<float>(2 * i + 1, 0) = match1_xy.at<float>(i, 1);//y}//如果是仿射變換if (model == string("affine")){Vec6f values;solve(A, B, values, DECOMP_QR);change = (Mat_<float>(3, 3) << values(0), values(1), values(4),values(2), values(3), values(5),+0.0f, +0.0f, 1.0f);Mat temp_1 = change(Range(0, 2), Range(0, 2)); //尺度和旋轉量Mat temp_2 = change(Range(0, 2), Range(2, 3)); //平移量Mat match2_xy_change = temp_1 * match2_xy_trans + repeat(temp_2, 1, N);Mat diff = match2_xy_change - match1_xy_trans; //求差pow(diff, 2.f, diff);rmse = (float)sqrt(sum(diff)(0) * 1.0) / N; //sum輸出是各個通道的和, / N 初始實在括號里面}//如果是透視變換else if (model == string("perspective")){/*透視變換模型[u'*w,v'*w, w]'=[u,v,w]' = [a1, a2, a5;a3, a4, a6;a7, a8, 1] * [x, y, 1]'[u',v']'=[x,y,0,0,1,0,-u'x, -u'y;0, 0, x, y, 0, 1, -v'x,-v'y] * [a1, a2, a3, a4, a5, a6, a7, a8]'即，Y = A*X *///構造 A 矩陣的后兩列Mat A2;A2.create(2 * N, 2, CV_32FC1);for (int i = 0; i < N; ++i){A2.at<float>(2 * i, 0) = match1_xy.at<float>(i, 0) * match2_xy.at<float>(i, 0) * (-1.f);A2.at<float>(2 * i, 1) = match1_xy.at<float>(i, 0) * match2_xy.at<float>(i, 1) * (-1.f);A2.at<float>(2 * i + 1, 0) = match1_xy.at<float>(i, 1) * match2_xy.at<float>(i, 0) * (-1.f);A2.at<float>(2 * i + 1, 1) = match1_xy.at<float>(i, 1) * match2_xy.at<float>(i, 1) * (-1.f);}Mat A1; //完成的 A 矩陣，(8,8)A1.create(2 * N, 8, CV_32FC1);A.copyTo(A1(Range::all(), Range(0, 6)));A2.copyTo(A1(Range::all(), Range(6, 8)));Mat values; //values中存放的是待求的8個參數solve(A1, B, values, DECOMP_QR); //(2N,8)*(8,1)=(2N,1)好像本身就有點超定矩陣的意思change.at<float>(0, 0) = values.at<float>(0);change.at<float>(0, 1) = values.at<float>(1);change.at<float>(0, 2) = values.at<float>(4);change.at<float>(1, 0) = values.at<float>(2);change.at<float>(1, 1) = values.at<float>(3);change.at<float>(1, 2) = values.at<float>(5);change.at<float>(2, 0) = values.at<float>(6);change.at<float>(2, 1) = values.at<float>(7);change.at<float>(2, 2) = 1.f;Mat temp1 = Mat::ones(1, N, CV_32FC1);Mat temp2;temp2.create(3, N, CV_32FC1);match2_xy_trans.copyTo(temp2(Range(0, 2), Range::all()));temp1.copyTo(temp2(Range(2, 3), Range::all()));Mat match2_xy_change = change * temp2; //待配準圖像中的特征點在參考圖中的映射結果Mat match2_xy_change_12 = match2_xy_change(Range(0, 2), Range::all());//float* temp_ptr = match2_xy_change.ptr<float>(2);float* temp_ptr = match2_xy_change.ptr<float>(2);for (int i = 0; i < N; ++i){float div_temp = temp_ptr[i];match2_xy_change_12.at<float>(0, i) = match2_xy_change_12.at<float>(0, i) / div_temp;match2_xy_change_12.at<float>(1, i) = match2_xy_change_12.at<float>(1, i) / div_temp;}Mat diff = match2_xy_change_12 - match1_xy_trans;pow(diff, 2, diff);rmse = (float)sqrt(sum(diff)(0) * 1.0) / N; //sum輸出是各個通道的和，rmse為輸入點的均方誤差}//如果是相似變換else if (model == string("similarity")){/*[x, y, 1, 0;y, -x, 0, 1] * [a, b, c, d]'=[u,v]*/Mat A3;A3.create(2 * N, 4, CV_32FC1);for (int i = 0; i < N; ++i){A3.at<float>(2 * i, 0) = match2_xy.at<float>(i, 0);A3.at<float>(2 * i, 1) = match2_xy.at<float>(i, 1);A3.at<float>(2 * i, 2) = 1.f;A3.at<float>(2 * i, 3) = 0.f;A3.at<float>(2 * i + 1, 0) = match2_xy.at<float>(i, 1);A3.at<float>(2 * i + 1, 1) = match2_xy.at<float>(i, 0) * (-1.f);A3.at<float>(2 * i + 1, 2) = 0.f;A3.at<float>(2 * i + 1, 3) = 1.f;}Vec4f values;solve(A3, B, values, DECOMP_QR);change = (Mat_<float>(3, 3) << values(0), values(1), values(2),values(1) * (-1.0f), values(0), values(3),+0.f, +0.f, 1.f);Mat temp_1 = change(Range(0, 2), Range(0, 2));//尺度和旋轉量Mat temp_2 = change(Range(0, 2), Range(2, 3));//平移量Mat match2_xy_change = temp_1 * match2_xy_trans + repeat(temp_2, 1, N);Mat diff = match2_xy_change - match1_xy_trans;//求差pow(diff, 2, diff);rmse = (float)sqrt(sum(diff)(0) * 1.0) / N;//sum輸出是各個通道的和}return change; }/********改進版LMS超定方程********/ Mat myMatch::improve_LMS(const Mat& match1_xy, const Mat& match2_xy, string model, float& rmse) {if (match1_xy.rows != match2_xy.rows)CV_Error(CV_StsBadArg, "LMS模塊輸入特征點對個數不一致！");if (!(model == string("affine") || model == string("similarity") ||model == string("perspective")))CV_Error(CV_StsBadArg, "LMS模塊圖像變換類型輸入錯誤！");const int N = match1_xy.rows; //特征點個數Mat match2_xy_trans, match1_xy_trans; //特征點坐標轉置transpose(match1_xy, match1_xy_trans); //矩陣轉置(2,M)transpose(match2_xy, match2_xy_trans);Mat change = Mat::zeros(3, 3, CV_32FC1); //變換矩陣//A*X=B,接下來部分仿射變換和透視變換一樣,如果特征點個數是M，則A=[2*M,6]矩陣//A=[x1,y1,0,0,1,0;0,0,x1,y1,0,1;.....xn,yn,0,0,1,0;0,0,xn,yn,0,1]，應該是改版過的Mat A = Mat::zeros(2 * N, 6, CV_32FC1);for (int i = 0; i < N; ++i){A.at<float>(2 * i, 0) = match2_xy.at<float>(i, 0);//xA.at<float>(2 * i, 1) = match2_xy.at<float>(i, 1);//yA.at<float>(2 * i, 4) = 1.f;//A.at<float>(2 * i, 0) = match2_xy.at<float>(i, 0);//x//A.at<float>(2 * i, 1) = match2_xy.at<float>(i, 1);//y//A.at<float>(2 * i, 2) = 1.f;A.at<float>(2 * i + 1, 2) = match2_xy.at<float>(i, 0);A.at<float>(2 * i + 1, 3) = match2_xy.at<float>(i, 1);A.at<float>(2 * i + 1, 5) = 1.f;/*A.at<float>(2 * i + 1, 3) = match2_xy.at<float>(i, 0);A.at<float>(2 * i + 1, 4) = match2_xy.at<float>(i, 1);A.at<float>(2 * i + 1, 5) = 1.f;*/}//如果特征點個數是M,那個B=[2*M,1]矩陣//B=[u1,v1,u2,v2,.....,un,vn]Mat B;B.create(2 * N, 1, CV_32FC1);for (int i = 0; i < N; ++i){B.at<float>(2 * i, 0) = match1_xy.at<float>(i, 0); //xB.at<float>(2 * i + 1, 0) = match1_xy.at<float>(i, 1);//y}//如果是仿射變換if (model == string("affine")){Vec6f values;solve(A, B, values, DECOMP_QR);change = (Mat_<float>(3, 3) << values(0), values(1), values(4),values(2), values(3), values(5),+0.0f, +0.0f, 1.0f);Mat temp_1 = change(Range(0, 2), Range(0, 2));//尺度和旋轉量Mat temp_2 = change(Range(0, 2), Range(2, 3));//平移量Mat match2_xy_change = temp_1 * match2_xy_trans + repeat(temp_2, 1, N);Mat diff = match2_xy_change - match1_xy_trans;//求差pow(diff, 2.f, diff);rmse = (float)sqrt(sum(diff)(0) * 1.0 / N);//sum輸出是各個通道的和}//如果是透視變換else if (model == string("perspective")){/*透視變換模型[u'*w,v'*w, w]'=[u,v,w]' = [a1, a2, a5;a3, a4, a6;a7, a8, 1] * [x, y, 1]'[u',v']'=[x,y,0,0,1,0,-u'x, -u'y;0, 0, x, y, 0, 1, -v'x,-v'y] * [a1, a2, a3, a4, a5, a6, a7, a8]'即，Y = A*X *///構造 A 矩陣的后兩列Mat A2;A2.create(2 * N, 2, CV_32FC1);for (int i = 0; i < N; ++i){A2.at<float>(2 * i, 0) = match1_xy.at<float>(i, 0) * match2_xy.at<float>(i, 0) * (-1.f);A2.at<float>(2 * i, 1) = match1_xy.at<float>(i, 0) * match2_xy.at<float>(i, 1) * (-1.f);A2.at<float>(2 * i + 1, 0) = match1_xy.at<float>(i, 1) * match2_xy.at<float>(i, 0) * (-1.f);A2.at<float>(2 * i + 1, 1) = match1_xy.at<float>(i, 1) * match2_xy.at<float>(i, 1) * (-1.f);}Mat A1; //完整的 A 矩陣，(8,8)A1.create(2 * N, 8, CV_32FC1);A.copyTo(A1(Range::all(), Range(0, 6)));A2.copyTo(A1(Range::all(), Range(6, 8)));Mat AA1, balance;Mat evects, evals;transpose(A1, AA1); //求矩陣 A1 的轉置balance = AA1 * A1;double a[8][8];for (int i = 0; i < 8; i++){for (int j = 0; j < 8; j++){a[i][j] = balance.at<float>(i, j);}}//構造輸入矩陣CvMat SrcMatrix = cvMat(8, 8, CV_32FC1, a);double b[8][8] ={{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};// 構造輸出特征向量矩陣,特征向量按行存儲,且與特征值相對應CvMat ProVector = cvMat(8, 8, CV_32FC1, b);// 構造輸出特征值矩陣,特征值按降序配列double c[8] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };CvMat ProValue = cvMat(8, 1, CV_32FC1, c);//求特征向量，最后一行特征向量即對應的 H 矩陣的參數cvEigenVV(&SrcMatrix, &ProVector, &ProValue, 1.0e-6F);輸出特征向量矩陣//for (int i = 0; i < 2; i++)//{// for (int j = 0; j < 2; j++)// printf("%f\t", cvmGet(&ProVector, i, j));// printf("\n");//}//把計算得到的最小特征值對應的特征變量賦給 H 矩陣change.at<float>(0, 0) = cvmGet(&ProVector, 7, 0);change.at<float>(0, 1) = cvmGet(&ProVector, 7, 1);change.at<float>(0, 2) = cvmGet(&ProVector, 7, 2);change.at<float>(1, 0) = cvmGet(&ProVector, 7, 3);change.at<float>(1, 1) = cvmGet(&ProVector, 7, 4);change.at<float>(1, 2) = cvmGet(&ProVector, 7, 5);change.at<float>(2, 0) = cvmGet(&ProVector, 7, 6);change.at<float>(2, 1) = cvmGet(&ProVector, 7, 7);change.at<float>(2, 2) = 1.f;Mat temp1 = Mat::ones(1, N, CV_32FC1);Mat temp2; //存放處理后的特征點(x,y,1)Ttemp2.create(3, N, CV_32FC1);match2_xy_trans.copyTo(temp2(Range(0, 2), Range::all()));temp1.copyTo(temp2(Range(2, 3), Range::all()));Mat match2_xy_change = change * temp2; //待配準圖像中的特征點在參考圖中的映射結果Mat match2_xy_change_12 = match2_xy_change(Range(0, 2), Range::all());//float* temp_ptr = match2_xy_change.ptr<float>(2);float* temp_ptr = match2_xy_change.ptr<float>(2);for (int i = 0; i < N; ++i){float div_temp = temp_ptr[i];match2_xy_change_12.at<float>(0, i) = match2_xy_change_12.at<float>(0, i) / div_temp;match2_xy_change_12.at<float>(1, i) = match2_xy_change_12.at<float>(1, i) / div_temp;}Mat diff = match2_xy_change_12 - match1_xy_trans;pow(diff, 2, diff);rmse = (float)sqrt(sum(diff)(0) * 1.0 / N); //sum輸出是各個通道的和，rmse為輸入點的均方誤差}//如果是相似變換else if (model == string("similarity")){/*[x, y, 1, 0;y, -x, 0, 1] * [a, b, c, d]'=[u,v]*/Mat A3;A3.create(2 * N, 4, CV_32FC1);for (int i = 0; i < N; ++i){A3.at<float>(2 * i, 0) = match2_xy.at<float>(i, 0);A3.at<float>(2 * i, 1) = match2_xy.at<float>(i, 1);A3.at<float>(2 * i, 2) = 1.f;A3.at<float>(2 * i, 3) = 0.f;A3.at<float>(2 * i + 1, 0) = match2_xy.at<float>(i, 1);A3.at<float>(2 * i + 1, 1) = match2_xy.at<float>(i, 0) * (-1.f);A3.at<float>(2 * i + 1, 2) = 0.f;A3.at<float>(2 * i + 1, 3) = 1.f;}Vec4f values;solve(A3, B, values, DECOMP_QR);change = (Mat_<float>(3, 3) << values(0), values(1), values(2),values(1) * (-1.0f), values(0), values(3),+0.f, +0.f, 1.f);Mat temp_1 = change(Range(0, 2), Range(0, 2));//尺度和旋轉量Mat temp_2 = change(Range(0, 2), Range(2, 3));//平移量Mat match2_xy_change = temp_1 * match2_xy_trans + repeat(temp_2, 1, N);Mat diff = match2_xy_change - match1_xy_trans;//求差pow(diff, 2, diff);rmse = (float)sqrt(sum(diff)(0) * 1.0) / N;//sum輸出是各個通道的和}return change; }/*********************該函數刪除錯誤的匹配點對****************************/ /*points_1表示參考圖像上匹配的特征點位置points_2表示待配準圖像上匹配的特征點位置model表示變換模型，“similarity”,"affine"，“perspective”threshold表示內點閾值inliers表示points_1和points_2中對應的點對是否是正確匹配，如果是，對應元素值為1，否則為0rmse表示最后所有正確匹配點對計算出來的誤差返回一個3 x 3矩陣，表示待配準圖像到參考圖像的變換矩陣*/ Mat myMatch::ransac(const vector<Point2f>& points_1, const vector<Point2f>& points_2, string model, float threshold, vector<bool>& inliers, float& rmse) {if (points_1.size() != points_2.size())CV_Error(CV_StsBadArg, "ransac模塊輸入特征點數量不一致！");if (!(model == string("affine") || model == string("similarity") ||model == string("perspective") || model == string("projective")))CV_Error(CV_StsBadArg, "ransac模塊圖像變換類型輸入錯誤！");const size_t N = points_1.size(); //特征點對數，size_t 類型常用來保存一個數據的大小，通常為整形，可移植性強int n; //相當于不同模型對應的標簽size_t max_iteration, iterations;//確定最大迭代次數，就目前來看此法過于簡單粗暴，可以使用自適應迭代次數法if (model == string("similarity")){n = 2;max_iteration = N * (N - 1) / 2;}else if (model == string("affine")){n = 3;max_iteration = N * (N - 1) * (N - 2) / (2 * 3);}else if (model == string("perspective")){n = 4;max_iteration = N * (N - 1) * (N - 2) / (2 * 3);}if (max_iteration > 800)iterations = 800;elseiterations = max_iteration;//取出保存在points_1和points_2中的點坐標，保存在Mat矩陣中，方便處理Mat arr_1, arr_2; //arr_1,和arr_2是一個[3 x N]的矩陣，每一列表示一個點坐標,第三行全是1arr_1.create(3, N, CV_32FC1);arr_2.create(3, N, CV_32FC1);//獲取矩陣每一行的首地址float* p10 = arr_1.ptr<float>(0), * p11 = arr_1.ptr<float>(1), * p12 = arr_1.ptr<float>(2);float* p20 = arr_2.ptr<float>(0), * p21 = arr_2.ptr<float>(1), * p22 = arr_2.ptr<float>(2);//把特征點放到矩陣中for (size_t i = 0; i < N; ++i){p10[i] = points_1[i].x;p11[i] = points_1[i].y;p12[i] = 1.f;p20[i] = points_2[i].x;p21[i] = points_2[i].y;p22[i] = 1.f;}Mat rand_mat; //特征點索引rand_mat.create(1, n, CV_32SC1);int* p = rand_mat.ptr<int>(0);Mat sub_arr1, sub_arr2; //存放隨機挑選的特征點sub_arr1.create(n, 2, CV_32FC1);sub_arr2.create(n, 2, CV_32FC1);Mat T; //待配準圖像到參考圖像的變換矩陣int most_consensus_num = 0; //當前最優一致集個數初始化為0vector<bool> right;right.resize(N);inliers.resize(N);for (size_t i = 0; i < iterations; ++i) //迭代次數{//隨機選擇n個不同的點對，不同的模型每次隨機選擇的個數不同while (1){randu(rand_mat, 0, N - 1); //隨機生成n個范圍在[0,N-1]之間的數，作為獲取特征點的索引//保證這n個點坐標不相同if (n == 2 && p[0] != p[1] &&(p10[p[0]] != p10[p[1]] || p11[p[0]] != p11[p[1]]) &&(p20[p[0]] != p20[p[1]] || p21[p[0]] != p21[p[1]]))break;if (n == 3 && p[0] != p[1] && p[0] != p[2] && p[1] != p[2] &&(p10[p[0]] != p10[p[1]] || p11[p[0]] != p11[p[1]]) &&(p10[p[0]] != p10[p[2]] || p11[p[0]] != p11[p[2]]) &&(p10[p[1]] != p10[p[2]] || p11[p[1]] != p11[p[2]]) &&(p20[p[0]] != p20[p[1]] || p21[p[0]] != p21[p[1]]) &&(p20[p[0]] != p20[p[2]] || p21[p[0]] != p21[p[2]]) &&(p20[p[1]] != p20[p[2]] || p21[p[1]] != p21[p[2]]))break;if (n == 4 && p[0] != p[1] && p[0] != p[2] && p[0] != p[3] &&p[1] != p[2] && p[1] != p[3] && p[2] != p[3] &&(p10[p[0]] != p10[p[1]] || p11[p[0]] != p11[p[1]]) &&(p10[p[0]] != p10[p[2]] || p11[p[0]] != p11[p[2]]) &&(p10[p[0]] != p10[p[3]] || p11[p[0]] != p11[p[3]]) &&(p10[p[1]] != p10[p[2]] || p11[p[1]] != p11[p[2]]) &&(p10[p[1]] != p10[p[3]] || p11[p[1]] != p11[p[3]]) &&(p10[p[2]] != p10[p[3]] || p11[p[2]] != p11[p[3]]) &&(p20[p[0]] != p20[p[1]] || p21[p[0]] != p21[p[1]]) &&(p20[p[0]] != p20[p[2]] || p21[p[0]] != p21[p[2]]) &&(p20[p[0]] != p20[p[3]] || p21[p[0]] != p21[p[3]]) &&(p20[p[1]] != p20[p[2]] || p21[p[1]] != p21[p[2]]) &&(p20[p[1]] != p20[p[3]] || p21[p[1]] != p21[p[3]]) &&(p20[p[2]] != p20[p[3]] || p21[p[2]] != p21[p[3]]))break;}//提取出n個點對for (int i = 0; i < n; ++i){sub_arr1.at<float>(i, 0) = p10[p[i]];sub_arr1.at<float>(i, 1) = p11[p[i]];sub_arr2.at<float>(i, 0) = p20[p[i]];sub_arr2.at<float>(i, 1) = p21[p[i]];}//根據這n個點對，計算初始變換矩陣 TT = LMS(sub_arr1, sub_arr2, model, rmse);int consensus_num = 0; //當前一致集(內點)個數if (model == string("perspective")){//變換矩陣計算待配準圖像特征點在參考圖像中的映射點Mat match2_xy_change = T * arr_2; //arr_2 中存放的是待配準圖像的特征點 (3,N),//match2_xy_change(Range(0, 2), Range::all())意思是提取 match2_xy_change 的 0、1 行，所有的列Mat match2_xy_change_12 = match2_xy_change(Range(0, 2), Range::all());//獲取 match2_xy_change 第二行首地址float* temp_ptr = match2_xy_change.ptr<float>(2);for (size_t i = 0; i < N; ++i){float div_temp = temp_ptr[i]; //match2_xy_change 第二行第 i 列值，除以 div_temp ，是為了保證第三行為 1，和原始坐標相對應match2_xy_change_12.at<float>(0, i) = match2_xy_change_12.at<float>(0, i) / div_temp;match2_xy_change_12.at<float>(1, i) = match2_xy_change_12.at<float>(1, i) / div_temp;}//計算待配準圖像特征點在參考圖像中的映射點與參考圖像中對應點的距離Mat diff = match2_xy_change_12 - arr_1(Range(0, 2), Range::all());pow(diff, 2, diff);//第一行和第二行求和，即兩點間距離的平方Mat add = diff(Range(0, 1), Range::all()) + diff(Range(1, 2), Range::all());float* p_add = add.ptr<float>(0);//遍歷所有距離，如果小于閾值，則認為是局內點for (size_t i = 0; i < N; ++i){if (p_add[i] < threshold) //初始 p_add[i] {right[i] = true;++consensus_num;}elseright[i] = false;}}else if (model == string("affine") || model == string("similarity")){Mat match2_xy_change = T * arr_2; //計算在參考圖像中的映射坐標Mat diff = match2_xy_change - arr_1;pow(diff, 2, diff);//第一行和第二行求和，計算特征點間的距離的平方Mat add = diff(Range(0, 1), Range::all()) + diff(Range(1, 2), Range::all());float* p_add = add.ptr<float>(0);for (size_t i = 0; i < N; ++i){//如果小于閾值if (p_add[i] < threshold){right[i] = true;++consensus_num;}elseright[i] = false;}}//判斷當前一致集是否是優于之前最優一致集，并更新當前最優一致集個數if (consensus_num > most_consensus_num){most_consensus_num = consensus_num;//把正確匹配的點賦予標簽 1 for (size_t i = 0; i < N; ++i)inliers[i] = right[i];}}//刪除重復點對for (size_t i = 0; i < N - 1; ++i){for (size_t j = i + 1; j < N; ++j){if (inliers[i] && inliers[j]){if (p10[i] == p10[j] && p11[i] == p11[j] && p20[i] == p20[j] && p21[i] == p21[j]){inliers[j] = false;--most_consensus_num;}}}}//迭代結束，獲得最優一致集合，根據這些最優一致集合計算出最終的變換關系 TMat consensus_arr1, consensus_arr2; //經過迭代后最終確認正確匹配的點consensus_arr1.create(most_consensus_num, 2, CV_32FC1);consensus_arr2.create(most_consensus_num, 2, CV_32FC1);int k = 0;for (size_t i = 0; i < N; ++i){if (inliers[i]){consensus_arr1.at<float>(k, 0) = p10[i];consensus_arr1.at<float>(k, 1) = p11[i];consensus_arr2.at<float>(k, 0) = p20[i];consensus_arr2.at<float>(k, 1) = p21[i];++k;}}int num_ransac = (model == string("similarity") ? 2 : (model == string("affine") ? 3 : 4));if (k < num_ransac)CV_Error(CV_StsBadArg, "ransac模塊刪除錯誤點對后剩下正確點對個數不足以計算出變換關系矩陣！");//利用迭代后正確匹配點計算變換矩陣，為什么不是挑選 n 個點計算變換矩陣T = LMS(consensus_arr1, consensus_arr2, model, rmse);return T; }/********************該函數生成兩幅圖的棋盤網格圖*************************/ /*image_1表示參考圖像image_2表示配準后的待配準圖像chessboard_1表示image_1的棋盤圖像chessboard_2表示image_2的棋盤圖像mosaic_image表示image_1和image_2的鑲嵌圖像width表示棋盤網格大小*/ void myMatch::mosaic_map(const Mat& image_1, const Mat& image_2, Mat& chessboard_1, Mat& chessboard_2, Mat& mosaic_image, int width) {if (image_1.size != image_2.size)CV_Error(CV_StsBadArg, "mosaic_map模塊輸入兩幅圖大小必須一致！");//生成image_1的棋盤網格圖chessboard_1 = image_1.clone();int rows_1 = chessboard_1.rows;int cols_1 = chessboard_1.cols;int row_grids_1 = cvFloor((double)rows_1 / width); //行方向網格個數int col_grids_1 = cvFloor((double)cols_1 / width); //列方向網格個數//指定區域像素賦值為零，便形成了棋盤圖//第一幅圖，第一行 2、4、6 像素值賦值零；第一幅圖與第二幅圖零像素位置交叉，以便兩幅圖交叉顯示for (int i = 0; i < row_grids_1; i = i + 2){for (int j = 1; j < col_grids_1; j = j + 2){Range range_x(j * width, (j + 1) * width);Range range_y(i * width, (i + 1) * width);chessboard_1(range_y, range_x) = 0;}}for (int i = 1; i < row_grids_1; i = i + 2){for (int j = 0; j < col_grids_1; j = j + 2){Range range_x(j * width, (j + 1) * width);Range range_y(i * width, (i + 1) * width);chessboard_1(range_y, range_x) = 0;}}//生成image_2的棋盤網格圖chessboard_2 = image_2.clone();int rows_2 = chessboard_2.rows;int cols_2 = chessboard_2.cols;int row_grids_2 = cvFloor((double)rows_2 / width);//行方向網格個數int col_grids_2 = cvFloor((double)cols_2 / width);//列方向網格個數//第二幅圖，第一行 1、3、5 像素值賦值零for (int i = 0; i < row_grids_2; i = i + 2){for (int j = 0; j < col_grids_2; j = j + 2){Range range_x(j * width, (j + 1) * width);Range range_y(i * width, (i + 1) * width);chessboard_2(range_y, range_x) = 0;}}for (int i = 1; i < row_grids_2; i = i + 2){for (int j = 1; j < col_grids_2; j = j + 2){Range range_x(j * width, (j + 1) * width);Range range_y(i * width, (i + 1) * width);chessboard_2(range_y, range_x) = 0;}}//兩個棋盤圖進行疊加，顯示配準效果mosaic_image = chessboard_1 + chessboard_2; }/*該函數對輸入圖像指定位置像素進行中值濾波，消除邊緣拼接陰影*/ /*image表示輸入的圖像position表示需要進行中值濾波的位置*/ inline void median_filter(Mat& image, const vector<vector<int>>& pos) {int channels = image.channels();switch (channels){case 1://單通道for (auto beg = pos.cbegin(); beg != pos.cend(); ++beg){int i = (*beg)[0];//yint j = (*beg)[1];//xuchar& pix_val = image.at<uchar>(i, j);vector<uchar> pixs;for (int row = -1; row <= 1; ++row){for (int col = -1; col <= 1; ++col){if (i + row >= 0 && i + row < image.rows && j + col >= 0 && j + col < image.cols){pixs.push_back(image.at<uchar>(i + row, j + col));}}}//排序std::sort(pixs.begin(), pixs.end());pix_val = pixs[pixs.size() / 2];}break;case 3://3通道for (auto beg = pos.cbegin(); beg != pos.cend(); ++beg){int i = (*beg)[0];//yint j = (*beg)[1];//xVec3b& pix_val = image.at<Vec3b>(i, j);vector<cv::Vec3b> pixs;for (int row = -1; row <= 1; ++row){for (int col = -1; col <= 1; ++col){if (i + row >= 0 && i + row < image.rows && j + col >= 0 && j + col < image.cols){pixs.push_back(image.at<Vec3b>(i + row, j + col));}}}//排序std::sort(pixs.begin(), pixs.end(),[pix_val](const Vec3b& a, const Vec3b& b)->bool {return sum((a).ddot(a))[0] < sum((b).ddot(b))[0];});pix_val = pixs[pixs.size() / 2];}break;default:break;} }/***************該函數把配準后的圖像進行融合*****************/ /*該函數功能主要是來對圖像進行融合，以顯示配準的效果*image_1表示參考圖像*image_2表示待配準圖像*T表示待配準圖像到參考圖像的轉換矩陣*fusion_image表示參考圖像和待配準圖像融合后的圖像*mosaic_image表示參考圖像和待配準圖像融合鑲嵌后的圖像，鑲嵌圖形是為了觀察匹配效果*matched_image表示把待配準圖像進行配準后的結果*/ void myMatch::image_fusion(const Mat& image_1, const Mat& image_2, const Mat T, Mat& fusion_image, Mat& matched_image) {//有關depth()的理解，詳解：https://blog.csdn.net/datouniao1/article/details/113524784if (!(image_1.depth() == CV_8U && image_2.depth() == CV_8U))CV_Error(CV_StsBadArg, "image_fusion模塊僅支持uchar類型圖像！");if (image_1.channels() == 4 || image_2.channels() == 4)CV_Error(CV_StsBadArg, "image_fusion模塊僅僅支持單通道或者3通道圖像");int rows_1 = image_1.rows, cols_1 = image_1.cols;int rows_2 = image_2.rows, cols_2 = image_2.cols;int channel_1 = image_1.channels();int channel_2 = image_2.channels();//可以對：彩色-彩色、彩色-灰色、灰色-彩色、灰色-灰色的配準Mat image_1_temp, image_2_temp;if (channel_1 == 3 && channel_2 == 3){image_1_temp = image_1;image_2_temp = image_2;}else if (channel_1 == 1 && channel_2 == 3){image_1_temp = image_1;cvtColor(image_2, image_2_temp, CV_RGB2GRAY); //顏色空間轉換，把彩色圖轉化為灰度圖}else if (channel_1 == 3 && channel_2 == 1){cvtColor(image_1, image_1_temp, CV_RGB2GRAY);image_2_temp = image_2;}else if (channel_1 == 1 && channel_2 == 1){image_1_temp = image_1;image_2_temp = image_2;}//創建一個（3，3）float 矩陣 Mat_ 是一個模板類Mat T_temp = (Mat_<float>(3, 3) << 1, 0, cols_1, 0, 1, rows_1, 0, 0, 1);Mat T_1 = T_temp * T;//對參考圖像和待配準圖像進行變換Mat trans_1, trans_2;//same type as image_2_temp trans_1 = Mat::zeros(3 * rows_1, 3 * cols_1, image_1_temp.type()); //創建擴大后的矩陣image_1_temp.copyTo(trans_1(Range(rows_1, 2 * rows_1), Range(cols_1, 2 * cols_1))); //把image_1_temp中的數據復制到擴大后矩陣的對應位置warpPerspective(image_2_temp, trans_2, T_1, Size(3 * cols_1, 3 * rows_1), INTER_LINEAR, 0, Scalar::all(0));/*功能：把image_2_temp投射到一個新的視平面，即變形*image_2_temp為輸入矩陣*trans_2為輸出矩陣，尺寸和輸入矩陣大小一致*T_1為變換矩陣，(3, 3)矩陣用于透視變換, (2, 2)用于線性變換, (1, 1)用于平移*warpPerspective函數功能是進行透視變換，T_1是(3,3)的透視變換矩陣*int flags=INTER_LINEAR為輸出圖像的插值方法*int borderMode=BORDER_CONSTANT，0 為圖像邊界的填充方式*const Scalar& borderValue=Scalar()：邊界的顏色設置，一般默認是0，Scalar::all(0)對影像邊界外進行填充*///使用簡單的均值法進行圖像融合Mat trans = trans_2.clone(); //把經過透視變換的image_2復制給transint nRows = rows_1;int nCols = cols_1;int len = nCols;bool flag_1 = false;bool flag_2 = false;vector<vector<int>> positions; //保存邊緣位置坐標switch (image_1_temp.channels()){case 1: //如果圖像1是單通道的for (int i = 0; i < nRows; ++i){uchar* ptr_1 = trans_1.ptr<uchar>(i + rows_1); //訪問trans_1中的指定行像素值uchar* ptr = trans.ptr<uchar>(i + rows_1); //訪問trans 中的指定行像素值for (int j = 0; j < nCols; ++j){if (ptr[j + len] == 0 && ptr_1[j + len] != 0) //非重合區域{flag_1 = true;if (flag_2) //表明從重合區域過度到了非重合區域{for (int p = -1; p <= 1; ++p) //保存邊界3x3區域像素{for (int q = -1; q <= 1; ++q){vector<int> pos;pos.push_back(i + rows_1 + p);pos.push_back(j + cols_1 + q);positions.push_back(pos);//保存邊緣位置坐標}}flag_2 = false;}ptr[j + len] = ptr_1[j + len];}else//對于重合區域{flag_2 = true;if (flag_1)//表明從非重合區域過度到了重合區域{for (int p = -1; p <= 1; ++p)//保存邊界3x3區域像素{for (int q = -1; q <= 1; ++q){vector<int> pos;pos.push_back(i + rows_1 + p);pos.push_back(j + cols_1 + q);positions.push_back(pos);//保存邊緣位置坐標}}flag_1 = false;}ptr[j + len] = saturate_cast<uchar>(((float)ptr[j + len] + (float)ptr_1[j + len]) / 2);}}}break;case 3: //如果圖像是三通道的len = len * image_1_temp.channels(); //3倍的列數for (int i = 0; i < nRows; ++i){uchar* ptr_1 = trans_1.ptr<uchar>(i + rows_1); //訪問trans_1中的指定行像素值uchar* ptr = trans.ptr<uchar>(i + rows_1); //訪問trans 中的指定行像素值for (int j = 0; j < nCols; ++j){int nj = j * image_1_temp.channels();//若兩張影像對應列像素值不同（3通道），則是非重合區，該過程僅僅是為了使配準后的影像進行融合，而非配準if (ptr[nj + len] == 0 && ptr[nj + len + 1] == 0 && ptr[nj + len + 2] == 0 &&ptr_1[nj + len] != 0 && ptr_1[nj + len + 1] != 0 && ptr_1[nj + len + 2] != 0){flag_1 = true;if (flag_2) //表明從重合區域過度到了非重合區域{for (int p = -1; p <= 1; ++p) //保存邊界3x3區域像素{for (int q = -1; q <= 1; ++q){vector<int> pos;pos.push_back(i + rows_1 + p);pos.push_back(j + cols_1 + q);positions.push_back(pos); //保存邊緣位置坐標}}flag_2 = false;}ptr[nj + len] = ptr_1[nj + len];ptr[nj + len + 1] = ptr_1[nj + len + 1];ptr[nj + len + 2] = ptr_1[nj + len + 2];}else{ //對于重合區域flag_2 = true;if (flag_1) //表明從非重合區域過度到了重合區域{for (int p = -1; p <= 1; ++p) //保存邊界3x3區域像素{for (int q = -1; q <= 1; ++q){vector<int> pos;pos.push_back(i + rows_1 + p);pos.push_back(j + cols_1 + q);positions.push_back(pos); //保存邊緣位置坐標}}flag_1 = false;}ptr[nj + len] = saturate_cast<uchar>(((float)ptr[nj + len] + (float)ptr_1[nj + len]) / 2);ptr[nj + len + 1] = saturate_cast<uchar>(((float)ptr[nj + len + 1] + (float)ptr_1[nj + len + 1]) / 2);ptr[nj + len + 2] = saturate_cast<uchar>(((float)ptr[nj + len + 2] + (float)ptr_1[nj + len + 2]) / 2);}}}break;default:break;}//根據獲取的邊緣區域的坐標，對邊緣像素進行中值濾波，消除邊緣效應median_filter(trans, positions);//刪除多余的區域Mat left_up = T_1 * (Mat_<float>(3, 1) << 0, 0, 1); //左上角Mat left_down = T_1 * (Mat_<float>(3, 1) << 0, rows_2 - 1, 1); //左下角Mat right_up = T_1 * (Mat_<float>(3, 1) << cols_2 - 1, 0, 1); //右上角Mat right_down = T_1 * (Mat_<float>(3, 1) << cols_2 - 1, rows_2 - 1, 1); //右下角//對于透視變換，需要除以一個因子left_up = left_up / left_up.at<float>(2, 0);left_down = left_down / left_down.at<float>(2, 0);right_up = right_up / right_up.at<float>(2, 0);right_down = right_down / right_down.at<float>(2, 0);//計算x,y坐標的范圍float temp_1 = min(left_up.at<float>(0, 0), left_down.at<float>(0, 0));float temp_2 = min(right_up.at<float>(0, 0), right_down.at<float>(0, 0));float min_x = min(temp_1, temp_2);temp_1 = max(left_up.at<float>(0, 0), left_down.at<float>(0, 0));temp_2 = max(right_up.at<float>(0, 0), right_down.at<float>(0, 0));float max_x = max(temp_1, temp_2);temp_1 = min(left_up.at<float>(1, 0), left_down.at<float>(1, 0));temp_2 = min(right_up.at<float>(1, 0), right_down.at<float>(1, 0));float min_y = min(temp_1, temp_2);temp_1 = max(left_up.at<float>(1, 0), left_down.at<float>(1, 0));temp_2 = max(right_up.at<float>(1, 0), right_down.at<float>(1, 0));float max_y = max(temp_1, temp_2);int X_min = max(cvFloor(min_x), 0);int X_max = min(cvCeil(max_x), 3 * cols_1 - 1);int Y_min = max(cvFloor(min_y), 0);int Y_max = min(cvCeil(max_y), 3 * rows_1 - 1);if (X_min > cols_1)X_min = cols_1;if (X_max < 2 * cols_1 - 1)X_max = 2 * cols_1 - 1;if (Y_min > rows_1)Y_min = rows_1;if (Y_max < 2 * rows_1 - 1)Y_max = 2 * rows_1 - 1;//提取有價值區域Range Y_range(Y_min, Y_max + 1);Range X_range(X_min, X_max + 1);fusion_image = trans(Y_range, X_range);matched_image = trans_2(Y_range, X_range);Mat ref_matched = trans_1(Y_range, X_range);//生成棋盤網格圖像/*Mat chessboard_1, chessboard_2;mosaic_map(trans_1(Y_range, X_range), trans_2(Y_range, X_range), chessboard_1, chessboard_2, mosaic_image, 100);*//*cv::imwrite(".\\image_save\\參考圖像棋盤圖像.jpg", chessboard_1);cv::imwrite(".\\image_save\\待配準圖像棋盤圖像.jpg", chessboard_2);*/cv::imwrite(".\\image_save\\配準后的參考圖像.jpg", ref_matched);cv::imwrite(".\\image_save\\配準后的待配準圖像.jpg", matched_image); }/******該函數計算參考圖像一個描述子和待配準圖像所有描述子的歐式距離，并獲得最近鄰和次近鄰距離，以及對應的索引*//*sub_des_1表示參考圖像的一個描述子*des_2表示待配準圖像描述子*num_des_2值待配準圖像描述子個數*dims_des指的是描述子維度*dis保存最近鄰和次近鄰距離*idx保存最近鄰和次近鄰索引*/ inline void min_dis_idx(const float* ptr_1, const Mat& des_2, int num_des2, int dims_des, float dis[2], int idx[2]) {float min_dis1 = 1000, min_dis2 = 2000;int min_idx1, min_idx2;for (int j = 0; j < num_des2; ++j){const float* ptr_des_2 = des_2.ptr<float>(j);float cur_dis = 0;for (int k = 0; k < dims_des; ++k){float diff = ptr_1[k] - ptr_des_2[k];cur_dis += diff * diff;}if (cur_dis < min_dis1) {min_dis1 = cur_dis;min_idx1 = j;}else if (cur_dis >= min_dis1 && cur_dis < min_dis2) {min_dis2 = cur_dis;min_idx2 = j;}}dis[0] = sqrt(min_dis1); dis[1] = sqrt(min_dis2);idx[0] = min_idx1; idx[1] = min_idx2; }/*加速版本的描述子匹配，返回匹配點候選集函數,該加速版本比前面版本速度提升了3倍*/ void myMatch::match_des(const Mat& des_1, const Mat& des_2, vector<vector<DMatch>>& dmatchs, DIS_CRIT dis_crite) {int num_des_1 = des_1.rows;int num_des_2 = des_2.rows;int dims_des = des_1.cols;vector<DMatch> match(2);//對于參考圖像上的每一點，和待配準圖像進行匹配if (dis_crite == 0) //歐幾里得距離{for (int i = 0; i < num_des_1; ++i) //對于參考圖像中的每個描述子{const float* ptr_des_1 = des_1.ptr<float>(i);float dis[2];int idx[2];min_dis_idx(ptr_des_1, des_2, num_des_2, dims_des, dis, idx);match[0].queryIdx = i;match[0].trainIdx = idx[0];match[0].distance = dis[0];match[1].queryIdx = i;match[1].trainIdx = idx[1];match[1].distance = dis[1];dmatchs.push_back(match);}}else if (dis_crite == 1)//cos距離{Mat trans_des2;transpose(des_2, trans_des2);double aa = (double)getTickCount();Mat mul_des = des_1 * trans_des2;//gemm(des_1, des_2, 1, Mat(), 0, mul_des, GEMM_2_T);double time1 = ((double)getTickCount() - aa) / getTickFrequency();cout << "cos距離中矩陣乘法花費時間： " << time1 << "s" << endl;for (int i = 0; i < num_des_1; ++i){float max_cos1 = -1000, max_cos2 = -2000;int max_idx1, max_idx2;float* ptr_1 = mul_des.ptr<float>(i);for (int j = 0; j < num_des_2; ++j){float cur_cos = ptr_1[j];if (cur_cos > max_cos1) {max_cos1 = cur_cos;max_idx1 = j;}else if (cur_cos <= max_cos1 && cur_cos > max_cos2) {max_cos2 = cur_cos;max_idx2 = j;}}match[0].queryIdx = i;match[0].trainIdx = max_idx1;match[0].distance = acosf(max_cos1);match[1].queryIdx = i;match[1].trainIdx = max_idx2;match[1].distance = acosf(max_cos2);dmatchs.push_back(match);}} }/*******************建立尺度直方圖、ROM 直方圖************************/ void myMatch::scale_ROM_Histogram(const vector<DMatch>& matches, float* scale_hist, float* ROM_hist, int n) {int len = matches.size();//使用AutoBuffer分配一段內存，這里多出4個空間的目的是為了方便后面平滑直方圖的需要AutoBuffer<float> buffer((4 * len) + n + 4);//X保存水平差分，Y保存數值差分，Mag保存梯度幅度，Ori保存梯度角度，W保存高斯權重float* X = buffer, * Y = buffer + len, * Mag = Y, * Ori = Y + len, * W = Ori + len;float* temp_hist = W + len + 2; //臨時保存直方圖數據for (int i = 0; i < n; ++i)temp_hist[i] = 0.f; //數據清零 }/*******************該函數刪除錯誤匹配點對,并完成匹配************************/ /*image_1表示參考圖像，image_2表示待配準圖像dmatchs表示最近鄰和次近鄰匹配點對keys_1表示參考圖像特征點集合keys_2表示待配準圖像特征點集合model表示變換模型right_matchs表示參考圖像和待配準圖像正確匹配點對matched_line表示在參考圖像和待配準圖像上繪制連接線該函數返回變換模型參數*/ Mat myMatch::match(const Mat& image_1, const Mat& image_2, const vector<vector<DMatch>>& dmatchs, vector<KeyPoint> keys_1,vector<KeyPoint> keys_2, string model, vector<DMatch>& right_matchs, Mat& matched_line, vector<DMatch>& init_matchs) {//獲取初始匹配的關鍵點的位置vector<Point2f> point_1, point_2;for (size_t i = 0; i < dmatchs.size(); ++i) //增加一個一個0.8，正確匹配點數增加2{double dis_1 = dmatchs[i][0].distance; //distance對應的是特征點描述符的歐式距離double dis_2 = dmatchs[i][1].distance;//比率測試篩選誤匹配點(初步篩選)，如果滿足則認為是有候選集中有正確匹配點if ((dis_1 / dis_2) < dis_ratio3) //最近鄰和次近鄰距離比閾值{//queryIdx、trainIdx和distance是DMatch類中的一些屬性 //pt是KeyPoint類中的成員，對應關鍵點的坐標point_1.push_back(keys_1[dmatchs[i][0].queryIdx].pt); //queryIdx對應的是特征描述子的下標，也是對應特征點的下標point_2.push_back(keys_2[dmatchs[i][0].trainIdx].pt); //trainIdx對應的是特征描述子的下標，也是對應特征點的下標init_matchs.push_back(dmatchs[i][0]); //保存正確的dmatchs}}cout << "距離比之后初始匹配點對個數是： " << init_matchs.size() << endl;int min_pairs = (model == string("similarity") ? 2 : (model == string("affine") ? 3 : 4));if (point_1.size() < min_pairs)CV_Error(CV_StsBadArg, "match模塊距離比階段匹配特征點個數不足！");//使用ransac算法再次對匹配點進行篩選，然后使用最后確定的匹配點計算變換矩陣的參數vector<bool> inliers; //存放的是 bool 類型的數據，對應特征點float rmse;//homography是一個(3,3)矩陣，是待配準影像到參考影像的變換矩陣，初始誤差閾值是 1.5Mat homography = ransac(point_1, point_2, model, 1.5, inliers, rmse);//提取出處正確匹配點對int i = 0;vector<Point2f> point_11, point_22;vector<DMatch>::iterator itt = init_matchs.begin();for (vector<bool>::iterator it = inliers.begin(); it != inliers.end(); ++it, ++itt){if (*it) //如果是正確匹配點對{right_matchs.push_back(*itt);// init_matchs 中匹配點對的存儲順序和 point_1 中特征點的存儲順序是一一對應的point_11.push_back(point_1.at(i));point_22.push_back(point_2.at(i));}++i;}cout << "使用RANSAC刪除錯誤點對,且返回正確匹配個數： " << right_matchs.size() << endl;cout << "誤差rmse: " << rmse << endl;//繪制初始匹配點對連線圖,此時初始匹配指的是經過 KnnMatch 篩選后的匹配Mat initial_matched; //輸出矩陣，類似畫布drawMatches(image_1, keys_1, image_2, keys_2, init_matchs, initial_matched,Scalar(255, 0, 255), Scalar(0, 255, 0), vector<char>()); //該函數用于繪制特征點并對匹配的特征點進行連線imwrite(".\\image_save\\初始匹配點對.jpg", initial_matched); //保存圖片,第一個顏色控制連線，第二個顏色控制特征點//繪制正確匹配點對連線圖drawMatches(image_1, keys_1, image_2, keys_2, right_matchs, matched_line,Scalar(255, 0, 255), Scalar(0, 255, 0), vector<char>());imwrite(".\\image_save\\正確匹配點對.jpg", matched_line);//保存和顯示檢測到的特征點Mat keys_image_1, keys_image_2; //輸出矩陣，類似畫布drawKeypoints(image_1, keys_1, keys_image_1, Scalar(0, 255, 0)); //該函數用于繪制圖像中的特征點drawKeypoints(image_2, keys_2, keys_image_2, Scalar(0, 255, 0));imwrite(".\\image_save\\參考圖像檢測到的特征點.jpg", keys_image_1);imwrite(".\\image_save\\待配準圖像檢測到的.jpg", keys_image_2);return homography; }myMatch::~myMatch() { }

5、 myDisplay.h

#pragma#include<iostream> #include<opencv2\core\core.hpp> #include<opencv2\features2d\features2d.hpp> #include<vector>using namespace cv; using namespace std;class myDisplay { public:myDisplay() {}void mosaic_pyramid(const vector<vector<Mat>>& pyramid, Mat& pyramid_image, int nOctaceLayers, string str);void write_mosaic_pyramid(const vector<vector<Mat>>& gauss_pyr_1, const vector<vector<Mat>>& dog_pyr_1,const vector<vector<Mat>>& gauss_pyr_2, const vector<vector<Mat>>& dog_pyr_2, int nOctaveLayers);//沒用到，后文沒有對其進行定義void write_keys_image(vector<KeyPoint>& keypoints_1, vector<KeyPoint>& keypoints_2,const Mat& image_1, const Mat& image_2, Mat& image_1_keys, Mat& image_2_keys, bool double_size); };

6、 myDisplay.cpp

#include "myDisplay.h"#include<opencv2\highgui\highgui.hpp> #include<vector> #include<sstream>/***************************該函數把金字塔放在一副圖像上(包絡高斯金字塔和DOG金字塔)******************/ /*pyramid表示高斯金字塔或者DOG金字塔pyramid_image表示生成的金字塔圖像nOctaveLayers表示每組中間層數，默認是3str表示是高斯金字塔還是DOG金字塔 */ void myDisplay::mosaic_pyramid(const vector<vector<Mat>>& pyramid, Mat& pyramid_image, int nOctaceLayers, string str) {//獲得每組金字塔圖像大小vector<vector<int>> temp_size;for (auto beg = pyramid.cbegin(); beg != pyramid.cend(); ++beg){vector<int> temp_vec;int rows = (*beg)[0].rows;int cols = (*beg)[0].cols;temp_vec.push_back(rows);temp_vec.push_back(cols);temp_size.push_back(temp_vec);}//計算最后生成的金字塔圖像pyramid_image的大小int total_rows = 0, total_cols = 0;for (auto beg = temp_size.begin(); beg != temp_size.end(); ++beg){total_rows = total_rows + (*beg)[0];//獲取行大小if (beg == temp_size.begin()) {if (str == string("高斯金字塔"))total_cols = (nOctaceLayers + 3) * ((*beg)[1]);//獲取列大小else if (str == string("DOG金字塔"))total_cols = (nOctaceLayers + 2) * ((*beg)[1]);//獲取列大小}}pyramid_image.create(total_rows, total_cols, CV_8UC1);int i = 0, accumulate_row = 0;for (auto beg = pyramid.cbegin(); beg != pyramid.cend(); ++beg, ++i){int accumulate_col = 0;accumulate_row = accumulate_row + temp_size[i][0];for (auto it = (*beg).cbegin(); it != (*beg).cend(); ++it){accumulate_col = accumulate_col + temp_size[i][1];Mat temp(pyramid_image, Rect(accumulate_col - temp_size[i][1], accumulate_row - temp_size[i][0], it->cols, it->rows));Mat temp_it;normalize(*it, temp_it, 0, 255, NORM_MINMAX, CV_32FC1);convertScaleAbs(temp_it, temp_it, 1, 0);temp_it.copyTo(temp);}} }/**************************該函數保存拼接后的高斯金字塔和DOG金字塔圖像**************************/ /*gauss_pyr_1表示參考高斯金字塔dog_pyr_1表示參考DOG金字塔gauss_pyr_2表示待配準高斯金字塔dog_pyr_2表示待配準DOG金字塔nOctaveLayers表示金字塔每組中間層數*/ void myDisplay::write_mosaic_pyramid(const vector<vector<Mat>>& gauss_pyr_1, const vector<vector<Mat>>& dog_pyr_1,const vector<vector<Mat>>& gauss_pyr_2, const vector<vector<Mat>>& dog_pyr_2, int nOctaveLayers) {//顯示參考和待配準高斯金字塔圖像Mat gauss_image_1, gauss_image_2;mosaic_pyramid(gauss_pyr_1, gauss_image_1, nOctaveLayers, string("高斯金字塔"));mosaic_pyramid(gauss_pyr_2, gauss_image_2, nOctaveLayers, string("高斯金字塔"));imwrite(".\\image_save\\參考圖像高斯金字塔.jpg", gauss_image_1);imwrite(".\\image_save\\待配準圖像高斯金字塔.jpg", gauss_image_2);//顯示參考和待配準DOG金字塔圖像Mat dog_image_1, dog_image_2;mosaic_pyramid(dog_pyr_1, dog_image_1, nOctaveLayers, string("DOG金字塔"));mosaic_pyramid(dog_pyr_2, dog_image_2, nOctaveLayers, string("DOG金字塔"));imwrite(".\\image_save\\參考圖像DOG金字塔.jpg", dog_image_1);imwrite(".\\image_save\\待配準圖像DOG金字塔.jpg", dog_image_2); }

7、?main.cpp

#include"mySift.h" #include"myDisplay.h" #include"myMatch.h" #include<direct.h>#include<opencv2\highgui\highgui.hpp> #include<opencv2\calib3d\calib3d.hpp> #include<opencv2\imgproc\imgproc.hpp>#include<fstream> #include<stdlib.h> #include<direct.h>int main(int argc, char* argv[]) {/********************** 1、讀入數據 **********************/Mat image_1 = imread("..\\test_data\\perspective_graf_1.ppm", -1);Mat image_2 = imread("..\\test_data\\perspective_graf_2.ppm", -1);string change_model = "perspective"; //affine為仿射變換，初始為perspectivestring folderPath = ".\\image_save";_mkdir(folderPath.c_str());double total_count_beg = (double)getTickCount(); //算法運行總時間開始計時mySift sift_1(0, 3, 0.04, 10, 1.6, true); //類對象/********************** 1、參考圖像特征點檢測和描述 **********************/vector<vector<Mat>> gauss_pyr_1, dog_pyr_1; //高斯金字塔和高斯差分金字塔vector<KeyPoint> keypoints_1; //特征點vector<vector<vector<float>>> all_cell_contrasts_1; //所有尺度層中所有單元格的對比度vector<vector<float>> average_contrast_1; //所有尺度層中多有單元格的平均對比度vector<vector<int>> n_cells_1; //所有尺度層，每一尺度層中所有單元格所需特征數vector<int> num_cell_1; //當前尺度層中所有單元格中的所需特征數vector<vector<int>> available_n_1; //所有尺度層，每一尺度層中所有單元格可得到特征數vector<int> available_num_1; //一個尺度層中所有單元格中可用特征數量vector<KeyPoint> final_keypoints1; //第一次篩選結果vector<KeyPoint> Final_keypoints1; //第二次篩選結果vector<KeyPoint> Final_Keypoints1; //第三次篩選結果Mat descriptors_1; //特征描述子double detect_1 = (double)getTickCount(); //特征點檢測運行時間開始計時sift_1.detect(image_1, gauss_pyr_1, dog_pyr_1, keypoints_1, all_cell_contrasts_1,average_contrast_1, n_cells_1, num_cell_1, available_n_1, available_num_1, final_keypoints1, Final_keypoints1, Final_Keypoints1); //特征點檢測cout << "參考圖像檢測出的總特征數 =" << keypoints_1.size() << endl;//getTickFrequency() 返回值為一秒的計時周期數，二者比值為特征點檢測時間double detect_time_1 = ((double)getTickCount() - detect_1) / getTickFrequency();cout << "參考圖像特征點檢測時間是： " << detect_time_1 << "s" << endl;double comput_1 = (double)getTickCount();vector<Mat> sar_harris_fun_1;vector<Mat> amplit_1;vector<Mat> orient_1;sift_1.comput_des(gauss_pyr_1, keypoints_1, amplit_1, orient_1, descriptors_1);double comput_time_1 = ((double)getTickCount() - comput_1) / getTickFrequency();cout << "參考圖像特征點描述時間是： " << comput_time_1 << "s" << endl;/********************** 1、待配準圖像特征點檢測和描述 **********************/vector<vector<Mat>> gauss_pyr_2, dog_pyr_2;vector<KeyPoint> keypoints_2;vector<vector<vector<float>>> all_cell_contrasts_2; //所有尺度層中所有單元格的對比度vector<vector<float>> average_contrast_2; //所有尺度層中多有單元格的平均對比度vector<vector<int>> n_cells_2; //所有尺度層，每一尺度層中所有單元格所需特征數vector<int> num_cell_2; //當前尺度層中所有單元格中的所需特征數vector<vector<int>> available_n_2; //所有尺度層，每一尺度層中的一個單元格可得到特征數vector<int> available_num_2; //一個尺度層中所有單元格中可用特征數量vector<KeyPoint> final_keypoints2; //第一次篩選結果vector<KeyPoint> Final_keypoints2; //第二次篩選結果vector<KeyPoint> Final_Keypoints2; //第三次篩選結果Mat descriptors_2;double detect_2 = (double)getTickCount();sift_1.detect(image_2, gauss_pyr_2, dog_pyr_2, keypoints_2, all_cell_contrasts_2,average_contrast_2, n_cells_2, num_cell_2, available_n_2, available_num_2, final_keypoints2, Final_keypoints2, Final_Keypoints2);cout << "待配準圖像檢測出的總特征數 =" << keypoints_2.size() << endl;double detect_time_2 = ((double)getTickCount() - detect_2) / getTickFrequency();cout << "待配準圖像特征點檢測時間是： " << detect_time_2 << "s" << endl;double comput_2 = (double)getTickCount();vector<Mat> sar_harris_fun_2;vector<Mat> amplit_2;vector<Mat> orient_2;sift_1.comput_des(gauss_pyr_2, keypoints_2, amplit_2, orient_2, descriptors_2);double comput_time_2 = ((double)getTickCount() - comput_2) / getTickFrequency();cout << "待配準特征點描述時間是： " << comput_time_2 << "s" << endl;/********************** 1、最近鄰與次近鄰距離比匹配，兩幅影像進行配準 **********************/myMatch mymatch;double match_time = (double)getTickCount(); //影像配準計時開始//knnMatch函數是DescriptorMatcher類的成員函數，FlannBasedMatcher是DescriptorMatcher的子類Ptr<DescriptorMatcher> matcher1 = new FlannBasedMatcher;Ptr<DescriptorMatcher> matcher2 = new FlannBasedMatcher;vector<vector<DMatch>> dmatchs; //vector<DMatch>中存放的是一個描述子可能匹配的候選集vector<vector<DMatch>> dmatch1;vector<vector<DMatch>> dmatch2;matcher1->knnMatch(descriptors_1, descriptors_2, dmatchs, 2);cout << "距離比之前初始匹配點對個數是：" << dmatchs.size() << endl;Mat matched_lines; //同名點連線vector<DMatch> init_matchs, right_matchs; //用于存放正確匹配的點//該函數返回的是空間映射模型參數Mat homography = mymatch.match(image_1, image_2, dmatchs, keypoints_1, keypoints_2, change_model, right_matchs, matched_lines, init_matchs);double match_time_2 = ((double)getTickCount() - match_time) / getTickFrequency();cout << "特征點匹配花費時間是： " << match_time_2 << "s" << endl;cout << change_model << "變換矩陣是：" << endl;cout << homography << endl;/********************** 1、把正確匹配點坐標寫入文件中 **********************/ofstream ofile;ofile.open(".\\position.txt"); //創建文件for (size_t i = 0; i < right_matchs.size(); ++i){ofile << keypoints_1[right_matchs[i].queryIdx].pt << " "<< keypoints_2[right_matchs[i].trainIdx].pt << endl;}/********************** 1、圖像融合 **********************/double fusion_beg = (double)getTickCount();Mat fusion_image, mosaic_image, regist_image;mymatch.image_fusion(image_1, image_2, homography, fusion_image, regist_image);imwrite(".\\image_save\\融合后的圖像.jpg", fusion_image);double fusion_time = ((double)getTickCount() - fusion_beg) / getTickFrequency();cout << "圖像融合花費時間是： " << fusion_time << "s" << endl;double total_time = ((double)getTickCount() - total_count_beg) / getTickFrequency();cout << "總花費時間是： " << total_time << "s" << endl;/********************** 1、生成棋盤圖，顯示配準結果 **********************///生成棋盤網格圖像Mat chessboard_1, chessboard_2, mosaic_images;Mat image1 = imread(".\\image_save\\配準后的參考圖像.jpg", -1);Mat image2 = imread(".\\image_save\\配準后的待配準圖像.jpg", -1);mymatch.mosaic_map(image1, image2, chessboard_1, chessboard_2, mosaic_images, 50);imwrite(".\\image_save\\參考圖像棋盤圖像.jpg", chessboard_1);imwrite(".\\image_save\\待配準圖像棋盤圖像.jpg", chessboard_2);imwrite(".\\image_save\\兩幅棋盤圖像疊加.jpg", mosaic_images);return 0; }

四、SIFT 算法實驗結果

1、實驗數據

說明：該數據是從網上下載的兩個時相影像，左側為舊時相影像，右側為新時相影像；?

2、實驗結果?

(1)特征點檢測結果

（2）最終匹配結果

??

（3）配準效果?

說明：把兩張圖像交叉重疊在一起，可以觀察其配準效果；?

?五、配置說明

????????該結果是在 VS2019 + opencv4.4.0 + Debug x64 環境下跑出來的；

??如果感覺復制麻煩，完整工程見：SIFT算法C++代碼實現-C/C++文檔類資源-CSDN下載

總結

以上是生活随笔為你收集整理的SIFT算法详解(附有完整代码)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。