SIFT算法实现及代码详解

#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算法的伪代码：
1. 尺度空间的构建:
a. 对图像进行高斯滤波，生成不同尺度的图像金字塔。

b. 使用拉普拉斯金字塔计算图像的尺度空间。

2. 关键点检测:
a. 对每个尺度空间进行极值点检测，找到图像中的候选关键点。

b. 对候选关键点进行精确定位，排除不稳定的关键点。

3. 方向计算:
a. 对每个关键点邻域内的像素计算梯度幅值和方向。

b. 构建梯度方向直方图，找到主要方向。

4. 描述子生成:
a. 将关键点邻域划分为若干个子区域。

b. 在每个子区域内计算梯度方向直方图，得到局部特征描述子。

5. 关键点匹配:
a. 使用特征描述子之间的距离进行关键点匹配。

b. 使用最近邻筛选和次近邻比率筛选进行匹配点的筛选。

其中，步骤2-4都是针对单个关键点进行操作的，可以进行并
行计算。

具体的代码实现可以参考开源的SIFT库，如OpenCV中的SIFT实现。

一、介绍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（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算法的前倆篇文章里头，已经对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算法的原理和步骤。

一、原理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算法（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中的实现方法。

前面们介绍了Harris和Shi-Tomasi角检测算法，这两种算法具有旋转不变性，但不具有尺度不变性，以下图为例，在左侧小图中可以检测到角，但图像被放后，在使用同样的窗口，就检测不到角了。

所以，们来介绍一种计算机视觉的算法，尺度不变特征转换即SIFT(Scale-invariantfeaturetransform)。

它用来侦测与描述影像中的局部性特征，它在空间尺度中寻找极值，并提取出其位置、尺度、旋转不变量，此算法由DavidLowe在1999年所发表，2004年完善总结。

应用范围包含物体辨识、机器人地图感知与导航、影像缝合、3D模型建立、手势辨识、影像追踪和动作比对等领域。

SIFT算法的实质在不同的尺度空间上查找关键(特征)，并计算出关键的方向。

SIFT 所查找到的关键一些十分突出，不会因光照，仿变换和噪音等因素而变化的，如角、边缘、暗区的亮及亮区的暗等。

1.1基本流程Lowe将SIFT算法分解为如下四步：尺度空间极值检测：搜索所有尺度上的图像位置。

通过高斯差分函数来识别潜在的对于尺度和旋转不变的关键。

关键定位：在每个候选的位置上，通过一个拟合精细的模型来确定位置和尺度。

关键的选择依据于它们的稳定程度。

关键方向确定：基于图像局部的梯度方向，分配给每个关键位置一个或多个方向。

所有后面的对图像数据的操作都相对于关键的方向、尺度和位置进行变换，从而保证了对于这些变换的不变性。

关键描述：在每个关键周围的邻域内，在选定的尺度上测量图像局部的梯度。

这些梯度作为关键的描述符，它允许比较的局部形状的变形或光照变化。

们就沿着Lowe的步骤，对SIFT算法的实现过程进行介绍：1.2尺度空间极值检测在不同的尺度空间不能使用相同的窗口检测极值，对小的关键使用小的窗口，对的关键使用的窗口，为了达到上述目的，们使用尺度空间滤波器。

高斯核可以产生多尺度空间的核函数。

-《Scale-spacetheory:Abasictoolforanalysingstructuresatdifferentscales》。

SIFT特征提取算法SIFT（Scale-Invariant Feature Transform）特征提取算法是一种用于图像的局部特征分析的算法。

它能够提取出图像中的关键点，并对这些关键点进行描述，从而可以用于图像匹配、物体识别等应用领域。

本文将详细介绍SIFT算法的原理和过程。

1.尺度空间构建SIFT算法首先通过使用高斯滤波器来构建图像的尺度空间，以便在不同尺度下检测关键点。

高斯滤波器可以通过一系列的高斯卷积操作实现，每次卷积之后对图像进行下采样（降低分辨率），得到不同尺度的图像。

2.关键点检测在尺度空间构建完成后，SIFT算法使用差分运算来检测关键点。

差分运算可以通过对图像进行高斯平滑操作来实现，然后计算相邻尺度之间的差分图像。

