SIFT算法實(shí)現(xiàn)及代碼詳解

經(jīng)典算法SIFT實(shí)現(xiàn)即代碼解釋：以下便是sift源碼庫編譯后的效果圖:14-.XL14-.XLTT訃：：為了給有興趣實(shí)現(xiàn)Sift算法的朋友提供個參考，特整理此文如下。要了解什么是sift算法，請參考：九、圖像特征提取與匹配之SIFT算法。ok,咱們下面,就來利用RobHess維護(hù)的sift庫來實(shí)現(xiàn)sift算法：首先，請下載RobHess維護(hù)的sift庫：/hess/code/sift/下載RobHess的這個壓縮包后，如果直接解壓縮，直接編譯，那么會出現(xiàn)下面的錯誤提示：編譯提示:errorC1083:Cannotopenincludefile:'cxcore.h':Nosuchfileordirectory，找不到這個頭文件。這個錯誤,是因為你還沒有安裝opencv，因為:cxcore.h和cv.h是開源的OPENCV頭文件，不是VC++的默認(rèn)安裝文件，所以你還得下載OpenCV并進(jìn)行安裝。然后，可以在OpenCV文件夾下找到你所需要的頭文件了。據(jù)網(wǎng)友稱，截止2010年4月4日，還沒有在VC6.0下成功使用opencv2.0的案例。所以，如果你是VC6.0的用戶請下載opencv1.0版本°vs的話,opencv2.0,1.0任意下載。以下，咱們就以vc6.0為平臺舉例，下載并安裝opencv1.0版本、gsl等。當(dāng)然,你也可以用vs編譯，同樣下載opencv（具體版本不受限制）、gsl等。請按以下步驟操作：一、下載opencv1.0/projects/opencvlibrary/files/opencv-win/1.O/OpenCV_1.O.exe/download二、安裝opencvl.O，配置Windows環(huán)境變量1、安裝注意：假如你是將OpenCV安裝到C:/ProgramFiles/OpenCV（如果你安裝的時候選擇不是安裝在C盤，則下面所有對應(yīng)的C盤都改為你所安裝在的那個“X盤”，即可），在安裝時選擇"將/OpenCV/bin加入系統(tǒng)變量"，打上“勾”。（Add/OpenCV/bintothesystermPATH。這一步確認(rèn)選上了之后，下面的檢查環(huán)境變量的步驟，便可免去）2、 檢查環(huán)境變量。為了確保上述步驟中，加入了系統(tǒng)變量，在安裝opencv1.0成功后，還得檢查C:/ProgramFiles/OpenCV/bin是否已經(jīng)被加入到環(huán)境變量PATH如果沒有，請加入。3、 最后是配置VisualC++6.0。全局設(shè)置菜單Tools->Options->Directories：先設(shè)置lib路徑，選擇Libraryfiles，在下方填入路徑：C:/ProgramFiles/OpenCV/lib然后選擇includefiles，在下方填入路徑（參考下圖）：C:/ProgramFiles/OpenCV/cxcore/includeC:/ProgramFiles/OpenCV/cv/includeC:/ProgramFiles/OpenCV/cvaux/includeC:/ProgramFiles/OpenCV/ml/includeC:/ProgramFiles/OpenCV/otherlibs/highguiC:/ProgramFiles/OpenCV/otherlibs/cvcam/include

