1.6 SURF特征提取與匹配

如果說SIFT算法中使用DoG對LoG進行了簡化，提高了搜索特征點的速度，那么，Bay等人提出的加速健壯特征（Speeded Up Robust Features, SURF）算法，則是對DoH的簡化和近似。雖然Lowe的SIFT算法已被認為是最有效的，也是最常用的特征點提取算法，但如果不借助于硬件的加速和專用圖形處理器的配合，SIFT算法在現有計算機性能的條件下仍然很難達到實時的程度。對于需要實時運算的場合，如導彈尋的制導時的目標識別，通常需要在毫秒級內完成特征點搜索、特征矢量生成、特征矢量匹配、目標鎖定等工作，這樣SIFT算法就很難適應這種需求了。SURF借鑒了SIFT中簡化近似的思想，將DoH中的高斯二階微分模板進行了近似簡化，使得模板對圖像的濾波只需要進行幾個簡單的加減法運算，并且，這種運算與濾波模板的尺寸無關。實驗證明，SURF算法較SIFT在運算速度上要快3倍左右，綜合性能要優于SIFT算法。

1.6.1 積分圖像

SURF算法中要用到積分圖像的概念。借助積分圖像，圖像與高斯二階微分模板的濾波轉化為對積分圖像的加減運算。積分圖像（Integral Image）的概念是由Viola和Jones提出來的，而將類似積分圖像用于盒子濾波（Box Filter）卻是Simard等人提出來的。

積分圖像中任意一點（i,j）的值ii（i,j），為原圖像左上角到任意點（i,j）相應的對角線區域灰度值的總和，即：式中，S（i,j）表示一列的積分，且S（i, －1）＝0,ii（－1,j）＝0。求積分圖像，只需對原圖像所有像素進行一遍掃描。

式中，p（i′,j′）表示原圖像中點（i′,j′）的灰度值，ii（i,j）可用下面兩式迭代計算得到：

如圖1.6.1所示，在求取窗口W內的像元灰度和時，不管窗口W 的大小如何，均可以用積分圖像的4個相應點（i1,j1）、（i2,j2）、（i3,j3）、（i4,j4）的值計算得到。也就是說，求取窗口W 內的像元灰度和與窗口的尺寸是無關的。窗口W 內像元的灰度和為：

假設有一灰度值均為1的圖像，那么該圖像中任一點（i,j）的積分圖像值實際就是圖像左上角點到任一點（i,j）構成的矩形區域面積（像元數）的大小。就是由W構成的矩形框包含的面積。

1.6.2 DoH近似

給定圖像I中一個點x（x,y），在點x處，尺度為σ的Hessian矩陣H（x,σ）定義如下：

式中，Lxx（x,σ）是高斯二階微分在點x處與圖像I的卷積，Lxy（x,σ）和Lyy（x,σ）具有同樣的含義。

由于二階高斯微分模板被離散化和剪裁的原因，導致了圖像在旋轉奇數倍的π/4時，即轉動到模板的對角線方向時，特征點檢測的重復性（Repeatability）降低；而在π/2時，特征點檢測的重復性最高。但這一小小的不足不影響我們使用Hessian矩陣進行特征點檢測。

為了將模板與圖像的卷積轉化成盒子濾波（Box Filter）運算，需要對高斯二階微分模板進行簡化，使得簡化后的模板只是由幾個矩形區域組成，矩形區域內填充同一值，如圖1.6.2所示，在簡化模板中白色區域的值為1，黑色區域的值為－1，灰色區域的值為0。

對于σ＝1.2的高斯二階微分濾波，設定模板的尺寸為9×9，并用它作為最小尺度空間值對圖像進行濾波和斑點檢測。使用Dxx、Dyy和Dxy表示模板與圖像進行卷積的結果，這樣，便可將Hessian矩陣的行列式作如下簡化：

式中，Y＝＝0.912≈0.9。理論上講，對于不同σ值和對應的模板尺寸，Y值是不同的，但為了簡化起見，可以認為它是一個常數。同樣，也可以認為C為常數，由于常數C不影響極大值求取，因此，便有

