结合OPENSIFT源码详解SIFT算法

一.算法介绍

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)》。

很多其他的特征类型也被用于图像识别中，例如图像轮廓、图像相位、颜色等。局部点特征与这些特征相结合可以增强识别的鲁棒性，未来的图像识别系统很可能采用多特征结合的方法。

二. 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是链接库的名字，把这一串”gtk-win32-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

7) 无法找到头文件dirent.h：在网上下载dirent.h - dirent API for Microsoft Visual Studio文件，放在include文件夹下，更改utils.h中<dirent.h>为”dirent.h”

上述修改之后，再生成解决方案应该就可以成功了。但程序运行还需要命令行参数。点击调试->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 ) + 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; reject feature */
  if( det <= 0 )
    return 1;
  if( tr * tr / det < ( curv_thr + 1.0 )*( curv_thr + 1.0 ) / curv_thr )
    return 0;
  return 1;

3.3 方向分配(Orientation assignment)

好了，我们已经确定了一些关键点。接下来需要一些参数来描述这个点，即关键点描述子(keypoint descriptor)。首先需要一个方向参数，以获取对图像旋转的不变性。

梯度方向是很容易想到的，但仅仅取这个点上的梯度方向，会不太稳定。Lowe采用取关键点周围一块区域内的梯度信息综合，作为关键点的方向。

方向计算是在高斯模糊后的图像L上进行的，且尺度与关键点所在尺度相同。二维空间上每一点的梯度（包括幅值m(x,y)和角度θ(x,y)）可以用以下方程表示，在calc_grad_mag_ori中实现:

将这个区域内所有点的梯度方向绘制成直方图，其中直方图将360度分成36份。每一个加入直方图中的方向需要根据其梯度幅值和到中央关键点距离进行权重（距离影响可用一个具有1.5σ的高斯函数实现，越远影响越小，σ为关键点所在尺度）。这部分在ori_hist中实现。根据高斯函数的3σ原理，距中心点3*1.5*σ时权值已经很小，可忽略不计，由此确定进行统计计算的区域大小。

最后得到的直方图反映了以关键点为中心的这部分图像区域的梯度方向分布。使用前先对直方图进行一次平滑操作，每相邻三个数据按[0.25,0.5,0.25]系数进行加权。具体见函数smooth_ori_hist。可以用直方图中最大值对应的角度来作为关键点的方向参数。对于大于最大值的80%的那些角度，也可以用来描述关键点的方向。据说，只有大约15%的关键点会有多重方向，但它们对于匹配的稳定性却意义重大。最后，直接使用直方图的角度将只有10度的角度分辨率，因此使用这个角度附近的共3个数据，采用抛物线插值方法来确定最终角度。这部分在calc_feature_oris中实现。

图6：方向直方图生成

一些参数：

/* default number of bins in histogram for orientation assignment 分配方向时直方图份数*/
#define SIFT_ORI_HIST_BINS 36
/* determines gaussian sigma for orientation assignment 分配方向时高斯权重参数因子*/
#define SIFT_ORI_SIG_FCTR 1.5
/* determines the radius of the region used in orientation assignment 分配方向时计算区域大小，根据3sigma外可忽略*/
#define SIFT_ORI_RADIUS 3.0 * SIFT_ORI_SIG_FCTR
/* orientation magnitude relative to max that results in new feature 分配方向时副方向相对主方向幅值的比例阈值*/
#define SIFT_ORI_PEAK_RATIO 0.8

主要代码：

hist= ori_hist( gauss_pyr[ddata->octv][ddata->intvl],  //生成方向分配所需直方图
       ddata->r, ddata->c, SIFT_ORI_HIST_BINS,
       cvRound( SIFT_ORI_RADIUS * ddata->scl_octv ),
       SIFT_ORI_SIG_FCTR * ddata->scl_octv );
