SIFT算法实现及代码详解
sift算法C代码详解
#ifndef SIFT_H#define SIFT_H#include "cxcore.h"/******************************** Structures ************************ *********///极值点检测中用到的结构//在SIFT特征提取过程中，此类型数据会被赋值给feature结构的feature_data成员struct detection_data{int r; //特征点所在的行int c; //特征点所在的列int octv; //高斯差分金字塔中，特征点所在的组int intvl; //高斯差分金字塔中，特征点所在的组中的层double subintvl; //特征点在层方向(σ方向,intvl方向)上的亚像素偏移量double scl_octv; //特征点所在的组的尺度};struct feature;/******************************* 一些默认参数 ************************* ****///高斯金字塔每组内的层数#define SIFT_INTVLS 3//第0层的初始尺度，即第0层高斯模糊所使用的参数#define SIFT_SIGMA 1.6//对比度阈值，针对归一化后的图像，用来去除不稳定特征/** default threshold on keypoint contrast |D(x)| */#define SIFT_CONTR_THR 0.04//主曲率比值的阈值，用来去除边缘特征/** default threshold on keypoint ratio of principle curvatures */#define SIFT_CURV_THR 10//是否将图像放大为之前的两倍/** double image size before pyramid construction? */#define SIFT_IMG_DBL 1//计算特征描述子过程中，计算方向直方图时，将特征点附近划分为d*d个区域，//每个区域生成一个直方图，SIFT_DESCR_WIDTH即d的默认值/** default width of descriptor histogram array */#define SIFT_DESCR_WIDTH 4//计算特征描述子过程中，每个方向直方图的bin个数/** default number of bins per histogram in descriptor array */#define SIFT_DESCR_HIST_BINS 8//输入图像的尺度为0.5/* assumed gaussian blur for input image */#define SIFT_INIT_SIGMA 0.5//边界的像素宽度，检测过程中将忽略边界线中的极值点，即只检测边界线以内是否存在极值点/* width of border in which to ignore keypoints */#define SIFT_IMG_BORDER 5//通过插值进行极值点精确定位时，最大差值次数，即关键点修正次数/* maximum steps of keypoint interpolation before failure */#define SIFT_MAX_INTERP_STEPS 5//特征点方向赋值过程中，梯度方向直方图中柱子(bin)的个数/* default number of bins in histogram for orientation assignment */#define SIFT_ORI_HIST_BINS 36//特征点方向赋值过程中，搜索邻域的半径为：3 * 1.5 * σ/* determines gaussian sigma for orientation assignment */#define SIFT_ORI_SIG_FCTR 1.5//特征点方向赋值过程中，搜索邻域的半径为：3 * 1.5 * σ/* determines the radius of the region used in orientation assignment * /#define SIFT_ORI_RADIUS 3.0 * SIFT_ORI_SIG_FCTR//特征点方向赋值过程中，梯度方向直方图的平滑次数，计算出梯度直方图后还要进行高斯平滑/* number of passes of orientation histogram smoothing */#define SIFT_ORI_SMOOTH_PASSES 2//特征点方向赋值过程中，梯度幅值达到最大值的80%则分裂为两个特征点/* orientation magnitude relative to max that results in new feature */#define SIFT_ORI_PEAK_RATIO 0.8//计算特征描述子过程中，特征点周围的d*d个区域中，每个区域的宽度为m*σ个像素，//SIFT_DESCR_SCL_FCTR即m的默认值，σ为特征点的尺度/* determines the size of a single descriptor orientation histogram */#define SIFT_DESCR_SCL_FCTR 3.0//计算特征描述子过程中，特征描述子向量中元素的阈值(最大值，并且是针对归一化后//的特征描述子)，超过此阈值的元素被强行赋值为此阈值/* threshold on magnitude of elements of descriptor vector */#define SIFT_DESCR_MAG_THR 0.2//计算特征描述子过程中，将浮点型的特征描述子变为整型时乘以的系数/* factor used to convert floating-point descriptor to unsigned char * /#define SIFT_INT_DESCR_FCTR 512.0//定义了一个带参数的函数宏，用来提取参数f中的feature_data成员并转换为detecti on_data格式的指针/* returns a feature's detection data */#define feat_detection_data(f) ( (struct detection_data*)(f->feature_ data) )/*************************** Function Prototypes ******************** *********//*使用默认参数在图像中提取SIFT特征点参数：img：图像指针feat：用来存储特征点的feature数组的指针static IplImage* downsample( IplImage* );//通过对高斯金字塔中每相邻两层图像相减来建立高斯差分金字塔static IplImage*** build_dog_pyr( IplImage***, int, int );//在尺度空间中检测极值点，通过插值精确定位，去除低对比度的点，去除边缘点，返回检测到的特征点序列static CvSeq* scale_space_extrema( IplImage***, int, int, double, int, CvMemStorage*);//通过在尺度空间中将一个像素点的值与其周围3*3*3邻域内的点比较来决定此点是否极值点(极大值或极小都行)static int is_extremum( IplImage***, int, int, int, int );//通过亚像素级插值进行极值点精确定位(修正极值点坐标)，并去除低对比度的极值点，将修正后的特征点组成feature结构返回static struct feature* interp_extremum( IplImage***, int, int, int, in t, int, double);//进行一次极值点差值，计算x，y，σ方向(层方向)上的子像素偏移量(增量)static void interp_step( IplImage***, int, int, int, int, double*, dou ble*, double* );//在DoG金字塔中计算某点的x方向、y方向以及尺度方向上的偏导数static CvMat* deriv_3D( IplImage***, int, int, int, int );//在DoG金字塔中计算某点的3*3海森矩阵static CvMat* hessian_3D( IplImage***, int, int, int, int );//计算被插值点的对比度：D + 0.5 * dD^T * Xstatic double interp_contr( IplImage***, int, int, int, int, double, d ouble, double );//为一个feature结构分配空间并初始化static struct feature* new_feature( void );//去除边缘响应，即通过计算主曲率比值判断某点是否边缘点static int is_too_edge_like( IplImage*, int, int, int );//计算特征点序列中每个特征点的尺度static void calc_feature_scales( CvSeq*, double, int );//将特征点序列中每个特征点的坐标减半(当设置了将图像放大为原图的2倍时，特征点检测完之后调用)static void adjust_for_img_dbl( CvSeq* );//计算每个特征点的梯度直方图，找出其主方向，若一个特征点有不止一个主方向，将其分为两个特征点static void calc_feature_oris( CvSeq*, IplImage*** );//计算指定像素点的梯度方向直方图，返回存放直方图的数组static double* ori_hist( IplImage*, int, int, int, int, double );//计算指定点的梯度的幅值magnitude和方向orientationstatic int calc_grad_mag_ori( IplImage*, int, int, double*, double* );//对梯度方向直方图进行高斯平滑，弥补因没有仿射不变性而产生的特征点不稳定的问题static void smooth_ori_hist( double*, int );//查找梯度直方图中主方向的梯度幅值，即查找直方图中最大bin的值static double dominant_ori( double*, int );//若当前特征点的直方图中某个bin的值大于给定的阈值，则新生成一个特征点并添加到特征点序列末尾static void add_good_ori_features( CvSeq*, double*, int, double, struc t feature* );//对输入的feature结构特征点做深拷贝，返回克隆生成的特征点的指针static struct feature* clone_feature( struct feature* );//计算特征点序列中每个特征点的特征描述子向量static void compute_descriptors( CvSeq*, IplImage***, int, int );//计算特征点附近区域的方向直方图，此直方图在计算特征描述子中要用到，返回值是一个d* d*n的三维数组static double*** descr_hist( IplImage*, int, int, double, double, int, int );static void interp_hist_entry( double***, double, double, double, doub le, int, int);//将某特征点的方向直方图转换为特征描述子向量，对特征描述子归一化并将所有元素转化为整型，存入指定特征点中static void hist_to_descr( double***, int, int, struct feature* );//归一化特征点的特征描述子，即将特征描述子数组中每个元素除以特征描述子的模static void normalize_descr( struct feature* );//比较函数，将特征点按尺度的降序排列，用在序列排序函数CvSeqSort中static int feature_cmp( void*, void*, void* );//释放计算特征描述子过程中用到的方向直方图的内存空间static void release_descr_hist( double****, int );//释放金字塔图像组的存储空间static void release_pyr( IplImage****, int, int );/*********************** Functions prototyped in sift.h ************** ********//**Finds SIFT features in an image using default parameter values. Alldetected features are stored in the array pointed to by \a feat.@param img the image in which to detect features@param feat a pointer to an array in which to store detected features@return Returns the number of features stored in \a feat or -1 on failu re@see _sift_features()*/int sift_features( IplImage* img, struct feature** feat ){return _sift_features( img, feat, SIFT_INTVLS, SIFT_SIGMA, SIF T_CONTR_THR,SIFT_CURV_THR, SIFT_IMG_DBL, SIFT_DESCR_WIDTH,SIFT_DESCR_HIST _BINS );}/*使用用户指定的参数在图像中提取SIFT特征点参数：img：输入图像feat：存储特征点的数组的指针，此数组的内存将在本函数中被分配，使用完后必须在调用出释放：free(*feat)intvls：每组的层数sigma：初始高斯平滑参数σcontr_thr：对比度阈值，针对归一化后的图像，用来去除不稳定特征curv_thr：去除边缘的特征的主曲率阈值img_dbl：是否将图像放大为之前的两倍descr_width：特征描述过程中，计算方向直方图时，将特征点附近划分为descr_width*d escr_width个区域，每个区域生成一个直方图descr_hist_bins：特征描述过程中，每个直方图中bin的个数返回值：提取的特征点个数，若返回-1表明提取失败*/int _sift_features( IplImage* img, struct feature** feat, int intvls,double sigma, double contr_thr, int c urv_thr,int img_dbl, int descr_width, int des cr_hist_bins ){IplImage* init_img;//原图经初始化后的图像IplImage*** gauss_pyr, *** dog_pyr;//三级指针，高斯金字塔图像组，Do G金字塔图像组CvMemStorage* storage;//存储器CvSeq* features;//存储特征点的序列，序列中存放的是struct feature类型的指针int octvs, i, n = 0;//输入参数检查if( ! img )fatal_error( "NULL pointer error, %s, line %d", __FILE __, __LINE__ );if( ! feat )fatal_error( "NULL pointer error, %s, line %d", __FILE __, __LINE__ );/*步骤一：建立尺度空间，即建立高斯差分(DoG)金字塔dog_pyr将原图转换为32位灰度图并归一化，然后进行一次高斯平滑，并根据参数img_dbl决定是否将图像尺寸放大为原图的2倍*/init_img = create_init_img( img, img_dbl, sigma );//计算高斯金字塔的组数octvsoctvs = log( MIN( init_img->width, init_img->height ) ) / log (2) - 2;//为了保证连续性，在每一层的顶层继续用高斯模糊生成3幅图像，所以高斯金字塔每组//有intvls+3层，DOG金字塔每组有intvls+2层//建立高斯金字塔gauss_pyr，是一个octvs*(intvls+3)的图像数组gauss_pyr = build_gauss_pyr( init_img, octvs, intvls, sigma );//建立高斯差分(DoG)金字塔dog_pyr，是一个octvs*(intvls+2)的图像数组dog_pyr = build_dog_pyr( gauss_pyr, octvs, intvls );/*步骤二：在尺度空间中检测极值点，并进行精确定位和筛选创建默认大小的内存存储器*/storage = cvCreateMemStorage( 0 );//在尺度空间中检测极值点，通过插值精确定位，去除低对比度的点，去除边缘点，//返回检测到的特征点序列features = scale_space_extrema( dog_pyr, octvs, intvls, contr_ thr,curv_thr, storage );//计算特征点序列features中每个特征点的尺度calc_feature_scales( features, sigma, intvls );if( img_dbl ) //若设置了将图像放大为原图的2倍adjust_for_img_dbl( features );//将特征点序列中每个特征点的坐标减半//(当设置了将图像放大为原图的2倍时，特征点检测完之后调用)/*步骤三：特征点方向赋值，完成此步骤后，每个特征点有三个信息：位置、尺度、方向*///计算每个特征点的梯度直方图，找出其主方向，若一个特征点有不止一个主方向，将其分为两个特征点calc_feature_oris( features, gauss_pyr );/*步骤四：计算特征描述子*///计算特征点序列中每个特征点的特征描述子向量compute_descriptors( features, gauss_pyr, descr_width, descr_h ist_bins );//按特征点尺度的降序排列序列中的元素的顺序，feature_cmp是自定义的比较函数cvSeqSort( features, (CvCmpFunc)feature_cmp, NULL );//将CvSeq类型的特征点序列features转换为通用的struct feature类型的数组featn = features->total;//特征点个数*feat = calloc( n, sizeof(struct feature) );//分配控件//将序列features中的元素拷贝到数组feat中，返回数组指针给feat*feat = cvCvtSeqToArray( features, *feat, CV_WHOLE_SEQ );//释放特征点数组feat中所有特征点的feature_data成员，因为此成员中的数据在检测完特征//点后就没用了for( i = 0; i < n; i++ ){free( (*feat)[i].feature_data );(*feat)[i].feature_data = NULL;}//释放各种临时数据的存储空间cvReleaseMemStorage( &storage );cvReleaseImage( &init_img );release_pyr( &gauss_pyr, octvs, intvls + 3 );release_pyr( &dog_pyr, octvs, intvls + 2 );//返回检测到的特征点的个数return n;}/*将原图转换为32位灰度图并归一化，然后进行一次高斯平滑，并根据参数img_dbl决定是否将图像尺寸放大为原图的2倍参数：img：输入的原图像img_dbl：是否将图像放大为之前的两倍sigma：初始高斯平滑参数σ返回值：初始化完成的图像*/static IplImage* create_init_img( IplImage* img, int img_dbl, double s igma ){IplImage* gray, * dbl;float sig_diff;//调用函数，将输入图像转换为32位灰度图，并归一化gray = convert_to_gray32( img );if( img_dbl ) //若设置了将图像放大为原图的2倍{//将图像长宽扩展一倍时，便有了底-1层，该层尺度为：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 );//高斯平滑，高斯核在x,y方向上的标准差都是sig_diffcvSmooth( dbl, dbl, CV_GAUSSIAN, 0, 0, sig_diff, sig_di ff );cvReleaseImage( &gray );return dbl;}else//不用放大为原图的2倍{//计算第0层的尺度sig_diff = sqrt( sigma * sigma - SIFT_INIT_SIGMA * SIFT _INIT_SIGMA );//高斯平滑，高斯核在x,y方向上的标准差都是sig_diffcvSmooth( gray, gray, CV_GAUSSIAN, 0, 0, sig_diff, sig_ diff );return gray;}}/*将输入图像转换为32位灰度图,并进行归一化参数：img：输入图像，3通道8位彩色图(BGR)或8位灰度图返回值：32位灰度图*/static IplImage* convert_to_gray32( IplImage* img ){IplImage* gray8, * gray32;gray32 = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );//首先将原图转换为8位单通道图像if( img->nChannels == 1 )//若原图本身就是单通道，直接克隆原图gray8 = cvClone( img );else//若原图是3通道图像{gray8 = cvCreateImage( cvGetSize(img), IPL_DEPTH_8U, 1 );//创建8位单通道图像cvCvtColor( img, gray8, CV_BGR2GRAY );//将原图转换为8为单通道图像}//然后将8为单通道图像gray8转换为32位单通道图像，并进行归一化处理(除以255)cvConvertScale( gray8, gray32, 1.0 / 255.0, 0 );//释放临时图像cvReleaseImage( &gray8 );//返回32位单通道图像return gray32;}/*根据输入参数建立高斯金字塔参数：base：输入图像，作为高斯金字塔的基图像octvs：高斯金字塔的组数intvls：每组的层数sigma：初始尺度返回值：高斯金字塔，是一个octvs*(intvls+3)的图像数组*/static IplImage*** build_gauss_pyr( IplImage* base, int octvs,int i ntvls, double sigma ){IplImage*** gauss_pyr;//为了保证连续性，在每一层的顶层继续用高斯模糊生成3幅图像，所以高斯金字塔每组有intvls+3层，DOG金字塔每组有intvls+2层double* sig = calloc( intvls + 3, sizeof(double));//每层的sigma 参数数组double sig_total, sig_prev, k;int i, o;//为高斯金字塔gauss_pyr分配空间，共octvs个元素，每个元素是一组图像的首指针gauss_pyr = calloc( octvs, sizeof( IplImage** ) );//为第i组图像gauss_pyr[i]分配空间，共intvls+3个元素，每个元素是一个图像指针for( i = 0; i < octvs; i++ )gauss_pyr[i] = calloc( intvls + 3, sizeof( IplImage* ) );/* precompute Gaussian sigmas using the following formul a:sigma_{total}^2 = sigma_{i}^2 + sigma_{i-1}^2 *///计算每次高斯模糊的sigma参数sig[0] = sigma;//初始尺度k = pow( 2.0, 1.0 / intvls );for( i = 1; i < intvls + 3; i++ ){sig_prev = pow( k, i - 1 ) * sigma;//i-1层的尺度sig_total = sig_prev * k;//i层的尺度sig[i] = sqrt( sig_total * sig_total - sig_prev * sig_p rev );}//逐组逐层生成高斯金字塔for( o = 0; o < octvs; o++ )//遍历组for( i = 0; i < intvls + 3; i++ )//遍历层{if( o == 0 && i == 0 )//第0组，第0层，就是原图像gauss_pyr[o][i] = cvCloneImage(base);else if( i == 0 )//新的一组的首层图像是由上一组最后一层图像下采样得到gauss_pyr[o][i] = downsample( gauss_pyr [o-1][intvls] );else//对上一层图像进行高斯平滑得到当前层图像{ //创建图像gauss_pyr[o][i] = cvCreateImage( cvGetS ize(gauss_pyr[o][i-1]),IPL_DEPTH_32F, 1 );//对上一层图像gauss_pyr[o][i-1]进行参数为sig[i]的高斯平滑，得到当前层图像gauss_pyr[o][i]cvSmooth( gauss_pyr[o][i-1], gauss_pyr [o][i], CV_GAUSSIAN, 0, 0, sig[i], sig[i] );}}free( sig );//释放sigma参数数组return gauss_pyr;//返回高斯金字塔}/*对输入图像做下采样生成其四分之一大小的图像(每个维度上减半)，使用最近邻差值方法参数：img：输入图像返回值：下采样后的图像*/static IplImage* downsample( IplImage* img ){//下采样图像IplImage* smaller = cvCreateImage( cvSize(img->width / 2, img ->height / 2), img->depth, img->nChannels );cvResize( img, smaller, CV_INTER_NN );//尺寸变换return smaller;}/*通过对高斯金字塔中每相邻两层图像相减来建立高斯差分金字塔参数：gauss_pyr：高斯金字塔octvs：组数intvls：每组的层数返回值：高斯差分金字塔，是一个octvs*(intvls+2)的图像数组*/static IplImage*** build_dog_pyr( IplImage*** gauss_pyr, int octvs, in t intvls ){IplImage*** dog_pyr;int i, o;//为高斯差分金字塔分配空间，共octvs个元素，每个元素是一组图像的首指针dog_pyr = calloc( octvs, sizeof( IplImage** ) );//为第i组图像dog_pyr[i]分配空间，共(intvls+2)个元素，每个元素是一个图像指针for( i = 0; i < octvs; i++ )dog_pyr[i] = calloc( intvls + 2, sizeof(IplImage*) );//逐组逐层计算差分图像for( o = 0; o < octvs; o++ )//遍历组for( i = 0; i < intvls + 2; i++ )//遍历层{ //创建DoG金字塔的第o组第i层的差分图像dog_pyr[o][i] = cvCreateImage( cvGetSize(gauss _pyr[o][i]), IPL_DEPTH_32F, 1 );//将高斯金字塔的第o组第i+1层图像和第i层图像相减来得到DoG金字塔的第o组第i层图像cvSub( gauss_pyr[o][i+1], gauss_pyr[o][i], dog _pyr[o][i], NULL );}return dog_pyr;//返回高斯差分金字塔}/*在尺度空间中检测极值点，通过插值精确定位，去除低对比度的点，去除边缘点，返回检测到的特征点序列参数：dog_pyr：高斯差分金字塔octvs：高斯差分金字塔的组数intvls：每组的层数contr_thr：对比度阈值，针对归一化后的图像，用来去除不稳定特征cur_thr：主曲率比值的阈值，用来去除边缘特征storage：存储器返回值：返回检测到的特征点的序列*/static CvSeq* scale_space_extrema( IplImage*** dog_pyr, int octvs, int intvls,double contr_thr, int curv_thr, CvMemStorage* storage ){CvSeq* features;//特征点序列double prelim_contr_thr = 0.5 * contr_thr / intvls;//像素的对比度阈值struct feature* feat;struct detection_data* ddata;int o, i, r, c;//在存储器storage上创建存储极值点的序列，其中存储feature结构类型的数据features = cvCreateSeq( 0, sizeof(CvSeq), sizeof(struct featur e), storage );/*遍历高斯差分金字塔，检测极值点*///SIFT_IMG_BORDER指明边界宽度，只检测边界线以内的极值点for( o = 0; o < octvs; o++ )//第o组for( i = 1; i <= intvls; i++ )//遍i层for(r = SIFT_IMG_BORDER; r < dog_pyr[o][0]->hei ght-SIFT_IMG_BORDER; r++)//第r行for(c = SIFT_IMG_BORDER; c < dog_pyr[o] [0]->width-SIFT_IMG_BORDER; c++)//第c列//进行初步的对比度检查，只有当归一化后的像素值大于对比度阈值prelim_contr_thr时才继续检测此像素点是//否可能是极值//调用函数pixval32f获取图像dog_p yr[o][i]的第r行第c列的点的坐标值，然后调用ABS宏求其绝对值if( ABS( pixval32f( dog_pyr[o] [i], r, c ) ) > prelim_contr_thr )//通过在尺度空间中将一个像素点的值与其周围3*3*3邻域内的点比较来决定此点是否极值点(极大值或//极小都行)if( is_extremum( dog_py r, o, i, r, c ) )//若是极值点{//由于极值点的检测是在离散空间中进行的，所以检测到的极值点并不一定是真正意义上的极//值点//因为真正的极值点可能位于两个像素之间，而在离散空间中只能精确到坐标点精度上//通过亚像素级插值进行极值点精确定位(修正极值点坐标)，并去除低对比度的极值点，将修//正后的特征点组成feature结构返回feat = interp_ex tremum(dog_pyr, o, i, r, c, intvls, contr_thr);//返回值非空，表明此点已被成功修正if( feat ){//调用宏feat_detection_data来提取参数feat中的feature_data成员并转换为detec // tion_data类型的指针ddata = feat_detection_data( feat );//去除边缘响应，即通过计算主曲率比值判断某点是否边缘点,返回值为0表示不是边//缘点，可做特征点if( ! is _too_edge_like( dog_pyr[ddata->octv][ddata->intvl], ddata->r, ddat a->c, cu rv_thr ) ){ cvSeqPush( features, feat );//向特征点序列features末尾插入新检测到 //的特征点feat}elsefree( ddata );free( fe at );}}return features;//返回特征点序列}/*通过在尺度空间中将一个像素点的值与其周围3*3*3邻域内的点比较来决定此点是否极值点(极大值或极小都行)参数：dog_pyr：高斯差分金字塔octv：像素点所在的组intvl：像素点所在的层r：像素点所在的行c：像素点所在的列返回值：若指定的像素点是极值点(极大值或极小值)，返回1；否则返回0*/static int is_extremum( IplImage*** dog_pyr, int octv, int intvl, int r, int c ){//调用函数pixval32f获取图像dog_pyr[octv][intvl]的第r行第c列的点的坐标值float val = pixval32f( dog_pyr[octv][intvl], r, c );int i, j, k;//检查是否最大值if( val > 0 ){for( i = -1; i <= 1; i++ )//层for( j = -1; j <= 1; j++ )//行for( k = -1; k <= 1; k++ )//列if( val < pixval32f( dog_pyr[oc tv][intvl+i], r + j, c + k ) )return 0;}//检查是否最小值else{for( i = -1; i <= 1; i++ )//层for( j = -1; j <= 1; j++ )//行for( k = -1; k <= 1; k++ )//列if( val > pixval32f( dog_pyr[oc tv][intvl+i], r + j, c + k ) )return 0;}return 1;}/*通过亚像素级插值进行极值点精确定位(修正极值点坐标)，并去除低对比度的极值点，将修正后的特征点组成feature结构返回参数：dog_pyr：高斯差分金字塔octv：像素点所在的组intvl：像素点所在的层r：像素点所在的行c：像素点所在的列intvls：每组的层数contr_thr：对比度阈值，针对归一化后的图像，用来去除不稳定特征返回值：返回经插值修正后的特征点(feature类型)；若经有限次插值依然无法精确到理想情况或者该点对比度过低，返回NULL*/static struct feature* interp_extremum( IplImage*** dog_pyr, int octv, int intvl,int r, int c, int intvls, double contr_thr ){struct feature* feat;//修正后的特征点struct detection_data* ddata;//与特征检测有关的结构，存在feature结构的feature_data成员中double xi, xr, xc, contr;//xi,xr,xc分别为亚像素的intvl(层),row (y),col(x)方向上的增量(偏移量)int i = 0;//插值次数//SIFT_MAX_INTERP_STEPS指定了关键点的最大插值次数，即最多修正多少次，默认是5while( i < SIFT_MAX_INTERP_STEPS ){//进行一次极值点差值，计算σ(层方向,intvl方向)，y，x方向上的子像素偏移量(增量)interp_step( dog_pyr, octv, intvl, r, c, &xi, &xr, &xc );//若在任意方向上的偏移量大于0.5时，意味着差值中心已经偏移到它的临近点上，所以必须改变当前关键点的位置坐标if( ABS( xi ) < 0.5 && ABS( xr ) < 0.5 && ABS( xc ) < 0.5 )//若三方向上偏移量都小于0.5，表示已经够精确，则不用继续插值break;//修正关键点的坐标，x,y,σ三方向上的原坐标加上偏移量取整(四舍五入)c += cvRound( xc );//x坐标修正r += cvRound( xr );//y坐标修正intvl += cvRound( xi );//σ方向，即层方向//若坐标修正后超出范围，则结束插值，返回NULLif( 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_BORDE R ){return NULL;}i++;}//若经过SIFT_MAX_INTERP_STEPS次插值后还没有修正到理想的精确位置，则返回NULL，即舍弃此极值点if( i >= SIFT_MAX_INTERP_STEPS )return NULL;//计算被插值点的对比度：D + 0.5 * dD^T * Xcontr = interp_contr( dog_pyr, octv, intvl, r, c, xi, xr, xc );if( ABS( contr ) < contr_thr / intvls )//若该点对比度过小，舍弃，返回NULLreturn NULL;//为一个特征点feature结构分配空间并初始化，返回特征点指针feat = new_feature();//调用宏feat_detection_data来提取参数feat中的feature_data成员并转换为detection_data类型的指针ddata = feat_detection_data( feat );//将修正后的坐标赋值给特征点feat//原图中特征点的x坐标，因为第octv组中的图的尺寸比原图小2^octv倍，所以坐标值要乘以2^octvfeat->img_pt.x = feat->x = ( c + xc ) * pow( 2.0, octv );//原图中特征点的y坐标，因为第octv组中的图的尺寸比原图小2^octv倍，所以坐标值要乘以2^octvfeat->img_pt.y = feat->y = ( r + xr ) * pow( 2.0, octv );ddata->r = r;//特征点所在的行ddata->c = c;//特征点所在的列ddata->octv = octv;//高斯差分金字塔中，特征点所在的组ddata->intvl = intvl;//高斯差分金字塔中，特征点所在的组中的层ddata->subintvl = xi;//特征点在层方向(σ方向,intvl方向)上的亚像素偏移量return feat;//返回特征点指针}/*进行一次极值点差值，计算x，y，σ方向(层方向)上的子像素偏移量(增量)参数：dog_pyr：高斯差分金字塔octv：像素点所在的组intvl：像素点所在的层r：像素点所在的行c：像素点所在的列xi：输出参数，层方向上的子像素增量(偏移)xr：输出参数，y方向上的子像素增量(偏移)xc：输出参数，x方向上的子像素增量(偏移)*/static void interp_step( IplImage*** dog_pyr, int octv, int intvl, int r, int c,double* xi, double* xr, double* xc ){CvMat* dD, * H, * H_inv, X;double x[3] = { 0 };//在DoG金字塔中计算某点的x方向、y方向以及尺度方向上的偏导数，结果存放在列向量dD中dD = deriv_3D( dog_pyr, octv, intvl, r, c );//在DoG金字塔中计算某点的3*3海森矩阵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 );//X = - H^(-1) * dD，H的三个元素分别是x,y,σ方向上的偏移量(具体见SIF T算法说明)cvGEMM( H_inv, dD, -1, NULL, 0, &X, 0 );cvReleaseMat( &dD );cvReleaseMat( &H );cvReleaseMat( &H_inv );*xi = x[2];//σ方向(层方向)偏移量*xr = x[1];//y方向上偏移量*xc = x[0];//x方向上偏移量}/*在DoG金字塔中计算某点的x方向、y方向以及尺度方向上的偏导数参数：dog_pyr：高斯差分金字塔octv：像素点所在的组intvl：像素点所在的层r：像素点所在的行c：像素点所在的列返回值：返回3个偏导数组成的列向量{ dI/dx, dI/dy, dI/ds }^T*/static CvMat* deriv_3D( IplImage*** dog_pyr, int octv, int intvl, int r, int c ){CvMat* dI;double dx, dy, ds;//求差分来代替偏导，这里是用的隔行求差取中值的梯度计算方法//求x方向上的差分来近似代替偏导数dx = ( pixval32f( dog_pyr[octv][intvl], r, c+1 ) -pixval32f( dog_pyr[octv][intvl], r, c-1 ) ) / 2.0;//求y方向上的差分来近似代替偏导数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;//组成列向量dI = cvCreateMat( 3, 1, CV_64FC1 );cvmSet( dI, 0, 0, dx );cvmSet( dI, 1, 0, dy );cvmSet( dI, 2, 0, ds );return dI;}/*在DoG金字塔中计算某点的3*3海森矩阵/ Ixx Ixy Ixs \| Ixy Iyy Iys |\ Ixs Iys Iss /参数：dog_pyr：高斯差分金字塔octv：像素点所在的组intvl：像素点所在的层r：像素点所在的行c：像素点所在的列返回值：返回3*3的海森矩阵*/static CvMat* hessian_3D( IplImage*** dog_pyr, int octv, int intvl, in t r, int c ){CvMat* H;double v, dxx, dyy, dss, dxy, dxs, dys;v = pixval32f( dog_pyr[octv][intvl], r, c );//该点的像素值//用差分近似代替倒数(具体公式见各种梯度的求法)//dxx = f(i+1,j) - 2f(i,j) + f(i-1,j)//dyy = f(i,j+1) - 2f(i,j) + f(i,j-1)dxx = ( pixval32f( dog_pyr[octv][intvl], r, c+1 ) +pixval32f( dog_pyr[octv][intvl], r, c-1 ) - 2 * v );dyy = ( pixval32f( dog_pyr[octv][intvl], r+1, c ) +pixval32f( dog_pyr[octv][intvl], r-1, c ) - 2 * v );dss = ( pixval32f( dog_pyr[octv][intvl+1], r, c ) +pixval32f( dog_pyr[octv][intvl-1], r, c ) - 2 * v );dxy = ( pixval32f( dog_pyr[octv][intvl], r+1, c+1 ) -pixval32f( dog_pyr[octv][intvl], r+1, c-1 ) -pixval32f( dog_pyr[octv][intvl], r-1, c+1 ) +pixval32f( dog_pyr[octv][intvl], r-1, c-1 ) ) / 4.0;dxs = ( pixval32f( dog_pyr[octv][intvl+1], r, c+1 ) -pixval32f( dog_pyr[octv][intvl+1], r, c-1 ) -pixval32f( dog_pyr[octv][intvl-1], r, c+1 ) +pixval32f( dog_pyr[octv][intvl-1], r, c-1 ) ) / 4.0;dys = ( pixval32f( dog_pyr[octv][intvl+1], r+1, c ) -pixval32f( dog_pyr[octv][intvl+1], r-1, c ) -pixval32f( dog_pyr[octv][intvl-1], r+1, c ) +pixval32f( dog_pyr[octv][intvl-1], r-1, c ) ) / 4.0;//组成海森矩阵H = cvCreateMat( 3, 3, CV_64FC1 );cvmSet( H, 0, 0, dxx );cvmSet( H, 0, 1, dxy );cvmSet( H, 0, 2, dxs );cvmSet( H, 1, 0, dxy );cvmSet( H, 1, 1, dyy );cvmSet( H, 1, 2, dys );cvmSet( H, 2, 0, dxs );cvmSet( H, 2, 1, dys );cvmSet( H, 2, 2, dss );return H;}/*计算被插值点的对比度：D + 0.5 * dD^T * X参数：dog_pyr：高斯差分金字塔octv：像素点所在的组intvl：像素点所在的层r：像素点所在的行c：像素点所在的列xi：层方向上的子像素增量xr：y方向上的子像素增量xc：x方向上的子像素增量返回值：插值点的对比度*/static double interp_contr( IplImage*** dog_pyr, int octv, int intvl, int r,int c, double xi, dou ble xr, double xc ){CvMat* dD, X, T;double t[1], x[3] = { xc, xr, xi };//偏移量组成的列向量X，其中是x，y，σ三方向上的偏移量cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );//矩阵乘法的结果T，是一个数值cvInitMatHeader( &T, 1, 1, CV_64FC1, t, CV_AUTOSTEP );//在DoG金字塔中计算某点的x方向、y方向以及尺度方向上的偏导数，结果存放在列向量dD中dD = deriv_3D( dog_pyr, octv, intvl, r, c );//矩阵乘法：T = dD^T * XcvGEMM( dD, &X, 1, NULL, 0, &T, CV_GEMM_A_T );cvReleaseMat( &dD );//返回计算出的对比度值：D + 0.5 * dD^T * X (具体公式推导见SIFT算法说明)return pixval32f( dog_pyr[octv][intvl], r, c ) + t[0] * 0.5;}/*为一个feature结构分配空间并初始化返回值：初始化完成的feature结构的指针*/static struct feature* new_feature( void ){struct feature* feat;//特征点指针struct detection_data* ddata;//与特征检测相关的结构feat = malloc( sizeof( struct feature ) );//分配空间memset( feat, 0, sizeof( struct feature ) );//清零ddata = malloc( sizeof( struct detection_data ) );memset( ddata, 0, sizeof( struct detection_data ) );feat->feature_data = ddata;//将特征检测相关的结构指针赋值给特征点的f eature_data成员feat->type = FEATURE_LOWE;//默认是LOWE类型的特征点return feat;}/*去除边缘响应，即通过计算主曲率比值判断某点是否边缘点参数：dog_img：此特征点所在的DoG图像r：特征点所在的行c：特征点所在的列cur_thr：主曲率比值的阈值，用来去除边缘特征返回值：0：此点是非边缘点；1：此点是边缘点*/static int is_too_edge_like( IplImage* dog_img, int r, int c, int curv _thr ){double d, dxx, dyy, dxy, tr, det;/*某点的主曲率与其海森矩阵的特征值成正比，为了避免直接计算特征值，这里只考虑特征值的比值可通过计算海森矩阵的迹tr(H)和行列式det(H)来计算特征值的比值设a是海森矩阵的较大特征值，b是较小的特征值，有a = r*b，r是大小特征值的比值tr(H) = a + b; det(H) = a*b;tr(H)^2 / det(H) = (a+b)^2 / ab = (r+1)^2/rr越大，越可能是边缘点；伴随r的增大，(r+1)^2/r 的值也增大，所以可通过(r +1)^2/r 判断主曲率比值是否满足条件*//* principal curvatures are computed using the trace and det of Hessian */d = pixval32f(dog_img, r, c);//调用函数pixval32f获取图像dog_img 的第r行第c列的点的坐标值//用差分近似代替偏导，求出海森矩阵的几个元素值/* / dxx dxy \\ dxy dyy / */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;//海森矩阵的行列式//若行列式为负，表明曲率有不同的符号，去除此点/* negative determinant -> curvatures have different signs; re ject feature */if( det <= 0 )return 1;//返回1表明此点是边缘点//通过式子：(r+1)^2/r 判断主曲率的比值是否满足条件，若小于阈值，表明不是边缘点if( tr * tr / det < ( curv_thr + 1.0 )*( curv_thr + 1.0 ) / curv _thr )return 0;//不是边缘点return 1;//是边缘点}/*计算特征点序列中每个特征点的尺度参数：features：特征点序列sigma：初始高斯平滑参数，即初始尺度intvls：尺度空间中每组的层数*/static void calc_feature_scales( CvSeq* features, double sigma, int in tvls ){struct feature* feat;struct detection_data* ddata;double intvl;int i, n;n = features->total;//总的特征点个数//遍历特征点for( i = 0; i < n; i++ ){//调用宏，获取序列features中的第i个元素，并强制转换为struct feature类型feat = CV_GET_SEQ_ELEM( struct feature, features, i );//调用宏feat_detection_data来提取参数feat中的feature_da ta成员并转换为detection_data类型的指针ddata = feat_detection_data( feat );//特征点所在的层数ddata->intvl加上特征点在层方向上的亚像素偏移量，得到特征点的较为精确的层数intvl = ddata->intvl + ddata->subintvl;//计算特征点的尺度(公式见SIFT算法说明)，并赋值给scl成员feat->scl = sigma * pow( 2.0, ddata->octv + intvl / int vls );//计算特征点所在的组的尺度，给detection_data的scl_octv成员赋值ddata->scl_octv = sigma * pow( 2.0, intvl / intvls );}}/*将特征点序列中每个特征点的坐标减半(当设置了将图像放大为原图的2倍时，特征点检测完之后调用)参数：features：特征点序列*/static void adjust_for_img_dbl( CvSeq* features ){struct feature* feat;int i, n;n = features->total;//总的特征点个数//遍历特征点for( i = 0; i < n; i++ )。
sift算法代码
sift算法代码
以下是SIFT算法的伪代码：
1. 尺度空间的构建:
a. 对图像进行高斯滤波，生成不同尺度的图像金字塔。
b. 使用拉普拉斯金字塔计算图像的尺度空间。
2. 关键点检测:
a. 对每个尺度空间进行极值点检测，找到图像中的候选关键点。
b. 对候选关键点进行精确定位，排除不稳定的关键点。
3. 方向计算:
a. 对每个关键点邻域内的像素计算梯度幅值和方向。
b. 构建梯度方向直方图，找到主要方向。
4. 描述子生成:
a. 将关键点邻域划分为若干个子区域。
b. 在每个子区域内计算梯度方向直方图，得到局部特征描述子。
5. 关键点匹配:
a. 使用特征描述子之间的距离进行关键点匹配。
b. 使用最近邻筛选和次近邻比率筛选进行匹配点的筛选。
其中，步骤2-4都是针对单个关键点进行操作的，可以进行并
行计算。
具体的代码实现可以参考开源的SIFT库，如OpenCV中的SIFT实现。
sift算法matlab复杂代码
一、介绍SIFT算法SIFT（Scale-Invariant Feature Transform）算法是一种用于图像处理和计算机视觉领域的特征提取算法，由David Lowe在1999年提出。
SIFT算法具有旋转、尺度、光照等方面的不变性，能够对图像进行稳健的特征点提取，被广泛应用于物体识别、图像匹配、图像拼接、三维重建等领域。
二、SIFT算法原理SIFT算法的主要原理包括尺度空间极值点检测、关键点定位、关键点方向确定、关键点描述等步骤。
其中，尺度空间极值点检测通过高斯差分金字塔来检测图像中的极值点，关键点定位则利用DoG响应函数进行关键点细化，关键点方向确定和关键点描述部分则通过梯度方向直方图和关键点周围区域的梯度幅度信息来完成。
三、使用Matlab实现SIFT算法在Matlab中实现SIFT算法，需要对SIFT算法的每个步骤进行详细的编程和调试。
需要编写代码进行图像的高斯金字塔和高斯差分金字塔的构建，计算尺度空间极值点，并进行关键点定位。
需要实现关键点的方向确定和描述子生成的算法。
将所有步骤整合在一起，完成SIFT算法的整体实现。
四、SIFT算法复杂代码的编写SIFT算法涉及到的步骤较多，需要编写复杂的代码来实现。
在编写SIFT算法的Matlab代码时，需要考虑到算法的高效性、可扩展性和稳定性。
具体来说，需要注意以下几点：1. 高斯差分金字塔和高斯金字塔的构建：在构建高斯差分金字塔时，需要编写代码实现图像的高斯滤波和图像的降采样操作，以得到不同尺度空间的图像。
还需要实现高斯差分金字塔的构建，以检测图像中的极值点。
2. 尺度空间极值点检测：在检测图像中的极值点时，需要编写代码实现对高斯差分金字塔的极值点检测算法，以找到图像中的潜在关键点。
3. 关键点的定位：关键点定位阶段需要编写代码实现对尺度空间极值点的精确定位，消除低对比度点和边缘响应点，并进行关键点的精细化操作。
4. 关键点的方向确定和描述子生成：在这一步骤中，需要编写代码实现对关键点周围区域的梯度幅度信息的计算和关键点方向的确定，以及生成关键点的描述子。
SIFT算法实现及代码详解
SIFT算法实现及代码详解SIFT（Scale-Invariant Feature Transform）是一种计算机视觉算法，用于在图像中寻找局部特征。
由于SIFT算法具有尺度不变性和旋转不变性等优点，因此广泛应用于目标检测、图像拼接、图像匹配等领域。
1.尺度空间构建：通过对原始图像进行多次高斯模糊，得到图像的尺度空间。
2.高斯差分金字塔构建：对尺度空间进行差分操作，得到高斯差分金字塔。
3.极值点检测：在高斯差分金字塔中检测局部极值点，用于确定关键点的位置。
4.关键点精确定位：通过对关键点进行插值，得到更精确的关键点位置。
5.方向分配：为每个关键点分配主方向，用于后续的旋转不变性。
6.关键点描述：对每个关键点周围的区域计算特征向量，用于描述关键点的特征。
7.特征匹配：通过比较特征向量的相似度，实现图像间的特征匹配。
下面是SIFT算法的Python实现（使用OpenCV库）：```pythonimport cv2def SIFT(image):#尺度空间构建gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)sift = cv2.xfeatures2d.SIFT_create# keypoints为检测到的关键点列表，descriptors为关键点的特征描述向量return keypoints, descriptors#读取图像image1 = cv2.imread('image1.png')image2 = cv2.imread('image2.png')#提取图像特征keypoints1, descriptors1 = SIFT(image1)keypoints2, descriptors2 = SIFT(image2)#特征匹配bf = cv2.BFMatcher(cv2.NORM_L2)matches = bf.match(descriptors1, descriptors2)matches = sorted(matches, key = lambda x:x.distance)#绘制匹配结果result = cv2.drawMatches(image1, keypoints1, image2, keypoints2, matches, None, flags=2)cv2.imshow('Result', result)cv2.waitKey(0)cv2.destroyAllWindows```上述代码中，首先定义了一个`SIFT`函数用于提取图像的特征。
SIFT算法C语言逐步实现详解
SIFT算法C语言逐步实现详解（上）引言：在我写的关于sift算法的前倆篇文章里头，已经对sift算法有了初步的介绍：九、图像特征提取与匹配之SIFT算法，而后在：九（续）、sift算法的编译与实现里，我也简单记录下了如何利用opencv，gsl等库编译运行sift程序。
但据一朋友表示，是否能用c语言实现sift算法，同时，尽量不用到opencv，gsl等第三方库之类的东西。
而且，Rob Hess维护的sift 库，也不好懂，有的人根本搞不懂是怎么一回事。
那么本文，就教你如何利用c语言一步一步实现sift算法，同时，你也就能真正明白sift算法到底是怎么一回事了。
ok，先看一下，本程序最终运行的效果图，sift 算法分为五个步骤（下文详述），对应以下第二--第六幅图：sift算法的步骤要实现一个算法，首先要完全理解这个算法的原理或思想。
咱们先来简单了解下，什么叫sift算法：sift，尺度不变特征转换，是一种电脑视觉的算法用来侦测与描述影像中的局部性特征，它在空间尺度中寻找极值点，并提取出其位置、尺度、旋转不变量，此算法由David Lowe 在1999年所发表，2004年完善总结。
所谓，Sift算法就是用不同尺度（标准差）的高斯函数对图像进行平滑，然后比较平滑后图像的差别，差别大的像素就是特征明显的点。
以下是sift算法的五个步骤：一、建立图像尺度空间(或高斯金字塔)，并检测极值点首先建立尺度空间，要使得图像具有尺度空间不变形，就要建立尺度空间，sift算法采用了高斯函数来建立尺度空间，高斯函数公式为：上述公式G(x,y,e)，即为尺度可变高斯函数。
而，一个图像的尺度空间L（x,y,e) ,定义为原始图像I(x,y)与上述的一个可变尺度的2维高斯函数G(x,y,e) 卷积运算。
即，原始影像I(x，y)在不同的尺度e下，与高斯函数G(x,y,e)进行卷积，得到L(x,y,e)，如下：以上的（x，y）是空间坐标，e，是尺度坐标，或尺度空间因子，e的大小决定平滑程度，大尺度对应图像的概貌特征，小尺度对应图像的细节特征。
sift算法的原理和步骤
sift算法的原理和步骤SIFT算法的原理和步骤SIFT算法是一种用于图像特征提取的算法，它能够从图像中提取出具有独特性、稳定性和可重复性的关键点，用于图像匹配、目标跟踪等任务。
本文将介绍SIFT算法的原理和步骤。
一、原理1. 尺度空间尺度空间是指同一物体在不同尺度下的表现形式。
SIFT算法采用高斯金字塔来实现尺度空间的构建，即将原始图像不断缩小并平滑处理，得到一系列模糊程度不同的图像。
2. 关键点检测在尺度空间中，SIFT算法采用DoG（Difference of Gaussian）来检测关键点。
DoG是指两个不同尺寸的高斯滤波器之间的差值，可以有效地提取出具有高斯拉普拉斯变换极值点（LoG）特征的区域。
3. 方向确定对于每个关键点，在其周围区域内计算梯度幅值和方向，并统计梯度直方图。
最终确定该关键点最显著的梯度方向作为其主方向。
4. 描述子生成以关键点为中心，生成一个16x16的方形区域，并将其分为4x4的小块。
对于每个小块，计算其内部像素的梯度方向直方图，并将其串联成一个128维的向量，作为该关键点的描述子。
5. 匹配通过计算不同图像之间的关键点描述子之间的距离来进行匹配。
采用最近邻法（Nearest Neighbor）和次近邻法（Second Nearest Neighbor）来进行筛选，从而得到最终的匹配结果。
二、步骤1. 构建高斯金字塔对于原始图像，采用高斯滤波器进行平滑处理，并将其缩小一定比例后再次平滑处理，得到一系列不同尺度下的图像。
这些图像构成了高斯金字塔。
2. 构建DoG金字塔在高斯金字塔中，相邻两层之间做差得到一组DoG金字塔。
通过在DoG金字塔上寻找局部极值点来检测关键点。
3. 确定关键点主方向对于每个关键点，在其周围区域内计算梯度幅值和方向，并统计梯度直方图。
最终确定该关键点最显著的梯度方向作为其主方向。
4. 生成描述子以关键点为中心，生成一个16x16的方形区域，并将其分为4x4的小块。
sift matlab 代码
sift matlab 代码Sift Matlab 代码Sift算法（Scale-invariant feature transform）是一种用于图像处理和计算机视觉中的特征提取算法。
Matlab是一种常用的编程语言和软件工具，用于实现各种算法和程序。
本文将介绍如何使用Matlab 编写Sift算法的代码，并对其原理和实现进行详细说明。
一、Sift算法原理Sift算法是一种基于局部特征的图像匹配算法，其主要思想是在图像中寻找关键点（keypoints），并对这些关键点进行描述，以便在不同图像之间进行匹配。
Sift算法具有尺度不变性和旋转不变性，能够稳定地提取图像的特征。
Sift算法的主要步骤包括尺度空间极值检测、关键点定位、关键点描述和关键点匹配等。
在尺度空间极值检测中，算法会在图像的不同尺度下检测局部极值点，这些点被认为是关键点的候选。
然后，通过对这些候选点进行精确定位和剔除不稳定点，最终得到真正的关键点。
接着，算法会对每个关键点周围的图像区域进行描述，生成描述子向量。
最后，通过比较不同图像的关键点描述子，实现图像匹配。
二、Matlab实现Sift算法在Matlab中实现Sift算法，可以使用现成的开源库或者自己编写代码。
一种常见的做法是使用vlfeat工具包，该工具包提供了Sift算法的Matlab接口，并包含了各种图像处理和特征提取的函数。
需要安装vlfeat工具包，并将其添加到Matlab的路径中。
然后，可以使用vl_sift函数来提取图像的Sift特征。
该函数会返回关键点的位置、尺度、方向以及描述子等信息，可以根据需要进行进一步处理和分析。
除了使用vlfeat工具包，也可以根据Sift算法的原理，编写自己的Sift代码。
在Matlab中，可以通过图像梯度计算、高斯金字塔构建和关键点描述等步骤，实现Sift算法的各个功能模块。
三、总结本文介绍了Sift算法的原理和在Matlab中的实现方法。
SIFT算法原理：SIFT算法详细介绍
前面们介绍了Harris和Shi-Tomasi角检测算法，这两种算法具有旋转不变性，但不具有尺度不变性，以下图为例，在左侧小图中可以检测到角，但图像被放后，在使用同样的窗口，就检测不到角了。
所以，们来介绍一种计算机视觉的算法，尺度不变特征转换即SIFT(Scale-invariantfeaturetransform)。
它用来侦测与描述影像中的局部性特征，它在空间尺度中寻找极值，并提取出其位置、尺度、旋转不变量，此算法由DavidLowe在1999年所发表，2004年完善总结。
应用范围包含物体辨识、机器人地图感知与导航、影像缝合、3D模型建立、手势辨识、影像追踪和动作比对等领域。
SIFT算法的实质在不同的尺度空间上查找关键(特征)，并计算出关键的方向。
SIFT 所查找到的关键一些十分突出，不会因光照，仿变换和噪音等因素而变化的，如角、边缘、暗区的亮及亮区的暗等。
1.1基本流程Lowe将SIFT算法分解为如下四步：尺度空间极值检测：搜索所有尺度上的图像位置。
通过高斯差分函数来识别潜在的对于尺度和旋转不变的关键。
关键定位：在每个候选的位置上，通过一个拟合精细的模型来确定位置和尺度。
关键的选择依据于它们的稳定程度。
关键方向确定：基于图像局部的梯度方向，分配给每个关键位置一个或多个方向。
所有后面的对图像数据的操作都相对于关键的方向、尺度和位置进行变换，从而保证了对于这些变换的不变性。
关键描述：在每个关键周围的邻域内，在选定的尺度上测量图像局部的梯度。
这些梯度作为关键的描述符，它允许比较的局部形状的变形或光照变化。
们就沿着Lowe的步骤，对SIFT算法的实现过程进行介绍：1.2尺度空间极值检测在不同的尺度空间不能使用相同的窗口检测极值，对小的关键使用小的窗口，对的关键使用的窗口，为了达到上述目的，们使用尺度空间滤波器。
高斯核可以产生多尺度空间的核函数。
-《Scale-spacetheory:Abasictoolforanalysingstructuresatdifferentscales》。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
经典算法SIFT实现即代码解释：以下便是sift源码库编译后的效果图：为了给有兴趣实现sift算法的朋友提供个参考，特整理此文如下。
要了解什么是sift算法，请参考：九、图像特征提取与匹配之SIFT算法。
ok，咱们下面，就来利用Rob Hess维护的sift 库来实现sift算法：首先，请下载Rob Hess维护的sift 库：/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/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加入系统变量"，打上“勾”。
（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/includeC:/Program Files/OpenCV/cv/includeC:/Program Files/OpenCV/cvaux/includeC:/Program Files/OpenCV/ml/includeC:/Program Files/OpenCV/otherlibs/highguiC:/Program Files/OpenCV/otherlibs/cvcam/include最后选择source files，在下方填入路径：C:/Program Files/OpenCV/cv/srcC:/Program Files/OpenCV/cxcore/srcC:/Program Files/OpenCV/cvaux/srcC:/Program Files/OpenCV/otherlibs/highguiC:/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也是一个库，也需要下载：/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目录复制到C:/VC6.0/Bin下；lib目录下的所有.lib文件全部复制到C:/VC6.0/Lib下。
然后，在tools－options－directories中，将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 file"odbccp32.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是直方图中的条数[cpp]view 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.//计算描述子，包括计算二维方向直方图和转换其为特征描述子45. compute_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<feature *>( calloc( n, sizeof(struct feature) ) );51. *feat = static_cast<feature *>( cvCvtSeqToArray( features, *feat, CV_WHOLE_SEQ ) );52.53.54.55.56.for( i = 0; i < n; i++ )57. {58. free( (*feat)[i].feature_data );59. (*feat)[i].feature_data = NULL;60. }61.62. cvReleaseMemStorage( &storage );63. cvReleaseImage( &init_img );64. release_pyr( &gauss_pyr, octvs, intvls + 3 );65. release_pyr( &dog_pyr, octvs, intvls + 2 );66.return n;67.}(1)初始化图像输入参数：这里不需要解释了该函数主要用来初始化图像，转换图像为32位灰度图以及进行高斯模糊。