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的运行结果。