SIFT算法实现C语言.pdf

经典算法 SIFT 实现即代码解释：以下便是 sift 源码库编译后的效果图：为了给有兴趣实现sift 算法的朋友提供个参考，特整理此文如下。要了解什么是 sift 算法，请参考：九、图像特征提取与匹配之SIFT 算法。ok，咱们下面，就来利用 Rob Hess维护的 sift 库来实现 sift 算法：首先，请下载 Rob Hess维护的 sift 库：http:/blogs.oregonstate.edu/hess/code/sift/下载 Rob Hess的这个压缩包后，如果直接解压缩，直接编译，那么会出现下面的错误提示：编 译提 示:error C1083:Cannot open include file:cxcore.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 等。当然，你也可以用 vs 编译，同样下载 opencv（具体版本不受限制）、gsl 等。请按以下步骤操作：一、下载 opencv1.0 http:/ opencv1.0，配置 Windows 环境变量 1、安装注意：假如你是将OpenCV 安装到 C:/Program Files/OpenCV（如果你安装的时候选择不是安装在C 盘，则下面所有对应的C 盘都改为你所安装在的那个“X盘”，即可），在安装时选择 将/OpenCV/bin 加入系统变量，打上“勾”。（Add/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，在下方填入路径(参考下图)：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/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 modules 附加上：cxcore.lib cv.lib ml.lib cvaux.lib highgui.lib cvcam.lib 当然，你不需要这么多lib，你可以只添加你需要的lib(见下图)三、下载 gsl，gsl 也是一个库，也需要下载：http:/ GSL 也是和 OpenCV 一样要把头文件和lib 的路径指定好。四、配置 gsl 将 C:/WinGsl/bin 中的 WinGsl.dll 和 WinGslD.dll 复制到 C:/VC6.0/Bin；将整个Gsl 目录 复制 到 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，因为这个应用程序未能启动。重新安装应用程序可能会修复此问题。”的错误。手动在我的电脑-属性-高级-环境变量-系统变量-path 添加 c:/program files/opencv/bin;添加完成后需要重启计算机。II、vc6.0 下配置了一下，可是编译程序时遇到如下一个错误：Linking.LINK:fatal error LNK1104:cannot open fileodbccp32.libcxcore.lib 可能是：在工程设置的时候添加连接库时没加空格或.来把两个文件名（odbccp32.lib cxcore.lib）分开。注意每一次操作后，记得保存。若经过以上所有的步骤之后，如果还不能正常编译，那就是还要稍微修改下你下载的 Rob Hess代码。ok，日后，若有空，再好好详细剖析下此sift 的源码。最后，祝你编译顺利。完。SIFT 代码详解：这是一个很强大的算法，主要用于图像配准和物体识别等领域，但是其计算量相比也比较大，性价比比较高的算法包括PCA-SIFT 和 SURF 其中 OpenCV 提供了 SURF 算法，但是为了方便理解。这里给出了Rob Hess 所实现的 SIFT 算法的实现以及注释，结合我自己的理解，如果，您有关于SIFT 算法不理解的地方咱们可以一起交流一下。或者您认为不详细的地方提出来。SIFT 算法的主要实现在sift.c 这个文件，其主要流程为：(1)首先创建初始图像，即通过将图像转换为32 位的灰度图，然后将图像使用三次插值来方大，之后通过高斯模糊处理(2)在此基础上进行高斯金字塔的构建以及高斯差分金字塔的构建(3)对图像进行极值点检测(4)计算特征向量的尺度(5)调整图像大小(6)计算特征的方向(7)计算描述子，其中包括计算二维方向直方图并转换直方图为特征描述子首先给出 sift 算法的整体框架代码：输入参数：img 为输入图像；feat 为所要提取的特征指针；intvl 指的是高斯金字塔和差分金字塔的层数；sigma 指的是图像初始化过程中高斯模糊所使用的参数；contr_thr是归一化之后的去除不稳定特征的阈值；curv_thr指的是去除边缘的特征的主曲率阈值；img_dbl是是否将图像放大为之前的两倍；descr_with用来计算特征描述子的方向直方图的宽度；descr_hist_bins是直方图中的条数cppview plaincopy1.int _sift_features(IplImage*img,struct feature*feat,int intvls,2.double sigma,double contr_thr,int curv_thr,3.int img_dbl,int descr_width,int descr_hist_bins)4.5.IplImage*init_img;6.IplImage*gauss_pyr,*dog_pyr;7.CvMemStorage*storage;8.CvSeq*features;9.int octvs,i,n=0;10.11./*check arguments*/12.if(!img)13.fatal_error(NULL pointer error,%s,line%d,_FILE_,_LINE_);14.15.if(!feat)16.fatal_error(NULL pointer error,%s,line%d,_FILE_,_LINE_);17.18./*build scale space pyramid;smallest dimension of top level is 4 pixels*/19./*构建高斯尺度空间金字塔，顶层最小的为4像素 */20.init_img=create_init_img(img,img_dbl,sigma);21.octvs=log(double MIN(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(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./计算特征的方向43.calc_feature_oris(features,gauss_pyr);44./计算描述子，包括计算二维方向直方图和转换其为特征描述子pute_descriptors(features,gauss_pyr,descr_width,descr_hist_bins);46.47./*sort features by decreasing scale and move from CvSeq to array*/48.cvSeqSort(features,(CvCmpFunc)feature_cmp,NULL);49.n=features-total;50.*feat=static_cast(calloc(n,sizeof(struct feature);51.*feat=static_cast(cvCvtSeqToArray(features,*feat,CV_WHOLE_SEQ);52.53.54.55.56.for(i=0;i width*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.return dbl;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.return gray;22.23.(2)构建高斯金字塔输入参数：octvs 是高斯金字塔的组invls 是高斯金字塔的层数sigma 是初始的高斯模糊参数，后续也通过它计算每一层所使用的sigmacppview plaincopy1.static IplImage*build_gauss_pyr(IplImage*base,intoctvs,int intvls,double sigma)2.3.IplImage*gauss_pyr;4.double*sig=static_cast(calloc(intvls+3,sizeof(double);5.double sig_total,sig_prev,k;6.int i,o;7.8.gauss_pyr=static_cast(calloc(octvs,sizeof(IplImage*);9.for(i=0;i octvs;i+)10.gauss_pyri=static_cast(calloc(intvls+3,sizeof(IplImage*);11.12./*13.precompute Gaussian sigmas using the following formula: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;i intvls+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_prev*sig_prev);25.26.27.28.for(o=0;o octvs;o+)29.for(i=0;i intvls+3;i+)30.31./对每一层进行降采样，形成高斯金字塔的每一层32.if(o=0&i=0)33.gauss_pyroi=cvCloneImage(base);34.35./*base of new octvave is halved image from end of previous octave*/36./每一组的第一层都是通过对前面一组的第一层降采样实现的37.elseif(i=0)38.gauss_pyroi=downsample(gauss_pyro-1intvls);39.40./*blur the current octaves last image to create the next one*/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.return gauss_pyr;53.降采样处理输入参数：不解释这就是降采样，其实就是将图像通过最近邻算法缩小为原来的一半cppview plaincopy1.static IplImage*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.return smaller;8.(3)构建高斯差分金字塔输入参数：不解释了参见上面的说明即可实际上差分金字塔的构成是通过对相邻层的图像进行相减获得的cppview plaincopy1.static IplImage*build_dog_pyr(IplImage*gauss_pyr,int octvs,int intvls)2.3.IplImage*dog_pyr;4.int i,o;5.6.dog_pyr=static_cast(calloc(octvs,sizeof(IplImage*);7.for(i=0;i octvs;i+)8.dog_pyri=static_cast(calloc(intvls+2,sizeof(IplImage*);9.10.for(o=0;o octvs;o+)11.for(i=0;i intvls+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.return dog_pyr;19.(4)极值点检测输入参数：contr_thr 是去除对比度低的点所采用的阈值curv_thr 是去除边缘特征的阈值cppview plaincopy1.static CvSeq*scale_space_extrema(IplImage*dog_pyr,int octvs,int intvls,2.double contr_thr,int curv_thr,3.CvMemStorage*storage)4.5.CvSeq*features;6.double prelim_contr_thr=0.5*contr_thr/intvls;7.struct feature*feat;8.struct detection_data*ddata;9.int o,i,r,c;10.11.features=cvCreateSeq(0,sizeof(CvSeq),sizeof(struct feature),storage);12.for(o=0;o octvs;o+)13.for(i=1;i=intvls;i+)14.for(r=SIFT_IMG_BORDER;r height-SIFT_IMG_BORDER;r+)15.for(c=SIFT_IMG_BORDER;c width-SIFT_IMG_BORDER;c+)16./*perform preliminary check on contrast*/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_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.return features;36.SIFT_IMG_BORDER是预定义的图像边缘；通过和对比度阈值比较去掉低对比度的点；而通过 is_extremum 来判断是否为极值点，如果是则通过极值点插值的方式获取亚像素的极值点的位置。然后通过 is_too_eage_like和所给的主曲率阈值判断是否为边缘点*判断是否为极值点其原理为:通过和高斯金字塔的上一层的9 个像素+本层的除了本像素自己的其他的 8 个像素和下一层的9 个像素进行比较看是否为这26 个像素中最小的一个或者是否为最大的一个，如果是则为极值点。cppview plaincopy1.staticint is_extremum(IplImage*dog_pyr,int octv,int intvl,int r,int c)2.3.float val=pixval32f(dog_pyroctvintvl,r,c);4.int i,j,k;5.6./*check for maximum*/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_pyroctvintvl+i,r+j,c+k)13.return 0;14.15.16./*check for minimum*/17.else18.19.for(i=-1;i=1;i+)20.for(j=-1;j=1;j+)21.for(k=-1;k pixval32f(dog_pyroctvintvl+i,r+j,c+k)23.return 0;24.25.26.return 1;27.*获取亚像素的极值点的位置cppview plaincopy1.staticstruct feature*interp_extremum(IplImage*dog_pyr,int octv,int intvl,2.int r,int c,int intvls,double contr_thr)3.4.struct feature*feat;5.struct detection_data*ddata;6.double xi,xr,xc,contr;/分别为亚像素的 intval,row,col的偏移 offset，和对比度7.int i=0;8.9.while(i SIFT_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+=cvRound(xi);18.19.if(intvl intvls|21.c SIFT_IMG_BORDER|22.r=dog_pyroctv0-width-SIFT_IMG_BORDER|24.r=dog_pyroctv0-height-SIFT_IMG_BORDER)25.26.return NULL;27.28.29.i+;30.31.32./确保极值点是经过最大5 步找到的33./*ensure convergence of interpolation*/34.if(i=SIFT_MAX_INTERP_STEPS)35.return NULL;36.37./获取找到的极值点的对比度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.return feat;54.*获取亚像素位置中所用到的函数cppview plaincopy1.staticvoid interp_step(IplImage*dog_pyr,int octv,int intvl,int r,int c,2.double*xi,double*xr,double*xc)3.4.CvMat*dD,*H,*H_inv,X;5.double x3=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=cvCreateMat(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 方向上的偏导数，高斯差分尺度空间金字塔中像素的尺度实际上在离散数据中计算偏导数是通过相邻像素的相减来计算的比如说计算x 方向的偏导数dx，则通过该向所的x 方向的后一个减去前一个然后除以 2 即可求的 dxcppview plaincopy1.static CvMat*deriv_3D(IplImage*dog_pyr,int octv,int intvl,int r,int c)2.3.CvMat*dI;4.double dx,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+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.return dI;19.*计算三维海森矩阵不需要讲什么，其实就是计算二次导数，计算方法也和一次导数的计算如出一辙。然后将结果放入到一个矩阵中去。cppview plaincopy1.static CvMat*hessian_3D(IplImage*dog_pyr,int octv,int intvl,int r,int c)2.3.CvMat*H;4.double v,dxx,dyy,dss,dxy,dxs,dys;5.6.v=pixval32f(dog_pyroctvintvl,r,c);7.dxx=(pixval32f(dog_pyroctvintvl,r,c+1)+8.pixval32f(dog_pyroctvintvl,r,c-1)-2*v);9.dyy=(pixval32f(dog_pyroctvintvl,r+1,c)+10.pixval32f(dog_pyroctvintvl,r-1,c)-2*v);11.dss=(pixval32f(dog_pyroctvintvl+1,r,c)+12.pixval32f(dog_pyroctvintvl-1,r,c)-2*v);13.dxy=(pixval32f(dog_pyroctvintvl,r+1,c+1)-14.pixval32f(dog_pyroctvintvl,r+1,c-1)-15.pixval32f(dog_pyroctvintvl,r-1,c+1)+16.pixval32f(dog_pyroctvintvl,r-1,c-1)/4.0;17.dxs=(pixval32f(dog_pyroctvintvl+1,r,c+1)-18.pixval32f(dog_pyroctvintvl+1,r,c-1)-19.pixval32f(dog_pyroctvintvl-1,r,c+1)+20.pixval32f(dog_pyroctvintvl-1,r,c-1)/4.0;21.dys=(pixval32f(dog_pyroctvintvl+1,r+1,c)-22.pixval32f(dog_pyroctvintvl+1,r-1,c)-23.pixval32f(dog_pyroctvintvl-1,r+1,c)+24.pixval32f(dog_pyroctvintvl-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.37.return H;38.*计算插入像素的对比度cppview plaincopy1.staticdouble interp_contr(IplImage*dog_pyr,int octv,int intvl,int r,2.int c,double xi,double xr,double xc)3.4.CvMat*dD,X,T;5.double t1,x3=xc,xr,xi;6.7.cvInitMatHeader(&X,3,1,CV_64FC1,x,CV_AUTOSTEP);8.cvInitMatHeader(&T,1,1,CV_64FC1,t,CV_AUTOSTEP);9.dD=deriv_3D(dog_pyr,octv,intvl,r,c);10.cvGEMM(dD,&X,1,NULL,0,&T,CV_GEMM_A_T);11.cvReleaseMat(&dD);12.13.return pixval32f(dog_pyroctvintvl,r,c)+t0*0.5;14.其中 cvGEMM 是矩阵的通用计算函数，至于CV_GEMM_A_T是计算 dD 的转置矩阵放入T 中*去除边缘相应通过计算所在特征向量的主曲率半径来判断特征是边缘的从而导致不稳定即去除边缘响应cppview plaincopy1.staticint is_too_edge_like(IplImage*dog_img,int r,int c,int curv_thr)2.3.double d,dxx,dyy,dxy,tr,det;4.5./*principal curvatures are computed using the trace and det of Hessian*/6.d=pixval32f(dog_img,r,c);7.dxx=pixval32f(dog_img,r,c+1)+pixval32f(dog_img,r,c-1)-2*d;8.dyy=pixval32f(dog_img,r+1,c)+pixval32f(dog_img,r-1,c)-2*d;9.dxy=(pixval32f(dog_img,r+1,c+1)-pixval32f(dog_img,r+1,c-1)-10.pixval32f(dog_img,r-1,c+1)+pixval32f(dog_img,r-1,c-1)/4.0;11.tr=dxx+dyy;12.det=dxx*dyy-dxy*dxy;13.14./*negative determinant-curvatures have different signs;reject feature*/15.if(det=0)16.return 1;17.18.if(tr*tr/det total;9.for(i=0;i intvl+ddata-subintvl;14.feat-scl=sigma*pow(2.0,ddata-octv+intvl/intvls);15.ddata-scl_octv=sigma*pow(2.0,intvl/intvls);16.17.(5)调整图像特征坐标、尺度、点的坐标大小为原来的一半javaview plaincopy1.staticvoid adjust_for_img_dbl(CvSeq*features)2.3.struct feature*feat;4.int i,n;5.6.n=features-total;7.for(i=0;i x/=2.0;11.feat-y/=2.0;12.feat-scl/=2.0;13.feat-img_pt.x/=2.0;14.feat-img_pt.y/=2.0;15.16.(6)给每一个图像特征向量计算规范化的方向cppview plaincopy1.staticvoid calc_feature_oris(CvSeq*features,IplImage*gauss_pyr)2.3.struct feature*feat;4.struct detection_data*ddata;5.double*hist;6.double omax;7.int i,j,n=features-total;8.9.10./遍历整个检测出来的特征点，计算每个特征点的直方图，然后平滑直方图去除突变，然后找到每一个特征点的主方向，并加入到好的方向特征数组中去11.for(i=0;i n;i+)12.13.feat=static_cast(malloc(sizeof(struct feature);14.cvSeqPopFront(features,feat);15.ddata=feat_detection_data(feat);16./计算给定的某个像素的灰度方向直方图17.hist=ori_hist(gauss_pyrddata-octvddata-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_HIST_BINS);24.25./描述子向量元素门限化26.add_good_ori_features(features,hist,SIFT_ORI_HIST_BINS,27.omax*SIFT_ORI_PEAK_RATIO,feat);28.free(ddata);29.free(feat);30.free(hist);31.32.*对所给像素计算灰度方向直方图以关键点为中心的邻域窗口内采样，并用直方图统计邻域像素的梯度方向。梯度直方图的范围是0360 度，其中每10 度一个柱，总共36 个柱cppview plaincopy1.staticdouble*ori_hist(IplImage*img,int r,int c,int n,int rad,do