本文大量參考了:http://blog.csdn.net/cxp2205455256/article/details/41311013;並且在此部落格的基礎之上，增加了一些新的內容。

一、SURF演算法的概述

       SURF，英文的全稱為Speed Up Robust Features,直譯為:加速版的具有魯棒特性的特徵演算法，是由Bay在2006年首次提出的。該演算法對經典的尺度不變特徵變換演算法(SIFT演算法)進行了改進，SIFT演算法最大的缺點就是如果不借助硬體或者專門的影象處理器進行加速的話，SIFT演算法很難達到實時處理的效果。SURF演算法最大的特徵在於採用了Haar特徵以及積分影象的概念，這大大的加速了程式的執行時間。SURF演算法不僅保持了SIFT演算法的尺度不變和旋轉不變的特性，而且對光照變化和放射變化同樣具有很強的魯棒性。SURF演算法一般應用於計算機視覺中的物體識別

二、積分影象
       積分影象的概念是由Viola和Jones提出的。積分影象中的任意一點（i，j）畫素的值為原影象左上角(影象原點)到任意點（i，j）相應的對焦區域的灰度值的總和，其數學公式如圖1所示:

       那麼，當我們想要計算影象中一個區域的積分，就只需計算這個區域的四個頂點在積分影象裡的值，便可以通過2步加法和2步減法計算得出，其數學公式如下：

三、SURF演算法的原理

    1、斑點檢測

       斑點:與周圍有著顏色和灰度差別的區域。在一個一維訊號中，讓它和高斯二階導數進行卷積，也就是拉普拉斯變換，那麼在訊號的邊緣處就會出現過零點的現象(其實就是:基於二階微分運算元---拉普拉斯運算元的邊緣檢測)。如下圖所示：

       高斯拉普拉斯Log邊緣檢測運算元的響應值就是在衡量影象的相似性，如下圖是一個影象的高斯拉普拉斯變換的三維圖和灰度圖的顯示，在影象中的斑點尺寸與高斯拉普拉斯函式的形狀趨於一致時，影象的拉普拉斯的響應值達到最大。

       Hession矩陣就是利用二階微分來進行斑點檢測，其矩陣如下：

2、構建Hessian矩陣構建高斯影象金字塔

      其實SURF構造的金字塔影象與SIFT構造的金字塔影象有很大的不同，正是因為這些不同才加快了SURF演算法關鍵點的檢測速度。SIFT採用的是DOG(高斯差分影象)影象，而SURF採用的是Hessian矩陣行列式近似值影象。Hessian矩陣式SURF演算法的核心。在數學中，Hessian矩陣是一個自變數為向量的實值函式的二階偏導陣列成的方塊矩陣。為了方便運算，假設函式f(x,y)，Hessian矩陣H是函式的二階偏導陣列成的方塊矩陣。首先來看看影象中某個畫素點的Hessian矩陣，如下圖所示:

即，我們在影象中的每一個畫素點，都可以求出一個Hessian矩陣，然後，我們再根據Hessian矩陣的判別式det(H)的正負號，來判別該點是不是極值點，Hessian矩陣的判別式如下所示:

       但是，由於特徵點需要具備尺度無關性，所在我們在Hessian矩陣的構造前，需要首先對影象進行高斯平滑濾波。這樣，在經過高斯平滑濾波之後，在對影象中的每一個畫素點進行Hessian矩陣的求解，具體的計算公式圖下所示:

                                                                  L(x,y) = G(t)*I(x,t)

       其中L(x,t)表示的是一幅影象在不同解析度下的表示，其本質表示的就是原始影象在不同尺度高斯核平滑之後的高斯平滑影象，可以利用高斯核G(t)和影象函式I(x)在點x的卷積來實現，其中的高斯核G(t)的計算公式為：高斯函式g(x)對變數x的二階導數。

       為了方便應用，Herbert Bay提出用近似值代替L(x,t)。為了平衡準確值與近似值之間的誤差引入了權值，權值隨著尺度的變化而變化，此時，Hessian矩陣的判別公式如下所示:

        由於求Hessian時，首先要先進行高斯平滑，然後求二階導數，這在離散的畫素點中應用模板的卷積就可以形成，那麼，我們就在想能不能用一個模板就可以將高斯平滑過程中的高斯平滑模板和在求影象二階導數過程中的拉普拉斯運算元的二階微分模板代替，答案是：這是可以的！然後，作者就提出了一種新的模板代替了這兩宗模板的組合。比如說Y方向上的模板就可以用下圖中的第二行的第二個模板代替；X方向上的模板就可以用第二行的第一個模板代替；X和Y方向上的二階混合偏導就可以用第一行中的第三個模板代替。其實，這就是Haar特徵中的:積分圖+邊緣模板、對角模板、中心環繞模板。至此，我們就用一個積分圖+模板的方式代替了高斯平滑和影象的二階求導的過程。具體的模板可以參考如下所示：

      上面講了這麼多，其實就只為了得到一張近似Hessian矩陣的行列式影象，SURF演算法中的Hessian矩陣行列式的影象就相當於SIFT演算法中的高斯差分影象，即DOG影象。但是，在金字塔影象中分為很多層，每一層叫一個Octave，每一個Octave中，又有幾張尺度不同的圖片(解析度不同，即高斯核函式的sigma不同)。在SIFT演算法中，同一個Octave層中的圖片尺寸相同，但是圖片的解析度(即圖片的模糊程度)不同。而不同的Octave層中的圖片尺寸大小也不相同，因為它是由上一層圖片降取樣得到的。在進行高斯模糊時，SIFT的高斯模板的大小是始終不變的，只是在不同的Octave之間，圖片的大小不一樣。而在SURF中，圖片的大小是一直不變的，不同Octave層待檢測圖片是改變了高斯模糊尺寸的大小得到的，當然了，同一個Octave中不同圖片用到的高斯模板的尺度(模糊程度)也是不同的。演算法允許尺度空間多層影象同時被處理，不需要對影象進行二次抽樣，從而提高了演算法的效能。當然在求解Hessian矩陣行列式影象的過程中，積分影象的使用，也是SURF演算法效能突飛猛進的原因之一。

