SIFT 特征點提取算法

SIFT特征點相對於ORB計算速度較慢，在沒有GPU加速情況下，無法滿足視覺里程計的實時性要求，或者無法運行在手機平台上，但是效果更好，精度更高。在應用時可以擇優選取，了解其本質原理的動機是為了自己使用時，可以對其進行修改，針對自己的應用場景優化算法。

有足夠的時間，可以去看D. Lowe的論文，理解起來更透徹．

1. 用高斯核構建尺度空間

對於構建的高斯金字塔，金字塔每層多張圖像合稱為一組（Octave），每組有多張（也叫層Interval）圖像。通常高斯金字塔最底層為原始圖像第0組，octave之間為降采樣過程，對應OpenCV里的函數為PryDown()。注意這里降采樣PryDown在直觀上是金字塔向上走，要注意區分。

金字塔的第i組octave通過對第i-1組降采樣獲得（通常降采樣的比例為2，也就是高斯模糊后去掉偶數行/列像素）。在每一組octave中還需要使用高斯核σ創建連續的尺度空間圖像: σ, kσ, k²σ,..., k^n-1σ，其中，k=2^1/s，s為每一組octave中再分的尺度層數。構建連續尺度空間的目的是為了提取特征點（角點）時，該特征點不僅在二維圖像空間中是梯度極值點，還需要在第三維的尺度空間也是極值點。

當然，檢測角點並不是直接在高斯卷積構建的尺度空間中，而是DoG（Difference of Gaussian）高斯差分尺度空間中，因此，為了保證DOG高斯差分尺度空間變化的連續性（是高斯卷積空間相鄰尺度的差），需要在每一層octave的高斯卷積尺度首尾兩邊多創建一個尺度。

這里注意理解連續性，不是微積分里面的連續性，指的是尺度空間連續，(σ, kσ, k²σ,..., k^s-1σ) ,(2σ, 2kσ, 2k²σ,..., 2k^s-1σ) 。從這里可以看出，尺度空間主要由σ主導．

2. 在DoG尺度空間中檢測極值點

　　雖然在高斯拉普拉斯LoG尺度空間中檢測極值點是最精確的，但是由於計算量比較大，通常使用DoG尺度空間對其進行近似，在DoG尺度空間中進行局部最值搜索。有關LoG可以參考http://fourier.eng.hmc.edu/e161/lectures/gradient/node8.html．LoG也就是將拉普拉斯算子作用於高斯核函數后形成的新的核函數，拉普拉斯算子的作用是求二階偏導，對噪點響應較敏感，因此需要先用高斯核對圖像進行平滑．

　　由於DoG尺度空間是離散的，我們只能通過比較DoG空間中某個點周圍26個點來獲取一個點是否可以作為極值(和ORB一樣，這里可以選較少的幾個先比較剔除)．

　　　　　　　　

但是為了進一步獲得精確的亞像素級別特征點位置，需要將DoG函數D(x,y,σ)在像素極值點附近二階展開，對亞像素極值求偏導后＝０求解．

　　　　

　　　　

第二式帶入第一式后可以得到亞像素極值點對應的DoG函數值，也就是對比度contrast：

　　　　

 

　　在檢測到DoG尺度空間中的極值點(extremum)后，使用兩個閾值來剔除質量不高（unstable）的點，contrastThreshold以及edgeThreshold．首先在DoG尺度空間中在potential的極值點附近二階泰勒展開，尋找更精確的極值點坐標，其中，contrastThreshold是為了剔除對比度不高的極值點，edgeThreshold是為了剔除邊緣點，因為邊緣點也滿足極值點的條件，但是不是需要的角點．

剔除對比度不高（不穩定）的點，|D(x)| < 0.03，這里假設DoG函數取值范圍[0, 1]．

剔除邊緣點借鑒了Harris Corner的檢測方法，計算特征點鄰域的Hessian矩陣。

　　　　

一次求導物理意義是變化率／梯度，二次求導物理意義是曲率，因此Hessian矩陣更能體現一個邊緣的特征程度。

該Hessian矩陣的特征值代表了邊緣特征的主曲率。兩個特征值之間的比例體現了鄰域兩個垂直方向的曲率差別，這里去掉比例較大（橫跨邊緣和沿着邊緣方向的曲率差別較大，>10）的點。注意這里不用特征值分解，而是用矩陣跡和行列式的性質，得到兩個特征值比例就行．

3. 計算特征點方向

在檢測到的特征點一個圓形鄰域像素集合上計算各自的梯度幅值和梯度方向：

　　　　　 

