SURF特征匹原理及源代碼

SURF是對SIFT的改進，相對於SIFT，主要優點是速度更快，更適合做實時特征檢查。

 

 

一、SURF原理：

　　相對於SIFT，SUFT采用Hssian算法檢測關鍵點，很大程度上提高了程序的運行速度，同時，在尺度空間的構建上，SIFT通過改變高斯卷積核的大小，構建不同的組，下文進行詳細介紹。

 

二、SURF的實現：

　　與SIFT步驟相同，可分為四步：尺度空間構建和極值檢測、特征點精確定位和過濾、特征方向賦值和特征點描述。

　　尺度空間構建和極值檢測：

　　　　SIFT通過高斯濾波后的同組圖像中的每層相減得到DOG空間來提取特征點，SURF通過Harris矩陣（極大提升速度的原因），求取X和Y的導數，進行角點的大致判斷。下面簡要介紹實現步驟：

　　　　　　　Harris矩陣的表示：　　　　

 

 　　上式中的L表示在經過高斯濾波后的圖像的像素值的大小，L_xx、L_xy、L_yy分別表示在x方向的二次偏導導、在x和y方向上的二次偏導和在y方向的二次偏導。在進行Harris矩陣求取時，需滿足尺度不變特性。即需先進行高斯濾波，得到多尺度的金字塔圖像（具體構造會在下面介紹）。根據上式可得到Harris矩陣表示為行列式的值的大小：

　　　　　　　　　　Det(H) = L_xxL_yy - L_xyL_xy

　　　　根據圖像像素點的離散性可知，圖像像素值的求導可通過相鄰像素點像素值的相減得到：

　　　　　　　　　　L_x = L(x+1,y) - L(x,y)　　　　L_xx = L(x+1,y) - L(x,y) - L(x,y)+L(x-1,y)  = L(x+1,y)+L(x-1,y)-2L(x,y)

　　　同理可求得L_xx、L_xy和L_yy的表達式，上式中的L均為經過高斯濾波后的像素值大小。由上面表達式可知，得到Harris響應值，需要先經過高斯平滑和二次求導兩步操作，SURF算法作者采用箱式濾波模板近似代替上面兩步，加快了求導的速度。下面暫時介紹SURF中箱式濾波原理（這個箱式濾波，在SURF中的使用方法，，，頭疼、自閉、脫發！！！！！！！）：

　　　SURF中箱式濾波：

　　　　　　箱式濾波可簡單理解為將像素點周圍划分為多個區域，再將每個局部區域內的像素值相加再乘以一個權重值后，進行相加。

　　　　　　　　

 

 

 　　　　上圖中，第一幅表示為經過高斯濾波后在y方向的二階求導，第二幅圖像為經過高斯濾波后，在xy方向的求導，第三幅圖像表示為箱式濾波在y方向的求導，第四幅圖表示箱式濾波為在xy方向上的求導。在采用箱式濾波的過程中，針對每個局部區域像素值的計算，可通過積分圖像進行快速計算，加快了計算的速度。（積分圖像在另一篇博客中已有介紹）

　　　　在采用箱式卷積核替換后，Harris響應值的計算可用下式表示：

　　　　　　　　　　　　Det(H) = D_xxD_yy -(wD_xy)²//其中D表示盒子模板與圖像的近似卷積  其中w為更好模擬而設的變量（原算法中建議取0.9）

　　　  通過判斷 Det(H)結果的符號，可簡要判斷是否為局部極值點，若為負值，則不是局部極值點，若為正值，則可能歸類為局部極值點。

尺寸空間構建：

　　　　相對於SIFT，SURF構建的圖像金字塔不同組圖像的尺寸大小相同，不同的是卷積核的大小（在此，，，查了很多博客，，，發現都寫成高斯卷積核大小，，，很想罵人，，，在開始寫為高斯卷積核，在計算Harris相應值時，又改為箱式濾波核，糾結了很久，看了不少博客，最終應該表示為箱式卷積核的大小）。下圖為一個簡要概括：

　　  　　　　

 　　　非極值抑制：

　　　　　　通過上面得到極值點，將其響應值大小與上下兩層響應值大小進行比較，進行判斷是否為極值。（原理和SIFT相同）

 

主方向的計算：

　　　　　　為確保得到的描述子具有旋轉不變性，需對每一個特征點分配一個主方向。具體方法為：以特征點為中心，以6s（s為特征點的尺度）為半徑的圓形內，對圖像進行Haar小波響應運算（相當於對圖像進行梯度運算，但采用積分圖像，更有效率。Haar小波響應運算暫不介紹，如有興趣，可網上百度查詢）。后設置一個以方向為中心，以pi/3為張角的扇形滑動窗口，通過一定的弧度（大概0.2）為步長，旋轉這個滑動窗口，對窗口內的圖像Haar小波的響應值進行累加。主方向最大的Harr響應累加值對應的方向即為主方向。

 

