SIFT代码V1版本（Matlab）实现以及思考的问题

    写在前面的话：根据网络上的资料整理学习了一下SIFT的原理，并且完成了代码实现。除了高斯滤波矩阵求逆这样的自带函数以外几乎所有的函数都是自己编写实现的。这里要非常感谢网络上的各位提供的资料，尤其是：

https://blog.csdn.net/zhaocj/article/details/42124473，这是他的原话“可以在线阅读，也可以免费下载（在这里，鄙视那些设置下载权限和积分的人！！！）。如果上述方法都不方便，可以留下email，向我索要。文章中错误的地方欢迎指正！”我把他的源码和注释理解消化以后自己编写了matlab版本的代码。我会附上我自己的代码，但是CSDN下载需要积分（实在是我个人积分太少，鄙视一下自己）。如果您积分不够可以私信我邮箱494431025@qq.com，我会免费提供的（有积分您还是送我一点吧，不麻烦邮件传递了）。分享的目的在于能够让大家共同学习，同时非常欢迎大家的指正。

1.什么是sift

    网络上的sift过程讲解一大把，估计很多人和我一样，一眼望去懵逼了。刚开始接触sift完全不明白这整个过程在干什么，看的人晕晕乎乎的。

    我个人的理解：SIFT就是找出图片中的二阶导数的极值点，然后任意一个极值点的周围邻域内的像素求描述符。该描述符就是邻域内像素点的一阶导数的某种组合。

2.sift过程

    正如前面所说，网络上的求解过程一大把。我这里就不啰嗦了，再说以我的语文水平也写不了别人那么好。偷个懒摘录人家的：

https://blog.csdn.net/ijuliet/article/details/4640624

    下面，先整理SIFT思路：
1. 输入图像，建议做double(width*=2, height*=2, size*=4), 并高斯过滤进行平滑。
2. 由图片size决定建几个塔，每塔几层图像(一般3-5层)。0塔的第0层是原始图像(或你double后的图像)，往上每一层是对其下一层进行Laplacian变换（高斯卷积，其中sigma值渐大，例如可以是sigma,k*sigma, k*k*sigma…），直观上看来越往上图片越模糊。塔间的图片是降采样关系，例如1塔的第0层可以由0塔的第3层down sample得到，然后进行与0塔类似的高斯卷积操作。
3.构建DoG金字塔。DoG金字塔由上一步生成的Gauss金字塔计算得到，塔数相同，每塔层数少1，因为DoG的每一层由Gauss的相邻两层相减得到。
4.在DoG塔里进行极值点检测，并根据用户预设的对比度阈值、主曲率阈值去除不合法特征点。极值点检测用的Non-Maximal Suppression，即在3*3*3个点中进行灰度值比较，最小或最大才过关。
5. 计算每个特征点的尺度。注意塔间尺度关系，sigma*2.0^(octvs+intvl/intvls)
6.计算每个特征点的梯度模值和方向。用特征点周围一个矩阵区域(patch)内的点来描述该特征点，用的直方图进行模值统计并寻找主方向，主方向可以不止一个。

7. 最后要生成64D或128D的特征描述符了。对齐主方向，计算方向直方图2D数组，假如每个直方图有8bin，那么64D(2*2*8bin)或128D(4*4*8bin)。

    如果你理解了第一部分的那句话，你就会发现整个过程都是围绕着这句话在做文章。

    在这里附上我自己程序中的程序框架：结合上述过程能够更好的理解

    本程序主要函数分为以下几大块
