經典演算法研究系列：九、影象特徵提取與匹配之SIFT演算法
經典演算法研究系列：九、SIFT演算法研究
作者:July、二零一一年二月十五日。
推薦閱讀：
David G. Lowe, "Distinctive image features from scale-invariant keypoints,"
International Journal of Computer Vision, 60, 2 (2004), pp. 91-110
---------------------------------------------
尺度不變特徵轉換(Scale-invariant feature transform 或 SIFT)是一種電腦視覺的演算法用來偵測與描述影像中的區域性性特徵，它在空間尺度中尋找極值點，並提取出其位置、尺度、旋轉不變數，此演算法由 David Lowe 在1999年所發表，2004年完善總結。
Sift演算法就是用不同尺度（標準差）的高斯函式對影象進行平滑，然後比較平滑後圖像的差別，
差別大的畫素就是特徵明顯的點。
一、Sift演算法的步驟
Sift（Scale Invariant Feature Transform）是一個很好的影象匹配演算法，
同時能處理亮度、平移、旋轉、尺度的變化，利用特徵點來提取特徵描述符，最後在特徵描述符之間尋找
匹配。
該演算法主要包括5個步驟進行匹配：
1、構建尺度空間，檢測極值點，獲得尺度不變性；
2、特徵點過濾並進行精確定位，剔除不穩定的特徵點；
3、在特徵點處提取特徵描述符，為特徵點分配方向值；
4、生成特徵描述子，利用特徵描述符尋找匹配點；
以特徵點為中心取16*16的鄰域作為取樣視窗，
將取樣點與特徵點的相對方向通過高斯加權後歸入包含8個bin的方向直方圖，
最後獲得4*4*8的128維特徵描述子。
示意圖如下：
5、計算變換引數。
當兩幅影象的Sift特徵向量生成以後，下一步就可以採用關鍵點特徵向量的歐式距離來作為兩幅影象中
關鍵點的相似性判定度量。
取圖1的某個關鍵點，通過遍歷找到影象2中的距離最近的兩個關鍵點。
在這兩個關鍵點中，如果次近距離除以最近距離小於某個闕值，則判定為一對匹配點。
最後，看下Sift 演算法效果圖：
下圖左邊部分Sift演算法匹配結果，右邊部分是其它演算法匹配結果：
二、Sift演算法的描述
在上述的Sift演算法步驟一中，提到了尺度空間，那麼什麼是尺度和尺度空間呢？
尺度就是受delta這個引數控制的表示。
而不同的L（x,y,delta)就構成了尺度空間，實際上，具體計算的時候，即使連續的高斯函式，都要被離
散為（一般為奇數大小）(2*k+1) *(2*k+1)矩陣，來和數字影象進行卷積運算。
David Lowe關於Sfit演算法，2004年發表在Int. Journal of Computer Vision的經典論文中，
對尺度空間（scal space)是這樣定義的 ：
It has been shown by Koenderink (1984) and Lindeberg (1994) that under a variety of
reasonable assumptions the only possible scale-space kernel is the Gaussian function.
Therefore,the scale space of an image is defined as a function, L(x; y; delta) that is
produced from the convolution of a variable-scale Gaussian, G(x; y; delta), with an input
image, I(x; y):
因此 ，一個影象的尺度空間，L（x,y,delta) ,
定義為原始影象I (x,y)與一個可變尺度的2維高斯函式G(x,y,delta) 卷積運算。
即，原始影像I(x，y)在不同的尺度e下，與高斯濾波器G(x,y,e)進行卷積，得到L(x,y,e)，如下：
L(x,y,e) = G(x,y,e)*I(x,y)
其中G(x,y,e)是尺度可變高斯函式，
G(x,y,e) = [1/2*pi*e2] * exp[ -(x2 + y2)/2e2]
（x，y）是空間座標， e是尺度座標。
為了更有效的在尺度空間檢測到穩定的關鍵點，提出了高斯差分尺度空間（DOG scale-space）。
利用不同尺度的高斯差分核與原始影象I(x,y) ，卷積生成。
D(x,y,e) = ((G(x,y,ke) - G(x,y,e)) * I(x,y)
= L(x,y,ke) - L(x,y,e)
DOG運算元計算簡單，是尺度歸一化的LoG運算元的近似。
Gaussian卷積是有尺寸大小的，使用同一尺寸的濾波器對兩幅包含有不同尺寸的同一物體的影象求區域性最值將有可能出現一方求得最值而另一方卻沒有的情況，但是容易知道假如物體的尺寸都一致的話它們的區域性最值將會相同。
SIFT的精妙之處在於採用影象金字塔的方法解決這一問題，我們可以把兩幅影象想象成是連續的，分別以它們作為底面作四稜錐，就像金字塔，那麼每一個 截面與原影象相似，那麼兩個金字塔中必然會有包含大小一致的物體的無窮個截面，但應用只能是離散的，所以我們只能構造有限層，層數越多當然越好，但處理時 間會相應增加，層數太少不行，因為向下取樣的截面中可能找不到尺寸大小一致的兩個物體的影象。
有了影象金字塔就可以對每一層求出區域性最值，但是這樣的穩定 點數目將會十分可觀，所以需要使用某種方法抑制去除一部分點，但又使得同一尺度下的穩定點得以儲存
影象金字塔的構建：影象金字塔共O組，每組有S層，下一組的影象由上一組影象降取樣得到。
如下圖：
三、Sift演算法的實現
作為一種匹配能力較強的區域性描述運算元，SIFT演算法的實現相當複雜，
不過David Lowe到底也還是用c++實現了它，標頭檔案裡就下述兩個關鍵函式，下面具體闡述下。
函式一：
int sift_features( IplImage* img, struct feature** feat )
這個函式就是用來提取影象中的特徵向量。
引數img為一個指向IplImage資料型別的指標，用來表示需要進行特徵提取的影象。
IplImage是opencv庫定義的影象基本型別（關於opencv是一個著名的影象處理類庫，詳細的介紹可以參見（http://www.opencv.org.cn）。
引數feat 是一個數組指標，用來儲存影象的特徵向量。
函式呼叫成功將返回特徵向量的數目，否則返回-1.
函式，完整表述如下：
int sift_features( IplImage* img, struct feature** feat )
{
return _sift_features( img, feat, SIFT_INTVLS, SIFT_SIGMA, SIFT_CONTR_THR,
SIFT_CURV_THR, SIFT_IMG_DBL, SIFT_DESCR_WIDTH,
SIFT_DESCR_HIST_BINS );
}如你所見，函式一sift_features呼叫下面的函式二。
函式二：
int _sift_features( IplImage* img, struct feature** feat, int intvls,double sigma, double
contr_thr, int curv_thr, int img_dbl, int descr_width, int descr_hist_bins )
稍微介紹下此函式的幾個引數：
intvls: 每個尺度空間的取樣間隔數，預設值為3.
sigma: 高斯平滑的數量，預設值1.6.
contr_thr:判定特徵點是否穩定，取值（0，1），預設為0.04，這個值越大，被剔除的特徵點就越多。
curv_thr:判定特徵點是否邊緣點，預設為6.
img_dbl:在建立尺度空間前如果影象被放大了1倍則取值為1，否則為0.
descr_width:計算特徵描述符時鄰域子塊的寬度，預設為4.
descr_hist_bins:計算特徵描述符時將特徵點鄰域進行投影的方向數，預設為8，分別是0，45，90，135
，180，215，270，315共8個方向。
以下是此函式二的完整表述：
int _sift_features( IplImage* img, struct feature** feat, int intvls,
double sigma, double contr_thr, int curv_thr,
int img_dbl, int descr_width, int descr_hist_bins )
{
IplImage* init_img;
IplImage*** gauss_pyr, *** dog_pyr;
CvMemStorage* storage;
CvSeq* features;
int octvs, i, n = 0;
/* check arguments */
if( ! img )
fatal_error( "NULL pointer error, %s, line %d", __FILE__, __LINE__ );
if( ! feat )
fatal_error( "NULL pointer error, %s, line %d", __FILE__, __LINE__ );
/* build scale space pyramid; smallest dimension of top level is ~4 pixels */
init_img = create_init_img( img, img_dbl, sigma );
octvs = log( MIN( init_img->width, init_img->height ) ) / log(2) - 2;
gauss_pyr = build_gauss_pyr( init_img, octvs, intvls, sigma );
dog_pyr = build_dog_pyr( gauss_pyr, octvs, intvls );
storage = cvCreateMemStorage( 0 );
features = scale_space_extrema( dog_pyr, octvs, intvls, contr_thr,
curv_thr, storage );
calc_feature_scales( features, sigma, intvls );
if( img_dbl )
adjust_for_img_dbl( features );
calc_feature_oris( features, gauss_pyr );
compute_descriptors( features, gauss_pyr, descr_width, descr_hist_bins );
/* sort features by decreasing scale and move from CvSeq to array */
cvSeqSort( features, (CvCmpFunc)feature_cmp, NULL );
n = features->total;
*feat = calloc( n, sizeof(struct feature) );
*feat = cvCvtSeqToArray( features, *feat, CV_WHOLE_SEQ );
for( i = 0; i < n; i++ )
{
free( (*feat)[i].feature_data );
(*feat)[i].feature_data = NULL;
}
cvReleaseMemStorage( &storage );
cvReleaseImage( &init_img );
release_pyr( &gauss_pyr, octvs, intvls + 3 );
release_pyr( &dog_pyr, octvs, intvls + 2 );
return n;
}這個函式是上述函式一的過載，作用是一樣的，實際上函式一隻不過是使用預設引數呼叫了函式二，
核心的程式碼都是在函式二中實現的。
sift創始人David Lowe的完整程式碼，包括他的論文，請到此處下載：
http://www.cs.ubc.ca/~lowe/keypoints
日後，本BLOG內，會具體剖析下上述David Lowe的Sift演算法程式碼。
完。July、二零一一年二月十五日。