SURF算法与源码分析、上

如果说SIFT算法中使用DOG对LOG进行了简化，提高了搜索特征点的速度，那么SURF算法则是对DoH的简化与近似。虽然SIFT算法已经被认为是最有效的，也是最常用的特征点提取的算法，但如果不借助于硬件的加速和专用图像处理器的配合，SIFT算法以现有的计算机仍然很难达到实时的程度。对于需要实时运算的场合，如基于特征点匹配的实时目标跟踪系统，每秒要处理8-24帧的图像，需要在毫秒级内完成特征点的搜索、特征矢量生成、特征矢量匹配、目标锁定等工作，这样SIFT算法就很难适应这种需求了。SURF借鉴了SIFT中简化近似的思想，把DoH中的高斯二阶微分模板进行了简化，使得模板对图像的滤波只需要进行几个简单的加减法运算，并且，这种运算与滤波器的尺度无关。实验证明，SURF算法较SIFT在运算速度上要快3倍左右。

1. 积分图像

SURF算法中要用到积分图像的概念。借助积分图像，图像与高斯二阶微分模板的滤波转化为对积分图像的加减运算。

积分图像中任意一点$(i,j)$的值$ii(i,j)$，为原图像左上角到点$(i,j)$相应的对角线区域灰度值的总和，即

$$ii(i,j) = sum_{rle i, cle j}p(r,c)$$

式中，$p(r,c)$表示图像中点$(r,c)$的灰度值，$ii(i,j)$可以用下面两式迭代计算得到

$$S(i,j) = S(i,j-1)+p(i,j)$$

$$ii(i,j)=ii(i-1,j)+S(i,j)$$

式中，$S(i,j)$表示一列的积分，且$S(i,-1)=0,ii(-1,j)=0$。求积分图像，只需要对原图像所有像素进行一遍扫描。

OpenCV中提供了用于计算积分图像的接口

/*
* src :输入图像，大小为M*N
* sum: 输出的积分图像，大小为（M+1）*（N+1）
* sdepth：用于指定sum的类型，-1表示与src类型一致
*/
void integral(InputArray src, OutputArray sum, int sdepth = -1);

值得注意的是OpenCV里的积分图大小比原图像多一行一列，那是因为OpenCV中积分图的计算公式为：

$$ii(i,j) = sum_{r< i, c< j}p(r,c)$$

一旦积分图计算好了，计算图像内任何矩形区域的像素值的和只需要三个加法，如上图所示。

2. DoH近似

在斑点检测这篇文章中已经提到过，我们可以利用Hessian矩阵行列式的极大值检测斑点。下面我们给出Hessian矩阵的定义。

给定图像$I$中的一个点$x(i,j)$，在点$x$处，尺度为$sigma$的Hessian矩阵$H(x,sigma)$定义如下：

$$H(x,sigma) = egin{bmatrix}L_{xx}(x,sigma) & L_{xy}(x,sigma) \ L_{xy}(x,sigma) & L_{yy}(x,sigma) end{bmatrix}$$

式中，$L_{xx}(x,sigma)$是高斯二阶微分$frac{partial ^2g(sigma)}{partial x^2}$在点$x$处与图像$I$的卷积，$L_{x,y}(x,sigma)$和$L_{yy}(x,sigma)$具有类似的含义。

下面显示的是上面三种高斯微分算子的图形。

但是利用Hessian行列式进行图像斑点检测时，有一个缺点。由于二阶高斯微分被离散化和裁剪的原因，导致了图像在旋转奇数倍的$pi/4$时，即转换到模板的对角线方向时，特征点检测的重复性降低（也就是说，原来特征点的地方，可能检测不到特征点了）。而在$pi/2$时，特征点检测的重现率真最高。但这一小小的不足不影响我们使用Hessian矩阵进行特征点的检测。

为了将模板与图产像的卷积转换为盒子滤波运算，我们需要对高斯二阶微分模板进行简化，使得简化后的模板只是由几个矩形区域组成，矩形区域内填充同一值，如下图所示，在简化模板中白色区域的值为正数，黑色区域的值为负数，灰度区域的值为0。

