FROM：?http://www.cnblogs.com/ronny/p/4048213.html

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

1. SURF特征點方向分配

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

與SIFT類似，使用σ=2s的高斯加權函數(shù)對Harr小波的響應值進行高斯加權。為了求取主方向值，需要設計一個以特征點為中心，張角為π/3的扇形滑動窗口。以步長為0.2弧度左右，旋轉(zhuǎn)這個滑動窗口，并對滑動窗口內(nèi)的圖像Harr小波響應值dx、dy進行累加，得到一個矢量(mw,θw)：

mw=∑wdx+∑wdy

θw=arctan(∑wdx/∑wdy)

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

θ=θw|max{mw}

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

圖 1? 求取主方向時扇形滑動窗口圍繞特征點轉(zhuǎn)動，統(tǒng)計Haar小波響應值，并計算方向角

2. 特征點特征矢量生成

生成特征點描述子與確定特征點方向有些類似，它需要計算圖像的Haar小波響應。不過，與主方向的確定不同的是，這次我們不是使用一個圓形區(qū)域，而是在一個矩形區(qū)域來計算Haar小波響應。以特征點為中心，沿上一節(jié)討論得到的主方向，沿主方向?qū)?span id="ze8trgl8bvbq" class="MathJax_Preview" style="margin:0px; padding:0px; color:rgb(136,136,136)">s20s×20s的圖像劃分為4×4個子塊，每個子塊利用尺寸2s的Harr模板進行響應值進行響應值計算，然后對響應值進行統(tǒng)計∑dx、∑|dx|、∑dy、∑|dy|形成特征矢量。如下圖2所示。圖中，以特征點為中心，以20s為邊長的矩形窗口為特征描述子計算使用的窗口，特征點到矩形邊框的線段表示特征點的主方向。

圖2 特征描述子表示

將20s的窗口劃分成4×4子窗口，每個子窗口有5s×5s個像素。使用尺寸為2s的Harr小波對子窗口圖像進行其響應值計算，共進行25次采樣，分別得到沿主方向的dy和垂直于主方向的dx。然后，以特征點為中心，對dy和dx進行高斯加權計算，高斯核的參數(shù)為σ=3.3s(即20s/6)。最后，分別對每個子塊的響應值進行統(tǒng)計，得到每個子塊的矢量：

V子塊=[∑dx,∑|dx|,∑dy,∑|dy|]

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

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

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

x=x0–j×scale×sin(θ)+i×scale×cos(θ)

y=y0–j×scale×cos(θ)+i×scale×sin(θ)

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

dx′=w(?dx×sin(θ)+dy×cos(θ))

dy′=w(?dx×cos(θ)+dy×sin(θ))

圖4 利用積分圖像進行Haar小波響應計算示意圖，左邊為旋轉(zhuǎn)后的圖像，右邊為旋轉(zhuǎn)前的圖像

3. 特征描述子的維數(shù)

一般而言，特征矢量的長度越長，特征矢量所承載的信息量就越大，特征描述子的獨特性就越好，但匹配時所付出的時間代價就越大。對于SURF描述子，可以將它擴展到用128維矢量來表示。具體方法是在求∑dx、∑|dx|時區(qū)分dy<0和dy≥0情況。同時，在求取∑dy、∑|dy|時區(qū)分dx<0和dx≥0情況。這樣，每個子塊就產(chǎn)生了8個梯度統(tǒng)計值，從而使描述子特征矢量的長度增加到8×4×4=128維。

為了實現(xiàn)快速匹配，SURF在特征矢量中增加了一個新的變量，即特征點的拉普拉斯響應正負號。在特征點檢測時，將Hessian矩陣的跡的正負號記錄下來，作為特征矢量中的一個變量。這樣做并不增加運算量，因為特征點檢測進已經(jīng)對Hessian矩陣的跡進行了計算。在特征匹配時，這個變量可以有效地節(jié)省搜索的時間，因為只有兩個具有相同正負號的特征點才有可能匹配，對于正負號不同的特征點就不進行相似性計算。

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

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

?

4. 源碼解析

特征點描述子的生成這一部分的代碼主要是通過SURFInvoker這個類來實現(xiàn)。在主流程中，通過一個parallel_for_()函數(shù)來并發(fā)計算。

struct SURFInvoker {enum{ORI_RADIUS = 6, ORI_WIN = 60, PATCH_SZ = 20};// Parametersconst Mat* img;const Mat* sum;vector<KeyPoint>* keypoints;Mat* descriptors;bool extended;bool upright;// Pre-calculated valuesint 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) 因為高斯加權的核的參數(shù)為2sigma// nOriSampleBound為 矩形框內(nèi)點的個數(shù)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為特征描述子的 區(qū)域大小 20s(s 這里初始為1了)/* 計算特征點方向用的 高斯分布 權值與坐標 */Mat G_ori = getGaussianKernel(2 * ORI_RADIUS + 1, SURF_ORI_SIGMA, CV_32F); // SURF_ORI_SIGMA = 1.2 *2 =2.5nOriSamples = 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) // 限制在圓形區(qū)域內(nèi) {apt[nOriSamples] = cvPoint(i, j);// 下面這里有個坐標轉(zhuǎn)換，因為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為圓形區(qū)域內(nèi)的點，nOriSampleBound是正方形區(qū)域的點/* 用于特征點描述子的高斯 權值 */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小波,參數(shù)為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區(qū)域的 梯度值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 表示最大的尺度 sint 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是當前層的尺度參數(shù) 1.2是第一層的參數(shù)，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模板放大，再根據(jù)積分圖計算，與前面求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;// 計算鄰域內(nèi)每個點的 梯度角度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 為扇形區(qū)域掃描的步長 {float sumx = 0, sumy = 0, temp_mod;for (j = 0; j < nangle; j++){// d是 分析到的那個點與 現(xiàn)在主方向的偏度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 dirfloat 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 easyresize(win, _patch, _patch.size(), 0, 0, INTER_AREA);// Calculate gradients in x and y with wavelets of size 2sfor (i = 0; i < PATCH_SZ; i++)for (j = 0; j < PATCH_SZ; j++){float dw = DW[i*PATCH_SZ + j]; // 高斯加權系數(shù)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 descriptorvec = 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++){// 每個方塊內(nèi)是一個5s * 5s的區(qū)域，每個方法由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則是依靠單個像素梯度的方向。

總結

以上是生活随笔為你收集整理的SURF算法与源码分析、下的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。