SIFT 算法实现及代码详解

safesky 贡献于2012-05-29

作者 User  创建于2012-04-17 13:38:00   修改者User  修改于2012-04-17 13:47:00字数27956

文档摘要: Sift是David Lowe于1999年提出的局部特征描述子,并于2004年进行了更深入的发展和完善。Sift特征匹配算法可以处理两幅图像之间发生平移、旋转、仿射变换情况下的匹配问题,具有很强的匹配能力。在Mikolajczyk对包括Sift算子在内的十种局部描述子所做的不变性对比实验中,Sift及其扩展算法已被证实在同类描述子中具有最强的健壮性。
关键词:

经典算法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://sourceforge.net/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/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://sourceforge.net/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 plaincopy 1. 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( calloc( n, sizeof(struct feature) ) );   51.     *feat = static_cast( 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位灰度图以及进行高斯模糊。 [cpp] view plaincopy 1. static IplImage* create_init_img( IplImage* img, int img_dbl, double sigma )   2. {   3.     IplImage* gray, * dbl;   4.     float sig_diff;   5.    6.     gray = convert_to_gray32( img );   7.     if( img_dbl )   8.     {   9.         sig_diff = sqrt( sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4 );   10.         dbl = cvCreateImage( cvSize( img->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.     else   18.     {   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是初始的高斯模糊参数,后续也通过它计算每一层所使用的sigma [cpp] view plaincopy 1. static IplImage*** build_gauss_pyr( IplImage* base, int octvs,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_pyr[i] = static_cast( calloc( intvls + 3, sizeof( IplImage* ) ) );   11.    12.     /*  13.         precompute Gaussian sigmas using the following formula:  14.         预计算每次高斯模糊的sigma  15.   16.         \sigma_{total}^2 = \sigma_{i}^2 + \sigma_{i-1}^2  17.     */   18.     sig[0] = 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.         sig[i] = 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_pyr[o][i] = cvCloneImage(base);   34.    35.             /* base of new octvave is halved image from end of previous octave */   36.             //每一组的第一层都是通过对前面一组的第一层降采样实现的   37.             else if( i == 0 )   38.                 gauss_pyr[o][i] = downsample( gauss_pyr[o-1][intvls] );   39.    40.             /* blur the current octave's last image to create the next one */   41.             //每一组的其他层则使通过使用不同sigma的高斯模糊来进行处理   42.             else   43.             {   44.                 gauss_pyr[o][i] = cvCreateImage( cvGetSize(gauss_pyr[o][i-1]),   45.                     IPL_DEPTH_32F, 1 );   46.                 cvSmooth( gauss_pyr[o][i-1], gauss_pyr[o][i],   47.                     CV_GAUSSIAN, 0, 0, sig[i], sig[i] );   48.             }   49.         }   50.    51.     free( sig );   52.     return gauss_pyr;   53. }   降采样处理 输入参数: 不解释 这就是降采样,其实就是将图像通过最近邻算法缩小为原来的一半 [cpp] view plaincopy 1. 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)构建高斯差分金字塔 输入参数: 不解释了参见上面的说明即可 实际上差分金字塔的构成是通过对相邻层的图像进行相减获得的 [cpp] view plaincopy 1. 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_pyr[i] = 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_pyr[o][i] = cvCreateImage( cvGetSize(gauss_pyr[o][i]),   14.                 IPL_DEPTH_32F, 1 );   15.             cvSub( gauss_pyr[o][i+1], gauss_pyr[o][i], dog_pyr[o][i], NULL );   16.         }   17.    18.     return dog_pyr;   19. }   (4)极值点检测 输入参数: contr_thr是去除对比度低的点所采用的阈值 curv_thr是去除边缘特征的阈值 [cpp] view plaincopy 1. 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 < dog_pyr[o][0]->height-SIFT_IMG_BORDER; r++)   15.                 for(c = SIFT_IMG_BORDER; c < dog_pyr[o][0]->width-SIFT_IMG_BORDER; c++)   16.                     /* perform preliminary check on contrast */   17.                     if( ABS( pixval32f( dog_pyr[o][i], 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_pyr[ddata->octv][ddata->intvl],   25.                                     ddata->r, ddata->c, curv_thr ) )   26.                                 {   27.                                     cvSeqPush( features, feat );   28.                                 }   29.                                 else   30.                                     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个像素中最小的一个或者是否为最大的一个,如果是则为极值点。 [cpp] view plaincopy 1. static int is_extremum( IplImage*** dog_pyr, int octv, int intvl, int r, int c )   2. {   3.     float val = pixval32f( dog_pyr[octv][intvl], 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_pyr[octv][intvl+i], r + j, c + k ) )   13.                         return 0;   14.     }   15.    16.     /* check for minimum */   17.     else   18.     {   19.         for( i = -1; i <= 1; i++ )   20.             for( j = -1; j <= 1; j++ )   21.                 for( k = -1; k <= 1; k++ )   22.                     if( val > pixval32f( dog_pyr[octv][intvl+i], r + j, c + k ) )   23.                         return 0;   24.     }   25.    26.     return 1;   27. }   *获取亚像素的极值点的位置 [cpp] view plaincopy 1. static struct 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 < 1  ||   20.             intvl > intvls  ||   21.             c < SIFT_IMG_BORDER  ||   22.             r < SIFT_IMG_BORDER  ||   23.             c >= dog_pyr[octv][0]->width - SIFT_IMG_BORDER  ||   24.             r >= dog_pyr[octv][0]->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 ) < contr_thr / intvls )   41.         return NULL;   42.     //若小于,则认为是极值点   43.     feat = new_feature();   44.     ddata = feat_detection_data( feat );   45.     feat->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. }   *获取亚像素位置中所用到的函数 [cpp] view plaincopy 1. static void 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 x[3] = { 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 = x[2];   23.     *xr = x[1];   24.     *xc = x[0];   25. }   *计算三维偏导数 计算在x和y方向上的偏导数,高斯差分尺度空间金字塔中像素的尺度 实际上在离散数据中计算偏导数是通过相邻像素的相减来计算的 比如说计算x方向的偏导数dx,则通过该向所的x方向的后一个减去前一个然后除以2即可求的dx [cpp] view plaincopy 1. 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_pyr[octv][intvl], r, c+1 ) -   7.         pixval32f( dog_pyr[octv][intvl], r, c-1 ) ) / 2.0;   8.     dy = ( pixval32f( dog_pyr[octv][intvl], r+1, c ) -   9.         pixval32f( dog_pyr[octv][intvl], r-1, c ) ) / 2.0;   10.     ds = ( pixval32f( dog_pyr[octv][intvl+1], r, c ) -   11.         pixval32f( dog_pyr[octv][intvl-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. }   *计算三维海森矩阵 不需要讲什么,其实就是计算二次导数,计算方法也和一次导数的计算如出一辙。 然后将结果放入到一个矩阵中去。 [cpp] view plaincopy 1. 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_pyr[octv][intvl], r, c );   7.     dxx = ( pixval32f( dog_pyr[octv][intvl], r, c+1 ) +    8.             pixval32f( dog_pyr[octv][intvl], r, c-1 ) - 2 * v );   9.     dyy = ( pixval32f( dog_pyr[octv][intvl], r+1, c ) +   10.             pixval32f( dog_pyr[octv][intvl], r-1, c ) - 2 * v );   11.     dss = ( pixval32f( dog_pyr[octv][intvl+1], r, c ) +   12.             pixval32f( dog_pyr[octv][intvl-1], r, c ) - 2 * v );   13.     dxy = ( pixval32f( dog_pyr[octv][intvl], r+1, c+1 ) -   14.             pixval32f( dog_pyr[octv][intvl], r+1, c-1 ) -   15.             pixval32f( dog_pyr[octv][intvl], r-1, c+1 ) +   16.             pixval32f( dog_pyr[octv][intvl], r-1, c-1 ) ) / 4.0;   17.     dxs = ( pixval32f( dog_pyr[octv][intvl+1], r, c+1 ) -   18.             pixval32f( dog_pyr[octv][intvl+1], r, c-1 ) -   19.             pixval32f( dog_pyr[octv][intvl-1], r, c+1 ) +   20.             pixval32f( dog_pyr[octv][intvl-1], r, c-1 ) ) / 4.0;   21.     dys = ( pixval32f( dog_pyr[octv][intvl+1], r+1, c ) -   22.             pixval32f( dog_pyr[octv][intvl+1], r-1, c ) -   23.             pixval32f( dog_pyr[octv][intvl-1], r+1, c ) +   24.             pixval32f( dog_pyr[octv][intvl-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. }   *计算插入像素的对比度 [cpp] view plaincopy 1. static double 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 t[1], x[3] = { 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_pyr[octv][intvl], r, c ) + t[0] * 0.5;   14. }   其中cvGEMM是矩阵的通用计算函数,至于CV_GEMM_A_T是计算dD的转置矩阵放入T中 *去除边缘相应 通过计算所在特征向量的主曲率半径来判断特征是边缘的从而导致不稳定 即去除边缘响应 [cpp] view plaincopy 1. static int 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 < ( curv_thr + 1.0 )*( curv_thr + 1.0 ) / curv_thr )   19.         return 0;   20.     return 1;   21. }   (4)计算特征向量的尺度 实际上是通过最初的sigma来获得每一层每一组的尺度 [cpp] view plaincopy 1. static void calc_feature_scales( CvSeq* features, double sigma, int intvls )   2. {   3.     struct feature* feat;   4.     struct detection_data* ddata;   5.     double intvl;   6.     int i, n;   7.    8.     n = features->total;   9.     for( i = 0; i < n; i++ )   10.     {   11.         feat = CV_GET_SEQ_ELEM( struct feature, features, i );   12.         ddata = feat_detection_data( feat );   13.         intvl = ddata->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)调整图像特征坐标、尺度、点的坐标大小为原来的一半 [java] view plaincopy 1. static void 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 < n; i++ )   8.     {   9.         feat = CV_GET_SEQ_ELEM( struct feature, features, i );   10.         feat->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)给每一个图像特征向量计算规范化的方向 [cpp] view plaincopy 1. static void 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_pyr[ddata->octv][ddata->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. }   *对所给像素计算灰度方向直方图 以关键点为中心的邻域窗口内采样,并用直方图统计邻域像素的梯度 方向。梯度直方图的范围是0~360度,其中每10度一个柱,总共36个柱 [cpp] view plaincopy 1. static double* ori_hist( IplImage* img, int r, int c, int n, int rad, double sigma)   2. {   3.     double* hist;   4.     double mag, ori, w, exp_denom, PI2 = CV_PI * 2.0;   5.     int bin, i, j;   6.    7.     hist = static_cast( calloc( n, sizeof( double ) ) );   8.     exp_denom = 2.0 * sigma * sigma;   9.     for( i = -rad; i <= rad; i++ )   10.         for( j = -rad; j <= rad; j++ )   11.             if( calc_grad_mag_ori( img, r + i, c + j, &mag, &ori ) )   12.             {   13.                 w = exp( -( i*i + j*j ) / exp_denom );   14.                 bin = cvRound( n * ( ori + CV_PI ) / PI2 );   15.                 bin = ( bin < n )? bin : 0;   16.                 hist[bin] += w * mag;   17.             }   18.    19.     return hist;   20. }   *计算所给像素的梯度大小和方向 每一个小格都代表了特征点邻域所在的尺度空间的一个像素 ,箭头方向代表了像素梯 度方向,箭头长度代表该像素的幅值也就是梯度的值 [cpp] view plaincopy 1. static int calc_grad_mag_ori( IplImage* img, int r, int c, double* mag, double* ori )   2. {   3.     double dx, dy;   4.    5.     if( r > 0  &&  r < img->height - 1  &&  c > 0  &&  c < img->width - 1 )   6.     {   7.         dx = pixval32f( img, r, c+1 ) - pixval32f( img, r, c-1 );   8.         dy = pixval32f( img, r-1, c ) - pixval32f( img, r+1, c );   9.         *mag = sqrt( dx*dx + dy*dy );   10.         *ori = atan2( dy, dx );   11.         return 1;   12.     }   13.    14.     else   15.         return 0;   16. }   *对方向直方图进行高斯模糊 使用高斯函数对直方图进行平滑,减少突变的影响。 [cpp] view plaincopy 1. static void smooth_ori_hist( double* hist, int n )   2. {   3.     double prev, tmp, h0 = hist[0];   4.     int i;   5.    6.     prev = hist[n-1];   7.     for( i = 0; i < n; i++ )   8.     {   9.         tmp = hist[i];   10.         hist[i] = 0.25 * prev + 0.5 * hist[i] +    11.             0.25 * ( ( i+1 == n )? h0 : hist[i+1] );   12.         prev = tmp;   13.     }   14. }   *在直方图中找到主方向的梯度 利用关键点邻域像素的梯度方向分布特性为每个关键点指定方向参数,使算子具备 旋转不变性。 [cpp] view plaincopy 1. static double dominant_ori( double* hist, int n )   2. {   3.     double omax;   4.     int maxbin, i;   5.    6.     omax = hist[0];   7.     maxbin = 0;   8.     for( i = 1; i < n; i++ )   9.         if( hist[i] > omax )   10.         {   11.             omax = hist[i];   12.             maxbin = i;   13.         }   14.     return omax;   15. }   *将大于某一个梯度大小阈值的特征向量加入到直方图中去 n为方向的个数 [cpp] view plaincopy 1. mag_thr描述子向量门限一般取0.2   [cpp] view plaincopy 1. static void add_good_ori_features( CvSeq* features, double* hist, int n,   2.                                    double mag_thr, struct feature* feat )   3. {   4.     struct feature* new_feat;   5.     double bin, PI2 = CV_PI * 2.0;   6.     int l, r, i;   7.    8.     for( i = 0; i < n; i++ )   9.     {   10.         l = ( i == 0 )? n - 1 : i-1;   11.         r = ( i + 1 ) % n;   12.            13.    14.         //描述子向量门限化,一般门限取0.2   15.         if( hist[i] > hist[l]  &&  hist[i] > hist[r]  &&  hist[i] >= mag_thr )   16.         {   17.             bin = i + interp_hist_peak( hist[l], hist[i], hist[r] );   18.             bin = ( bin < 0 )? n + bin : ( bin >= n )? bin - n : bin;   19.             new_feat = clone_feature( feat );   20.             new_feat->ori = ( ( PI2 * bin ) / n ) - CV_PI;   21.             cvSeqPush( features, new_feat );   22.             free( new_feat );   23.         }   24.     }   25. }   (7)计算特征描述子 [cpp] view plaincopy 1. static void compute_descriptors( CvSeq* features, IplImage*** gauss_pyr, int d, int n)   2. {   3.     struct feature* feat;   4.     struct detection_data* ddata;   5.     double*** hist;   6.     int i, k = features->total;   7.    8.     for( i = 0; i < k; i++ )   9.     {   10.         feat = CV_GET_SEQ_ELEM( struct feature, features, i );   11.         ddata = feat_detection_data( feat );   12.         //计算二维方向直方图   13.         hist = descr_hist( gauss_pyr[ddata->octv][ddata->intvl], ddata->r,   14.             ddata->c, feat->ori, ddata->scl_octv, d, n );   15.         //将二维方向直方图转换为特征描述子   16.         hist_to_descr( hist, d, n, feat );   17.         release_descr_hist( &hist, d );   18.     }   19. }   *计算二维方向直方图 [cpp] view plaincopy 1. static double*** descr_hist( IplImage* img, int r, int c, double ori,   2.                              double scl, int d, int n )   3. {   4.     double*** hist;   5.     double cos_t, sin_t, hist_width, exp_denom, r_rot, c_rot, grad_mag,   6.         grad_ori, w, rbin, cbin, obin, bins_per_rad, PI2 = 2.0 * CV_PI;   7.     int radius, i, j;   8.    9.     hist = static_cast( calloc( d, sizeof( double** ) ) );   10.     for( i = 0; i < d; i++ )   11.     {   12.         hist[i] =static_cast( calloc( d, sizeof( double* ) ) );   13.         for( j = 0; j < d; j++ )   14.             hist[i][j] = static_cast( calloc( n, sizeof( double ) ) );   15.     }   16.    17.     cos_t = cos( ori );   18.     sin_t = sin( ori );   19.     bins_per_rad = n / PI2;   20.     exp_denom = d * d * 0.5;   21.     hist_width = SIFT_DESCR_SCL_FCTR * scl;   22.     radius = hist_width * sqrt(2.0) * ( d + 1.0 ) * 0.5 + 0.5;   23.     for( i = -radius; i <= radius; i++ )   24.         for( j = -radius; j <= radius; j++ )   25.         {   26.             /*  27.             即将坐标移至关键点主方向  28.             计算采用的直方图数组中相对于方向旋转的坐标  29.             Calculate sample's histogram array coords rotated relative to ori.  30.             Subtract 0.5 so samples that fall e.g. in the center of row 1 (i.e.  31.             r_rot = 1.5) have full weight placed in row 1 after interpolation.  32.             */   33.             c_rot = ( j * cos_t - i * sin_t ) / hist_width;   34.             r_rot = ( j * sin_t + i * cos_t ) / hist_width;   35.             rbin = r_rot + d / 2 - 0.5;   36.             cbin = c_rot + d / 2 - 0.5;   37.    38.             if( rbin > -1.0  &&  rbin < d  &&  cbin > -1.0  &&  cbin < d )   39.                 if( calc_grad_mag_ori( img, r + i, c + j, &grad_mag, &grad_ori ))   40.                 {   41.                     grad_ori -= ori;   42.                     while( grad_ori < 0.0 )   43.                         grad_ori += PI2;   44.                     while( grad_ori >= PI2 )   45.                         grad_ori -= PI2;   46.    47.                     obin = grad_ori * bins_per_rad;   48.                     w = exp( -(c_rot * c_rot + r_rot * r_rot) / exp_denom );   49.                     interp_hist_entry( hist, rbin, cbin, obin, grad_mag * w, d, n );   50.                 }   51.         }   52.    53.     return hist;   54. }   *插入一个entry进入到方向直方图中从而形成特征描述子 这个,我也不怎么明白。。。 [cpp] view plaincopy 1. static void interp_hist_entry( double*** hist, double rbin, double cbin,   2.                                double obin, double mag, int d, int n )   3. {   4.     double d_r, d_c, d_o, v_r, v_c, v_o;   5.     double** row, * h;   6.     int r0, c0, o0, rb, cb, ob, r, c, o;   7.    8.     r0 = cvFloor( rbin );   9.     c0 = cvFloor( cbin );   10.     o0 = cvFloor( obin );   11.     d_r = rbin - r0;   12.     d_c = cbin - c0;   13.     d_o = obin - o0;   14.    15.     /*  16.     The entry is distributed into up to 8 bins.  Each entry into a bin  17.     is multiplied by a weight of 1 - d for each dimension, where d is the  18.     distance from the center value of the bin measured in bin units.  19.     */   20.     for( r = 0; r <= 1; r++ )   21.     {   22.         rb = r0 + r;   23.         if( rb >= 0  &&  rb < d )   24.         {   25.             v_r = mag * ( ( r == 0 )? 1.0 - d_r : d_r );   26.             row = hist[rb];   27.             for( c = 0; c <= 1; c++ )   28.             {   29.                 cb = c0 + c;   30.                 if( cb >= 0  &&  cb < d )   31.                 {   32.                     v_c = v_r * ( ( c == 0 )? 1.0 - d_c : d_c );   33.                     h = row[cb];   34.                     for( o = 0; o <= 1; o++ )   35.                     {   36.                         ob = ( o0 + o ) % n;   37.                         v_o = v_c * ( ( o == 0 )? 1.0 - d_o : d_o );   38.                         h[ob] += v_o;   39.                     }   40.                 }   41.             }   42.         }   43.     }   44. }   *将二维直方图转换为特征描述子 实际上是归一化描述子和转换为整数 [cpp] view plaincopy 1. static void hist_to_descr( double*** hist, int d, int n, struct feature* feat )   2. {   3.     int int_val, i, r, c, o, k = 0;   4.    5.     for( r = 0; r < d; r++ )   6.         for( c = 0; c < d; c++ )   7.             for( o = 0; o < n; o++ )   8.                 feat->descr[k++] = hist[r][c][o];   9.    10.     feat->d = k;   11.     normalize_descr( feat );   12.     for( i = 0; i < k; i++ )   13.         if( feat->descr[i] > SIFT_DESCR_MAG_THR )   14.             feat->descr[i] = SIFT_DESCR_MAG_THR;   15.     normalize_descr( feat );   16.    17.     /* convert floating-point descriptor to integer valued descriptor */   18.     for( i = 0; i < k; i++ )   19.     {   20.         int_val = SIFT_INT_DESCR_FCTR * feat->descr[i];   21.         feat->descr[i] = MIN( 255, int_val );   22.     }   23. }   *归一化描述子 [cpp] view plaincopy 1. static void normalize_descr( struct feature* feat )   2. {   3.     double cur, len_inv, len_sq = 0.0;   4.     int i, d = feat->d;//为描述子长度128维   5.        6.    7.     //如何进行归一化特征描述子来降低对光照的影响   8.     //主要就是将每一个特征的平方求和,然后开方,然后去其倒数,然后乘以每一个特征描述子的梯度值,这样就得到了归一化的特征描述子   9.     for( i = 0; i < d; i++ )   10.     {   11.         cur = feat->descr[i];   12.         len_sq += cur*cur;   13.     }   14.     len_inv = 1.0 / sqrt( len_sq );   15.     for( i = 0; i < d; i++ )   16.         feat->descr[i] *= len_inv;   17. }   下面给出存储在文件中的SIFT特征分析: 下面给出特征文件的部分特征 [plain] view plaincopy 1. 114  128     2. 101.350424   136.130888   40.169873   0.771085  orientation   3.  0 0 0 0 3 1 0 0 2 23 46 15 18 3 0 0 6 20 13 1   4.  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 88 36 0 0   5.  81 95 57 47 185 114 2 7 185 155 19 6 19 6 1 22 22 0 0 0   6.  0 0 0 1 0 0 0 0 37 8 0 0 91 12 0 1 185 144 11 35   7.  185 50 0 0 23 28 8 95 40 1 0 0 0 0 0 4 0 0 0 0   8.  0 0 0 0 11 5 0 0 4 2 0 0 49 20 0 0 1 0 0 1   9.  0 0 0 0 0 0 0 0   10. 127.871534 71.100559 15.768594 -2.024589   11.  1 2 2 72 63 12 1 1 133 93 1 4 2 7 4 44 133 115 0 0   12.  0 0 0 20 9 4 0 0 0 0 0 0 23 0 1 9 107 20 1 8   13.  133 5 0 0 0 1 5 133 132 14 0 0 0 0 8 133 14 1 0 0   14.  0 0 0 8 26 0 0 0 126 37 8 22 133 47 0 0 0 0 3 52   15.  131 41 0 0 0 0 2 36 1 0 0 0 0 0 0 2 2 0 0 0   16.  34 105 80 24 111 15 0 0 0 1 55 66 79 21 0 0 0 0 0 5   17.  0 0 0 0 0 0 0 0   下面给出说明: [plain] view plaincopy 1. 114 特征数目 128  向量维度   2.    3. 关键点坐标   4. 101.350424  y坐标   5. 136.130888  x坐标   6. 40.169873  scale 尺度   7. 0.771085  orientation  关键点的梯度方向   8.    9.    10. 16个种子点的8个方向向量的信息共128个信息   11.    12.  0 0 0 0 3 1 0 0 2 23 46 15 18 3 0 0 6 20 13 1   13.  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 88 36 0 0   14.  81 95 57 47 185 114 2 7 185 155 19 6 19 6 1 22 22 0 0 0   15.  0 0 0 1 0 0 0 0 37 8 0 0 91 12 0 1 185 144 11 35   16.  185 50 0 0 23 28 8 95 40 1 0 0 0 0 0 4 0 0 0 0   17.  0 0 0 0 11 5 0 0 4 2 0 0 49 20 0 0 1 0 0 1   18.  0 0 0 0 0 0 0 0   19.    20. 下一组关键点向量   21. 127.871534 y坐标   22. 71.100559 x坐标   23. 15.768594 尺度=sigma*2^(高斯模糊)   24. -2.024589 梯度方向m(x,y)。。。。等等   最后附上一个Rob Hess的源码下载地址吧 http://blogs.oregonstate.edu/hess/code/sift/ k-d tree的优化查找算法BBF   BBF(Best Bin First)是一种改进的k-d树最近邻查询算法。从前两篇标准的k-d树查询过程可以看出其搜索过程中的“回溯”是由“查询路径”来决定的,并没有考虑查询路径上数据点本身的一些性质。BBF的查询思路就是将“查询路径”上的节点进行排序,如按各自分割超平面(称为Bin)与查询点的距离排序。回溯检查总是从优先级最高的(Best Bin)的树节点开始。另外BBF还设置了一个运行超时限制,当优先级队列中的所有节点都经过检查或者超出时间限制时,算法返回当前找到的最好结果作为近似的最近邻。采用了best-bin-first search方法就可以将k-d树扩展到高维数据集上。   下面我们通过大牛Rob Hess基于OpenCV的SIFT实现中的相关代码来具体学习下BBF算法。 /* Finds an image feature's approximate k nearest neighbors in a kd tree using Best Bin First search. @param kd_root root of an image feature kd tree @param feat image feature for whose neighbors to search @param k number of neighbors to find @param nbrs pointer to an array in which to store pointers to neighbors in order of increasing descriptor distance @param max_nn_chks search is cut off after examining this many tree entries @return Returns the number of neighbors found and stored in nbrs, or -1 on error. */ //参数和返回值参看以上注释 //基于k-d tree + bbf的k近邻查找函数 int kdtree_bbf_knn( struct kd_node* kd_root, struct feature* feat, int k, struct feature*** nbrs, int max_nn_chks ) { struct kd_node* expl;      //expl是特征k-d tree中的一个节点 struct min_pq* min_pq; //min_pq是优先级队列 struct feature* tree_feat, ** _nbrs;//tree_feat是一个SIFT特征,_nbrs中存放着查找出来的近邻特征节点 struct bbf_data* bbf_data;   //bbf_data是一个用来存放临时特征数据和特征间距离的缓存结构 int i, t = 0, n = 0;       //t是运行时限,n是查找出来的近邻个数 if( ! nbrs || ! feat || ! kd_root ) { fprintf( stderr, "Warning: NULL pointer error, %s, line %d\n", __FILE__, __LINE__ ); return -1; } _nbrs = calloc( k, sizeof( struct feature* ) );  //给查找结果分配相应大小的内存 min_pq = minpq_init();                 //min_pq队列初始化 minpq_insert( min_pq, kd_root, 0 );         //将根节点先插入到min_pq优先级队列中 while( min_pq->n > 0 && t < max_nn_chks ) //min_pq队列没有回溯完且未达到时限 { expl = ( struct kd_node*)minpq_extract_min( min_pq );//从min_pq中提取优先级最高的节点(并移除) if( ! expl ) { fprintf( stderr, "Warning: PQ unexpectedly empty, %s line %d\n", __FILE__, __LINE__ ); goto fail; } expl = explore_to_leaf( expl, feat, min_pq );    //从expl节点开始查找到叶子节点(下详)  if( ! expl ) { fprintf( stderr, "Warning: PQ unexpectedly empty, %s line %d\n", __FILE__, __LINE__ ); goto fail; } for( i = 0; i < expl->n; i++ )  //开始比较查找最近邻 { tree_feat = &expl->features[i]; bbf_data = malloc( sizeof( struct bbf_data ) ); if( ! bbf_data ) { fprintf( stderr, "Warning: unable to allocate memory," " %s line %d\n", __FILE__, __LINE__ ); goto fail; } bbf_data->old_data = tree_feat->feature_data; bbf_data->d = descr_dist_sq(feat, tree_feat);  //计算叶子节点特征和目标特征的距离 tree_feat->feature_data = bbf_data; n += insert_into_nbr_array( tree_feat, _nbrs, n, k );//判断并插入符合条件的近邻到_nbrs中 } t++; }    //释放内存并返回结果 minpq_release( &min_pq ); for( i = 0; i < n; i++ ) { bbf_data = _nbrs[i]->feature_data; _nbrs[i]->feature_data = bbf_data->old_data; free( bbf_data ); } *nbrs = _nbrs; return n; fail: minpq_release( &min_pq ); for( i = 0; i < n; i++ ) { bbf_data = _nbrs[i]->feature_data; _nbrs[i]->feature_data = bbf_data->old_data; free( bbf_data ); } free( _nbrs ); *nbrs = NULL; return -1; }    整个kdtree_bbf_knn函数包括了优先级队列的建立和k邻近查找两个过程。其中最关键的两个数据结构就是min_pq优先级队列和_nbrs存放k邻近结果的队列。min_pq优先级队列是按照各节点的分割超平面和目标查询特征点之间的距离升序排列的,第一个节点就是最小距离(优先级最高)的节点。另外_nbrs中也是按照与目标特征的距离升序排列,直接取结果的前k个特征就是对应的k近邻。注意:上述代码中的一些数据结构的定义以及一些对应的函数如:minpq_insert,minpq_extract_min, insert_into_nbr_array, descr_dist_sq等在这里就不贴了,详细代码可参看Rob Hess主页(http://blogs.oregonstate.edu/hess/)中代码参考文档。   下面来详细看看函数explore_to_leaf是如何实现的。 /* Explores a kd tree from a given node to a leaf. Branching decisions are made at each node based on the descriptor of a given feature. Each node examined but not explored is put into a priority queue to be explored later, keyed based on the distance from its partition key value to the given feature's desctiptor. @param kd_node root of the subtree to be explored @param feat feature upon which branching decisions are based @param min_pq a minimizing priority queue into which tree nodes are placed as described above @return Returns a pointer to the leaf node at which exploration ends or NULL on error. */ //参数和返回值参看以上注释 //搜索路径和优先级队列的生成函数 static struct kd_node* explore_to_leaf( struct kd_node* kd_node, struct feature* feat, struct min_pq* min_pq ) { struct kd_node* unexpl, * expl = kd_node;  //unexpl中存放着优先级队列的候选特征点                              //expl为开始搜索节点 double kv;                     //kv是分割维度的数据 int ki;                     //ki是分割维度序号 while( expl && ! expl->leaf ) { ki = expl->ki;                //获得分割节点的ki,kv数据 kv = expl->kv; if( ki >= feat->d ) { fprintf( stderr, "Warning: comparing imcompatible descriptors, %s" \ " line %d\n", __FILE__, __LINE__ ); return NULL; } if( feat->descr[ki] <= kv )        //目标特征和分割节点分割维上的数据比较 { unexpl = expl->kd_right;       //小于右子树根节点成为候选节点  expl = expl->kd_left; //并进入左子树搜索 } else { unexpl = expl->kd_left;        //大于左子树根节点成为候选节点 expl = expl->kd_right;        //并进入右子树搜索 }      //将候选节点unexpl根据目标与分割超平面的距离插入到优先级队列中 if( minpq_insert( min_pq, unexpl, ABS( kv - feat->descr[ki] ) ) )   { fprintf( stderr, "Warning: unable to insert into PQ, %s, line %d\n", __FILE__, __LINE__ ); return NULL; } } return expl;  //返回搜索路径中最后的叶子节点 }   从explore_to_leaf函数的实现中可以看到,优先级队列和搜索路径是同时生成的,这也是BBF算法的精髓所在:在二叉搜索的时候将搜索路径另一侧的分支加入到优先级队列中,供回溯时查找。而优先级队列的排序就是根据目标特征与分割超平面的距离ABS( kv - feat->descr[ki] ) 注意:是目标特征和分割超平面间的距离,不是候选节点和分割超平面的距离。如还是上两篇例子中的数据,查找(2,4.5)的k近邻,当搜索到(5,4)节点时,应将(4,7)节点加入搜索路径而将(2,3)节点选为优先级队列的候选节点,优先级的计算是:abs(4 - 4.5) = 0.5。

下载文档到电脑,查找使用更方便

文档的实际排版效果,会与网站的显示效果略有不同!!

需要 10 金币 [ 分享文档获得金币 ] 3 人已下载

下载文档