当前位置：文档之家› SIFT算法实现C语言

SIFT算法实现C语言

经典算法SIFT实现即代码解释：

以下便是sift源码库编译后的效果图：

为了给有兴趣实现sift算法的朋友提供个参考，特整理此文如下。要了解什么是sift算法，请参考：九、图像特征提取与匹配之SIFT算法。ok，咱们下面，就来利用Rob Hess维护的sift 库来实现sift算法：

首先，请下载Rob Hess维护的sift 库：

https://www.doczj.com/doc/0314012183.html,/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

https://www.doczj.com/doc/0314012183.html,/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也是一个库，也需要下载：

https://www.doczj.com/doc/0314012183.html,/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;

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;

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 );

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;

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*) )

);

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->in

tvl],

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;

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;

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 };

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;

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;

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 };

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;

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;

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;

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;

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.}

*对所给像素计算灰度方向直方图

RC4加密算法C语言实现

RC4 加密算法 C 语言实现代码文件名 RC4.cpp Encrypt.h (代码详见后文) 备注：将以上两个文件放在相同的路径(建议不要放在中文路径下)编译执行！编译环境 Microsoft Visual C++ 6.0 C-Free 5.0 代码解释 RC4 加密算法是大名鼎鼎的RSA 三人组中的头号人物Ron Rivest 在1987 年设计的密钥长度可变的流加密算法簇。之所以称其为簇，是由于其核心部分的S-box 长度可为任意，但一般为256字节。该算法的速度可以达到DES 加密的10倍左右。 RC4 算法的原理很简单，包括初始化算法和伪随机子密码生成算法两大部分。假设S-box 长度和密钥长度均为为n。先来看看算法的初始化部分(用类C伪代码表示)： for (i=0; i

} 得到的子密码sub_k用以和明文进行xor运算，得到密文，解密过程也完全相同。 RC4加密算法在C++中的实现： RC4函数(加密/解密)：其实RC4只有加密，将密文再加密一次，就是解密了。 GetKey函数：随机字符串产生器。 ByteToHex函数：把字节码转为十六进制码，一个字节两个十六进制。十六进制字符串非常适合在HTTP中传输。 HexToByte函数：把十六进制字符串，转为字节码。。 Encrypt函数：把字符串经RC4加密后，再把密文转为十六进制字符串返回，可直接用于传输。 Decrypt函数：直接密码十六进制字符串密文，再解密，返回字符串明文。源代码以下为Encrypt.h文件代码 #ifndef _ENCRYPT_RC4_ #defi ne _ENCRYPT_RC4_ #in clude #defi ne BOX_LEN 256 int GetKey(c onst un sig ned char* pass, int pass_le n, un sig ned char *out); int RC4(c onst un sig ned char* data, int data_le n, const un sig ned char* key, int key_le n, un sig ned char* out, i nt* out_le n); static void swap_byte( un sig ned char* a, un sig ned char* b); char* En crypt(co nst char* szSource, const char* szPassWord); // 加密，返回加密结果 char* Decrypt(co nst char* szSource, con st char* szPassWord); // 解密，返回解密结果 char* ByteToHex(c onst un sig ned char* vByte, const int vLe n); // 把字节码pbBuffer转为十六进制字符串，方便传输 unsigned char* HexToByte(const char* szHex); // 把十六进制字符串转为字节码pbBuffer，解码 #e ndif // #ifndef _ENCRYPT_RC4_

SIFT算法原理

