一.算法介绍 SIFT算法（Scale-Invariant feature transform，尺度不变特征变换）通过在图像中提取独特性不变特征，可以实现物体或场景在不同视角下的可靠匹配。其提取的特征对于图像缩放、旋转和一定范围内的三维仿射变换、噪声叠加、光照变化均具有不变性。由于特征的高度独特性，场景中的每一个特征都有很大的可能在由多幅图像提取的特征数据库中得到正确的匹配结果。因此使用这些特征可以用于物体识别。识别首先需将每个特征与由已知物体组成的特征数据库进行快速最邻匹配，然后通过霍夫变换确定属于同一个物体的特征点对，最后通过对一致性姿态参数的最小二乘法进行验证。这种方法可以在噪声和光照影响下获得稳定的识别效果，而不失实时性。 使用一系列局部关键点进行图像匹配最早可追溯至Moravec(1981)，他在立体匹配中使用了角点检测。Moravec detector 在1988年被Harris和Stephens改进，Harris又在1992年改进用于动作追踪，自此Harris corner detector被广泛应用。1997年Schmid和Mohr做出了开拓性的工作，证明局部的不变特征可以应用于图像识别。但是Harris corner detector对图像尺度改变非常敏感，所以对不同大小的图像匹配效果不佳。 Lowe,D.G.1999年发表论文《Object recognition from local scale-invariant features》，提出了一种具有尺度不变性的特征。之后，又有一系列研究试图将这种不变性扩展至对全部的仿射变换有效。在2004年Lowe发表的论文《Distinctive Image Features from Scale-Invariant Keypoints》中，其SIFT特征可以适用于一定范围内的仿射变换。Y.Ke和R.Sukthankar 2004年在《PCA-SIFT：A More Distinctive Representation for Local Image Descriptors》中提出PCA-SIFT，将原算法中的直方图法换成主元分析法(Principle Component Analysis)，实现描述子的降维。Herbert Bay 于2006年提出SURF加速稳健特征算法，实现算法加速，并于2008年发表论文《Speeded-Up Robust Features(SURF)》。 很多其他的特征类型也被用于图像识别中，例如图像轮廓、图像相位、颜 {MOD}等。局部点特征与这些特征相结合可以增强识别的鲁棒性，未来的图像识别系统很可能采用多特征结合的方法。

二. OPENSIFT编译运行

David Lowe，SIFT算法的提出者，在论文之外提供了其算法的示例程序——SIFT demo program (Version 4, July 2005)，可惜其关键部分没有给出源码。Rob Hess用C语言写了开源的SIFT实现——OpenSIFT，且具有与原作者示例相差不多的运行性能。由于SIFT算法实现比较复杂，局囿于本人能力与时间问题，本节算法实现细节将主要结合Lowe论文和OpenSIFT源码进行，本人没有重新编写代码，仅对代码进行注释。 虽然有开源代码可直接下载，但作者用gcc编译，代码还依赖于另外两个开源项目OpenCV和GTK+，不能在我熟悉Windows环境下直接使用，我将其加入到Microsoft Visual Studio 2012的项目中稍作改动并编译通过。

2.1 新建VS项目

首先在VS2012中新建项目，模板选择Visual C++->Win32->Win32控制台应用程序，名称可填写SIFT，因只有一个项目，可以将为解决方案创建目录取消勾选。将下载的OpenSIFT中的include和 src目录拷贝至项目目录下，在VS的解决方案资源管理器下，分别在头文件和源文件上右击选择添加->现有项，添加include中的h文件和src中的c文件。注意src中的match.c、dspfeat.c和siftfeat.c中都含有main函数，分别实现两张图片匹配、读取存在文件中的图像特征并显示和计算输入图像的SIFT特征，源文件中只能添加三个中的一个。我这里选择添加match.c 这时直接点击本地Windows调试器会有一大堆错误，还需要进行下面的配置。

2.2 OpenCV配置

