SIFT(Scale-Invariant Feature Transform)

尺度不变特征转换(Scale-invariant feature transform或SIFT)是一种电脑视觉的算法用来侦测与描述影像中的局部性特征，它在空间尺度中寻找极值点，并提取出其位置、尺度、旋转不变量，此算法由 David Lowe在1999年所发表，2004年完善总结。在CNN用来图像特征提取之前，SIFT算是手工设计特征提取算法的巅峰了。关于sift算法，有很多大佬讲的很好了，参考：

1. SIFT原始论文翻译：

   https://www.cnblogs.com/cuteshongshong/archive/2012/05/25/2506374.html

2. 大佬讲解文章：

   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

3. 源码学习：

   vlFeat： vision Lab Feature library

   ​ https://www.vlfeat.org/

   ​ https://github.com/vlfeat/vlfeat

   ​ https://www.vlfeat.org/overview/sift.html

   opencv

   https://github.com/opencv/opencv/blob/68d15fc62edad980f1ffa15ee478438335f39cc3/modules/features2d/src/sift.dispatch.cpp

   python：

   https://github.com/rmislam/PythonSIFT

首先看一下sift特征提取和匹配算法的步骤：

1. 构建尺度空间， 生成高斯差分金字塔（DOG金字塔）
2. 空间极值点检测(关键点的初步查探)
3. 稳健关键点的精确定位和筛选
   • 精确定位：离散空间的局部极值，通过迭代法寻找连续空间极值
   • 阈值筛选：消除小于阈值的极值
   • 消除边缘效应：去除掉可能是图片卷积过程中边缘效应产生的极值点
4. 稳定关键点的方向信息分配
5. 关键点描述
6. 特征点匹配

sift详细解释，大佬们已经讲的很好了，这里只记录下自己的一些困惑和理解，希望能帮助到和我一样的同学。(个人理解不一定完全对，欢迎纠错)

1. 高斯差分金字塔

sift为什么要建立高斯差分金字塔呢？弄清楚这个，得注意两个关键词：高斯， 差分。

高斯

熟悉图像处理的同学，都知道高斯卷积核能用来平滑图像，使图像变得模糊，如opencv中的GaussianBlur函数。而sift算法是为了找到尺度不变的特征点，就是图片中物体的特征点不随尺度变化，所以要得到物体不同尺度下的图片，对图片进行高斯卷积可以模拟物体在不同尺度下的图片(想象下，我们把图片缩小，图片中的物体会变小，同时也变得模糊)

差分

得到物体在不同尺度下的图片后，很明显下一步要做的就是特征提取，所谓图片特征就是特别突出的像素点，如比局部背景像素点都大的像素点，就是图像处理中常说的边缘，角点等。在图像处理中，laplacian算子能用来提取物体边缘，观察下laplacian卷积核，就能发现laplacian算子找到的这些边缘都是局部极值，如下是一个3x3的laplacian卷积核：(比周围四个像素点都大或都小才会被提取出来)

\[\left[ \begin{array}{ll} 0& 1& 0& \\ 1& -4& 1& \\ 0& 1& 0& \\ \end{array} \right] \]

所以总结下，先对图片进行高斯卷积，再进行拉普拉斯卷积，就能得到不同尺度下物体的特征，这一步骤常称为高斯拉普拉斯(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这里采用的就是类似的牛顿迭代法。牛顿迭代法的公式如下：

\[x^{(t+1)} = x^{(t)} - H_t^{-1}g_t\\ 其中g_t为x_t点的梯度，H_t为海森矩阵 \]

上述公式含义就是, 对于函数\(f(x)\)，在坐标点\(x^{(t)}\)，沿着\(- H_t^{-1}g_t\)方向移动，就能得到函数值比其更小的值，只要我们进行不断的迭代，坐标点就会收敛于极值点处。

所以回到我们的原始问题中，在图像上某一点处，我们是能计算该像素点处的一阶导数和二阶导数，也就能得到梯度\(g\)和海森矩阵\(H\), 那么在离散空间的极值点处，通过牛顿迭代法，就能得到连续空间的极值点。可以看下sift算法源码中是怎么计算的：

有x, y, s三个维度，所以梯度为三个维度的一阶偏导数，海森矩阵是三个维度的二阶偏导数， 即：

\[梯度：g = [dx, dy, ds]\\ 海森矩阵：H = \left[ \begin{array}{ll} dxx& dxy& dxs& \\ dyx& dyy& dys& \\ dsx& dsy& dss& \\ \end{array} \right] \]

看下代码实现：

梯度计算

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正是我们需要的值，而矩阵方程的解可以采用最小二乘法求解。

\[Hx = g\\ x = H^{-1}g \]

计算处\(- 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

迭代结束后，精确定位了极值点位置， 根据泰勒公式可以得到该位置对应极值点的近似值，泰勒公式的二次展开如下：

\[f(\vec X) = f(\vec X_0) + (\vec X-\vec X_0)^T g + \frac{1}{2}(\vec X-\vec X_0)^TH (\vec X-\vec X_0)^T\\ \]

上面式子中，\(\vec X = (x, y, s)\)表示尺度空间中某一点，即从\(\vec X_0\)移动到$\vec X \(处，根据我们迭代公式知道，移动值\)\vec X-\vec X_0\(为\)- H_t^{-1}g_t$，带入上述公式得到:

\[f(\vec X) = f(\vec X_0) + (\vec X-\vec X_0)^T g + \frac{1}{2}(\vec X-\vec X_0)^TH (\vec X-\vec X_0)^T\\ =f(\vec X_0) + (- H^{-1}g)^Tg + \frac{1}{2}(- H^{-1}g)^TH(- H^{-1}g)\\ =f(\vec X_0) + g^T(-H^{-1})g + \frac{1}{2}g^T(-H^{-1})H(- H^{-1}g)\\ = f(\vec X_0) + g^T(-H^{-1})g + \frac{1}{2}g^T(H^{-1}g)\\ =f(\vec X_0) + \frac{1}{2}g^T(-H^{-1}g) \]

对应代码中如下：

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方向的二阶导数，

\[海森矩阵：H = \left[ \begin{array}{ll} dxx& dxy& dxs& \\ dyx& dyy& dys& \\ dsx& dsy& dss& \\ \end{array} \right] \]

假设\(dxx=\alpha, dyy=\beta, \alpha = r \beta\), 则可以构建下面的等式， 很明显r=1时，这个表达式取得最小值，r大于1，或者小于1时，表达式值增大。所以只需要设置一个阈值，计算下面表达式的值，是否超过阈值即可。

\[\frac{(\alpha + \beta)^2}{\alpha \beta} = \frac{(r\beta + \beta)^2}{r\beta^2}=\frac{(r + 1)^2}{r} \]

代码中如下：

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中还有很多细节和参数值计算方法，多看几遍代码和讲解文章应该就清楚了。