SIFT特征提取和代碼分析

本文隻記錄sift特征提取過程和sift的擴充應用，并分析了opensift的代碼。如果想詳細了解sift的理論知識請參見Rachel-Zhang的文章。這裡沒分析OpenCV的代碼，是因為相比之下opensift代碼結構更加清楚，可讀性更好。

一、SIFT提取過程

1. 對圖像寬高放大1倍，并假定圖像已被0.5高斯濾波過，為了達到初始為1.6高斯的效果，再用 1.62−0.52−−−−−−−−−√ 高斯濾波一次，公式參考wiki-高斯模糊；
2. 建立高斯金字塔，塔數 o=log2min(w,h) ，每塔含6層，高斯方差為 σ,kσ,k2σ,k3σ,k4σ,k5σ ，下一個塔時将圖像寬高縮小1倍再高斯濾過；
3. 高斯差分，對每個塔内各層做差，形成DoG圖像；
4. 對DoG圖像求局部極值點并亞像素化，并用 X^=−∂2D−1∂X2∂D∂X 去疊代 X ，這裡的理論本質是利用牛頓疊代法來最優化D，可參見李航的《統計學習方法》講牛頓疊代法那部分；
5. 過濾過 D 值太小的點；
6. 去除邊緣響應的點，利用Hessian矩陣模和迹的關系設計濾波器，H=[DxxDxyDxyDyy]，這裡區分一下harris角點中的自相關矩陣。圖像二階導數為 Dxx=Dx+1+Dx−1−2Dx ， Dxy=(Dx+1,y+1−Dx+1,y−1−Dx−1,y+1+Dx−1,y−1)/4 ；
7. 計算關鍵點視窗下的梯度直方圖，視窗大小為 3×1.5σoct ，高斯權重中方差為 1.5σoct ，其中 oct 為關鍵點所在塔内的尺度，注意這裡的尺度已亞像素化了，并且累加梯度強度時未采用線性插值法，計算比較粗略，但在後面計算描述子時就非常精細了；
8. 計算梯度直方圖的峰值，取峰值的80%作為門檻值，直方圖中大于門檻值的極值點會成為新關鍵點，并附上該方向加後到特征點清單中，是以一個關鍵點可以會多次加到關鍵點清單中，但關鍵點的附屬方向不一樣；
9. 計算關鍵點旋轉視窗下的梯度直方圖，視窗大小為 3σoct×2√×(d+1)+12 ，其中 d 表示描述子的寬度為4。計算時以關鍵點為中點，從(−r,−r)、 (−r+1,−r) 、…、 (r,r) 坐标進行旋轉 [cosθsinθ−sinθcosθ]×[xy] ，然後 [x′y′] 映射到 4×4 的網格内，角度值為 tan−1dydx− 主方向角，累加梯度強度時采用了線性插值法，最後形成 4×4×8=128 維向量；
10. 對特征向量歸一化，可以用 L2−Hys 歸一化。此時所有關鍵點的描述子計算完成；
11. 對所有關鍵點按照所在尺度（ σ×2oct+kS ，其中 k 已被亞像素化過）進行降序排序，這是為了友善以後關鍵點比對。

二、擴充知識

• SIFT算法的擴充
  • PCA-SIFT，其它步驟都一樣，隻在特征描述階段有所差別，具體為：對每個特征點周圍選擇一個大小為41*41像素的像塊（主方向旋轉後的），計算垂直和水準的梯度，形成39*39*2=3042的向量，對所有的特征點做成一個矩陣k*3042做PCA，選取n個特征向量做成投影矩陣。一般為36維，比對速度更快，但區分度下降，并且延長了特征的計算時間。
  • GLOH（Gradient Location-Orientation Histogram），也隻是在特征描述階段有所差別，具體為：把原來SIFT中4*4棋盤格的子塊改成放射狀的3個同心圓，對大的兩個同心圓等分8份，共形成8+8+1=17的區域，直方圖16維，再做PCA。區分度更高但是資料壓縮花銷時間太長。
• SIFT的一些應用
  • sparse sift，做特征比對，有圖像拼接的應用，圖像拷貝檢測；
  • dense sift，用BoW模型，做圖像分類，目辨別别（連結以後加上）；

三、代碼分析

SIFT特征檢測的資料結構

struct detection_data
{
  int r;            // 所在圖像中的行号
  int c;            // 所在圖像中的列号
  int octv;         // 所在尺度空間的金字塔數，octv越大圖像越小
  int intvl;        // 所在金字塔中的層數
  double subintvl;  // 亞像素化時金字塔中的層數的小數位
  double scl_octv;  // sigma * pow(2, (intvl + subintval) / intvls)
};
           

SIFT關鍵點檢測與描述的代碼結構

int _sift_features(...)
{
    init_img = create_init_img(...);    // 圖像寬高放大倍，假定圖像初始化尺度為
    gauss_pyr = build_gauss_pyr(...);   // 建立高斯金字塔圖像，存在gauss_pyr[octvs][intval+]中，函數中為了減少使用pow()，做了優化，見補充
    dog_pyr = build_dog_pyr(...);       // 建立高斯差分金字塔
    features = scale_space_extrema(...) // 尋找高斯差分空間極大值點，完成了亞像素化和邊緣響應
    calc_feature_scales(...);           // 計算關鍵點的各尺度，即關鍵點的尺度值=sigma*pow(, octv + (intvl + subintvl) / intvls)
    if(img_dbl)                         // 如果初始化時圖像double過，則需要把關鍵點的坐标縮小
        adjust_for_img_dbl( features );
    calc_feature_oris(...);             // 計算關鍵點主方向，此時一個關鍵點有多個主方向，則會複制産生新的關鍵點，這裡計算的直方圖未插值，但最後對直方圖做了平滑
    compute_descriptors(...);           // 對關鍵點進行描述
    cvSeqSort(...);                     // 對關鍵點按尺度降序排序
}
           

補充

1、計算高斯金字塔時的trick，普通情況下是利用σk=sigma×2ks核去高斯模糊，而代碼中用了遞推式 sig[i]=(kiσ)2−(ki−1σ)2−−−−−−−−−−−−−−√=ki−1σk2−1−−−−−√ ，是以這個樣遞推 sig[0]=σ 、 sig[1]=σk2−1−−−−−√ 、 sig[i]=k×sig[i−1]