3.1.1尺度空间极值检测尺度空间理论最早出现于计算机视觉领域，当时其目的是模拟图像数据的多尺度特征。随后Koendetink 利用扩散方程来描述尺度空间滤波过程，并由此证明高斯核是实现尺度变换的唯一变换核。Lindeberg ，Babaud 等人通过不同的推导进一步证明高斯核是唯一的线性核。因此，尺度空间理论的主要思想是利用高斯核对原始图像进行尺度变换，获得图像多尺度下的尺度空间表示序列，对这些序列进行尺度空间特征提取。二维高斯函数定义如下： 222()/221 (,,)2x y G x y e σσπσ-+= （5）一幅二维图像，在不同尺度下的尺度空间表示可由图像与高斯核卷积得到： (,,(,,)*(,)L x y G x y I x y σσ)= （6）其中(x,y )为图像点的像素坐标，I(x,y )为图像数据, L 代表了图像的尺度空间。σ称为尺度空间因子，它也是高斯正态分布的方差，其反映了图像被平滑的程度，其值越小表征图像被平滑程度越小，相应尺度越小。大尺度对应于图像的概貌特征，小尺度对应于图像的细节特征。因此，选择合适的尺度因子平滑是建立尺度空间的关键。在这一步里面，主要是建立高斯金字塔和DOG(Difference of Gaussian)金字塔，然后在DOG 金字塔里面进行极值检测，以初步确定特征点的位置和所在尺度。 (1)建立高斯金字塔为了得到在不同尺度空间下的稳定特征点，将图像(,)I x y 与不同尺度因子下的高斯核(,,)G x y σ进行卷积操作，构成高斯金字塔。高斯金字塔有o 阶，一般选择4阶，每一阶有s 层尺度图像，s 一般选择5层。在高斯金字塔的构成中要注意，第1阶的第l 层是放大2倍的原始图像，其目的是为了得到更多的特征点；在同一阶中相邻两层的尺度因子比例系数是k ，则第1阶第2层的尺度因子是k σ，然后其它层以此类推则可；第2阶的第l 层由第一阶的中间层尺度图像进行子抽样获得，其尺度因子是2k σ，然后第2阶的第2层的尺度因子是第1层的k 倍即3 k σ。第3阶的第1层由第2阶的中间层尺度图像进行子抽样获得。其它阶的构成以此类推。 (2)建立DOG 金字塔 DOG 即相邻两尺度空间函数之差，用(,,)D x y σ来表示，如公式(3)所示： (,,)((,,)(,,))*(,)(,,)(,,)D x y G x y k G x y I x y L x y k L x y σσσσσ=-=- (7) DOG 金字塔通过高斯金字塔中相邻尺度空间函数相减即可，如图1所示。在图中，DOG 金字塔的第l 层的尺度因子与高斯金字塔的第l 层是一致的，其它阶也一样。

SIFT算法实现及代码详解

经典算法SIFT实现即代码解释：以下便是sift源码库编译后的效果图：

为了给有兴趣实现sift算法的朋友提供个参考，特整理此文如下。要了解什么是sift算法，请参考：九、图像特征提取与匹配之SIFT算法。ok，咱们下面，就来利用Rob Hess维护的sift 库来实现sift算法：首先，请下载Rob Hess维护的sift 库： https://www.doczj.com/doc/0314012183.html,/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 https://www.doczj.com/doc/0314012183.html,/projects/opencvlibrary/files/opencv-win/1.0/OpenCV_1.0.exe

SIFT算法英文详解

