Java教程

sift算法

本文主要是介绍sift算法,对大家解决编程问题具有一定的参考价值,需要的程序猿们随着小编来一起学习吧!

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

这篇关于sift算法的文章就介绍到这儿,希望我们推荐的文章对大家有所帮助,也希望大家多多支持为之网!