四、利用非極大值抑制初步確定關鍵點的位置(特徵點)

   1、尺度金字塔構造

    在SURF中，採用不斷增大盒子濾波器模板尺寸與積分影象求取Hession矩陣響應，然後在響應影象上採用3D非極大值抑制，求取各種不同尺度

的斑點，以下是兩種不同的金字塔，SURF的金字塔屬於第二種：

    SURF中採用9X9尺寸的濾波器作為起始濾波器，之後的濾波器尺寸可由以下公式計算得出：

octave、interval在公式中都是從1開始，也就是當第0組第0層時，在公式中octave = 1， interval = 1。

    採用這種方式來定義濾波器尺寸的理由如下;

    濾波器響應長度、濾波器尺寸、組索引O、層索引S

尺度sigma之間的關係如下：

    採用類似的方法來處理其他幾組的模板序列。其方法是將濾波器尺寸增加量翻倍（6，12，24，38）。這樣，可以得到第二組的濾波器尺寸，它們分別為15，27，39，51。第三組的濾波器尺寸為27，51，75，99。如果原始影象的尺寸仍然大於對應的濾波器尺寸，尺度空間的分析還可以進行第四組，其對應的模板尺寸分別為51，99，147，195。下圖顯示了第一組至第三組的濾波器尺寸變化： 

      在通常尺度分析情況下，隨著尺度的增大，被檢測到的斑點數量迅速衰減。所以一般進行3-4組就可以了，與此同時，為了減少運算量，提高計算的速度，可以考慮在濾波時，將取樣間隔設為2。
    2、為了在影象及不同尺寸中定位興趣點，我們用了3×3×3鄰域非最大值抑制：

      此步驟和SIFT演算法類似，將經過Hessian矩陣處理過的每個畫素點與其三維領域的26個點的大小進行比較，如果它是這26個點中的最大值或者最小值，則保留下來，當做初步確定的關鍵點(特徵點)。檢測過程中，使用與該尺度影象解析度相對應大小的濾波器進行檢測。以3*3的濾波器為例，該尺度層影象中的9個畫素點之一的檢測特徵點與自身的8個特徵點和其上下兩個尺度層中的各9個點進行比較，共計與26個點比較。

    3、區域性極大值精確定位

      這裡也和SIFT演算法類似，採用三維線性插值演算法得到亞畫素級別的特徵點(關鍵點)，同時也去掉小於一定閾值的點，增加極值是檢測到的特徵點數量減少，最終只有幾個特徵最強的點被檢測了出來。

五、特徵點描述符

1、特徵點方向分配

    以特徵點為中心，計算半徑為6s(S為特徵點所在的尺度值)的鄰域內的點在x、y方向的Haar小波(Haar小波邊長取4s)響應，Harr小波

模板如圖所示：

    計算出影象在哈爾小波的x和y方向上的響應值之後，對兩個值進行因子為2S的高斯加權，加權後的值分別表示在水平和垂直方向上的方向分量。