对于$sigma = 1.2$的高斯二阶微分滤波器，我们设定模板的尺寸为$9 imes9$的大小，并用它作为最小尺度空间值对图像进行滤波和斑点检测。我们使用$D_{xx}、D_{xy}$和$D_{yy}$表示模板与图像进行卷积的结果。这样，便可以将Hessian矩阵的行列式作如下的简化。

$$Det(H) = L_{xx}L_{yy} – L_{xy}L_{xy} = D_{xx}frac{L_{xx}}{D_{xx}}D_{yy}frac{L_{yy}}{D_{yy}} - D_{xy}frac{L_{xy}}{D_{xy}}D_{xy}frac{L_{xy}}{D_{xy}}  approx D_{xx}D_{yy} – (wD_{xy})^2$$

滤波器响应的相关权重$w$是为了平衡Hessian行列式的表示式。这是为了保持高斯核与近似高斯核的一致性。

$$w = frac{|L_{xy}(sigma)|_F|D_{xx}(sigma)_F|}{|L_{xx}(sigma)|_F|D_{xy}(sigma)_F|} = 0.912$$

其中$|X|_F$为Frobenius范数。理论上来说对于不同的$sigma$的值和对应尺寸的模板尺寸，$w$值是不同的，但为了简化起见，可以认为它是同一个常数。

使用近似的Hessian矩阵行列式来表示图像中某一点$x$处的斑点响应值，遍历图像中所有的像元点，便形成了在某一尺度下琉璃点检测的响应图像。使用不同的模板尺寸，便形成了多尺度斑点响应的金字塔图像，利用这一金字塔图像，就可以进行斑点响应极值点的搜索，其过程完全与SIFT算法类同。

3. 尺度空间表示

通常想要获取不同尺度的斑点，必须建立图像的尺度空间金字塔。一般的方法是通过不同$sigma$的高斯函数，对图像进行平滑滤波，然后重采样图像以获得更高一层的金字塔图像。SIFT特征检测算法中就是通过相邻两层图像金字塔相减得到DoG图像，然后再在DoG图像上进行斑点和边缘检测工作的。

由于采用了盒子滤波和积分图像，所以，我们并不需要像SIFT算法那样去直接建立图像金字塔，而是采用不断增大盒子滤波模板的尺寸的间接方法。通过不同尺寸盒子滤波模板与积分图像求取Hessian矩阵行列式的响应图像。然后在响应图像上采用3D非最大值抑制，求取各种不同尺度的斑点。

如前所述，我们使用$9 imes9$的模板对图像进行滤波，其结果作为最初始的尺度空间层（此时，尺度值为s=1.2，近似$sigma=1.2$的高斯微分），后续的层将通过逐步放大滤波模板尺寸，以及放大后的模板不断与图像进行滤波得到。由于采用盒子滤波和积分图像，滤波过程并不随着滤波模板尺寸的增加而使运算工作量增加。

与SIFT算法类似，我们需要将尺度空间划分为若干组（Octaves）。一个组代表了逐步放大的滤波模板对同一输入图像进行滤波的一系列响应图。每个组又由若干固定的层组成。由于积分图像离散化的原因，两个层之间的最小尺度变化量是由高斯二阶微分滤波器在微分方向上对正负斑点响应长度$l_0$决定的，它是盒子滤波器模板尺寸的$1/3$。对于$9 imes9$的模板，它的$l_0=3$。一下层的响应长度至少应该在$l_0$的基础上增加2个像素，以保证一边一个像素，即$l_0 = 5$。这样模板的尺寸就为$15 imes15$。以此类推，我们可以得到一个尺寸增大模板序列，它们的尺寸分别为：$9 imes9，15 imes15，21 imes21，27 imes27$，黑色、白色区域的长度增加偶数个像素，以保证一个中心像素的存在。

       

采用类似的方法来处理其他几组的模板序列。其方法是将滤波器尺寸增加量翻倍（6，12，24，38）。这样，可以得到第二组的滤波器尺寸，它们分别为15，27，39，51。第三组的滤波器尺寸为27，51，75，99。如果原始图像的尺寸仍然大于对应的滤波器尺寸，尺度空间的分析还可以进行第四组，其对应的模板尺寸分别为51，99，147和195。下图显示了第一组至第三组的滤波器尺寸变化。