对差分图像进行极值检测，即寻找局部最大和最小值的像素点，这些像素点就是图像中的关键点。

3.关键点精确定位关键点的精确定位是通过拟合关键点周围的局部图像来实现的。

SIFT算法使用了一种高度鲁棒的方法，即利用关键点周围梯度的方向和大小来进行拟合。

具体来说，SIFT算法在关键点周围计算图像的梯度幅值和方向，并构建梯度直方图。

然后通过在梯度直方图中寻找局部极值来确定关键点的方向。

4.关键点描述关键点的描述是为了提取关键点周围的特征向量，用于后续的匹配和识别。

SIFT算法使用了一种局部特征描述算法，即将关键点周围的图像区域划分为小的子区域，并计算每个子区域的梯度方向直方图。

然后将这些直方图组合起来，构成一个维度较高的特征向量。

5.特征向量匹配在完成关键点描述之后，SIFT算法使用一种近似的最近邻方法来进行特征向量的匹配。

具体来说，使用KD树或者暴力匹配的方法来寻找两幅图像中最相似的特征向量。

通过计算特征向量之间的距离，可以找到最相似的匹配对。

6.尺度不变性SIFT算法具有尺度不变性的特点，即对于图像的缩放、旋转和视角变化等变换具有较好的鲁棒性。

这是因为在特征提取的过程中，SIFT算法对图像进行了多尺度的分析，并利用了关键点周围的梯度信息进行描述。

SIFT算法实现C语言SIFT算法（Scale-Invariant Feature Transform）是一种用于图像特征提取和匹配的计算机视觉算法。

它的主要思想是通过检测和描述图像中的关键点，提取出具有尺度不变性和旋转不变性的特征点，从而适应不同尺度和旋转的图像匹配。

本文将介绍如何用C语言实现SIFT算法。

1.数据结构首先，我们需要定义一些用于存储图像信息和特征点信息的数据结构。

以下是一个简单的数据结构定义：```ctypedef structint x;int y;float scale;float orientation;float descriptor[128];} Keypoint;typedef structunsigned char *data;int width;int height;} Image;```在这个数据结构中，Keypoint用于存储关键点的位置、尺度、方向和描述子，Image用于存储图像的像素数据、宽度和高度。

2.关键点检测关键点检测是SIFT算法的第一步。

我们可以使用一个简单的角点检测算法，如Harris角点检测，来找到图像中具有显著度较高的角点。

以下是一个简单的关键点检测函数的实现：```cvoid detectKeypoints(Image *input, Keypoint *output, int*numKeypoints)// Harris corner detection implementation//...```这个函数接受一个图像作为输入，并在图像上检测并输出关键点。

输出参数包括关键点的数组output和关键点的数量numKeypoints。

3.关键点描述关键点描述是SIFT算法的第二步。

对于每个关键点，我们需要计算一个描述子，用于描述关键点周围的图像区域。

以下是一个简单的关键点描述函数的实现：```c//...```这个函数接受一个图像和关键点数组作为输入，并为每个关键点计算描述子。

一、介绍SIFT特征提取算法SIFT（Scale-Invariant Feature Transform）是一种用于图像处理和计算机视觉领域的特征提取算法，由David Lowe在1999年提出。

SIFT算法能够在不同尺度、旋转和光照条件下检测和描述图像中的关键特征点。

二、SIFT特征提取的原理1. 尺度空间极值检测SIFT算法首先通过高斯滤波器构建图像的尺度空间金字塔，在不同尺度上检测图像中的极值点，用于定位关键特征点。

2. 关键点定位在尺度空间金字塔上定位极值点，并通过尺度空间的拟合插值定位关键点的位置和尺度。

3. 方向确定对关键点周围的梯度方向进行统计，选择主要梯度方向作为关键点的方向。

4. 关键点描述基于关键点周围的梯度幅值和方向构建特征向量，用于描述关键点的外观特征。

三、MATLAB实现SIFT特征提取代码MATLAB提供了丰富的图像处理工具包，其中也包括对SIFT算法的支持。

