Rob Hess sift源码详解（一）

网上已经有许多关于sift算法详细解释。本文适合于已经对sift算法有一定了解的读者，写本文的目的在于对Rob Hess的sift算法的源代码进行注释，从代码的角度对提升对sift算法的理解。闲话少说，进入正文。

首先从网址http://robwhess.github.io/opensift/上下载sift源代码。解压之后不能直接运行，因为作者是在linux下编译的，所以需要自己在windows下调试，具体过程网上搜一下都有，本文在此不赘述了。

查看文件中的match.c文件包含了sift算法匹配的整个过程。

extern int sift_features( IplImage* img, struct feature** feat );

这个函数包含整个sift算法的流程，函数返回值为特征点的数目，如果没有找到特征点则返回-1。img定义为输入图像，feat定义为指向结构体feature的指针。

int sift_features( IplImage* img, struct feature** feat ){return _sift_features( img,                       feat,                    SIFT_INTVLS,                    SIFT_SIGMA,                    SIFT_CONTR_THR,                    SIFT_CURV_THR,                    SIFT_IMG_DBL,                    SIFT_DESCR_WIDTH,                    SIFT_DESCR_HIST_BINS );}int _sift_features( IplImage* img, struct feature** feat, intintvls,        doublesigma, doublecontr_thr, intcurv_thr,        intimg_dbl,intdescr_width, intdescr_hist_bins){IplImage* init_img;IplImage***gauss_pyr, *** dog_pyr;CvMemStorage*storage;CvSeq* features;intoctvs, i, n = 0; /* checkarguments *///判断是否正确读入图像if( ! img )fatal_error("NULL pointer error, %s, line %d",  __FILE__, __LINE__ );if( ! feat )fatal_error("NULL pointer error, %s, line %d",  __FILE__, __LINE__ ); /* buildscale space pyramid; smallest dimension of top level is ~4 pixels */init_img = create_init_img(img, img_dbl, sigma );//初始化图像octvs = log( MIN(init_img->width, init_img->height ) ) / log(2) - 2;//定义金字塔的组数，n=log2(min(M,N))-t，由于塔顶图像最小维数为4 pixels所以t=log2(4)=2gauss_pyr =build_gauss_pyr( init_img, octvs, intvls, sigma );//建立高斯金字塔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 );//关键点定位calc_feature_scales(features, sigma, intvls );//计算关键点尺度if(img_dbl )//adjust_for_img_dbl(features );//由于初始图像尺度加倍，最终缩小图像尺度calc_feature_oris( features,gauss_pyr );//关键点方向分配compute_descriptors(features, gauss_pyr, descr_width, descr_hist_bins );//定义特征点描述符//释放资源/* sortfeatures by decreasing scale and move from CvSeq to array */cvSeqSort( features,(CvCmpFunc)feature_cmp, NULL );//将特征向量按序排列n = features->total;*feat = calloc( n, sizeof(struct feature) );*feat = cvCvtSeqToArray( features, *feat,CV_WHOLE_SEQ );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;}//初始图像尺度扩大staticIplImage* create_init_img( IplImage* img, intimg_dbl, doublesigma ){IplImage* gray, *dbl;doublesig_diff; gray = convert_to_gray32( img );//将输入图像转化为32比特位灰度值if(img_dbl )//如果输入图像作为第0组第0层图像，则丢弃最高空域的采样率，则原始图像在构建高斯金字塔之前尺寸加倍生成-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 );//采用立方插值法扩大一倍       cvSmooth(dbl, dbl,CV_GAUSSIAN, 0, 0, sig_diff, sig_diff );//高斯模糊       cvReleaseImage(&gray);       return dbl;}else//{       sig_diff = sqrt(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA );       cvSmooth( gray,gray, CV_GAUSSIAN, 0, 0, sig_diff, sig_diff );       return gray;}} staticIplImage* convert_to_gray32( IplImage* img ){  IplImage* gray8, *gray32;   gray32 = cvCreateImage(cvGetSize(img),IPL_DEPTH_32F, 1 );  if( img->nChannels == 1 )//输入图像为1通道      gray8 = cvClone(img );  else//输入图像为3通道    {      gray8 = cvCreateImage(cvGetSize(img),IPL_DEPTH_8U, 1 );      cvCvtColor(img,gray8, CV_BGR2GRAY );//转换为灰度图    }      cvConvertScale(gray8, gray32, 1.0 / 255.0, 0 );//gray32=1.0/255*gray8，gray32的范围为[0，1]      cvReleaseImage(&gray8);      return gray32;}static IplImage*** build_gauss_pyr(IplImage* base, int octvs,int intvls, double sigma ){  IplImage*** gauss_pyr;  const int _intvls = intvls;  double sig[_intvls+3], sig_total, sig_prev, k;  inti, o;   gauss_pyr = calloc( octvs, sizeof( IplImage** ) );  for( i = 0; i < octvs; i++ )     gauss_pyr[i] = calloc( intvls + 3, sizeof( IplImage *) );//若每组中检测intvls个极值点，则高斯差分（DOG）金字塔需intvls+2层图像，而DOG金字塔由高斯金字塔相邻两层相减得到，则高斯金字塔需intvls+3层图像，本源码中intvls=3   /*   precompute Gaussian sigmas using the following formula:   \sigma_{total}^2 = \sigma_{i}^2 + \sigma_{i-1}^2  */  sig[0] = sigma;  k =pow( 2.0, 1.0 / intvls );//k=2^(1/intvls)  for( i = 1; i < intvls + 3; i++ )    {      sig_prev = pow( k, i - 1 ) * sigma;      sig_total = sig_prev * k;      sig[i] = sqrt( sig_total * sig_total - sig_prev * sig_prev );//sigma迭代公式    }   for( o = 0; o < octvs; o++ )      for( i = 0; i < intvls + 3; i++ )      {         if(o == 0  &&  i == 0 )            gauss_pyr[o][i] = cvCloneImage(base);          /*base of new octvave is halved image from end of previous octave */         elseif( i == 0 )           gauss_pyr[o][i] = downsample(gauss_pyr[o-1][intvls] );//高斯金字塔上一组图像的初始图像由前一组图像的倒数第三张图像隔点采样得到。                   /*blur the current octave's last image to create the next one */         else           {             gauss_pyr[o][i] = cvCreateImage(cvGetSize(gauss_pyr[o][i-1]),IPL_DEPTH_32F, 1 );             cvSmooth( gauss_pyr[o][i-1],gauss_pyr[o][i],CV_GAUSSIAN, 0, 0, sig[i], sig[i] );//           }     }      return gauss_pyr;} static IplImage* downsample( IplImage* img)//图像下采样{ IplImage* smaller = cvCreateImage( cvSize(img->width / 2,img->height / 2),img->depth, img->nChannels );//图像尺寸变为原来的1/2 cvResize( img, smaller, CV_INTER_NN );//采用最近邻插值  return smaller;}//建立DOG金字塔static IplImage*** build_dog_pyr(IplImage*** gauss_pyr, int octvs, int intvls ){  IplImage*** dog_pyr;  inti, o;   dog_pyr = calloc( octvs, sizeof( IplImage** ) );  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_pyr[o][i]= cvCreateImage( cvGetSize(gauss_pyr[o][i]),IPL_DEPTH_32F, 1 );         cvSub(gauss_pyr[o][i+1], gauss_pyr[o][i], dog_pyr[o][i], NULL );//高斯金字塔图像相减     }  return dog_pyr;}

未完待续