1.主函数match_yhb.m,读入两张图片并输出匹配的结果
2.函数 SIFT_yhb.m为生成SIFT的特征点以及描述符
3.函数Define（）被SIFT_yhb.m调用，主要是生成全局变量。有关的参数定义全部在这里
4.函数sigma_num（）被Define调用，生成sigma数组。均是相对于本组内的尺度，并作为全局变量。
5.buildgausspyr()被SIFT_yhb.m调用生成高斯金字塔
6.buildDogpyr(）被SIFT_yhb.m调用生成Dog金字塔
7.findExtrma(）被SIFT_yhb.m调用，初步寻找极值点
8.adjustExtrPoint(）被findExtrma(）调用。精确点的进一步判断选择，并得出精确位置
9.oritHist(）被findExtrma(）调用。返回极值点的第一主方向，在 SIFT_yhb.m中进一步得出多个主方向

10.calcDescriptors(）被SIFT_yhb.m调用，生成128维描述符。

3.sift思考

    当你隐隐约约看明白了这个过程的时候，你肯定有很多的疑问。至少我是有很多疑问的。在这里我将整个过程中思考的一些问题总结一下，我估计这里面会有很多不准确的地方，非常欢迎批评指正（但不是说脏话，艹效果糟糕嘛，辣鸡），如果能附上您个人的见解那是最美不过。

3.1为什么要分组，分了组以后为啥还要分层？

    直观的理解，拍摄的第1张图片是近距离拍摄的，第2张图片是远距离拍摄的。当然分组（相当于图片缩放）以后两张图片能够对物体本身能够有更好的匹配。

    分层就是利用不同的尺度模糊，这里贴上网络上的原话：尺度空间理论的基本思想是：在图像信息处理模型中引入一个被视为尺度的参数，通过连续变化尺度参数获得多尺度下的尺度空间表示序列，对这些序列进行尺度空间主轮廓的提取，并以该主轮廓作为一种特征向量，实现边缘、角点检测和不同分辨率上的特征提取等。

    大概意思是图片越模糊主要就显示主轮廓。通过不同的尺度下寻找二阶极值点能够模拟灵长类动物的眼睛实际上看物体，尺度越大图像越模糊。我没法猜测LOWE当时提出尺度连续的意义，我个人认为意义不大，为了求得Dog塔极值点构造必要的4层差不多就够了（Dog上下需要一层做比较，即需要3层，高斯比Dog多1层即4层）。实际上程序运行时，在大尺度的图像中寻找到的极值点也确实明显少很多。当然因为Dog是Log的近似，尺寸连续的高斯4层还是必要的。

3.2 如何理解尺度，尺寸，以及为啥要连续的尺度？

    首先一定要明白这一点尺度相对于尺寸（图片大小，像素点的数量）才有意义。比如一张图片扩大两倍以后用3.2的尺度去高斯模糊，和原图片用1.6的尺度模糊，是相同的效果。我还是尽量少帖公式，这样子理解吧，显然当图片某点距离原点距离n，放大k倍以后的图片距离为k*n。正态分布公式中，距离的平方除以方差（这里指的是尺度）的平方，因此放大以后图片为了相似的权重效果也需要放大尺度为k倍。到这里你会发现如果所有分组分层的图片放大到同一尺寸的图片后尺度还真正是连续的。即：

因此组内用模糊即可，但是这两个公式其实是等价的。一个是相当于在组内的尺寸尺度，一个是同一个尺寸下的尺度。理解了尺度与尺寸的关系就明白了倒数第3层降采样的尺度关系，我这里不贴公式啰嗦了。

         上文也讲了，我个人认为连续尺度没有多大必要（当然你认为尺度大的模糊图能够定位到大的轮廓那就相当有必要了），但是还是必须做。那凭啥连续的尺度是等比数列还不是等差数列呢？这里不得不搬出：

         这个公式没有仔细研究，但是不妨碍高尺度下(大致)理解这个公式。意思是两个不同的高斯层相减得到Log的近似。很明显两个高斯是一个等比为k的尺度关系。其实这个公式是值得深挖的，因为涉及到后面极值点寻找的问题。不同的Dog层找极值点，难道你就不担心不同的层之间其实不是二阶极值点Log的等比例近似。这么说吧，比如上面一层是Log的8倍近似，下面一层是Log的4倍近似，这样当下面层内的极值点与上面一层做极大极小比较时，下面的值被上面的放大的极值点给Pass掉了。希望大家提意见，如果仅仅是Log的（k-1）倍近似，那当然没问题，大家都是相同的等比k，这样Dog图中不用担心相邻层被不同的系数放大（也就是大家都是Log的（k-1）倍近似）。

结论：尺度等比例连续本身就有必要，构造Log的等比例近似。如果你觉得尺度大的图里面有检测到大轮廓，那就多构造几层，至少Lowe建议是这样的，构造s=3层。

3.3 冥冥之中偏偏选中了用二阶极值点作为判断依据？

         2002年Mikolajczyk在详细的实验比较中发现尺度归一化的高斯拉普拉斯函数的极大值和极小值同其它的特征提取函数（如：梯度，Hessian或Harris角特征）比较，能够产生最稳定的图像特征。

         这一段话是摘录的，如果有人研究了该论文请给出建设性的意见。至少我知道比如一个白色矩形，一阶导数的极大极小值是边界。相反角点处的一阶导数并不大，但是角点的二阶导数比边界大。而sift就是想定位比如角点、边缘点、暗区域的亮点以及亮区域的暗点，既然两幅图像中有相同的景物，那么使用某种方法分别提取各自的稳定点。我当时是不理解的，不过做了个实验看看这些点的一阶导数二阶导数慢慢的就明白了。实际上定位边界不是没有用而是错误率很高，人类的眼睛是对边界敏感，但是机器做匹配不好做。后面会接着讲这个问题。

3.4 不稳定特征点判据如何归一化？

    LOWE指出即使是极值点但是其值本身比较小或者精确定位泰勒拟合以后还是比较小（阈值0.04，相较于1.0），要去除。但是归一化是个麻烦事，我们也知道程序中很多时候都是做了等比例的放大或者缩小的。找极值点没问题，但是极值点排除就出问题了。因为极值点排除是一个相对于1.0的绝对阈值0.04，比如Dog中最小是-0.2最大是0.1（为啥给出这样的-2倍关系是有依据的，不然你自己在图上画一个阶跃信号，二阶导数的最小值为最大值的-2倍）。你该如何归一化。归一到（0-1）显然不对，最小的二阶导数归一到0附近反而被0.04干掉了。归一到（-1~1）之间，貌似也不对，因为原来的零点被抬高了，新的零点也不是实际上二阶导数的零点。我自己的思路是不归一化，直接用层内极小值极大值的绝对值的最大者乘以0.04作为阈值，这样子将Dog中零点附近的不稳定点排除。

3.5 Hessian矩阵去除边界点？

    Lowe说DOG算子会产生较强的边缘响应，需要剔除不稳定的边缘响应点。请问这句话你能看出啥意思不？至少我是完全不知道为啥好不容易检测出了极值点而且处于人类敏感的边界位置处，好端端的还是要被干掉。要不大神怎么是大神呢。不相信人家是可以的，但是吃过苦头你就明白了。我现在的程序就有错误匹配，这些错误点往往就是在边界位置处。再来看这句话似乎理解了一点点，边界位置很容易产生极值点，而边界是很不好判断的。其实打开程序运行结果，看到错误点你就明白了，可能这个边界的点会被误定位到另外一个边界上的点。即使这两个边界方向不同（因为sift的特征描述符对主方向旋转，所以特征描述符是具有旋转不变性质的）。

    这也是我下一步可能思考的问题，如何降低错误匹配率，暂时想到的是既然我已经通过hessian矩阵删除掉边界特征点（其实是做不到完全删除边界周围的特征点的）还是有错误匹配，我打算考虑主方向。如果主方向五花八门那就是最好。如果我继续优化程序降低了错误率我会贴上结果以及总结。

    还有一点，OpenCV中Hessian矩阵是在Dog中求，我是在高斯中求，因为我认为既然判断边界当然应该在源图中找，而不是Dog中找，而且Dog已经是二阶导数，继续Hessian矩阵相当于4阶导数。这样的高阶导数对数据是很敏感的。

4百思不得解的小问题（第3部分的问题都能理解这里只是小菜）

4.1找极值点是连续的上下左右找，不同组之间怎么办？

    我刚看这个过程的时候也有这些疑问，不过你看了代码就不会有疑惑了，答案是当然不会在不同组之间去找极值点。

4.2主方向要干吗？

    主方向就是为了匹配旋转不变性，特征点描述符都是以这个主方向为基准方向做描述的，而且有可能不止一个主方向喔，不相信你就用一个白色矩形做实验，而且不同主方向正如你所想的相差90°。在这里我也提一下程序中，主方向统计完36个柱子后，平滑以后，求主方向的确切值不是采用参考代码中的抛物线插值，而是用左右两个柱子和本身柱子做重心得到。还是采用的白色矩形做实验，发现重心法得到的结果比抛物线插值更接近实际主方向。

4.3特征点描述符的三线性插值？

    在本文一开头的那篇参考文献中，我一开始真没明白到底要干啥。我在这里插一下我个人的理解。相当于一个点的梯度，对它领域的不同方块d，和不同的方向柱n的贡献，一共被八个顶点加权。那么该点被切了三刀，共产生8个立方体，对某一个顶点的贡献就是与顶点最远的那个立方体的体积。

4.4 领域的范围？

    这个问题是我到目前都没有很好的理解的，或许这个问题本身就不大，比如3σ，这个σ到底是指的组内的尺度，还是相较于基础组的尺度，或者说相较于上一层高斯滤波采用的尺度，个人认为最合适的应该还是组内的尺度。不过我程序中好像不是这么写的。

    最后贴上效果图，匹配点不够多，有一定的错误率。LOWE主页的程序同样也有错误率。或许SIFT的错误率是一个不可避免的话题。请各位大神给出建设性的意见。

代码下载链接：https://download.csdn.net/download/csuyhb/10330738