sift算法理解
SIFT(Scale-Invariant Feature Transform)
尺度不变特征转换(Scale-invariant feature transform或SIFT)是一种电脑视觉的算法用来侦测与描述影像中的局部性特征，它在空间尺度中寻找极值点，并提取出其位置、尺度、旋转不变量，此算法由 David Lowe在1999年所发表，2004年完善总结。在CNN用来图像特征提取之前，SIFT算是手工设计特征提取算法的巅峰了。关于sift算法，有很多大佬讲的很好了，参考：
-
SIFT原始论文翻译：
https://www.cnblogs.com/cuteshongshong/archive/2012/05/25/2506374.html
-
大佬讲解文章：
https://blog.csdn.net/zddblog/article/details/7521424
https://blog.csdn.net/qq_37174526/article/details/100052271
https://www.cnblogs.com/Alliswell-WP/p/SIFT.html
-
源码学习：
vlFeat： vision Lab Feature library
https://github.com/vlfeat/vlfeat
https://www.vlfeat.org/overview/sift.html
opencv
python：
首先看一下sift特征提取和匹配算法的步骤：
- 构建尺度空间， 生成高斯差分金字塔（DOG金字塔）
- 空间极值点检测(关键点的初步查探)
- 稳健关键点的精确定位和筛选
- 精确定位：离散空间的局部极值，通过迭代法寻找连续空间极值
- 阈值筛选：消除小于阈值的极值
- 消除边缘效应：去除掉可能是图片卷积过程中边缘效应产生的极值点
- 稳定关键点的方向信息分配
- 关键点描述
- 特征点匹配
sift详细解释，大佬们已经讲的很好了，这里只记录下自己的一些困惑和理解，希望能帮助到和我一样的同学。(个人理解不一定完全对，欢迎纠错)
1. 高斯差分金字塔
sift为什么要建立高斯差分金字塔呢？弄清楚这个，得注意两个关键词：高斯， 差分。
高斯
熟悉图像处理的同学，都知道高斯卷积核能用来平滑图像，使图像变得模糊，如opencv中的GaussianBlur函数。而sift算法是为了找到尺度不变的特征点，就是图片中物体的特征点不随尺度变化，所以要得到物体不同尺度下的图片，对图片进行高斯卷积可以模拟物体在不同尺度下的图片(想象下，我们把图片缩小，图片中的物体会变小，同时也变得模糊)
差分
得到物体在不同尺度下的图片后，很明显下一步要做的就是特征提取，所谓图片特征就是特别突出的像素点，如比局部背景像素点都大的像素点，就是图像处理中常说的边缘，角点等。在图像处理中，laplacian算子能用来提取物体边缘，观察下laplacian卷积核，就能发现laplacian算子找到的这些边缘都是局部极值，如下是一个3x3的laplacian卷积核：(比周围四个像素点都大或都小才会被提取出来)
所以总结下，先对图片进行高斯卷积，再进行拉普拉斯卷积，就能得到不同尺度下物体的特征，这一步骤常称为高斯拉普拉斯(LoG, Laplacian of Gaussian)，但是这个计算比较耗时，大佬们发现了一种近似LoG的算法，但是计算更加简单，就是高斯差分(DoG, Difference of Gaussian)。 高斯差分DoG的步骤，就是采用不同sigma的高斯核分别对两幅图片进行卷积，然后求两幅图片相减，即差分。
采用LoG提取特征，LoG具有旋转不变性：https://www.cnblogs.com/Alliswell-WP/p/Laplace_Gaussian.html
关于Dog是LoG的近似，推导过程：https://www.cnblogs.com/silence-cho/p/13621975.html
2. 稳健关键点的精确定位和筛选
sift算法在找到尺度空间中的极值点后，还进行了一步精确定位和筛选，从而获得稳定的关键点，主要包括三步：
- 精确定位：离散空间的局部极值，通过迭代法寻找连续空间极值
- 阈值筛选：消除小于阈值的极值
- 消除边缘效应：去除掉可能是图片卷积过程中边缘效应产生的极值点
2.1 精确定位
在确定极值点时，是通过与其邻近的26个点作比较，如下图所示，中间尺度图片上的点，与同一尺度的8个点，上下两个尺度各9个点，共26个点进行大小比较，确定其是不是极值点。如果我们建立一个三维坐标系，包括x， y， s三个维度，分别表示图片的x轴方向，y轴方向，以及图片的尺度，对于三个维度下的每个整数坐标点与其邻域26个整数坐标点进行比较，很明显我们得到是一个离散空间极值点，如果想要得到连续空间的极值点，则需要对这个极值点的位置进行修正，即精确定位极值点。
如何确定连续空间的极值点呢？一提到极值点，熟悉深度学习的同学很快就会想到梯度下降法和牛顿迭代法，sift这里采用的就是类似的牛顿迭代法。牛顿迭代法的公式如下：
上述公式含义就是, 对于函数\(f(x)\)，在坐标点\(x^{(t)}\)，沿着\(- H_t^{-1}g_t\)方向移动，就能得到函数值比其更小的值，只要我们进行不断的迭代，坐标点就会收敛于极值点处。
所以回到我们的原始问题中，在图像上某一点处，我们是能计算该像素点处的一阶导数和二阶导数，也就能得到梯度\(g\)和海森矩阵\(H\), 那么在离散空间的极值点处，通过牛顿迭代法，就能得到连续空间的极值点。可以看下sift算法源码中是怎么计算的：
有x, y, s三个维度，所以梯度为三个维度的一阶偏导数，海森矩阵是三个维度的二阶偏导数， 即：
看下代码实现：
梯度计算
def computeGradientAtCenterPixel(pixel_array):
dx = 0.5 * (pixel_array[1, 1, 2] - pixel_array[1, 1, 0])
dy = 0.5 * (pixel_array[1, 2, 1] - pixel_array[1, 0, 1])
ds = 0.5 * (pixel_array[2, 1, 1] - pixel_array[0, 1, 1])
return array([dx, dy, ds])
函数说明：
"""Approximate gradient at center pixel [1, 1, 1] of 3x3x3 array using central difference formula of order O(h^2), where h is the step size
# With step size h, the central difference formula of order O(h^2) for f'(x) is (f(x + h) - f(x - h)) / (2 * h)
# Here h = 1, so the formula simplifies to f'(x) = (f(x + 1) - f(x - 1)) / 2
# NOTE: x corresponds to second array axis, y corresponds to first array axis, and s (scale) corresponds to third array axis
"""
海森矩阵计算
def computeHessianAtCenterPixel(pixel_array):
center_pixel_value = pixel_array[1, 1, 1]
dxx = pixel_array[1, 1, 2] - 2 * center_pixel_value + pixel_array[1, 1, 0]
dyy = pixel_array[1, 2, 1] - 2 * center_pixel_value + pixel_array[1, 0, 1]
dss = pixel_array[2, 1, 1] - 2 * center_pixel_value + pixel_array[0, 1, 1]
dxy = 0.25 * (pixel_array[1, 2, 2] - pixel_array[1, 2, 0] - pixel_array[1, 0, 2] + pixel_array[1, 0, 0])
dxs = 0.25 * (pixel_array[2, 1, 2] - pixel_array[2, 1, 0] - pixel_array[0, 1, 2] + pixel_array[0, 1, 0])
dys = 0.25 * (pixel_array[2, 2, 1] - pixel_array[2, 0, 1] - pixel_array[0, 2, 1] + pixel_array[0, 0, 1])
return array([[dxx, dxy, dxs],
[dxy, dyy, dys],
[dxs, dys, dss]])
函数说明：
"""Approximate Hessian at center pixel [1, 1, 1] of 3x3x3 array using central difference formula of order O(h^2), where h is the step size
"""
# With step size h, the central difference formula of order O(h^2) for f''(x) is (f(x + h) - 2 * f(x) + f(x - h)) / (h ^ 2)
# Here h = 1, so the formula simplifies to f''(x) = f(x + 1) - 2 * f(x) + f(x - 1)
# With step size h, the central difference formula of order O(h^2) for (d^2) f(x, y) / (dx dy) = (f(x + h, y + h) - f(x + h, y - h) - f(x - h, y + h) + f(x - h, y - h)) / (4 * h ^ 2)
# Here h = 1, so the formula simplifies to (d^2) f(x, y) / (dx dy) = (f(x + 1, y + 1) - f(x + 1, y - 1) - f(x - 1, y + 1) + f(x - 1, y - 1)) / 4
# NOTE: x corresponds to second array axis, y corresponds to first array axis, and s (scale) corresponds to third array axis
然后计算\(- H_t^{-1}g_t\)时，一般不直接计算H的逆矩阵，而是采用最小二乘法求解矩阵方程的方式，如下所示，矩阵方程\(Hx = g\)的解x正是我们需要的值，而矩阵方程的解可以采用最小二乘法求解。
计算处\(- H_t^{-1}g_t\)后，分别对x,y,s三个维度坐标值进行更新，sift算法中只对大于0.5的值进行更新，太小就不更新(小于0.5,说明此处就是极值点，不需要更新)，而且只进行了5次迭代,看下代码如下：
for attempt_index in range(num_attempts_until_convergence): # num_attempts_until_convergence=5，表示5次迭代
# need to convert from uint8 to float32 to compute derivatives and need to rescale pixel values to [0, 1] to apply Lowe's thresholds
first_image, second_image, third_image = dog_images_in_octave[image_index-1:image_index+2]
pixel_cube = stack([first_image[i-1:i+2, j-1:j+2],
second_image[i-1:i+2, j-1:j+2],
third_image[i-1:i+2, j-1:j+2]]).astype('float32') / 255.
gradient = computeGradientAtCenterPixel(pixel_cube) # 梯度值g
hessian = computeHessianAtCenterPixel(pixel_cube) # 海森矩阵H
extremum_update = -lstsq(hessian, gradient, rcond=None)[0] # 最小二乘法求解H^-1 *g
if abs(extremum_update[0]) < 0.5 and abs(extremum_update[1]) < 0.5 and abs(extremum_update[2]) < 0.5:
break
j += int(round(extremum_update[0])) #更新x
i += int(round(extremum_update[1])) #更新y
image_index += int(round(extremum_update[2])) #更新s
# make sure the new pixel_cube will lie entirely within the image
if i < image_border_width or i >= image_shape[0] - image_border_width or j < image_border_width or j >= image_shape[1] - image_border_width or image_index < 1 or image_index > num_intervals:
extremum_is_outside_image = True
break
迭代结束后，精确定位了极值点位置， 根据泰勒公式可以得到该位置对应极值点的近似值，泰勒公式的二次展开如下：
上面式子中，\(\vec X = (x, y, s)\)表示尺度空间中某一点，即从 \(\vec X_0\) 移动到 \(\vec X\) 处，根据我们迭代公式知道，移动值\(\vec X-\vec X_0\)为\(- H_t^{-1}g_t\)，带入上述公式得到:
对应代码中如下：
functionValueAtUpdatedExtremum = pixel_cube[1, 1, 1] + 0.5 * dot(gradient, extremum_update) # extremum_update即为-H^-1g
牛顿迭代法参考： https://zhuanlan.zhihu.com/p/37588590
https://www.cnblogs.com/happylion/p/4172632.html
2.2 阈值筛选
这一步比较简单，在上述精确定位极值点位置后，根据泰勒公式得到该位置的对应极值，sift中设置了一个阈值(contrast_threshold)，如果小于阈值就删除掉这个极值点。
2.3 消除边缘效应
边缘效应指对图片进行卷积时，对于图片上下左右边界处的像素点，常需要对图片补0，方便进行卷积运算，会导致边界处产生较强的效应，得到错误的极值点。对于图片边界处错误的极值点，其x方向和y方向上二阶导数，会出现一大一小。如下图中，在图片的上边边界，x方向上的二阶导数dxx比较小(两侧都位于图像中)，y方向上的二阶导数dyy较大(一侧位于图像中，对应值为像素值，另一侧位于图像外面，对应值为填充的0)
所以边缘效应产生的极值点，一个方向上二阶导数大，另一个方向上二阶导数小，那么只需过滤掉两者比值较大的极值点。上面计算的海森矩阵包括了x方向和y方向的二阶导数，
假设\(dxx=\alpha, dyy=\beta, \alpha = r \beta\), 则可以构建下面的等式， 很明显r=1时，这个表达式取得最小值，r大于1，或者小于1时，表达式值增大。所以只需要设置一个阈值，计算下面表达式的值，是否超过阈值即可。
代码中如下：
functionValueAtUpdatedExtremum = pixel_cube[1, 1, 1] + 0.5 * dot(gradient, extremum_update)
if abs(functionValueAtUpdatedExtremum) * num_intervals >= contrast_threshold:
xy_hessian = hessian[:2, :2] # 海森矩阵
xy_hessian_trace = trace(xy_hessian) # α+β
xy_hessian_det = det(xy_hessian) # αβ
if xy_hessian_det > 0 and eigenvalue_ratio * (xy_hessian_trace ** 2) < ((eigenvalue_ratio + 1) ** 2) * xy_hessian_det: # eigenvalue_ratio 即为r
# Contrast check passed -- construct and return OpenCV KeyPoint object
keypoint = KeyPoint()
keypoint.pt = ((j + extremum_update[0]) * (2 ** octave_index), (i + extremum_update[1]) * (2 ** octave_index))
keypoint.octave = octave_index + image_index * (2 ** 8) + int(round((extremum_update[2] + 0.5) * 255)) * (2 ** 16)
keypoint.size = sigma * (2 ** ((image_index + extremum_update[2]) / float32(num_intervals))) * (2 ** (octave_index + 1)) # octave_index + 1 because the input image was doubled
keypoint.response = abs(functionValueAtUpdatedExtremum)
return keypoint, image_index
边缘效应参考：https://www.cnblogs.com/pegasus/archive/2011/05/19/2051416.html
https://blog.csdn.net/qq_34886403/article/details/83589108
https://zhuanlan.zhihu.com/p/366778678
3. 关键点描述
在生成关键点的描述子(128维向量)时，有个地方也值得注意下，sift算法采用了三线性插值(逆运算)。
如下图所示，这里也有三个维度：x，y，angle，分别表示关键点的x值，y值，幅角，而对应的值f(x, y, angle)表示幅值。下面图中红色点是我们计算得到一个坐标点(x, y, angle), 这个坐标可能是小数，如(1.5, 1.5 1.5), 我们需要将其幅值分配到与其相邻的八个整数坐标点上, 就需要采用三次插值逆运算，即分别在x维度，y维度，angle维度进行幅值分配。
(注意：sift中将360度分成了8个区间，每个区间45度，所以angle维度中：1表示45度，2表示90度，。。。7表示360度)
最后，sift中还有很多细节和参数值计算方法，多看几遍代码和讲解文章应该就清楚了。