Editor|Tabs|Debug|Compatibility|BuildOptiv&sDirectoriesWorks]I~&Platform;|Win32Platform;|Win32|lncludefilesDirectories：C:\ProqramFiles^MicrosoftVisualStudio\VC98\INCLUDECAProqramFiles^MicrosoftVisualStudio\VC98\MFC\INCLUDECAProgramFiles^MicrosoftVisualStudio\VC90\ATmNCLUDEC:\PROGRAMFILES\OPENCV\CXCORE\INCLUDEC:\PROGRAMFILES\OPENCV\CV\INCLUDEC:\PROGRAMFILES\OPENCV\MmNCLUDEC:\PROGRAMFILES\OPENCV\CVAUX\INCLUDEC:\PROGRAMFILES\OPENCV\OTHERLIBS\HIGHGUIC:\PROGRAMFILESPFENCWQTHERL舊S'CYCA兇INCLUDEOK Cancel最后選擇sourcefiles，在下方填入路徑：C:/ProgramFiles/OpenCV/cv/srcC:/ProgramFiles/OpenCV/cxcore/srcC:/ProgramFiles/OpenCV/cvaux/srcC:/ProgramFiles/OpenCV/otherlibs/highguiC:/ProgramFiles/OpenCV/otherlibs/cvcam/src/windows項目設(shè)置每創(chuàng)建一個將要使用OpenCV的VCProject,都需要給它指定需要的lib。菜單：Project->Settings，然后將Settingfor選為AllConfigurations，然后選擇右邊的link標(biāo)簽，在Object/librarymodules附加上：cxcore.libcv.libml.libcvaux.libhighgui.libcvcam.lib當(dāng)然，你不需要這么多l(xiāng)ib,你可以只添加你需要的lib（見下圖）三、 下載gsl,gsl也是一個庫，也需要下載：/projects/gnuwin32/files/gsl/1.8/gsl-1.8.exe/download。在編譯時候GSL也是和OpenCV—樣要把頭文件和lib的路徑指定好。四、 配置gsl將C:/WinGsl/bin中的WinGsl.dll和WinGslD.dll復(fù)制到C:/VC6.0/Bin；將整個Gsl目錄復(fù)制到C:/VC6.0/Bin下；lib目錄下的所有.lib文件全部復(fù)制到C:/VC6.0/Lib下。然后，在tools—options—directories中，將C:/WinGsl下的lib,gsl分別加入到庫文件和頭文件的搜索路徑中。以下是可能會出現(xiàn)的錯誤情況處理：I、 OpenCV安裝后“沒有找到cxcorelOO.dll，啲錯誤處理在安裝時選擇“將/OpenCV/bin加入系統(tǒng)變量"(Add/OpenCV/bintothesystermPATH)。但該選項并不一定能成功添加到系統(tǒng)變量，如果編寫的程序在運(yùn)行時出現(xiàn)“沒有找到cxcore100.dll,因為這個應(yīng)用程序未能啟動。重新安裝應(yīng)用程序可能會修復(fù)此問題?！钡腻e誤。手動在我的電腦-＞屬性-＞高級-＞環(huán)境變量-＞系統(tǒng)變量-＞path添加c:/programfiles/opencv/bin;添加完成后需要重啟計算機(jī)。II、 vc6.0下配置了一下，可是編譯程序時遇到如下一個錯誤：Linking...LINK:fatalerrorLNK1104:cannotopenfile"odbccp32.libcxcore.lib"可能是：在工程設(shè)置的時候添加連接庫時沒加空格或 ?來把兩個文件名(odbccp32.libcxcore.lib)分開。注意每一次操作后，記得保存。若經(jīng)過以上所有的步驟之后，如果還不能正常編譯，那就是還要稍微修改下你下載的RobHess代碼。ok，日后，若有空，再好好詳細(xì)剖析下此sift的源碼。最后，祝你編譯順利。宀完。SIFT代碼詳解：這是一個很強(qiáng)大的算法，主要用于圖像配準(zhǔn)和物體識別等領(lǐng)域，但是其計算量相比也比較大，性價比比較高的算法包括PCA-SIFT和SURF其中OpenCV提供了SURF算法，但是為了方便理解。這里給出了RobHess所實(shí)現(xiàn)的SIFT算法的實(shí)現(xiàn)以及注釋，結(jié)合我自己的理解，如果，您有關(guān)于SIFT算法不理解的地方咱們可以一起交流一下?；蛘吣J(rèn)為不詳細(xì)的地方提出來。SIFT算法的主要實(shí)現(xiàn)在sift.c這個文件，其主要流程為：首先創(chuàng)建初始圖像，即通過將圖像轉(zhuǎn)換為32位的灰度圖，然后將圖像使用三次插值來方大，之后通過高斯模糊處理在此基礎(chǔ)上進(jìn)行高斯金字塔的構(gòu)建以及高斯差分金字塔的構(gòu)建對圖像進(jìn)行極值點(diǎn)檢測計算特征向量的尺度調(diào)整圖像大小計算特征的方向計算描述子，其中包括計算二維方向直方圖并轉(zhuǎn)換直方圖為特征描述子首先給出sift算法的整體框架代碼：輸入?yún)?shù)：img為輸入圖像；feat為所要提取的特征指針；intvl指的是高斯金字塔和差分金字塔的層數(shù)；sigma指的是圖像初始化過程中高斯模糊所使用的參數(shù)；contr_thr是歸一化之后的去除不穩(wěn)定特征的閾值；curv_thr指的是去除邊緣的特征的主曲率閾值；img_dbl是是否將圖像放大為之前的兩倍；descr_with用來計算特征描述子的方向直方圖的寬度；descr_hist_bins是直方圖中的條數(shù)[cpp]viewplaincopyint_sift_features(IplImage*img,structfeature**feat,intintvls,doublesigma,doublecontr_thr,intcurv_thr,intimg_dbl,intdescr_width,intdescr_hist_bins){IplImage*init_img;IplImage***gauss_pyr,***dog_pyr;CvMemStorage*storage;CvSeq*features;intoctvs,i,n=0;10.11. /*checkarguments*/4.if(!img)fatal_error("NULLpointererror,%s,line%d",__FILE__,__LINE__);if(!feat)fatal_error("NULLpointererror,%s,line%d",__FILE__,__LINE__);/*buildscalespacepyramid;smallestdimensionoftoplevelis~4pixels*//*構(gòu)建高斯尺度空間金字塔，頂層最小的為4像素*/init_img=create_init_img(img,img_dbl,sigma);octvs=log(doubleMIN(init_img->width,init_img->height))/log(2.0)-2;//構(gòu)建高斯金字塔和高斯差分金字塔gauss_pyr=build_gauss_pyr(init_img,octvs,intvls,sigma);dog_pyr=build_dog_pyr(gauss_pyr,octvs,intvls);storage=cvCreateMemStorage(0);//尺度空間極值點(diǎn)檢測features=scale_space_extrema(dog_pyr,octvs,intvls,contr_thr,curv_thr,storage);//畫出去除低對比度的極值點(diǎn)//draw_extrempoint(img,features);//計算特征向量的尺度calc_feature_scales(features,sigma,intvls);if(img_dbl)adjust_for_img_dbl(features);//計算特征的方向calc_feature_oris(features,gauss_pyr);//計算描述子，包括計算二維方向直方圖和轉(zhuǎn)換其為特征描述子compute_descriptors(features,gauss_pyr,descr_width,descr_hist_bins);/*sortfeaturesbydecreasingscaleandmovefromCvSeqtoarray*/cvSeqSort(features,(CvCmpFunc)feature_cmp,NULL);n=features->total;*feat=static_cast<feature*>(calloc(n,sizeof(structfeature)));*feat=static_cast<feature*>(cvCvtSeqToArray(features,*feat,CV_WHOLE_SEQ))8.59.60.for(i=0;i<n;i++)8.59.60.for(i=0;i<n;i++){free((*feat)[i].feature_data);(*feat)[i].feature_data=NULL;}2.{.2.61.cvReleaseMemStorage(&storage);cvReleaseImage(&init_img);release_pyr(&gauss_pyr,octvs,intvls+3);release_pyr(&dog_pyr,octvs,intvls+2);returnn;}(1)初始化圖像輸入?yún)?shù)：這里不需要解釋了該函數(shù)主要用來初始化圖像，轉(zhuǎn)換圖像為32位灰度圖以及進(jìn)行高斯模糊。[cpp]viewplaincopystaticIplImage*create_init_img(IplImage*img,intimg_dbl,doublesigma)IplImage*gray,*dbl;floatsig_diff;gray=convert_to_gray32(img);if(img_dbl){sig_diff=sqrt(sigma*sigma-SIFT_INIT_SIGMA*SIFT_INIT_SIGMA*4);dbl=cvCreateImage(cvSize(img->width*2,img->height*2),IPL_DEPTH_32F,1);cvResize(gray,dbl,CV_INTER_CUBIC);cvSmooth(dbl,dbl,CV_GAUSSIAN,0,0,sig_diff,sig_diff);cvReleaseImage(&gray);returndbl;}else{sig_diff=sqrt(sigma*sigma-SIFT_INIT_SIGMA*SIFT_INIT_SIGMA);cvSmooth(gray,gray,CV_GAUSSIAN,0,0,sig_diff,sig_diff);returngray;}

