平臺：win10 x64 +VS 2015專業版 +opencv-2.4.11 + gtk_-bundle_2.24.10_win32html
主要參考：1.代碼：RobHess的SIFT源碼c++
2.書：王永明 王貴錦 《圖像局部不變性特徵與描述》算法
SIFT四步驟和特徵匹配及篩選：數組
步驟一：創建尺度空間，即創建高斯差分(DoG)金字塔dog_pyr數據結構
步驟二：在尺度空間中檢測極值點，並進行精肯定位和篩選建立默認大小的內存存儲器ide
步驟三：特徵點方向賦值，完成此步驟後，每一個特徵點有三個信息：位置、尺度、方向函數
步驟二：在尺度空間中檢測極值點，並進行精肯定位和篩選建立默認大小的內存存儲器
問題及解答：
（1）問題描述：如何找尋關鍵點？在幾層之間對比？思路是什麼？
答：極值的檢測是在DoG空間進行的，檢測是之前點爲中心，3pixel*3pixel*3pixel的立方體爲鄰域，判斷當前點是否爲局部最大或最小。以下圖所示，橘黃色爲當前檢測點，綠色點爲其鄰域。由於要比較當前點的上下層圖像，因此極值檢測從DoG每層的第2幅圖像開始，終止於每層的倒數第2幅圖像（第1幅沒有下層，最後1幅沒有上層，沒法比較）。
參考書P84-P85及圖4-6
（2）問題描述：爲何還要對找尋的關鍵點進行插值？爲何要精肯定位？
答：以上極值點的搜索時在離散空間中進行的，檢測到的極值點並非真正意義上的極值點。以下圖所示，連續空間中極值與離散空間的區別。一般經過插值的方式，利用離散的值來插值，求取接近真正的極值的點。
對於一維函數，利用泰勒級數，將其展開爲二次函數：
f(x) ≈ f(0) + f'(0)x + f''(0)x2
對於二維函數，泰勒展開爲：
矩陣表示爲：
矩陣表示爲：
矢量表示爲：
當矢量爲n維時，有：
書中P86頁有錯誤，我已經更正，並附在：https://www.cnblogs.com/Alliswell-WP/p/SIFT.html
求取f(x)的極值，只需求取∂f/∂x = 0。對於極值，x,y,σ三個變量，即爲三維空間。利用三維子像元插值，設其函數爲D(x, y, σ)，令x = (x, y, σ)T，那麼在第一節中找到的極值點進行泰勒展開爲（式1）以下：
（式1）
其中D爲極值點的值，∂DT/∂x爲在極值點各自變量的倒數，∂2D/∂x2爲其在展開點相應的矩陣。對上式求導，另∂D(x)/∂x = 0，結果以下式，對應的（爲了表示方便，ˆx代替）向量即爲真正極值點偏離插值點的量。求解得（式2）以下：
（式2）
最終極值點的位置即爲插值點x+ˆx，且屢次迭代能夠提升精度（通常爲5次迭代）。
參考書P85-P87及圖4-7和圖4-8
（3）問題描述：既然知道了檢測到的極值點並非真正意義上的極值點，那怎麼怎麼插值呢？
答：圖像離散的像素點組成的，因此要差分代替偏導，計算一階偏導，二階偏導構建Hessian矩陣。
（4）問題描述：對關鍵點定位後，須要剔除一些很差的KeyPoint，那什麼是很差的KeyPoint的呢？如何剔除？
答：
1.DoG響應較低的點，即極值較小的點。
2.響應較強的點也不是穩定的特徵點。DoG對圖像中的邊緣有較強的響應值，因此落在圖像邊緣的點也不是穩定的特徵點。一方面圖像邊緣上的點是很難定位的，具備定位的歧義性；另外一方面這樣的點很容易受到噪聲的干擾變得不穩定。
對於第一種，只需計算矯正後的點的響應值D(ˆx)，響應值小於必定閾值，即認爲該點效應較小，將其剔除。將（式2）帶入（式1），求解得：
在Lowe文章中，將|D(ˆx)|<0.03（圖像灰度歸一化爲[0,1]）的特徵點剔除。
對於第二種，利用Hessian矩陣來剔除。一個平坦的DoG響應峯值在橫跨邊緣的地方有較大的主曲率，而在垂直邊緣的地方有較小的主曲率。主曲率能夠經過2×2的Hessian矩陣H求出：
D值能夠經過求臨近點差分獲得。H的特徵值與D的主曲率成正比，具體可參見Harris角點檢測算法。爲了不求具體的值，咱們能夠經過H將特徵值的比例表示出來。令爲最大特徵值，爲最小特徵值，那麼：
Tr(H)表示矩陣H的跡，Det(H)表示H的行列式。令表示最大特徵值與最小特徵值的比值，則有：
上式與兩個特徵值的比例有關。隨着主曲率比值的增長，也會增長。咱們只須要去掉比率大於必定值的特徵點。Lowe論文中去掉r=10的點。
參考書P87-P88
（5）問題描述：RobHess的SIFT源碼如何實現步驟二，大概思路是這樣的？
答：（5.1）代碼及說明：
/*步驟二：在尺度空間中檢測極值點，並進行精肯定位和篩選建立默認大小的內存存儲器*/
storage = cvCreateMemStorage( 0 ); //調用opencv的cvCreateMemStorage函數，函數功能：用來建立一個內存存儲器，來統一管理各類動態對象的內存。函數返回一個新建立的內存存儲器指針。
//在尺度空間中檢測極值點，經過插值精肯定位，去除低對比度的點，去除邊緣點，
//返回檢測到的特徵點序列
features = scale_space_extrema( dog_pyr, octvs, intvls, contr_thr, curv_thr, storage );
//計算特徵點序列features中每一個特徵點的尺度
calc_feature_scales( features, sigma, intvls );
if( img_dbl ) //若設置了將圖像放大爲原圖的2倍
adjust_for_img_dbl( features );//將特徵點序列中每一個特徵點的座標減半
//(當設置了將圖像放大爲原圖的2倍時，特徵點檢測完以後調用)
（5.2）scale_space_extrema代碼及說明：
/*在尺度空間中檢測極值點，經過插值（三維二次函數）精肯定位，去除低對比度的點，去除邊緣點，返回檢測到的特徵點序列
參數：
dog_pyr：高斯差分金字塔
octvs：高斯差分金字塔的組數
intvls：每組的層數
contr_thr：對比度閾值，針對歸一化後的圖像，用來去除不穩定特徵
cur_thr：主曲率比值的閾值，用來去除邊緣特徵
storage：存儲器
返回值：返回檢測到的特徵點的序列*/
static CvSeq* scale_space_extrema( IplImage*** dog_pyr, int octvs, int intvls,
double contr_thr, int curv_thr, CvMemStorage* storage ) //octvs=O(log(min(length,width))/log(2)-2)；intvls=3；curv_thr=0.04；curv_thr=10
{
CvSeq* features;//特徵點序列
double prelim_contr_thr = 0.5 * contr_thr / intvls;//像素的對比度閾值，此處與LOWE論文不一樣，，|D(x)|<0.03，書P87
struct feature* feat;
struct detection_data* ddata;
int o, i, r, c;
//在存儲器storage上建立存儲極值點的序列，其中存儲feature結構類型的數據
features = cvCreateSeq( 0, sizeof(CvSeq), sizeof(struct feature), storage );
/*遍歷高斯差分金字塔，檢測極值點*/
//SIFT_IMG_BORDER指明邊界寬度，只檢測邊界線之內的極值點
for( o = 0; o < octvs; o++ )//第o組
for( i = 1; i <= intvls; i++ )//遍i層
for(r = SIFT_IMG_BORDER; r < dog_pyr[o][0]->height-SIFT_IMG_BORDER; r++)//第r行
for(c = SIFT_IMG_BORDER; c < dog_pyr[o][0]->width-SIFT_IMG_BORDER; c++)//第c列
//進行初步的對比度檢查，只有當歸一化後的像素值大於對比度
//閾值prelim_contr_thr時才繼續檢測此像素點是否多是極值
//調用函數pixval32f獲取圖像dog_pyr[o][i]的第r行第c列的點的座標值，
//而後調用ABS宏求其絕對值
if( ABS( pixval32f( dog_pyr[o][i], r, c ) ) > prelim_contr_thr )
//經過在尺度空間中將一個像素點的值與其周圍3*3*3鄰域內的點比較來
//決定此點是否極值點(極大值或極小都行)
if( is_extremum( dog_pyr, o, i, r, c ) )//如果極值點
{
//因爲極值點的檢測是在離散空間中進行的，因此檢測到的極值點並
//不必定是真正意義上的極值點
//由於真正的極值點可能位於兩個像素之間，而在離散空間中只能精
//確到座標點精度上
//經過亞像素級插值進行極值點精肯定位(修正極值點座標)，並去除
//低對比度的極值點，將修正後的特徵點組成feature結構返回
feat = interp_extremum(dog_pyr, o, i, r, c, intvls, contr_thr);
//返回值非空，代表此點已被成功修正
if( feat )
{
//調用宏feat_detection_data來提取參數feat中的feature_data成員
//並轉換爲detection_data類型的指針
ddata = feat_detection_data( feat );
//去除邊緣響應，即經過計算主曲率比值判斷某點是否邊緣點,
//返回值爲0表示不是邊緣點，可作特徵點
if( ! is_too_edge_like( dog_pyr[ddata->octv][ddata->intvl],
ddata->r, ddata->c, curv_thr ) )
{
cvSeqPush( features, feat );//向特徵點序列features末尾插
//入新檢測到的特徵點feat
}
else
free( ddata );
free( feat );
}
}
return features;//返回特徵點序列
}
（5.2.1）CvSeq的理解：
說明：動態結構序列CvSeq是全部OpenCv動態數據結構的基礎。有兩種類型的序列：稠密序列，稀疏序列：
其中一種是——稠密序列都派生自CvSeq，他們用來表明可擴展的一維數組 — 向量、棧、隊列和雙端隊列。數據間不存在空隙（連續存儲）。若是元素元素從序列中間被刪除或插入新的元素到序列，那麼此元素後邊的相關元素所有被移動。
total表示稠密序列的元素個數，或者稀疏序列被分配的節點數。elem_size表示序列中每一個元素佔用的字節數。block_max是最近一個內存的最大邊界指針。ptr表示當寫指針。delta_elems表示序列間隔尺寸。storage指向序列存儲的內存塊的指針。free_blocks表示空的塊列表。first指向第一個序列塊。
參看博客：1）OpenCV——CvSeq動態結構序列：https://www.cnblogs.com/farewell-farewell/p/5999908.html
2）CvSeq的理解：https://blog.csdn.net/weijianmeng/article/details/7173560
（5.2.2）cvCreateSeq的理解：
說明：CvSeq *cvCreateSeq(int seq_flags, int header_size, int elem_size, CvMemStorage *storage);
header_size：序列頭的大小，一般爲sizeof(CvSeq)。
elem_size：存儲元素的大小。
storage：內存存儲器，添加元素的時候，就會從內存存儲器申請空間。
參看博客：opencv建立序列cvCreateSeq與插入元素cvSeqPush的運用：https://www.cnblogs.com/farewell-farewell/p/5999908.html
（5.2.3）for(r = SIFT_IMG_BORDER; r < dog_pyr[o][0]->height-SIFT_IMG_BORDER; r++)//第r行
for(c = SIFT_IMG_BORDER; c < dog_pyr[o][0]->width-SIFT_IMG_BORDER; c++)//第c列
中爲何須要設置SIFT_IMG_BORDER指明邊界寬度：
答：只檢測邊界線之內的極值點，由於邊界容易受到噪聲的干擾而不穩定。
（5.2.4）pixval32f(dog_pyr[o][i])函數功能？ABS(pixval32f(dog_pyr[o][i]))爲何加絕對值？爲何須要判斷高斯差分金字塔（DOG）中某個像素小於0的狀況？
答：utils.c中定義了這個內部函數：
static inline float pixval32f( IplImage* img, int r, int c )
{
return ( (float*)(img->imageData + img->widthStep*r) )[c];
}
具體怎麼訪問圖像元素？
參看博客：1） opencv——訪問圖像元素(imagedata widthstep)：https://blog.csdn.net/xiaofeilong321/article/details/13287697
2）opencv學習筆記（八）IplImage* 訪問圖像像素的值：https://www.cnblogs.com/codingmengmeng/p/6559724.html
（5.2.5）is_extremum代碼及說明：
/*經過在尺度空間中將一個像素點的值與其周圍3*3*3鄰域內的點比較來決定此點是否極值點
(極大值或極小都行)
參數：
dog_pyr：高斯差分金字塔
octv：像素點所在的組
intvl：像素點所在的層
r：像素點所在的行
c：像素點所在的列
返回值：若指定的像素點是極值點(極大值或極小值)，返回1；不然返回0*/
static int is_extremum( IplImage*** dog_pyr, int octv, int intvl, int r, int c )
{
//調用函數pixval32f獲取圖像dog_pyr[octv][intvl]的第r行第c列的點的座標值
float val = pixval32f( dog_pyr[octv][intvl], r, c );
int i, j, k;
//檢查是否最大值
if( val > 0 )
{
for( i = -1; i <= 1; i++ )//層
for( j = -1; j <= 1; j++ )//行
for( k = -1; k <= 1; k++ )//列
if( val < pixval32f( dog_pyr[octv][intvl+i], r + j, c + k ) )
return 0; //不是極值點，退出此函數
}
//檢查是否最小值
else
{
for( i = -1; i <= 1; i++ )//層
for( j = -1; j <= 1; j++ )//行
for( k = -1; k <= 1; k++ )//列
if( val > pixval32f( dog_pyr[octv][intvl+i], r + j, c + k ) )
return 0;
}
return 1; //是極值點，返回1
}
（5.2.6）interp_extremum代碼及說明：
剔除不穩定點，精肯定位關鍵點位置
/*經過亞像素級插值進行極值點精肯定位(修正極值點座標)，並去除低對比度的極值點，
將修正後的特徵點組成feature結構返回
參數：
dog_pyr：高斯差分金字塔
octv：像素點所在的組
intvl：像素點所在的層
r：像素點所在的行
c：像素點所在的列
intvls：每組的層數
contr_thr：對比度閾值，針對歸一化後的圖像，用來去除不穩定特徵
返回值：返回經插值修正後的特徵點(feature類型)；若經有限次插值依然沒法精確到理想狀況
或者該點對比度太低，返回NULL*/
static struct feature* interp_extremum( IplImage*** dog_pyr, int octv, int intvl,
int r, int c, int intvls, double contr_thr ) //intvls=3; contr_thr=0.04
{
struct feature* feat;//修正後的特徵點
struct detection_data* ddata;//與特徵檢測有關的結構，存在feature結構的feature_data成員中
double xi, xr, xc, contr;//xi,xr,xc分別爲亞像素的intvl(層),row(y),col(x)方向上的
//增量(偏移量)
int i = 0;//插值次數
//SIFT_MAX_INTERP_STEPS指定了關鍵點的最大插值次數，即最多修正多少次，默認是5
while( i < SIFT_MAX_INTERP_STEPS )
{
//進行一次極值點差值，計算σ(層方向,intvl方向)，y，x方向上的子像素偏移量(增量)
interp_step( dog_pyr, octv, intvl, r, c, &xi, &xr, &xc );
//若在任意方向上的偏移量大於0.5時，意味着差值中心已經偏移到它的臨近點上，
//因此必須改變當前關鍵點的位置座標
if( ABS( xi ) < 0.5 && ABS( xr ) < 0.5 && ABS( xc ) < 0.5 )//若三方向上偏移量
//都小於0.5，表示已經夠精確，則不用繼續插值
break;
//修正關鍵點的座標，x,y,σ三方向上的原座標加上偏移量取整(四捨五入)
c += cvRound( xc );//x座標修正
r += cvRound( xr );//y座標修正
intvl += cvRound( xi );//σ方向，即層方向
//若座標修正後超出範圍，則結束插值，返回NULL
if( intvl < 1 || //層座標插以後越界
intvl > intvls ||
c < SIFT_IMG_BORDER || //行列座標插以後到邊界線內
r < SIFT_IMG_BORDER ||
c >= dog_pyr[octv][0]->width - SIFT_IMG_BORDER ||
r >= dog_pyr[octv][0]->height - SIFT_IMG_BORDER )
{
return NULL;
}
i++;
}
//若通過SIFT_MAX_INTERP_STEPS次插值後尚未修正到理想的精確位置，則返回NULL，
//即捨棄此極值點
if( i >= SIFT_MAX_INTERP_STEPS )
return NULL;
//計算被插值點的對比度：D + 0.5 * dD^T * X
contr = interp_contr( dog_pyr, octv, intvl, r, c, xi, xr, xc );
//此處與書上的不太同樣，RobHess此處的閾值使用的是0.04/S
if( ABS( contr ) < contr_thr / intvls )//若該點對比度太小，捨棄，返回NULL
return NULL;
//爲一個特徵點feature結構分配空間並初始化，返回特徵點指針
feat = new_feature();
//調用宏feat_detection_data來提取參數feat中的feature_data成員並轉換爲
//detection_data類型的指針
ddata = feat_detection_data( feat );
//將修正後的座標賦值給特徵點feat
//原圖中特徵點的x座標，由於第octv組中的圖的尺寸比原圖小2^octv倍，
//因此座標值要乘以2^octv,最後針對的是原圖的匹配
feat->img_pt.x = feat->x = ( c + xc ) * pow( 2.0, octv );
//原圖中特徵點的y座標，由於第octv組中的圖的尺寸比原圖小2^octv倍，
//因此座標值要乘以2^octv
feat->img_pt.y = feat->y = ( r + xr ) * pow( 2.0, octv );
ddata->r = r;//特徵點所在的行
ddata->c = c;//特徵點所在的列
ddata->octv = octv;//高斯差分金字塔中，特徵點所在的組
ddata->intvl = intvl;//高斯差分金字塔中，特徵點所在的組中的層
ddata->subintvl = xi;//特徵點在層方向(σ方向,intvl方向)上的亞像素偏移量
return feat;//返回特徵點指針
}
（5.2.6.1）cvRound()說明：
答：函數cvRound，cvFloor，cvCeil 都是用一種舍入的方法將輸入浮點數轉換成整數：
cvRound()：返回跟參數最接近的整數值，即四捨五入；
cvFloor()：返回不大於參數的最大整數值，即向下取整；
cvCeil()：返回不小於參數的最小整數值，即向上取整；
參看:【雜談opencv】OpenCV中的cvRound()、cvFloor()、 cvCeil()函數講解——https://blog.csdn.net/sinat_36264666/article/details/78849125
（5.2.6.2）interp_step代碼及說明：
答：
/*進行一次極值點差值，計算x，y，σ方向(層方向)上的子像素偏移量(增量)
參數：
dog_pyr：高斯差分金字塔
octv：像素點所在的組
intvl：像素點所在的層
r：像素點所在的行
c：像素點所在的列
xi：輸出參數，層方向上的子像素增量(偏移)
xr：輸出參數，y方向上的子像素增量(偏移)
xc：輸出參數，x方向上的子像素增量(偏移)*/
static void interp_step( IplImage*** dog_pyr, int octv, int intvl, int r, int c,
double* xi, double* xr, double* xc )
{
CvMat* dD, * H, * H_inv, X;
double x[3] = { 0 };
//在DoG金字塔中計算某點的x方向、y方向以及尺度方向上的偏導數，結果存放在列向量dD中
dD = deriv_3D( dog_pyr, octv, intvl, r, c );
//在DoG金字塔中計算某點的3*3海森矩陣
H = hessian_3D( dog_pyr, octv, intvl, r, c );
H_inv = cvCreateMat( 3, 3, CV_64FC1 );//海森矩陣的逆陣
cvInvert( H, H_inv, CV_SVD );
cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );
//X = - H^(-1) * dD，H的三個元素分別是x,y,σ方向上的偏移量(具體見SIFT算法說明)
cvGEMM( H_inv, dD, -1, NULL, 0, &X, 0 );
cvReleaseMat( &dD );
cvReleaseMat( &H );
cvReleaseMat( &H_inv );
*xi = x[2];//σ方向(層方向)偏移量
*xr = x[1];//y方向上偏移量
*xc = x[0];//x方向上偏移量
}
（5.2.6.2.1）deriv_3D代碼及說明：
答：
/*在DoG金字塔中計算某點的x方向、y方向以及尺度方向上的偏導數
參數：
dog_pyr：高斯差分金字塔
octv：像素點所在的組
intvl：像素點所在的層
r：像素點所在的行
c：像素點所在的列
返回值：返回3個偏導數組成的列向量{ dI/dx, dI/dy, dI/ds }^T*/
static CvMat* deriv_3D( IplImage*** dog_pyr, int octv, int intvl, int r, int c )
{
CvMat* dI;
double dx, dy, ds;
//求差分來代替偏導，這裏是用的隔行求差取中值的梯度計算方法
//求x方向上的差分來近似代替偏導數
dx = ( pixval32f( dog_pyr[octv][intvl], r, c+1 ) -
pixval32f( dog_pyr[octv][intvl], r, c-1 ) ) / 2.0;
//求y方向上的差分來近似代替偏導數
dy = ( pixval32f( dog_pyr[octv][intvl], r+1, c ) -
pixval32f( dog_pyr[octv][intvl], r-1, c ) ) / 2.0;
//求層間的差分來近似代替尺度方向上的偏導數
ds = ( pixval32f( dog_pyr[octv][intvl+1], r, c )
pixval32f( dog_pyr[octv][intvl-1], r, c ) ) / 2.0;
//組成列向量
dI = cvCreateMat( 3, 1, CV_64FC1 );
cvmSet( dI, 0, 0, dx );
cvmSet( dI, 1, 0, dy );
cvmSet( dI, 2, 0, ds );
return dI;
}
（5.2.6.2.1.1）問題issue：差分代替偏導？一階二階都使用嗎？怎麼操做？
差分代替偏導！一階狀況：
一階差分代替一階偏導，即梯度，[-1，0，1]
行方向[-1，0，1]
-1
列方向[ 0 ]
1
二階狀況：
參看：圖像處理中的一階偏導數和二階偏導數——https://www.docin.com/p-749102193.html
（5.2.6.2.1.2）cvmSet說明
答：opencv中cvmSet爲逐點賦值
cvmSet( dI, 0, 0, dx ); dx
cvmSet( dI, 1, 0, dy ); dI = [ dy ]
cvmSet( dI, 2, 0, ds ); ds
（5.2.6.2.2）hessian_3D代碼及說明：
答：//差分代替偏導，二階狀況
/*在DoG金字塔中計算某點的3*3海森矩陣
/ Ixx Ixy Ixs \
| Ixy Iyy Iys |
\ Ixs Iys Iss /
參數：
dog_pyr：高斯差分金字塔
octv：像素點所在的組
intvl：像素點所在的層
r：像素點所在的行
c：像素點所在的列
返回值：返回3*3的海森矩陣
*/
static CvMat* hessian_3D( IplImage*** dog_pyr, int octv, int intvl, int r, int c )
{
CvMat* H;
double v, dxx, dyy, dss, dxy, dxs, dys;
v = pixval32f( dog_pyr[octv][intvl], r, c );//該點的像素值
//用差分近似代替倒數(具體公式見各類梯度的求法)
//dxx = f(i+1,j) - 2f(i,j) + f(i-1,j)
//dyy = f(i,j+1) - 2f(i,j) + f(i,j-1)
dxx = ( pixval32f( dog_pyr[octv][intvl], r, c+1 ) +
pixval32f( dog_pyr[octv][intvl], r, c-1 ) - 2 * v );
dyy = ( pixval32f( dog_pyr[octv][intvl], r+1, c ) +
pixval32f( dog_pyr[octv][intvl], r-1, c ) - 2 * v );
dss = ( pixval32f( dog_pyr[octv][intvl+1], r, c ) +
pixval32f( dog_pyr[octv][intvl-1], r, c ) - 2 * v );
dxy = ( pixval32f( dog_pyr[octv][intvl], r+1, c+1 ) -
pixval32f( dog_pyr[octv][intvl], r+1, c-1 ) -
pixval32f( dog_pyr[octv][intvl], r-1, c+1 ) +
pixval32f( dog_pyr[octv][intvl], r-1, c-1 ) ) / 4.0;
dxs = ( pixval32f( dog_pyr[octv][intvl+1], r, c+1 ) -
pixval32f( dog_pyr[octv][intvl+1], r, c-1 ) -
pixval32f( dog_pyr[octv][intvl-1], r, c+1 ) +
pixval32f( dog_pyr[octv][intvl-1], r, c-1 ) ) / 4.0;
dys = ( pixval32f( dog_pyr[octv][intvl+1], r+1, c ) -
pixval32f( dog_pyr[octv][intvl+1], r-1, c ) -
pixval32f( dog_pyr[octv][intvl-1], r+1, c ) +
pixval32f( dog_pyr[octv][intvl-1], r-1, c ) ) / 4.0;
//組成海森矩陣
H = cvCreateMat( 3, 3, CV_64FC1 );
cvmSet( H, 0, 0, dxx );
cvmSet( H, 0, 1, dxy );
cvmSet( H, 0, 2, dxs );
cvmSet( H, 1, 0, dxy );
cvmSet( H, 1, 1, dyy );
cvmSet( H, 1, 2, dys );
cvmSet( H, 2, 0, dxs );
cvmSet( H, 2, 1, dys );
cvmSet( H, 2, 2, dss );
return H;
}
（5.2.6.2.3）cvInvert函數說明：
答：
double cvInvert(//矩陣取逆
const CvArr* src,//目標矩陣
CvArr* dst,//結果矩陣
int method = CV_LU//逆運算方法
);
其中method有
方法的參數值 含義
CV_LU 高斯消去法
CV_SVD 奇異值分解
CV_SVD_SYM 對稱矩陣的SVD
參看：《學習opencv》筆記——矩陣和圖像操做——cvInRange,cvInRangeS,cvInvert and cvMahalonobis——https://blog.csdn.net/zhurui_idea/article/details/28630589
（5.2.6.2.4）問題issue：cvInitMatHeader的一些問題：
答：
1. cvInitMatHeader的用法很古怪。第一個參數必須是CvMat格式的，官方的文檔是CvMat *mat，可是這裏要注意，mat必須是初始化過的，單獨定義一個指針如CvMat *data1；把data1帶入cvInitMatHeader函數後，編譯器會報錯，顯示沒有初始化！下面的代碼纔是正確的
2.下面的代碼我在使用的時候想生成兩個mat，而後被另外一個函數f調用，發現f調用時失敗的，用下面代碼生成的數據被釋放掉了。因此這種用法只能在一個函數體內使用，跨函數是不行的。
CvMat data1; CvMat responses1; //將數組中數據轉化成opencv支持的mat格式 cvInitMatHeader(&data1, total, var_count, CV_32FC1, pData); cvInitMatHeader(&responses1, total, 1, CV_32SC1, pRespon);
參看：1）【OpenCV學習】矩陣cvInitMatHeader和cvCreateMat：https://blog.csdn.net/cc1949/article/details/22476251
2）cvInitMatHeader的一些問題：https://www.xuebuyuan.com/1050480.html
3）【OpenCV學習】矩陣cvInitMatHeader和cvCreateMat：舉例說明看——https://blog.csdn.net/cc1949/article/details/22476251
（5.2.6.2.5）cvGEMM函數說明：
答：void cvGEMM( const CvArr* src1, const CvArr* src2, double alpha, const CvArr* src3, double beta, CvArr* dst, int tABC=0 );
這是通用矩陣乘法，其中各個參數表示：
src1：第一輸入數組
src2：第二輸入數組
alpha：係數
src3「第三輸入數組（偏移量），若是沒有偏移量，能夠爲空（NULL）
beta：表示偏移量的係數
dst：輸出數組
tABC：轉置操做標誌，能夠是0。當爲0時，沒有轉置。
參看博客：cvGEMM、cvMatMul和cvMatMulAdd的定義——https://blog.csdn.net/liulianfanjianshi/article/details/11737921
（5.2.6.3）interp_contr代碼及說明：
答：
/*計算被插值點的對比度：D + 0.5 * dD^T * X
參數：
dog_pyr：高斯差分金字塔
octv：像素點所在的組
intvl：像素點所在的層
r：像素點所在的行
c：像素點所在的列
xi：層方向上的子像素增量
xr：y方向上的子像素增量
xc：x方向上的子像素增量
返回值：插值點的對比度*/
static double interp_contr( IplImage*** dog_pyr, int octv, int intvl, int r,
int c, double xi, double xr, double xc )
{
CvMat* dD, X, T;
double t[1], x[3] = { xc, xr, xi };
//偏移量組成的列向量X，其中是x，y，σ三方向上的偏移量
cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );
//矩陣乘法的結果T，是一個數值
cvInitMatHeader( &T, 1, 1, CV_64FC1, t, CV_AUTOSTEP );
//在DoG金字塔中計算某點的x方向、y方向以及尺度方向上的偏導數，結果存放在列向量dD中
dD = deriv_3D( dog_pyr, octv, intvl, r, c );
//矩陣乘法：T = dD^T * X
cvGEMM( dD, &X, 1, NULL, 0, &T, CV_GEMM_A_T );
cvReleaseMat( &dD );
//返回計算出的對比度值：D + 0.5 * dD^T * X (具體公式推導見SIFT算法說明)
return pixval32f( dog_pyr[octv][intvl], r, c ) + t[0] * 0.5;
}
（5.2.6.4）new_feature()代碼及說明：
答：
/*爲一個feature結構分配空間並初始化
返回值：初始化完成的feature結構的指針*/
static struct feature* new_feature( void )
{
struct feature* feat;//特徵點指針
struct detection_data* ddata;//與特徵檢測相關的結構
feat = malloc( sizeof( struct feature ) );//分配空間
memset( feat, 0, sizeof( struct feature ) );//清零
ddata = malloc( sizeof( struct detection_data ) );
memset( ddata, 0, sizeof( struct detection_data ) );
feat->feature_data = ddata;//將特徵檢測相關的結構指針賦值給特徵點的feature_data成員
feat->type = FEATURE_LOWE;//默認是LOWE類型的特徵點
return feat;
}
（5.2.6.5）pow()說明：
答：C 庫函數 double pow(double x, double y) 返回 x 的 y 次冪，即 x^y。
參看:C 庫函數 - pow()——https://www.runoob.com/cprogramming/c-function-pow.html
（5.2.7）is_too_edge_like代碼及說明：
/*去除邊緣響應，即經過計算主曲率比值判斷某點是否邊緣點
參數：
dog_img：此特徵點所在的DoG圖像
r：特徵點所在的行
c：特徵點所在的列
cur_thr：主曲率比值的閾值，用來去除邊緣特徵
返回值：0：此點是非邊緣點；1：此點是邊緣點*/
static int is_too_edge_like( IplImage* dog_img, int r, int c, int curv_thr )
{
double d, dxx, dyy, dxy, tr, det;
/*某點的主曲率與其海森矩陣的特徵值成正比，爲了不直接計算特徵值，這裏只考慮
特徵值的比值可經過計算海森矩陣的跡tr(H)和行列式det(H)來計算特徵值的比值
設a是海森矩陣的較大特徵值，b是較小的特徵值，有a = r*b，r是大小特徵值的比值
tr(H) = a + b; det(H) = a*b;
tr(H)^2 / det(H) = (a+b)^2 / ab = (r+1)^2/r
r越大，越多是邊緣點；伴隨r的增大，(r+1)^2/r 的值也增大，因此可經過(r+1)^2/r 判斷
主曲率比值是否知足條件*/
/* principal curvatures are computed using the trace and det of Hessian */
d = pixval32f(dog_img, r, c);//調用函數pixval32f獲取圖像dog_img的第r行第c列的點的座標值
//用差分近似代替偏導，求出海森矩陣的幾個元素值
/* / dxx dxy \
\ dxy dyy / */
dxx = pixval32f( dog_img, r, c+1 ) + pixval32f( dog_img, r, c-1 ) - 2 * d;
dyy = pixval32f( dog_img, r+1, c ) + pixval32f( dog_img, r-1, c ) - 2 * d;
dxy = ( pixval32f(dog_img, r+1, c+1) - pixval32f(dog_img, r+1, c-1) -
pixval32f(dog_img, r-1, c+1) + pixval32f(dog_img, r-1, c-1) ) / 4.0;
tr = dxx + dyy;//海森矩陣的跡
det = dxx * dyy - dxy * dxy;//海森矩陣的行列式
//若行列式爲負，代表曲率有不一樣的符號，去除此點
/* negative determinant -> curvatures have different signs; reject feature */
if( det <= 0 )
return 1;//返回1代表此點是邊緣點
//經過式子：(r+1)^2/r 判斷主曲率的比值是否知足條件，若小於閾值，代表不是邊緣點
if( tr * tr / det < ( curv_thr + 1.0 )*( curv_thr + 1.0 ) / curv_thr )
return 0;//不是邊緣點
return 1;//是邊緣點
}
（5.2.8）cvSeqPush代碼及說明：
函數原型： char * cvSeqPush(CvSeq *seq, void *element = NULL);//添加元素到序列的尾部
參看：opencv建立序列cvCreateSeq與插入元素cvSeqPush的運用：https://blog.csdn.net/gdut2015go/article/details/46494677
（5.3）calc_feature_scales代碼及說明：
/*計算特徵點序列中每一個特徵點的尺度
參數：
features：特徵點序列
sigma：初始高斯平滑參數，即初始尺度
intvls：尺度空間中每組的層數*/
static void calc_feature_scales( CvSeq* features, double sigma, int intvls )
{
struct feature* feat;
struct detection_data* ddata;
double intvl;
int i, n;
n = features->total;//總的特徵點個數
//遍歷特徵點
for( i = 0; i < n; i++ )
{
//調用宏，獲取序列features中的第i個元素，並強制轉換爲struct feature類型
feat = CV_GET_SEQ_ELEM( struct feature, features, i );
//調用宏feat_detection_data來提取參數feat中的feature_data成員並轉換爲
//detection_data類型的指針
ddata = feat_detection_data( feat );
//特徵點所在的層數ddata->intvl加上特徵點在層方向上的亞像素偏移量，獲得
//特徵點的較爲精確的層數
intvl = ddata->intvl + ddata->subintvl;
//計算特徵點的尺度(公式見SIFT算法說明)，並賦值給scl成員
feat->scl = sigma * pow( 2.0, ddata->octv + intvl / intvls );
//計算特徵點所在的組的尺度，給detection_data的scl_octv成員賦值
ddata->scl_octv = sigma * pow( 2.0, intvl / intvls );
}
}
（5.3.1）CV_GET_SEQ_ELEM函數說明：
答案：用法：從所給序列中取出元素的地址，注意：獲得的是地址，即指針
因此，關鍵在於序列中存放的是那種類型的數據，若存放的爲地址，那用這個宏獲得的就是指針的指針。
例程：
1.序列中存放的爲CvPoint
CvPoint pt=*CV_GET_SEQ_ELEM(CvPoint,hull,i);
2.序列中存放的爲CvPoint*，即指針
CvPoint* pt=*CV_GET_SEQ_ELEM(CvPoint*,hull,i);
或者
CvPoint pt=**CV_GET_SEQ_ELEM(CvPoint*,hull,i);
參看：關於宏CV_GET_SEQ_ELEM——https://blog.csdn.net/u013089125/article/details/20126215
（5.4）adjust_for_img_dbl代碼及說明：
/*將特徵點序列中每一個特徵點的座標減半(當設置了將圖像放大爲原圖的2倍時，特徵點檢測
完以後調用)
參數：
features：特徵點序列*/
static void adjust_for_img_dbl( CvSeq* features )
{
struct feature* feat;
int i, n;
n = features->total;//總的特徵點個數
//遍歷特徵點
for( i = 0; i < n; i++ )
{
//調用宏，獲取序列features中的第i個元素，並強制轉換爲struct feature類型
feat = CV_GET_SEQ_ELEM( struct feature, features, i );
//將特徵點的x,y座標和尺度都減半
feat->x /= 2.0;
feat->y /= 2.0;
feat->scl /= 2.0;
feat->img_pt.x /= 2.0;
feat->img_pt.y /= 2.0;
}
}
附：
1.問題描述：C版本的cvGetTickFrequency()函數和C++版本的getTickFrequency()計算時間的公式不同?
解決方案：參看博客：https://blog.csdn.net/chaipp0607/article/details/71056580
注意：OpenCV C版本的cvGetTickFrequency()函數和C++版本的getTickFrequency()的單位不同