SIFT算法实现C语言

1、经典算法SIFT实现即代码解释：以下便是sift源码库编译后的效果图： 为了给有兴趣实现sift算法的朋友提供个参考，特整理此文如下。要了解什么是sift算法，请参考：九、图像特征提取与匹配之SIFT算法。ok，咱们下面，就来利用Rob Hess维护的sift 库来实现sift算法： 首先，请下载Rob Hess维护的sift 库：/hess/code/sift/ 下载Rob Hess的这个压缩包后，如果直接解压缩，直接编译，那么会出现下面的错误提示：编译提示:error C1083: Cannot open include file: c

2、xcore.h: No such file or directory，找不到这个头文件。 这个错误，是因为你还没有安装opencv，因为：cxcore.h和cv.h是开源的OPEN CV头文件,不是VC+的默认安装文件,所以你还得下载OpenCV并进行安装。然后，可以在OpenCV文件夹下找到你所需要的头文件了。 据网友称，截止2010年4月4日，还没有在VC6.0下成功使用opencv2.0的案例。所以，如果你是VC6.0的用户请下载opencv1.0版本。vs的话，opencv2.0,1.0任意下载。 以下，咱们就以vc6.0为平台举例，下载并安装opencv1.0版本、gsl等。当然，你

3、也可以用vs编译，同样下载opencv（具体版本不受限制）、gsl等。 请按以下步骤操作：一、下载opencv1.0/projects/opencvlibrary/files/opencv-win/1.0/OpenCV_1.0.exe/download二、安装opencv1.0，配置Windows环境变量1、安装注意：假如你是将OpenCV安装到C:/Program Files/OpenCV（如果你安装的时候选择不是安装在C盘，则下面所有对应的C盘都改为你所安装在的那个“X盘”，即可），在安装时选择将/OpenCV/bin加入系统变量，打上“勾”。（A

4、dd/OpenCV/bin to the systerm PATH。这一步确认选上了之后，下面的检查环境变量的步骤，便可免去）2、检查环境变量。为了确保上述步骤中，加入了系统变量，在安装opencv1.0成功后，还得检查C:/Program Files/OpenCV/bin是否已经被加入到环境变量PATH，如果没有，请加入。3、最后是配置Visual C+ 6.0。全局设置 菜单Tools-Options-Directories：先设置lib路径，选择Library files，在下方填入路径：C:/Program Files/OpenCV/lib 然后选择include files，在下方填

5、入路径(参考下图)： C:/Program Files/OpenCV/cxcore/include C:/Program Files/OpenCV/cv/include C:/Program Files/OpenCV/cvaux/include C:/Program Files/OpenCV/ml/include C:/Program Files/OpenCV/otherlibs/highgui C:/Program Files/OpenCV/otherlibs/cvcam/include 最后选择source files，在下方填入路径： C:/Program Files/OpenCV/cv

6、/src C:/Program Files/OpenCV/cxcore/src C:/Program Files/OpenCV/cvaux/src C:/Program Files/OpenCV/otherlibs/highgui C:/Program Files/OpenCV/otherlibs/cvcam/src/windows项目设置 每创建一个将要使用OpenCV的VC Project，都需要给它指定需要的lib。菜单：Project-Settings，然后将Setting for选为All Configurations，然后选择右边的link标签，在Object/library mo

7、dules附加上： cxcore.lib cv.lib ml.lib cvaux.lib highgui.lib cvcam.lib 当然，你不需要这么多lib，你可以只添加你需要的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复制到C:/VC6.0/Bin；将整个Gsl目

8、录复制到C:/VC6.0/Bin下；lib目录下的所有.lib文件全部复制到C:/VC6.0/Lib下。 然后，在toolsoptionsdirectories中，将C:/WinGsl下的lib，gsl分别加入到库文件和头文件的搜索路径中。以下是可能会出现的错误情况处理：I、OpenCV安装后“没有找到cxcore100.dll”的错误处理在安装时选择“将/OpenCV/bin加入系统变量”（Add/OpenCV/bin to the systerm PATH）。 但该选项并不一定能成功添加到系统变量，如果编写的程序在运行时出现“没有找到cxcore100.dll，因为这个应用程序未能启动。重