構(gòu)建高斯金字塔輸入?yún)?shù)：octvs是高斯金字塔的組invls是高斯金字塔的層數(shù)sigma是初始的高斯模糊參數(shù)，后續(xù)也通過它計算每一層所使用的sigma[cpp]viewplaincopy<spanstyle="font-size:13px;">staticIplImage***build_gauss_pyr(IplImage*base,intoctvs,intintvls,doublesigma){IplImage***gauss_pyr;double*sig=static_cast<double*>(calloc(intvls+3,sizeof(double)));doublesig_total,sig_prev, k;inti,o;7.gauss_pyr=static_cast<IplImage***>(calloc(octvs,sizeof(IplImage**)));for(i=0;i<octvs;i++)gauss_pyr[i]=static_cast<IplImage**>(calloc(intvls+3,sizeof(IplImage*)));9.30.precomputeGaussiansigmasusingthefollowingformula:預(yù)計算每次高斯模糊的sigma\sigma_{total}A2=\sigma_{i}A2+\sigma_{i-1}A2*/sig[0]=sigma;k=pow(2.0,1.0/intvls);for(i=1;i<intvls+3;i++){sig_prev=pow(k,i-1)*sigma;sig_total=sig_prev*k;sig[i]=sqrt(sig_total*sig_total-sig_prev*sig_prev);}for(o=0;o<octvs;o++)for(i=0;i<intvls+3;i++){//對每一層進(jìn)行降采樣，形成高斯金字塔的每一層if(o==0 &&i==0)

34.35./*baseofnewoctvaveishalvedimagefromendofpreviousoctave*/36.//每一組的第一層都是通過對前面一組的第一層降采樣實(shí)現(xiàn)的37.elseif(i==0)38.gauss_pyr[o][i]=downsample(gauss_pyr[o-1][intvls]);39.40./*blurthecurrentoctave'slastimagetocreatethenextone*/41.//每一組的其他層則使通過使用不同sigma的高斯模糊來進(jìn)行處理42.else43.44.gauss_pyr[o][i]=cvCreateImage(cvGetSize(gauss_pyr[o][i-1]),45.IPL_DEPTH_32F,1);46.cvSmooth(gauss_pyr[o][i-1],gauss_pyr[o][i],47.CV_GAUSSIAN,0,0,sig[i],sig[i]);48.34.35./*baseofnewoctvaveishalvedimagefromendofpreviousoctave*/36.//每一組的第一層都是通過對前面一組的第一層降采樣實(shí)現(xiàn)的37.elseif(i==0)38.gauss_pyr[o][i]=downsample(gauss_pyr[o-1][intvls]);39.40./*blurthecurrentoctave'slastimagetocreatethenextone*/41.//每一組的其他層則使通過使用不同sigma的高斯模糊來進(jìn)行處理42.else43.44.gauss_pyr[o][i]=cvCreateImage(cvGetSize(gauss_pyr[o][i-1]),45.IPL_DEPTH_32F,1);46.cvSmooth(gauss_pyr[o][i-1],gauss_pyr[o][i],47.CV_GAUSSIAN,0,0,sig[i],sig[i]);1.free(sig);52.returngauss_pyr;53.}</span>降采樣處理輸入?yún)?shù)：不解釋這就是降采樣，其實(shí)就是將圖像通過最近鄰算法縮小為原來的一半[cpp]viewplaincopystaticIplImage*downsample(IplImage*img){IplImage*smaller=cvCreateImage(cvSize(img->width/2,img->height/2),img->depth,img->nChannels);cvResize(img,smaller,CV_INTER_NN);6.returnsmaller;}構(gòu)建高斯差分金字塔輸入?yún)?shù):不解釋了參見上面的說明即可實(shí)際上差分金字塔的構(gòu)成是通過對相鄰層的圖像進(jìn)行相減獲得的[cpp]viewplaincopy<spanstyle="font-size:16px;">staticIplImage***build_dog_pyr(IplImage***gauss_pyr,intoctvs,intintvls){IplImage***dog_pyr;inti,o;5.dog_pyr=static_cast<IplImage***>(calloc(octvs,sizeof(IplImage**)));for(i=0;i<octvs;i++)dog_pyr[i]=static_cast<IplImage**>(calloc(intvls+2,sizeof(IplImage*)));9.for(o=0;o<octvs;o++)for(i=0;i<intvls+2;i++){dog_pyr[o][i]=cvCreateImage(cvGetSize(gauss_pyr[o][i]),IPL_DEPTH_32F,1);cvSub(gauss_pyr[o][i+1],gauss_pyr[o][i],dog_pyr[o][i],NULL);}17.returndog_pyr;}</span>(4)極值點(diǎn)檢測輸入?yún)?shù)：contr_thr是去除對比度低的點(diǎn)所采用的閾值curv_thr是去除邊緣特征的閾值[cpp]viewplaincopy1?staticCvSeq*scale_space_extrema(IplImage***dog_pyr,intoctvs,intintvls,2?doublecontr_thr,intcurv_thr,3?CvMemStorage*storage)4.{5?CvSeq*features;

.10.doubleprelim_contr_thr=0.5*contr_thr/intvls;structfeature*feat;structdetection_data*ddata;into,i,r,c;11.features=cvCreateSeq(0,sizeof(CvSeq),sizeof(structfeature),storage);12.for(o=0;o<octvs;o++)13.for(i=1;i<=intvls;i++)14.for(r=SIFT_IMG_BORDER;r<dog_pyr[o][0]->height-SIFT_IMG_BORDER;r++)15.for(c=SIFT_IMG_BORDER;c<dog_pyr[o][0]->width-SIFT_IMG_BORDER；c++16./*performpreliminarycheckoncontrast*/17.if(ABS(pixval32f(dog_pyr[o][i],r,c))>prelim_contr_thr)18.if.10.doubleprelim_contr_thr=0.5*contr_thr/intvls;structfeature*feat;structdetection_data*ddata;into,i,r,c;11.features=cvCreateSeq(0,sizeof(CvSeq),sizeof(structfeature),storage);12.for(o=0;o<octvs;o++)13.for(i=1;i<=intvls;i++)14.for(r=SIFT_IMG_BORDER;r<dog_pyr[o][0]->height-SIFT_IMG_BORDER;r++)15.for(c=SIFT_IMG_BORDER;c<dog_pyr[o][0]->width-SIFT_IMG_BORDER；c++16./*performpreliminarycheckoncontrast*/17.if(ABS(pixval32f(dog_pyr[o][i],r,c))>prelim_contr_thr)18.if(is_extremum(dog_pyr,o,i,r,c))19.20.feat=interp_extremum(dog_pyr,o,i,r,c,intvls,contr_0.31.thr);tvl],if(feat)ddata=feat_detection_data(feat);if(!is_too_edge_like(dog_pyr[ddata->octv][ddata->inddata->r,ddata->c,curv_thr))cvSeqPush(features,feat);elsefree(ddata);free(feat);32.33.34.returnfeatures;}SIFT_IMG_BORDER是預(yù)定義的圖像邊緣；通過和對比度閾值比較去掉低對比度的點(diǎn)；而通過is_extremum來判斷是否為極值點(diǎn)，如果是則通過極值點(diǎn)插值的方式獲取亞像素的極值點(diǎn)的位置。然后通過is_too_eage_like和所給的主曲率閾值判斷是否為邊緣點(diǎn)判斷是否為極值點(diǎn)

其原理為:通過和高斯金字塔的上一層的9個像素+本層的除了本像素自己的其他的8個像素和下一層的9個像素進(jìn)行比較看是否為這26個像素中最小的一個或者是否為最大的一個，如果是則為極值點(diǎn)。[cpp]viewplaincopy1.staticintis_extremum(IplImage***dog_pyr,intoctv,intintvl,intr,intc)2.{3.floatval=pixval32f(dog_pyr[octv][intvl],r,c);4.inti,j,k;5.6./*checkformaximum*/7.if(val>0)8.{9.for(i=-1;i<=1;i++)10.for(j=-1;j<=1;j++)11.for(k=-1;k<=1;k++)12.if(val<pixval32f(dog_pyr[octv][intvl+i],r+j,c+k))13.return0;14.}15.16./*checkforminimum*/17.else18.{19.for(i=-1;i<=1;i++)20.for(j=-1;j<=1;j++)21.for(k=-1;k<=1;k++)22.if(val>pixval32f(dog_pyr[octv][intvl+i],r+j,c+k))23.return0;24.}25.26.return1;27.}*獲取亞像素的極值點(diǎn)的位置[cpp]viewplaincopyoctv,intintvl,staticstructoctv,intintvl,2.intr,intc,intintvls,doublecontr_thr)

2.{4.5?4.5?.3.44.45.structdetection_data*ddata;doublexi,xr,xc,contr;//分別為亞像素的intval,row,col的偏移offset,和對比度inti=0;while(i<SIFT_MAX_INTERP_STEPS)//重新確定極值點(diǎn)并重新定位的操作只能循環(huán)5次{interp_step(dog_pyr,octv,intvl,r,c,&xi,&xr,&xc);if(ABS(xi)<0.5&&ABS(xr)<0.5&&ABS(xc)<0.5)//如果滿足條件就停止尋找break;//否則繼續(xù)尋找極值點(diǎn)c+=cvRound(xc);r+=cvRound(xr);intvl+=cvRound(xi);if(intvl<1 ||intvl>intvls||c<SIFT_IMG_BORDER||r<SIFT_IMG_BORDER||c>=dog_pyr[octv][0]->width-SIFT_IMG_BORDER||r>=dog_pyr[octv][0]->height-SIFT_IMG_BORDER){returnNULL;}i++;}//確保極值點(diǎn)是經(jīng)過最大5步找到的/*ensureconvergenceofinterpolation*/if(i>=SIFT_MAX_INTERP_STEPS)returnNULL;//獲取找到的極值點(diǎn)的對比度contr=interp_contr(dog_pyr,octv,intvl,r,c,xi,xr,xc);//判斷極值點(diǎn)是否小于某一個閾值if(ABS(contr)<contr_thr/intvls)returnNULL;//若小于，則認(rèn)為是極值點(diǎn)feat=new_feature();ddata=feat_detection_data(feat);feat->img_pt.x=feat->x=(c+xc)*pow(2.0,octv);

46.feat->img_pt.y=feat->y=(r+xr)*pow(2.0,octv);47.ddata->r=r;48.ddata->c=c;49.ddata->octv=octv;50.ddata->intvl=intvl;51.ddata->subintvl=xi;52.53.returnfeat;54.}*獲取亞像素位置中所用到的函數(shù)[cpp]viewplaincopy1.staticvoidinterp_step(IplImage***dog_pyr,intoctv,intintvl,intr,intc,2?double*xi,double*xr,double*xc)3.{4.CvMat*dD,*H,*H_inv,X;5?doublex[3]={0};6?7.8?//計算三維偏導(dǎo)數(shù)dD=deriv_3D(dog_pyr,octv,intvl,r,c);//計算三維海森矩陣H=hessian_3D(dog_pyr,octv,intvl,r,c);H_inv=cvCreateMat(3,3,CV_64FC1);cvInvert(H,H_inv,CV_SVD);cvInitMatHeader(&X,3,1,CV_64FC1,x,CV_AUTOSTEP);15.cvGEMM(H_inv,dD,-1,NULL,0,&X,0)；17.cvReleaseMat(&dD);cvReleaseMat(&H);cvReleaseMat( &H_inv);21.TOC\o"1-5"\h\z*xi = x[2];*xr = x[1];*xc = x[0];}計算三維偏導(dǎo)數(shù)

計算在x和y方向上的偏導(dǎo)數(shù)，高斯差分尺度空間金字塔中像素的尺度實(shí)際上在離散數(shù)據(jù)中計算偏導(dǎo)數(shù)是通過相鄰像素的相減來計算的比如說計算x方向的偏導(dǎo)數(shù)dx,則通過該向所的x方向的后一個減去前一個然后除以2即可求的dx[cpp]viewplaincopystaticCvMat*deriv_3D(IplImage***dog_pyr,intoctv,intintvl,intr,intc){CvMat*dI;doubledx,dy,ds;5.dx=(pixval32f(dog_pyr[octv][intvl],r,c+1)-pixval32f(dog_pyr[octv][intvl],r,c-1))/2.0;dy=(pixval32f(dog_pyr[octv][intvl],r+1,c)-pixval32f(dog_pyr[octv][intvl],r-1,c))/2.0;ds=(pixval32f(dog_pyr[octv][intvl+1],r,c)-pixval32f(dog_pyr[octv][intvl-1],r,c))/2.0;12.13?dI=cvCreateMat(3,1,CV_64FC1);14.cvmSet(di,0,0,dx);15.cvmSet(di,1,0,dy);16.cvmSet(di,2,0,ds);17.18.returndi;19.}*計算三維海森矩陣不需要講什么，其實(shí)就是計算二次導(dǎo)數(shù)，計算方法也和一次導(dǎo)數(shù)的計算如出一轍。然后將結(jié)果放入到一個矩陣中去。[cpp]viewplaincopystaticCvMat*hessian_3D(邛lImage***dog_pyr,intoctv,intintvl,intr,intc){CvMat*H;doublev,dxx,dyy,dss,dxy,dxs,dys;

5.v=pixval32f(dog_pyr[octv][intvl],r,c);7.dxx=(pixval32f(dog_pyr[octv][intvl],r,c+1)+8.pixval32f(dog_pyr[octv][intvl],r,c-1)-2*v);9.dyy=(pixval32f(dog_pyr[octv][intvl],r+1,c)+10.pixval32f(dog_pyr[octv][intvl],r-1,c)-2*v);11.dss=(pixval32f(dog_pyr[octv][intvl+1],r,c)+12.pixval32f(dog_pyr[octv][intvl-1],r,c)-2*v);13.dxy=(pixval32f(dog_pyr[octv][intvl],r+1,c+1)-14.pixval32f(dog_pyr[octv][intvl],r+1,c-1)-15.pixval32f(dog_pyr[octv][intvl],r-1,c+1)+16.pixval32f(dog_pyr[octv][intvl],r-1,c-1))/4.0;17.dxs=(pixval32f(dog_pyr[octv][intvl+1],r,c+1)-18.pixval32f(dog_pyr[octv][intvl+1],r,c-1)-19.pixval32f(dog_pyr[octv][intvl-1],r,c+1)+20.pixval32f(dog_pyr[octv][intvl-1],r,c-1))/4.0;21.dys=(pixval32f(dog_pyr[octv][intvl+1],r+1,c)-22.pixval32f(dog_pyr[octv][intvl+1],r-1,c)-23.pixval32f(dog_pyr[octv][intvl-1],r+1,c)+24.pixval32f(dog_pyr[octv][intvl-1],r-1,c))/4.0;25.26.H=cvCreateMat(3,3,CV_64FC1);27.cvmSet(H,0,0,dxx);28.cvmSet(H,0,1,dxy);29.cvmSet(H,0,2,dxs);30.cvmSet(H,1,0,dxy);31.cvmSet(H,1,1,dyy);32.cvmSet(H,1,2,dys);33.cvmSet(H,2,0,dxs);34.cvmSet(H,2,1,dys);35.cvmSet(H,2,2,dss);36.returnH;}*計算插入像素的對比度[cpp]viewplaincopy..5.staticdoubleinterp_contr(IplImage***dog_pyr,intoctv,intintvl,intr,intc,doublexi,doublexr,doublexc){CvMat*dD,X,T;doublet[1],x[3]={xc,xr,xi};6.TOC\o"1-5"\h\zcvInitMatHeader(&X,3, 1,CV_64FC1, x, CV_AUTOSTEP );cvInitMatHeader(&T,1, 1,CV_64FC1, t, CV_AUTOSTEP );dD=deriv_3D(dog_pyr, octv,intvl, r, c);cvGEMM(dD,&X,1,NULL,0,&T,CV_GEMM_A_T)；cvReleaseMat(&dD);12.returnpixval32f(dog_pyr[octv][intvl],r,c)+t[0]*0.5;}其中cvGEMM是矩陣的通用計算函數(shù)，至于CV_GEMM_A_T是計算dD的轉(zhuǎn)置矩陣放入T中*去除邊緣相應(yīng)通過計算所在特征向量的主曲率半徑來判斷特征是邊緣的從而導(dǎo)致不穩(wěn)定即去除邊緣響應(yīng)[cpp]viewplaincopystaticintis_too_edge_like(IplImage*dog_img,intr,intc,intcurv_thr)TOC\o"1-5"\h\z{doubled,dxx,dyy,dxy,tr,det;4./*principalcurvaturesarecomputedusingthetraceanddetofHessian*/d= pixval32f(dog_img,r, c);dxx = pixval32f(dog_img, r,c+1 ) + pixval32f( dog_img, r,c-1 ) - 2 * d;dyy = pixval32f(dog_img, r+1,c ) + pixval32f( dog_img, r-1,c ) - 2 * d;dxy = (pixval32f(dog_img,r+1,c+1) -pixval32f(dog_img,r+1,c-1)-pixval32f(dog_img,r-1,c+1)+pixval32f(dog_img,r-1,c-1))/4.0;tr=dxx+dyy;det=dxx*dyy-dxy*dxy;13./*negativedeterminant->curvatureshavedifferentsigns;rejectfeature*/if(det<=0)return1;17.if(tr*tr/det<(curv_thr+1.0)*(curv_thr+1.0)/curv_thr)return0;return1;}計算特征向量的尺度

實(shí)際上是通過最初的sigma來獲得每一層每一組的尺度[cpp]viewplaincopystaticvoidcalc_feature_scales(CvSeq*features,doublesigma,intintvls)TOC\o"1-5"\h\z{struct feature* feat;struct detection_data*ddata;double intvl;inti, n;7.n=features->total;for(i=0;i<n;i++){feat=CV_GET_SEQ_ELEM(structfeature,features,i);ddata=feat_detection_data(feat);intvl=ddata->intvl+ddata->subintvl;feat->scl=sigma*pow(2.0,ddata->octv+intvl/intvls);ddata->scl_octv=sigma*pow(2.0,intvl/intvls);}}調(diào)整圖像特征坐標(biāo)、尺度、點(diǎn)的坐標(biāo)大小為原來的一半[java]viewplaincopystaticvoidadjust_for_img_dbl(CvSeq*features){structfeature*feat;inti,n;.3.14.15.}for(i=0;i<n;i++){feat=CV_GET_SEQ_ELEM(structfeature,features,i);feat->x/=2.0;feat->y/=2.0;feat->scl/=2.0;feat->img_pt.x/=2.0;feat->img_pt.y/=2.0;}給每一個圖像特征向量計算規(guī)范化的方向[cpp]viewplaincopystaticvoidcalc_feature_oris(CvSeq*features,IplImage***gauss_pyr){struct feature* feat;struct detection_data*ddata;double*hist;double omax;inti,j,n=features->total;8.9.//遍歷整個檢測出來的特征點(diǎn)，計算每個特征點(diǎn)的直方圖，然后平滑直方圖去除突變，然后找到每一個特征點(diǎn)的主方向，并加入到好的方向特征數(shù)組中去11.for(:i=0;i<n;i++)12.{13.feat=static_cast<feature*>(malloc(sizeof(structfeature)));14.cvSeqPopFront(features,feat);15.ddata=feat_detection_data(feat);16.//計算給定的某個像素的灰度方向直方圖17.hist=ori_hist(gauss_pyr[ddata->octv][ddata->intvl],18.ddata->r,ddata->c,SIFT_ORI_HIST_BINS,19.cvRound(SIFT_ORI_RADIUS*ddata->scl_octv),20.SIFT_ORI_SIG_FCTR*ddata->scl_octv);21.for(j=0;j<SIFT_ORI_SMOOTH_PASSES；j++)22.smooth_ori_hist(hist,SIFT_ORI_HIST_BINS);23.omax=dominant_ori(hist,SIFT_ORI_