在通常尺度分析情况下，随着尺度的增大，被检测到的斑点数量迅速衰减。所以一般进行3-4组就可以了，与此同时，为了减少运算量，提高计算的速度，可以考虑在滤波时，将采样间隔设为2。

对于尺寸为L的模板，当用它与积分图运算来近似二维高斯核的滤波时，对应的二维高斯核的参数$sigma = 1.2 imes frac{L}{9}$，这一点至关重要，尤其是在后面计算描述子时，用于计算邻域的半径时。

4. 兴趣点的定位

为了在图像及不同尺寸中定位兴趣点，我们用了$3 imes3 imes3$邻域非最大值抑制。具体的步骤基本与SIFT一致，而且Hessian矩阵行列式的最大值在尺度和图像空间被插值。

下面显示了我们用的快速Hessian检测子检测到的兴趣点。

5. SURF源码解析

这份源码来自OpenCV nonfree模块。

这里先介绍SURF特征点定位这一块，关于特征点的描述下一篇文章再介绍。

5.1 主干函数 fastHessianDetector

特征点定位的主干函数为fastHessianDetector，该函数接受一个积分图像，以及尺寸相关的参数，组数与每组的层数，检测到的特征点保存在vector<KeyPoint>类型的结构中。

static void fastHessianDetector(const Mat& sum, const Mat& msum, vector<KeyPoint>& keypoints,
    int nOctaves, int nOctaveLayers, float hessianThreshold)
{
    /*first Octave图像采样的步长，第二组的时候加倍，以此内推
    增加这个值，将会加快特征点检测的速度，但是会让特征点的提取变得不稳定*/
    const int SAMPLE_STEP0 = 1;

    int nTotalLayers = (nOctaveLayers + 2)*nOctaves; // 尺度空间的总图像数
    int nMiddleLayers = nOctaveLayers*nOctaves; // 用于检测特征点的层的 总数，也就是中间层的总数

    vector<Mat> dets(nTotalLayers); // 每一层图像 对应的 Hessian行列式的值
    vector<Mat> traces(nTotalLayers); // 每一层图像 对应的 Hessian矩阵的迹的值
    vector<int> sizes(nTotalLayers); // 每一层用的 Harr模板的大小
    vector<int> sampleSteps(nTotalLayers); // 每一层用的采样步长 
    vector<int> middleIndices(nMiddleLayers); // 中间层的索引值

    keypoints.clear();

    // 为上面的对象分配空间，并赋予合适的值
    int index = 0, middleIndex = 0, step = SAMPLE_STEP0;

    for (int octave = 0; octave < nOctaves; octave++)
    {
        for (int layer = 0; layer < nOctaveLayers + 2; layer++)
        {
            /*这里sum.rows - 1是因为 sum是积分图，它的大小是原图像大小加1*/
            dets[index].create((sum.rows - 1) / step, (sum.cols - 1) / step, CV_32F); // 这里面有除以遍历图像用的步长
            traces[index].create((sum.rows - 1) / step, (sum.cols - 1) / step, CV_32F);
            sizes[index] = (SURF_HAAR_SIZE0 + SURF_HAAR_SIZE_INC*layer) << octave;
            sampleSteps[index] = step;

            if (0 < layer && layer <= nOctaveLayers)
                middleIndices[middleIndex++] = index;
            index++;
        }
        step *= 2;
    }
    // Calculate hessian determinant and trace samples in each layer
    for (int i = 0; i < nTotalLayers; i++)
    {
        calcLayerDetAndTrace(sum, sizes[i], sampleSteps[i], dets[i], traces[i]);
    }

    // Find maxima in the determinant of the hessian
    for (int i = 0; i < nMiddleLayers; i++)
    {
        int layer = middleIndices[i];
        int octave = i / nOctaveLayers;
        findMaximaInLayer(sum, msum, dets, traces, sizes, keypoints, octave, layer, hessianThreshold, sampleSteps[layer]);
    }

    std::sort(keypoints.begin(), keypoints.end(), KeypointGreater());
}

5.2 计算Hessian矩阵的行列式与迹calcLayerDetAndTrace

这个函数首先定义了尺寸为9的第一层图像的三个模板。模板分别为一个$3 imes5$、$3 imes5$、$4 imes5$的二维数组表示，数组的每一行表示一个黑白块的位置参数。函数里只初始化了第一层图像的模板参数，后面其他组其他层的Harr模板参数都是用resizeHaarPattern这个函数来计算的。这个函数返回的是一个SurfHF的结构体，这个结构体由两个点及一个权重构成。