下面给出MATLAB实现SIFT特征提取的简单示例代码：```matlab读入图像I = imread('image.jpg');I = single(rgb2gray(I));提取SIFT特征点[f, d] = vl_sift(I);显示特征点imshow(I);h1 = vl_plotframe(f);h2 = vl_plotsiftdescriptor(d, f);set(h1,'color','k','linewidth',3);set(h2,'color','g')；```上述代码中，首先使用`imread`读入一张图像，并将其转化为灰度图像。

然后使用`vl_sift`函数提取图像中的SIFT特征点，其中`f`为特征点的位置和尺度，`d`为特征点的描述子。

最后使用`vl_plotframe`和`vl_plotsiftdescriptor`函数将特征点和特征描述子可视化出来。

SIFT算法详解SIFT（Scale-Invariant Feature Transform）算法是一种用于计算图像特征的算法，最早由David Lowe于1999年提出，是一种用于在不同尺度和旋转下具有不变性的特征描述算法。

1.关键点检测：SIFT算法首先对图像进行多尺度空间的检测，即在图像中寻找尺度空间极值点。

为了实现尺度不变性，SIFT算法采用了高斯金字塔的方法来对图像进行多尺度模糊处理。

金字塔的每一层可以通过对上一层图像进行下采样和卷积来生成。

然后，在每一层金字塔上使用一种拉普拉斯算子来计算尺度空间的极值点，这些极值点包括图像的边缘点、角点和斑点等。

极值点的检测是通过对每个像素点的周围像素点进行梯度计算和方向估计来实现的。

2.关键点描述：在关键点检测之后，SIFT算法通过计算每个关键点的局部图像梯度直方图来对关键点进行描述。

首先，将关键点周围的领域分为若干个子区域，一般为16个，每个子区域可以通过将关键点周围的区域划分为4×4个小区域来实现。

然后，对每个小区域中的像素点计算其梯度幅值和方向，并将其分到相应的方向直方图中。

通过对所有子区域的梯度方向直方图进行拼接，就可以得到整个关键点的描述子。

描述子的维度一般为128维，特征向量的每个维度包含了关键点的领域中的梯度方向信息。

3.特征匹配：在关键点描述之后，SIFT算法使用一种基于特征向量距离的匹配方法来实现图像特征点的匹配。

常用的方法是计算两个特征向量之间的欧式距离，并使用一个阈值来筛选出较为相似的特征点。

为了提高匹配的准确性和稳定性，SIFT算法使用了一种自适应的阈值筛选方法，通过计算两个特征向量距离的比值来判断特征点的匹配性。

总结来说，SIFT算法通过对图像进行多尺度的检测、关键点的描述和特征的匹配，实现了对图像特征的提取和匹配。

这种算法对于图像的尺度变化、旋转变换和部分遮挡等具有一定的不变性，因此被广泛应用于图像处理、目标识别和三维重建等领域。

SIFT特征提取算法详解SIFT（Scale-Invariant Feature Transform）特征提取算法是一种用于在图像中寻找关键点的方法，该算法具有尺度不变性，能够提取出能够对图像进行描述的稳定的局部特征。

本文将从算法原理、关键步骤以及优缺点等方面进行详细解析。

一、算法原理SIFT算法的主要思想是通过寻找关键点和计算这些关键点的特征向量，来描述图像中的局部特征，并实现对尺度、旋转和亮度的不变性。

1.尺度空间极值检测SIFT算法首先通过DoG（Difference of Gaussian）来检测图像中的关键点。

DoG是指在不同尺度下的高斯平滑图像之间的差异，这样可以有效地检测出图像中的尺度空间极值点。

2.关键点精确定位通过在DoG金字塔中精确定位关键点，SIFT算法可以找到图像中的关键点。

该算法使用了Hessian矩阵的近似来计算特征点的位置和尺度，并用Taylor展开来精确计算特征点的位置。

3.方向分配为了使特征具有旋转不变性，SIFT算法通过计算特征点周围的梯度方向直方图来分配特征点的主方向。

具体步骤是确定特征点周围像素的梯度幅值和方向，并将方向分配到主方向上。

4.特征描述在计算了特征点的尺度和方向之后，SIFT算法通过计算特征点周围的局部图像块的梯度直方图来描述特征。