for( j = 0; j < SIFT_ORI_SMOOTH_PASSES; j++ )
    smooth_ori_hist( hist, SIFT_ORI_HIST_BINS );    //直方图平滑
omax = dominant_ori( hist, SIFT_ORI_HIST_BINS );    //直方图峰值
add_good_ori_features( features, hist, SIFT_ORI_HIST_BINS,  //增加其他方向并进行二次插值
         omax * SIFT_ORI_PEAK_RATIO, feat );

3.4 局部图像描述子(The local image descriptor)

前面已经获取了关键点的位置、尺度和方向。如果想要匹配，我们可以直接将两幅图像关键点附近的区域进行相关性分析计算，根据相似度来判断是否匹配。

这样做有两个问题，一个是计算量太大，一个是这种方法对于仿射变换很敏感，因为仿射变换后图像由于扭曲，相关性已经降低了。

直接进行图像相关分析计算量大，我们可以折中只取几个参数来描述这个小区域，只比较这几个参数的差别。比如说前面分配方向时用到的梯度方向分布直方图其实就可以作为一种描述方法。我们可以将图像旋转至插值计算后的关键点主方向后再进行一次直方图的计算，这样可以得到具有旋转不变性的36个描述参数。只进行36个参数的比较，计算量就大大降低了。

对于第2个问题的解决，Lowe说也是从别人的方法中得到了启示。我们从不同角度观察物体时，物体上的同一个点在物体投影中的相对位置可能会发生变化。也就是说，如果我们以关键点附近的多个点为中心分别计算梯度方向分布直方图，那么即使两幅以不同角度投影的图像在它们的精确关键点位置处不能很好地匹配，它们也极有可能在关键点附近点上得到较好的匹配。这样也增加了描述参数的数量，减少了误匹配的可能。

论文最后采用了4×4个方向直方图的方案，相当于分别以16个点为中心计算，而每个直方图分为8个大的方向，最后形成一个4×4×8=128维的描述向量。

/** default width of descriptor histogram array 描述子直方图矩阵默认尺寸*/
#define SIFT_DESCR_WIDTH 4
/** default number of bins per histogram in descriptor array 描述子直方图默认方向份数*/
#define SIFT_DESCR_HIST_BINS 8

论文中说这个描述向量是由16×16采样像素点计算得来的，但程序中和网上下的PPT似乎并没有严格按照论文。一直没看明白，后来通过图片搜索找到这篇文章——【OpenCV】SIFT原理与源码分析：关键点描述，这位同学解释得比较清楚。

图7：描述子采样区域大小示意

这里描述子的采样区域是与关键点所在尺度有关的。与上一节一样，各点梯度的计算是在与关键点对应的高斯模糊图像上进行的。将关键点附近划分成d×d个子区域（默认d=4），每个子区域宽度为hist_width=3σ个像元。考虑到实际计算时需要双线性插值，故计算的图像区域宽度为hist_width*(d+1)，再考虑旋转，则实际运算时需要进行判断的图像区域宽度为hist_width*(d+1)*sqrt(2).

程序中需要对所有图7最外面实线大方框中的点进行判断，但最后直方图区域仅包含在里面实线小框（可以旋转）中的点。虚线框与小实线框间的点对于位于实线框边缘的的直方图也是有贡献的。这样理解之后再来看descr_hist中的这一段程序就比较清楚了。

图8：区域坐标轴旋转

首先需要旋转坐标轴，使旋转后x方向与关键点方向相同。点在旋转后坐标系下的新坐标表示为:

图9：按子区域生成直方图矩阵

