SIFT特征點提取

本文轉載自查看原文 2019-10-16 11:09 1994 VS2015_opencv_C++提高

一、 SIFT算法

1、算法簡介

尺度不變特征轉換即SIFT (Scale-invariant feature transform)是一種計算機視覺的算法。它用來偵測與描述影像中的局部性特征，

它在空間尺度中尋找極值點，並提取出其位置、尺度、旋轉不變量，此算法由 David Lowe在1999年所發表，2004年完善總結。

其應用范圍包含物體辨識、機器人地圖感知與導航、影像縫合、3D模型建立、手勢辨識、影像追蹤和動作比對。

局部影像特征的描述與偵測可以幫助辨識物體，SIFT特征是基於物體上的一些局部外觀的興趣點而與影像的大小和旋轉無關。對於光線、噪聲、些微視角改變的容忍度也相當高。基於這些特性，它們是高度顯著而且相對容易擷取，在母數龐大的特征數據庫中，很容易辨識物體而且鮮有誤認。使用 SIFT特征描述對於部分物體遮蔽的偵測率也相當高，甚至只需要3個以上的SIFT物體特征就足以計算出位置與方位。在現今的電腦硬件速度下和小型的特征數據庫條件下，辨識速度可接近即時運算。SIFT特征的信息量大，適合在海量數據庫中快速准確匹配。

SIFT算法的實質是在不同的尺度空間上查找關鍵點(特征點)，並計算出關鍵點的方向。SIFT所查找到的關鍵點是一些十分突出，不會因光照，仿射變換和噪音等因素而變化的點，如角點、邊緣點、暗區的亮點及亮區的暗點等。

2、SIFT算法流程圖

二、SIFT算法操作步驟

1、圖像金字塔

1.1、高斯金字塔

圖像高斯金字塔(Gaussian Pyramid)是采用高斯函數對圖像進行模糊以及降采樣處理得到。其形成過程可如下圖所示：

其中高斯模糊系數計算公式如下：

1.1.1、高斯函數與圖像卷積

根據3σ原則，使用NxN的模板在圖像每一個像素點處操作，其中N=[（6σ+1）]且向上取最鄰近奇數。

其操作如下圖：

1.1.2、分離高斯卷積

上面這樣直接與圖像卷積，速度比較慢，同時圖像邊緣信息也會損失嚴重。后來，后來、、、，不知哪位學者發現，可以使用分離的高斯卷積（即先用1xN的模板沿着X方向對圖像卷積一次，然后用Nx1的模板沿着Y方向對圖像再卷積一次，其中N=[（6σ+1）]且向上取最鄰近奇數），這樣既省時也減小了直接卷積對圖像邊緣信息的嚴重損失。

1.1.3、高斯金子塔源碼分析

