SIFT算法

SIFT（Scale Invariant Feature Transform），尺度不变特征变换。具有旋转不变性、尺度不变性、亮度变化保持不变性，是一种非常稳定的局部特征。在目标检测和特征提取方向占据着重要的地位。

SIFT算法所查找到的关键点是一些很突出，不因光照、仿射变换和噪声等因素而变化的点，比如角点、边缘点、暗区亮点或亮区暗点等。

SIFT算法步骤

1. 尺度空间极值检测：在各个尺度的图像上搜索极值位置。通过高斯差分函数识别潜在的对于尺度和旋转不变的关键点。

2. 确定关键点位置：在每个潜在的位置上，通过一个拟合的精细模型确定位置和尺度。关键点的选择依据于其稳定程度。

3. 确定关键点方向：基于局部图像的梯度方向，给每个关键点位置分配了一个或多个方向。所有后面对图像数据的操作都相对于关键点的方向、尺度和位置进行变换，从而保证此类变换的不变性。

4. 关键点描述：在每个关键点周围邻域内，在给定尺度上测量图像局部的梯度，这些梯度作为关键点的描述符，允许比较大的局部形状的变化或光照变化。

尺度空间极值检测

当我们近距离和远距离看一个物体时，模糊程度是不一样的，从近距离到远距离图像细节越来越模糊。在SIFT中，图像的尺度的意思也是类似的，注意这里的尺度空间指的不是图像的大小，而是图像的模糊程度。

从频域的角度看，图像的模糊程度代表的频域信息的低频成分和高频成分。模糊图像代表信息大部分都集中低频，仅含少量的高频信息。图像清晰代表信息成分丰富，高低频的信息都有。

高斯核是唯一可以产生多尺度空间的核函数。这一点在T Lindeber的文献《Scale-space theory: a basic tool for analyzing structures at different scales》中被证明。

图像的尺度空间 $L(x,y,\sigma)$ 定义为原始图像(x,y)与一个可变尺度的二维高斯函数 $G(x,y,\sigma)$ 的卷积运算，即：

$L(x,y,\sigma) = G(x,y,\sigma) * I(x,y)$

其中：

$G(x,y,\sigma) = \frac{1}{2\pi\sigma^2}e^{-(\frac{x^2 + y ^2}{2\sigma^2})}$

$\sigma$ 是尺度空间因子，它决定了图像的模糊程度。当 $\sigma$ 较大时，尺度较大，图像模糊，显示的是概况为主。 $\sigma$ 较小时，尺度较小，显示的是更多细节。

到这里，相信大家对尺度空间的概念应该有了一定认识。如果还不清楚，我这里举一个不太恰当的例子：

假设有两个人A和B站在同一窗户旁边看向远处的一颗树。A是马赛人，视力极其牛逼；B是一个300度的近视眼。A和B所看到的树的图像就属于不同尺度空间。 $\sigma$ 在这里就可以理解为近视度数。

在计算高斯函数离散近似时，距离 $3\sigma$ 之外的像素所起的作用几乎为0，因此可以忽略掉这些像素，一般在实际应用时，高斯卷积核只要计算 $(6\sigma+1) *(6\sigma+1)$ 就可以了。

在SIFT做尺度空间极值检测时，需要构建高斯金字塔，高斯金字塔的细节，请在“参考”中的链接"高斯金字塔"中去查看。得到高斯金字塔后，我们就得到了不同大小的图片组（octave），每个octave由多个不同尺度的图像组成。

原图地址：高斯金字塔_小鹿焕焕的博客-CSDN博客

要找到尺度空间下的极值，可以使用LoG（高斯拉普拉斯方法）来求得，这在边缘检测笔记中有提到，其本质就是使用图像的二阶导数来计算。但LoG的计算量较大，在SIFT算法原始论文中，提到了使用高斯擦还分函数来做归一化拉普拉斯算子的近似，通过两个相邻高斯尺度空间图像直接相减即可：

$D(x, y, \sigma)=(G(x, y, k\sigma)-G(x, y,\sigma))* I(x, y)$

$= L(x, y, k\sigma) - L(x, y, \sigma)$

D(DoG)函数和LoG之间的关系：

具体证明可以参考下面的链接：