將360度方向划分為方向直方圖中的36個bins，獲取直方圖的peak，作為特征點的方向。如果存在大於80%*peak幅值的方向，則在同一個尺度特征點(x,y,σ)上新建一個方向不同的特征點，有多少個新建多少個特征點．論文中說，相同位置的不會超過15%, 但是可以顯著增強匹配的穩定性．注意這里的鄰域是有區分度的，距離特征點近的像素在計算直方圖時權重較大，使用高斯分布加權處理。

獲得關鍵點主方向后，每個關鍵點有三個信息(x,y,σ,θ)：位置、尺度、方向。

 

４. 計算128維描述子

在對鄰域像素點采樣前，要旋轉xy坐標系和關鍵點的方向對齊，使得特征點具備方向不變的特征（旋轉相機后，特征點描述子不變）。

將旋轉后區域划分為4×4=16個子區域（每個區域間隔為mσ像元），在每個子區域內計算8個方向的梯度直方圖，繪制每個方向梯度方向的累加值，形成一個種子點。

與求主方向不同的是，此時，每個子區域梯度方向直方圖將0°~360°划分為8個方向區間，每個區間為45°。即每個種子點有8個方向區間的梯度強度信息。由於存在d×d，即4×4個子區域，所以最終共有4×4×8=128個數據，形成128維SIFT特征矢量。

計算128維空間中描述子的距離，SIFT的實現中計算的是歐氏距離.

 

練習：

OpenCV中提供的goodFeatureToTrack()可以使用Harris和Shi-Tomasi角點檢測算法對圖像進行角點提取，注意沒有構建尺寸空間，且沒有提供描述子計算，是最簡單的角點提取算法，但是通過對其參數的調整，可以了解一些實現特點：

#include "common.h"
#include <opencv2/nonfree/features2d.hpp>

using namespace std;
using namespace cv;

Mat origin, img;
vector<Point2f> keyPoints;
string title("FeatureDetector");
int maxCorners(10);
int qualityLevel(2);
int minDistance(10);
int blockSize(3);
int useHarrisDetector(0);
int k(4);
int method(1);

void detectCorner(int, void*);

int main(){
  origin = imread("/home/shang/dataset/opencv/building.jpg",CV_LOAD_IMAGE_COLOR);
  if(!origin.data){
    cerr << "No data!"<< endl;
    return -1;
  }
  cvtColor(origin, img, CV_BGR2GRAY);
  namedWindow(title);
  createTrackbar("maxCorners", title, &maxCorners, 200, detectCorner);
  createTrackbar("qualityLevel (%)", title, &qualityLevel, 100, detectCorner);
  createTrackbar("minDistance", title, &minDistance, 400, detectCorner);
  createTrackbar("blockSize", title, &blockSize, 10, detectCorner);
  createTrackbar("useHarrisDetector", title, &useHarrisDetector, 1, detectCorner);
  createTrackbar("k when in HarrisDetector", title, &k, 100, detectCorner);
  createTrackbar("Method", title, &method, 2, detectCorner);
  detectCorner(0,0);
  imshow(title, img);
  
  while(true){
    if(waitKey(0)==27)
      break;
  }
}

void detectCorner(int, void*){
  Mat result ;
  img.copyTo(result);
  switch (method){
    case 1:
      if(qualityLevel==0) {
    qualityLevel = 1;
    setTrackbarPos("qualityLevel (%)", title, 1);
      }
      if(blockSize==0) {
    blockSize = 1;
    setTrackbarPos("blockSize", title, 1);
      }
      goodFeaturesToTrack(result, 
              keyPoints, 
              maxCorners, 
              qualityLevel/100.0, 
              minDistance, 
              noArray(), 
              blockSize,
              useHarrisDetector,
              k/100.0
            );
    
      cvtColor(result, result, CV_GRAY2RGB);
      for(vector<Point2f>::const_iterator it=keyPoints.begin(); it!=keyPoints.end(); it++){
    circle(result, *it, 8, Scalar(0,0,255),2);
      }
      imshow(title, result);
      break;
    case 2:
      SiftFeatureDetector detector(100);
      vector<KeyPoint> keyPoints;
      detector.detect(result, keyPoints);
      cvtColor(result, result, CV_GRAY2RGB);
      for(vector<KeyPoint>::const_iterator it=keyPoints.begin(); it!=keyPoints.end(); it++){
    rectangle(result, Rect(it->pt.x-8,it->pt.y-8,16,16), Scalar(0,255,0),2);
      }
      imshow(title, result);
    break;
  }
  
}

 

參考：

SIFT特征提取分析

【OpenCV】SIFT原理與源碼分析