for (o = 0; o < octvs; o++)//金字塔組數為octvs，
    for (i = 0; i < intvls + 3; i++)//每一組有intvls + 3 層，intvls一般為3
        {
        if (o == 0 && i == 0)//如果是第一組第1層
            gauss_pyr[o][i] = cvCloneImage(base);//base 為原始灰度圖像經過升采樣或降采樣得到的圖像
    /* base of new octvave is halved image from end of previous octave */
        else if (i == 0)//建立非第一組的第1層
            gauss_pyr[o][i] = downsample(gauss_pyr[o - 1][intvls]);//降采樣圖像
        /* blur the current octave's last image to create the next one */
        else//建立非第一組的非第1層
            {
            gauss_pyr[o][i] = cvCreateImage(cvGetSize(gauss_pyr[o][i - 1]),IPL_DEPTH_32F, 1);
            cvSmooth(gauss_pyr[o][i - 1], gauss_pyr[o][i],CV_GAUSSIAN, 0, 0, sig[i], sig[i]);// sig[i]為模糊系數
            }//cvSmooth 為平滑處理函數，也即模糊處理。CV_GAUSSIAN 為選用高斯函數對圖像模糊

	return gauss_pyr;//返回建好的金字塔

1.2、高斯差分金字塔

2002年Mikolajczyk在詳細的實驗比較中發現尺度歸一化的高斯拉普拉斯函數的極大值和極小值同其它的特征提取函數，例如：梯度，Hessian或Harris角特征比較，能夠產生最穩定的圖像特征。而Lindeberg早在1994年就發現高斯差分函數（簡稱DOG算子）與尺度歸一化的高斯拉普拉斯函數非常近似。如下式：

其中k-1是個常數，並不影響極值點位置的求取。

1.2.1、差分金字塔的建立

差分金字塔的是在高斯金字塔的基礎上操作的，其建立過程是：在高斯金子塔中的每組中相鄰兩層相減（下一層減上一層）就生成高斯差分金字塔. 高斯差分金字塔其操作如下圖：

1.2.2、差分金字塔源碼分析

for (o = 0; o < octvs; o++)//octvs為高斯金字塔組數
    for (i = 0; i < intvls + 2; i++)//因為相減，故高斯金字塔中每組有（intvls + 2）層圖像
        {
        dog_pyr[o][i] = cvCreateImage(cvGetSize(gauss_pyr[o][i]),IPL_DEPTH_32F, 1);
        cvSub(gauss_pyr[o][i + 1], gauss_pyr[o][i], dog_pyr[o][i], NULL);//cvSub為opencv內置相減函數
        }
    return dog_pyr;//返回高斯差分金字塔

2、空間極值點（即關鍵點）檢測

關鍵點是由DOG空間的局部極值點組成的，關鍵點的初步探查是通過同一組內各DoG相鄰兩層圖像之間比較完成的。為了尋找DoG函數的極值點，每一個像素點要和它所有的相鄰點比較，看其是否比它的圖像域和尺度域的相鄰點大或者小。如圖下圖所示，中間的檢測點和它同尺度的8個相鄰點和上下相鄰尺度對應的9×2個點共26個點比較，以確保在尺度空間和二維圖像空間都檢測到極值點。

2.1、極值點檢測過程

2.1.1、極值點檢測示意

2.1.2、極值點檢測源碼分析

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))//pixval32f為提取圖像像素位置上的灰度值
            return 0;}
        else /* check for minimum */
        {
        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))//r c為圖像的行數和列數，dog_pyr為高斯差分圖
            return 0;

2.2、關鍵點定位

以上方法檢測到的極值點是離散空間的極值點，以下通過擬合三維二次函數來精確確定關鍵點的位置和尺度，同時去除低對比度的關鍵點和不穩定的邊緣響應點(因為DoG算子會產生較強的邊緣響應)，以增強匹配穩定性、提高抗噪聲能力。

2.2.1、關鍵點精確定位

離散空間的極值點並不是真正的極值點，下圖顯示了二維函數離散空間得到的極值點與連續空間極值點的差別。利用已知的離散空間點插值得到的連續空間極值點的方法叫做子像素插值。

為了提高關鍵點的穩定性，需要對尺度空間DoG函數進行曲線插值。利用DoG函數在尺度空間的Taylor展開式(插值函數)為：

上面算式的矩陣表示如下：

其中，X求導並讓方程等於零，可以得到極值點的偏移量為：

對應極值點，方程的值為：

其中, X^代表相對插值中心的偏移量，當它在任一維度上的偏移量大於0.5時（即x或y或 σ），意味着插值中心已經偏移到它的鄰近點上，所以必須改變當前關鍵點的位置。同時在新的位置上反復插值直到收斂；也有可能超出所設定的迭代次數或者超出圖像邊界的范圍，此時這樣的點應該刪除，在Lowe中進行了5次迭代。另外，過小的點易受噪聲的干擾而變得不穩定，所以將小於某個經驗值(Lowe論文中使用0.03，Rob Hess等人實現時使用0.04/S)的極值點刪除。同時，在此過程中獲取特征點的精確位置(原位置加上擬合的偏移量)以及尺度(σ)。