不過，在實際計算濾波響應值時需要將模板盒子尺寸（面積）進行歸一化處理，以保證使用一個統一的Frobenius范數，能適應所有的濾波尺寸。如對于9×9模板的Lx x和L yy盒子的面積為15,Lxy的盒子面積為9。一般而言，如果盒子內部填充值為vn∈｛,－1,1－｝2，盒子對應的四個角點積分圖像值為，盒子面積分別為sxx、syy（sxx＝syy）和sxy，那么，盒子濾波響應值為：

從Dxx、Dyy和Dxy的計算公式可以看出，它們的運算量與模板的尺寸是無關的。計算Dxx和Dyy只有12次加減法和4次乘法，計算Dxy只有16次加減法和5次乘法。

使用近似的Hessian矩陣行列式來表示圖像中某一點x處的斑點響應值，遍歷圖像中所有的像元點，便形成了在某一尺度下斑點檢測的響應圖像。使用不同的模板尺寸，便形成了多尺度斑點響應的金字塔圖像，利用這一金字塔圖像，就可以進行斑點響應極值點的搜索，其過程完全與SIFT算法類同。

1.6.3 尺度空間表示

通常要想獲取不同尺度的斑點，必須建立圖像的尺度空間金字塔。一般的方法是通過采用不同σ的高斯函數，對圖像進行平滑濾波，然后，重采樣圖像以獲得更高一層的金字塔圖像。Lowe在SIFT方法中就是通過相鄰兩層圖像金字塔相減得到DoG圖像，然后再在DoG圖像上進行斑點和邊緣檢測工作。

由于采用了盒子濾波和積分圖像，所以，并不需要像SIFT算法那樣去直接建立圖像金字塔，而是采用不斷增大盒子濾波模板尺寸的間接方法。通過不同尺寸盒子濾波模板與積分圖像求取Hessian矩陣行列式的響應圖像，然后，在響應圖像上采用3D非最大值抑制，求取各種不同尺度的斑點。

如前所述，使用9×9的模板對圖像進行濾波，其結果作為最初始的尺度空間層（此時，尺度值s＝1.2，近似σ＝1.2的高斯微分），后續的層將通過逐步放大濾波模板尺寸，以及放大后的模板不斷與圖像進行濾波得到。由于采用盒子濾波和積分圖像，濾波過程并不隨著濾波模板尺寸的增加而使運算工作量增加。

與SIFT算法類似，需要將尺度空間劃分成若干組（Octaves）。一個組代表了逐步放大的濾波模板對同一輸入圖像進行濾波的一系列響應圖，每個組又由若干固定的層組成。由于積分圖像離散化的原因，兩個層之間的最小尺度變化量是由高斯二階微分濾波器在微分方向上對正負斑點響應長度l0決定的，它是盒子濾波模板尺寸的1/3。對于9×9的盒子濾波模板而言，l0為3。下一個層的響應長度至少應該在l0的基礎上增加2個像元，以保證一邊一個像元，即：l0＝5，這樣，模板的尺寸就為15×15，如圖1.6.3所示。以此類推，可以得到一個尺寸逐漸增大的模板序列，它們的尺寸分別為：9×9、15× 15、21×21、27×27、39×39，黑色、白色區域的長度增加偶數個像元，以保證一個中心像元的存在。

采用類似的方法來處理其他組的模板序列。其方法是將濾波器尺寸增加量翻倍（6、12、24、48）。這樣，可以得到第二組的濾波器尺寸，它們分別為15、27、39、51；第三組的濾波器尺寸為27、51、75、99。如果原始圖像尺寸仍然大于對應的濾波器尺寸，尺度空間的分析還可進行第四組，其對應的模板尺寸分別是51、99、147和195。圖1.6.4給出了第一到第三組的濾波器尺寸變化的圖形表示。對數水平軸代表尺度、組之間有相互重疊，其目的是為了覆蓋所有可能的尺度。在通常尺度分析情況下，隨著尺度的增大，被檢測到的斑點數量迅速衰減。與此同時，為了減少運算量，提高計算的速度，可以考慮在濾波時，將采樣間隔設為2°。