OpenCV是一个计算机视觉库。在OpenCV中文网站下载最新的OpenCV for Windows，当前是Version 2.4.4，并参考安装文档安装，可参照其中的VC 2010 Express下安装OpenCV2.4.3，基本上都是一样的，配置时直接对之前建立的SIFT项目进行即可。有一些配置对于32位系统和64位系统是不一样的，个人建议跟我一样非专业的，64位系统也还是按照32位系统来配置，这样兼容性好一些，后面也少一些麻烦。另外文中说要将TBB目录加入环境变量Path中，我没有找到，不加好像也没什么问题。然后要在配置属性->链接器->输入->附加依赖项中增加lib文件，我的版本是2.4.4，要将文中的243改成244。可以直接粘贴下面的内容在编辑框中。 Debug配置下： opencv_calib3d244d.lib;opencv_contrib244d.lib;opencv_core244d.lib;opencv_features2d244d.lib;opencv_flann244d.lib;opencv_gpu244d.lib;opencv_highgui244d.lib;opencv_imgproc244d.lib;opencv_legacy244d.lib;opencv_ml244d.lib;opencv_objdetect244d.lib;opencv_ts244d.lib;opencv_video244d.lib; Release配置下： opencv_contrib244.lib;opencv_core244.lib;opencv_features2d244.lib;opencv_flann244.lib;opencv_gpu244.lib;opencv_highgui244.lib;opencv_imgproc244.lib;opencv_legacy244.lib;opencv_ml244.lib;opencv_objdetect244.lib;opencv_ts244.lib;opencv_video244.lib;

2.3 GTK+配置

GTK+是一个图形工具包，在OpenSIFT主要用它来实现窗口和画一些线。在The GTK+ Project这个网站上下载GTK+，我选择Windows(32-bit)版本，当前是Version 2.24，里面有很多包可以下载，像我一样不知道怎么选的可以点all-in-one bundle.下载下来是一个zip压缩包，里面有一个txt说明。把压缩包解压到一个空目录下，比如我放在C:Program Files (x86)GTK 下，然后将bin的路径添加进系统环境变量PATH中。之后按说明验证，Win+R输入cmd运行，在cmd中输入“pkg-config --cflags gtk+-2.0” ，会有一些输出，输入 “gtk-demo” ，会出现一个示例，演示GTK+的一些功能控件。 接下来就跟OpenCV一样，要在VS2012的项目中进行一番配置了。在CMD中输入运行“pkg-config --cflags --libs gtk+-2.0”,可以看到需要包含的目录和链接库。可以将这些输出导入txt文件中，运行“pkg-config --cflags --libs gtk+-2.0 > G:gtk.txt”,打开G:gtk.txt，内容如下： Files (x86)/GTK/include/gtk-2.0 -mms-bitfields (x86)/GTK/lib/gtk-2.0/include (x86)/GTK/include/atk-1.0 (x86)/GTK/include/cairo (x86)/GTK/include/gdk-pixbuf-2.0 (x86)/GTK/include/pango-1.0 (x86)/GTK/include/glib-2.0 (x86)/GTK/lib/glib-2.0/include (x86)/GTK/include (x86)/GTK/include/freetype2 (x86)/GTK/include/libpng14 -IC:/Program Files (x86)/GTK/lib -LC:/Program -lgtk-win32-2.0 -lgdk-win32-2.0 -latk-1.0 -lgio-2.0 -lpangowin32-1.0 -lgdi32 -lpangocairo-1.0 -lgdk_pixbuf-2.0 -lpango-1.0 -lcairo -lgobject-2.0 -lgmodule-2.0 -lgthread-2.0 -lglib-2.0 -lintl 然后根据这个来添加配置。在VS项目属性的VC++目录->包含目录中添加/GTK/include/gtk-2.0 到 /GTK/include/libpng14的这些路径，注意要用带盘符的完整路径，中间的那个-mms-bitfields不用管它。-I后面的C:/Program Files (x86)/GTK/lib要添加在库目录中。再后面的-l是链接库的名字，把这一串”in32-2.0.lib;gdk-win32-2.0.lib;atk-1.0.lib;gio-2.0.lib;pangowin32-1.0.lib;gdi32.lib;pangocairo-1.0.lib;gdk_pixbuf-2.0.lib;pango-1.0.lib;cairo.lib;gobject-2.0.lib;gmodule-2.0.lib;gthread-2.0.lib;glib-2.0.lib;intl.lib;”添加进配置属性->链接器->输入->附加依赖项 中就行了。

2.4 代码修改