这些梯度直方图能够描述特征点周围的领域特征，并且具有尺度和旋转不变性。

二、关键步骤详解1.构建高斯金字塔SIFT算法首先对输入的图像进行平滑处理，然后通过不断降采样生成高斯金字塔。

高斯金字塔能够提供多尺度图像。

2.构建差分金字塔利用高斯金字塔计算DoG金字塔。

通过对相邻两层高斯图像进行差分操作，可以得到相应的两层DoG图像。

3.检测尺度空间极值点在DoG金字塔上，通过比较每一个像素点与其8邻域和相邻的26个像素点的像素值来确定极值点。

具体方法是通过比较每一层图像的像素值与周围像素值进行判定。

4.精确定位关键点通过利用Taylor展开公式来精确定位关键点的位置，这可以减小特征点在尺度和位置上的误差。

特征点匹配——SIFT算法详解SIFT（Scale-Invariant Feature Transform）是一种用于在图像中寻找关键点并进行匹配的算法。

该算法由David Lowe在1999年发布，并且一直被广泛应用于计算机视觉领域。

SIFT算法具有尺度不变性和旋转不变性，可以在不同的图像尺度和旋转角度下进行特征点的匹配。

SIFT算法的主要步骤包括关键点检测、关键点描述和特征点匹配。

关键点检测：在一张图像中，关键点通常是指存在于不同尺度和方向上的局部最大值或局部最小值。

SIFT使用高斯差分金字塔来检测关键点。

首先，通过对原始图像进行高斯模糊，创建一个金字塔，然后在每一组金字塔中计算高斯差分图像。

接着，通过比较每个像素周围的8个像素和自身像素的差值，找到局部极值点。

最后，使用尺度空间极大值抑制来进一步过滤出稳定的关键点。

关键点描述：在关键点检测后，需要对每个关键点进行描述。

SIFT使用局部图像梯度的直方图来描述关键点。

首先，在每个关键点周围的16x16像素块上计算梯度的幅值和方向。

然后将这个块分成16个4x4的子块，并在每个子块上计算一个8方向的直方图。

最后，将这些直方图连接起来形成一个128维的向量，用来表示该关键点。

在对两幅图像提取出关键点并进行描述后，需要对这些特征点进行匹配。

SIFT使用欧式距离来计算特征向量之间的相似性。

对于每个特征点，将其描述子与另一幅图像中的所有描述子进行比较，选择最佳匹配的特征点对。

SIFT算法在匹配过程中还引入了RANSAC算法来排除错误的匹配。

RANSAC（Random Sample Consensus）是一种迭代的鲁棒性估计方法，可以通过随机选择一个小子集来估计模型参数，并通过计算剩余误差和阈值来确定最终的模型。

总结一下，SIFT算法通过关键点检测、关键点描述和特征点匹配三个步骤来实现对图像中的特征点进行匹配。

该算法具有尺度不变性和旋转不变性，可以在不同的图像尺度和旋转角度下进行特征点的匹配。

SIFT算法原理SIFT算法详细介绍SIFT（Scale-Invariant Feature Transform）算法是一种用于图像特征提取和匹配的算法，被广泛应用于计算机视觉领域。

SIFT算法的核心思想是在不同尺度空间中检测和描述图像的局部特征，通过这些特征可实现图像的匹配、识别和定位。

1.尺度空间极值点检测：首先，SIFT算法通过高斯平滑滤波器构建高斯金字塔，每一层都是由上一层进行降采样得到的。

在不同尺度空间的图像上，通过比较每个像素点与其相邻像素点的灰度值，找出尺度空间极值点，作为关键点候选。

2.关键点定位：通过对尺度空间的极值点进行精确定位，以获取更加准确的关键点。

在关键点候选集合中，利用尺度空间的极值点与梯度方向直方图等局部特征进行过滤，剔除低对比度和边缘响应的关键点。

3.方向分配：为了使SIFT算法具有旋转不变性，需要为每个关键点分配主方向。

采用梯度直方图对关键点周围的图像区域进行统计，找到梯度方向直方图的峰值，作为该关键点的主方向。

4.关键点描述：在每个关键点周围的图像区域内，构建128维的特征向量，用于描述该关键点的局部特征。

