Surf(Speed Up Robust Feature)
Surf算法的原理 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??
1.構建Hessian矩陣構造高斯金字塔尺度空間
其實surf構造的金字塔圖像與sift有很大不同，就是因為這些不同才加快了其檢測的速度。Sift采用的是DOG圖像，而surf采用的是Hessian矩陣行列式近似值圖像。Hessian矩陣是Surf算法的核心，為了方便運算，假設函數f(z，y)，Hessian矩陣H是由函數，偏導數組成。首先來看看圖像中某個像素點的Hessian矩陣，如下：
即每一個像素點都可以求出一個Hessian矩陣。
H矩陣判別式為：
判別式的值是H矩陣的特征值，可以利用判定結果的符號將所有點分類，根據判別式取值正負，來判別該點是或不是極值點。
在SURF算法中，用圖像像素l(x，y)即為函數值f(x，y)，選用二階標準高斯函數作為濾波器，通過特定核間的卷積計算二階偏導數，這樣便能計算出H矩陣的三個矩陣元素L_xx,L_xy,L_yy 從而計算出H矩陣：
但是由于我們的特征點需要具備尺度無關性，所以在進行Hessian矩陣構造前，需要對其進行高斯濾波。
這樣，經過濾波后在進行Hessian的計算，其公式如下：
L(x,t)是一幅圖像在不同解析度下的表示，可以利用高斯核G(t)與圖像函數 I(x) 在點x的卷積來實現，其中高斯核G(t)：
g(x)為高斯函數，t為高斯方差。通過這種方法可以為圖像中每個像素計算出其H行列式的決定值，并用這個值來判別特征點。為方便應用，Herbert Bay提出用近似值現代替L(x,t)。為平衡準確值與近似值間的誤差引入權值叫，權值隨尺度變化，則H矩陣判別式可表示為：
其中0.9是作者給出的一個經驗值，其實它是有一套理論計算的，具體去看surf的英文論文。
由于求Hessian時要先高斯平滑，然后求二階導數，這在離散的像素點是用模板卷積形成的，這2中操作合在一起用一個模板代替就可以了，比如說y方向上的模板如下：
該圖的左邊即用高斯平滑然后在y方向上求二階導數的模板，為了加快運算用了近似處理，其處理結果如右圖所示，這樣就簡化了很多。并且右圖可以采用積分圖來運算，大大的加快了速度，關于積分圖的介紹，可以去查閱相關的資料。
同理，x和y方向的二階混合偏導模板如下所示：
上面講的這么多只是得到了一張近似hessian行列式圖，這類似sift中的DOG圖，但是在金字塔圖像中分為很多層，每一層叫做一個octave，每一個octave中又有幾張尺度不同的圖片。在sift算法中，同一個octave層中的圖片尺寸(即大小)相同，但是尺度(即模糊程度)不同，而不同的octave層中的圖片尺寸大小也不相同，因為它是由上一層圖片降采樣得到的。在進行高斯模糊時，sift的高斯模板大小是始終不變的，只是在不同的octave之間改變圖片的大小。而在surf中，圖片的大小是一直不變的，不同的octave層得到的待檢測圖片是改變高斯模糊尺寸大小得到的，當然了，同一個octave中個的圖片用到的高斯模板尺度也不同。算法允許尺度空間多層圖像同時被處理，不需對圖像進行二次抽樣，從而提高算法性能。左圖是傳統方式建立一個如圖所示的金字塔結構，圖像的寸是變化的，并且運 算會反復使用高斯函數對子層進行平滑處理，右圖說明Surf算法使原始圖像保持不變而只改變濾波器大小。Surf采用這種方法節省了降采樣過程，其處理速度自然也就提上去了。其金字塔圖像如下所示：
2. 利用非極大值抑制初步確定特征點
? ? ? ? 此步驟和sift類似，將經過hessian矩陣處理過的每個像素點與其3維領域的26個點進行大小比較，如果它是這26個點中的最大值或者最小值，則保留下來，當做初步的特征點。檢測過程中使用與該尺度層圖像解析度相對應大小的濾波器進行檢測，以3×3的濾波器為例，該尺度層圖像中9個像素點之一圖2檢測特征點與自身尺度層中其余8個點和在其之上及之下的兩個尺度層9個點進行比較，共26個點，圖2中標記‘x’的像素點的特征值若大于周圍像素則可確定該點為該區域的特征點。
?
3. 精確定位極值點
? ? ? ? ?這里也和sift算法中的類似，采用3維線性插值法得到亞像素級的特征點，同時也去掉那些值小于一定閾值的點，增加極值使檢測到的特征點數量減少，最終只有幾個特征最強點會被檢測出來。
?
4. 選取特征點的主方向。
? ? ? ? 這一步與sift也大有不同。Sift選取特征點主方向是采用在特征點領域內統計其梯度直方圖，取直方圖bin值最大的以及超過最大bin值80%的那些方向做為特征點的主方向。
? ? ? ?而在surf中，不統計其梯度直方圖，而是統計特征點領域內的harr小波特征。即在特征點的領域(比如說，半徑為6s的圓內，s為該點所在的尺度)內，統計60度扇形內所有點的水平haar小波特征和垂直haar小波特征總和，haar小波的尺寸變長為4s，這樣一個扇形得到了一個值。然后60度扇形以一定間隔進行旋轉，最后將最大值那個扇形的方向作為該特征點的主方向。該過程的示意圖如下：
5. 構造surf特征點描述算子
? ? ? ? 在sift中，是在特征點周圍取16*16的鄰域，并把該領域化為4*4個的小區域，每個小區域統計8個方向梯度，最后得到4*4*8=128維的向量，該向量作為該點的sift描述子。
在surf中，也是在特征點周圍取一個正方形框，框的邊長為20s(s是所檢測到該特征點所在的尺度)。該框帶方向，方向當然就是第4步檢測出來的主方向了。然后把該框分為16個子區域，每個子區域統計25個像素的水平方向和垂直方向的haar小波特征，這里的水平和垂直方向都是相對主方向而言的。該haar小波特征為水平方向值之和，水平方向絕對值之和，垂直方向之和，垂直方向絕對值之和。該過程的示意圖如下所示：
這樣每個小區域就有4個值，所以每個特征點就是16*4=64維的向量，相比sift而言，少了一半，這在特征匹配過程中會大大加快匹配速度。
6.結束語
Surf采用Henssian矩陣獲取圖像局部最值還是十分穩定的，但是在求主方向階段太過于依賴局部區域像素的梯度方向，有可能使得找到的主方向不準確，后面的特征向量提取以及匹配都嚴重依賴于主方向，即使不大偏差角度也可以造成后面特征匹配的放大誤差，從而匹配不成功；另外圖像金字塔的層取 得不足夠緊密也會使得尺度有誤差，后面的特征向量提取同樣依賴相應的尺度，發明者在這個問題上的折中解決方法是取適量的層然后進行插值。
Sift是一種只 利用到灰度性質的算法，忽略了色彩信息，后面又出現了幾種據說比Surf更穩定的描述器其中一些利用到了色彩信息，讓我們拭目以待吧。
代碼： ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??
來源：OpenCV/sample/c中的find_obj.cpp代碼
?
需仔細注意：
1.定位部分：通過透視變換，畫出了目標在圖像中的位置，但是這么做會浪費很多時間，可以改進：
2.flann尋找最近的臨近Keypoints：
首先，利用圖像，構建多維查找樹，然后，利用Knn算法找到最近的Keypoints (KNN算法：http://blog.csdn.net/sangni007/article/details/7482890）
[cpp]view plaincopyprint?
?? ?? cv::flann::Index?flann_index(m_image,?cv::flann::KDTreeIndexParams(4));???? ?? flann_index.knnSearch(m_object,?m_indices,?m_dists,?2,?cv::flann::SearchParams(64)?);???
flann算法有很多功能，
文檔：http://opencv.itseez.com/modules/flann/doc/flann_fast_approximate_nearest_neighbor_search.html?highlight=flann#fast-approximate-nearest-neighbor-search
?
[cpp]view plaincopyprint?
? ? ? ? ? ?? #include?"opencv2/objdetect/objdetect.hpp"?? #include?"opencv2/features2d/features2d.hpp"?? #include?"opencv2/highgui/highgui.hpp"?? #include?"opencv2/calib3d/calib3d.hpp"?? #include?"opencv2/imgproc/imgproc_c.h"#include?<iostream>?? #include?<vector>?? #include?<stdio.h>using?namespace?std;?? void?help()?? {?? ????printf(?? ????????"This?program?demonstrated?the?use?of?the?SURF?Detector?and?Descriptor?using\n"?? ????????"either?FLANN?(fast?approx?nearst?neighbor?classification)?or?brute?force?matching\n"?? ????????"on?planar?objects.\n"?? ????????"Usage:\n"?? ????????"./find_obj?<object_filename>?<scene_filename>,?default?is?box.png??and?box_in_scene.png\n\n");?? ????return;?? }?? #define?USE_FLANN?? IplImage*?image?=?0;double?compareSURFDescriptors(?const?float*?d1,?const?float*?d2,?double?best,?int?length?)?? {?? ????double?total_cost?=?0;?? ????assert(?length?%?4?==?0?);?? ????for(?int?i?=?0;?i?<?length;?i?+=?4?)?? ????{?? ????????double?t0?=?d1[i??]?-?d2[i??];?? ????????double?t1?=?d1[i+1]?-?d2[i+1];?? ????????double?t2?=?d1[i+2]?-?d2[i+2];?? ????????double?t3?=?d1[i+3]?-?d2[i+3];?? ????????total_cost?+=?t0*t0?+?t1*t1?+?t2*t2?+?t3*t3;?? ????????if(?total_cost?>?best?)?? ????????????break;?? ????}?? ????return?total_cost;?? }?? int?naiveNearestNeighbor(?const?float*?vec,?int?laplacian,?? ??????????????????????const?CvSeq*?model_keypoints,?? ??????????????????????const?CvSeq*?model_descriptors?)?? {?? ????int?length?=?(int)(model_descriptors->elem_size/sizeof(float));?? ????int?i,?neighbor?=?-1;?? ????double?d,?dist1?=?1e6,?dist2?=?1e6;?? ????CvSeqReader?reader,?kreader;?? ????cvStartReadSeq(?model_keypoints,?&kreader,?0?);?? ????cvStartReadSeq(?model_descriptors,?&reader,?0?);????for(?i?=?0;?i?<?model_descriptors->total;?i++?)?? ????{?? ????????const?CvSURFPoint*?kp?=?(const?CvSURFPoint*)kreader.ptr;?? ????????const?float*?mvec?=?(const?float*)reader.ptr;?? ?????CV_NEXT_SEQ_ELEM(?kreader.seq->elem_size,?kreader?);?? ????????CV_NEXT_SEQ_ELEM(?reader.seq->elem_size,?reader?);?? ????????if(?laplacian?!=?kp->laplacian?)?? ????????????continue;?? ????????d?=?compareSURFDescriptors(?vec,?mvec,?dist2,?length?);?? ????????if(?d?<?dist1?)?? ????????{?? ????????????dist2?=?dist1;?? ????????????dist1?=?d;?? ????????????neighbor?=?i;?? ????????}?? ????????else?if?(?d?<?dist2?)?? ????????????dist2?=?d;?? ????}?? ????if?(?dist1?<?0.6*dist2?)?? ????????return?neighbor;?? ????return?-1;?? }?? ?? ?? void?findPairs(?const?CvSeq*?objectKeypoints,?const?CvSeq*?objectDescriptors,?? ???????????const?CvSeq*?imageKeypoints,?const?CvSeq*?imageDescriptors,?vector<int>&?ptpairs?)?? {?? ????int?i;?? ????CvSeqReader?reader,?kreader;?? ????cvStartReadSeq(?objectKeypoints,?&kreader?);?? ????cvStartReadSeq(?objectDescriptors,?&reader?);?? ????ptpairs.clear();????for(?i?=?0;?i?<?objectDescriptors->total;?i++?)?? ????{?? ????????const?CvSURFPoint*?kp?=?(const?CvSURFPoint*)kreader.ptr;?? ????????const?float*?descriptor?=?(const?float*)reader.ptr;?? ????????CV_NEXT_SEQ_ELEM(?kreader.seq->elem_size,?kreader?);?? ????????CV_NEXT_SEQ_ELEM(?reader.seq->elem_size,?reader?);?? ????????int?nearest_neighbor?=?naiveNearestNeighbor(?descriptor,?kp->laplacian,?imageKeypoints,?imageDescriptors?);?? ????????if(?nearest_neighbor?>=?0?)?? ????????{?? ????????????ptpairs.push_back(i);?? ????????????ptpairs.push_back(nearest_neighbor);?? ????????}?? ????}?? }?? void?flannFindPairs(?const?CvSeq*,?const?CvSeq*?objectDescriptors,?? ???????????const?CvSeq*,?const?CvSeq*?imageDescriptors,?vector<int>&?ptpairs?)?? {?? ?int?length?=?(int)(objectDescriptors->elem_size/sizeof(float));????cv::Mat?m_object(objectDescriptors->total,?length,?CV_32F);?? ?cv::Mat?m_image(imageDescriptors->total,?length,?CV_32F);?? ??? ????CvSeqReader?obj_reader;?? ?float*?obj_ptr?=?m_object.ptr<float>(0);?? ????cvStartReadSeq(?objectDescriptors,?&obj_reader?);?? ??? ????for(int?i?=?0;?i?<?objectDescriptors->total;?i++?)?? ????{?? ????????const?float*?descriptor?=?(const?float*)obj_reader.ptr;?? ????????CV_NEXT_SEQ_ELEM(?obj_reader.seq->elem_size,?obj_reader?);?? ????????memcpy(obj_ptr,?descriptor,?length*sizeof(float));?? ????????obj_ptr?+=?length;?? ????}?? ??? ????CvSeqReader?img_reader;?? ?float*?img_ptr?=?m_image.ptr<float>(0);?? ????cvStartReadSeq(?imageDescriptors,?&img_reader?);?? ????for(int?i?=?0;?i?<?imageDescriptors->total;?i++?)?? ????{?? ????????const?float*?descriptor?=?(const?float*)img_reader.ptr;?? ????????CV_NEXT_SEQ_ELEM(?img_reader.seq->elem_size,?img_reader?);?? ????????memcpy(img_ptr,?descriptor,?length*sizeof(float));?? ????????img_ptr?+=?length;?? ????}?????? ????cv::Mat?m_indices(objectDescriptors->total,?2,?CV_32S);?? ????cv::Mat?m_dists(objectDescriptors->total,?2,?CV_32F);?? ??? ??? ????cv::flann::Index?flann_index(m_image,?cv::flann::KDTreeIndexParams(4));???? ??? ????flann_index.knnSearch(m_object,?m_indices,?m_dists,?2,?cv::flann::SearchParams(64)?);??? ????float*?dists_ptr?=?m_dists.ptr<float>(0);?? ????for?(int?i=0;i<m_indices.rows;++i)??? ?{?? ?????if?(dists_ptr[2*i]<0.6*dists_ptr[2*i+1])?? ??{?? ??????ptpairs.push_back(i);?? ??????ptpairs.push_back(indices_ptr[2*i]);?? ?????}?? ????}?? }?? ?? ?? int?locatePlanarObject(?const?CvSeq*?objectKeypoints,?const?CvSeq*?objectDescriptors,?? ?const?CvSeq*?imageKeypoints,?const?CvSeq*?imageDescriptors,?? ?const?CvPoint?src_corners[4],?CvPoint?dst_corners[4]?)?? {?? ????double?h[9];?? ????CvMat?_h?=?cvMat(3,?3,?CV_64F,?h);?? ????vector<int>?ptpairs;?? ????vector<CvPoint2D32f>?pt1,?pt2;?? ????CvMat?_pt1,?_pt2;?? ????int?i,?n;#ifdef?USE_FLANN?? ????flannFindPairs(?objectKeypoints,?objectDescriptors,?imageKeypoints,?imageDescriptors,?ptpairs?);?? #else?? ????findPairs(?objectKeypoints,?objectDescriptors,?imageKeypoints,?imageDescriptors,?ptpairs?);?? #endif????n?=?(int)(ptpairs.size()/2);?? ????if(?n?<?4?)?? ????????return?0;????pt1.resize(n);?? ????pt2.resize(n);?? ????for(?i?=?0;?i?<?n;?i++?)?? ????{?? ????????pt1[i]?=?((CvSURFPoint*)cvGetSeqElem(objectKeypoints,ptpairs[i*2]))->pt;?? ????????pt2[i]?=?((CvSURFPoint*)cvGetSeqElem(imageKeypoints,ptpairs[i*2+1]))->pt;?? ????}????_pt1?=?cvMat(1,?n,?CV_32FC2,?&pt1[0]?);?? ????_pt2?=?cvMat(1,?n,?CV_32FC2,?&pt2[0]?);?? ????if(?!cvFindHomography(?&_pt1,?&_pt2,?&_h,?CV_RANSAC,?5?))?? ????????return?0;????for(?i?=?0;?i?<?4;?i++?)?? ????{?? ????????double?x?=?src_corners[i].x,?y?=?src_corners[i].y;?? ????????double?Z?=?1./(h[6]*x?+?h[7]*y?+?h[8]);?? ????????double?X?=?(h[0]*x?+?h[1]*y?+?h[2])*Z;?? ????????double?Y?=?(h[3]*x?+?h[4]*y?+?h[5])*Z;?? ????????dst_corners[i]?=?cvPoint(cvRound(X),?cvRound(Y));?? ????}????return?1;?? }?? ?? int?main(int?argc,?char**?argv)?? {?? ??? ????const?char*?object_filename?=?argc?==?3???argv[1]?:?"D:/src.jpg";?? ????const?char*?scene_filename?=?argc?==?3???argv[2]?:?"D:/Demo.jpg";????help();????IplImage*?object?=?cvLoadImage(?object_filename,?CV_LOAD_IMAGE_GRAYSCALE?);?? ????IplImage*?image?=?cvLoadImage(?scene_filename,?CV_LOAD_IMAGE_GRAYSCALE?);?? ????if(?!object?||?!image?)?? ????{?? ????????fprintf(?stderr,?"Can?not?load?%s?and/or?%s\n",?? ????????????object_filename,?scene_filename?);?? ????????exit(-1);?? ????}?? ??? ????CvMemStorage*?storage?=?cvCreateMemStorage(0);????cvNamedWindow("Object",?1);?? ????cvNamedWindow("Object?Correspond",?1);????static?CvScalar?colors[]?=??? ????{?? ????????{{0,0,255}},?? ????????{{0,128,255}},?? ????????{{0,255,255}},?? ????????{{0,255,0}},?? ????????{{255,128,0}},?? ????????{{255,255,0}},?? ????????{{255,0,0}},?? ????????{{255,0,255}},?? ????????{{255,255,255}}?? ????};?? ???? ????IplImage*?object_color?=?cvCreateImage(cvGetSize(object),?8,?3);?? ????cvCvtColor(?object,?object_color,?CV_GRAY2BGR?);????CvSeq*?objectKeypoints?=?0,?*objectDescriptors?=?0;?? ????CvSeq*?imageKeypoints?=?0,?*imageDescriptors?=?0;?? ????int?i;?? ?? ? ? ? ? ? ? ? ? ? ? ? ?? ?CvSURFParams?params?=?cvSURFParams(500,?1);????double?tt?=?(double)cvGetTickCount();?? ?? ? ? ? ? ? ?? ????cvExtractSURF(?object,?0,?&objectKeypoints,?&objectDescriptors,?storage,?params?);?? ????printf("Object?Descriptors:?%d\n",?objectDescriptors->total);????cvExtractSURF(?image,?0,?&imageKeypoints,?&imageDescriptors,?storage,?params?);?? ????printf("Image?Descriptors:?%d\n",?imageDescriptors->total);?? ????tt?=?(double)cvGetTickCount()?-?tt;????printf(?"Extraction?time?=?%gms\n",?tt/(cvGetTickFrequency()*1000.));?? ????CvPoint?src_corners[4]?=?{{0,0},?{object->width,0},?{object->width,?object->height},?{0,?object->height}};?? ?????? ?CvPoint?dst_corners[4];?? ????IplImage*?correspond?=?cvCreateImage(?cvSize(image->width,?object->height+image->height),?8,?1?);?? ??? ??? ????cvSetImageROI(?correspond,?cvRect(?0,?0,?object->width,?object->height?)?);?? ????cvCopy(?object,?correspond?);?? ????cvSetImageROI(?correspond,?cvRect(?0,?object->height,?correspond->width,?correspond->height?)?);?? ????cvCopy(?image,?correspond?);?? ????cvResetImageROI(?correspond?);#ifdef?USE_FLANN?? ????printf("Using?approximate?nearest?neighbor?search\n");?? #endif?? ??? ????if(?locatePlanarObject(?objectKeypoints,?objectDescriptors,?imageKeypoints,?? ????????imageDescriptors,?src_corners,?dst_corners?))?? ????{?? ????????for(?i?=?0;?i?<?4;?i++?)?? ????????{?? ????????????CvPoint?r1?=?dst_corners[i%4];?? ????????????CvPoint?r2?=?dst_corners[(i+1)%4];?? ???cvLine(?correspond,?cvPoint(r1.x,?r1.y+object->height?),?? ????cvPoint(r2.x,?r2.y+object->height?),?colors[8]?);?? ????????}?? ????}?? ??? ??? ????vector<int>?ptpairs;?? #ifdef?USE_FLANN?? ????flannFindPairs(?objectKeypoints,?objectDescriptors,?imageKeypoints,?imageDescriptors,?ptpairs?);?? #else?? ????findPairs(?objectKeypoints,?objectDescriptors,?imageKeypoints,?imageDescriptors,?ptpairs?);?? #endif?? ??? ????for(?i?=?0;?i?<?(int)ptpairs.size();?i?+=?2?)?? ????{?? ????????CvSURFPoint*?r1?=?(CvSURFPoint*)cvGetSeqElem(?objectKeypoints,?ptpairs[i]?);?? ????????CvSURFPoint*?r2?=?(CvSURFPoint*)cvGetSeqElem(?imageKeypoints,?ptpairs[i+1]?);?? ????????cvLine(?correspond,?cvPointFrom32f(r1->pt),?? ????????????cvPoint(cvRound(r2->pt.x),?cvRound(r2->pt.y+object->height)),?colors[8]?);?? ????}????cvShowImage(?"Object?Correspond",?correspond?);?? ??? ????for(?i?=?0;?i?<?objectKeypoints->total;?i++?)?? ????{?? ????????CvSURFPoint*?r?=?(CvSURFPoint*)cvGetSeqElem(?objectKeypoints,?i?);?? ????????CvPoint?center;?? ????????int?radius;?? ????????center.x?=?cvRound(r->pt.x);?? ????????center.y?=?cvRound(r->pt.y);?? ????????radius?=?cvRound(r->size*1.2/9.*2);?? ????????cvCircle(?object_color,?center,?radius,?colors[0],?1,?8,?0?);?? ????}?? ????cvShowImage(?"Object",?object_color?);????cvWaitKey(0);??? ????cvDestroyWindow("Object");?? ????cvDestroyWindow("Object?Correspond");????return?0;?? }?? ????
?
?
?
from: http://blog.csdn.net/yangtrees/article/details/7482960
《新程序員》：云原生和全面數字化實踐50位技術專家共同創作，文字、視頻、音頻交互閱讀
總結
以上是生活随笔為你收集整理的学习OpenCV——Surf（特征点篇）flann快速最近邻搜索算法的全部內容，希望文章能夠幫你解決所遇到的問題。
如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。