9、新安装应用程序可能会修复此问题。”的错误。 手动在我的电脑-属性-高级-环境变量-系统变量-path添加c:/program files/opencv/bin;添加完成后需要重启计算机。II、vc6.0下配置了一下，可是编译程序时遇到如下一个错误：Linking. LINK : fatal error LNK1104: cannot open fileodbccp32.libcxcore.lib 可能是：在工程设置的时候添加连接库时没加空格或.来把两个文件名（odbccp32.lib cxcore.lib）分开。注意每一次操作后，记得保存。 若经过以上所有的步骤之后，如果还不能正常编译，那就是

10、还要稍微修改下你下载的Rob Hess代码。ok，日后，若有空，再好好详细剖析下此sift的源码。最后，祝你编译顺利。 完。SIFT代码详解：这是一个很强大的算法，主要用于图像配准和物体识别等领域，但是其计算量相比也比较大，性价比比较高的算法包括PCA-SIFT和SURF其中OpenCV提供了SURF算法，但是为了方便理解。这里给出了Rob Hess所实现的SIFT算法的实现以及注释，结合我自己的理解，如果，您有关于SIFT算法不理解的地方咱们可以一起交流一下。或者您认为不详细的地方提出来。 SIFT算法的主要实现在sift.c这个文件，其主要流程为：(1)首先创建初始图像，即通过将图像转换为

11、32位的灰度图，然后将图像使用三次插值来方大，之后通过高斯模糊处理(2)在此基础上进行高斯金字塔的构建以及高斯差分金字塔的构建(3)对图像进行极值点检测(4)计算特征向量的尺度(5)调整图像大小(6)计算特征的方向(7)计算描述子，其中包括计算二维方向直方图并转换直方图为特征描述子首先给出sift算法的整体框架代码：输入参数：img为输入图像；feat为所要提取的特征指针；intvl指的是高斯金字塔和差分金字塔的层数；sigma指的是图像初始化过程中高斯模糊所使用的参数； contr_thr是归一化之后的去除不稳定特征的阈值；curv_thr指的是去除边缘的特征的主曲率阈值；img_dbl是是

12、否将图像放大为之前的两倍；descr_with用来计算特征描述子的方向直方图的宽度；descr_hist_bins是直方图中的条数cppview plaincopy1. int_sift_features(IplImage*img,structfeature*feat,intintvls,2. doublesigma,doublecontr_thr,intcurv_thr,3. intimg_dbl,intdescr_width,intdescr_hist_bins)4. 5. IplImage*init_img;6. IplImage*gauss_pyr,*dog_pyr;7. CvMemS