通过计算每个像素点周围的梯度幅值和方向，并将其转换为特征向量的模式，构成关键点的描述符。

5.特征匹配：对于待匹配的图像，通过计算两个图像的特征向量之间的欧氏距离，进行特征匹配。

采用距离比值方法，选取最匹配和次匹配之间的距离比值小于预先设定的阈值的特征点，进行最终匹配。

6.特征筛选：通过应用RANSAC（随机抽样一致性）算法，对匹配结果进行筛选，剔除外点，保留内点，从而得到更准确的匹配结果。

SIFT算法的应用非常广泛，常用于目标识别、图像拼接、图像检索、三维重建和运动估计等领域。

在目标识别中，通过提取图像的SIFT特征，并与模板图像进行匹配，可以实现目标的识别和定位。

在图像拼接中，通过寻找多幅图像之间的共同特征点，并根据特征点的位置关系进行配准，可以实现图像的自动拼接。

SIFT算法原理SIFT算法详细介绍SIFT（Scale-Invariant Feature Transform）是一种用于计算机视觉领域的算法，用于检测和描述图像中的局部特征。

它被广泛应用于目标识别、图像匹配和三维重建等任务中。

SIFT算法通过提取图像的稳定特征点来实现图像的尺度和旋转不变性。

下面详细介绍SIFT算法的原理和步骤。

1.尺度空间极值检测：SIFT算法首先在不同的尺度空间中通过高斯差分金字塔（Difference of Gaussian，DoG）寻找稳定的特征点。

通过对输入图像进行高斯滤波，得到一系列不同尺度的图像。

然后，通过对相邻的两个不同尺度的图像进行差分操作，得到高斯差分图像，即DoG金字塔。

接着，在DoG金字塔中寻找局部极值点，即该点的像素值在其周围的3×3×3邻域内最大或最小。

2.生成关键点：在尺度空间极值点检测后，通过插值计算亚像素精度的关键点位置，以获得更精确的特征点位置。

对比邻域像素的梯度幅值和方向，重新定位关键点位置。

3.消除边缘响应：排除低对比度的稳定特征点和位于边缘的特征点，以提高匹配的准确性。

通过计算Hessian矩阵的迹和行列式来判断是否为边缘响应。

4.计算主方向：为了使SIFT算法对旋转具有不变性，对每个关键点计算该点的主方向。

在关键点周围的邻域内，计算梯度幅值和方向直方图，选取主方向作为该特征点的方向描述符。

5.生成特征描述子：在关键点检测和主方向计算后，利用关键点附近的图像区域创建描述子。

以关键点为中心，将图像区域分为若干个子区域，并在每个子区域内计算局部特征。

对每个子区域，计算梯度幅值和方向直方图，形成一个向量。

最后将这些向量串联形成一个特征向量，作为该特征点的描述子。

6.特征点匹配：使用描述子来匹配不同图像中的特征点。

通过计算两个特征点描述子之间的距离来判断它们的相似性。

通常使用欧氏距离或余弦相似度来度量特征点之间的差异。

然后，根据距离进行特征点匹配，通过选取最佳匹配对的阈值来过滤不准确的匹配。

Python语⾔实现SIFT算法⽬录⼀、什么是SIFT算法⼆、准备⼯作2.1 实验设备2.2 OpenCV安装三、实验⼯作3.1 图像选择3.2 程序实现3.3 程序结果本⽂侧重于如何使⽤Python语⾔实现SIFT算法所有程序已打包：⼀、什么是SIFT算法 SIFT，即尺度不变特征变换（Scale-invariant feature transform，SIFT），是⽤于图像处理领域的⼀种描述。

这种描述具有尺度不变性，可在图像中检测出关键点，是⼀种局部特征描述⼦。

⼆、准备⼯作2.1 实验设备 本⽂在Windows10系统上，使⽤pycharm软件完成所有实验。

2.2 OpenCV安装 我们可以使⽤OpenCV库中的cv2.xfeatures2d.SIFT_create()函数实现SIFT，但由于专利保护，很多版本的OpenCV库已⽆法提供该函数，⽬前仅3.4.2.16版本的OpenCV库可使⽤此函数。