SIFT: Scale Invariant Feature Transform The algorithm SIFT is quite an involved algorithm. It has a lot going on and can be come confusing, So I’ve split up the entire algorithm into multiple parts. Here’s an outline of what happens in SIFT. Constructing a scale space This is the initial preparation. You create internal representations of the original image to ensure scale invariance. This is done by generating a “scale space”. LoG Approximation The Laplacian of Gaussian is great for finding interesting points (or key points) in an image. But it’s computationally expensive. So we cheat and approximate it using the representation created earlier. Finding keypoints With the super fast approximation, we now try to find key points. These are maxima and minima in the Difference of Gaussian image we calculate in step 2 Get rid of bad key points Edges and low contrast regions are bad keypoints. Eliminating these makes the algorithm efficient and robust. A technique similar to the Harris Corner Detector is used here. Assigning an orientation to the keypoints An orientation is calculated for each key point. Any further calculations are done relative to this orientation. This effectively cancels out the effect of orientation, making it rotation invariant. Generate SIFT features Finally, with scale and rotation invariance in place, one more representation is generated. This helps uniquely identify features. Lets say you have 50,000 features. With this representation, you can easily identify the feature you’re looking for (sa y, a particular eye, or a sign board). That was an overview of the entire algorithm. Over the next few days, I’ll go through each step in detail. Finally, I’ll show you how to implement SIFT in OpenCV! What do I do with SIFT features? After you run through the algorithm, you’ll have SIFT features for your image. Once you have these, you can do whatever you want. Track images, detect and identify objects (which can be partly hidden as well), or whatever you can think of. We’ll get into this later as well. But the catch is, this algorithm is patented. >.< So, it’s good enough for academic purposes. But if you’re looking to make something commercial, look for something else! [Thanks to aLu for pointing out SURF is patented too] 1. Constructing a scale space Real world objects are meaningful only at a certain scale. You might see a sugar cube perfectly on a table. But if looking at the entire milky way, then it simply does not exist. This multi-scale nature of objects is quite common in nature. And a scale space attempts to replicate this concept

SIFT 特征提取算法详解

SIFT 特征提取算法总结主要步骤 1)、尺度空间的生成； 2)、检测尺度空间极值点； 3)、精确定位极值点； 4)、为每个关键点指定方向参数； 5)、关键点描述子的生成。 L(x,y,σ), σ= 1.6 a good tradeoff

D(x,y,σ), σ= 1.6 a good tradeoff

关于尺度空间的理解说明：图中的2是必须的，尺度空间是连续的。在 Lowe 的论文中，将第0层的初始尺度定为1.6，图片的初始尺度定为0.5. 在检测极值点前对原始图像的高斯平滑以致图像丢失高频信息，所以Lowe 建议在建立尺度空间前首先对原始图像长宽扩展一倍，以保留原始图像信息，增加特征点数量。尺度越大图像越模糊。 next octave 是由first octave 降采样得到（如2），尺度空间的所有取值，s为每组层数，一般为3~5 在DOG尺度空间下的极值点同一组中的相邻尺度（由于k的取值关系，肯定是上下层）之间进行寻找

在极值比较的过程中，每一组图像的首末两层是无法进行极值比较的，为了满足尺度变化的连续性，我们在每一组图像的顶层继续用高斯模糊生成了 3 幅图像，高斯金字塔有每组S+3层图像。DOG金字塔每组有S+2层图像.

If ratio > (r+1)2/(r), throw it out (SIFT uses r=10) 表示DOG金字塔中某一尺度的图像x方向求导两次通过拟和三维二次函数以精确确定关键点的位置和尺度（达到亚像素精度）？

直方图中的峰值就是主方向，其他的达到最大值80%的方向可作为辅助方向 Identify peak and assign orientation and sum of magnitude to key point The user may choose a threshold to exclude key points based on their assigned sum of magnitudes. 利用关键点邻域像素的梯度方向分布特性为每个关键点指定方向参数，使算子具备旋转不变性。以关键点为中心的邻域窗口内采样，并用直方图统计邻域像素的梯度方向。梯度直方图的范围是0～360度，其中每10度一个柱，总共36个柱。随着距中心点越远的领域其对直方图的贡献也响应减小.Lowe论文中还提到要使用高斯函数对直方图进行平滑，减少突变的影响。

补充全排列算法C语言实现

