上一篇文章 SURF演算法與原始碼分析、上 中主要分析的是SURF特徵點定位的演算法原理與相關OpenCV中的原始碼分析，這篇文章接著上篇文章對已經定位到的SURF特徵點進行特徵描述。這一步至關重要，這是SURF特徵點匹配的基礎。總體來說演算法思路和SIFT相似，只是每一步都做了不同程度的近似與簡化，提高了效率。

1. SURF特徵點方向分配

為了保證特徵向量具有旋轉不變性，與SIFT特徵一樣，需要對每個特徵點分配一個主方向。為些，我們需要以特徵點為中心，以$6s$（$s = 1.2 *L /9$為特徵點的尺度）為半徑的圓形區域，對影象進行Haar小波響應運算。這樣做實際就是對影象進行梯度運算只不過是我們需要利用積分影象，提高計算影象梯度的效率。在SIFT特徵描述子中我們在求取特徵點主方向時，以是特徵點為中心，在以4.5$\sigma$為半徑的鄰域內計算梯度方向直方圖。事實上，兩種方法在求取特徵點主方向時，考慮到Haar小波的模板頻寬，實際計算梯度的影象區域是相同的。用於計算梯度的Harr小波的尺度為4s。

與SIFT類似，使用$\sigma = 2s$的高斯加權函式對Harr小波的響應值進行高斯加權。為了求取主方向值，需要設計一個以特徵點為中心，張角為$\pi/3$的扇形滑動視窗。以步長為0.2弧度左右，旋轉這個滑動視窗，並對滑動視窗內的影象Harr小波響應值dx、dy進行累加，得到一個向量$(m_w,\theta_w)$：

$$m_w = \sum _w dx + \sum_w dy$$

$$\theta_w = arctan(\sum_w dx / \sum_w dy)$$

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

$$\theta = \theta_w|max\{m_w\}$$

可以依照SIFT求方方向時策略，當存在另一個相當於主峰值80%能量的峰值時，則將這個方向認為是該特徵點的輔方向。一個特徵點可能會被指定具有多個方向（一個主方向，一個以上輔方向），這可以增強匹配的魯棒性。和SIFT的描述子類似，如果在$m_w$中出現另一個大於主峰能量$max\{m_w\} 80%$時的次峰，可以將該特徵點複製成兩個特徵點。一個主的方向為最大響應能量所對應的方向，另一個主方向為次大響應能量所對應的方向。

圖 1  求取主方向時扇形滑動視窗圍繞特徵點轉動，統計Haar小波響應值，並計算方向角

2. 特徵點特徵向量生成

生成特徵點描述子與確定特徵點方向有些類似，它需要計算影象的Haar小波響應。不過，與主方向的確定不同的是，這次我們不是使用一個圓形區域，而是在一個矩形區域來計算Haar小波響應。以特徵點為中心，沿上一節討論得到的主方向，沿主方向將$s20s\times20s$的影象劃分為$4\times4$個子塊，每個子塊利用尺寸$2s$的Harr模板進行響應值進行響應值計算，然後對響應值進行統計$\sum dx$、$\sum |dx|$、$\sum dy$、$\sum |dy|$形成特徵向量。如下圖2所示。圖中，以特徵點為中心，以20s為邊長的矩形視窗為特徵描述子計算使用的視窗，特徵點到矩形邊框的線段表示特徵點的主方向。

圖2 特徵描述子表示

將$20s$的視窗劃分成$4\times4$子視窗，每個子視窗有$5s\times5s$個畫素。使用尺寸為$2s$的Harr小波對子視窗影象進行其響應值計算，共進行25次取樣，分別得到沿主方向的dy和垂直於主方向的dx。然後，以特徵點為中心，對dy和dx進行高斯加權計算，高斯核的引數為$\sigma = 3.3s (即20s/6)$。最後，分別對每個子塊的響應值進行統計，得到每個子塊的向量：

$$V_{子塊}=\left[\sum dx,\sum |dx|,\sum dy,\sum |dy|\right]$$

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

圖3 不同的影象密度模式得到的不同的描述子結果

為了充分利用積分影象進行Haar小波的響應計算，我們並不直接旋轉Haar小波模板求得其響應值，而是在積影象上先使用水平和垂直的Haar模板求得響應值dy和dx，然後根據主方向旋轉dx和dy與主方向操持一致，如下圖4所示。為了求得旋轉後Haar小波響應值，首先要得到旋轉前影象的位置。旋轉前後圖偈的位置關係，可以通過點的旋轉公式得到：