安装教程： （1）查看当前版本opencv：进⼊cmd（组合键win+R，输⼊cmd），输⼊conda list，查看当前pycharm所有库并找到opencv-python，若找不到库，说明没有安装。

（2）卸载原版本（在cmd中输⼊：pip uninstall opencv） （3）安装新版本（在cmd中输⼊：pip install opencv-python==3.4.2.16 -i "https:///simple/" （4）安装附属库（在cmd中输⼊：pip install opencv-contrib-python==3.4.2.16 -i "https:///simple/"）三、实验⼯作3.1 图像选择 这⾥选择经典的lena图像作为实验对象，为了选择⼀个待匹配图像，本⽂使⽤如下代码对lena图像进⾏逆时针45°旋转。

from PIL import Imageimg = Image.open('lena.png')img2 = img.rotate(45) # 逆时针旋转45°img2.save("lena_rot45.png")img2.show()参考图像与待匹配图像（即旋转图像）如下图所⽰：3.2 程序实现"""图像匹配——SIFT点特征匹配实现步骤：（1）读取图像；（2）定义sift算⼦；（3）通过sift算⼦对需要匹配的图像进⾏特征点获取；a.可获取各匹配图像经过sift算⼦的特征点数⽬（4）可视化特征点（在原图中标记为圆圈）；a.为⽅便观察，可将匹配图像横向拼接（5）图像匹配（特征点匹配）；a.通过调整ratio获取需要进⾏图像匹配的特征点数量（ratio值越⼤，匹配的线条越密集，但错误匹配点也会增多）b.通过索引ratio选择固定的特征点进⾏图像匹配（6）将待匹配图像通过旋转、变换等⽅式将其与⽬标图像对齐"""import cv2 # opencv版本需为3.4.2.16import numpy as np # 矩阵运算库import time # 时间库original_lena = cv2.imread('lena.png') # 读取lena原图lena_rot45 = cv2.imread('lena_rot45.png') # 读取lena旋转45°图sift = cv2.xfeatures2d.SIFT_create()# 获取各个图像的特征点及sift特征向量# 返回值kp包含sift特征的⽅向、位置、⼤⼩等信息；des的shape为（sift_num， 128）， sift_num表⽰图像检测到的sift特征数量(kp1, des1) = sift.detectAndCompute(original_lena, None)(kp2, des2) = sift.detectAndCompute(lena_rot45, None)# 特征点数⽬显⽰print("=========================================")print("=========================================")print('lena 原图特征点数⽬：', des1.shape[0])print('lena 旋转图特征点数⽬：', des2.shape[0])print("=========================================")print("=========================================")# 举例说明kp中的参数信息for i in range(2):print("关键点", i)print("数据类型:", type(kp1[i]))print("关键点坐标:", kp1[i].pt)print("邻域直径:", kp1[i].size)print("⽅向:", kp1[i].angle)print("所在的图像⾦字塔的组:", kp1[i].octave)print("=========================================")print("=========================================")"""⾸先对原图和旋转图进⾏特征匹配，即图original_lena和图lena_rot45"""# 绘制特征点，并显⽰为红⾊圆圈sift_original_lena = cv2.drawKeypoints(original_lena, kp1, original_lena, color=(255, 0, 255))sift_lena_rot45 = cv2.drawKeypoints(lena_rot45, kp2, lena_rot45, color=(255, 0, 255))sift_cat1 = np.hstack((sift_original_lena, sift_lena_rot45)) # 对提取特征点后的图像进⾏横向拼接cv2.imwrite("sift_cat1.png", sift_cat1)print('原图与旋转图特征点绘制图像已保存')cv2.imshow("sift_point1", sift_cat1)cv2.waitKey()# 特征点匹配# K近邻算法求取在空间中距离最近的K个数据点，并将这些数据点归为⼀类start = time.time() # 计算匹配点匹配时间bf = cv2.BFMatcher()matches1 = bf.knnMatch(des1, des2, k=2)print('⽤于原图和旋转图图像匹配的所有特征点数⽬：', len(matches1))# 调整ratio# ratio=0.4：对于准确度要求⾼的匹配；# ratio=0.6：对于匹配点数⽬要求⽐较多的匹配；# ratio=0.5：⼀般情况下。