13、torage*storage;8. CvSeq*features;9. intoctvs,i,n=0;10. 11. /*checkarguments*/12. if(!img)13. fatal_error(NULLpointererror,%s,line%d,_FILE_,_LINE_);14. 15. if(!feat)16. fatal_error(NULLpointererror,%s,line%d,_FILE_,_LINE_);17. 18. /*buildscalespacepyramid;smallestdimensionoftoplevelis4pixels*/19. /*构

14、建高斯尺度空间金字塔，顶层最小的为4像素*/20. init_img=create_init_img(img,img_dbl,sigma);21. octvs=log(doubleMIN(init_img-width,init_img-height)/log(2.0)-2;22. /构建高斯金字塔和高斯差分金字塔23. gauss_pyr=build_gauss_pyr(init_img,octvs,intvls,sigma);24. dog_pyr=build_dog_pyr(gauss_pyr,octvs,intvls);25. 26. storage=cvCreateMemStorage

15、(0);27. 28. /尺度空间极值点检测29. features=scale_space_extrema(dog_pyr,octvs,intvls,contr_thr,30. curv_thr,storage);31. 32. /画出去除低对比度的极值点33. /draw_extrempoint(img,features);34. 35. 36. 37. 38. /计算特征向量的尺度39. calc_feature_scales(features,sigma,intvls);40. if(img_dbl)41. adjust_for_img_dbl(features);42. /计算特征的

16、方向43. calc_feature_oris(features,gauss_pyr);44. /计算描述子，包括计算二维方向直方图和转换其为特征描述子45. compute_descriptors(features,gauss_pyr,descr_width,descr_hist_bins);46. 47. /*sortfeaturesbydecreasingscaleandmovefromCvSeqtoarray*/48. cvSeqSort(features,(CvCmpFunc)feature_cmp,NULL);49. n=features-total;50. *feat=stati

17、c_cast(calloc(n,sizeof(structfeature);51. *feat=static_cast(cvCvtSeqToArray(features,*feat,CV_WHOLE_SEQ);52. 53. 54. 55. 56. for(i=0;iwidth*2,img-height*2),11. IPL_DEPTH_32F,1);12. cvResize(gray,dbl,CV_INTER_CUBIC);13. cvSmooth(dbl,dbl,CV_GAUSSIAN,0,0,sig_diff,sig_diff);14. cvReleaseImage(&gray);15.

18、 returndbl;16. 17. else18. 19. sig_diff=sqrt(sigma*sigma-SIFT_INIT_SIGMA*SIFT_INIT_SIGMA);20. cvSmooth(gray,gray,CV_GAUSSIAN,0,0,sig_diff,sig_diff);21. returngray;22. 23. (2)构建高斯金字塔输入参数：octvs是高斯金字塔的组invls是高斯金字塔的层数sigma是初始的高斯模糊参数，后续也通过它计算每一层所使用的sigmacppview plaincopy1. staticIplImage*build_gauss_pyr(

19、IplImage*base,intoctvs,intintvls,doublesigma)2. 3. IplImage*gauss_pyr;4. double*sig=static_cast(calloc(intvls+3,sizeof(double);5. doublesig_total,sig_prev,k;6. inti,o;7. 8. gauss_pyr=static_cast(calloc(octvs,sizeof(IplImage*);9. for(i=0;ioctvs;i+)10. gauss_pyri=static_cast(calloc(intvls+3,sizeof(Ipl

20、Image*);11. 12. /*13. precomputeGaussiansigmasusingthefollowingformula:14. 预计算每次高斯模糊的sigma15. 16. sigma_total2=sigma_i2+sigma_i-1217. */18. sig0=sigma;19. k=pow(2.0,1.0/intvls);20. for(i=1;iintvls+3;i+)21. 22. sig_prev=pow(k,i-1)*sigma;23. sig_total=sig_prev*k;24. sigi=sqrt(sig_total*sig_total-sig_p

21、rev*sig_prev);25. 26. 27. 28. for(o=0;ooctvs;o+)29. for(i=0;iintvls+3;i+)30. 31. /对每一层进行降采样，形成高斯金字塔的每一层32. if(o=0&i=0)33. gauss_pyroi=cvCloneImage(base);34. 35. /*baseofnewoctvaveishalvedimagefromendofpreviousoctave*/36. /每一组的第一层都是通过对前面一组的第一层降采样实现的37. elseif(i=0)38. gauss_pyroi=downsample(gauss_pyro

22、-1intvls);39. 40. /*blurthecurrentoctaveslastimagetocreatethenextone*/41. /每一组的其他层则使通过使用不同sigma的高斯模糊来进行处理42. else43. 44. gauss_pyroi=cvCreateImage(cvGetSize(gauss_pyroi-1),45. IPL_DEPTH_32F,1);46. cvSmooth(gauss_pyroi-1,gauss_pyroi,47. CV_GAUSSIAN,0,0,sigi,sigi);48. 49. 50. 51. free(sig);52. returng

23、auss_pyr;53. 降采样处理输入参数：不解释这就是降采样，其实就是将图像通过最近邻算法缩小为原来的一半cppview plaincopy1. staticIplImage*downsample(IplImage*img)2. 3. IplImage*smaller=cvCreateImage(cvSize(img-width/2,img-height/2),4. img-depth,img-nChannels);5. cvResize(img,smaller,CV_INTER_NN);6. 7. returnsmaller;8. (3)构建高斯差分金字塔输入参数：不解释了参见上面的说明

24、即可实际上差分金字塔的构成是通过对相邻层的图像进行相减获得的cppview plaincopy1. staticIplImage*build_dog_pyr(IplImage*gauss_pyr,intoctvs,intintvls)2. 3. IplImage*dog_pyr;4. inti,o;5. 6. dog_pyr=static_cast(calloc(octvs,sizeof(IplImage*);7. for(i=0;ioctvs;i+)8. dog_pyri=static_cast(calloc(intvls+2,sizeof(IplImage*);9. 10. for(o=0

25、;ooctvs;o+)11. for(i=0;iintvls+2;i+)12. 13. dog_pyroi=cvCreateImage(cvGetSize(gauss_pyroi),14. IPL_DEPTH_32F,1);15. cvSub(gauss_pyroi+1,gauss_pyroi,dog_pyroi,NULL);16. 17. 18. returndog_pyr;19. (4)极值点检测输入参数：contr_thr是去除对比度低的点所采用的阈值curv_thr是去除边缘特征的阈值cppview plaincopy1. staticCvSeq*scale_space_extrema

26、(IplImage*dog_pyr,intoctvs,intintvls,2. doublecontr_thr,intcurv_thr,3. CvMemStorage*storage)4. 5. CvSeq*features;6. doubleprelim_contr_thr=0.5*contr_thr/intvls;7. structfeature*feat;8. structdetection_data*ddata;9. into,i,r,c;10. 11. features=cvCreateSeq(0,sizeof(CvSeq),sizeof(structfeature),storage

27、);12. for(o=0;ooctvs;o+)13. for(i=1;i=intvls;i+)14. for(r=SIFT_IMG_BORDER;rheight-SIFT_IMG_BORDER;r+)15. for(c=SIFT_IMG_BORDER;cwidth-SIFT_IMG_BORDER;c+)16. /*performpreliminarycheckoncontrast*/17. if(ABS(pixval32f(dog_pyroi,r,c)prelim_contr_thr)18. if(is_extremum(dog_pyr,o,i,r,c)19. 20. feat=interp

28、_extremum(dog_pyr,o,i,r,c,intvls,contr_thr);21. if(feat)22. 23. ddata=feat_detection_data(feat);24. if(!is_too_edge_like(dog_pyrddata-octvddata-intvl,25. ddata-r,ddata-c,curv_thr)26. 27. cvSeqPush(features,feat);28. 29. else30. free(ddata);31. free(feat);32. 33. 34. 35. returnfeatures;36. SIFT_IMG_B

29、ORDER是预定义的图像边缘；通过和对比度阈值比较去掉低对比度的点；而通过is_extremum来判断是否为极值点，如果是则通过极值点插值的方式获取亚像素的极值点的位置。然后通过is_too_eage_like和所给的主曲率阈值判断是否为边缘点*判断是否为极值点其原理为:通过和高斯金字塔的上一层的9个像素+本层的除了本像素自己的其他的8个像素和下一层的9个像素进行比较看是否为这26个像素中最小的一个或者是否为最大的一个，如果是则为极值点。cppview plaincopy1. staticintis_extremum(IplImage*dog_pyr,intoctv,intintvl,intr

30、,intc)2. 3. floatval=pixval32f(dog_pyroctvintvl,r,c);4. inti,j,k;5. 6. /*checkformaximum*/7. if(val0)8. 9. for(i=-1;i=1;i+)10. for(j=-1;j=1;j+)11. for(k=-1;k=1;k+)12. if(valpixval32f(dog_pyroctvintvl+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=

31、1;j+)21. for(k=-1;kpixval32f(dog_pyroctvintvl+i,r+j,c+k)23. return0;24. 25. 26. return1;27. *获取亚像素的极值点的位置cppview plaincopy1. staticstructfeature*interp_extremum(IplImage*dog_pyr,intoctv,intintvl,2. intr,intc,intintvls,doublecontr_thr)3. 4. structfeature*feat;5. structdetection_data*ddata;6. doublexi

32、,xr,xc,contr;/分别为亚像素的intval,row,col的偏移offset，和对比度7. inti=0;8. 9. while(iSIFT_MAX_INTERP_STEPS)/重新确定极值点并重新定位的操作只能循环5次10. 11. interp_step(dog_pyr,octv,intvl,r,c,&xi,&xr,&xc);12. if(ABS(xi)0.5&ABS(xr)0.5&ABS(xc)0.5)/如果满足条件就停止寻找13. break;14. /否则继续寻找极值点15. c+=cvRound(xc);16. r+=cvRound(xr);17. intvl+=cvR

33、ound(xi);18. 19. if(intvlintvls|21. cSIFT_IMG_BORDER|22. r=dog_pyroctv0-width-SIFT_IMG_BORDER|24. r=dog_pyroctv0-height-SIFT_IMG_BORDER)25. 26. returnNULL;27. 28. 29. i+;30. 31. 32. /确保极值点是经过最大5步找到的33. /*ensureconvergenceofinterpolation*/34. if(i=SIFT_MAX_INTERP_STEPS)35. returnNULL;36. 37. /获取找到的极值

34、点的对比度38. contr=interp_contr(dog_pyr,octv,intvl,r,c,xi,xr,xc);39. /判断极值点是否小于某一个阈值40. if(ABS(contr)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. *获取亚像

35、素位置中所用到的函数cppview plaincopy1. 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. doublex3=0;6. 7. 8. /计算三维偏导数9. dD=deriv_3D(dog_pyr,octv,intvl,r,c);10. /计算三维海森矩阵11. H=hessian_3D(dog_pyr,octv,intvl,r,c);12. H_inv=cvCreateMa

36、t(3,3,CV_64FC1);13. cvInvert(H,H_inv,CV_SVD);14. cvInitMatHeader(&X,3,1,CV_64FC1,x,CV_AUTOSTEP);15. 16. cvGEMM(H_inv,dD,-1,NULL,0,&X,0);17. 18. cvReleaseMat(&dD);19. cvReleaseMat(&H);20. cvReleaseMat(&H_inv);21. 22. *xi=x2;23. *xr=x1;24. *xc=x0;25. *计算三维偏导数计算在x和y方向上的偏导数，高斯差分尺度空间金字塔中像素的尺度实际上在离散数据中计算偏

37、导数是通过相邻像素的相减来计算的比如说计算x方向的偏导数dx，则通过该向所的x方向的后一个减去前一个然后除以2即可求的dxcppview plaincopy1. staticCvMat*deriv_3D(IplImage*dog_pyr,intoctv,intintvl,intr,intc)2. 3. CvMat*dI;4. doubledx,dy,ds;5. 6. dx=(pixval32f(dog_pyroctvintvl,r,c+1)-7. pixval32f(dog_pyroctvintvl,r,c-1)/2.0;8. dy=(pixval32f(dog_pyroctvintvl,r+

38、1,c)-9. pixval32f(dog_pyroctvintvl,r-1,c)/2.0;10. ds=(pixval32f(dog_pyroctvintvl+1,r,c)-11. pixval32f(dog_pyroctvintvl-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. *计算三维海森矩阵不需要讲什么，其实就是计算二次导数，计算方法也和一次导数的计算如出一辙。然后将结果放入到一个矩阵中去。cppview plaincopy1. staticCvMat*hessian_3D(IplImage*dog_pyr,intoctv,inti