下面考虑每一点对于直方图的贡献大小，即权重。首先要给离关键点较近的点以较高的权值。这可以通过在梯度幅值上乘以一个高斯函数来实现，高斯函数的σ参数取为描述子区域宽度的一半，即为0.5*d*hist_width. 除此之外，这次要生成多个相邻区域的直方图，我们希望这些相邻直方图不要变化太大。所以采用双线性插值，每一个点对离它最近的4个子区域的直方图都会有影响，影响程度视其离该区域中心点（即图9中网格的交点，除去边缘的交点，共有16个交点，即16个子区域）在x、y两个方向上的距离。可以在每一维上简单地乘以1-dis,dis是以hist_width为单位长度的距离，注意当dis>1时要乘的是0。另外，角度不是像方向分配时那样直接四舍五入，而是分配到距它最近的两个方向上。综上，可得

图9中每个子区域用8个箭头代表直方图的8个方向，箭头的长度对应直方图在这个方向上的值。

以下程序生成描述子直方图：

static double*** descr_hist( IplImage* img, int r, int c, double ori,
                 double scl, int d, int n )
{
  double*** hist;
  double cos_t, sin_t, hist_width, exp_denom, r_rot, c_rot, grad_mag,
    grad_ori, w, rbin, cbin, obin, bins_per_rad, PI2 = 2.0 * CV_PI;
  int radius, i, j;

  hist = calloc( d, sizeof( double** ) );   //共有d*d个直方图，每个直方图包含n个数据
  for( i = 0; i < d; i++ )
    {
      hist[i] = calloc( d, sizeof( double* ) );
      for( j = 0; j < d; j++ )
    hist[i][j] = calloc( n, sizeof( double ) );
    }
  
  cos_t = cos( ori );
  sin_t = sin( ori );
  bins_per_rad = n / PI2;   //直方图每份覆盖角度范围
  exp_denom = d * d * 0.5;  //高斯系数取为描述子宽度一半。本来应该为(0.5*d*hist_width)^2更好理解，这是为了后面w计算方便
  hist_width = SIFT_DESCR_SCL_FCTR * scl;   //3sigma，scl为关键点的尺度.
  radius = hist_width * sqrt(2) * ( d + 1.0 ) * 0.5 + 0.5;  //描述子覆盖区域半径，最后加0.5以实现四舍五入
  for( i = -radius; i <= radius; i++ )  //包含在大实线框中的点都需要判断
    for( j = -radius; j <= radius; j++ )
      {
    /*
      Calculate sample's histogram array coords rotated relative to ori.
      Subtract 0.5 so samples that fall e.g. in the center of row 1 (i.e.
      r_rot = 1.5) have full weight placed in row 1 after interpolation.
    */
    c_rot = ( j * cos_t - i * sin_t ) / hist_width;//计算在旋转坐标系下的坐标表示（点不动，坐标轴转），以小区域宽度为步长单位
    r_rot = ( j * sin_t + i * cos_t ) / hist_width;
    rbin = r_rot + d / 2 - 0.5; //rbin是在新坐标系下的坐标，这里单位长度为hist_width，这里又进行平移，使原点变到（1.5，1.5)
    cbin = c_rot + d / 2 - 0.5; //加d/2后将小实线框的左上角平移至原点，再减0.5又回去一点
        //下面是看这个点是否在虚线框内。虚线框宽度在这个坐标系下为d+1，如果rbin不减0.5就更容易理解了，这样rbin>-0.5 && rbin<d+0.5，rbin的取值范围恰是d+1
    if( rbin > -1.0  &&  rbin < d  &&  cbin > -1.0  &&  cbin < d )  
      if( calc_grad_mag_ori( img, r + i, c + j, &grad_mag, &grad_ori )) //计算此点的梯度方向，并化成相对于关键点方向的角度值。grad_mag梯度幅值
        {
          grad_ori -= ori;
          while( grad_ori < 0.0 )
        grad_ori += PI2;
          while( grad_ori >= PI2 )
        grad_ori -= PI2;
          
          obin = grad_ori * bins_per_rad;   //角度值映射至0~7之间
          w = exp( -(c_rot * c_rot + r_rot * r_rot) / exp_denom ); //计算以关键点为中心的高斯权重系数,旋转对称，旋转不影响其值。
          //由此公式得来w=exp(-(i^2+j^2)/(2*(0.5*d*hist_width)^2))=exp(-(c_rot^2+r_cot^2)/(0.5*d*d))
          interp_hist_entry( hist, rbin, cbin, obin, grad_mag * w, d, n ); //通过双线性插值计算直方图
        }
      }
  return hist;

static void interp_hist_entry( double*** hist, double rbin, double cbin,
                   double obin, double mag, int d, int n )
{
  double d_r, d_c, d_o, v_r, v_c, v_o;
  double** row, * h;
  int r0, c0, o0, rb, cb, ob, r, c, o;

  r0 = cvFloor( rbin ); //返回不大于参数的最大整数值
  c0 = cvFloor( cbin );
  o0 = cvFloor( obin );
  d_r = rbin - r0;  //得到小数部分
  d_c = cbin - c0;
  d_o = obin - o0;

  /*
    The entry is distributed into up to 8 bins.  Each entry into a bin
    is multiplied by a weight of 1 - d for each dimension, where d is the
    distance from the center value of the bin measured in bin units.
  */
  //这里也可以看出前面rbin、cbin为何减0.5。这样原点上这点的d_r、d_c均为0.5，即原点上这点的方向将被平均分配叠加在它
  //周围4个直方图上面
  for( r = 0; r <= 1; r++ )
    {
      rb = r0 + r;
      if( rb >= 0  &&  rb < d )
    {
      v_r = mag * ( ( r == 0 )? 1.0 - d_r : d_r );  //公式中的1-d_r
      row = hist[rb];   //第rb行上的hist
      for( c = 0; c <= 1; c++ )
        {
          cb = c0 + c;
          if( cb >= 0  &&  cb < d )
        {
          v_c = v_r * ( ( c == 0 )? 1.0 - d_c : d_c );
          h = row[cb]; //h表示第rb行cb列上的hist
          for( o = 0; o <= 1; o++ )
            {
              ob = ( o0 + o ) % n;
              v_o = v_c * ( ( o == 0 )? 1.0 - d_o : d_o );  //角度不是像ori_hist中那样直接四舍五入，而是分配到距它最近的两个方向上
              h[ob] += v_o;
            }
        }
        }
    }
    }
}

下面我们考虑这个描述特征对于光照变化的不变性。首先，将这个128维向量化为单位长度，通过normalize_descr。通过这个归一化操作，将消除等比例光强变化的影响。对于所有光强增加一个常数的情况，由于我们选取的特征是通过差值操作获得的，所以也不会有影响。但是，对于非线性的光照变化，对于某些点的梯度幅值可能会有较大影响，但是对梯度方向的影响可能会比较小。因此，我们限制在向量单位化后，每个值都不能大于0.2，即把大于0.2的值改为0.2，然后重新进行一次单位化。

/* threshold on magnitude of elements of descriptor vector 描述向量元素最大阈值*/
#define SIFT_DESCR_MAG_THR 0.2

3.5 应用于物体识别(Application to object recognition)

首先我们需要一副参考图像，通过上面描述的方法计算这幅图像的所有关键点。然后来了另外一张欲匹配的图像，也计算出这幅图像的所有特征点。对于后者的特征点，我们可以利用欧氏距离在参考图像中找到一个最近点。欧氏距离即128维差值的平方和再开平方。

易知，最近点总是存在的，但不一定是一个正确的匹配点。通过定义一个最近距离的阈值来判断并不是一个好方法。一个更有效的方法是观察最近点距离与次最近点距离的比值。其思想如下，如果正确的匹配点真的存在，那么它应该是唯一的，而且就是最近点，且这个距离应该比它到其它点的距离要小得多，因为其它点都是错误的匹配点，理因有更大的距离。但如果这个点在参考图像中确实没有匹配点，即所有的点都是错误匹配点，则最近点和次最近点的距离值相近的可能性很大，因为所有点差值都不过是高维的随机数。

论文中，作者将所有最近点与次最近点距离之比大于0.8的匹配视为错误的匹配，这样可以减少90%的误匹配，而对正确匹配数量的减少只有不到5%。OpenSIFT默认用了一个更小的阈值，估计是错匹配太多了。

/* threshold on squared ratio of distances between NN and 2nd NN 最近点与次最近点距离之比*/
#define NN_SQ_DIST_RATIO_THR 0.49

在高维空间寻找最近点是一个难题。算法例如k-d树对于超过10维的空间就不能提供有效的加速。因此，作者在k-d树算法基础上得到一种寻找近似最近点的算法——Best-Bin-First(BBF)算法。这部分代码实现在kdtree.c中，使用kdtree_build构建特征kd树，用kdtree_bbf_knn返回2个近似的最近点，再根据这两个值之比判断是否将其作为正确匹配。算法实现细节还未曾研究，以后有时间再看。网上有人写了这部分的注释，有兴趣可以参看——k-d tree的优化查找算法BBF

对于物体识别来说，我们希望通过尽可能少的特征匹配来实现。实际上，如果两幅图像间是仿射变换的关系，且能保证所有的匹配点都是正确的，那么只需要3个匹配对就可以确定这个仿射变换的6个参数。实际应用时，会有大量的匹配对，我们需要从中筛选出属于同一个物体的匹配对，并综合得到仿射变换参数。论文中说如果有效匹配占全部匹配比例小于50%，常用的算法比如RANSAC（随机抽取一致）和Least Median of Squares（最小二乘）表现很差，而实际上这个比例通常小于1%（比例如此小是因为作者考虑的是实际应用的情况，即特征数据库是由很多只含单一物体的小图片的特征组成的，而来匹配的是一副包含很多物体的大图片，对每一个待识别物体来说，与其有关的匹配只占所有匹配的一小部分，如图10），所以作者提出用Hough变换来实现甄别。

图10：复杂场景中的物体识别定位

如果只考虑仅包含一个物体的场景的话，就要简单一些了。由于有效匹配占的比例较高，可以采用RANSAC算法。OpenSIFT在match.c中提供了RANSAC的实现，默认是被注释的。图11是对图1中匹配进一步计算变换参数后的效果。

RANSAC我也不懂，无法讲解。以后有机会再看。

图11：匹配后得到变换参数

3.6 其它图片效果

对于给定的3张船的照片，仅需修改程序的图片名称输入参数，即可获得匹配效果。

图12：img1与img5匹配效果

图13：img1与img6匹配效果

四. SIFT算法的应用及改进

SIFT算法目前已经在很多领域得到应用。包括但不限于：

Ø 物体识别

Ø 机器人定位与导航

Ø 全景图像拼接

Ø 三维场景建模

Ø 动作识别

Ø 指纹与人脸识别

Ø 视频跟踪

Ø ……

SIFT算法也有一些不足，后来者提出了很多改进，在算法介绍一节里也有提及。

首先是实时性不高。这点如果大家实际运行了前面几幅船的图片的匹配应该会有体会。SURF据说要快一些。

然后是关键点主要依附于边缘拐角，因此对于边缘模糊和边缘光滑的图像，提取出的关键点就很少。

SIFT是在灰度图像上进行，不能充分利用图像的彩色信息，所以有CSIFT的扩展(Colored SIFT)。

还有抗仿射SIFT变换ASIFT(Affine-SIFT)。

备注：

本文原为课程作业，文字内容系本人参阅Lowe论文及网上资料后之个人理解，其中偏误，还请自甄，若能指教，不胜感激。

另外若本文有幸被转载，请保留出处：http://cql.cnblogs.com

附：本文对应的OPENSIFT源码及我加入VS项目后修改注释的版本下载https://files.cnblogs.com/cql/OPENSIFT.zip