DoG和LoG算子 - 知乎DoG(Difference of Gaussian)算子和LoG(Laplacian of Gaussian)算子是常用的极值点检测(Blob Detection)两种方法，高斯卷积是为了进行尺度变换，那么LapLacian呢。因此这里首先引入LapLacian算子。图像边缘检测…https://zhuanlan.zhihu.com/p/49447503 通过DoG算子，我们可以得到高斯差分金字塔，具体计算很简单，对每一个octave中相邻两个尺度的图像做减法即可：

可以看到，DoG金字塔每个Octave的层数是原始高斯金字塔层数 - 1。

接下来，通过DoG金字塔来查找不同尺度空间中的局部最大值或最小值。如下图：

对于某个像素点（图中以X标记）来说，它需要与自己相邻的8个邻域，以及尺度空间中相邻的上下两层的18个临近点（2x9）做比较。如果该点的值在所有这些点中是最大或最小的，则该点可能是一个关键点，这个点会作为候选点。

局部最大值的搜索过程从每组的第二层开始，以第二层为当前层，对当前层的DoG图像中的每一个点都和26个（当前层8个+上下临近层各9个）临近点进行比较。搜索得到的极值点包含的信息既有位置坐标（DoG图像坐标），也有尺度空间信息（第几层）。当第二层搜索完成后，以第三层为当前层重复局部最大值搜索过程。当S = 3时，表示每组里面要搜索3层，因此DoG里需要有S+2层，最初构建的高斯金字塔就需要有S+3层。

确定关键点位置（精确定位关键点）

前一步我们找到的候选极值点，是在离散空间上的候选极值点，这样检测出的极值点，和真实极值点位置有一定差异，为了更精确地定位到真正的极值点，SIFT算法中使用了尺度空间函数 $D(x,y,\sigma)$ 的三元二阶泰勒展开来定位。

假设我们检测到的极值点为 $X_0 = (x_0,y_0,\sigma_0)$ ，则该点对应的泰勒展开为：

$f(\begin{bmatrix} x\\ y\\ \sigma \end{bmatrix}) = f(\begin{bmatrix} x_0\\ y_0\\ \sigma_0 \end{bmatrix}) + \begin{bmatrix} \frac{\partial f}{\partial x}, & \frac{\partial f}{\partial y} ,& \frac{\partial f}{\partial \sigma} \end{bmatrix}\begin{bmatrix} \begin{bmatrix} x\\ y\\ \sigma \end{bmatrix} - \begin{bmatrix} x_0\\ y_0\\ \sigma_0 \end{bmatrix} \end{bmatrix}$

$+ \frac{1}{2}\begin{bmatrix} \begin{bmatrix} x\\ y\\ \sigma \end{bmatrix} - \begin{bmatrix} x_0\\ y_0\\ \sigma_0 \end{bmatrix} \end{bmatrix}^T\begin{bmatrix} \frac{\partial ^2f}{\partial x \partial x}, & \frac{\partial ^2f}{\partial x \partial y}, & \frac{\partial ^2f}{\partial x \partial \sigma} \\ \frac{\partial ^2f}{\partial x \partial y}, & \frac{\partial ^2f}{\partial y \partial y}, & \frac{\partial ^2f}{\partial y \partial \sigma}\\ \frac{\partial ^2f}{\partial x \partial \sigma}, &\frac{\partial ^2f}{\partial x \partial \sigma} & \frac{\partial ^2f}{\partial \sigma \partial \sigma} \end{bmatrix}\begin{bmatrix} \begin{bmatrix} x\\ y\\ \sigma \end{bmatrix} - \begin{bmatrix} x_0\\ y_0\\ \sigma_0 \end{bmatrix} \end{bmatrix}$

我们用向量 $\widehat{X}$ 替换 $\begin{bmatrix} x\\ y\\ \sigma \end{bmatrix} - \begin{bmatrix} x_0\\ y_0\\ \sigma_0 \end{bmatrix}$ ，它表示曲线上某个点相对于候选极值点 $X_0$ 的一个偏移量，用 $X$ 表示向量 $\begin{bmatrix} x\\ y\\ \sigma \end{bmatrix}$ 。则会得到：

$f(X) = f(X_0) + \frac{\partial f}{\partial X}^T\widehat{X} + \frac{1}{2}\widehat{X}^T\frac{\partial^2f}{\partial X^2}\widehat{X}$

这里的 $f(X)$ 就是我们说的尺度空间函数 $D(X)$ 的近似多项式形式,因此我们可以简单认为 $f(X)$ 和 $D(X)$ 是一样的，后续推导会继续使用 $f(X)$ 这个符号。