$$x = x_0 –j\times scale\times sin(\theta)+i\times scale\times cos(\theta)$$

$$y = y_0 –j\times scale\times cos(\theta)+i\times scale\times sin(\theta)$$

在得到點$(j,i)$在旋轉前對應積分影象的位置$(x,y)$後，利用積分影象與水平、垂直Harr小波，求得水平與垂直兩個方向的響應值dx和dy。對dx和dy進行高斯加權處理，並根據主方向的角度，對dx和dy進行旋轉變換，從而，得到旋轉後的dx’和dy’。其計算公式如下：

$$dx’ = w(-dx\times sin(\theta)+dy\times cos(\theta))$$

$$dy’ = w(-dx\times cos(\theta)+dy\times sin(\theta))$$

圖4 利用積分影象進行Haar小波響應計算示意圖，左邊為旋轉後的影象，右邊為旋轉前的影象

3. 特徵描述子的維數

一般而言，特徵向量的長度越長，特徵向量所承載的資訊量就越大，特徵描述子的獨特性就越好，但匹配時所付出的時間代價就越大。對於SURF描述子，可以將它擴充套件到用128維向量來表示。具體方法是在求$\sum dx$、$\sum |dx|$時區分$dy<0$和$dy\ge 0$情況。同時，在求取$\sum dy$、$\sum |dy|$時區分$dx<0$和$dx\ge 0$情況。這樣，每個子塊就產生了8個梯度統計值，從而使描述子特徵向量的長度增加到$8\times4\times4=128$維。

為了實現快速匹配，SURF在特徵向量中增加了一個新的變數，即特徵點的拉普拉斯響應正負號。在特徵點檢測時，將Hessian矩陣的跡的正負號記錄下來，作為特徵向量中的一個變數。這樣做並不增加運算量，因為特徵點檢測進已經對Hessian矩陣的跡進行了計算。在特徵匹配時，這個變數可以有效地節省搜尋的時間，因為只有兩個具有相同正負號的特徵點才有可能匹配，對於正負號不同的特徵點就不進行相似性計算。

簡單地說，我們可以根據特徵點的響應值符號，將特徵點分成兩組，一組是具有拉普拉斯正響應的特徵點，一組是具有拉普拉斯負響應的特徵點，匹配時，只有符號相同組中的特徵點才能進行相互匹配。顯然，這樣可以節省特徵點匹配的時間。如下圖5所示。

圖5 黑背景下的亮斑和白背景下的黑斑 因為它們的拉普拉斯響應正負號不同，不會對它們進行匹配

4. 原始碼解析

特徵點描述子的生成這一部分的程式碼主要是通過SURFInvoker這個類來實現。在主流程中，通過一個parallel_for_()函式來併發計算。

struct SURFInvoker
{
    enum{ORI_RADIUS = 6, ORI_WIN = 60, PATCH_SZ = 20};
    // Parameters
    const Mat* img;
    const Mat* sum;
    vector<KeyPoint>* keypoints;
    Mat* descriptors;
    bool extended;
    bool upright;

    // Pre-calculated values
    int nOriSamples;
    vector<Point> apt; // 特徵點周圍用於描述方向的鄰域的點
    vector<float> aptw; // 描述 方向時的 高斯 權
    vector<float> DW;