Harr特徵值反應了影象灰度變化的情況，那麼這個主方向就是描述那些灰度變化特別劇烈的區域方向。

    接著，以特徵點為中心，張角為π/3的扇形滑動，計算視窗內的Harr小波響應值dx、dy的累加：

    扇形視窗的滑動如圖所示：

    在OpenSURF的C++程式碼實現如下：

1. for(int i = -6; i <= 6; ++i)   
2. {  
3.   for(int j = -6; j <= 6; ++j)   
4.   {  
5.     if(i*i + j*j < 36)   
6.     {  
7.       gauss = static_cast<float>(gauss25[id[i+6]][id[j+6]]);  
8.       resX[idx] = gauss * haarX(r+j*s, c+i*s, 4*s);  
9.       resY[idx] = gauss * haarY(r+j*s, c+i*s, 4*s);  
10.       Ang[idx] = getAngle(resX[idx], resY[idx]);  
11.       ++idx;  
12.     }  
13.   }  
14. }  

    通過i，j來控制以特徵點為中心的6X6的範圍，(i*i + j*j < 36)則篩選落在以特徵點為中心，半徑為6s的圓形區域內的點，然後計算HaarX與HarrY並與通過事先計算好的Gauss

濾波。並計算出每個點的角度。

    最後將最大值那個扇形的方向作為該特徵點的主方向。在OpenSURF尋找扇形中具有最大值得方向程式碼如下：

1. for(ang1 = 0; ang1 < 2*pi;  ang1+=0.15f) {  
2.    ang2 = ( ang1+pi/3.0f > 2*pi ? ang1-5.0f*pi/3.0f : ang1+pi/3.0f);  
3.    sumX = sumY = 0.f;   
4.    for(unsigned int k = 0; k < Ang.size(); ++k)   
5.    {  
6.      // get angle from the x-axis of the sample point
7.      constfloat & ang = Ang[k];  
8.      // determine whether the point is within the window
9.      if (ang1 < ang2 && ang1 < ang && ang < ang2)   
10.      {  
11.        sumX+=resX[k];    
12.        sumY+=resY[k];  
13.      }   
14.      elseif (ang2 < ang1 &&   
15.        ((ang > 0 && ang < ang2) || (ang > ang1 && ang < 2*pi) ))   
16.      {  
17.        sumX+=resX[k];    
18.        sumY+=resY[k];  
19.      }  
20.    }  
21.    // if the vector produced from this window is longer than all 
22.    // previous vectors then this forms the new dominant direction
23.    if (sumX*sumX + sumY*sumY > max)   
24.    {  
25.      // store largest orientation
26.      max = sumX*sumX + sumY*sumY;  
27.      orientation = getAngle(sumX, sumY);  
28.    }  
29.  }  

六、構造SURF特徵點的特徵描述子

       在SIFT演算法中，是在特徵點周圍取16*16的鄰域，並且把該鄰域化為4*4的小區域，每個小區域統計8個方向的梯度，最後得到4*4*8=128維的向量，並且把該向量作為SIFT特徵點的特徵描述子。

       在SURF演算法中，也是在特徵點周圍取一個正方形框，框的邊長為20s(s是所檢測到該特徵點所在的尺度)。該框帶方向，方向當然就是上面一步檢測出來的主方向。然後把該框分為16個子區域，每個子區域統計25個畫素的水平方向和垂直方向的haar小波特徵，這裡的水平方向和垂直方向是相對於主方向而言的。該haar小波特徵為水平方向值之和、水平方向絕對值之和和垂直方向之和、垂直方向絕對值之和。這樣每個小區域就有4個值，所以每個特徵點就是16*4維的向量，相比SIFT而言，少了一半，這在特徵匹配的過程中會大大的加速匹配的速度。

    以特徵點為中心，沿主方向將20SX20S的影象劃分為4X4個子塊，每個子塊用尺寸2S的Harr模板進行響應值計算，並統計每個子塊中

在計算這個矩形區域時並不是先把它旋轉到主方向，而是先計算出每一個點的Harr響應值dx、dy並高斯加權處理後，把dx、dy進行旋轉變換，計算

公式如下：

    在OpenSURF的實現原始碼中採用的是另外一種方式，通過點旋轉公式，把點旋轉到主方向上並進行最近鄰插值的對應點，公式如下：

七、匹配

     為了加速匹配過程，SURF藉助Laplacian（在之前計算Hessian是可以順便得出，不佔用太多的時間）的符號使匹配過程索引加快。這樣可以將下面的情況區分開，然後在進行描述符匹配：