1、創建高斯金字塔第-1組
1.1、將源圖片轉成灰度圖
void ConvertToGray(const Mat& src, Mat& dst) { cv::Size size = src.size(); if (dst.empty()) dst.create(size, CV_64F); //[1]利用Mat類的成員函數創建Mat容器 uchar* srcData = src.data; //[2]指向存儲所有像素值的矩陣的數據區域 pixel_t* dstData = (pixel_t*)dst.data; //[3]指向dst的矩陣體數據區域 int dstStep = dst.step / sizeof(dstData[0]); //[4]一行圖像含有多少字節數 for (int j = 0; j<src.cols; j++) { for (int i = 0; i<src.rows; i++) { double b = (srcData + src.step*i + src.channels()*j)[0] / 255.0; double g = (srcData + src.step*i + src.channels()*j)[1] / 255.0; double r = (srcData + src.step*i + src.channels()*j)[2] / 255.0; *(dstData + dstStep*i + dst.channels()*j) = (r + g + b) / 3.0; } } }
1.2、進行上采樣
/************************************************************************************************************************* *模塊說明： * 線性插值放大，將圖像放大2倍 **************************************************************************************************************************/ void UpSample(const Mat &src, Mat &dst) { if (src.channels() != 1) return; dst.create(src.rows * 2, src.cols * 2, src.type()); //dst的行列是原圖2倍大小的行列 pixel_t* srcData = (pixel_t*)src.data; pixel_t* dstData = (pixel_t*)dst.data; int srcStep = src.step / sizeof(srcData[0]); //一行圖像含有多少字節數 int dstStep = dst.step / sizeof(dstData[0]); int m = 0, n = 0; for (int j = 0; j < src.cols - 1; j++, n += 2) { m = 0; for (int i = 0; i < src.rows - 1; i++, m += 2) { /*第1步：n=0 m=0，n=2,m=0; 第2步：n=0 m=1，n=2,m=1; 第3步：n=1 m=0，n=3,m=0; 第4步：n=1 m=1，n=3,m=1; 第5步：最后2行2列 */ /*第1步：n=0 m=2，n=2,m=2; 第2步：n=0 m=3，n=2,m=3; 第3步：n=1 m=2，n=3,m=2; 第4步：n=1 m=3，n=3,m=3; 假設5*5的矩陣 */ double sample = *(srcData + srcStep * i + src.channels() * j); *(dstData + dstStep * m + dst.channels() * n) = sample; double rs = *(srcData + srcStep * (i)+src.channels()*j) + (*(srcData + srcStep * (i + 1) + src.channels()*j)); *(dstData + dstStep * (m + 1) + dst.channels() * n) = rs / 2; double cs = *(srcData + srcStep * i + src.channels()*(j)) + (*(srcData + srcStep * i + src.channels()*(j + 1))); *(dstData + dstStep * m + dst.channels() * (n + 1)) = cs / 2; double center = (*(srcData + srcStep * (i + 1) + src.channels() * j)) + (*(srcData + srcStep * i + src.channels() * j)) + (*(srcData + srcStep * (i + 1) + src.channels() * (j + 1))) + (*(srcData + srcStep * i + src.channels() * (j + 1))); *(dstData + dstStep * (m + 1) + dst.channels() * (n + 1)) = center / 4; } } if (dst.rows < 3 || dst.cols < 3) return; //最后兩行兩列 for (int k = dst.rows - 1; k >= 0; k--) { *(dstData + dstStep *(k)+dst.channels()*(dst.cols - 2)) = *(dstData + dstStep *(k)+dst.channels()*(dst.cols - 3)); *(dstData + dstStep *(k)+dst.channels()*(dst.cols - 1)) = *(dstData + dstStep *(k)+dst.channels()*(dst.cols - 3)); } for (int k = dst.cols - 1; k >= 0; k--) { *(dstData + dstStep *(dst.rows - 2) + dst.channels()*(k)) = *(dstData + dstStep *(dst.rows - 3) + dst.channels()*(k)); *(dstData + dstStep *(dst.rows - 1) + dst.channels()*(k)) = *(dstData + dstStep *(dst.rows - 3) + dst.channels()*(k)); }
}
1.3、求該層σ（sigma）
1.4、進行高斯模糊
計算金字塔層數（根據圖像長、寬中最小的那個）
這里呢？low論文中高斯金字塔負一層創建,具體的可參考：https://www.cnblogs.com/xujianqing/p/5988627.html
/*************************************************************************************************************************
*模塊說明：
* OpenCv中的高斯平滑函數
**************************************************************************************************************************/
void GaussianSmooth(const Mat &src, Mat &dst, double sigma)
{
GaussianBlur(src, dst, Size(0, 0), sigma, sigma);
}
/*************************************************************************************************************************
*模塊說明：
* 模塊1--創建初始灰度圖像------初始圖像先將原圖像灰度化，再擴大一倍后，使用高斯模糊平滑
*功能說明：
* 在最開始建立高斯金字塔的時候，要預先模糊輸入的圖像來作為第0個組的第0層圖像，這時相當於丟棄了最高的空域的采樣率。因
* 此，通常的做法是先將圖像的尺度擴大一倍來生成第-1組。我們假定初始的輸入圖像為了抗擊混淆現象，已經對其進行了sigma=0.5
* 的高斯模糊，如果輸入圖像的尺寸用雙線性插值擴大一倍，那么相當於sigma=1.0
*這樣做的原因:
* 在檢測極值點前，對原始圖像的高斯平滑以致圖像丟失高頻信息，所以Lowe建議在建立尺度空間前，首先對原始圖像長寬擴展一倍，
* 以便保留原始圖像的信息(特別是圖像的高頻信息，比如邊緣，角點)，增加特征點的數量。
*函數功能:
* 這個函數的作用是用於創建高斯金字塔的-1層圖像，在主函數中定義 sigma =1.6，INIT_SIGMA=O.5;
在lowe的論文中，他定義圖片的尺度是1.6，被攝像頭模糊的尺度是0.5，但是這個1.52對我們來說並沒有啥意思，
lowe為了獲得更多的特征點（姑且這么解釋吧）他先對原始圖像進行擴大一倍。然后呢再對它進行高斯平滑作為第一階第一層，
高斯模糊尺度為=(1.6* 1.6 - (0.5 * 2) * (0.5 * 2))
**************************************************************************************************************************/
void CreateInitSmoothGray(const Mat &src, Mat &dst, double sigma = SIGMA)
{
cv::Mat gray; //[1]保存原始圖像的灰度圖像
cv::Mat up; //[2]保存原始圖像的上采樣圖像
//std::cout << "CreateInitSmoothGray sigma = " << sigma << std::endl; //查看sigma初始值
ConvertToGray(src, gray); //灰度圖像
UpSample(gray, up); //[3]圖像的上采樣
//[4]高斯金字塔的-1組的sigma的計算公式
double sigma_init = sqrt(sigma * sigma - (INIT_SIGMA * 2) * (INIT_SIGMA * 2));// INIT_SIGM初始sigma=0.5; -1層的sigma
GaussianSmooth(up, dst, sigma_init); //[5]高斯平滑
}
2、創建高斯金字塔
2.1、高斯模糊卷積
2.2、降采樣
隔點采用：
/************************************************************************************************************************* *模塊說明： * 隔點采樣 **************************************************************************************************************************/ void DownSample(const Mat& src, Mat& dst) { if (src.channels() != 1) return; if (src.cols <= 1 || src.rows <= 1) { src.copyTo(dst); return; } dst.create((int)(src.rows / 2), (int)(src.cols / 2), src.type()); //dst的行列=src/2 pixel_t* srcData = (pixel_t*)src.data; pixel_t* dstData = (pixel_t*)dst.data; int srcStep = src.step / sizeof(srcData[0]); //src一行圖像含有多少字節數 int dstStep = dst.step / sizeof(dstData[0]); int m = 0, n = 0; for (int j = 0; j < src.cols; j += 2, n++) { m = 0; for (int i = 0; i < src.rows; i += 2, m++) { pixel_t sample = *(srcData + srcStep * i + src.channels() * j); if (m < dst.rows && n < dst.cols) { *(dstData + dstStep * m + dst.channels() * n) = sample; //將src隔點采樣到dst中去。剛好src行列/2 } } } } /************************************
上述2步原理參考我的博文：SIFT解析（一）建立高斯金字塔
;
也就是說o代表哪一組，S表示哪一層，s表示那一層中第幾幅圖片，也就是第幾層。
S=3,第0（o=0）組的（S+3）6幅圖像（及第0組內有6層圖像）分別表示為（s第幾層的圖片）
又
大家如果覺得還是不清楚，可以讀一下百度文庫中這篇文章中的高斯金字塔構建原理，講解的很清楚：
https://wenku.baidu.com/view/d7edd2464b73f242336c5ffa.html
這部分代碼可以提現為
/************************************************************************************************************************* *模塊說明： * OpenCv中的高斯平滑函數 **************************************************************************************************************************/ void GaussianSmooth(const Mat &src, Mat &dst, double sigma) { GaussianBlur(src, dst, Size(0, 0), sigma, sigma); } /************************************************************************************************************************* *模塊說明： * 模塊1--創建初始灰度圖像------初始圖像先將原圖像灰度化，再擴大一倍后，使用高斯模糊平滑 *功能說明： * 在最開始建立高斯金字塔的時候，要預先模糊輸入的圖像來作為第0個組的第0層圖像，這時相當於丟棄了最高的空域的采樣率。因 * 此，通常的做法是先將圖像的尺度擴大一倍來生成第-1組。我們假定初始的輸入圖像為了抗擊混淆現象，已經對其進行了sigma=0.5 * 的高斯模糊，如果輸入圖像的尺寸用雙線性插值擴大一倍，那么相當於sigma=1.0 *這樣做的原因: * 在檢測極值點前，對原始圖像的高斯平滑以致圖像丟失高頻信息，所以Lowe建議在建立尺度空間前，首先對原始圖像長寬擴展一倍， * 以便保留原始圖像的信息(特別是圖像的高頻信息，比如邊緣，角點)，增加特征點的數量。 *函數功能: * 這個函數的作用是用於創建高斯金字塔的-1層圖像，在主函數中定義 sigma =1.6 在lowe的論文中，他定義圖片的尺度是1.6，被攝像頭模糊的尺度是0.5，但是這個1.52對我們來說並沒有啥意思， lowe為了獲得更多的特征點（姑且這么解釋吧）他先對原始圖像進行擴大一倍。然后呢再對它進行高斯平滑作為第一階第一層， 高斯模糊尺度為=(1.6* 1.6 - (0.5 * 2) * (0.5 * 2)) **************************************************************************************************************************/ void CreateInitSmoothGray(const Mat &src, Mat &dst, double sigma = SIGMA) { cv::Mat gray; //[1]保存原始圖像的灰度圖像 cv::Mat up; //[2]保存原始圖像的上采樣圖像 //std::cout << "CreateInitSmoothGray sigma = " << sigma << std::endl; //查看sigma初始值 ConvertToGray(src, gray); //灰度圖像 UpSample(gray, up); //[3]圖像的上采樣 //[4]高斯金字塔的-1組的sigma的計算公式 double sigma_init = sqrt(sigma * sigma - (INIT_SIGMA * 2) * (INIT_SIGMA * 2));// INIT_SIGM初始sigma=0.5; -1層的sigma GaussianSmooth(up, dst, sigma_init); //[5]高斯平滑 } /************************************************************************************************************************* *模塊說明： * 模塊二：3.3 圖像高斯金字塔的構建 在主函數中：int octaves = log((double)min(init_gray.rows, init_gray.cols)) / log(2.0) - 2; //計算高斯金字塔的層數 高斯金字塔中每組圖像中有三層/張圖片 INTERVALS=3；SIGMA=1.6； **************************************************************************************************************************/ void GaussianPyramid(const Mat &src, vector<Mat>&gauss_pyr, int octaves, int intervals = INTERVALS, double sigma = SIGMA) { //std::cout << "sigma is = " << sigma << std::endl; double *sigmas = new double[intervals + 3]; double k = pow(2.0, 1.0 / intervals); //pow（x,y）,x的y次方；k=2的（1/s)次方 // std::cout << "k is = " << k << std::endl; sigmas[0] = sigma; // std::cout << "sigma0 is = " << sigmas[0] << std::endl; double sig_prev; double sig_total; //用循環來表示第0組內各層圖片的尺度,這個可以表示為0組的各層圖像尺度； for (int i = 1; i < intervals + 3; i++) { sig_prev = pow(k, i - 1) * sigma; //第i-1層尺度 sig_total = sig_prev * k; //第i層尺度 sigmas[i] = sqrt(sig_total * sig_total - sig_prev * sig_prev); //組內每層尺度坐標 } for (int o = 0; o < octaves; o++) { //每組多三層 for (int i = 0; i < intervals + 3; i++) { Mat mat; if (o == 0 && i == 0) { src.copyTo(mat); //第0組，第0層，就是原圖像 } else if (i == 0) //新的一組的首層圖像是由上一組最后一層圖像下采樣得到 { DownSample(gauss_pyr[(o - 1)*(intervals + 3) + intervals], mat); } else //對上一層圖像gauss_pyr[i-1]進行參數為sig[i]的高斯平滑，得到當前層圖像 { GaussianSmooth(gauss_pyr[o * (intervals + 3) + i - 1], mat, sigmas[i]); } gauss_pyr.push_back(mat); //將mat中的像素添加到gauss_pyr數組里面；也就是說每組從第0層到第5層的圖片像素裝進vector數組中； } } delete[] sigmas; //釋放sigmas參數數組 }
但是在原理中，我發現一個很費解的問題
而在論文代碼中，並沒有看出在哪體現出來每一處的2^o次方這個概念。后來查找文章才發現？
在極值比較的過程中，每一組圖像的首末兩層是無法進行極值比較的，為了滿足尺度變化的連續性；
==============================================================================================================================
這里有的童鞋不理解什么叫“為了滿足尺度變化的連續性”，現在做仔細闡述：
假設s=3，也就是每個塔里有3層，則k=2^1/s=2^1/3，那么按照上圖可得Gauss Space和DoG space 分別有3個（s個）和2個（s-1個）分量，
在DoG space中，1st-octave兩項分別是σ,kσ; 2nd-octave兩項分別是2σ,2kσ;由於無法比較極值，我們必須在高斯空間繼續添加高斯模糊項，
使得形成σ,kσ,k2σ,k3σ,k4σ這樣就可以選擇DoG space中的中間三項kσ,k2σ,k3σ（只有左右都有才能有極值），
那么下一octave中（由上一層降采樣獲得）所得三項即為2kσ,2k2σ,2k3σ，其首項2kσ=2^4/3。剛好與上一octave末項k3σ=2^3/3尺度變化連續起來，
所以每次要在Gaussian space添加3項，每組（塔）共S+3層圖像，相應的DoG金字塔有S+2層圖像。
==============================================================================================================================
這個尺度因子都是在原圖上進行的。而在算法實現過程中，采用高斯平滑是在上一層圖像上再疊加高斯平滑。即我們在程序中看到的平滑因子為（4-5）
我們在源代碼上同一時候也沒有看到組間的2倍的關系，實際在對每組的平滑因子都是一樣的，2倍的關系是因為在降採樣的過程中產生的，第二層的第一張圖是由第一層的平滑因子為2σ的圖像上（即倒數第三張）降採樣得到，此時圖像平滑因子為σ，所以繼續採用以上的平滑因子。但在原圖上看。形成了所有的空間尺度。
這里一定要寫下參考大神的博格：
https://www.cnblogs.com/Alliswell-WP/p/SIFT_code1.html //這個是很多問題的解析,包括各種問題； https://blog.csdn.net/xw20084898/article/details/16832755 //這個參考了主要是尺度變換的連續性
3、求高斯差分金字塔（用DoG）
兩個高斯金字塔層相減生成一個差分金字塔。G（x, y, σ）已由上面算出
G（x, y, σ1）- G（x, y, σ2）
/************************************************************************************************************************* *模塊說明： * 圖像的差分 **************************************************************************************************************************/ void Sub(const Mat& a, const Mat& b, Mat & c) { if (a.rows != b.rows || a.cols != b.cols || a.type() != b.type()) return; if (!c.empty()) return; c.create(a.size(), a.type()); pixel_t* ap = (pixel_t*)a.data; pixel_t* ap_end = (pixel_t*)a.dataend; pixel_t* bp = (pixel_t*)b.data; pixel_t* cp = (pixel_t*)c.data; int step = a.step / sizeof(pixel_t); while (ap != ap_end) { *cp++ = *ap++ - *bp++; } } /************************************************************************************************************************* *模塊說明： * 模塊三：3.4 高斯差分金字塔 *功能說明： * 1--2002年，Mikolajczyk在詳細的實驗比較中發現尺度歸一化的高斯拉普拉斯函數的極大值和極小值同其他的特征提取函數，例如： * 梯度，Hessian或者Harris角點特征比較，能夠產生穩定的圖像特征。 * 2--而Lindberg早在1994年發現高斯差分函數（Difference of Gaussian,簡稱DOG算子）與尺度歸一化的拉普拉斯函數非常相似，因此， * 利用圖像間的差分代替微分。 **************************************************************************************************************************/ void DogPyramid(const vector<Mat>& gauss_pyr, vector<Mat>& dog_pyr, int octaves, int intervals = INTERVALS) { for (int o = 0; o < octaves; o++) { for (int i = 1; i < intervals + 3; i++) { Mat mat; Sub(gauss_pyr[o*(intervals + 3) + i], gauss_pyr[o*(intervals + 3) + i - 1], mat); dog_pyr.push_back(mat); } } }
總結：
高斯金字塔共分o組(Octave), 每組又分S層(Layer)。 組內各層圖像的分辨率是相同的，即長和寬相同，但尺度逐漸增加，即越往塔頂圖像越模糊。而下-組的圖像是由上一組圖像按照隔點降采樣得到的，即圖像的長和寬分別減半。
octaves = log((double)min(init_gray.rows, init_gray.cols)) / log(2.0) - 2; //計算高斯金字塔的層數
本實驗中o=6，S=s+3=6,小寫的s是高斯差分金字塔的層數；
但是我們在代碼中可以看到，每組中由6層，sub函數差分圖像后每組還剩5層；
那為什么我們說小寫的s（s=3）是高斯差分金字塔的層數；這就體現在后面的極值點檢測中；
4、極值點檢測
4.1、排除閾值小的點
4.2、判斷是否是極值
4.3、修正極值點，刪除不穩定的點
上述2步原理參考我的博文：SIFT算法原理（2）-極值點的精確定位
首先，第一個問題：
問題描述：如何找尋關鍵點？在幾層之間對比？思路是什么？
答：極值的檢測是在DoG空間進行的，檢測是以前點為中心，3pixel*3pixel*3pixel的立方體為鄰域，判斷當前點是否為局部最大或最小。如下圖所示，橘黃色為當前檢測點，綠色點為其鄰域。因為要比較當前點的上下層圖像，所以極值檢測從DoG每層的第2幅圖像開始，終止於每層的倒數第2幅圖像（第1幅沒有下層，最后1幅沒有上層，無法比較）。所以實際上我們極值點只在中間3層進行檢測
極值點初步檢測：
/*************************************************************************************************************************
*模塊說明：
* 模塊四：3.5 空間極值點的檢測(關鍵點的初步探查)
*功能說明：
* 1--關鍵點是由DOG空間的局部極值點組成的，關鍵點的初步探查是通過同一組內各DoG相鄰兩層圖像之間的比較完成的。為了尋找DoG
* 函數的極值點，每一個像素點都要和它所有相鄰的點比較，看其是否比它的圖像域和尺度域相鄰的點大還是小。
* 2--當然這樣產生的極值點並不全都是穩定的特征點，因為某些極值點相應較弱，而且DOG算子會產生較強的邊緣響應。
**************************************************************************************************************************/
void DetectionLocalExtrema(const vector<Mat>& dog_pyr, vector<Keypoint>& extrema, int octaves, int intervals = INTERVALS)
{
double thresh = 0.5 * DXTHRESHOLD / intervals; //|D(x^)| < 0.03 ; DXTHRESHOLD=0.03極值點太多
for (int o = 0; o < octaves; o++)
{
//第一層和最后一層極值忽略
for (int i = 1; i < (intervals + 2) - 1; i++) //Dog圖像每一組5層（差分6-1），又第一層和最后一程忽略；
{
int index = o*(intervals + 2) + i; //[1]圖片索引的定位
pixel_t *data = (pixel_t *)dog_pyr[index].data; //[2]獲取圖片的矩陣體的首地址
// std::cout << "首地址為 = " << sizeof(data[0]) << std::endl;
int step = dog_pyr[index].step / sizeof(data[0]); //[3]說明矩陣在存儲空間中的存儲是以線性空間的方式存放的
for (int y = IMG_BORDER; y < dog_pyr[index].rows - IMG_BORDER; y++) //IMG_BORDER=5 ，這里我查找資料說的是指檢測邊界以內的點，因為邊界易受干擾
{
for (int x = IMG_BORDER; x < dog_pyr[index].cols - IMG_BORDER; x++)
{
// std::cout << "x is = " << x << std::endl;
pixel_t val = *(data + y*step + x); //遍歷像素點
if (std::fabs(val) > thresh) //[4]排除閾值過小的點，fabs（）取絕對值
{
if (isExtremum(x, y, dog_pyr, index)) //[5]判斷是否是極值
{
Keypoint *extrmum = InterploationExtremum(x, y, dog_pyr, index, o, i);//極值點x,y；當前組當前層；
if (extrmum)
{
if (passEdgeResponse(extrmum->x, extrmum->y, dog_pyr, index))
{
extrmum->val = *(data + extrmum->y*step + extrmum->x);
extrema.push_back(*extrmum);
}
delete extrmum;
}//extrmum
}//isExtremum
}//std::fabs
}//for x
}//for y
}
}
}
isExtremum這個函數主要是判斷極值點，主要原理參考圖2；代碼如下：
/*************************************************************************************************************************
*模塊說明：
* 模塊四的第一步:3.4-空間極值點的初步檢測(關鍵點的初步探查)
*功能說明:
* 1--在高斯差分函數之后，得到一系列的關鍵點的疑似點，我們需要對這些關鍵點的疑似點初步進行檢測和篩選
* 檢測是以前點為中心，3pixel*3pixel*3pixel的立方體為鄰域，判斷當前點是否為局部最大或最小。i,j,k=-1,0,1表示當前點的鄰域
**************************************************************************************************************************/
bool isExtremum(int x, int y, const vector<Mat>& dog_pyr, int index)
{
pixel_t * data = (pixel_t *)dog_pyr[index].data;
int step = dog_pyr[index].step / sizeof(data[0]);
pixel_t val = *(data + y*step + x);
if (val > 0)
{
for (int i = -1; i <= 1; i++) //層
{
int stp = dog_pyr[index + i].step / sizeof(data[0]);
for (int j = -1; j <= 1; j++) //行
{
for (int k = -1; k <= 1; k++) //列
{
if (val < *((pixel_t*)dog_pyr[index + i].data + stp*(y + j) + (x + k)))
{
return false; //不是極值點退出
}
}
}
}
}
else
{
for (int i = -1; i <= 1; i++)
{
int stp = dog_pyr[index + i].step / sizeof(pixel_t);
for (int j = -1; j <= 1; j++)
{
for (int k = -1; k <= 1; k++)
{
//檢查最小極值
if (val > *((pixel_t*)dog_pyr[index + i].data + stp*(y + j) + (x + k)))
{
return false; //不是極值點退出
}
}
}
}
}
return true; //是級值點，返回
}
一個定義不好的高斯差分算子的極值在橫跨邊緣的地方有較大的主曲率，而在垂直邊緣的方向有較小的主曲率。DOG算子會產生較強的邊緣響應，需要剔除不穩定的邊緣響應點。獲取特征點處的Hessian矩陣，主曲率通過一個2x2 的Hessian矩陣H求出（D的主曲率和H的特征值成正比）；
原理參考：SIFT算法原理（2）-極值點的精確定位
代碼如下：
/*************************************************************************************************************************
*模塊說明：
* 模塊四的第三步:4.2--消除邊緣響應點
*功能說明:
* 1）一個定義不好的高斯差分算子的極值在橫跨邊緣的地方有較大的住主曲率，在垂直邊緣的方向有較小的主曲率。
* 2）DOG算子會產生較強的邊緣響應，需要剔除不穩定的邊緣響應點，獲取特征點處的Hessian矩陣，主曲率通過一個2*2的Hessian矩
* 陣H求出
* 3）主曲率D和Hessian矩陣的特征值成正比，公式(r+1)*(r+1)/r的值在兩個特征值相等時最小；這個值越大，說明兩個特征值的比值
* 越大，即在某一個方向的梯度值越大，而在另一個方向的梯度值越小，而邊緣恰恰就是這種情況。所以，為了剔除邊緣響應點，
* 需要讓該比值小於一定的閾值，因此，為了檢測主曲率是否在某閾值r下，只需檢測。CSDN論文中的公式(4-7成立)，成立，將關
* 鍵點保留，反之，剔除。
**************************************************************************************************************************/
#define DAt(x, y) (*(data+(y)*step+(x)))
bool passEdgeResponse(int x, int y, const vector<Mat>& dog_pyr, int index, double r = RATIO)
{
pixel_t *data = (pixel_t *)dog_pyr[index].data; //極值點的圖片索引定位
int step = dog_pyr[index].step / sizeof(data[0]);
pixel_t val = *(data + y*step + x); //極值點
double Dxx;
double Dyy;
double Dxy;
double Tr_h; //[1]Hessian矩陣的跡
double Det_h; //[2]Hessian矩陣所對應的行列式的值
Dxx = DAt(x + 1, y) + DAt(x - 1, y) - 2 * val;
Dyy = DAt(x, y + 1) + DAt(x, y - 1) - 2 * val;
Dxy = (DAt(x + 1, y + 1) + DAt(x - 1, y - 1) - DAt(x - 1, y + 1) - DAt(x + 1, y - 1)) / 4.0;
Tr_h = Dxx + Dyy;
Det_h = Dxx * Dyy - Dxy * Dxy;
if (Det_h <= 0)return false;
if (Tr_h * Tr_h / Det_h < (r + 1) * (r + 1) / r) return true;
return false;
}
5、計算極值點尺度，圖像特征縮放
/*************************************************************************************************************************
*模塊說明：
* 模塊五：計算特征點對應的尺度
*
**************************************************************************************************************************/
void CalculateScale(vector<Keypoint>& features, double sigma = SIGMA, int intervals = INTERVALS)
{
double intvl = 0;
for (int i = 0; i < features.size(); i++)
{
intvl = features[i].interval + features[i].offset_interval;
features[i].scale = sigma * pow(2.0, features[i].octave + intvl / intervals);
features[i].octave_scale = sigma * pow(2.0, intvl / intervals);
}
}
//對擴大的圖像特征縮放
/*將特征點序列中每個特征點的坐標減半(圖像金字塔設置了將圖像放大為原圖的2倍時，特征點檢測完之后調用)
/features：特征點序列*/
void HalfFeatures(vector<Keypoint>& features)
{
for (int i = 0; i < features.size(); i++)
{
features[i].dx /= 2;
features[i].dy /= 2;
features[i].scale /= 2;
}
}
關鍵點的x代表cols,列；關鍵點的y代表rows,行
features[i].scale = sigma * pow(2.0, features[i].octave + intvl / intervals);//尺度
features[i].octave_scale = sigma * pow(2.0, intvl / intervals);
6、關鍵點方向分配
6.1、計算梯度直方圖
6.2、對直方圖做兩次高斯平滑
6.3、求直方圖中的主方向
6.4、使求出的主方向更加精確
原理可參考：RobHess的SIFT代碼解析步驟三
主代碼：
/********************************************************************************************************************************
*模塊說明:
* 模塊六：5 關鍵點方向分配
*功能說明：
* 1）為了使描述符具有旋轉不變性，需要利用圖像的局部特征為每一個關鍵點分配一個基准方向。使用圖像梯度的方法求取局部結構的穩定
* 方向。
* 2）對於在DOG金字塔中檢測出來的關鍵點，采集其所在高斯金字塔圖像3sigma鄰域窗口內像素的梯度和方向梯度和方向特征。
* 3）梯度的模和方向如下所示:
* 4) 在完成關鍵點的梯度計算后，使用直方圖統計鄰域內像素的梯度和方向。梯度直方圖將0~360度的方向范圍分為36個柱，其中每柱10度，
* 如圖5.1所示，直方圖的峰值方向代表了關鍵點的主方向
cvRound()：返回跟參數最接近的整數值，即四舍五入；
*********************************************************************************************************************************/
void OrientationAssignment(vector<Keypoint>& extrema, vector<Keypoint>& features, const vector<Mat>& gauss_pyr)
{
int n = extrema.size();
double *hist;
for (int i = 0; i < n; i++)
{
//ORI_HIST_BINS=36;ORI_WINDOW_RADIUS=3.0*1.5,特征點方向賦值過程中，搜索鄰域的半徑為：3*1.5*σ
//ORI_SIGMA_TIMES=1.5,
hist = CalculateOrientationHistogram(gauss_pyr[extrema[i].octave*(INTERVALS + 3) + extrema[i].interval],
extrema[i].x, extrema[i].y, ORI_HIST_BINS, cvRound(ORI_WINDOW_RADIUS*extrema[i].octave_scale),
ORI_SIGMA_TIMES*extrema[i].octave_scale); //[1]計算梯度的方向直方圖
for (int j = 0; j < ORI_SMOOTH_TIMES; j++)
GaussSmoothOriHist(hist, ORI_HIST_BINS); //[2]對方向直方圖進行高斯平滑
double highest_peak = DominantDirection(hist, ORI_HIST_BINS); //[3]求取方向直方圖中的峰值
//[4]計算更加精確的關鍵點主方向
CalcOriFeatures(extrema[i], features, hist, ORI_HIST_BINS, highest_peak*ORI_PEAK_RATIO); //ORI_PEAK_RATIO=0.8
//精確關鍵點的主方向，拋物插值
delete[] hist;
}
}
調用代碼相關代碼如下：
問題描述：為什么要計算特征梯度直方圖？怎么計算？
答：在完成特征點鄰域的高斯圖像的梯度計算后，使用直方圖統計鄰域內像素的梯度方向和幅值。梯度方向直方圖的橫軸是梯度方向角，縱軸是梯度方向角對應的（帶高斯權重）梯度幅值累加值。梯度方向直方圖將。0°~360°的范圍，分為36個柱，每10°為一個柱。直方圖的峰值代表了該特征點處鄰域內圖像梯度的主方向，也即該特征點的主方向，
關鍵點的直方圖計算
/********************************************************************************************************************************
*模塊說明:
* 模塊六---步驟1：計算梯度的方向直方圖
*功能說明：
* 1）直方圖以每10度為一個柱，共36個柱，柱代表的方向為為像素點的梯度方向，柱的長短代表了梯度幅值。
* 2）根據Lowe的建議，直方圖統計采用3*1.5*sigma
* 3）在直方圖統計時，每相鄰三個像素點采用高斯加權，根據Lowe的建議，模板采用[0.25,0.5,0.25],並且連續加權兩次
*結 論：
* 圖像的關鍵點檢測完畢后，每個關鍵點就擁有三個信息：位置、尺度、方向；同時也就使關鍵點具備平移、縮放和旋轉不變性
bin=36.
*********************************************************************************************************************************/
double* CalculateOrientationHistogram(const Mat& gauss, int x, int y, int bins, int radius, double sigma)
{
//std::cout << "首地址1為 = " << radius << std::endl; //radius=3*1.5?*?σ(extrema[i].octave_scale)
//std::cout << "首地址2為 = " << bins << std::endl; // sigma=1.5*σ(extrema[i].octave_scale)
double* hist = new double[bins]; //[1]動態分配一個double類型的數組，用來存放梯度直方圖
for (int i = 0; i < bins; i++) //[2]給這個數組初始化
*(hist + i) = 0.0;
double mag; //[3]關鍵點的梯度幅值
double ori; //[4]關鍵點的梯度方向
double weight;
int bin;
const double PI2 = 2.0*CV_PI; //2Π
double econs = -1.0 / (2.0*sigma*sigma);
for (int i = -radius; i <= radius; i++)
{
for (int j = -radius; j <= radius; j++)
{
if (CalcGradMagOri(gauss, x + i, y + j, mag, ori)) //[5]計算該關鍵點的梯度幅值和方向
{
// std::cout << "ori為 = " << ori << std::endl;
weight = exp((i*i + j*j)*econs); // 該點的梯度幅值權重
bin = cvRound(bins * (CV_PI - ori) / PI2); //[6]對一個double行的數進行四舍五入，返回一個整形的數;
// bin = cvRound(bins * (CV_PI + ori) / PI2); //[6]對一個double行的數進行四舍五入，返回一個整形的數;
//計算梯度的方向對應的直方圖中的bin下標,ori范圍(-Π，Π]，+Π后范圍(0，2Π]
// std::cout << "bin1為 = " << bin << std::endl;
bin = bin < bins ? bin : 0; //把bin的下標歸到0~36；bin<bins，bin=bin;否則bin=0；
//std::cout << "bin2為 = " << bin << std::endl;
hist[bin] += mag * weight; //[7]統計梯度的方向直方圖
}
}
}
return hist;
}
這個公式主要實在計算特征值的方向和幅值，體現在CalcGradMagOri（）函數上，代碼如下：
/********************************************************************************************************************************
*模塊說明:
* 模塊六---步驟2：計算關鍵點的梯度和梯度方向
*功能說明：
* 1）計算關鍵點(x,y)處的梯度幅值和梯度方向
* 2）將所計算出來的梯度幅值和梯度方向保存在變量mag和ori中
*********************************************************************************************************************************/
bool CalcGradMagOri(const Mat& gauss, int x, int y, double& mag, double& ori)
{
if (x > 0 && x < gauss.cols - 1 && y > 0 && y < gauss.rows - 1)
{
pixel_t *data = (pixel_t*)gauss.data;
int step = gauss.step / sizeof(*data);
double dx = *(data + step*y + (x + 1)) - (*(data + step*y + (x - 1))); //[1]利用X方向上的差分代替微分dx
double dy = *(data + step*(y + 1) + x) - (*(data + step*(y - 1) + x)); //[2]利用Y方向上的差分代替微分dy
mag = sqrt(dx*dx + dy*dy); //[3]計算該關鍵點的梯度幅值
ori = atan2(dy, dx); //[4]計算該關鍵點的梯度方向
return true;
}
else
return false;
}
atan2相關知識匯總——https://blog.csdn.net/weixin_42142612/article/details/80972768
對直方圖做兩次高斯平滑
LOWE建議對直方圖進行平滑，降低突變的影響，彌補因沒有仿射不變性而產生的特征點不穩定的問題
hist[0]=0.25*hist[35]+0.5*hist[0]+0.25*hist[1];
/********************************************************************************************************************************
*模塊說明:
* 模塊六---步驟3：對梯度方向直方圖進行連續兩次的高斯平滑
*功能說明：
* 1）在直方圖統計時，每相鄰三個像素點采用高斯加權，根據Lowe的建議，模板采用[0.25,0.5,0.25],並且連續加權兩次
* 2）對直方圖進行兩次平滑
*********************************************************************************************************************************/
void GaussSmoothOriHist(double *hist, int n)
{
double prev = hist[n - 1];
double temp;
double h0 = hist[0];
for (int i = 0; i < n; i++)
{
temp = hist[i];
hist[i] = 0.25 * prev + 0.5 * hist[i] + 0.25 * (i + 1 >= n ? h0 : hist[i + 1]);//對方向直方圖進行高斯平滑
prev = temp;
}
}
6.3、求直方圖中的主方向
6.4、使求出的主方向更加精確
方向直方圖的峰值則代表了該特征點處鄰域梯度的方向，以直方圖中最大值作為該關鍵點的主方向。為了增強匹配的魯棒性，只保留峰值大於主方向峰值80％的方向（另一個方向的直方圖的幅值大於主方向峰值80％的方向）作為該關鍵點的輔方向。
/********************************************************************************************************************************
*模塊說明:
* 模塊六---步驟5：計算更加精確的關鍵點主方向----拋物插值
*功能說明：
* 1）方向直方圖的峰值則代表了該特征點的方向，以直方圖中的最大值作為該關鍵點的主方向。為了增強匹配的魯棒性，只保留峰值大於主
* 方向峰值80%的方向作為改關鍵點的輔方向。因此，對於同一梯度值得多個峰值的關鍵點位置，在相同位置和尺度將會有多個關鍵點被
* 創建但方向不同。僅有15%的關鍵點被賦予多個方向，但是可以明顯的提高關鍵點的穩定性。
* 2）在實際編程中，就是把該關鍵點復制成多份關鍵點，並將方向值分別賦給這些復制后的關鍵點
* 3）並且，離散的梯度直方圖要進行【插值擬合處理】，來求得更加精確的方向角度值
*********************************************************************************************************************************/
#define Parabola_Interpolate(l, c, r) (0.5*((l)-(r))/((l)-2.0*(c)+(r))) //求小柱子在直方圖中的索引號
void CalcOriFeatures(const Keypoint& keypoint, vector<Keypoint>& features, const double *hist, int n, double mag_thr)
{
double bin;
double PI2 = CV_PI * 2.0;
int l;
int r;
//遍歷直方圖
for (int i = 0; i < n; i++)
{
l = (i == 0) ? n - 1 : i - 1; //前一個(左邊的)bin的下標，
r = (i + 1) % n; //后一個(右邊的)bin的下標
//i=0,l=35,r=1;
//hist[i]是極值
//若當前的bin是局部極值(比前一個和后一個bin都大)，並且值大於給定的幅值閾值，則新生成一個特征點並添加到特征點序列末尾
if (hist[i] > hist[l] && hist[i] > hist[r] && hist[i] >= mag_thr) //mag_thr=0.8*關鍵點最大值
{
//根據左、中、右三個bin的值對當前bin進行直方圖插值
bin = i + Parabola_Interpolate(hist[l], hist[i], hist[r]);
bin = (bin < 0) ? (bin + n) : (bin >= n ? (bin - n) : bin); //將插值結果規范到[0,n]內
Keypoint new_key;
CopyKeypoint(keypoint, new_key); //復制特征點
new_key.ori = ((PI2 * bin) / n) - CV_PI; //方向還原為（-Π，Π]，bin最大值36，2Π*36/36-Π=Π，bin最小值0，2Π*0/36-Π=-Π
features.push_back(new_key); //插入到特征點序列末尾
}
}
}
//復制關鍵點
void CopyKeypoint(const Keypoint& src, Keypoint& dst)
{
dst.dx = src.dx;
dst.dy = src.dy;
dst.interval = src.interval;
dst.octave = src.octave;
dst.octave_scale = src.octave_scale;
dst.offset_interval = src.offset_interval;
dst.offset_x = src.offset_x;
dst.offset_y = src.offset_y;
dst.ori = src.ori;
dst.scale = src.scale;
dst.val = src.val;
dst.x = src.x;
dst.y = src.y;
}
再次討論一下features
關鍵點的x代表cols,列；關鍵點的y代表rows,行 features[i].scale = sigma * pow(2.0, features[i].octave + intvl / intervals);//尺度 features[i].octave_scale = sigma * pow(2.0, intvl / intervals); 梯度幅值和梯度方向保存在變量mag和ori中
7、關鍵點描述
7.1、確定描述子所需的鄰域區域
問題1：所需的圖像區域半徑radius怎么計算？
答：將關鍵點附近的區域划分為d*d(Lowe建議d=4)個子區域，每個子區域作為一個種子點，每個種子點有8個方向。考慮到實際計算時，需要采用三線性插值，所需圖像窗口邊長為3x3xσ_oct x(d+1) 。在考慮到旋轉因素(方便下一步將坐標軸旋轉到關鍵點的方向)，如下圖所示，實際計算所需的圖像區域半徑為：
部分代碼：
//為了保證特征描述子具有旋轉不變性，要以特征點為中心，在附近鄰域內旋轉θ角，即旋轉為特征點的方向 double cos_ori = cos(ori); double sin_ori = sin(ori); //6.1高斯權值，sigma等於描述字窗口寬度的一半 double sigma = 0.5 * width; double conste = -1.0 / (2 * sigma*sigma); double PI2 = CV_PI * 2; //計算特征描述子過程中，特征點周圍的width*width個區域中，每個區域的寬度為3*σ個像素， double sub_hist_width = DESCR_SCALE_ADJUST * octave_scale; //DESCR_SCALE_ADJUST=3 //【1】計算描述子所需的圖像領域區域的半徑 //考慮到要進行雙線性插值，每個區域的寬度應為:DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 ) //在考慮到旋轉因素，每個區域的寬度應為：DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )* sqrt(2) //所以搜索的半徑是：DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )* sqrt(2) / 2 int radius = (sub_hist_width*sqrt(2.0)*(width + 1)) / 2.0 + 0.5; //[1]0.5取四舍五入
7.2、（依據6求出的方向）旋轉坐標軸，使具有旋轉不變性
將坐標軸旋轉為關鍵點的方向，以確保旋轉不變性。
/遍歷每個區域的像素
for (int i = -radius; i <= radius; i++)
{
for (int j = -radius; j <= radius; j++)
{
//坐標旋轉為主方向,以確保旋轉不變性。
double rot_x = (cos_ori * j - sin_ori * i) / sub_hist_width;
double rot_y = (sin_ori * j + cos_ori * i) / sub_hist_width;
//將鄰域內的采樣點分配到對應的子區域內，將子區域內的梯度值分配到8個方向上，計算其權值
double xbin = rot_x + width / 2 - 0.5; //[2]xbin,ybin為落在4*4窗口中的下標值
double ybin = rot_y + width / 2 - 0.5;
7.3、將鄰域內采樣點分配到旋轉后的對應子區域，確認種子點
if (xbin > -1.0 && xbin < width && ybin > -1.0 && ybin < width)
{
if (CalcGradMagOri(gauss, x + j, y + i, grad_mag, grad_ori)) //[3]計算關鍵點的梯度方向
{
grad_ori = (CV_PI - grad_ori) - ori;
while (grad_ori < 0.0)
grad_ori += PI2;
while (grad_ori >= PI2)
grad_ori -= PI2; //把梯度方向歸入[0，2Π]
double obin = grad_ori * (bins / PI2);
//梯度方向*(8/2Π)歸到[0，8]，8為直方圖有8個bin，2Π為角度取值范圍的長度，把方向的索引歸於0~8.
//conste = -1.0 / (2 * sigma*sigma);
//計算權重
double weight = exp(conste*(rot_x*rot_x + rot_y * rot_y));
7.4、計算種子點的8個方向的梯度信息
三線性差值：
/*
功能說明：
將一個條目插入到形成特征描述符的方向直方圖數組中，該條目分配最多8個方向，對於每個維度每個方向乘以（1-d） 的權重，
其中d是以bin測量的bin中心值距離。
//hist為方向直方圖的三維數組
//xbin子行
//ybin子列
//obin子方向
//mag權重
//d方向直方圖的寬度
//bins每個直方圖中bin的個數
三線性插值
*/
void InterpHistEntry(double ***hist, double xbin, double ybin, double obin, double mag, int bins, int d)
{
double d_r, d_c, d_o, v_r, v_c, v_o;
double** row, *h;
int r0, c0, o0, rb, cb, ob, r, c, o;
r0 = cvFloor(ybin); //向下取整
c0 = cvFloor(xbin);
o0 = cvFloor(obin);
d_r = ybin - r0; //小數余項
d_c = xbin - c0;
d_o = obin - o0;
//這里也可以看出前面rbin，cbin為何減0.5，這樣原點上這點d_r,d_c均為0.5，即原點上這點的方向將被平均分配疊加在它周圍4個直方圖上面。
/*
做插值：
xbin,ybin,obin:種子點所在子窗口的位置和方向
所有種子點都將落在4*4的窗口中
r0,c0取不大於xbin，ybin的正整數
r0,c0只能取到0,1,2
xbin,ybin在(-1, 2)
r0取不大於xbin的正整數時。
r0+0 <= xbin <= r0+1
mag在區間[r0,r1]上做插值
obin同理
*/
for (r = 0; r <= 1; r++) // 沿着行方向不超過1個單位長度
{
rb = r0 + r;
if (rb >= 0 && rb < d) //如果此刻還在真正的描述子區間內
{
v_r = mag * ((r == 0) ? 1.0 - d_r : d_r); //對近鄰兩行的貢獻因子
row = hist[rb]; //第rb行上的hist
for (c = 0; c <= 1; c++) //沿着列方向不超過1個單位長度
{
cb = c0 + c;
if (cb >= 0 && cb < d) //如果此刻還在真正的描述子區間內
{
v_c = v_r * ((c == 0) ? 1.0 - d_c : d_c); //對近鄰兩行的貢獻因子
h = row[cb]; //h表示第rb行cb列上的hist
for (o = 0; o <= 1; o++) //沿着直方圖方向不超過1個單位長度
{
ob = (o0 + o) % bins; //bins=8，8個柱子
v_o = v_c * ((o == 0) ? 1.0 - d_o : d_o); //對近鄰兩個方向的貢獻因子
h[ob] += v_o;
}
}
}
}
}
}
這里三線性差值還不是很懂，后續再繼續學習
前面幾個步驟的總代碼：
/********************************************************************************************************************************
*模塊說明:
* 模塊七--步驟1:計算描述子的直方圖
*功能說明：
gauss：圖像指針
x：特征點所在的行
y：特征點所在的列
ori：特征點的主方向
octave_scale：特征點的尺度//特征點所在組的尺度
width：計算方向直方圖時，將特征點附近划分為width*width個區域，每個區域生成一個直方圖，默認width為4
bins：每個直方圖中bins的個數
返回值：double類型的三維數組，即一個d*d的二維數組，數組中每個元素是一個有n個bins的直方圖數組
*********************************************************************************************************************************/
double*** CalculateDescrHist(const Mat& gauss, int x, int y, double octave_scale, double ori, int bins, int width)
{
double ***hist = new double**[width]; //width*width*bins的三維直方圖數組 //為第一維分配空間
for (int i = 0; i < width; i++)
{
hist[i] = new double*[width]; //為第二維分配空間
for (int j = 0; j < width; j++)
{
hist[i][j] = new double[bins]; //為第三維分配空間
}
}
for (int r = 0; r < width; r++)
for (int c = 0; c < width; c++)
for (int o = 0; o < bins; o++)
hist[r][c][o] = 0.0;
//為了保證特征描述子具有旋轉不變性，要以特征點為中心，在附近鄰域內旋轉θ角，即旋轉為特征點的方向
double cos_ori = cos(ori);
double sin_ori = sin(ori);
//6.1高斯權值，sigma等於描述字窗口寬度的一半
double sigma = 0.5 * width;
double conste = -1.0 / (2 * sigma*sigma);
double PI2 = CV_PI * 2;
//計算特征描述子過程中，特征點周圍的width*width個區域中，每個區域的寬度為3*σ個像素，
double sub_hist_width = DESCR_SCALE_ADJUST * octave_scale; //DESCR_SCALE_ADJUST=3
//【1】計算描述子所需的圖像領域區域的半徑
//考慮到要進行雙線性插值，每個區域的寬度應為:DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )
//在考慮到旋轉因素，每個區域的寬度應為：DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )* sqrt(2)
//所以搜索的半徑是：DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )* sqrt(2) / 2
int radius = (sub_hist_width*sqrt(2.0)*(width + 1)) / 2.0 + 0.5; //[1]0.5取四舍五入
double grad_ori;
double grad_mag;
//遍歷每個區域的像素
for (int i = -radius; i <= radius; i++)
{
for (int j = -radius; j <= radius; j++)
{
//坐標旋轉為主方向,以確保旋轉不變性。
double rot_x = (cos_ori * j - sin_ori * i) / sub_hist_width;
double rot_y = (sin_ori * j + cos_ori * i) / sub_hist_width;
//將鄰域內的采樣點分配到對應的子區域內，將子區域內的梯度值分配到8個方向上，計算其權值
double xbin = rot_x + width / 2 - 0.5; //[2]xbin,ybin為落在4*4窗口中的下標值
double ybin = rot_y + width / 2 - 0.5;
if (xbin > -1.0 && xbin < width && ybin > -1.0 && ybin < width)
{
if (CalcGradMagOri(gauss, x + j, y + i, grad_mag, grad_ori)) //[3]計算關鍵點的梯度方向
{
grad_ori = (CV_PI - grad_ori) - ori;
while (grad_ori < 0.0)
grad_ori += PI2;
while (grad_ori >= PI2)
grad_ori -= PI2; //把梯度方向歸入[0，2Π]
double obin = grad_ori * (bins / PI2);
//梯度方向*(8/2Π)歸到[0，8]，8為直方圖有8個bin，2Π為角度取值范圍的長度，把方向的索引歸於0~8.
//conste = -1.0 / (2 * sigma*sigma);
//計算權重
double weight = exp(conste*(rot_x*rot_x + rot_y * rot_y));
InterpHistEntry(hist, xbin, ybin, obin, grad_mag*weight, bins, width); //線性差值求每個種子點八個方向的梯度。
}
}
}
}
return hist;
}
7.5、將8個方向的梯度變為特征向量
7.6、歸一化消除光照影響
7.7、描述子門限，再次歸一化
特征描述子
如上統計的4*4*8=128個梯度信息即為該關鍵點的特征向量。
特征向量形成后，為了去除光照變化的影響，需要對它們進行歸一化處理，對於圖像灰度值整體漂移，圖像各點的梯度是鄰域像素相減得到，所以也能去除。得到的描述子向量為H=(h1,h2,.......,h128)，歸一化后的特征向量為L=(L1,L2,......,L128)，則
描述子的門限化
非線性光照，相機飽和度變化對造成某些方向的梯度值過大，而對方向的影響微弱。因此設置門限值(向量歸一化后，一般取0.2)截斷較大的梯度值(大於0.2的則就令它等於0.2，小於0.2的則保持不變)。然后再進行一次歸一化處理，提高特征的鑒別性。
//歸一化特征點的特征描述子，即將特征描述子數組中每個元素除以特征描述子的模
void NormalizeDescr(Keypoint& feat)
{
double cur, len_inv, len_sq = 0.0;
int i, d = feat.descr_length; //特征描述子的維數
//求特征描述子的模
for (i = 0; i < d; i++)
{
cur = feat.descriptor[i];
len_sq += cur*cur;
}
len_inv = 1.0 / sqrt(len_sq); //sqrt( len_sq )為特征描述子的模
for (i = 0; i < d; i++)
feat.descriptor[i] *= len_inv; //特征描述子中每個元素除以特征描述子的模，完成歸一化
}
/********************************************************************************************************************************
*模塊說明:
* 模塊七--步驟2:直方圖到描述子的轉換
*功能說明：
將某特征點的方向直方圖轉換為特征描述子向量，對特征描述子歸一化並將所有元素轉化為整型，存入指定特征點中
參數：
hist：width*width*bins的三維直方圖數組
width：計算方向直方圖時，將特征點附近划分為width*width個區域，每個區域生成一個直方圖
bins：每個直方圖的bin個數
feature：特征點指針，將計算好的特征描述子存入其中
*********************************************************************************************************************************/
void HistToDescriptor(double ***hist, int width, int bins, Keypoint& feature)
{
int int_val, i, r, c, o, k = 0;
//遍歷width*width*bins的三維直方圖數組，將其中的所有數據(一般是128個)都存入feature結構的descriptor成員中
for (r = 0; r < width; r++)
for (c = 0; c < width; c++)
for (o = 0; o < bins; o++)
{
feature.descriptor[k++] = hist[r][c][o];
}
feature.descr_length = k; //特征描述子的維數，一般是128
NormalizeDescr(feature); //[1]描述子特征向量歸一化
//歸一化特征點的特征描述子，即將特征描述子數組中每個元素除以特征描述子的模
// DESCR_MAG_THR 0.2
//遍歷特征描述子向量，將超過閾值DESCR_MAG_THR的元素強行賦值DESCR_MAG_THR
for (i = 0; i < k; i++) //[2]描述子向量門限
if (feature.descriptor[i] > DESCR_MAG_THR)
feature.descriptor[i] = DESCR_MAG_THR;
NormalizeDescr(feature); //[3]描述子進行最后一次的歸一化操作
//遍歷特征描述子向量，每個元素乘以系數INT_DESCR_FCTR( 512.0)來變為整型，並且最大值不能超過255
for (i = 0; i < k; i++) //[4]將單精度浮點型的描述子轉換為整形的描述子
{
int_val = INT_DESCR_FCTR * feature.descriptor[i];
feature.descriptor[i] = min(255, int_val);
}
}
上述3步原理參考我的博文：SIFT算法原理（3）-確定關鍵點的主方位，構建關鍵點描述符
如何畫出已找到的關鍵點：
一些畫出點的代碼信息：
opencv數據結構CvScalar，【雜談opencv】OpenCV中的cvRound()、cvFloor()、 cvCeil()函數講解，
opencv學習中——CvPoint、CvSize、CvRect、CV_RGB、cvRectangle
代碼實現：
/********************************************************************************************************************************
*模塊說明:
* 模塊七:6 關鍵點描述
*功能說明：
* 1）通過以上步驟，對於一個關鍵點，擁有三個信息：位置、尺度、方向
* 2）接下來就是為每個關鍵點建立一個描述符，用一組向量來將這個關鍵點描述出來，使其不隨各種變化而變化，比如光照、視角變化等等
* 3）這個描述子不但包括關鍵點，也包含關鍵點周圍對其貢獻的像素點，並且描述符應該有較高的獨特性，以便於特征點正確的匹配概率
* 1）SIFT描述子----是關鍵點鄰域高斯圖像梯度統計結果的一種表示。
* 2）通過對關鍵點周圍圖像區域分塊，計算塊內梯度直方圖，生成具有獨特性的向量
* 3）這個向量是該區域圖像信息的一種表述和抽象，具有唯一性。
*Lowe論文：
* Lowe建議描述子使用在關鍵點尺度空間內4*4的窗口中計算的8個方向的梯度信息，共4*4*8=128維向量來表征。具體的步驟如下所示:
* 1)確定計算描述子所需的圖像區域
* 2）將坐標軸旋轉為關鍵點的方向，以確保旋轉不變性，如CSDN博文中的圖6.2所示；旋轉后鄰域采樣點的新坐標可以通過公式(6-2)計算
* 3）將鄰域內的采樣點分配到對應的子區域，將子區域內的梯度值分配到8個方向上，計算其權值
* 4）插值計算每個種子點八個方向的梯度
* 5）如上統計的4*4*8=128個梯度信息即為該關鍵點的特征向量。特征向量形成后，為了去除光照變化的影響，需要對它們進行歸一化處理，
* 對於圖像灰度值整體漂移，圖像各點的梯度是鄰域像素相減得到的，所以也能去除。得到的描述子向量為H，歸一化之后的向量為L
* 6）描述子向量門限。非線性光照，相機飽和度變化對造成某些方向的梯度值過大，而對方向的影響微弱。因此，設置門限值（向量歸一化
* 后，一般取0.2）截斷較大的梯度值。然后，在進行一次歸一化處理，提高特征的鑒別性。
* 7）按特征點的尺度對特征描述向量進行排序
* 8）至此，SIFT特征描述向量生成。
將關鍵點附近的區域划分為d*d(Lowe建議d=4)個子區域，每個子區域作為一個種子點，每個種子點有8個方向。
*********************************************************************************************************************************/
void DescriptorRepresentation(vector<Keypoint>& features, const vector<Mat>& gauss_pyr, int bins, int width) //bins=8,width=4
{
double ***hist;
for (int i = 0; i < features.size(); i++)
{ //[1]計算描述子的直方圖
//gauss_pyr[features[i].octave*(INTERVALS + 3) + features[i].interval 計算特征點所在矩陣中的位置；
hist = CalculateDescrHist(gauss_pyr[features[i].octave*(INTERVALS + 3) + features[i].interval],
features[i].x, features[i].y, features[i].octave_scale, features[i].ori, bins, width);
HistToDescriptor(hist, width, bins, features[i]); //[2]直方圖到描述子的轉換
/*釋放計算特征描述子過程中用到的方向直方圖的內存空間
hist：方向直方圖的指針，是一個width*width*bins的三維直方圖數組
*/
for (int j = 0; j < width; j++)
{
for (int k = 0; k < width; k++)
{
delete[] hist[j][k]; //釋放第三維的內存
}
delete[] hist[j]; //釋放第二維的內存
}
delete[] hist; //釋放第一維的內存
}
}
上面多維數組可參考：圖解c/c＋＋多級指針與“多維”數組
排序函數用法：
sort函數的用法(C++排序庫函數的調用)
sift代碼調用總：
/*******************************************************************************************************************
*函數說明:
* 最大的模塊1：SIFT算法模塊
*函數參數說明:
* 1---const Mat &src---------------准備進行特征點檢測的原始圖片
* 2---Vector<Keypoint>& features---用來存儲檢測出來的關鍵點
* 3---double sigma-----------------sigma值
* 4---int intervals----------------關鍵點所在的層數
********************************************************************************************************************/
void Sift(const Mat &src, vector<Keypoint>& features, double sigma, int intervals)
{
std::cout << "【Step_one】Create -1 octave gaussian pyramid image" << std::endl;
std::cout << "sigma = " << sigma << std::endl;
std::cout << "intervals = " << intervals << std::endl;
cv::Mat init_gray;
CreateInitSmoothGray(src, init_gray, sigma); //灰度圖
int octaves = log((double)min(init_gray.rows, init_gray.cols)) / log(2.0) - 2; //計算高斯金字塔的層數
std::cout << "【1】The height and width of init_gray_img = " << init_gray.rows << "*" << init_gray.cols << std::endl;
std::cout << "【2】The octaves of the gauss pyramid = " << octaves << std::endl;
std::cout << "【Step_two】Building gaussian pyramid ..." << std::endl;
std::vector<Mat> gauss_pyr;
GaussianPyramid(init_gray, gauss_pyr, octaves, intervals, sigma); //高斯金字塔
write_pyr(gauss_pyr, "gausspyramid");
std::cout << "【Step_three】Building difference of gaussian pyramid..." << std::endl;
vector<Mat> dog_pyr;
DogPyramid(gauss_pyr, dog_pyr, octaves, intervals); //差分金字塔
write_pyr(dog_pyr, "dogpyramid");
std::cout << "【Step_four】Deatecting local extrema..." << std::endl;
vector<Keypoint> extrema;
DetectionLocalExtrema(dog_pyr, extrema, octaves, intervals); //極值點初探（排除較小點，找到極值點，消除邊緣響應）
std::cout << "【3】keypoints cout: " << extrema.size() << " --" << std::endl;
std::cout << "【4】extrema detection finished." << std::endl;
std::cout << "【5】please look dir gausspyramid, dogpyramid and extrema.txt.--" << endl;
std::cout << "【Step_five】CalculateScale..." << std::endl;
CalculateScale(extrema, sigma, intervals); //計算特征點的尺度
HalfFeatures(extrema); //圖像縮放
std::cout << "【Step_six】Orientation assignment..." << std::endl;
OrientationAssignment(extrema, features, gauss_pyr); //關鍵點直方圖計算，使用的extrema是縮放后的圖像
std::cout << "【6】features count: " << features.size() << std::endl;
std::cout << "【Step_seven】DescriptorRepresentation..." << std::endl;
//DESCR_HIST_BINS=8;DESCR_WINDOW_WIDTH=4;
DescriptorRepresentation(features, gauss_pyr, DESCR_HIST_BINS, DESCR_WINDOW_WIDTH); //關鍵點主方向的描述
sort(features.begin(), features.end(), FeatureCmp); //排序， FeatureCmp是排序方法
cout << "finished." << endl;
}
最后放一下總的代碼
vs2015+0pencv3.2.0實現
sift.cpp
#include"sift.h"
#include<fstream>
#include<iostream>
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <opencv2/core/core.hpp>
#include <stdio.h>
#include <stdlib.h>
using namespace std;
using namespace cv;
/*************************************************************************************************************************
*模塊說明：
* 主函數；
**************************************************************************************************************************/
int main(int argc, char **argv)
{
cv::Mat src = imread("E:\\VS2015Opencv\\vs2015\\project\\picture\\05.jpg");
if (src.empty())
{
cout << "jpg open error! " << endl;
getchar();
return -1;
}
if (src.channels() != 3) return -2;
vector<Keypoint> features;
Sift(src, features, 1.6); //【1】SIFT特征點檢測和特征點描述
DrawKeyPoints(src, features); //【2】畫出關鍵點(特征點)
DrawSiftFeatures(src, features); //【3】畫出SIFT特征點
write_features(features, "descriptor.txt"); //【4】將特征點寫入文本
cv::imshow("src", src);
cv::waitKey();
return 0;
}
/*************************************************************************************************************************
*模塊說明：
* 圖像灰度化函數----將彩色圖像轉為灰度圖像
**************************************************************************************************************************/
void ConvertToGray(const Mat& src, Mat& dst)
{
cv::Size size = src.size();
if (dst.empty())
dst.create(size, CV_64F); //[1]利用Mat類的成員函數創建Mat容器
uchar* srcData = src.data; //[2]指向存儲所有像素值的矩陣的數據區域
pixel_t* dstData = (pixel_t*)dst.data; //[3]指向dst的矩陣體數據區域
int dstStep = dst.step / sizeof(dstData[0]); //[4]一行圖像含有多少字節數
for (int j = 0; j<src.cols; j++)
{
for (int i = 0; i<src.rows; i++)
{
double b = (srcData + src.step*i + src.channels()*j)[0] / 255.0;
double g = (srcData + src.step*i + src.channels()*j)[1] / 255.0;
double r = (srcData + src.step*i + src.channels()*j)[2] / 255.0;
*(dstData + dstStep*i + dst.channels()*j) = (r + g + b) / 3.0;
}
}
}
/*************************************************************************************************************************
*模塊說明：
* 隔點采樣
**************************************************************************************************************************/
void DownSample(const Mat& src, Mat& dst)
{
if (src.channels() != 1)
return;
if (src.cols <= 1 || src.rows <= 1)
{
src.copyTo(dst);
return;
}
dst.create((int)(src.rows / 2), (int)(src.cols / 2), src.type()); //dst的行列=src/2
pixel_t* srcData = (pixel_t*)src.data;
pixel_t* dstData = (pixel_t*)dst.data;
int srcStep = src.step / sizeof(srcData[0]); //src一行圖像含有多少字節數
int dstStep = dst.step / sizeof(dstData[0]);
int m = 0, n = 0;
for (int j = 0; j < src.cols; j += 2, n++)
{
m = 0;
for (int i = 0; i < src.rows; i += 2, m++)
{
pixel_t sample = *(srcData + srcStep * i + src.channels() * j);
if (m < dst.rows && n < dst.cols)
{
*(dstData + dstStep * m + dst.channels() * n) = sample; //將src隔點采樣到dst中去。剛好src行列/2
}
}
}
}
/*************************************************************************************************************************
*模塊說明：
* 線性插值放大，將圖像放大2倍
**************************************************************************************************************************/
void UpSample(const Mat &src, Mat &dst)
{
if (src.channels() != 1)
return;
dst.create(src.rows * 2, src.cols * 2, src.type()); //dst的行列是原圖2倍大小的行列
pixel_t* srcData = (pixel_t*)src.data;
pixel_t* dstData = (pixel_t*)dst.data;
int srcStep = src.step / sizeof(srcData[0]); //一行圖像含有多少字節數
int dstStep = dst.step / sizeof(dstData[0]);
int m = 0, n = 0;
for (int j = 0; j < src.cols - 1; j++, n += 2)
{
m = 0;
for (int i = 0; i < src.rows - 1; i++, m += 2)
{
/*第1步：n=0 m=0，n=2,m=0; 第2步：n=0 m=1，n=2,m=1; 第3步：n=1 m=0，n=3,m=0; 第4步：n=1 m=1，n=3,m=1; 第5步：最后2行2列 */
/*第1步：n=0 m=2，n=2,m=2; 第2步：n=0 m=3，n=2,m=3; 第3步：n=1 m=2，n=3,m=2; 第4步：n=1 m=3，n=3,m=3; 假設5*5的矩陣 */
double sample = *(srcData + srcStep * i + src.channels() * j);
*(dstData + dstStep * m + dst.channels() * n) = sample;
double rs = *(srcData + srcStep * (i)+src.channels()*j) + (*(srcData + srcStep * (i + 1) + src.channels()*j));
*(dstData + dstStep * (m + 1) + dst.channels() * n) = rs / 2;
double cs = *(srcData + srcStep * i + src.channels()*(j)) + (*(srcData + srcStep * i + src.channels()*(j + 1)));
*(dstData + dstStep * m + dst.channels() * (n + 1)) = cs / 2;
double center = (*(srcData + srcStep * (i + 1) + src.channels() * j))
+ (*(srcData + srcStep * i + src.channels() * j))
+ (*(srcData + srcStep * (i + 1) + src.channels() * (j + 1)))
+ (*(srcData + srcStep * i + src.channels() * (j + 1)));
*(dstData + dstStep * (m + 1) + dst.channels() * (n + 1)) = center / 4;
}
}
if (dst.rows < 3 || dst.cols < 3)
return;
//最后兩行兩列
for (int k = dst.rows - 1; k >= 0; k--)
{
*(dstData + dstStep *(k)+dst.channels()*(dst.cols - 2)) = *(dstData + dstStep *(k)+dst.channels()*(dst.cols - 3));
*(dstData + dstStep *(k)+dst.channels()*(dst.cols - 1)) = *(dstData + dstStep *(k)+dst.channels()*(dst.cols - 3));
}
for (int k = dst.cols - 1; k >= 0; k--)
{
*(dstData + dstStep *(dst.rows - 2) + dst.channels()*(k)) = *(dstData + dstStep *(dst.rows - 3) + dst.channels()*(k));
*(dstData + dstStep *(dst.rows - 1) + dst.channels()*(k)) = *(dstData + dstStep *(dst.rows - 3) + dst.channels()*(k));
}
}
/*************************************************************************************************************************
*模塊說明：
* OpenCv中的高斯平滑函數
**************************************************************************************************************************/
void GaussianSmooth(const Mat &src, Mat &dst, double sigma)
{
GaussianBlur(src, dst, Size(0, 0), sigma, sigma);
}
/*************************************************************************************************************************
*模塊說明：
* 模塊1--創建初始灰度圖像------初始圖像先將原圖像灰度化，再擴大一倍后，使用高斯模糊平滑
*功能說明：
* 在最開始建立高斯金字塔的時候，要預先模糊輸入的圖像來作為第0個組的第0層圖像，這時相當於丟棄了最高的空域的采樣率。因
* 此，通常的做法是先將圖像的尺度擴大一倍來生成第-1組。我們假定初始的輸入圖像為了抗擊混淆現象，已經對其進行了sigma=0.5
* 的高斯模糊，如果輸入圖像的尺寸用雙線性插值擴大一倍，那么相當於sigma=1.0
*這樣做的原因:
* 在檢測極值點前，對原始圖像的高斯平滑以致圖像丟失高頻信息，所以Lowe建議在建立尺度空間前，首先對原始圖像長寬擴展一倍，
* 以便保留原始圖像的信息(特別是圖像的高頻信息，比如邊緣，角點)，增加特征點的數量。
*函數功能:
* 這個函數的作用是用於創建高斯金字塔的-1層圖像，在主函數中定義 sigma =1.6
在lowe的論文中，他定義圖片的尺度是1.6，被攝像頭模糊的尺度是0.5，但是這個1.52對我們來說並沒有啥意思，
lowe為了獲得更多的特征點（姑且這么解釋吧）他先對原始圖像進行擴大一倍。然后呢再對它進行高斯平滑作為第一階第一層，
高斯模糊尺度為=(1.6* 1.6 - (0.5 * 2) * (0.5 * 2))
**************************************************************************************************************************/
void CreateInitSmoothGray(const Mat &src, Mat &dst, double sigma = SIGMA)
{
cv::Mat gray; //[1]保存原始圖像的灰度圖像
cv::Mat up; //[2]保存原始圖像的上采樣圖像
//std::cout << "CreateInitSmoothGray sigma = " << sigma << std::endl; //查看sigma初始值
ConvertToGray(src, gray); //灰度圖像
UpSample(gray, up); //[3]圖像的上采樣
//[4]高斯金字塔的-1組的sigma的計算公式
double sigma_init = sqrt(sigma * sigma - (INIT_SIGMA * 2) * (INIT_SIGMA * 2));// INIT_SIGM初始sigma=0.5; -1層的sigma
GaussianSmooth(up, dst, sigma_init); //[5]高斯平滑
}
/*************************************************************************************************************************
*模塊說明：
* 模塊二：3.3 圖像高斯金字塔的構建
在主函數中：int octaves = log((double)min(init_gray.rows, init_gray.cols)) / log(2.0) - 2; //計算高斯金字塔的層數
高斯金字塔中每組圖像中有三層/張圖片 INTERVALS=3；SIGMA=1.6；
**************************************************************************************************************************/
void GaussianPyramid(const Mat &src, vector<Mat>&gauss_pyr, int octaves, int intervals = INTERVALS, double sigma = SIGMA)
{
//std::cout << "sigma is = " << sigma << std::endl;
double *sigmas = new double[intervals + 3];
double k = pow(2.0, 1.0 / intervals); //pow（x,y）,x的y次方；k=2的（1/s)次方
// std::cout << "k is = " << k << std::endl;
sigmas[0] = sigma;
// std::cout << "sigma0 is = " << sigmas[0] << std::endl;
double sig_prev;
double sig_total;
//用循環來表示第0組內各層圖片的尺度，這個可以表示為0組的各層圖像尺度；
for (int i = 1; i < intervals + 3; i++)
{
sig_prev = pow(k, i - 1) * sigma; //第i-1層尺度
sig_total = sig_prev * k; //第i層尺度
sigmas[i] = sqrt(sig_total * sig_total - sig_prev * sig_prev); //組內每層尺度坐標
}
for (int o = 0; o < octaves; o++)
{
//每組多三層
for (int i = 0; i < intervals + 3; i++)
{
Mat mat;
if (o == 0 && i == 0) //第0組，第0層，就是原圖像
{
src.copyTo(mat);
}
else if (i == 0) //新的一組的首層圖像是由上一組最后一層圖像下采樣得到
{
DownSample(gauss_pyr[(o - 1)*(intervals + 3) + intervals], mat);
}
else //對上一層圖像gauss_pyr[i-1]進行參數為sig[i]的高斯平滑，得到當前層圖像
{
GaussianSmooth(gauss_pyr[o * (intervals + 3) + i - 1], mat, sigmas[i]);
// std::cout << "sigmas[i] is = " << sigmas[i] << std::endl;
}
gauss_pyr.push_back(mat); //將mat中的像素添加到gauss_pyr數組里面；也就是說每組從第0層到第5層的圖片像素裝進vector數組中；
}
}
delete[] sigmas;
}
/*************************************************************************************************************************
*模塊說明：
* 圖像的差分
**************************************************************************************************************************/
void Sub(const Mat& a, const Mat& b, Mat & c)
{
if (a.rows != b.rows || a.cols != b.cols || a.type() != b.type())
return;
if (!c.empty())
return;
c.create(a.size(), a.type());
pixel_t* ap = (pixel_t*)a.data;
pixel_t* ap_end = (pixel_t*)a.dataend;
pixel_t* bp = (pixel_t*)b.data;
pixel_t* cp = (pixel_t*)c.data;
int step = a.step / sizeof(pixel_t);
while (ap != ap_end)
{
*cp++ = *ap++ - *bp++;
}
}
/*************************************************************************************************************************
*模塊說明：
* 模塊三：3.4 高斯差分金字塔
*功能說明：
* 1--2002年，Mikolajczyk在詳細的實驗比較中發現尺度歸一化的高斯拉普拉斯函數的極大值和極小值同其他的特征提取函數，例如：
* 梯度，Hessian或者Harris角點特征比較，能夠產生穩定的圖像特征。
* 2--而Lindberg早在1994年發現高斯差分函數（Difference of Gaussian,簡稱DOG算子）與尺度歸一化的拉普拉斯函數非常相似，因此，
* 利用圖像間的差分代替微分。
**************************************************************************************************************************/
void DogPyramid(const vector<Mat>& gauss_pyr, vector<Mat>& dog_pyr, int octaves, int intervals = INTERVALS)
{
for (int o = 0; o < octaves; o++)
{
for (int i = 1; i < intervals + 3; i++)
{
Mat mat;
Sub(gauss_pyr[o*(intervals + 3) + i], gauss_pyr[o*(intervals + 3) + i - 1], mat);
dog_pyr.push_back(mat);
}
}
}
/*************************************************************************************************************************
*模塊說明：
* 模塊四的第一步:3.4-空間極值點的初步檢測(關鍵點的初步探查)
*功能說明:
* 1--在高斯差分函數之后，得到一系列的關鍵點的疑似點，我們需要對這些關鍵點的疑似點初步進行檢測和篩選
* 檢測是以前點為中心，3pixel*3pixel*3pixel的立方體為鄰域，判斷當前點是否為局部最大或最小。i,j,k=-1,0,1表示當前點的鄰域
**************************************************************************************************************************/
bool isExtremum(int x, int y, const vector<Mat>& dog_pyr, int index)
{
pixel_t * data = (pixel_t *)dog_pyr[index].data;
int step = dog_pyr[index].step / sizeof(data[0]);
pixel_t val = *(data + y*step + x);
if (val > 0)
{
for (int i = -1; i <= 1; i++) //層
{
int stp = dog_pyr[index + i].step / sizeof(data[0]);
for (int j = -1; j <= 1; j++) //行
{
for (int k = -1; k <= 1; k++) //列
{
if (val < *((pixel_t*)dog_pyr[index + i].data + stp*(y + j) + (x + k)))
{
return false; //不是極值點退出
}
}
}
}
}
else
{
for (int i = -1; i <= 1; i++)
{
int stp = dog_pyr[index + i].step / sizeof(pixel_t);
for (int j = -1; j <= 1; j++)
{
for (int k = -1; k <= 1; k++)
{
//檢查最小極值
if (val > *((pixel_t*)dog_pyr[index + i].data + stp*(y + j) + (x + k)))
{
return false; //不是極值點退出
}
}
}
}
}
return true; //是級值點，返回
}
/*************************************************************************************************************************
*模塊說明：
* 模塊四的第三步:4.2--消除邊緣響應點
*功能說明:
* 1）一個定義不好的高斯差分算子的極值在橫跨邊緣的地方有較大的住主曲率，在垂直邊緣的方向有較小的主曲率。
* 2）DOG算子會產生較強的邊緣響應，需要剔除不穩定的邊緣響應點，獲取特征點處的Hessian矩陣，主曲率通過一個2*2的Hessian矩
* 陣H求出
* 3）主曲率D和Hessian矩陣的特征值成正比，公式(r+1)*(r+1)/r的值在兩個特征值相等時最小；這個值越大，說明兩個特征值的比值
* 越大，即在某一個方向的梯度值越大，而在另一個方向的梯度值越小，而邊緣恰恰就是這種情況。所以，為了剔除邊緣響應點，
* 需要讓該比值小於一定的閾值，因此，為了檢測主曲率是否在某閾值r下，只需檢測。CSDN論文中的公式(4-7成立)，成立，將關
* 鍵點保留，反之，剔除。
**************************************************************************************************************************/
#define DAt(x, y) (*(data+(y)*step+(x)))
bool passEdgeResponse(int x, int y, const vector<Mat>& dog_pyr, int index, double r = RATIO)
{
pixel_t *data = (pixel_t *)dog_pyr[index].data; //極值點的圖片索引定位
int step = dog_pyr[index].step / sizeof(data[0]);
pixel_t val = *(data + y*step + x); //極值點
double Dxx;
double Dyy;
double Dxy;
double Tr_h; //[1]Hessian矩陣的跡
double Det_h; //[2]Hessian矩陣所對應的行列式的值
Dxx = DAt(x + 1, y) + DAt(x - 1, y) - 2 * val;
Dyy = DAt(x, y + 1) + DAt(x, y - 1) - 2 * val;
Dxy = (DAt(x + 1, y + 1) + DAt(x - 1, y - 1) - DAt(x - 1, y + 1) - DAt(x + 1, y - 1)) / 4.0;
Tr_h = Dxx + Dyy;
Det_h = Dxx * Dyy - Dxy * Dxy;
if (Det_h <= 0)return false;
if (Tr_h * Tr_h / Det_h < (r + 1) * (r + 1) / r) return true;
return false;
}
/*************************************************************************************************************************
*模塊說明：
* 有限差分求導？
**************************************************************************************************************************/
#define Hat(i, j) (*(H+(i)*3 + (j)))
double PyrAt(const vector<Mat>& pyr, int index, int x, int y)
{
pixel_t *data = (pixel_t*)pyr[index].data;
int step = pyr[index].step / sizeof(data[0]);
pixel_t val = *(data + y*step + x);
return val;
}
/*************************************************************************************************************************
*模塊說明：
* 有限差分求導？
**************************************************************************************************************************/
#define At(index, x, y) (PyrAt(dog_pyr, (index), (x), (y)))
//3維D(x)一階偏導,dx列向量
void DerivativeOf3D(int x, int y, const vector<Mat>& dog_pyr, int index, double *dx)
{
double Dx = (At(index, x + 1, y) - At(index, x - 1, y)) / 2.0;
double Dy = (At(index, x, y + 1) - At(index, x, y - 1)) / 2.0;
double Ds = (At(index + 1, x, y) - At(index - 1, x, y)) / 2.0;
dx[0] = Dx;
dx[1] = Dy;
dx[2] = Ds;
}
//3維D(x)二階偏導，即Hessian矩陣
void Hessian3D(int x, int y, const vector<Mat>& dog_pyr, int index, double *H)
{
double val, Dxx, Dyy, Dss, Dxy, Dxs, Dys;
val = At(index, x, y);
Dxx = At(index, x + 1, y) + At(index, x - 1, y) - 2 * val;
Dyy = At(index, x, y + 1) + At(index, x, y - 1) - 2 * val;
Dss = At(index + 1, x, y) + At(index - 1, x, y) - 2 * val;
Dxy = (At(index, x + 1, y + 1) + At(index, x - 1, y - 1)
- At(index, x + 1, y - 1) - At(index, x - 1, y + 1)) / 4.0;
Dxs = (At(index + 1, x + 1, y) + At(index - 1, x - 1, y)
- At(index - 1, x + 1, y) - At(index + 1, x - 1, y)) / 4.0;
Dys = (At(index + 1, x, y + 1) + At(index - 1, x, y - 1)
- At(index + 1, x, y - 1) - At(index - 1, x, y + 1)) / 4.0;
Hat(0, 0) = Dxx;
Hat(1, 1) = Dyy;
Hat(2, 2) = Dss;
Hat(1, 0) = Hat(0, 1) = Dxy;
Hat(2, 0) = Hat(0, 2) = Dxs;
Hat(2, 1) = Hat(1, 2) = Dys;
}
/*************************************************************************************************************************
*模塊說明：
* 4.4 三階矩陣求逆
**************************************************************************************************************************/
#define HIat(i, j) (*(H_inve+(i)*3 + (j)))
//3*3階矩陣求逆
bool Inverse3D(const double *H, double *H_inve)
{
double A = Hat(0, 0)*Hat(1, 1)*Hat(2, 2)
+ Hat(0, 1)*Hat(1, 2)*Hat(2, 0)
+ Hat(0, 2)*Hat(1, 0)*Hat(2, 1)
- Hat(0, 0)*Hat(1, 2)*Hat(2, 1)
- Hat(0, 1)*Hat(1, 0)*Hat(2, 2)
- Hat(0, 2)*Hat(1, 1)*Hat(2, 0);
if (fabs(A) < 1e-10) return false;
HIat(0, 0) = Hat(1, 1) * Hat(2, 2) - Hat(2, 1)*Hat(1, 2);
HIat(0, 1) = -(Hat(0, 1) * Hat(2, 2) - Hat(2, 1) * Hat(0, 2));
HIat(0, 2) = Hat(0, 1) * Hat(1, 2) - Hat(0, 2)*Hat(1, 1);
HIat(1, 0) = Hat(1, 2) * Hat(2, 0) - Hat(2, 2)*Hat(1, 0);
HIat(1, 1) = -(Hat(0, 2) * Hat(2, 0) - Hat(0, 0) * Hat(2, 2));
HIat(1, 2) = Hat(0, 2) * Hat(1, 0) - Hat(0, 0)*Hat(1, 2);
HIat(2, 0) = Hat(1, 0) * Hat(2, 1) - Hat(1, 1)*Hat(2, 0);
HIat(2, 1) = -(Hat(0, 0) * Hat(2, 1) - Hat(0, 1) * Hat(2, 0));
HIat(2, 2) = Hat(0, 0) * Hat(1, 1) - Hat(0, 1)*Hat(1, 0);
for (int i = 0; i < 9; i++)
{
*(H_inve + i) /= A;
}
return true;
}
/*************************************************************************************************************************
*模塊說明：
*
**************************************************************************************************************************/
//計算x^
void GetOffsetX(int x, int y, const vector<Mat>& dog_pyr, int index, double *offset_x)
{
//x^ = -H^(-1) * dx; dx = (Dx, Dy, Ds)^T
double H[9], H_inve[9] = { 0 };
Hessian3D(x, y, dog_pyr, index, H);
Inverse3D(H, H_inve);
double dx[3];
DerivativeOf3D(x, y, dog_pyr, index, dx);
for (int i = 0; i < 3; i++)
{
offset_x[i] = 0.0;
for (int j = 0; j < 3; j++)
{
offset_x[i] += H_inve[i * 3 + j] * dx[j];
}
offset_x[i] = -offset_x[i];
}
}
//計算|D(x^)|
double GetFabsDx(int x, int y, const vector<Mat>& dog_pyr, int index, const double* offset_x)
{
//|D(x^)|=D + 0.5 * dx * offset_x; dx=(Dx, Dy, Ds)^T
double dx[3];
DerivativeOf3D(x, y, dog_pyr, index, dx);
double term = 0.0;
for (int i = 0; i < 3; i++)
term += dx[i] * offset_x[i];
pixel_t *data = (pixel_t *)dog_pyr[index].data;
int step = dog_pyr[index].step / sizeof(data[0]);
pixel_t val = *(data + y*step + x);
return fabs(val + 0.5 * term);
}
/*************************************************************************************************************************
*模塊說明：
* 模塊四的第二步:修正極值點，刪除不穩定的點
*功能說明:
* 1--根據高斯差分函數產生的極值點並不全都是穩定的特征點，因為某些極值點的響應較弱，而且DOG算子會產生較強的邊緣響應
* 2--以上方法檢測到的極值點是離散空間的極值點，下面通過擬合三維二次函數來精確定位關鍵點的位置和尺度，同時去除對比度
* 低和不穩定的邊緣響應點(因為DOG算子會產生較強的邊緣響應)，以增強匹配的穩定性、提高抗噪聲的能力。
* 3--修正極值點，刪除不穩定點，|D(x)| < 0.03 Lowe 2004
**************************************************************************************************************************/
Keypoint* InterploationExtremum(int x, int y, const vector<Mat>& dog_pyr, int index, int octave, int interval, double dxthreshold = DXTHRESHOLD)
{
//計算x=(x,y,sigma)^T
//x^ = -H^(-1) * dx; dx = (Dx, Dy, Ds)^T
double offset_x[3] = { 0 };
const Mat &mat = dog_pyr[index];
int idx = index;
int intvl = interval;
int i = 0;
while (i < MAX_INTERPOLATION_STEPS)
{
GetOffsetX(x, y, dog_pyr, idx, offset_x);
//4. Accurate keypoint localization. Lowe
//如果offset_x 的任一維度大於0.5，it means that the extremum lies closer to a different sample point.
if (fabs(offset_x[0]) < 0.5 && fabs(offset_x[1]) < 0.5 && fabs(offset_x[2]) < 0.5)
break;
//用周圍的點代替
x += cvRound(offset_x[0]);
y += cvRound(offset_x[1]);
interval += cvRound(offset_x[2]);
idx = index - intvl + interval;
//此處保證檢測邊時 x+1,y+1和x-1, y-1有效
if (interval < 1 || interval > INTERVALS || x >= mat.cols - 1 || x < 2 || y >= mat.rows - 1 || y < 2)
{
return NULL;
}
i++;
}
//竄改失敗
if (i >= MAX_INTERPOLATION_STEPS)
return NULL;
//rejecting unstable extrema
//|D(x^)| < 0.03取經驗值
if (GetFabsDx(x, y, dog_pyr, idx, offset_x) < dxthreshold / INTERVALS)
{
return NULL;
}
Keypoint *keypoint = new Keypoint;
keypoint->x = x;
keypoint->y = y;
keypoint->offset_x = offset_x[0];
keypoint->offset_y = offset_x[1];
keypoint->interval = interval;
keypoint->offset_interval = offset_x[2];
keypoint->octave = octave;
keypoint->dx = (x + offset_x[0])*pow(2.0, octave);
keypoint->dy = (y + offset_x[1])*pow(2.0, octave);
return keypoint;
}
/*************************************************************************************************************************
*模塊說明：
* 模塊四：3.5 空間極值點的檢測(關鍵點的初步探查)
*功能說明：
* 1--關鍵點是由DOG空間的局部極值點組成的，關鍵點的初步探查是通過同一組內各DoG相鄰兩層圖像之間的比較完成的。為了尋找DoG
* 函數的極值點，每一個像素點都要和它所有相鄰的點比較，看其是否比它的圖像域和尺度域相鄰的點大還是小。
* 2--當然這樣產生的極值點並不全都是穩定的特征點，因為某些極值點相應較弱，而且DOG算子會產生較強的邊緣響應。
**************************************************************************************************************************/
void DetectionLocalExtrema(const vector<Mat>& dog_pyr, vector<Keypoint>& extrema, int octaves, int intervals = INTERVALS)
{
double thresh = 0.5 * DXTHRESHOLD / intervals; //|D(x^)| < 0.03 ; DXTHRESHOLD=0.03極值點太多
for (int o = 0; o < octaves; o++)
{
//第一層和最后一層極值忽略
for (int i = 1; i < (intervals + 2) - 1; i++) //Dog圖像每一組5層（差分6-1），又第一層和最后一程忽略；
{
int index = o*(intervals + 2) + i; //[1]圖片索引的定位
pixel_t *data = (pixel_t *)dog_pyr[index].data; //[2]獲取圖片的矩陣體的首地址
// std::cout << "首地址為 = " << sizeof(data[0]) << std::endl;
int step = dog_pyr[index].step / sizeof(data[0]); //[3]說明矩陣在存儲空間中的存儲是以線性空間的方式存放的
for (int y = IMG_BORDER; y < dog_pyr[index].rows - IMG_BORDER; y++) //IMG_BORDER=5 ，這里我查找資料說的是指檢測邊界以內的點，因為邊界易受干擾
{
for (int x = IMG_BORDER; x < dog_pyr[index].cols - IMG_BORDER; x++)
{
// std::cout << "x is = " << x << std::endl;
pixel_t val = *(data + y*step + x); //遍歷像素點
if (std::fabs(val) > thresh) //[4]排除閾值過小的點，fabs（）取絕對值
{
if (isExtremum(x, y, dog_pyr, index)) //[5]判斷是否是極值
{
Keypoint *extrmum = InterploationExtremum(x, y, dog_pyr, index, o, i);//極值點x,y；當前組當前層；
if (extrmum)
{
if (passEdgeResponse(extrmum->x, extrmum->y, dog_pyr, index))
{
extrmum->val = *(data + extrmum->y*step + extrmum->x);
extrema.push_back(*extrmum);
}
delete extrmum;
}//extrmum
}//isExtremum
}//std::fabs
}//for x
}//for y
}
}
}
/*************************************************************************************************************************
*模塊說明：
* 模塊五：計算特征點對應的尺度
*
**************************************************************************************************************************/
void CalculateScale(vector<Keypoint>& features, double sigma = SIGMA, int intervals = INTERVALS)
{
double intvl = 0;
for (int i = 0; i < features.size(); i++)
{
intvl = features[i].interval + features[i].offset_interval;
features[i].scale = sigma * pow(2.0, features[i].octave + intvl / intervals);
features[i].octave_scale = sigma * pow(2.0, intvl / intervals);
}
}
//對擴大的圖像特征縮放
/*將特征點序列中每個特征點的坐標減半(圖像金字塔設置了將圖像放大為原圖的2倍時，特征點檢測完之后調用)
/features：特征點序列*/
void HalfFeatures(vector<Keypoint>& features)
{
for (int i = 0; i < features.size(); i++)
{
features[i].dx /= 2;
features[i].dy /= 2;
features[i].scale /= 2;
}
}
/********************************************************************************************************************************
*模塊說明:
* 模塊六---步驟2：計算關鍵點的梯度和梯度方向
*功能說明：
* 1）計算關鍵點(x,y)處的梯度幅值和梯度方向
* 2）將所計算出來的梯度幅值和梯度方向保存在變量mag和ori中
*********************************************************************************************************************************/
bool CalcGradMagOri(const Mat& gauss, int x, int y, double& mag, double& ori)
{
if (x > 0 && x < gauss.cols - 1 && y > 0 && y < gauss.rows - 1)
{
pixel_t *data = (pixel_t*)gauss.data;
int step = gauss.step / sizeof(*data);
double dx = *(data + step*y + (x + 1)) - (*(data + step*y + (x - 1))); //[1]利用X方向上的差分代替微分dx
double dy = *(data + step*(y + 1) + x) - (*(data + step*(y - 1) + x)); //[2]利用Y方向上的差分代替微分dy
mag = sqrt(dx*dx + dy*dy); //[3]計算該關鍵點的梯度幅值
ori = atan2(dy, dx); //[4]計算該關鍵點的梯度方向
return true;
}
else
return false;
}
/********************************************************************************************************************************
*模塊說明:
* 模塊六---步驟1：計算梯度的方向直方圖
*功能說明：
* 1）直方圖以每10度為一個柱，共36個柱，柱代表的方向為為像素點的梯度方向，柱的長短代表了梯度幅值。
* 2）根據Lowe的建議，直方圖統計采用3*1.5*sigma
* 3）在直方圖統計時，每相鄰三個像素點采用高斯加權，根據Lowe的建議，模板采用[0.25,0.5,0.25],並且連續加權兩次
*結 論：
* 圖像的關鍵點檢測完畢后，每個關鍵點就擁有三個信息：位置、尺度、方向；同時也就使關鍵點具備平移、縮放和旋轉不變性
bin=36.
*********************************************************************************************************************************/
double* CalculateOrientationHistogram(const Mat& gauss, int x, int y, int bins, int radius, double sigma)
{
//std::cout << "首地址1為 = " << radius << std::endl; //radius=3*1.5?*?σ(extrema[i].octave_scale)
//std::cout << "首地址2為 = " << bins << std::endl; // sigma=1.5*σ(extrema[i].octave_scale)
double* hist = new double[bins]; //[1]動態分配一個double類型的數組，用來存放梯度直方圖
for (int i = 0; i < bins; i++) //[2]給這個數組初始化
*(hist + i) = 0.0;
double mag; //[3]關鍵點的梯度幅值
double ori; //[4]關鍵點的梯度方向
double weight;
int bin;
const double PI2 = 2.0*CV_PI; //2Π
double econs = -1.0 / (2.0*sigma*sigma);
for (int i = -radius; i <= radius; i++)
{
for (int j = -radius; j <= radius; j++)
{
if (CalcGradMagOri(gauss, x + i, y + j, mag, ori)) //[5]計算該關鍵點的梯度幅值和方向
{
// std::cout << "ori為 = " << ori << std::endl;
weight = exp((i*i + j*j)*econs); // 該點的梯度幅值權重
bin = cvRound(bins * (CV_PI - ori) / PI2); //[6]對一個double行的數進行四舍五入，返回一個整形的數;
// bin = cvRound(bins * (CV_PI + ori) / PI2); //[6]對一個double行的數進行四舍五入，返回一個整形的數;
//計算梯度的方向對應的直方圖中的bin下標,ori范圍(-Π，Π]，+Π后范圍(0，2Π]
// std::cout << "bin1為 = " << bin << std::endl;
bin = bin < bins ? bin : 0; //把bin的下標歸到0~36；bin<bins，bin=bin;否則bin=0；
//std::cout << "bin2為 = " << bin << std::endl;
hist[bin] += mag * weight; //[7]統計梯度的方向直方圖
}
}
}
return hist;
}
/********************************************************************************************************************************
*模塊說明:
* 模塊六---步驟3：對梯度方向直方圖進行連續兩次的高斯平滑
*功能說明：
* 1）在直方圖統計時，每相鄰三個像素點采用高斯加權，根據Lowe的建議，模板采用[0.25,0.5,0.25],並且連續加權兩次
* 2）對直方圖進行兩次平滑
*********************************************************************************************************************************/
void GaussSmoothOriHist(double *hist, int n)
{
double prev = hist[n - 1];
double temp;
double h0 = hist[0];
for (int i = 0; i < n; i++)
{
temp = hist[i];
hist[i] = 0.25 * prev + 0.5 * hist[i] + 0.25 * (i + 1 >= n ? h0 : hist[i + 1]);//對方向直方圖進行高斯平滑
prev = temp;
}
}
/********************************************************************************************************************************
*模塊說明:
* 模塊六---步驟4：計算方向直方圖中的主方向
*********************************************************************************************************************************/
double DominantDirection(double *hist, int n)
{
double maxd = hist[0];
for (int i = 1; i < n; i++)
{
if (hist[i] > maxd) //求取36個柱中的最大峰值
maxd = hist[i];
}
return maxd;
}
//復制關鍵點
void CopyKeypoint(const Keypoint& src, Keypoint& dst)
{
dst.dx = src.dx;
dst.dy = src.dy;
dst.interval = src.interval;
dst.octave = src.octave;
dst.octave_scale = src.octave_scale;
dst.offset_interval = src.offset_interval;
dst.offset_x = src.offset_x;
dst.offset_y = src.offset_y;
dst.ori = src.ori;
dst.scale = src.scale;
dst.val = src.val;
dst.x = src.x;
dst.y = src.y;
}
/********************************************************************************************************************************
*模塊說明:
* 模塊六---步驟5：計算更加精確的關鍵點主方向----拋物插值
*功能說明：
* 1）方向直方圖的峰值則代表了該特征點的方向，以直方圖中的最大值作為該關鍵點的主方向。為了增強匹配的魯棒性，只保留峰值大於主
* 方向峰值80%的方向作為改關鍵點的輔方向。因此，對於同一梯度值得多個峰值的關鍵點位置，在相同位置和尺度將會有多個關鍵點被
* 創建但方向不同。僅有15%的關鍵點被賦予多個方向，但是可以明顯的提高關鍵點的穩定性。
* 2）在實際編程中，就是把該關鍵點復制成多份關鍵點，並將方向值分別賦給這些復制后的關鍵點
* 3）並且，離散的梯度直方圖要進行【插值擬合處理】，來求得更加精確的方向角度值
*********************************************************************************************************************************/
#define Parabola_Interpolate(l, c, r) (0.5*((l)-(r))/((l)-2.0*(c)+(r))) //求小柱子在直方圖中的索引號
void CalcOriFeatures(const Keypoint& keypoint, vector<Keypoint>& features, const double *hist, int n, double mag_thr)
{
double bin;
double PI2 = CV_PI * 2.0;
int l;
int r;
//遍歷直方圖
for (int i = 0; i < n; i++)
{
l = (i == 0) ? n - 1 : i - 1; //前一個(左邊的)bin的下標，
r = (i + 1) % n; //后一個(右邊的)bin的下標
//i=0,l=35,r=1;
//hist[i]是極值
//若當前的bin是局部極值(比前一個和后一個bin都大)，並且值大於給定的幅值閾值，則新生成一個特征點並添加到特征點序列末尾
if (hist[i] > hist[l] && hist[i] > hist[r] && hist[i] >= mag_thr) //mag_thr=0.8*關鍵點最大值
{
//根據左、中、右三個bin的值對當前bin進行直方圖插值
bin = i + Parabola_Interpolate(hist[l], hist[i], hist[r]);
bin = (bin < 0) ? (bin + n) : (bin >= n ? (bin - n) : bin); //將插值結果規范到[0,n]內
Keypoint new_key;
CopyKeypoint(keypoint, new_key); //復制特征點
new_key.ori = ((PI2 * bin) / n) - CV_PI; //方向還原為（-Π，Π]，bin最大值36，2Π*36/36-Π=Π，bin最小值0，2Π*0/36-Π=-Π
features.push_back(new_key); //插入到特征點序列末尾
}
}
}
/********************************************************************************************************************************
*模塊說明:
* 模塊六：5 關鍵點方向分配
*功能說明：
* 1）為了使描述符具有旋轉不變性，需要利用圖像的局部特征為每一個關鍵點分配一個基准方向。使用圖像梯度的方法求取局部結構的穩定
* 方向。
* 2）對於在DOG金字塔中檢測出來的關鍵點，采集其所在高斯金字塔圖像3sigma鄰域窗口內像素的梯度和方向梯度和方向特征。
* 3）梯度的模和方向如下所示:
* 4) 在完成關鍵點的梯度計算后，使用直方圖統計鄰域內像素的梯度和方向。梯度直方圖將0~360度的方向范圍分為36個柱，其中每柱10度，
* 如圖5.1所示，直方圖的峰值方向代表了關鍵點的主方向
cvRound()：返回跟參數最接近的整數值，即四舍五入；
*********************************************************************************************************************************/
void OrientationAssignment(vector<Keypoint>& extrema, vector<Keypoint>& features, const vector<Mat>& gauss_pyr)
{
int n = extrema.size();
double *hist;
for (int i = 0; i < n; i++)
{
//ORI_HIST_BINS=36;ORI_WINDOW_RADIUS=3.0*1.5,特征點方向賦值過程中，搜索鄰域的半徑為：3*1.5*σ
//ORI_SIGMA_TIMES=1.5,
//gauss_pyr[extrema[i].octave*(INTERVALS + 3) + extrema[i].interval是在計算關鍵點在數組中的位置；
hist = CalculateOrientationHistogram(gauss_pyr[extrema[i].octave*(INTERVALS + 3) + extrema[i].interval],
extrema[i].x, extrema[i].y, ORI_HIST_BINS, cvRound(ORI_WINDOW_RADIUS*extrema[i].octave_scale),
ORI_SIGMA_TIMES*extrema[i].octave_scale); //[1]計算梯度的方向直方圖
for (int j = 0; j < ORI_SMOOTH_TIMES; j++)
GaussSmoothOriHist(hist, ORI_HIST_BINS); //[2]對方向直方圖進行高斯平滑
double highest_peak = DominantDirection(hist, ORI_HIST_BINS); //[3]求取方向直方圖中的峰值
//[4]計算更加精確的關鍵點主方向
CalcOriFeatures(extrema[i], features, hist, ORI_HIST_BINS, highest_peak*ORI_PEAK_RATIO); //ORI_PEAK_RATIO=0.8
//精確關鍵點的主方向，拋物插值
delete[] hist;
}
}
/*
功能說明：
將一個條目插入到形成特征描述符的方向直方圖數組中，該條目分配最多8個方向，對於每個維度每個方向乘以（1-d） 的權重，
其中d是以bin測量的bin中心值距離。
//hist為方向直方圖的三維數組
//xbin子行
//ybin子列
//obin子方向
//mag權重
//d方向直方圖的寬度
//bins每個直方圖中bin的個數
三線性插值
*/
void InterpHistEntry(double ***hist, double xbin, double ybin, double obin, double mag, int bins, int d)
{
double d_r, d_c, d_o, v_r, v_c, v_o;
double** row, *h;
int r0, c0, o0, rb, cb, ob, r, c, o;
r0 = cvFloor(ybin); //向下取整
c0 = cvFloor(xbin);
o0 = cvFloor(obin);
d_r = ybin - r0; //小數余項
d_c = xbin - c0;
d_o = obin - o0;
//這里也可以看出前面rbin，cbin為何減0.5，這樣原點上這點d_r,d_c均為0.5，即原點上這點的方向將被平均分配疊加在它周圍4個直方圖上面。
/*
做插值：
xbin,ybin,obin:種子點所在子窗口的位置和方向
所有種子點都將落在4*4的窗口中
r0,c0取不大於xbin，ybin的正整數
r0,c0只能取到0,1,2
xbin,ybin在(-1, 2)
r0取不大於xbin的正整數時。
r0+0 <= xbin <= r0+1
mag在區間[r0,r1]上做插值
obin同理
*/
for (r = 0; r <= 1; r++) // 沿着行方向不超過1個單位長度
{
rb = r0 + r;
if (rb >= 0 && rb < d) //如果此刻還在真正的描述子區間內
{
v_r = mag * ((r == 0) ? 1.0 - d_r : d_r); //對近鄰兩行的貢獻因子
row = hist[rb]; //第rb行上的hist
for (c = 0; c <= 1; c++) //沿着列方向不超過1個單位長度
{
cb = c0 + c;
if (cb >= 0 && cb < d) //如果此刻還在真正的描述子區間內
{
v_c = v_r * ((c == 0) ? 1.0 - d_c : d_c); //對近鄰兩行的貢獻因子
h = row[cb]; //h表示第rb行cb列上的hist
for (o = 0; o <= 1; o++) //沿着直方圖方向不超過1個單位長度
{
ob = (o0 + o) % bins; //bins=8，8個柱子
v_o = v_c * ((o == 0) ? 1.0 - d_o : d_o); //對近鄰兩個方向的貢獻因子
h[ob] += v_o;
}
}
}
}
}
}
/********************************************************************************************************************************
*模塊說明:
* 模塊七--步驟1:計算描述子的直方圖
*功能說明：
gauss：圖像指針
x：特征點所在的行
y：特征點所在的列
ori：特征點的主方向
octave_scale：特征點的尺度//特征點所在組的尺度
width：計算方向直方圖時，將特征點附近划分為width*width個區域，每個區域生成一個直方圖，默認width為4
bins：每個直方圖中bins的個數
返回值：double類型的三維數組，即一個d*d的二維數組，數組中每個元素是一個有n個bins的直方圖數組
*********************************************************************************************************************************/
double*** CalculateDescrHist(const Mat& gauss, int x, int y, double octave_scale, double ori, int bins, int width)
{
double ***hist = new double**[width]; //width*width*bins的三維直方圖數組 //為第一維分配空間
for (int i = 0; i < width; i++)
{
hist[i] = new double*[width]; //為第二維分配空間
for (int j = 0; j < width; j++)
{
hist[i][j] = new double[bins]; //為第三維分配空間
}
}
for (int r = 0; r < width; r++)
for (int c = 0; c < width; c++)
for (int o = 0; o < bins; o++)
hist[r][c][o] = 0.0;
//為了保證特征描述子具有旋轉不變性，要以特征點為中心，在附近鄰域內旋轉θ角，即旋轉為特征點的方向
double cos_ori = cos(ori);
double sin_ori = sin(ori);
//6.1高斯權值，sigma等於描述字窗口寬度的一半
double sigma = 0.5 * width;
double conste = -1.0 / (2 * sigma*sigma);
double PI2 = CV_PI * 2;
//計算特征描述子過程中，特征點周圍的width*width個區域中，每個區域的寬度為3*σ個像素，
double sub_hist_width = DESCR_SCALE_ADJUST * octave_scale; //DESCR_SCALE_ADJUST=3
//【1】計算描述子所需的圖像領域區域的半徑
//考慮到要進行雙線性插值，每個區域的寬度應為:DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )
//在考慮到旋轉因素，每個區域的寬度應為：DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )* sqrt(2)
//所以搜索的半徑是：DESCR_SCALE_ADJUST * octave_scale * ( width + 1.0 )* sqrt(2) / 2
int radius = (sub_hist_width*sqrt(2.0)*(width + 1)) / 2.0 + 0.5; //[1]0.5取四舍五入
double grad_ori;
double grad_mag;
//遍歷每個區域的像素
for (int i = -radius; i <= radius; i++)
{
for (int j = -radius; j <= radius; j++)
{
//坐標旋轉為主方向,以確保旋轉不變性。
double rot_x = (cos_ori * j - sin_ori * i) / sub_hist_width;
double rot_y = (sin_ori * j + cos_ori * i) / sub_hist_width;
//將鄰域內的采樣點分配到對應的子區域內，將子區域內的梯度值分配到8個方向上，計算其權值
double xbin = rot_x + width / 2 - 0.5; //[2]xbin,ybin為落在4*4窗口中的下標值
double ybin = rot_y + width / 2 - 0.5;
if (xbin > -1.0 && xbin < width && ybin > -1.0 && ybin < width)
{
if (CalcGradMagOri(gauss, x + j, y + i, grad_mag, grad_ori)) //[3]計算關鍵點的梯度方向
{
grad_ori = (CV_PI - grad_ori) - ori;
while (grad_ori < 0.0)
grad_ori += PI2;
while (grad_ori >= PI2)
grad_ori -= PI2; //把梯度方向歸入[0，2Π]
double obin = grad_ori * (bins / PI2);
//梯度方向*(8/2Π)歸到[0，8]，8為直方圖有8個bin，2Π為角度取值范圍的長度，把方向的索引歸於0~8.
//conste = -1.0 / (2 * sigma*sigma);
//計算權重
double weight = exp(conste*(rot_x*rot_x + rot_y * rot_y));
InterpHistEntry(hist, xbin, ybin, obin, grad_mag*weight, bins, width); //線性差值求每個種子點八個方向的梯度。
}
}
}
}
return hist;
}
//歸一化特征點的特征描述子，即將特征描述子數組中每個元素除以特征描述子的模
void NormalizeDescr(Keypoint& feat)
{
double cur, len_inv, len_sq = 0.0;
int i, d = feat.descr_length; //特征描述子的維數
//求特征描述子的模
for (i = 0; i < d; i++)
{
cur = feat.descriptor[i];
len_sq += cur*cur;
}
len_inv = 1.0 / sqrt(len_sq); //sqrt( len_sq )為特征描述子的模
for (i = 0; i < d; i++)
feat.descriptor[i] *= len_inv; //特征描述子中每個元素除以特征描述子的模，完成歸一化
}
/********************************************************************************************************************************
*模塊說明:
* 模塊七--步驟2:直方圖到描述子的轉換
*功能說明：
將某特征點的方向直方圖轉換為特征描述子向量，對特征描述子歸一化並將所有元素轉化為整型，存入指定特征點中
參數：
hist：width*width*bins的三維直方圖數組
width：計算方向直方圖時，將特征點附近划分為width*width個區域，每個區域生成一個直方圖
bins：每個直方圖的bin個數
feature：特征點指針，將計算好的特征描述子存入其中
*********************************************************************************************************************************/
void HistToDescriptor(double ***hist, int width, int bins, Keypoint& feature)
{
int int_val, i, r, c, o, k = 0;
//遍歷width*width*bins的三維直方圖數組，將其中的所有數據(一般是128個)都存入feature結構的descriptor成員中
for (r = 0; r < width; r++)
for (c = 0; c < width; c++)
for (o = 0; o < bins; o++)
{
feature.descriptor[k++] = hist[r][c][o];
}
feature.descr_length = k; //特征描述子的維數，一般是128
NormalizeDescr(feature); //[1]描述子特征向量歸一化
//歸一化特征點的特征描述子，即將特征描述子數組中每個元素除以特征描述子的模
// DESCR_MAG_THR 0.2
//遍歷特征描述子向量，將超過閾值DESCR_MAG_THR的元素強行賦值DESCR_MAG_THR
for (i = 0; i < k; i++) //[2]描述子向量門限
if (feature.descriptor[i] > DESCR_MAG_THR)
feature.descriptor[i] = DESCR_MAG_THR;
NormalizeDescr(feature); //[3]描述子進行最后一次的歸一化操作
//遍歷特征描述子向量，每個元素乘以系數INT_DESCR_FCTR( 512.0)來變為整型，並且最大值不能超過255
for (i = 0; i < k; i++) //[4]將單精度浮點型的描述子轉換為整形的描述子
{
int_val = INT_DESCR_FCTR * feature.descriptor[i];
feature.descriptor[i] = min(255, int_val);
}
}
/********************************************************************************************************************************
*模塊說明:
* 模塊七:6 關鍵點描述
*功能說明：
* 1）通過以上步驟，對於一個關鍵點，擁有三個信息：位置、尺度、方向
* 2）接下來就是為每個關鍵點建立一個描述符，用一組向量來將這個關鍵點描述出來，使其不隨各種變化而變化，比如光照、視角變化等等
* 3）這個描述子不但包括關鍵點，也包含關鍵點周圍對其貢獻的像素點，並且描述符應該有較高的獨特性，以便於特征點正確的匹配概率
* 1）SIFT描述子----是關鍵點鄰域高斯圖像梯度統計結果的一種表示。
* 2）通過對關鍵點周圍圖像區域分塊，計算塊內梯度直方圖，生成具有獨特性的向量
* 3）這個向量是該區域圖像信息的一種表述和抽象，具有唯一性。
*Lowe論文：
* Lowe建議描述子使用在關鍵點尺度空間內4*4的窗口中計算的8個方向的梯度信息，共4*4*8=128維向量來表征。具體的步驟如下所示:
* 1)確定計算描述子所需的圖像區域
* 2）將坐標軸旋轉為關鍵點的方向，以確保旋轉不變性，如CSDN博文中的圖6.2所示；旋轉后鄰域采樣點的新坐標可以通過公式(6-2)計算
* 3）將鄰域內的采樣點分配到對應的子區域，將子區域內的梯度值分配到8個方向上，計算其權值
* 4）插值計算每個種子點八個方向的梯度
* 5）如上統計的4*4*8=128個梯度信息即為該關鍵點的特征向量。特征向量形成后，為了去除光照變化的影響，需要對它們進行歸一化處理，
* 對於圖像灰度值整體漂移，圖像各點的梯度是鄰域像素相減得到的，所以也能去除。得到的描述子向量為H，歸一化之后的向量為L
* 6）描述子向量門限。非線性光照，相機飽和度變化對造成某些方向的梯度值過大，而對方向的影響微弱。因此，設置門限值（向量歸一化
* 后，一般取0.2）截斷較大的梯度值。然后，在進行一次歸一化處理，提高特征的鑒別性。
* 7）按特征點的尺度對特征描述向量進行排序
* 8）至此，SIFT特征描述向量生成。
將關鍵點附近的區域划分為d*d(Lowe建議d=4)個子區域，每個子區域作為一個種子點，每個種子點有8個方向。
*********************************************************************************************************************************/
void DescriptorRepresentation(vector<Keypoint>& features, const vector<Mat>& gauss_pyr, int bins, int width) //bins=8,width=4
{
double ***hist;
for (int i = 0; i < features.size(); i++)
{ //[1]計算描述子的直方圖
//gauss_pyr[features[i].octave*(INTERVALS + 3) + features[i].interval 計算特征點所在矩陣中的位置；
hist = CalculateDescrHist(gauss_pyr[features[i].octave*(INTERVALS + 3) + features[i].interval],
features[i].x, features[i].y, features[i].octave_scale, features[i].ori, bins, width);
HistToDescriptor(hist, width, bins, features[i]); //[2]直方圖到描述子的轉換
/*釋放計算特征描述子過程中用到的方向直方圖的內存空間
hist：方向直方圖的指針，是一個width*width*bins的三維直方圖數組
*/
for (int j = 0; j < width; j++)
{
for (int k = 0; k < width; k++)
{
delete[] hist[j][k]; //釋放第三維的內存
}
delete[] hist[j]; //釋放第二維的內存
}
delete[] hist; //釋放第一維的內存
}
}
/*比較函數，將特征點按尺度的降序排列
參數：
f1：第一個特征點的指針
f2：第二個特征點的指針
返回值：如果f1的尺度小於f2的尺度，返回1；否則返回-1；若相等返回0
*/
bool FeatureCmp(Keypoint& f1, Keypoint& f2)
{
return f1.scale < f2.scale;
}
/*******************************************************************************************************************
*函數說明:
* 最大的模塊1：SIFT算法模塊
*函數參數說明:
* 1---const Mat &src---------------准備進行特征點檢測的原始圖片
* 2---Vector<Keypoint>& features---用來存儲檢測出來的關鍵點
* 3---double sigma-----------------sigma值
* 4---int intervals----------------關鍵點所在的層數
********************************************************************************************************************/
void Sift(const Mat &src, vector<Keypoint>& features, double sigma, int intervals)
{
std::cout << "【Step_one】Create -1 octave gaussian pyramid image" << std::endl;
std::cout << "sigma = " << sigma << std::endl;
std::cout << "intervals = " << intervals << std::endl;
cv::Mat init_gray;
CreateInitSmoothGray(src, init_gray, sigma); //灰度圖
int octaves = log((double)min(init_gray.rows, init_gray.cols)) / log(2.0) - 2; //計算高斯金字塔的層數
std::cout << "【1】The height and width of init_gray_img = " << init_gray.rows << "*" << init_gray.cols << std::endl;
std::cout << "【2】The octaves of the gauss pyramid = " << octaves << std::endl;
std::cout << "【Step_two】Building gaussian pyramid ..." << std::endl;
std::vector<Mat> gauss_pyr;
GaussianPyramid(init_gray, gauss_pyr, octaves, intervals, sigma); //高斯金字塔
write_pyr(gauss_pyr, "gausspyramid");
std::cout << "【Step_three】Building difference of gaussian pyramid..." << std::endl;
vector<Mat> dog_pyr;
DogPyramid(gauss_pyr, dog_pyr, octaves, intervals); //差分金字塔
write_pyr(dog_pyr, "dogpyramid");
std::cout << "【Step_four】Deatecting local extrema..." << std::endl;
vector<Keypoint> extrema;
DetectionLocalExtrema(dog_pyr, extrema, octaves, intervals); //極值點初探（排除較小點，找到極值點，消除邊緣響應）
std::cout << "【3】keypoints cout: " << extrema.size() << " --" << std::endl;
std::cout << "【4】extrema detection finished." << std::endl;
std::cout << "【5】please look dir gausspyramid, dogpyramid and extrema.txt.--" << endl;
std::cout << "【Step_five】CalculateScale..." << std::endl;
CalculateScale(extrema, sigma, intervals); //計算特征點的尺度
HalfFeatures(extrema); //圖像縮放
std::cout << "【Step_six】Orientation assignment..." << std::endl;
OrientationAssignment(extrema, features, gauss_pyr); //關鍵點直方圖計算，使用的extrema是縮放后的圖像
std::cout << "【6】features count: " << features.size() << std::endl;
std::cout << "【Step_seven】DescriptorRepresentation..." << std::endl;
//DESCR_HIST_BINS=8;DESCR_WINDOW_WIDTH=4;
DescriptorRepresentation(features, gauss_pyr, DESCR_HIST_BINS, DESCR_WINDOW_WIDTH); //關鍵點主方向的描述
sort(features.begin(), features.end(), FeatureCmp); //排序， FeatureCmp是排序方法
cout << "finished." << endl;
}
/*******************************************************************************************************************
*函數說明:
* 畫出SIFT特征點的具體函數
********************************************************************************************************************/
void DrawSiftFeature(Mat& src, Keypoint& feat, CvScalar color)
{
int len, hlen, blen, start_x, start_y, end_x, end_y, h1_x, h1_y, h2_x, h2_y;
double scl, ori;
double scale = 5.0;
double hscale = 0.75;
CvPoint start, end, h1, h2;
/* compute points for an arrow scaled and rotated by feat's scl and ori */
start_x = cvRound(feat.dx);
start_y = cvRound(feat.dy);
scl = feat.scale;
ori = feat.ori;
len = cvRound(scl * scale);
hlen = cvRound(scl * hscale);
blen = len - hlen;
end_x = cvRound(len * cos(ori)) + start_x;
end_y = cvRound(len * -sin(ori)) + start_y;
h1_x = cvRound(blen * cos(ori + CV_PI / 18.0)) + start_x;
h1_y = cvRound(blen * -sin(ori + CV_PI / 18.0)) + start_y;
h2_x = cvRound(blen * cos(ori - CV_PI / 18.0)) + start_x;
h2_y = cvRound(blen * -sin(ori - CV_PI / 18.0)) + start_y;
start = cvPoint(start_x, start_y);
end = cvPoint(end_x, end_y);
h1 = cvPoint(h1_x, h1_y);
h2 = cvPoint(h2_x, h2_y);
line(src, start, end, color, 1, 8, 0);
line(src, end, h1, color, 1, 8, 0);
line(src, end, h2, color, 1, 8, 0);
}
/*******************************************************************************************************************
*函數說明:
* 最大的模塊3：畫出SIFT特征點
********************************************************************************************************************/
void DrawSiftFeatures(Mat& src, vector<Keypoint>& features)
{
CvScalar color = CV_RGB(0, 255, 0);
for (int i = 0; i < features.size(); i++)
{
DrawSiftFeature(src, features[i], color);
}
}
/*******************************************************************************************************************
*函數說明:
* 最大的模塊2：畫出關鍵點KeyPoints
********************************************************************************************************************/
void DrawKeyPoints(Mat &src, vector<Keypoint>& keypoints)
{
int j = 0;
for (int i = 0; i < keypoints.size(); i++)
{
CvScalar color = { 255, 0 ,0 };
circle(src, Point(keypoints[i].dx, keypoints[i].dy), 3, color);
j++;
}
}
const char* GetFileName(const char* dir, int i)
{
char *name = new char[50];
sprintf(name, "%s\\%d\.jpg", dir, i);
return name;
}
void cv64f_to_cv8U(const Mat& src, Mat& dst)
{
double* data = (double*)src.data;
int step = src.step / sizeof(*data);
if (!dst.empty())
return;
dst.create(src.size(), CV_8U);
uchar* dst_data = dst.data;
for (int i = 0, m = 0; i < src.cols; i++, m++)
{
for (int j = 0, n = 0; j < src.rows; j++, n++)
{
*(dst_data + dst.step*j + i) = (uchar)(*(data + step*j + i) * 255);
}
}
}
//通過轉換后保存的圖像，會失真,和imshow顯示出的圖像相差很大
void writecv64f(const char* filename, const Mat& mat)
{
Mat dst;
cv64f_to_cv8U(mat, dst);
imwrite(filename, dst);
}
void write_pyr(const vector<Mat>& pyr, const char* dir)
{
for (int i = 0; i < pyr.size(); i++)
{
writecv64f(GetFileName(dir, i), pyr[i]);
}
}
void read_features(vector<Keypoint>&features, const char*file)
{
ifstream in(file);
int n = 0, dims = 0;
in >> n >> dims;
cout << n << " " << dims << endl;
for (int i = 0; i < n; i++)
{
Keypoint key;
in >> key.dy >> key.dx >> key.scale >> key.ori;
for (int j = 0; j < dims; j++)
{
in >> key.descriptor[j];
}
features.push_back(key);
}
in.close();
}
/*******************************************************************************************************************
*函數說明:
* 最大的模塊4：將特征點寫入文本文件
********************************************************************************************************************/
void write_features(const vector<Keypoint>&features, const char*file)
{
ofstream dout(file);
dout << features.size() << " " << FEATURE_ELEMENT_LENGTH << endl;
for (int i = 0; i < features.size(); i++)
{
dout << features[i].dy << " " << features[i].dx << " " << features[i].scale << " " << features[i].ori << endl;
for (int j = 0; j < FEATURE_ELEMENT_LENGTH; j++)
{
if (j % 20 == 0)
dout << endl;
dout << features[i].descriptor[j] << " ";
}
dout << endl;
}
dout.close();
}
sift.h
#ifndef SIFT_H
#define SIFT_H
#include <opencv2/core/core.hpp>
#include "opencv2/highgui/highgui.hpp"
using namespace cv;
using namespace std;
typedef double pixel_t; //【1】像素類型
#define INIT_SIGMA 0.5 //【2】初始sigma
#define SIGMA 1.6
#define INTERVALS 3 //【3】高斯金字塔中每組圖像中有三層/張圖片
#define RATIO 10 //【4】半徑r
#define MAX_INTERPOLATION_STEPS 5 //【5】最大空間間隔
#define DXTHRESHOLD 0.03 //【6】|D(x^)| < 0.03 0.03極值點太多
#define ORI_HIST_BINS 36 //【7】bings=36
#define ORI_SIGMA_TIMES 1.5
#define ORI_WINDOW_RADIUS 3.0 * ORI_SIGMA_TIMES
#define ORI_SMOOTH_TIMES 2
#define ORI_PEAK_RATIO 0.8
#define FEATURE_ELEMENT_LENGTH 128
#define DESCR_HIST_BINS 8
#define IMG_BORDER 5
#define DESCR_WINDOW_WIDTH 4
#define DESCR_SCALE_ADJUST 3
#define DESCR_MAG_THR 0.2
#define INT_DESCR_FCTR 512.0
/*********************************************************************************************
*模塊說明：
* 關鍵點/特征點的結構體聲明
*注意點1：
* 在高斯金字塔構建的過程中，一幅圖像可以產生好幾組(octave)圖像，一組圖像包含幾層(inteval)
* 圖像
*********************************************************************************************/
struct Keypoint
{
int octave; //【1】關鍵點所在組
int interval; //【2】關鍵點所在層
double offset_interval; //【3】調整后的層的增量
int x; //【4】x,y坐標,根據octave和interval可取的層內圖像
int y;
double scale; //【5】空間尺度坐標scale = sigma0*pow(2.0, o+s/S)
double dx; //【6】特征點坐標，該坐標被縮放成原圖像大小
double dy;
double offset_x;
double offset_y;
//============================================================
//1---高斯金字塔組內各層尺度坐標，不同組的相同層的sigma值相同
//2---關鍵點所在組的組內尺度
//============================================================
double octave_scale; //【1】offset_i;
double ori; //【2】方向
int descr_length;
double descriptor[FEATURE_ELEMENT_LENGTH]; //【3】特征點描述符
double val; //【4】極值
};
/*********************************************************************************
*模塊說明：
* SIFT算法內,所有成員函數的聲明
*********************************************************************************/
void read_features(vector<Keypoint>&features, const char*file);
void write_features(const vector<Keypoint>&features, const char*file);
void testInverse3D();
void write_pyr(const vector<Mat>& pyr, const char* dir);
void DrawKeyPoints(Mat &src, vector<Keypoint>& keypoints);
const char* GetFileName(const char* dir, int i);
void ConvertToGray(const Mat& src, Mat& dst);
void DownSample(const Mat& src, Mat& dst);
void UpSample(const Mat& src, Mat& dst);
void GaussianTemplateSmooth(const Mat &src, Mat &dst, double);
void GaussianSmooth2D(const Mat &src, Mat &dst, double sigma);
void GaussianSmooth(const Mat &src, Mat &dst, double sigma);
void Sift(const Mat &src, vector<Keypoint>& features, double sigma = SIGMA, int intervals = INTERVALS);
void CreateInitSmoothGray(const Mat &src, Mat &dst, double);
void GaussianPyramid(const Mat &src, vector<Mat>&gauss_pyr, int octaves, int intervals, double sigma);
void Sub(const Mat& a, const Mat& b, Mat & c);
void DogPyramid(const vector<Mat>& gauss_pyr, vector<Mat>& dog_pyr, int octaves, int intervals);
void DetectionLocalExtrema(const vector<Mat>& dog_pyr, vector<Keypoint>& extrema, int octaves, int intervals);
void DrawSiftFeatures(Mat& src, vector<Keypoint>& features);
#endif
參考：
鏈接：https://blog.csdn.net/jkldlnx/article/details/72388669
OpenCV成長之路(9)：特征點檢測與圖像匹配
最后c代碼解釋：
https://wenku.baidu.com/view/d7edd2464b73f242336c5ffa.html
https://www.cnblogs.com/Alliswell-WP/tag/SIFT/