我们知道，二阶泰勒展开的目的是为了得到一个多项式来近似出原函数。得到这个展开式后，要求真实的极值点，其实就是求这个函数一阶导数为0的位置。

首先，求出这个函数的一阶导数：

$f(X_0)$ 是常数项，求导后为0。

$\frac{\partial f}{\partial X}^T\widehat{X}$ 包含X的一次方表达，求导后只剩下左边的常数项 $\frac{\partial f}{\partial X}^T$ 。

最右边的 $\frac{1}{2}\widehat{X}^T\frac{\partial^2f}{\partial X^2}\widehat{X}$ 这种形式的矩阵是一个二次型，对其求导可得 $\frac{1}{2}(\frac{\partial^2 f}{\partial X^2} + \frac{\partial ^2f}{\partial X}^T)\widehat{X}$

由于二次型矩阵的系数矩阵的转置和其本身相等，因此 $\frac{1}{2}(\frac{\partial^2 f}{\partial X^2} + \frac{\partial ^2f}{\partial X}^T)\widehat{X} =\frac{\partial^2 f}{\partial X^2} \widehat{X}$

综上，可得函数的一阶导为：

$\frac{\partial f(X)}{\partial X} =\frac{\partial f}{\partial X}^T +\frac{\partial^2 f}{\partial X^2} \widehat{X}$

让导数为零，求解：

$\frac{\partial f}{\partial X}^T +\frac{\partial^2 f}{\partial X^2} \widehat{X} = 0 \Rightarrow \widehat{X} = -\frac{\partial^2 f}{\partial X^2} ^{-1}\frac{\partial f}{\partial X}$

将 $\widehat{X}$ 代回到泰勒展开式中，有:

$f(X) = f(X_0) + \frac{\partial f}{\partial X}^T\widehat{X} + \frac{1}{2}(-\frac{\partial^2f}{\partial X^2}^{-1}\frac{\partial f}{\partial X})^T\frac{\partial^2f}{\partial X^2}(-\frac{\partial^2f}{\partial X^2}^{-1}\frac{\partial f}{\partial X})$

$= f(X_0) + \frac{\partial f}{\partial X}^T\widehat{X} + \frac{1}{2}\frac{\partial f}{\partial X}^T\frac{\partial^2f}{\partial X^2}^{-T}\frac{\partial^2f}{\partial X^2}(-\frac{\partial^2f}{\partial X^2}^{-1}\frac{\partial f}{\partial X})$

$= f(X_0) + \frac{\partial f}{\partial X}^T\widehat{X} + \frac{1}{2}\frac{\partial f}{\partial X}^T(\frac{\partial^2f}{\partial X^2}^{-1}\frac{\partial f}{\partial X})$

$= f(X_0) + \frac{\partial f}{\partial X}^T\widehat{X} + \frac{1}{2}\frac{\partial f}{\partial X}^T(-\widehat{X})$

$= f(X_0) + \frac{1}{2}\frac{\partial f}{\partial X}^T(\widehat{X})$

定位极值点的过程，本质是通过拟合出的曲线来求极值点位置并不断迭代：

迭代的意思是找到一个更精确的极值点位置后，再以它为基础，去拟合出一条曲线定位这条曲线上的极值点。因此需要注意两个问题：迭代次数要加以限制（多次计算后都没有收敛，三个分量都没有小于0.5）；如果求出的极值点位置超出了一定范围（比如0.5表示点的位置太远，因为泰勒展开式在特定点位置附近能较好地拟合出曲线，太远了就不准确了）。

舍去低对比度的点：