2.2.2、消除邊緣響應

一個定義不好的高斯差分算子的極值在橫跨邊緣的地方有較大的主曲率，而在垂直邊緣的方向有較小的主曲率。DOG算子會產生較強的邊緣響應，需要剔除不穩定的邊緣響應點。獲取特征點處的Hessian矩陣，主曲率通過一個2x2 的Hessian矩陣H求出（D的主曲率和H的特征值成正比）：

假設H的特征值為α和β（α、β代表x和y方向的梯度）且α>β。令α=rβ則有：

其中Tr(H)求取H的對角元素和；Det(H)為求H的行列式值。

則公式(r+1)^2/r的值在兩個特征值相等時最小，隨着的增大而增大。值越大，說明兩個特征值的比值越大，即在某一個方向的梯度值越大，而在另一個方向的梯度值越小，而邊緣恰恰就是這種情況。所以為了剔除邊緣響應點，需要讓該比值小於一定的閾值，因此，為了檢測主曲率是否在某域值r下，只需檢測:

論文建議r=10,OpenCv也采用r=10

2.2.3、精確定位中的泰勒插值源碼分析

while (i < SIFT_MAX_INTERP_STEPS)//SIFT_MAX_INTERP_STEPS=5為最大迭代次數，避免長時迭代
        {
        interp_step(dog_pyr, octv, intvl, r, c, &xi, &xr, &xc);// 泰勒展開擬合，xi,xr,xc依次為x、y、σ方向偏移量, 
        if (ABS(xi) < 0.5  &&  ABS(xr) < 0.5  &&  ABS(xc) < 0.5)//如果當前偏移量絕對值中的每個值均小於0.5，退出迭代
            break;
        c += cvRound(xc);//計算行坐標，cvRound 為四舍五入。
        r += cvRound(xr);
        intvl += cvRound(xi);
        if (intvl < 1 ||//不在計算的圖像層中
            intvl > intvls ||//高斯差分每組的層數為intvls 
            c < SIFT_IMG_BORDER ||//靠近圖像邊緣5個像素的區域不做檢測，SIFT_IMG_BORDER=5,
            r < SIFT_IMG_BORDER ||
            c >= dog_pyr[octv][0]->width - SIFT_IMG_BORDER ||//靠近圖像邊緣5個像素的區域不做檢測
            r >= dog_pyr[octv][0]->height - SIFT_IMG_BORDER)
            {
            return NULL;
            }
        i++;//迭代計數
        }

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 };
    dD = deriv_3D(dog_pyr, octv, intvl, r, c);//一階偏導數
    H = hessian_3D(dog_pyr, octv, intvl, r, c);//Hessian 矩陣即二階導數組成的矩陣
    H_inv = cvCreateMat(3, 3, CV_64FC1);
    cvInvert(H, H_inv, CV_SVD);//求Hessian矩陣的逆矩陣
    cvInitMatHeader(&X, 3, 1, CV_64FC1, x, CV_AUTOSTEP);
    cvGEMM(H_inv, dD, -1, NULL, 0, &X, 0); //cvGEMM為矩陣乘法，//第一個矩陣的系數；//H_inv、dD第一二個矩陣//-1矩陣前的常數//X結果矩陣
    cvReleaseMat(&dD); 
    cvReleaseMat(&H); 
    cvReleaseMat(&H_inv);
    *xi = x[2]; 
    *xr = x[1];
    *xc = x[0];
    }

3、關鍵點方向分配

為了使描述符具有旋轉不變性，需要利用圖像的局部特征為給每一個關鍵點分配一個基准方向。使用圖像梯度的方法求取局部結構的穩定方向。

3.1、特征點的梯度

3.1.1、梯度的計算