struct SurfHF
{
    int p0, p1, p2, p3;
    float w;

    SurfHF() : p0(0), p1(0), p2(0), p3(0), w(0) {}
};

resizeHaarPattern这个函数非常的巧妙，它把模板中的点坐标。转换到在积分图中的相对（模板左上角点）坐标。

static void
resizeHaarPattern(const int src[][5], SurfHF* dst, int n, int oldSize, int newSize, int widthStep)
{
    float ratio = (float)newSize / oldSize;
    for (int k = 0; k < n; k++)
    {
        int dx1 = cvRound(ratio*src[k][0]);
        int dy1 = cvRound(ratio*src[k][1]);
        int dx2 = cvRound(ratio*src[k][2]);
        int dy2 = cvRound(ratio*src[k][3]);
        /*巧妙的坐标转换*/
        dst[k].p0 = dy1*widthStep + dx1; // 转换为一个相对距离，距离模板左上角点的  在积分图中的距离 !!important!!
        dst[k].p1 = dy2*widthStep + dx1; 
        dst[k].p2 = dy1*widthStep + dx2;
        dst[k].p3 = dy2*widthStep + dx2;
        dst[k].w = src[k][4] / ((float)(dx2 - dx1)*(dy2 - dy1));// 原来的+1,+2用 覆盖的所有像素点平均。
    }
}

在用积分图计算近似卷积时，用的是calcHaarPattern函数。这个函数比较简单，只用知道左上与右下角坐标即可。

inline float calcHaarPattern(const int* origin, const SurfHF* f, int n)
{
    /*orgin即为积分图，n为模板中 黑白 块的个数 */
    double d = 0;
    for (int k = 0; k < n; k++)
        d += (origin[f[k].p0] + origin[f[k].p3] - origin[f[k].p1] - origin[f[k].p2])*f[k].w;
    return (float)d;
}

最终我们可以看到了整个calcLayerDetAndTrack的代码

static void calcLayerDetAndTrace(const Mat& sum, int size, int sampleStep,
    Mat& det, Mat& trace)
{
    const int NX = 3, NY = 3, NXY = 4;
    const int dx_s[NX][5] = { { 0, 2, 3, 7, 1 }, { 3, 2, 6, 7, -2 }, { 6, 2, 9, 7, 1 } };
    const int dy_s[NY][5] = { { 2, 0, 7, 3, 1 }, { 2, 3, 7, 6, -2 }, { 2, 6, 7, 9, 1 } };
    const int dxy_s[NXY][5] = { { 1, 1, 4, 4, 1 }, { 5, 1, 8, 4, -1 }, { 1, 5, 4, 8, -1 }, { 5, 5, 8, 8, 1 } };

    SurfHF Dx[NX], Dy[NY], Dxy[NXY];

    if (size > sum.rows - 1 || size > sum.cols - 1)
        return;
    resizeHaarPattern(dx_s, Dx, NX, 9, size, sum.cols);
    resizeHaarPattern(dy_s, Dy, NY, 9, size, sum.cols);
    resizeHaarPattern(dxy_s, Dxy, NXY, 9, size, sum.cols);

    /* The integral image 'sum' is one pixel bigger than the source image */
    int samples_i = 1 + (sum.rows - 1 - size) / sampleStep; // 最大能遍历到的 行坐标，因为要减掉一个模板的尺寸
    int samples_j = 1 + (sum.cols - 1 - size) / sampleStep; // 最大能遍历到的 列坐标

    /* Ignore pixels where some of the kernel is outside the image */
    int margin = (size / 2) / sampleStep;

    for (int i = 0; i < samples_i; i++)
    {
        /*坐标为(i,j)的点是模板左上角的点，所以实际现在模板分析是的i+margin,j+margin点处的响应*/
        const int* sum_ptr = sum.ptr<int>(i*sampleStep);
        float* det_ptr = &det.at<float>(i + margin, margin); // 左边空隙为 margin
        float* trace_ptr = &trace.at<float>(i + margin, margin);
        for (int j = 0; j < samples_j; j++)
        {
            float dx = calcHaarPattern(sum_ptr, Dx, 3);
            float dy = calcHaarPattern(sum_ptr, Dy, 3);
            float dxy = calcHaarPattern(sum_ptr, Dxy, 4);
            sum_ptr += sampleStep;
            det_ptr[j] = dx*dy - 0.81f*dxy*dxy;
            trace_ptr[j] = dx + dy;
        }
    }
}

