九之再續：教你一步一步用c語言實現sift算法、上

                     教你一步一步用c語言實現sift算法、上

作者：July、二零一一年三月十二日
出處：http://blog.csdn.net/v_JULY_v
參考：Rob Hess維護的sift 庫
環境：windows xp+vc6.0
條件：c語言實現。
說明：本BLOG內會陸續一一實現所有經典算法。
------------------------

引言：
    在我寫的關于sift算法的前倆篇文章里頭，已經對sift算法有了初步的介紹：九、圖像特征提取與匹配之SIFT算法，而后在：九（續）、sift算法的編譯與實現里，我也簡單記錄下了如何利用opencv，gsl等庫編譯運行sift程序。
    但據一朋友表示，是否能用c語言實現sift算法，同時，盡量不用到opencv，gsl等第三方庫之類的東西。而且，Rob Hess維護的sift 庫，也不好懂，有的人根本搞不懂是怎么一回事。
    那么本文，就教你如何利用c語言一步一步實現sift算法，同時，你也就能真正明白sift算法到底是怎么一回事了。

    ok，先看一下，本程序最終運行的效果圖，sift 算法分為五個步驟（下文詳述），對應以下第二--第六幅圖：

sift算法的步驟
    要實現一個算法，首先要完全理解這個算法的原理或思想。咱們先來簡單了解下，什么叫sift算法：
    sift，尺度不變特征轉換，是一種電腦視覺的算法用來偵測與描述影像中的局部性特征，它在空間尺度中尋找極值點，并提取出其位置、尺度、旋轉不變量，此算法由 David Lowe 在1999年所發表，2004年完善總結。

    所謂，Sift算法就是用不同尺度（標準差）的高斯函數對圖像進行平滑，然后比較平滑后圖像的差別，