仿照SIFT算法，可以給出濾波響應長度l、濾波器的尺寸L、組索引ο、層索引s、尺度σ之間的相互關系：

濾波器可以采用矢量數據結構來表示。數據結構中分別包含濾波器中每個盒子的坐標、盒子中的填充值、盒子的面積等信息。坐標既可以濾波器中心像元為原點，也可以左上角像元為原點。

圖1.6.5所示為不同尺寸時的濾波器模板的圖形化表示。

1.6.4 SURF特征描述算子

SIFT特征描述算子在生成特征矢量時使用的是高斯圖像，而SURF特征描述子在生成特征矢量時使用的則是積分圖像。這樣做的目的就是要充分利用在特征點檢測時形成的中間結果（積分圖像），避免在特征矢量生成時對圖像進行重復運算。

為了保證特征矢量具有旋轉不變性，與SIFT算法一樣，需要對每個特征點分配一個主方向。為此，需要在以特征點為中心，以6s（s為特征點的尺度）為半徑的圓形區域內，對圖像進行Haar小波響應運算。

在SIFT特征描述算子中，在求取特征點主方向時，是以特征點為中心，在以1.6σ為半徑的鄰域內計算梯度方向直方圖的。事實上，兩種方法在求取特征點主方向時，考慮到Haar小波的模板帶寬，實際計算梯度的圖像區域是相同的，如圖1.6.6所示。

與SIFT算法類似，使用σ＝2s的高斯加權函數對Haar小波的響應值進行高斯加權。為了求取主方向，需要設計一個以特征點為中心，張角為的扇形滑動窗口。如圖1.6.7所示，以步長約0.2弧度旋轉這個滑動窗口，并對滑動窗口內圖像Haar小波響應值dx、dy進行累加，得到一個矢量（mw,θw）：

主方向為最大Haar響應累加值所對應的方向，也就是最長矢量所對應的方向，即：

可以仿照SIFT算法求主方向時的策略，當存在另一個相當于主峰值80%能量的峰值時，則將這個方向認為是該特征點的輔方向。一個特征點可能會被指定具有多個方向（一個主方向，一個以上輔方向），這可以增強匹配的魯棒性。和SIFT的描述子類似，如果當在mw中出現另一個大于主峰能量max｛mw｝80 %時的次峰，可以將該特征點復制成兩個特征點：一個主方向為最大響應能量所對應的方向，另一個的主方向為次最大響應能量所對應的方向。

生成特征點描述算子與確定特征點的方向有些類似，它需要計算圖像的Haar小波響應。不過，與主方向的確定不同的是，這次不是使用一個圓形區域，而是在一個矩形區域來計算Haar小波響應。以特征點為中心，沿特征點的主方向將20×20的圖像劃分成4×4個子塊，每個子塊利用尺寸2s的Haar小波模板進行響應值計算，然后對響應值進行統計∑dx、∑|dx|、∑dy、∑|dy| 形成特征矢量，如圖1.6.8所示。圖中，以特征點為中心，以20s為邊長的矩形窗口為特征描述算子計算使用的窗口，特征點到矩形邊框的線段表示特征點的主方向。

將20s的窗口劃分成4×4子窗口，每個子窗口中有5×5個像元，使用尺度為2s的Haar小波對子窗口圖像進行其響應值計算，共進行25次采樣，分別得到沿主方向的dy和垂直于主方向的dx。然后，以特征點為中心，對dx、dy進行高斯加權計算，其中σ＝3.3s。最后，分別對每個子塊的響應值進行統計，得到每個子塊的矢量：

由于共有4×4個子塊，因此，特征描述算子共由4×4×4＝64維特征矢量組成。SURF描述算子不僅具有尺度和旋轉不變性，而且對光照的變化也具有不變性。使用小波響應本身就具有亮度不變性，而對比度的不變性則是通過將特征矢量進行歸一化來實現。圖1.6.9給出了三種不同圖像模式的子塊得到的不同結果。對于實際圖像描述算子，可以認為它們是由這三種不同模式圖像的描述子組合而成的。