對於在DOG金字塔中檢測出的關鍵點點，采集其所在高斯金字塔圖像3σ領域窗口內像素的梯度和方向分布特征。梯度的模值和方向如下：

L為關鍵點所在的尺度空間值，按Lowe的建議，梯度的模值m(x,y)按 σ=1.5σ_oct 的高斯分布加成，按尺度采樣的3σ原則，領域窗口半徑為 3x1.5σ_oct。

3.1.1、梯度直方圖

在完成關鍵點的梯度計算后，使用直方圖統計領域內像素的梯度和方向。梯度直方圖將0~360度的方向范圍分為36個柱(bins)，其中每柱10度。如圖5.1所示，直方圖的峰值方向代表了關鍵點的主方向，(為簡化，圖中只畫了八個方向的直方圖)。

3.2、特征點主方向的確定

方向直方圖的峰值則代表了該特征點處鄰域梯度的方向，以直方圖中最大值作為該關鍵點的主方向。為了增強匹配的魯棒性，只保留峰值大於主方向峰值80％的方向作為該關鍵點的輔方向。因此，對於同一梯度值的多個峰值的關鍵點位置，在相同位置和尺度將會有多個關鍵點被創建但方向不同。僅有15％的關鍵點被賦予多個方向，但可以明顯的提高關鍵點匹配的穩定性。實際編程實現中，就是把該關鍵點復制成多份關鍵點，並將方向值分別賦給這些復制后的關鍵點，並且，離散的梯度方向直方圖要進行插值擬合處理，來求得更精確的方向角度值。

3.2.1、梯度圖像的平滑處理

為了防止某個梯度方向角度因受到噪聲的干擾而突變，我們還需要對梯度方向直方圖進行平滑處理。Opencv 所使用的平滑公式為：

其中i∈[0,35],h 和H 分別表示平滑前和平滑后的直方圖。由於角度是循環的，即0⁰=360⁰，如果出現h(j)，j超出了(0,…,35)的范圍，那么可以通過圓周循環的方法找到它所對應的、在0⁰=360⁰之間的值，如h(-1) = h(35)。

3.2.2、梯度直方圖拋物線插值

假設我們在第i個小柱子要找一個精確的方向，那么由上面分析知道：

設插值拋物線方程為h(t)=at²+bt+c,其中a、b、c為拋物線的系數，t為自變量，t∈[-1,1],此拋物線求導並令它等於0。

即h(t)´=0 得t_max=-b/(2a)

現在把這三個插值點帶入方程可得：

3.2.3、拋物線插值源碼分析