5.3 局部最大值搜索findMaximaInLayer

这里算法思路很简单，值得注意的是里面的一些坐标的转换很巧妙，里面比较重的函数就是interpolateKeypoint函数，通过插值计算最大值点。

/*
* Maxima location interpolation as described in "Invariant Features from
* Interest Point Groups" by Matthew Brown and David Lowe. This is performed by
* fitting a 3D quadratic to a set of neighbouring samples.
*
* The gradient vector and Hessian matrix at the initial keypoint location are
* approximated using central differences. The linear system Ax = b is then
* solved, where A is the Hessian, b is the negative gradient, and x is the
* offset of the interpolated maxima coordinates from the initial estimate.
* This is equivalent to an iteration of Netwon's optimisation algorithm.
*
* N9 contains the samples in the 3x3x3 neighbourhood of the maxima
* dx is the sampling step in x
* dy is the sampling step in y
* ds is the sampling step in size
* point contains the keypoint coordinates and scale to be modified
*
* Return value is 1 if interpolation was successful, 0 on failure.
*/

static int
interpolateKeypoint(float N9[3][9], int dx, int dy, int ds, KeyPoint& kpt)
{
    Vec3f b(-(N9[1][5] - N9[1][3]) / 2,  // Negative 1st deriv with respect to x
        -(N9[1][7] - N9[1][1]) / 2,  // Negative 1st deriv with respect to y
        -(N9[2][4] - N9[0][4]) / 2); // Negative 1st deriv with respect to s

    Matx33f A(
        N9[1][3] - 2 * N9[1][4] + N9[1][5],            // 2nd deriv x, x
        (N9[1][8] - N9[1][6] - N9[1][2] + N9[1][0]) / 4, // 2nd deriv x, y
        (N9[2][5] - N9[2][3] - N9[0][5] + N9[0][3]) / 4, // 2nd deriv x, s
        (N9[1][8] - N9[1][6] - N9[1][2] + N9[1][0]) / 4, // 2nd deriv x, y
        N9[1][1] - 2 * N9[1][4] + N9[1][7],            // 2nd deriv y, y
        (N9[2][7] - N9[2][1] - N9[0][7] + N9[0][1]) / 4, // 2nd deriv y, s
        (N9[2][5] - N9[2][3] - N9[0][5] + N9[0][3]) / 4, // 2nd deriv x, s
        (N9[2][7] - N9[2][1] - N9[0][7] + N9[0][1]) / 4, // 2nd deriv y, s
        N9[0][4] - 2 * N9[1][4] + N9[2][4]);           // 2nd deriv s, s

    Vec3f x = A.solve(b, DECOMP_LU);

    bool ok = (x[0] != 0 || x[1] != 0 || x[2] != 0) &&
        std::abs(x[0]) <= 1 && std::abs(x[1]) <= 1 && std::abs(x[2]) <= 1;

    if (ok)
    {
        kpt.pt.x += x[0] * dx;
        kpt.pt.y += x[1] * dy;
        kpt.size = (float)cvRound(kpt.size + x[2] * ds);
    }
    return ok;
}