為了充分利用積分圖像進行Haar小波的響應計算，不直接通過旋轉Haar小波模板求得其響應值，而是在積分圖像上先使用水平和垂直的Haar小波模板求其響應值，在求得響應值dx和dy后，然后根據主方向旋轉dx和dy，使其與主方向保持一致。為了求旋轉后的Haar小波響應值，首先要得到旋轉前圖像的位置。旋轉前后圖像的位置關系，可以通過點的旋轉公式得到：

在得到點（j,i）在旋轉前對應積分圖像的位置（x,y）后，利用積分圖像與水平、垂直Haar小波求得水平和垂直兩個方向的響應值dx和dy。對dx和dy進行高斯加權處理，并根據主方向的角度，對dx和dy進行旋轉變換，從而，得到旋轉后的dx′和dy′。其計算公式如下：

圖1.6.10說明在有噪聲干擾下，SURF特征描述算子與沒有噪聲干擾時具有相同的特征矢量。

一般而言，特征矢量的長度越長，特征矢量所承載的信息量就越大，特征描述子的獨特性就越好，但匹配所付出的時間代價就越大。對于SURF描述算子，可以將它擴展用到128維矢量來表示。具體的方法是在求∑dx、∑|dx|時，區分dy＜0和dy≥0情況。同樣，在求取∑dy、∑|dy|時，區分dx＜0和dx≥0情況。這樣，每個子塊就產生了8個梯度統計值，從而使描述子特征矢量的長度增加到8×4×4＝128維。

為了實現快速匹配，SURF算法在特征矢量中增加了一個新的變量，即特征點的拉普拉斯響應正負號。在特征點檢測時，將Hessian矩陣的跡（Trace）的正負號記錄下來，作為特征矢量的一個變量。這樣做并不增加運算量，因為特征點檢測時已經對Hessian矩陣的跡進行計算了。在特征匹配時，這個變量可以有效地節省搜索時間，因為只有兩個正負號相同的特征點才可能匹配，對于不同正負號的特征點就不再進行相似性計算了。簡單地說，可以根據特征點的響應值符號，將特征點分成兩組，一組是具有拉普拉斯正響應的特征點，另一組是具有拉普拉斯負響應的特征點。匹配時，只有符號相同組中的特征點才能進行相互匹配。

1.6.5 程序實現

在最新的MATLAB 8.0中，可以調用計算機視覺系統工具箱中（Computer Vision System Toolbox）的detectSURFFeatures（）函數來檢測輸入灰度圖像的SURF特征，其具體使用方法如下：

POINTS＝detectSURFFeatures（I, Name, Value）

功能：用于檢測灰度圖像的SURF特征。

輸入：I——待檢測的灰度圖像；

Name, Value——MetricThreshold：矩陣閾值，只有大于該閾值的才能確定

為SURF特征點，默認值為1000；

——NumOctaves：組數（Octaves），默認值為3；

——NumScaleLevels：每組的尺度層數，默認值為4。

輸出：POINTS——返回的一個對象，其中包含著SURF特征點的信息。

當獲得一幅圖像的SURF特征點的信息后，可以調用extractFeatures（）函數獲取SURF特征向量。

extractFeatures（）函數的具體使用方法如下：

［FEATURES, VALID＿POINTS］＝extractFeatures（I, POINTS）

功能：提取特征點的特征向量。

輸入：I——灰度圖像；

POINTS——特征點信息。

輸出：FEATURES——特征描述向量；

VALID＿POINTS——有效特征點坐標。