#define interp_hist_peak( l, c, r ) ( 0.5 * ((l)-(r)) / ((l) - 2.0*(c) + (r)) )//插值計算式，l為左側柱子值，r為左側柱子值
static void add_good_ori_features(CvSeq* features, double* hist, int n,
                                  double mag_thr, struct feature* feat)//精確主方向及輔方向
    {
    struct feature* new_feat;
    double bin, PI2 = CV_PI * 2.0;//CV_PI=pi
    int l, r, i;
    for (i = 0; i < n; i++)// 直方圖有n=36個小柱子
        {
        l = (i == 0) ? n - 1 : i - 1;//把小柱子看成是循環的，角度的取值為0-360即一個圓周
        r = (i + 1) % n;
        //只對小柱子的值大於等於主峰80%且此小柱子比左右兩邊小柱子都高的柱子進行拋物線插值
        if (hist[i] > hist[l] && hist[i] > hist[r] && hist[i] >= mag_thr)// mag_thr為>=80%的最高峰值
            {
            bin = i + interp_hist_peak(hist[l], hist[i], hist[r]);//interp_hist_peak 插值函數
            bin = (bin < 0) ? n + bin : (bin >= n) ? bin - n : bin;//角度取值約束在0-360之間，且是連續循環的
            new_feat = clone_feature(feat);//幅值特征點
            new_feat->ori = ((PI2 * bin) / n) - CV_PI;//？
            cvSeqPush(features, new_feat);
            free(new_feat);
            }

至此，圖像的關鍵點已檢測完畢，每個關鍵點有三個信息：位置、所處尺度、方向。由此可以確定一個SIFT特征區域。

4、特征點描述符
通過以上步驟，對於每一個關鍵點，擁有三個信息：位置、尺度以及方向。接下來就是為每個關鍵點建立一個描述符，使其不隨各種變化而改變，比如光照變化、視角變化等等。並且描述符應該有較高的獨特性，以便於提高特征點正確匹配的概率。

4.1、特征的生成過程
4.1.1、確定計算描述子所需的區域

將關鍵點附近的區域划分為d*d(Lowe建議d=4)個子區域，每個子區域作為一個種子點，每個種子點有8個方向。考慮到實際計算時，需要采用三線性插值，所需圖像窗口邊長為3x3xσ_oct x(d+1) 。在考慮到旋轉因素(方便下一步將坐標軸旋轉到關鍵點的方向)，如下圖6.1所示，實際計算所需的圖像區域半徑為：

4.1.2、坐標軸旋轉至主方向

將坐標軸旋轉為關鍵點的方向，以確保旋轉不變性。

4.1.3、梯度直方圖的生成

將鄰域內的采樣點分配到對應的子區域內，將子區域內的梯度值分配到8個方向上，計算其權值。

旋轉后的采樣點落在子區域的下標為

4.1.4、三線性插值

插值計算每個種子點八個方向的梯度。

采樣點在子區域中的下標(x'',y'') (圖中藍色窗口內紅色點)線性插值，計算其對每個種子點的貢獻。如圖中的紅色點，落在第0行和第1行之間，對這兩行都有貢獻。對第0行第3列種子點的貢獻因子為dr，對第1行第3列的貢獻因子為1-dr，同理，對鄰近兩列的貢獻因子為dc和1-dc，對鄰近兩個方向的貢獻因子為do和1-do。則最終累加在每個方向上的梯度大小為：

其中k，m，n為0。（像素點超出了對要插值區間的四個臨近子區間所在范圍）或為1（像素點處在對要插值區間的四個鄰近子區間之一所在范圍）。

4.1.5、特征描述子

如上統計的4*4*8=128個梯度信息即為該關鍵點的特征向量。

特征向量形成后，為了去除光照變化的影響，需要對它們進行歸一化處理，對於圖像灰度值整體漂移，圖像各點的梯度是鄰域像素相減得到，所以也能去除。得到的描述子向量為H=(h1,h2,.......,h128)，歸一化后的特征向量為L=(L1,L2,......,L128)，則

4.1.6、描述子的門限化

非線性光照，相機飽和度變化對造成某些方向的梯度值過大，而對方向的影響微弱。因此設置門限值(向量歸一化后，一般取0.2)截斷較大的梯度值(大於0.2的則就令它等於0.2，小於0.2的則保持不變)。然后再進行一次歸一化處理，提高特征的鑒別性。

4.2、描述子相關分析

用一組圖來概括描述子的生成過程

4.2.1、描述子生成總括

4.2.3、描述子三線性插值源碼分析

static void interp_hist_entry(double*** hist, double rbin, double cbin,
                                  double obin, double mag, int d, int n)
        {
        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(rbin);//向下取整
        c0 = cvFloor(cbin);
        o0 = cvFloor(obin);
        d_r = rbin - r0;//小數余項
        d_c = cbin - c0;
        d_o = obin - o0;
        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);//d_r = rbin - r0;
                row = hist[rb];
                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];
                        for (o = 0; o <= 1; o++)//沿着直方圖方向不超過1個單位角度 
                            { 
                            ob = (o0 + o) % n；//n=8，8個小柱子
                            v_o = v_c * ((o == 0) ? 1.0 - d_o : d_o);
                            h[ob] += v_o;
                            }
                        }
                    }
                }