差別大的像素就是特征明顯的點。

    以下是sift算法的五個步驟：
    一、建立圖像尺度空間(或高斯金字塔)，并檢測極值點

    首先建立尺度空間，要使得圖像具有尺度空間不變形，就要建立尺度空間，sift算法采用了高斯函數來建立尺度空間，高斯函數公式為：
        G(x,y,e) = [1/2*pi*e^2] * exp[ -(x^2 + y^2)/2e^2]
    上述公式G(x,y,e)，即為尺度可變高斯函數。

    而，一個圖像的尺度空間L（x,y,e) ,定義為原始圖像I(x,y)與上述的一個可變尺度的2維高斯函數G(x,y,e) 卷積運算。
    即，原始影像I(x，y)在不同的尺度e下，與高斯函數G(x,y,e)進行卷積，得到L(x,y,e)，如下：
        L(x,y,e) = G(x,y,e)*I(x,y)

    以上的（x，y）是空間坐標， e，是尺度坐標，或尺度空間因子，e的大小決定平滑程度，大尺度對應圖像的概貌特征，小尺度對應圖像的細節特征。大的e值對應粗糙尺度(低分辨率)，反之，對應精細尺度(高分辨率)。

    尺度，受e這個參數控制的表示。而不同的L（x,y,e）就構成了尺度空間，具體計算的時候，即使連續的高斯函數，都被離散為（一般為奇數大小）(2*k+1) *(2*k+1)矩陣，來和數字圖像進行卷積運算。
    隨著e的變化，建立起不同的尺度空間，或稱之為建立起圖像的高斯金字塔。

    但，像上述L(x,y,e) = G(x,y,e)*I(x,y)的操作，在進行高斯卷積時，整個圖像就要遍歷所有的像素進行卷積(邊界點除外)，于此，就造成了時間和空間上的很大浪費。
    為了更有效的在尺度空間檢測到穩定的關鍵點，也為了縮小時間和空間復雜度，對上述的操作作了一個改建：即，提出了高斯差分尺度空間（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算子的近似。

    ok，耐心點，咱們再來總結一下上述內容：
    1、高斯卷積
    在組建一組尺度空間后，再組建下一組尺度空間，對上一組尺度空間的最后一幅圖像進行二分之一采樣，得到下一組尺度空間的第一幅圖像，然后進行像建立第一組尺度空間那樣的操作，得到第二組尺度空間，公式定義為
         L(x,y,e) = G(x,y,e)*I(x,y)

    圖像金字塔的構建：圖像金字塔共O組，每組有S層，下一組的圖像由上一組圖像降采樣得到，效果圖，圖A如下(左為上一組，右為下一組)：

   
    2、高斯差分
    在尺度空間建立完畢后，為了能夠找到穩定的關鍵點，采用高斯差分的方法來檢測那些在局部位置的極值點，即采用倆個相鄰的尺度中的圖像相減，即公式定義為：
        D(x,y,e) = ((G(x,y,ke) - G(x,y,e)) * I(x,y)
                 = L(x,y,ke) - L(x,y,e)

    效果圖，圖B：

    SIFT的精妙之處在于采用圖像金字塔的方法解決這一問題，我們可以把兩幅圖像想象成是連續的，分別以它們作為底面作四棱錐，就像金字塔，那么每一個 截面與原圖像相似，那么兩個金字塔中必然會有包含大小一致的物體的無窮個截面，但應用只能是離散的，所以我們只能構造有限層，層數越多當然越好，但處理時 間會相應增加，層數太少不行，因為向下采樣的截面中可能找不到尺寸大小一致的兩個物體的圖像。

    咱們再來具體闡述下構造D(x,y,e)的詳細步驟：
    1、首先采用不同尺度因子的高斯核對圖像進行卷積以得到圖像的不同尺度空間，將這一組圖像作為金子塔圖像的第一層。
    2、接著對第一層圖像中的2倍尺度圖像（相對于該層第一幅圖像的2倍尺度）以2倍像素距離進行下采樣來得到金子塔圖像的第二層中的第一幅圖像，對該圖像采用不同尺度因子的高斯核進行卷積，以獲得金字塔圖像中第二層的一組圖像。
    3、再以金字塔圖像中第二層中的2倍尺度圖像（相對于該層第一幅圖像的2倍尺度）以2倍像素距離進行下采樣來得到金字塔圖像的第三層中的第一幅圖像，對該圖像采用不同尺度因子的高斯核進行卷積，以獲得金字塔圖像中第三層的一組圖像。這樣依次類推，從而獲得了金字塔圖像的每一層中的一組圖像，如下圖所示：

    4、對上圖得到的每一層相鄰的高斯圖像相減，就得到了高斯差分圖像，如下述第一幅圖所示。下述第二幅圖中的右列顯示了將每組中相鄰圖像相減所生成的高斯差分圖像的結果，限于篇幅，圖中只給出了第一層和第二層高斯差分圖像的計算（下述倆幅圖統稱為圖2）：

    5、因為高斯差分函數是歸一化的高斯拉普拉斯函數的近似，所以可以從高斯差分金字塔分層結構提取出圖像中的極值點作為候選的特征點。對DOG 尺度空間每個點與相鄰尺度和相鄰位置的點逐個進行比較，得到的局部極值位置即為特征點所處的位置和對應的尺度。

    二、檢測關鍵點
    為了尋找尺度空間的極值點，每一個采樣點要和它所有的相鄰點比較，看其是否比它的圖像域和尺度域的相鄰點大或者小。如下圖，圖3所示，中間的檢測點和它同尺度的8個相鄰點和上下相鄰尺度對應的9×2個點共26個點比較，以確保在尺度空間和二維圖像空間都檢測到極值點。

    因為需要同相鄰尺度進行比較，所以在一組高斯差分圖像中只能檢測到兩個尺度的極值點（如上述第二幅圖中右圖的五角星標識），而其它尺度的極值點檢測則需要在圖像金字塔的上一層高斯差分圖像中進行。依次類推，最終在圖像金字塔中不同層的高斯差分圖像中完成不同尺度極值的檢測。
    當然這樣產生的極值點并不都是穩定的特征點，因為某些極值點響應較弱，而且DOG算子會產生較強的邊緣響應。

    三、關鍵點方向的分配
    為了使描述符具有旋轉不變性，需要利用圖像的局部特征為給每一個關鍵點分配一個方向。利用關鍵點鄰域像素的梯度及方向分布的特性，可以得到梯度模值和方向如下：

   
    其中，尺度為每個關鍵點各自所在的尺度。
    在以關鍵點為中心的鄰域窗口內采樣，并用直方圖統計鄰域像素的梯度方向。梯度直方圖的范圍是0～360度，其中每10度一個方向，總共36個方向。
直方圖的峰值則代表了該關鍵點處鄰域梯度的主方向，即作為該關鍵點的方向。

    在計算方向直方圖時，需要用一個參數等于關鍵點所在尺度1.5倍的高斯權重窗對方向直方圖進行加權，上圖中用藍色的圓形表示，中心處的藍色較重，表示權值最大，邊緣處顏色潛，表示權值小。如下圖所示，該示例中為了簡化給出了8方向的方向直方圖計算結果，實際sift創始人David Lowe的原論文中采用36方向的直方圖。

    方向直方圖的峰值則代表了該特征點處鄰域梯度的方向，以直方圖中最大值作為該關鍵點的主方向。為了增強匹配的魯棒性，只保留峰值大于主方向峰值80％的方向作為該關鍵點的輔方向。因此，對于同一梯度值的多個峰值的關鍵點位置，在相同位置和尺度將會有多個關鍵點被創建但方向不同。僅有15％的關鍵點被賦予多個方向，但可以明顯的提高關鍵點匹配的穩定性。
       
    至此，圖像的關鍵點已檢測完畢，每個關鍵點有三個信息：位置、所處尺度、方向。由此可以確定一個SIFT特征區域。

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

    接下來以關鍵點為中心取8×8的窗口。
    上圖，圖5中左部分的中央黑點為當前關鍵點的位置，每個小格代表關鍵點鄰域所在尺度空間的一個像素，箭頭方向代表該像素的梯度方向，箭頭長度代表梯度模值，圖中藍色的圈代表高斯加權的范圍（越靠近關鍵點的像素梯度方向信息貢獻越大）。
    然后在每4×4的小塊上計算8個方向的梯度方向直方圖，繪制每個梯度方向的累加值，即可形成一個種子點，如圖5右部分所示。此圖中一個關鍵點由2×2共4個種子點組成，每個種子點有8個方向向量信息。這種鄰域方向性信息聯合的思想增強了算法抗噪聲的能力，同時對于含有定位誤差的特征匹配也提供了較好的容錯性。

    實際計算過程中，為了增強匹配的穩健性，Lowe建議對每個關鍵點使用4×4共16個種子點來描述，這樣對于一個關鍵點就可以產生128個數據，即最終形成128維的SIFT特征向量。此時SIFT特征向量已經去除了尺度變化、旋轉等幾何變形因素的影響，再繼續將特征向量的長度歸一化，則可以進一步去除光照變化的影響。

    五、最后一步：當兩幅圖像的SIFT特征向量生成后，下一步我們采用關鍵點特征向量的歐式距離來作為兩幅圖像中關鍵點的相似性判定度量。取上圖中，圖像A中的某個關鍵點，并找出其與圖像B中歐式距離最近的前兩個關鍵點，在這兩個關鍵點中，如果最近的距離除以次近的距離少于某個比例閾值，則接受這一對匹配點。降低這個比例閾值，SIFT匹配點數目會減少，但更加穩定。

sift算法的逐步c實現
    ok，上文攪了那么多的理論，如果你沒有看懂它，咋辦列?沒關系，下面，咱們來一步一步實現此sift算法，即使你沒有看到上述的理論，慢慢的，你也會明白sift算法到底是怎么一回事，sift算法到底是怎么實現的...。
   yeah，請看：

前期工作：
在具體編寫核心函數之前，得先做幾個前期的準備工作：

0、頭文件：

3、聲明幾個圖像的基本處理函數：
    圖像處理基本函數，其實也可以用OPENCV的函數代替，但本文，咱們選擇了用c語言實現，盡量不用第三方庫的東西，所以，還得自己編寫這些函數：

主體實現
    ok，以上所有的工作都就緒以后，那么接下來，咱們就先來編寫main函數，因為你一看主函數之后，你就立馬能發現sift算法的工作流程及其原理了。
    （主函數中涉及到的函數，下一篇文章：一、教你一步一步用c語言實現sift算法、下，咱們自會一個一個編寫）：