如果說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_{r\le i,\ c\le 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) = \begin{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\times9$的大小，並用它作為最小尺度空間值對影象進行濾波和斑點檢測。我們使用$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\times9$的模板對影象進行濾波，其結果作為最初始的尺度空間層（此時，尺度值為s=1.2，近似$\sigma=1.2$的高斯微分），後續的層將通過逐步放大濾波模板尺寸，以及放大後的模板不斷與影象進行濾波得到。由於採用盒子濾波和積分影象，濾波過程並不隨著濾波模板尺寸的增加而使運算工作量增加。

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

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

在通常尺度分析情況下，隨著尺度的增大，被檢測到的斑點數量迅速衰減。所以一般進行3-4組就可以了，與此同時，為了減少運算量，提高計算的速度，可以考慮在濾波時，將取樣間隔設為2。

對於尺寸為L的模板，當用它與積分圖運算來近似二維高斯核的濾波時，對應的二維高斯核的引數$\sigma = 1.2 \times \frac{L}{9}$，這一點至關重要，尤其是在後面計算描述子時，用於計算鄰域的半徑時。

4. 興趣點的定位

為了在影象及不同尺寸中定位興趣點，我們用了$3\times3\times3$鄰域非最大值抑制。具體的步驟基本與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\times5$、$3\times5$、$4\times5$的二維陣列表示，陣列的每一行表示一個黑白塊的位置引數。函式裡只初始化了第一層影象的模板引數，後面其他組其他層的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\n", point.pt.x, point.pt.y, point.size );*/
                        keypoints.push_back(kpt);
                    }
                }
            }
        }
    }
}

6. 總結

總體來說，如果理解了SIFT演算法，再來看SURF演算法會發現思路非常簡單。尤其是區域性最大值查詢方面，基本一致。關鍵還是一個用積分圖來簡化卷積的思路，以及怎麼用不同的模板來近似原來尺度空間中的高斯濾波器。

這一篇主要討論分析的是SURF的定位問題，下面還有SURF特徵點的方向計算與描述子的生成，將在下一篇文章中詳細描述。