特征匹配部分由ORB篇已介紹OPENCV中特征匹配需要用到的一些函數和類的封裝完成，本篇不再介紹。SIFT和SURF由於版權問題，（SIFT在2020年（今年）3月6日專利有限期20年過期，OPENCV后續的版本中可能會有相應接口。）在opencv4.1+中沒有函數接口，可通過安裝對應版本的opencv-contrib進行使用。
1 背景
SIFT（Scale Invariant Feature Transform，尺度不變特征變換匹配算法）由David G. Lowe教授在1999年（《Object Recognition from Local Scale-Invariant Features》）提出，在2004年（《Distinctive Image Features from Scale-Invariant Keypoints》）得以完善。SIFT可以應用到物體辨識、機器人地圖感知與導航、影像縫合、3D模型建立、手勢辨識、影像追蹤和動作比對等方向。SIFT 特征是基於物體上的一些局部外觀的興趣點而與影像的大小和旋轉無關。對於光線、噪聲、些微視角改變的容忍度也相當高。基於這些特性，它們是高度顯著而且相對容易擷取，在母數龐大的特征數據庫中，很容易辨識物體而且鮮有誤認。使用 SIFT特征描述對於部分物體遮蔽的偵測率也相當高，甚至只需要3個以上的SIFT物體特征就足以計算出位置與方位。在現今的電腦硬件速度下和小型的特征數據庫條件下，辨識速度可接近即時運算。SIFT特征的信息量大，適合在海量數據庫中快速准確匹配，最大的缺點是計算量大、運行時間相對較長。
2 SIFT特征原理
2.1概述
SIFT算法通過金字塔出現的構造進行不同尺度的特征點的提取和描述子的計算。整體可分為四個步驟：尺度空間極值檢測、關鍵點定位、方向計算和關鍵點描述子計算。下面按步驟介紹。
2.2尺度空間極值檢測：
尺度空間方法將傳統的單尺度圖像信息處理技術納入尺度不斷變化的動態分析框架中，更容易獲取圖像的本質特征。可用於表示輪廓在不同尺度下的形狀、大小和位置等。通過構建尺度空間，滿足檢測出的像素點的尺度不變特性。
尺度空間的構建可采用高斯模糊卷積核進行構建，SIFT算法中通過高斯核與圖像的卷積操作，達到圖像不同尺度下的表示效果。詳細如下。
高斯卷積核中每個值的大小如下：
其中σ表示方差，通過改變σ值的大小，獲取不同尺度對應的高斯卷積核。高斯卷積核與圖像的公式如下：
L表示輸出圖像，G表示二維高斯卷積核，I表示輸入圖像的像素值。*表示進行卷積運算。σ 可表示為尺度空間因子。
在SIFT算法中，采用圖像金字塔的方式表示尺度空間因子。整個金字塔可分為兩部分：從上到下分為多組和每組包含多層。
如上圖，最下面一組的第一層為原圖像。第二組的第一層圖像通過對第一組的倒數第三層圖像下采樣獲取。依次類推直到得到所需的組數。對每一組，通過采用不同σ大小的高斯卷積核進行卷積得到多層圖像。其中σ的公式如下：
其中o表示為組數，s表示為層數,S表示為層總數。
2002年，研究者發現尺度歸一化的高斯拉普拉斯函數相對於其它的特征提取函數，可產生更穩定的特征。其中高斯拉普拉斯函數:
而早在1994年，有研究者發現高斯差分函數（簡稱DOG算子）與尺度歸一化的高斯拉普拉斯函數非常近似。其中高斯拉普拉斯函數與D（x,y, σ）的關系可通過如下推導：
SIFT算法作者使用高效的高斯差分替換拉普拉斯算子進行角點檢測。公式如下：
D（x,y,σ) = (G(x,y,kσ)-G(x,y,σ))*I(x,y) = L(x,y,kσ)-L(x,y,σ)
在構建好尺度空間后，為尋找空間中的極值點，SIFT算法通過將像素點與周圍所有相鄰點的像素值的大小進行比較，判斷其是否為極值點。如下圖
將中間點與周圍8個點和上下兩層共9*2共26點進行比較，確保在尺度空間和二維空間均檢測到極值點。（判斷是否在相鄰中，像素值最大或最小）。
2.3 關鍵點定位：
通過 2.2可得到離散空間的極值點，SIFT算法通過擬合二次函數，精確關鍵點的位置，同時去除低對比度的關鍵點和不穩定的邊緣相應點以增強匹配穩定性和提高抗噪能力。
通過上述計算，得到的像素點坐標均以整數為單位，在像素點之間仍有無數的點未考慮。下圖為離散空間極值點與連續空間極值點之間的區別。
為求取連續空間上的極值點，可通過對尺度空間DOG函數進行泰勒展開求解：
其中，X=（x,y,σ)T。對上式求導，並讓方程等於零，得到極值點的偏移量和方程的值分別為：
代表相對插值中心的偏移量，當它的任一維度上的偏移量大於0.5時，以為插值中心偏移到相鄰的像素點上，應改變當前的關鍵點的位置，在新的位置上反復插值直到收斂。同時。|D(x)|小於某個值的點時，易受噪聲的干擾而變得不穩定，SIFT算法將|D(x)|小於某個經驗值的極值點刪除。
通過DOG算子計算得到極值點，可產生較強的邊緣響應，需要剔除不穩定的邊緣響應點。SIFT算法通過采用Harris檢測算法進一步削弱邊緣不穩定性，詳細如下。
首先得到一個Hessian矩陣：
H的特征值α和β可表示為x和y方向的梯度，得下列公式：
Trace(H) = Dxx+Dyy = α+β,
Det(H) = DxxDyy -(Dxy)2 =αβ
對α和β,假設α=rβ，得到：
當r=1時，上式有最小值。由Harris篇可知，當α和β的比例越大時（即像素值一個方向變化很大，另一方向變化小），極可能為邊緣上的點。故為剔除邊緣響應點，需讓該比值小於一定的閾值，可通過檢測在某閾值r下，是否符合下列等式即可(SIFT發表文章中,r取值為10)：
2.4 關鍵點方向計算：
通過上述算法得到的關鍵點仍不具備旋轉不變形，SIFT通過采集其所在高斯金字塔圖像的3σ領域窗口像素的幅度值和方向分布特征獲取角度信息，計算如下：
其中，m(x,y)表示梯度值，θ表示該像素點的角度。通過上述公式計算像素點周圍鄰域窗口每個點的像素的幅度值和方向大小，再通過統計直方圖的方式進行方向的確定，詳細如下。
首先將0~360度划分為36個部分，每個部分的范圍為10°，進行直方圖統計，圖像表示如下：
其中，峰值代表關鍵點的主方向。同時為增強特征匹配的魯棒性，保留峰值大於主方向峰值80%的方向作為關鍵點的輔方向。
2.5 描述子計算：
通過上述步驟，計算出的關鍵點擁有三個信息：位置、尺度和方向。在進行特征匹配時，仍需要一個描述符用於特征匹配。SIFT算法描述符的計算可分為4部分，下面按步驟進行介紹。
2.5.1 確定計算描述子所需的圖像區域：
特征描述子的計算與點所在尺度有關，因此對關鍵點的梯度的計算應在對應的高斯圖像上進行。SIFT算法將關鍵點附近的鄰域划分為d*d（一般為4）個子區域，每個子區域做為一個種子點，每個種子點有8個方向。每個子區域的大小與2.4計算關鍵點的方向相同，均為3σ。
2.5.2 旋轉圖像：
對於2.5.1的圖像區域，根據2.4計算得到的主方向，進行旋轉，旋轉公式如下：
2.5.3 子區域種子點的角度分配：
對圖像進行旋轉后，計算每個子區域內的梯度值和方向，通過加權函數分配到對應種子點的8個方向上。
如上圖所示，其中右右圖可表示為每個種子點的方向信息。
2.5.4 描述子計算：
通過2.5.3，得到的4*4*8=128個梯度信息即為該關鍵點的特征向量。在形成特征向量后，為去除光照變化的影響，需進行歸一化處理。歸一化的數據即為該關鍵點的描述子。
3 代碼實現
3.1 API和demo：
API:
/* 在opencv中，暫時沒有SIFT的接口，需要通過安裝opencv的擴展庫opencv-contrib使用
*/
cv::Ptr<cv::xfeatures2d::SIFT> sift = cv::xfeatures2d::SIFT::create(int nfeatures = 0, int nOctaveLayers = 3, double contrastThreshold = 0.04, double edgeThreshold = 10, double sigma = 1.6);
//特征的數量，0表示計算每個特征方向，１表示不計算方向（速度更快）;表示金子塔的組數；用於過濾區域中的弱特征，閾值越大，特征越少；用於過濾類似邊緣的閾值；高斯輸入層級。
Demo:
#include"opencv2/opencv.hpp" #include"opencv2/xfeatures2d.hpp" using namespace std; using namespace cv; int main() { Mat src1 = imread("T-.png"); Mat src2 = imread("DT.png"); //特征點提取方法 cv::Ptr<cv::xfeatures2d::SIFT> sift = cv::xfeatures2d::SIFT::create(); //特征點提取 vector<KeyPoint> kp1,kp2; sift->detect(src1,kp1); sift->detect(src2,kp2); //畫出特征點 Mat keypointImage1,key_pointImage2,descriptor1,descriptor2; drawKeypoints(src1,kp1,keypointImage1,cv::Scalar::all(-1),DrawMatchesFlags::DEFAULT); imshow("SIFTKP1",keypointImage1); drawKeypoints(src2,kp2,key_pointImage2,cv::Scalar(-1),DrawMatchesFlags::DEFAULT); imshow("SIFTKP2",key_pointImage2); cv::Mat match_pointL,match_pointR; sift->detectAndCompute(src1,cv::Mat(),kp1,match_pointL); sift->detectAndCompute(src2,cv::Mat(),kp2,match_pointR); if(match_pointL.type()!=CV_32F||match_pointR.type()!=CV_32F) { match_pointL.convertTo(match_pointL,CV_32F); match_pointL.convertTo(match_pointR,CV_32F); } vector<DMatch> matches; cv::Ptr<cv::DescriptorMatcher> matcher = cv::DescriptorMatcher::create("FlannBased"); matcher->match(match_pointL,match_pointR,matches); double maxDist = 10; for(int i =0;i<match_pointL.rows;i++) { double dist = matches[i].distance; if(dist>maxDist) maxDist= dist; } vector<DMatch> good_matches,need_matches; vector<float>distance; for(int i =0;i<match_pointL.rows;i++) { if(matches[i].distance<0.1*maxDist) ////調參褚 0.1越小 越精確 官方推薦0.5 如果確定點 可改變 { good_matches.push_back(matches[i]); distance.push_back(matches[i].distance); } } Mat imageOutput; cv::drawMatches(src1,kp1,src2,kp2,good_matches,imageOutput); vector<Point>first_image,second_image; for(int i =0;i<need_matches.size();i++) { first_image.push_back(kp1[need_matches[i].queryIdx].pt); second_image.push_back(kp2[need_matches[i].trainIdx].pt); } cv::imshow("匹配圖片",imageOutput); waitKey(0); return 0; }
源代碼：
/* 假設nOctaveLayers = 3 第0族第0幅： sigma=1.6 第0族第1幅： sigma=1.6 2^(1./3)*sigma 第0族第2幅： sigma=1.6 2^(2./3)*sigma 第0族第3幅： sigma=1.6 2^(3./3)*sigma 第0族第4幅： sigma=1.6 2^(4./3)*sigma 第0族第5幅： sigma=1.6 2^(5./3)*sigma 第0族-->第1族 大小變為一半 尺度為2*sigma 用的是第三幅圖像的sigma 第1族第0幅： sigma=1.6 2^(3./3)*sigma 第1族第1幅： sigma=1.6 2^(4./3)*sigma 第1族第2幅： sigma=1.6 2^(5./3)*sigma 第1族第3幅： sigma=1.6 2^(6./3)*sigma 第1族第4幅： sigma=1.6 2^(7./3)*sigma 第1族第5幅： sigma=1.6 2^(8./3)*sigma ..... 第o族第i幅： sigma=1.6 2^(o + i/3))*sigma gaussian pyramid 一共有nOctaves * （nOctaveLayers + 3）幅圖像 */ void SIFT::buildGaussianPyramid( const Mat& base, vector<Mat>& pyr, int nOctaves ) const { vector<double> sig(nOctaveLayers + 3); pyr.resize(nOctaves*(nOctaveLayers + 3)); // precompute Gaussian sigmas using the following formula: // \sigma_{total}^2 = \sigma_{i}^2 + \sigma_{i-1}^2 // 之所以做 \sigma_{i}^2 = \sigma_{total}^2 - \sigma_{i-1}^2 // 是因為每一張圖像都是在上一張圖像的基礎上做的高斯模糊，累計起來其實就是對原始圖像做了2^(o-1)*2^(n-1/nOctaveLayers)尺度的高斯模糊 sig[0] = sigma; double k = pow( 2., 1. / nOctaveLayers ); for( int i = 1; i < nOctaveLayers + 3; i++ ) { double sig_prev = pow(k, (double)(i-1))*sigma; double sig_total = sig_prev*k; sig[i] = std::sqrt(sig_total*sig_total - sig_prev*sig_prev); } for( int o = 0; o < nOctaves; o++ ) { for( int i = 0; i < nOctaveLayers + 3; i++ ) { Mat& dst = pyr[o*(nOctaveLayers + 3) + i]; if( o == 0 && i == 0 ) dst = base; // base of new octave is halved image from end of previous octave else if( i == 0 ) { const Mat& src = pyr[(o-1)*(nOctaveLayers + 3) + nOctaveLayers]; resize(src, dst, Size(src.cols/2, src.rows/2), 0, 0, INTER_NEAREST); } else { const Mat& src = pyr[o*(nOctaveLayers + 3) + i-1]; GaussianBlur(src, dst, Size(), sig[i], sig[i]); } } } } /* DOG Pyramid 一共有nOctaves * （nOctaveLayers + 2）幅圖像 */ void SIFT::buildDoGPyramid( const vector<Mat>& gpyr, vector<Mat>& dogpyr ) const { int nOctaves = (int)gpyr.size()/(nOctaveLayers + 3); dogpyr.resize( nOctaves*(nOctaveLayers + 2) ); for( int o = 0; o < nOctaves; o++ ) { for( int i = 0; i < nOctaveLayers + 2; i++ ) { const Mat& src1 = gpyr[o*(nOctaveLayers + 3) + i]; const Mat& src2 = gpyr[o*(nOctaveLayers + 3) + i + 1]; Mat& dst = dogpyr[o*(nOctaveLayers + 2) + i]; subtract(src2, src1, dst, noArray(), DataType<sift_wt>::type); } } } // // Interpolates a scale-space extremum's location and scale to subpixel // accuracy to form an image feature. Rejects features with low contrast. // Based on Section 4 of Lowe's paper. static bool adjustLocalExtrema( const vector<Mat>& dog_pyr, KeyPoint& kpt, int octv, int& layer, int& r, int& c, int nOctaveLayers, float contrastThreshold, float edgeThreshold, float sigma ) { const float img_scale = 1.f/(255*SIFT_FIXPT_SCALE); //圖像的尺度，將uchar型的圖像數轉化為0~1之間的圖像 const 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, contr=0; int i = 0; for( ; i < SIFT_MAX_INTERP_STEPS; i++ ) //SIFT_MAX_INTERP_STEPS為迭代次數 { int idx = octv*(nOctaveLayers+2) + layer; //取該特征點所在差分金字塔的位置索引 const Mat& img = dog_pyr[idx]; const Mat& prev = dog_pyr[idx-1]; const Mat& next = dog_pyr[idx+1]; Vec3f dD((img.at<sift_wt>(r, c+1) - img.at<sift_wt>(r, c-1))*deriv_scale, (img.at<sift_wt>(r+1, c) - img.at<sift_wt>(r-1, c))*deriv_scale, (next.at<sift_wt>(r, c) - prev.at<sift_wt>(r, c))*deriv_scale); //分別是對x,y,z軸的差分組成的差分（微分）向量 float v2 = (float)img.at<sift_wt>(r, c)*2; float dxx = (img.at<sift_wt>(r, c+1) + img.at<sift_wt>(r, c-1) - v2)*second_deriv_scale;//對x軸的二階微分 float dyy = (img.at<sift_wt>(r+1, c) + img.at<sift_wt>(r-1, c) - v2)*second_deriv_scale;//對y軸的二階微分 float dss = (next.at<sift_wt>(r, c) + prev.at<sift_wt>(r, c) - v2)*second_deriv_scale; //對z軸的二階微分 float dxy = (img.at<sift_wt>(r+1, c+1) - img.at<sift_wt>(r+1, c-1) - img.at<sift_wt>(r-1, c+1) + img.at<sift_wt>(r-1, c-1))*cross_deriv_scale; //對x，y軸的分別微分 float dxs = (next.at<sift_wt>(r, c+1) - next.at<sift_wt>(r, c-1) - prev.at<sift_wt>(r, c+1) + prev.at<sift_wt>(r, c-1))*cross_deriv_scale; //對x，z軸的分別微分 float dys = (next.at<sift_wt>(r+1, c) - next.at<sift_wt>(r-1, c) - prev.at<sift_wt>(r+1, c) + prev.at<sift_wt>(r-1, c))*cross_deriv_scale; //對y，z軸的分別微分 Matx33f H(dxx, dxy, dxs, dxy, dyy, dys, dxs, dys, dss); //組成的海森矩陣 //求解Hx=dD 負梯度方向（二階梯度下降法） Vec3f X = H.solve(dD, DECOMP_LU); //梯度下降方向 xi = -X[2]; xr = -X[1]; xc = -X[0]; //設置終止條件，證明該點已經為極值點 if( std::abs(xi) < 0.5f && std::abs(xr) < 0.5f && std::abs(xc) < 0.5f ) break; //如果該點為發散的，則返回錯誤 if( std::abs(xi) > (float)(INT_MAX/3) || std::abs(xr) > (float)(INT_MAX/3) || std::abs(xc) > (float)(INT_MAX/3) ) return false; //進行梯度下降 c += cvRound(xc); r += cvRound(xr); layer += cvRound(xi); //如果點坐標（x，y，z）越界 則同樣返回錯誤 if( layer < 1 || layer > nOctaveLayers || c < SIFT_IMG_BORDER || c >= img.cols - SIFT_IMG_BORDER || r < SIFT_IMG_BORDER || r >= img.rows - SIFT_IMG_BORDER ) return false; } // ensure convergence of interpolation 驗證迭代次數沒有大於最大迭代次數 if( i >= SIFT_MAX_INTERP_STEPS ) return false; { int idx = octv*(nOctaveLayers+2) + layer; //取該特征點所在差分金字塔的位置索引 const Mat& img = dog_pyr[idx]; const Mat& prev = dog_pyr[idx-1]; const Mat& next = dog_pyr[idx+1]; Matx31f dD((img.at<sift_wt>(r, c+1) - img.at<sift_wt>(r, c-1))*deriv_scale, (img.at<sift_wt>(r+1, c) - img.at<sift_wt>(r-1, c))*deriv_scale, (next.at<sift_wt>(r, c) - prev.at<sift_wt>(r, c))*deriv_scale); //微分向量 float t = dD.dot(Matx31f(xc, xr, xi)); //內積 將微分向量與梯度下降的方向做內積 //？？？？？？？？？？？？？？？？？？？？？？？？ contr = img.at<sift_wt>(r, c)*img_scale + t * 0.5f; if( std::abs( contr ) * nOctaveLayers < contrastThreshold ) return false; // principal curvatures are computed using the trace and det of Hessian // 避免計算特征值的麻煩，海森矩陣的特征值與主曲率成正比 ，這里用Tr(H)^2/Det(H)>edgeThreshold來近似代替主曲率 // 設a,b分別為海森矩陣特征值，則令r=a/b Tr(H)=a+b Det(H)=a*b Tr(H)^2/Det(H)=(r+1)^2/r float v2 = img.at<sift_wt>(r, c)*2.f; float dxx = (img.at<sift_wt>(r, c+1) + img.at<sift_wt>(r, c-1) - v2)*second_deriv_scale; float dyy = (img.at<sift_wt>(r+1, c) + img.at<sift_wt>(r-1, c) - v2)*second_deriv_scale; float dxy = (img.at<sift_wt>(r+1, c+1) - img.at<sift_wt>(r+1, c-1) - img.at<sift_wt>(r-1, c+1) + img.at<sift_wt>(r-1, c-1)) * cross_deriv_scale; float tr = dxx + dyy; float det = dxx * dyy - dxy * dxy; if( det <= 0 || tr*tr*edgeThreshold >= (edgeThreshold + 1)*(edgeThreshold + 1)*det ) return false; } // 特征點的精確位置（對應原始圖像的位置（x,y）） kpt.pt.x = (c + xc) * (1 << octv); kpt.pt.y = (r + xr) * (1 << octv); // 特征點所在的族 kpt.octave = octv + (layer << 8) + (cvRound((xi + 0.5)*255) << 16); // 該特征點所在的尺度大小為：該特征點所在尺度 2^(octv+layer/nOctaveLayers)*sigma * 2 因為平時所說的尺度為半徑，這里是長度 kpt.size = sigma*powf(2.f, (layer + xi) / nOctaveLayers)*(1 << octv)*2; kpt.response = std::abs(contr); return true; } // // Detects features at extrema in DoG scale space. Bad features are discarded // based on contrast and ratio of principal curvatures. void SIFT::findScaleSpaceExtrema( const vector<Mat>& gauss_pyr, const vector<Mat>& dog_pyr, vector<KeyPoint>& keypoints ) const { //由於gauss_pyr一共存在nOctaves*(nOctaveLayers + 3)幅圖像 所以族數nOctaves = gauss_pyr.size()/(nOctaveLayers + 3); int nOctaves = (int)gauss_pyr.size()/(nOctaveLayers + 3); int threshold = cvFloor(0.5 * contrastThreshold / nOctaveLayers * 255 * SIFT_FIXPT_SCALE); const int n = SIFT_ORI_HIST_BINS; //將方向分為36個bin float hist[n]; KeyPoint kpt; keypoints.clear(); for( int o = 0; o < nOctaves; o++ ) for( int i = 1; i <= nOctaveLayers; i++ ) //注意這里是從第1張圖像到第nOctaveLayers圖像 而第0張圖像和第nOctaveLayers+1張圖像將不作處理，也是為什么我們要取nOctaveLayers+3張高斯圖像的原因 { int idx = o*(nOctaveLayers+2)+i; //取差分金字塔當前族當前層的圖像及其上一張圖像和下一張圖像 const Mat& img = dog_pyr[idx]; const Mat& prev = dog_pyr[idx-1]; const Mat& next = dog_pyr[idx+1]; int step = (int)img.step1(); int rows = img.rows, cols = img.cols; for( int r = SIFT_IMG_BORDER; r < rows-SIFT_IMG_BORDER; r++) { const sift_wt* currptr = img.ptr<sift_wt>(r); const sift_wt* prevptr = prev.ptr<sift_wt>(r); const sift_wt* nextptr = next.ptr<sift_wt>(r); for( int c = SIFT_IMG_BORDER; c < cols-SIFT_IMG_BORDER; c++) { sift_wt val = currptr[c]; // find local extrema with pixel accuracy 如果該點大於三維坐標系下周邊所有臨近的26個點，才證明該點有可能是極值點 if( std::abs(val) > threshold && ((val > 0 && val >= currptr[c-1] && val >= currptr[c+1] && val >= currptr[c-step-1] && val >= currptr[c-step] && val >= currptr[c-step+1] && val >= currptr[c+step-1] && val >= currptr[c+step] && val >= currptr[c+step+1] && val >= nextptr[c] && val >= nextptr[c-1] && val >= nextptr[c+1] && val >= nextptr[c-step-1] && val >= nextptr[c-step] && val >= nextptr[c-step+1] && val >= nextptr[c+step-1] && val >= nextptr[c+step] && val >= nextptr[c+step+1] && val >= prevptr[c] && val >= prevptr[c-1] && val >= prevptr[c+1] && val >= prevptr[c-step-1] && val >= prevptr[c-step] && val >= prevptr[c-step+1] && val >= prevptr[c+step-1] && val >= prevptr[c+step] && val >= prevptr[c+step+1]) || (val < 0 && val <= currptr[c-1] && val <= currptr[c+1] && val <= currptr[c-step-1] && val <= currptr[c-step] && val <= currptr[c-step+1] && val <= currptr[c+step-1] && val <= currptr[c+step] && val <= currptr[c+step+1] && val <= nextptr[c] && val <= nextptr[c-1] && val <= nextptr[c+1] && val <= nextptr[c-step-1] && val <= nextptr[c-step] && val <= nextptr[c-step+1] && val <= nextptr[c+step-1] && val <= nextptr[c+step] && val <= nextptr[c+step+1] && val <= prevptr[c] && val <= prevptr[c-1] && val <= prevptr[c+1] && val <= prevptr[c-step-1] && val <= prevptr[c-step] && val <= prevptr[c-step+1] && val <= prevptr[c+step-1] && val <= prevptr[c+step] && val <= prevptr[c+step+1]))) { int r1 = r, c1 = c, layer = i; //將極值點進行差值精確找到極值點的位置，提高精度的位極值點的x,y坐標以及該極值點所處的層數layer 注意並不優化極值點所在的族數 //並且該函數將剔除差分金字塔的邊緣效應 if( !adjustLocalExtrema(dog_pyr, kpt, o, layer, r1, c1, nOctaveLayers, (float)contrastThreshold, (float)edgeThreshold, (float)sigma) ) continue; //計算該特征點所在層數相對於該族第0層圖像的尺度比例 sigma float scl_octv = kpt.size*0.5f/(1 << o); //計算該特征點在高斯金字塔下的方向梯度直方圖 float omax = calcOrientationHist(gauss_pyr[o*(nOctaveLayers+3) + layer], Point(c1, r1), cvRound(SIFT_ORI_RADIUS * scl_octv), //// 9sigma/2 SIFT_ORI_SIG_FCTR * scl_octv, // 3sigma/2 hist, n); float mag_thr = (float)(omax * SIFT_ORI_PEAK_RATIO); //將所有hist值大於主方向對應hist值得所有方向都作為特征點的方向加入特征點vector容器中 找尋輔助方向 for( int j = 0; j < n; j++ ) { int l = j > 0 ? j - 1 : n - 1; int r2 = j < n-1 ? j + 1 : 0; if( hist[j] > hist[l] && hist[j] > hist[r2] && hist[j] >= mag_thr ) { float bin = j + 0.5f * (hist[l]-hist[r2]) / (hist[l] - 2*hist[j] + hist[r2]); bin = bin < 0 ? n + bin : bin >= n ? bin - n : bin; kpt.angle = 360.f - (float)((360.f/n) * bin); if(std::abs(kpt.angle - 360.f) < FLT_EPSILON) kpt.angle = 0.f; keypoints.push_back(kpt); } } } } } } } // Computes a gradient orientation histogram at a specified pixel //在高斯金字塔下計算特征點的梯度直方圖 static float calcOrientationHist( const Mat& img, Point pt, int radius, float sigma, float* hist, int n ) { int i, j, k, len = (radius*2+1)*(radius*2+1); float expf_scale = -1.f/(2.f * sigma * sigma); AutoBuffer<float> buf(len*4 + n+4); //為下面的變量開辟地址空間 float *X = buf, *Y = X + len, *Mag = X, *Ori = Y + len, *W = Ori + len; float* temphist = W + len + 2; for( i = 0; i < n; i++ ) //36個bin temphist數組存儲該bin的權值 temphist[i] = 0.f; for( i = -radius, k = 0; i <= radius; i++ ) { int y = pt.y + i; if( y <= 0 || y >= img.rows - 1 ) continue; for( j = -radius; j <= radius; j++ ) { int x = pt.x + j; if( x <= 0 || x >= img.cols - 1 ) continue; float dx = (float)(img.at<sift_wt>(y, x+1) - img.at<sift_wt>(y, x-1)); float dy = (float)(img.at<sift_wt>(y-1, x) - img.at<sift_wt>(y+1, x)); // X為每個特征點對應范圍內x方向的梯度 Y對應的是每個特征點對應范圍內y方向的梯度 W是每個特征點對應范圍內的權值 X[k] = dx; Y[k] = dy; W[k] = (i*i + j*j)*expf_scale; k++; } } len = k; // compute gradient values, orientations and the weights over the pixel neighborhood // 梯度的權值 根據高斯核的尺度sigma決定（半徑大小） exp(W, W, len); // 梯度的幅角 fastAtan2(Y, X, Ori, len, true); // 梯度的幅值 magnitude(X, Y, Mag, len); for( k = 0; k < len; k++ ) { int bin = cvRound((n/360.f)*Ori[k]); //計算當前角度屬於哪個bin if( bin >= n ) bin -= n; if( bin < 0 ) bin += n; temphist[bin] += W[k]*Mag[k]; //根據幅值和權重加入到該bin的方向直方圖中 } // smooth the histogram 對直方圖進行濾波（高斯濾波） temphist[-1] = temphist[n-1]; temphist[-2] = temphist[n-2]; temphist[n] = temphist[0]; temphist[n+1] = temphist[1]; for( i = 0; i < n; i++ ) { hist[i] = (temphist[i-2] + temphist[i+2])*(1.f/16.f) + (temphist[i-1] + temphist[i+1])*(4.f/16.f) + temphist[i]*(6.f/16.f); } //得到主方向 對應最大的值 float maxval = hist[0]; for( i = 1; i < n; i++ ) maxval = std::max(maxval, hist[i]); return maxval; } //計算SIFT描述子 static void calcSIFTDescriptor( const Mat& img, Point2f ptf, float ori, float scl, int d, int n, float* dst ) { Point pt(cvRound(ptf.x), cvRound(ptf.y)); float cos_t = cosf(ori*(float)(CV_PI/180)); float sin_t = sinf(ori*(float)(CV_PI/180)); float bins_per_rad = n / 360.f; float exp_scale = -1.f/(d * d * 0.5f); float hist_width = SIFT_DESCR_SCL_FCTR * scl; //3*sigma 3sigma之外的點對當前點的影響不大，關聯性不大（3sigma原則） int radius = cvRound(hist_width * 1.4142135623730951f * (d + 1) * 0.5f); //滿足旋轉型和三次線性插值 將半徑變為此 // Clip the radius to the diagonal of the image to avoid autobuffer too large exception radius = std::min(radius, (int) sqrt((double) img.cols*img.cols + img.rows*img.rows)); cos_t /= hist_width; //為了計算時y''= 1/3sigma*(j * cos_t - i * sin_t) sin_t /= hist_width; //為了計算時x''= 1/3sigma*(j * sin_t + i * cos_t) int i, j, k, len = (radius*2+1)*(radius*2+1), histlen = (d+2)*(d+2)*(n+2); int rows = img.rows, cols = img.cols; AutoBuffer<float> buf(len*6 + histlen); float *X = buf, *Y = X + len, *Mag = Y, *Ori = Mag + len, *W = Ori + len; float *RBin = W + len, *CBin = RBin + len, *hist = CBin + len; for( i = 0; i < d+2; i++ ) { for( j = 0; j < d+2; j++ ) for( k = 0; k < n+2; k++ ) hist[(i*(d+2) + j)*(n+2) + k] = 0.; } for( i = -radius, k = 0; i <= radius; i++ ) for( j = -radius; j <= radius; j++ ) { // Calculate sample's histogram array coords rotated relative to ori. // Subtract 0.5 so samples that fall e.g. in the center of row 1 (i.e. // r_rot = 1.5) have full weight placed in row 1 after interpolation. //c_rot r_rot范圍在-d/2~d/2 之間 cbin rbin范圍在0~d之間 r,c是對應點的坐標 float c_rot = j * cos_t - i * sin_t; float r_rot = j * sin_t + i * cos_t; float rbin = r_rot + d/2 - 0.5f; float cbin = c_rot + d/2 - 0.5f; int r = pt.y + i, c = pt.x + j; if( rbin > -1 && rbin < d && cbin > -1 && cbin < d && r > 0 && r < rows - 1 && c > 0 && c < cols - 1 ) { float dx = (float)(img.at<sift_wt>(r, c+1) - img.at<sift_wt>(r, c-1)); float dy = (float)(img.at<sift_wt>(r-1, c) - img.at<sift_wt>(r+1, c)); //X[k] Y[k]存儲的是對應點相對於x軸y軸的梯度 RBin[k] CBin[k]存儲的是對應點轉化成0-d的坐標 X[k] = dx; Y[k] = dy; RBin[k] = rbin; CBin[k] = cbin; //權重 W[k] = (c_rot * c_rot + r_rot * r_rot)*exp_scale; k++; } } // radius*radius的范圍內點的個數 len = k; //計算梯度的方向，幅值，以及權重 fastAtan2(Y, X, Ori, len, true); magnitude(X, Y, Mag, len); exp(W, W, len); //遍歷radius*radius的范圍內所有的點 for( k = 0; k < len; k++ ) { float rbin = RBin[k], cbin = CBin[k]; float obin = (Ori[k] - ori)*bins_per_rad; //相對於主方向的角度 float mag = Mag[k]*W[k]; //當前點的梯度權值 int r0 = cvFloor( rbin ); int c0 = cvFloor( cbin ); int o0 = cvFloor( obin ); rbin -= r0; //dr 對於鄰近兩行的貢獻因子dr和1-dr cbin -= c0; //dc 對於鄰近兩列的貢獻因子dc和1-dc obin -= o0; //do 對於鄰近兩個方向的貢獻因子為do和 1-do if( o0 < 0 ) o0 += n; if( o0 >= n ) o0 -= n; // histogram update using tri-linear interpolation 三次線性插值 計算每個點對於所處的bin以及相鄰bin的權重 float v_r1 = mag*rbin, v_r0 = mag - v_r1; float v_rc11 = v_r1*cbin, v_rc10 = v_r1 - v_rc11; float v_rc01 = v_r0*cbin, v_rc00 = v_r0 - v_rc01; float v_rco111 = v_rc11*obin, v_rco110 = v_rc11 - v_rco111; float v_rco101 = v_rc10*obin, v_rco100 = v_rc10 - v_rco101; float v_rco011 = v_rc01*obin, v_rco010 = v_rc01 - v_rco011; float v_rco001 = v_rc00*obin, 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; } // finalize histogram, since the orientation histograms are circular 最終的梯度直方圖 for( i = 0; i < d; i++ ) for( 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]; for( k = 0; k < n; k++ ) dst[(i*d + j)*n + k] = hist[idx+k]; } // copy histogram to the descriptor, // apply hysteresis thresholding // and scale the result, so that it can be easily converted // to byte array float nrm2 = 0; // 4*4*8 描述子向量的個數 //描述子中128維向量存儲的是4*4個小方格里8個方向分別的權重和 （先歸一化然后在將0-1float數轉化為 255的uchar數） len = d*d*n; for( k = 0; k < len; k++ ) nrm2 += dst[k]*dst[k]; float thr = std::sqrt(nrm2)*SIFT_DESCR_MAG_THR;//添加閾值，向量二范數大於向量歸一化后0.2，則認為是由於非線性光照造成的 for( i = 0, nrm2 = 0; i < k; i++ ) //重新計算h^2的和 { float val = std::min(dst[i], thr); dst[i] = val; nrm2 += val*val; } nrm2 = SIFT_INT_DESCR_FCTR/std::max(std::sqrt(nrm2), FLT_EPSILON); #if 1 for( k = 0; k < len; k++ ) { dst[k] = saturate_cast<uchar>(dst[k]*nrm2); //對向量進行歸一化並轉化為uchar型 } #else float nrm1 = 0; for( k = 0; k < len; k++ ) { dst[k] *= nrm2; nrm1 += dst[k]; } nrm1 = 1.f/std::max(nrm1, FLT_EPSILON); for( k = 0; k < len; k++ ) { dst[k] = std::sqrt(dst[k] * nrm1);//saturate_cast<uchar>(std::sqrt(dst[k] * nrm1)*SIFT_INT_DESCR_FCTR); } #endif } static void calcDescriptors(const vector<Mat>& gpyr, const vector<KeyPoint>& keypoints, Mat& descriptors, int nOctaveLayers, int firstOctave ) { int d = SIFT_DESCR_WIDTH, n = SIFT_DESCR_HIST_BINS; //找尋檢測到的沒各特征點的 for( size_t i = 0; i < keypoints.size(); i++ ) { KeyPoint kpt = keypoints[i]; int octave, layer; float scale; unpackOctave(kpt, octave, layer, scale); CV_Assert(octave >= firstOctave && layer <= nOctaveLayers+2); float size=kpt.size*scale; Point2f ptf(kpt.pt.x*scale, kpt.pt.y*scale); //在高斯金字塔下計算描述子 const Mat& img = gpyr[(octave - firstOctave)*(nOctaveLayers + 3) + layer]; float angle = 360.f - kpt.angle; if(std::abs(angle - 360.f) < FLT_EPSILON) angle = 0.f; //img為高斯金字塔下的圖像 ptf也是高斯金字塔下對應層數的坐標 calcSIFTDescriptor(img, ptf, angle, size*0.5f, d, n, descriptors.ptr<float>((int)i)); } } void SIFT::operator()(InputArray _image, InputArray _mask, vector<KeyPoint>& keypoints, OutputArray _descriptors, bool useProvidedKeypoints) const { int firstOctave = -1, actualNOctaves = 0, actualNLayers = 0; Mat image = _image.getMat(), mask = _mask.getMat(); if( image.empty() || image.depth() != CV_8U ) CV_Error( CV_StsBadArg, "image is empty or has incorrect depth (!=CV_8U)" ); if( !mask.empty() && mask.type() != CV_8UC1 ) CV_Error( CV_StsBadArg, "mask has incorrect type (!=CV_8UC1)" ); if( useProvidedKeypoints ) { firstOctave = 0; int maxOctave = INT_MIN; for( size_t i = 0; i < keypoints.size(); i++ ) { int octave, layer; float scale; unpackOctave(keypoints[i], octave, layer, scale); firstOctave = std::min(firstOctave, octave); maxOctave = std::max(maxOctave, octave); actualNLayers = std::max(actualNLayers, layer-2); } firstOctave = std::min(firstOctave, 0); CV_Assert( firstOctave >= -1 && actualNLayers <= nOctaveLayers ); actualNOctaves = maxOctave - firstOctave + 1; } Mat base = createInitialImage(image, firstOctave < 0, (float)sigma); vector<Mat> gpyr, dogpyr; int nOctaves = actualNOctaves > 0 ? actualNOctaves : cvRound(log( (double)std::min( base.cols, base.rows ) ) / log(2.) - 2) - firstOctave; //double t, tf = getTickFrequency(); //t = (double)getTickCount(); buildGaussianPyramid(base, gpyr, nOctaves); buildDoGPyramid(gpyr, dogpyr); //t = (double)getTickCount() - t; //printf("pyramid construction time: %g\n", t*1000./tf); if( !useProvidedKeypoints ) { //t = (double)getTickCount(); findScaleSpaceExtrema(gpyr, dogpyr, keypoints); //去除重復點 KeyPointsFilter::removeDuplicated( keypoints ); //從中選取最優的前nfeatures個特征點 if( nfeatures > 0 ) KeyPointsFilter::retainBest(keypoints, nfeatures); //t = (double)getTickCount() - t; //printf("keypoint detection time: %g\n", t*1000./tf); if( firstOctave < 0 ) for( size_t i = 0; i < keypoints.size(); i++ ) { KeyPoint& kpt = keypoints[i]; float scale = 1.f/(float)(1 << -firstOctave); kpt.octave = (kpt.octave & ~255) | ((kpt.octave + firstOctave) & 255); kpt.pt *= scale; kpt.size *= scale; } if( !mask.empty() ) KeyPointsFilter::runByPixelsMask( keypoints, mask ); } else { // filter keypoints by mask //KeyPointsFilter::runByPixelsMask( keypoints, mask ); } if( _descriptors.needed() ) { //t = (double)getTickCount(); int dsize = descriptorSize(); _descriptors.create((int)keypoints.size(), dsize, CV_32F); Mat descriptors = _descriptors.getMat(); calcDescriptors(gpyr, keypoints, descriptors, nOctaveLayers, firstOctave); //t = (double)getTickCount() - t; //printf("descriptor extraction time: %g\n", t*1000./tf); } }
參考：
1.基於SIFT特征目標跟蹤算法研究 藺海峰 馬宇峰 宋濤
2.https://blog.csdn.net/u010440456/article/details/81483145
3.https://blog.csdn.net/qq_30356613/article/details/78534101