    SURFInvoker(const Mat& _img, const Mat& _sum,
        vector<KeyPoint>& _keypoints, Mat& _descriptors,
        bool _extended, bool _upright)
    {
        keypoints = &_keypoints;
        descriptors = &_descriptors;
        img = &_img;
        sum = &_sum;
        extended = _extended;
        upright = _upright;

        // 用於描述特徵點的 方向的 鄰域大小： 12*sigma+1 (sigma =1.2) 因為高斯加權的核的引數為2sigma
        // nOriSampleBound為 矩形框內點的個數
        const int nOriSampleBound = (2 * ORI_RADIUS + 1)*(2 * ORI_RADIUS + 1); // 這裡把s近似為1 ORI_DADIUS = 6s

        // 分配大小 
        apt.resize(nOriSampleBound);
        aptw.resize(nOriSampleBound);
        DW.resize(PATCH_SZ*PATCH_SZ); // PATHC_SZ為特徵描述子的 區域大小 20s(s 這裡初始為1了)

        /* 計算特徵點方向用的 高斯分佈 權值與座標 */
        Mat G_ori = getGaussianKernel(2 * ORI_RADIUS + 1, SURF_ORI_SIGMA, CV_32F); // SURF_ORI_SIGMA = 1.2 *2 =2.5
        nOriSamples = 0;
        for (int i = -ORI_RADIUS; i <= ORI_RADIUS; i++)
        {
            for (int j = -ORI_RADIUS; j <= ORI_RADIUS; j++)
            {
                if (i*i + j*j <= ORI_RADIUS*ORI_RADIUS) // 限制在圓形區域內
                {
                    apt[nOriSamples] = cvPoint(i, j);
                    // 下面這裡有個座標轉換，因為i,j都是從-ORI_RADIUS開始的。
                    aptw[nOriSamples++] = G_ori.at<float>(i + ORI_RADIUS, 0) * G_ori.at<float>(j + ORI_RADIUS, 0);
                }
            }
        }
        CV_Assert(nOriSamples <= nOriSampleBound); // nOriSamples為圓形區域內的點，nOriSampleBound是正方形區域的點

        /* 用於特徵點描述子的高斯 權值 */
        Mat G_desc = getGaussianKernel(PATCH_SZ, SURF_DESC_SIGMA, CV_32F); // 用於生成特徵描述子的 高斯加權 sigma = 3.3s (s初取1)
        for (int i = 0; i < PATCH_SZ; i++)
        {
            for (int j = 0; j < PATCH_SZ; j++)
                DW[i*PATCH_SZ + j] = G_desc.at<float>(i, 0) * G_desc.at<float>(j, 0);
        }

        /* x與y方向上的 Harr小波,引數為4s */
        const int NX = 2, NY = 2;
        const int dx_s[NX][5] = { { 0, 0, 2, 4, -1 }, { 2, 0, 4, 4, 1 } };
        const int dy_s[NY][5] = { { 0, 0, 4, 2, 1 }, { 0, 2, 4, 4, -1 } };

        float X[nOriSampleBound], Y[nOriSampleBound], angle[nOriSampleBound]; // 用於計算特生點主方向
        uchar PATCH[PATCH_SZ + 1][PATCH_SZ + 1];
        float DX[PATCH_SZ][PATCH_SZ], DY[PATCH_SZ][PATCH_SZ]; // 20s * 20s區域的 梯度值
        CvMat matX = cvMat(1, nOriSampleBound, CV_32F, X);
        CvMat matY = cvMat(1, nOriSampleBound, CV_32F, Y);
        CvMat _angle = cvMat(1, nOriSampleBound, CV_32F, angle);
        Mat _patch(PATCH_SZ + 1, PATCH_SZ + 1, CV_8U, PATCH);

        int dsize = extended ? 128 : 64;

        int k, k1 = 0, k2 = (int)(*keypoints).size();// k2為Harr小波的 模板尺寸
        float maxSize = 0;
        for (k = k1; k < k2; k++)
        {
            maxSize = std::max(maxSize, (*keypoints)[k].size);
        }
        // maxSize*1.2/9 表示最大的尺度 s
        int imaxSize = std::max(cvCeil((PATCH_SZ + 1)*maxSize*1.2f / 9.0f), 1);
        Ptr<CvMat> winbuf = cvCreateMat(1, imaxSize*imaxSize, CV_8U);
        for (k = k1; k < k2; k++)
        {
            int i, j, kk, nangle;
            float* vec;
            SurfHF dx_t[NX], dy_t[NY];
            KeyPoint& kp = (*keypoints)[k];
            float size = kp.size;
            Point2f center = kp.pt;
            /* s是當前層的尺度引數 1.2是第一層的引數，9是第一層的模板大小*/
            float s = size*1.2f / 9.0f;
            /* grad_wav_size是 harr梯度模板的大小 邊長為 4s */
            int grad_wav_size = 2 * cvRound(2 * s);
            if (sum->rows < grad_wav_size || sum->cols < grad_wav_size)
            {
                /* when grad_wav_size is too big,
                * the sampling of gradient will be meaningless
                * mark keypoint for deletion. */
                kp.size = -1;
                continue;
            }

            float descriptor_dir = 360.f - 90.f;
            if (upright == 0)
            {
                // 這一步 是計算梯度值，先將harr模板放大，再根據積分圖計算，與前面求D_x，D_y一致類似
                resizeHaarPattern(dx_s, dx_t, NX, 4, grad_wav_size, sum->cols);
                resizeHaarPattern(dy_s, dy_t, NY, 4, grad_wav_size, sum->cols);
                for (kk = 0, nangle = 0; kk < nOriSamples; kk++)
                {
                    int x = cvRound(center.x + apt[kk].x*s - (float)(grad_wav_size - 1) / 2);
                    int y = cvRound(center.y + apt[kk].y*s - (float)(grad_wav_size - 1) / 2);
                    if (y < 0 || y >= sum->rows - grad_wav_size ||
                        x < 0 || x >= sum->cols - grad_wav_size)
                        continue;
                    const int* ptr = &sum->at<int>(y, x);
                    float vx = calcHaarPattern(ptr, dx_t, 2);
                    float vy = calcHaarPattern(ptr, dy_t, 2);
                    X[nangle] = vx*aptw[kk];
                    Y[nangle] = vy*aptw[kk];
                    nangle++;
                }
                if (nangle == 0)
                {
                    // No gradient could be sampled because the keypoint is too
                    // near too one or more of the sides of the image. As we
                    // therefore cannot find a dominant direction, we skip this
                    // keypoint and mark it for later deletion from the sequence.
                    kp.size = -1;
                    continue;
                }
                matX.cols = matY.cols = _angle.cols = nangle;
                // 計算鄰域內每個點的 梯度角度
                cvCartToPolar(&matX, &matY, 0, &_angle, 1);

                float bestx = 0, besty = 0, descriptor_mod = 0;
                for (i = 0; i < 360; i += SURF_ORI_SEARCH_INC) // SURF_ORI_SEARCH_INC 為扇形區域掃描的步長
                {
                    float sumx = 0, sumy = 0, temp_mod;
                    for (j = 0; j < nangle; j++)
                    {
                        // d是 分析到的那個點與 現在主方向的偏度
                        int d = std::abs(cvRound(angle[j]) - i);
                        if (d < ORI_WIN / 2 || d > 360 - ORI_WIN / 2)
                        {
                            sumx += X[j];
                            sumy += Y[j];
                        }
                    }
                    temp_mod = sumx*sumx + sumy*sumy;
                    // descriptor_mod 是最大峰值
                    if (temp_mod > descriptor_mod)
                    {
                        descriptor_mod = temp_mod;
                        bestx = sumx;
                        besty = sumy;
                    }
                }
                descriptor_dir = fastAtan2(-besty, bestx);
            }
            kp.angle = descriptor_dir;
            if (!descriptors || !descriptors->data)
                continue;

            /* 用特徵點周圍20*s為邊長的鄰域 計算特徵描述子 */
            int win_size = (int)((PATCH_SZ + 1)*s);
            CV_Assert(winbuf->cols >= win_size*win_size);
            Mat win(win_size, win_size, CV_8U, winbuf->data.ptr);

            if (!upright)
            {
                descriptor_dir *= (float)(CV_PI / 180); // 特徵點的主方向 弧度值
                float sin_dir = -std::sin(descriptor_dir); //  - sin dir
                float cos_dir = std::cos(descriptor_dir);

                float win_offset = -(float)(win_size - 1) / 2;
                float start_x = center.x + win_offset*cos_dir + win_offset*sin_dir;
                float start_y = center.y - win_offset*sin_dir + win_offset*cos_dir;
                uchar* WIN = win.data;

                int ncols1 = img->cols - 1, nrows1 = img->rows - 1;
                size_t imgstep = img->step;
                for (i = 0; i < win_size; i++, start_x += sin_dir, start_y += cos_dir)
                {
                    double pixel_x = start_x;
                    double pixel_y = start_y;
                    for (j = 0; j < win_size; j++, pixel_x += cos_dir, pixel_y -= sin_dir)
                    {
                        int ix = cvFloor(pixel_x), iy = cvFloor(pixel_y);
                        if ((unsigned)ix < (unsigned)ncols1 &&
                            (unsigned)iy < (unsigned)nrows1)
                        {
                            float a = (float)(pixel_x - ix), b = (float)(pixel_y - iy);
                            const uchar* imgptr = &img->at<uchar>(iy, ix);
                            WIN[i*win_size + j] = (uchar)
                                cvRound(imgptr[0] * (1.f - a)*(1.f - b) +
                                imgptr[1] * a*(1.f - b) +
                                imgptr[imgstep] * (1.f - a)*b +
                                imgptr[imgstep + 1] * a*b);
                        }
                        else
                        {
                            int x = std::min(std::max(cvRound(pixel_x), 0), ncols1);
                            int y = std::min(std::max(cvRound(pixel_y), 0), nrows1);
                            WIN[i*win_size + j] = img->at<uchar>(y, x);
                        }
                    }
                }
            }
            else
            {

                float win_offset = -(float)(win_size - 1) / 2;
                int start_x = cvRound(center.x + win_offset);
                int start_y = cvRound(center.y - win_offset);
                uchar* WIN = win.data;
                for (i = 0; i < win_size; i++, start_x++)
                {
                    int pixel_x = start_x;
                    int pixel_y = start_y;
                    for (j = 0; j < win_size; j++, pixel_y--)
                    {
                        int x = MAX(pixel_x, 0);
                        int y = MAX(pixel_y, 0);
                        x = MIN(x, img->cols - 1);
                        y = MIN(y, img->rows - 1);
                        WIN[i*win_size + j] = img->at<uchar>(y, x);
                    }
                }
            }
            // Scale the window to size PATCH_SZ so each pixel's size is s. This
            // makes calculating the gradients with wavelets of size 2s easy
            resize(win, _patch, _patch.size(), 0, 0, INTER_AREA);

            // Calculate gradients in x and y with wavelets of size 2s
            for (i = 0; i < PATCH_SZ; i++)
            for (j = 0; j < PATCH_SZ; j++)
            {
                float dw = DW[i*PATCH_SZ + j]; // 高斯加權係數
                float vx = (PATCH[i][j + 1] - PATCH[i][j] + PATCH[i + 1][j + 1] - PATCH[i + 1][j])*dw;
                float vy = (PATCH[i + 1][j] - PATCH[i][j] + PATCH[i + 1][j + 1] - PATCH[i][j + 1])*dw;
                DX[i][j] = vx;
                DY[i][j] = vy;
            }

            // Construct the descriptor
            vec = descriptors->ptr<float>(k);
            for (kk = 0; kk < dsize; kk++)
                vec[kk] = 0;
            double square_mag = 0;
            if (extended)
            {
                // 128維描述子，考慮dx與dy的正負號
                for (i = 0; i < 4; i++)
                for (j = 0; j < 4; j++)
                {
                    // 每個方塊內是一個5s * 5s的區域，每個方法由8個特徵描述
                    for (int y = i * 5; y < i * 5 + 5; y++)
                    {
                        for (int x = j * 5; x < j * 5 + 5; x++)
                        {
                            float tx = DX[y][x], ty = DY[y][x];
                            if (ty >= 0)
                            {
                                vec[0] += tx;
                                vec[1] += (float)fabs(tx);
                            }
                            else {
                                vec[2] += tx;
                                vec[3] += (float)fabs(tx);
                            }
                            if (tx >= 0)
                            {
                                vec[4] += ty;
                                vec[5] += (float)fabs(ty);
                            }
                            else {
                                vec[6] += ty;
                                vec[7] += (float)fabs(ty);
                            }
                        }
                    }
                    for (kk = 0; kk < 8; kk++)
                        square_mag += vec[kk] * vec[kk];
                    vec += 8;
                }
            }
            else
            {
                // 64位描述子
                for (i = 0; i < 4; i++)
                for (j = 0; j < 4; j++)
                {
                    for (int y = i * 5; y < i * 5 + 5; y++)
                    {
                        for (int x = j * 5; x < j * 5 + 5; x++)
                        {
                            float tx = DX[y][x], ty = DY[y][x];
                            vec[0] += tx; vec[1] += ty;
                            vec[2] += (float)fabs(tx); vec[3] += (float)fabs(ty);
                        }
                    }
                    for (kk = 0; kk < 4; kk++)
                        square_mag += vec[kk] * vec[kk];
                    vec += 4;
                }
            }
            // 歸一化 描述子 以滿足 光照不變性
            vec = descriptors->ptr<float>(k);
            float scale = (float)(1. / (sqrt(square_mag) + DBL_EPSILON));
            for (kk = 0; kk < dsize; kk++)
                vec[kk] *= scale;
        }
    }
};

5. 總結

實際上有文獻指出，SURF比SIFT工作更出色。他們認為主要是因為SURF在求取描述子特徵向量時，是對一個子塊的梯度資訊進行求和，而SIFT則是依靠單個畫素梯度的方向。