SURF特征點描述子計算：

 　　　　　　　　 在SURF算法中，也是在特征點周圍取一個正方形框，框的邊長為20s(s是所檢測到該特征點所在的尺度)。該框帶方向，方向當然就是上面一步檢測出來的主方向。然后把該框分為16個子區域，每個子區域統計25個像素的水平方向和垂直方向的haar小波特征，這里的水平方向和垂直方向是相對於主方向而言的。該haar小波特征為水平方向值之和、水平方向絕對值之和和垂直方向之和、垂直方向絕對值之和。這樣每個小區域就有4個值，所以每個特征點就是16*4維的向量，相比SIFT而言，少了一半，這在特征匹配的過程中會大大的加速匹配的速度。

 　在OpenSURF的實現源碼中采用的方式，通過點旋轉公式，把點旋轉到主方向上並進行最近鄰插值的對應點，公式如下：

 

 

 

 

 

OPENCV-contrib內接口：

    cv::Ptr<cv::xfeatures2d::SURF> surf = cv::xfeatures2d::SURF::create();

 

完整實現：

#include <opencv2/opencv.hpp>
#include <opencv2/xfeatures2d.hpp>

int main()
{
    cv::Mat imageL = cv::imread("Ｔ.png");
    cv::Mat imageR = cv::imread("-DT.png");


    //提取特征點方法


    cv::Ptr<cv::xfeatures2d::SURF> surf = cv::xfeatures2d::SURF::create();

    //特征點
    std::vector<cv::KeyPoint> keyPointL, keyPointR;
    //單獨提取特征點
    surf->detect(imageL, keyPointL);
    surf->detect(imageR, keyPointR);

    //畫特征點
    cv::Mat keyPointImageL;
    cv::Mat keyPointImageR;
    drawKeypoints(imageL, keyPointL, keyPointImageL, cv::Scalar::all(-1), cv::DrawMatchesFlags::DRAW_RICH_KEYPOINTS);
    drawKeypoints(imageR, keyPointR, keyPointImageR, cv::Scalar::all(-1), cv::DrawMatchesFlags::DRAW_RICH_KEYPOINTS);

    //顯示窗口
    cv::namedWindow("KeyPoints of imageL");
    cv::namedWindow("KeyPoints of imageR");

    //顯示特征點
    cv::imshow("KeyPoints of imageL", keyPointImageL);
    cv::imshow("KeyPoints of imageR", keyPointImageR);

    //特征點匹配
    cv::Mat despL, despR;
    //提取特征點並計算特征描述子
    surf->detectAndCompute(imageL, cv::Mat(), keyPointL, despL);
    surf->detectAndCompute(imageR, cv::Mat(), keyPointR, despR);

 
    std::vector<cv::DMatch> matches;

    //如果采用flannBased方法 那么 desp通過orb的到的類型不同需要先轉換類型
    if (despL.type() != CV_32F || despR.type() != CV_32F)
    {
        despL.convertTo(despL, CV_32F);
        despR.convertTo(despR, CV_32F);
    }

    cv::Ptr<cv::DescriptorMatcher> matcher = cv::DescriptorMatcher::create("FlannBased");
    matcher->match(despL, despR, matches);

    //計算特征點距離的最大值
    double maxDist = 0;
    for (int i = 0; i < despL.rows; i++)
    {
        double dist = matches[i].distance;
        if (dist > maxDist)
            maxDist = dist;
    }

    //挑選好的匹配點
    std::vector< cv::DMatch > good_matches;
    for (int i = 0; i < despL.rows; i++)
    {
        if (matches[i].distance < 0.03*maxDist)
        {
            good_matches.push_back(matches[i]);
        }
    }



    cv::Mat imageOutput;
    cv::drawMatches(imageL, keyPointL, imageR, keyPointR, good_matches, imageOutput);

    cv::namedWindow("picture of matching");
    cv::imshow("picture of matching", imageOutput);
    cv::waitKey(0);
    return 0;
}

 

源代碼：

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;
        }
    }
};

 

 

 

參考：

 https://blog.csdn.net/songzitea/article/details/16986423

https://blog.csdn.net/msq19895070/article/details/24290599

https://wenku.baidu.com/view/d71c9d2e905f804d2b160b4e767f5acfa1c783b4.html?rec_flag=default

https://blog.csdn.net/qq_37764129/article/details/80969515