static void findMaximaInLayer(const Mat& sum, const Mat& mask_sum,
    const vector<Mat>& dets, const vector<Mat>& traces,
    const vector<int>& sizes, vector<KeyPoint>& keypoints,
    int octave, int layer, float hessianThreshold, int sampleStep)
{
    // Wavelet Data
    const int NM = 1;
    const int dm[NM][5] = { { 0, 0, 9, 9, 1 } };
    SurfHF Dm;

    int size = sizes[layer];

    // 当前层图像的大小
    int layer_rows = (sum.rows - 1) / sampleStep;
    int layer_cols = (sum.cols - 1) / sampleStep;

    // 边界区域大小，考虑的下一层的模板大小
    int margin = (sizes[layer + 1] / 2) / sampleStep + 1;

    if (!mask_sum.empty())
        resizeHaarPattern(dm, &Dm, NM, 9, size, mask_sum.cols);

    int step = (int)(dets[layer].step / dets[layer].elemSize());

    for (int i = margin; i < layer_rows - margin; i++)
    {
        const float* det_ptr = dets[layer].ptr<float>(i);
        const float* trace_ptr = traces[layer].ptr<float>(i);
        for (int j = margin; j < layer_cols - margin; j++)
        {
            float val0 = det_ptr[j]; // 中心点的值
            if (val0 > hessianThreshold)
            {
                // 模板左上角的坐标
                int sum_i = sampleStep*(i - (size / 2) / sampleStep);
                int sum_j = sampleStep*(j - (size / 2) / sampleStep);

                /* The 3x3x3 neighbouring samples around the maxima.
                The maxima is included at N9[1][4] */

                const float *det1 = &dets[layer - 1].at<float>(i, j);
                const float *det2 = &dets[layer].at<float>(i, j);
                const float *det3 = &dets[layer + 1].at<float>(i, j);
                float N9[3][9] = { { det1[-step - 1], det1[-step], det1[-step + 1],
                    det1[-1], det1[0], det1[1],
                    det1[step - 1], det1[step], det1[step + 1] },
                    { det2[-step - 1], det2[-step], det2[-step + 1],
                    det2[-1], det2[0], det2[1],
                    det2[step - 1], det2[step], det2[step + 1] },
                    { det3[-step - 1], det3[-step], det3[-step + 1],
                    det3[-1], det3[0], det3[1],
                    det3[step - 1], det3[step], det3[step + 1] } };

                /* Check the mask - why not just check the mask at the center of the wavelet? */
                if (!mask_sum.empty())
                {
                    const int* mask_ptr = &mask_sum.at<int>(sum_i, sum_j);
                    float mval = calcHaarPattern(mask_ptr, &Dm, 1);
                    if (mval < 0.5)
                        continue;
                }

                /* 检测val0，是否在N9里极大值，？？为什么不检测极小值呢*/
                if (val0 > N9[0][0] && val0 > N9[0][1] && val0 > N9[0][2] &&
                    val0 > N9[0][3] && val0 > N9[0][4] && val0 > N9[0][5] &&
                    val0 > N9[0][6] && val0 > N9[0][7] && val0 > N9[0][8] &&
                    val0 > N9[1][0] && val0 > N9[1][1] && val0 > N9[1][2] &&
                    val0 > N9[1][3] && val0 > N9[1][5] &&
                    val0 > N9[1][6] && val0 > N9[1][7] && val0 > N9[1][8] &&
                    val0 > N9[2][0] && val0 > N9[2][1] && val0 > N9[2][2] &&
                    val0 > N9[2][3] && val0 > N9[2][4] && val0 > N9[2][5] &&
                    val0 > N9[2][6] && val0 > N9[2][7] && val0 > N9[2][8])
                {
                    /* Calculate the wavelet center coordinates for the maxima */
                    float center_i = sum_i + (size - 1)*0.5f;
                    float center_j = sum_j + (size - 1)*0.5f;

                    KeyPoint kpt(center_j, center_i, (float)sizes[layer],
                        -1, val0, octave, CV_SIGN(trace_ptr[j]));

                    /* 局部极大值插值，用Hessian，类似于SIFT里的插值，里面没有迭代5次，只进行了一次查找，why?  */
                    int ds = size - sizes[layer - 1];
                    int interp_ok = interpolateKeypoint(N9, sampleStep, sampleStep, ds, kpt);

                    /* Sometimes the interpolation step gives a negative size etc. */
                    if (interp_ok)
                    {
                        /*printf( "KeyPoint %f %f %d
", point.pt.x, point.pt.y, point.size );*/
                        keypoints.push_back(kpt);
                    }
                }
            }
        }
    }
}

6. 总结

总体来说，如果理解了SIFT算法，再来看SURF算法会发现思路非常简单。尤其是局部最大值查找方面，基本一致。关键还是一个用积分图来简化卷积的思路，以及怎么用不同的模板来近似原来尺度空间中的高斯滤波器。

这一篇主要讨论分析的是SURF的定位问题，下面还有SURF特征点的方向计算与描述子的生成，将在下一篇文章中详细描述。