字符串全排列算法C语言实现问题描述：输入一个字符串，打印出该字符串中字符的所有排列。输入: 123 输出: 123 132 213 231 312 321 问题分析：现象分析：这种问题，从直观感觉就是用递归方法来实现即：把复杂问题逐渐简单化，进而得出具体代码实现。（如何发现一个问题可以使用递归方式来解决？经分析可以发现：M 个数的排列方法与N(N>M)个数的排列方法没有区别，处理方法与数据的个数没有关系。处理方法的一致性，就可以采用递归）。 3个数（123）排列，第一位1不动，剩下两个数（23）的排列，只要相互颠倒一下，就可以出现关于1开头的所有排列 123 132 把2换到第一位，保持不动，剩下的两个数（13）的排列，只要相互颠倒一下，就可以出现关于2开头的所有排列 213 231 同理，把3换到第一位，可得到 312 321 扩展：把3个数的所有排列，前面加一个4，就可以得到关于4开头的所有的排列4123 4132 4213 4231 4312 4321 若把4与后续数据中的任意一个数据交换，通过完成对后续三个数的全排列，就可以得到相应的数开头的四数的排列。总结：对4个数的排列，可以转换成首位不动，完成对3个数的排列对3个数的排列，可以转换成首位不动，完成对2个数的排列对2个数的排列，可以转换成首位不动，完成对1个数的排列对于1个数，无排列，直接输出结果算法实现说明：

对n个数的排列，分成两步：（1）首位不动，完成对n-1个数的排列，（2）循环将后续其他数换到首位，再次进行n-1个数的排列注：排列完成后，最终的串要与原串相同 C语言代码实现： #include #include //将串左移一位，首位存到末尾。 void shift( char *s ) { if ( !s||!s[0] ) return ; //security code . null string char ch=s[0]; int i=0; while( s[++i] ) s[i-1]=s[i] ; s[i-1]=ch; } //本函数对于一个已排序好的数据进行全排列 void permutation(char list[], int head ) { int i,len; len=strlen(list) ; if ( len-head == 1 ) //后续没有再排列的，则输出排列数据 { printf( "%s\n", list ); } else { for (i = k; i

SIFT算法分析

SIFT算法分析 1 SIFT 主要思想 SIFT算法是一种提取局部特征的算法，在尺度空间寻找极值点，提取位置，尺度，旋转不变量。 2 SIFT 算法的主要特点： a)SIFT特征是图像的局部特征，其对旋转、尺度缩放、亮度变化保持不变性，对视角变化、仿射变换、噪声也保持一定程度的稳定性。 b)独特性(Distinctiveness)好，信息量丰富，适用于在海量特征数据库中进行快速、准确的匹配。 c)多量性，即使少数的几个物体也可以产生大量SIFT特征向量。 d)高速性，经优化的SIFT匹配算法甚至可以达到实时的要求。 e)可扩展性，可以很方便的与其他形式的特征向量进行联合。 3 SIFT 算法流程图：

4 SIFT 算法详细 1）尺度空间的生成尺度空间理论目的是模拟图像数据的多尺度特征。高斯卷积核是实现尺度变换的唯一线性核，于是一副二维图像的尺度空间定义为： L( x, y, ) G( x, y, ) I (x, y) 其中G(x, y, ) 是尺度可变高斯函数，G( x, y, ) 2 1 2 y2 (x ) 2 e / 2 2 （x，y）是空间坐标，是尺度坐标。大小决定图像的平滑程度，大尺度对应图像的概貌特征，小尺度对应图像的细节特征。大的值对应粗糙尺度(低分辨率)，反之，对应精细尺度(高分辨率)。为了有效的在尺度空间检测到稳定的关键点，提出了高斯差分尺度空间（DOG scale-space）。利用不同尺度的高斯差分核与图像卷积生成。 D( x, y, ) (G( x, y,k ) G( x, y, )) I ( x, y) L( x, y,k ) L( x, y, ) DOG算子计算简单，是尺度归一化的LoG算子的近似。图像金字塔的构建：图像金字塔共O组，每组有S层，下一组的图像由上一组图像降采样得到。图1由两组高斯尺度空间图像示例金字塔的构建，第二组的第一副图像由第一组的第一副到最后一副图像由一个因子2降采样得到。图2 DoG算子的构建：图1 Two octaves of a Gaussian scale-space image pyramid with s =2 intervals. The first image in the second octave is created by down sampling to last image in the previous