这时再重新生成解决方案，还是会有很多错误。需要逐一解决。 1) 无法打开包括文件：因为头文件都放在include文件夹中，编译器无法定位。在项目属性->配置属性->VC++目录->包含目录中加入$(ProjectDir)include 2) 在”inline”之后应输入’(’：在utils.h中加 #define inline __inline 3) “srandom”,"random”未定义：将xform.c中的srandom改为srand，random改为rand 4) 应输入常量表达式：sift.c中，应该是不支持这样实现动态数组。 const int _intvls = intvls; double sig[_intvls+3], sig_total, sig_prev, k; 这两行改为 double *sig = (double *)calloc(intvls + 3,sizeof(double)); double sig_total, sig_prev, k; 再在最后return前加上free(sig); 5) “M_PI”未声明：发生在imgfeatures.c中。在imgfeatures.h中加入 #define M_PI 3.14159265358979323846 6) “snprintf”未定义：在utils.c中，将snprintf改为_snprintf 上述修改之后，再生成解决方案应该就可以成功了。但程序运行还需要命令行参数。点击调试->SIFT属性，在配置属性->调试->命令参数栏中输入”beaver.png beaver_xform.png”，确定。然后再点启动调试就可以看到运行效果了，当然之前要将这两幅图像放在工程目录下。程序运行匹配效果如图1： 图1：OpenSIFT图像匹配效果 为了直观显示每幅图的关键点信息，我们把siftfeat.c中的一部分代码复制到match.c的main函数中，并适当修改。 if( display ) { draw_features( img1, feat1, n1 ); display_big_img( img1, argv[1] ); cvWaitKey( 0 ); } 运行效果如图2：   图2：SIFT关键点 图中射线端点为关键点位置，此外还直观表示了关键点的尺度和方向信息。

三. 算法详解

大致按照Lowe论文顺序讲解

3.1 尺度空间极值点检测(Detection of scale-space extrema)

要实现图像匹配，我们想在构成图像的所有点中通过某种方法选取出一些具有代表性的点，重要的是这些特殊点在图像经过伸缩、旋转、叠加噪声和其他一些变化后仍能被提取出来——即具有不变性，这样我们就可以利用这些具有不变性的特殊点来进行图像的匹配。Lowe提出了一种寻找特殊点的方法，用这种方法找到的点对于图像的尺度变化具有不变性。 前面提到了，最早的图像匹配使用的是角点检测（corner detector）。应该说，图像中物体的边缘是一个比较好的特征。比较常用的边缘检测方法有一阶导数边缘检测、二阶导数边缘检测。拉普拉斯边缘检测就属于二阶导数边缘检测。拉普拉斯算子为 Marr和Hildrith提出了高斯拉普拉斯边缘检测算子（LOG算子，表示为 ），在拉普拉斯操作之前先进行高斯平滑操作，以抑制高频噪声。高斯函数表示为G(x,y,σ)，与原始图像做卷积后，相当于对图像进行了低通滤波，且σ越大，则滤波器截止频率越低，图像也因失去更多的高频信息而更模糊。Lindeberg又提出将上述结果乘以系数后可获得尺度不变性，即尺度归一化的高斯拉普拉斯算子，。为什么这样可以获得尺度不变性，我也不明白，需查看论文。但大致的原理可以从后面的高斯差分中看出，即σ的增大和图像尺寸的减小是有某种联系的，构建图像金字塔的过程即是在所有尺度特征中离散采样，这样可以保证两幅不同尺寸的图像的金字塔在某些层是可以相匹配的。）） 这样处理后得到边缘检测图像，我们可以在其中寻找局部极值点来作为这幅图像的特征点。如果令中的σ不断变大（为了计算机实现，只能取离散的σ值，对应于图3的Scale），就可以得到很多输出图像，这样局部极值点的寻找就扩展至三维空间(x,y,σ)。正如图3所示，每个像素点X需要与上下及周围共26个邻近点O比较来确定局部最大最小极值。 图3：极值点检测 应该说，这样已经可以得到所谓的尺度空间极值点了，只不过直接运用的计算量比较大。Lowe采用了一些方法来减少计算量。 首先是用高斯差分(DOG)来近似。Lowe说可以参考热扩散方程得到，这个公式应该比较关键，但我尚不清楚为何     左边的高斯偏导数可以用差分来近似： 其中k趋近于1就可以取等号了。所以， 从中可以看出，高斯差分与尺度归一化的高斯拉普拉斯算子仅相差一个比例因子(k-1)。如果我们取k为定值，则在高斯差分中找到的极值点位置大致相当于在中的位置。 下一个减少计算量的方法是分组（Octave），即通过缩小图片尺寸来减少卷积计算。后面会详细解释。 图4：分组高斯差分 上面提到了高斯差分，要实际应用有几个参数要确定。初始的σ值，尺度扩大的比例系数k值，还有到底算多少层结束。后面的说明用σ容易混乱，以后的σ就当作一个常数，令，高斯函数为G(x,y,Scale) 然后这样分组，Scale=σ~2σ的为一组，Scale=2σ~4σ的为一组……。我们知道，随着Scale的变大，高斯卷积核矩阵也需要取大，这样计算量就会增加。分组时每一组图像的尺寸取为上组的一半，这样就避免了卷积核矩阵的扩大。例如原始图像与2σ高斯核卷积结果隔行隔列抽取出的数据，与原始图像隔行隔列抽取的数据与σ高斯核卷积得到的结果是一样的，但后者的计算量明显要少得多。当然这样也会造成一些信息的丢失，但高Scale下的图像本来就很模糊了，所以影响不大。 每一组中取几个离散值呢？Lowe表示取s=3个比较经济，当然，如果不考虑计算速度的话，取多点可以得到更多的特征点。这样可以得到k，即k=2^(1/s)=2^(1/3).因为不同组的高斯差分尺寸大小不同，不能比较极值，所以每一组需额外几张图片来保证不同组高斯差分间的连续性。每组共需生成s+3=6张高斯处理过的图像。从图5中可以看得比较清楚 图5：DOG空间生成示意 从图中看出，第1组由σ~(k^5)σ的6幅高斯化图像，产生了5幅高斯差分（DOG）图像，去除首尾不能比较极值的，可以在中间3幅DOG图像（对应kσ~(k^3)σ）中找到特征极值点。第2组则由(k^3)σ~(k^8)σ的6幅高斯化图像最终得到有极值点的(k^4)σ~(k^6)σ与上组相连的3幅DOG图像。注意图像下采样导致尺寸减半后，实际计算时用的Scale参数也要相应减半。因此每一组第1幅高斯图像可直接由上层第4幅图像尺寸减半下采样获取（因为k^3=2）。 初始的σ值，Lowe在综合考虑后取为1.6。但他说这样就不能充分利用图像中的高频信息，所以将原始图像以线性插值的方法扩大两倍作为第一层，这一举措大大增加了关键点的数量。另外假设原始图像由于噪声等影响已经相当于被σ=0.5的高斯函数模糊化过了，因此为了达到σ=1.6的高斯模糊效果，只需再使用一个σ=√((1.6)^2-(0.5)^2 )的高斯模糊就行了。（由于大的高斯模糊的半径是所用多个高斯模糊半径平方和的平方根，见维基百科：高斯模糊） 计算终止条件，OpenSIFT自己注释是说图像尺寸每一维不能小于4像素，但程序中这样计算貌似是不小于8像素，不过这个参数应该是无伤大雅的。如下： octvs = log( MIN( init_img->width, init_img->height ) ) / log(2) - 2; sift.h 中定义的参数 /** default number of sampled intervals per octave 每一组分的层数 */ #define SIFT_INTVLS 3 /** default sigma for initial gaussian smoothing sigma初始值*/ #define SIFT_SIGMA 1.6 /** double image size before pyramid construction? 是否扩大两倍图像作为初始层*/ #define SIFT_IMG_DBL 1 /* assumed gaussian blur for input image 假定的初始模糊度*/ #define SIFT_INIT_SIGMA 0.5 sift.c 中sift_features 默认使用了这些参数。 /* build scale 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; //看一共能分几组(Octave) gauss_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 ); //在高斯差分金字塔中寻找尺度空间极值点 中间的一些函数我都在源码中增添了中文注释。其实Rob Hess的代码很清晰，只要明白了原理就不难读懂。 寻找极值点时用到了一些限制参数，将在下文中讲解。

3.2 关键点精确定位(Accurate keypoint licalization)

前面已经在尺度中间中找到了一些极值关键点，但这些点并不是都能使用，要添加限制条件去除一些不好的。另外作者还提供了将关键点坐标提高至亚像素级精度的方法。 前面我们已经在三维尺度空间D(x,y,σ)中找到了极值点，但由于空间是离散化的，极值点的位置精度也被限制于离散点上。我们通过在离散点进行泰勒二次展开来提高精度。 上述公式来自论文。注意是一个向量。这种多维的泰勒展开我不太常用，而且作者的公式我觉得表达的不是很准确，全部用X，比较混乱。可以对比普通一维的来理解 总之，极值点应该出现在导数为0时。对D(X)求导，令其为0 得准确极值点与中间离散点的偏差 相关程序如下，偏差计算在sift.c中interp_step函数内实现： dD = deriv_3D( dog_pyr, octv, intvl, r, c ); //计算了dD/dx,dD/dy,dD/d(sigma),返回为列向量dD/dX H =hessian_3D( dog_pyr, octv, intvl, r, c ); //用离散值近似计算出三维hessian矩阵,即公式中d2D/dX2 H_inv =cvCreateMat( 3, 3, CV_64FC1 ); cvInvert( H, H_inv, CV_SVD ); //矩阵求逆 cvInitMatHeader( &X, 3, 1,CV_64FC1, x, CV_AUTOSTEP ); cvGEMM( H_inv, dD, -1, NULL, 0, &X, 0 ); //计算结果为-(d2D/dX2)^(-1)*(dD/dX),即off_X 如果得到的偏差X ̂在任一维上大于0.5，则说明准确极值点离另一个采样离散点更近。在这种情况下，需要通过在那个采样点上进行泰勒展开再进行一次计算。最后将偏差加在采样点坐标上以获得更为精确的极值点估计坐标。 最后极值点处的具体数值为 这部分在interp_contr中 dD = deriv_3D( dog_pyr, octv, intvl, r, c ); //dD/dX cvGEMM( dD, &X, 1, NULL, 0, &T, CV_GEMM_A_T );//(dD^T/dX)*off_X return pixval32f( dog_pyr[octv][intvl], r, c ) + t[0] * 0.5; //D+(dD^T/dX)*off_X/2 前面说了，D(X)近似于归一化的高斯拉普拉斯，较低的值代表着较低的对比度。而低对比度的点对噪声比较敏感，因此论文中将|D(x)| 小于0.03的极值点都抛弃了（假定图像像素值变化范围为[0,1]）。程序中取的0.04，并在interp_extremum中实现了上述功能 /** default threshold on keypoint contrast |D(x)| 对比度最小值限制*/ #define SIFT_CONTR_THR 0.04 我们已经知道，高斯差分在图像边缘处会有比较大的值。但并非所有的边缘点都可以作为稳定的好的特征关键点。我们需要排除这样一些点，在穿过边缘的方向上变化较大，而在相垂直的方向上变化很小。（剩下的应该就是包括角点这样的点了） 这里需要用到2×2的海森(Hessian)矩阵H： 通过The Hessian matrix and its eigenvalues可知，矩阵的两个特征值分别表征了两个D在两个相垂直方向上的变化率。设α为幅值较大的特征值，β为幅值较小的特征值。取r为两个特征值之比，α=rβ 根据上面所述，我们需要排除r较大的点。由于我们只关心r，可以用以下方法来避免特征值的计算。 根据矩阵知识，有： 特征值为负直接抛弃，其他有， 易知r≥1，此时((r+1)^2)/r是单调递增的。因此欲限制r小于某个阈值，只需满足 在论文中，限制r不能大于10.这部分实现在is_too_edge_like中 /** default threshold on keypoint ratio of principle curvatures 关键点处曲率比值最大值限制*/ #define SIFT_CURV_THR 10 /* principal curvatures are computed using the trace and det of Hessian */ d =pixval32f(dog_img, r, c); dxx = pixval32f( dog_img, r, c+1 ) + pixval32f( dog_img, r, c-1 ) - 2 * d;dyy = pixval32f( dog_img, r+1, c ) +<