【例1.6.1】 基于SURF尺度不變特征點在復雜環境下進行目標探測的MATLAB程序。

    ***********************************************************
    %  步驟 1：讀入圖像
    boxImage ＝ imread(＇stapleRemover.jpg＇);
    figure; imshow(boxImage);
    title(＇Image of a Box＇);
    sceneImage ＝ imread(＇clutteredDesk.jpg＇);
    figure; imshow(sceneImage);
    title(＇Image of a Cluttered Scene＇);
    %  步驟2：檢測SURF特征點
    boxPoints ＝ detectSURFFeatures(boxImage);
    scenePoints ＝ detectSURFFeatures(sceneImage);
    figure; imshow(boxImage);
    title(＇100 Strongest Feature Points from Box Image＇);
    hold on;
    plot(boxPoints.selectStrongest(100));
    figure; imshow(sceneImage);
    title(＇300 Strongest Feature Points from Scene Image＇);
    hold on;
    plot(scenePoints.selectStrongest(300));
    %  步驟3：提取特征向量
    [boxFeatures, boxPoints]＝ extractFeatures(boxImage, boxPoints);
    [sceneFeatures, scenePoints]＝ extractFeatures(sceneImage, scenePoints);
    %  步驟4：進行特征點粗匹配
    boxPairs ＝ matchFeatures(boxFeatures, sceneFeatures);
    matchedBoxPoints ＝ boxPoints(boxPairs(:,1), :);
    matchedScenePoints ＝ scenePoints(boxPairs(:,2), :);
    figure;
    showMatchedFeatures(boxImage, sceneImage, matchedBoxPoints, ...
        matchedScenePoints, ＇montage＇);
    title(＇PutativelyMatched Points(Including Outliers)＇);
    %  步驟5：根據匹配的特征點確定目標在場景圖像中的位置
    [tform, inlierBoxPoints, inlierScenePoints]＝ ...
        estimateGeometricTransform(matchedBoxPoints, matchedScenePoints, ＇affine＇);
    figure;
    showMatchedFeatures(boxImage, sceneImage, inlierBoxPoints, ...
        inlierScenePoints, ＇montage＇);
    %  步驟6：將檢測到的物體在場景中框出
    title(＇Matched Points(Inliers Only)＇);
    boxPolygon ＝[1,1; ...                             % top－left
            size(boxImage,2),1; ...                    % top－right
            size(boxImage,2), size(boxImage,1); ...    % bottom－right
            1, size(boxImage,1); ...                   % bottom－left
            1,1];                                      % top－left again to close the polygon
    newBoxPolygon ＝ transformPointsForward(tform, boxPolygon);
    figure; imshow(sceneImage);
    hold on;
    line(newBoxPolygon(:,1), newBoxPolygon(:,2), ＇Color＇, ＇y＇);
    title(＇Detected Box＇);
    %  檢測另外一個目標
    elephantImage ＝ imread(＇elephant.jpg＇);
    figure; imshow(elephantImage);
    title(＇Image of an Elephant＇);
    elephantPoints ＝ detectSURFFeatures(elephantImage);
    figure; imshow(elephantImage);
    hold on;
    plot(elephantPoints.selectStrongest(100));
    title(＇100 Strongest Feature Points from Elephant Image＇);
    [elephantFeatures, elephantPoints]＝ extractFeatures(elephantImage, elephantPoints);
    elephantPairs ＝ matchFeatures(elephantFeatures, sceneFeatures, ＇MaxRatio＇,0.9);
    [tform, inlierElephantPoints, inlierScenePoints]＝ ...
        estimateGeometricTransform(matchedElephantPoints, matchedScenePoints, ＇affine＇);
    figure;
    showMatchedFeatures(elephantImage, sceneImage, inlierElephantPoints, ...
        inlierScenePoints, ＇montage＇);
    title(＇Matched Points(Inliers Only)＇);
    elephantPolygon ＝[1,1; ...                                 % top－left
            size(elephantImage,2),1; ...                        % top－right
            size(elephantImage,2), size(elephantImage,1); ...   % bottom－right
            1, size(elephantImage,1); ...                       % bottom－left
            1,1];                                               % top－left again to close the polygon
    newElephantPolygon ＝ transformPointsForward(tform, elephantPolygon);
    figure;
    imshow(sceneImage);
    hold on;
    line(newBoxPolygon(:,1), newBoxPolygon(:,2), ＇Color＇, ＇y＇);
    line(newElephantPolygon(:,1), newElephantPolygon(:,2), ＇Color＇, ＇g＇);
    title(＇Detected Elephant and Box＇);
    ***********************************************************

圖1.6.11所示為例1.6.1的運行結果。