前面通过泰勒展开得到了比较粗的关键点位置，此时的关键点可能会包含一些低对比度的像素噪声。通过对比极值点的灰度值和一个阈值(一般是0.03左右，论文里假定像素值的范围是[0,1]）来决定是否要舍去极值点。小于阈值算法认为该点是一个噪声，会舍去。

去除边缘响应：

DoG在边缘上会有很强的响应，沿着边缘上的点如果不做处理也会被判断为关键点。为了去除边缘响应，需要引入黑森（海森）矩阵来做处理。

$H(x,y) = \begin{bmatrix} D_{xx}(x,y) & D_{xy}(x,y) \\ D_{xy}(x,y)& D_{yy}(x,y) \end{bmatrix}$

引入黑森矩阵H的目的是因为H描述了曲线的局部曲率，黑森矩阵的特征值标定了函数在相应特征向量方向上变化的快慢程度。如果函数在该点的不同方向上的曲率变化差不多，则其为角点。黑森矩阵的特征值和局部曲率是成正比的 (可以类比一下Harris角点检测算法，如果检测出边，则其M矩阵的特征值中是某一个大，某一个小；如果检测到角点，则两个特征值都大)。

SIFT算法中，也没有直接去求H矩阵的特征值，而是通过矩阵的迹和行列式的比值来间接得知两者的关系:

$Trace(H) = D_{xx} + D_{yy} = \alpha + \beta$

$Det(H) = D_{xx}D_{yy} - (D_{xy})^2 = \alpha\beta$

其中 $\alpha$ 和 $\beta$ 为两个特征值，并且 $\alpha > \beta, \alpha = r\beta$ 。表示 $\alpha$ 为两个特征值的最大的一个，且为 $\beta$ 的r倍。

当Det(H)小于0时，去掉点X，因此此时两个特征值异号，变化不均匀，为边缘。

当 $\alpha$ 和 $\beta$ 同号并且 $\alpha > \beta$ ，那么r大于1，通过迹的平方和行列式比值来做检测：

$\frac{Trace(H)^2}{Det(H)} = \frac{(\alpha + \beta)^2}{\alpha\beta} = \frac{(r+1)^2\beta^2}{r\beta^2} = \frac{(r+1)^2}{r}$

从这个表达式可以看出，只有两个特征值几乎相等时，这个表达式值最小，表达式的值随着r变大而变大。函数 $\frac{(r+1)^2}{r} \ \ \{r \geqslant 1\}$ 是一个单调递增的函数，最小值为4：

因此，在SIFT算法中，会设定一个阈值r（论文中取的是10），然后判断下面这个不等式是否成立：

$\frac{Trace(H)^2}{Det(H)} < \frac{(r+1)^2}{r}$

如果不满足此条件，则表示该点不符合要求，会被丢弃。

关键点检测过程概览总结

(a) 原始233x189图像；

(b) DoG函数极值检测出的最初832个关键点位置；

(d) 去除边缘响应后，最终剩下536个关键点。

确定关键点的方向

为了让关键点具有旋转不变性，SIFT算法会将特征点区域统一旋转到特定方向，这一方向即为特征点区域的主方向。如何确定这个方向呢？

对于每个关键点，SIFT统计以该点所在的高斯尺度图像中，以该点为圆心、半径为3*1.5 $\sigma$ （此半径并未从我看到的论文中查到，参考的是网络中的资料，请注意）的圆内，所有的像素的梯度方向及其梯度幅值。每个点的梯度和幅值计算方法如下：

从前面精确定位关键点的小节我们知道，关键点的尺度（亚像素级）大概率和原始高斯金字塔中的特定 $\sigma$ 尺度不是同一个值。这里说所的高斯图像的尺度，是通过找到和亚像素级的关键点尺度值最接近的那一层高斯图像。（注意：此时计算梯度和幅值是在高斯金字塔中的图像中进行，并非DoG金字塔中的图像）。

将计算得到的信息用直方图表示（论文中以每10度为一个范围，也就是有36个方向，下面的图为示意图，没有画出36个方向）：

直方图的每一个bin（柱）为10度，在圆内梯度方向在这个角度的像素都找出来，然后将他们的幅值相加作为这个bin（柱）的高度。需要注意的是，离圆心越远的像素对幅值的贡献程度就越小，实际计算幅值相加时，并非简单地加和，而是会通过高斯函数加权而来，高斯函数的方差为1.5 $\sigma$ 。

直方图峰值所对应的方向为主方向，任何大于峰值80%的方向为特征点的辅方向。这里得到的方向角度值是一个范围，要更精确地知道具体的角度，可以通过拟合直方图中主方向和辅方向的抛物线来获取。

现在，我们就得到了关键点的三个关键信息：位置、尺度和方向 $(x,y,\sigma,\theta)$ 。

关键点特征描述符

接下来，我们要对关键点建立特征描述符，描述符不仅包含了关键点的信息，也包含了关键点周围对其有贡献的像素点信息。主要思路是将关键点周围的像素区域进行分块，计算块内的梯度直方图，然后得到特征向量。

描述符和特征点所在的尺度是相关的，因此描述符是在关键点对应的高斯尺度图像上进行生成的。以关键点为中心，将其附近邻域划分为d*d个子区域（d一般为4）,每个区域是一个正方形。统计每个正方形区域内的梯度信息。对于要统计的关键点的区域，为了保证旋转不变性，需要先将图像的x轴旋转到关键点的主方向上。

然后，在每个子块内，统计像素在8个方向上的梯度值，统计梯度值的时候一般会进行高斯加权。

d一般为4，表示有16个子块，每个子块会得到8个值（8个方向），将这8个梯度信息记录下来就得到了关键点的特征描述符。因此，一个特征描述符的维度是：4*4*8 = 128。

参考资料

高斯金字塔

高斯金字塔_小鹿焕焕的博客-CSDN博客1、为什么要构建高斯金字塔高斯金字塔模仿的是图像的不同的尺度，尺度应该怎样理解？对于一副图像，你近距离观察图像，与你在一米之外观察，看到的图像效果是不同的，前者比较清晰，后者比较模糊，前者比较大，后者比较小，通过前者能看到图像的一些细节信息，通过后者能看到图像的一些轮廓的信息，这就是图像的尺度，图像的尺度是自然存在的，并不是人为创造的。好了，到这里我们明白了，其实以前对一幅图像的处理还是比较单调的，因为我们的关注点只落在二维空间，并没有考虑到“图像的纵深”这样一个概念，如果将这些内容考虑进去我们是不是会_高斯金字塔https://blog.csdn.net/qq_44771627/article/details/123840316 SIFT算法原始论文

https://www.cs.ubc.ca/~lowe/papers/ijcv04.pdfhttps://www.cs.ubc.ca/~lowe/papers/ijcv04.pdf

图像算法：Difference of Gaussian(DOG) 高斯函数差分_高斯差分_SSS_369的博客-CSDN博客概念Difference of Gaussian(DOG)是高斯函数的差分。它是可以通过将图像与高斯函数进行卷积得到一幅图像的低通滤波结果，即去噪过程，这里的Gaussian和高斯低通滤波器的高斯一样，是一个函数，即为正态分布函数。同时，它对高斯拉普拉斯LoG的近似，在某一尺度上的特征检测可以通过对两个相邻高斯尺度空间的图像相减，得到DoG的响应值图像。基本理论首先，高斯函数..._高斯差分https://blog.csdn.net/sss_369/article/details/84674639SIFT算法详解_夜半罟霖的博客-CSDN博客大纲引言一、高斯金字塔二、高斯差分金字塔三、特征点处理1.阈值化2.非极大值抑制3. 二阶泰勒修正待续引言 SIFT算法是为了解决图片的匹配问题，想要从图像中提取一种对图像的大小和旋转变化保持鲁棒的特征，从而实现匹配。这一算法的灵感也十分的直观：人眼观测两张图片是否匹配时会注意到其中的典型区域（特征点部分），如果我们能够实现这一特征点区域提取过程，再对所提取到的区域进行描述就可以实现特征匹配了。于是问题就演变成了以下几个子问题：应该选取什么样的点作为特征点呢？：人眼对图像中的高频区域更加的敏感，由此_sift算法https://blog.csdn.net/Dr_maker/article/details/121442210

关于二次型矩阵和其求导的相关内容，可参考这里：

二次型 - AcWing高数https://www.acwing.com/blog/content/35231/

二次型求导的推导过程? - 知乎直接给你看我们上课的PPT吧下图中的a对应xT，b对应Ax，过程就是直接套下面这个公式，最后一步，求导要…https://www.zhihu.com/question/22455493/answer/131136246?utm_id=0

极值点的精确定位_哔哩哔哩_bilibili极值点的精确定位是6.SIFT(尺度不变特征变换)的第2集视频，该合集共计5集，视频收藏或关注UP主，及时了解更多相关视频内容。https://www.bilibili.com/video/BV1Qb411W7cK/?p=2&vd_source=474bff49614e62744eb84e9f8340d91a

多元函数的泰勒展开

多元函数的泰勒展开 - 知乎类比一元函数的泰勒公式 f(x) = f(x_k)+(x-x_k)f'(x_k)+\frac{1}{2!}(x-x_k)^2f''(x_k)+\cdot\cdot\cdot+\frac{1}{n!}(x-x_k)^nf^{(n)}(x_k)+o^n 根据多元函数求偏导的法则，我们可以推出多元函数的泰…https://zhuanlan.zhihu.com/p/490766872

黑森矩阵

https://www.cnblogs.com/zywnnblog/p/16605038.htmlhttps://www.cnblogs.com/zywnnblog/p/16605038.html 有限差分方法