角点检测与SIFT

关注【CV算法恩仇录】

角点检测

图像特征的好坏取决于信息量，信息量越丰富特征越好。接下来，举例说明图像的信息量，如下图中的矩形框，称为窗口，哪个窗口框住的图像具有信息量。从左往右看，第一张图片没有找到边缘，第二张图片找到边缘，第三张图片找到了一个两条边缘的交点，称为角点，窗口在这个角点任意方向稍微移动，像素值会巨大变化，第三张图片中窗口的信息量最大。

检测图像角点的经典的算法是 Harris 角点检测（Harris Corner Detector），它的思想是使用一个固定尺寸的窗口，在图中滑动，若窗口内的灰度值变化巨大，此窗口内区域存在角点，主要步骤如下：

• 计算窗口往 $x$ 和 $y$ 方向移动时像素值的变化量。
• 在矩形框中计算角点响应函数。
• 阈值处理，当函数值大于阈值，保留窗口特征，否则舍弃。

接下来详细介绍每个步骤，第一步计算窗口内像素值的变化量，给定一个窗口函数 $W$，是给窗口框住的像素值复制权重，常见是权重都为 $1$。定义一个函数 $E(u,v)$ 计算变化量，$u,v$ 是窗口整体的平移量，窗口函数的权重为 $1$，公式细节如下：

使用二阶泰勒公式对公式进行推导，结果如下

其中，$u,v$ 是常量， $r$ 是旋转因子，不影响公式的数值，${\lambda }_1$ 和 ${\lambda }_2$ 是矩阵的特征值，决定了公式的数值。由于特征值的计算复杂，Harris 把计算特征值转变为计算行列式和矩阵迹，称为角点响应函数 $R$，公式如下，$k$ 是一个常数，在 $0.04$ 与 $0.06$ 之间。

两个特征值与角点的关系是：当两个特征值都很小时，窗口滑动的是平坦区域，灰度值基本没有变化；当其中一个特征值远大于另一个特征值时，窗口滑动到边缘；当两个特征值都比较大，且差别不大，窗口滑动到角点，最后通过阈值选择最终的角点。

OpenCV里cornerHarris函数实现了角点检测，以一个棋牌的图片为例子，函数的代码示例如下，检测的角点展示在棋盘格上的红点。

import cv2 as cv

filename = './书籍的图片/棋牌格.png' 
img = cv.imread(filename)
img_gray = cv.cvtColor(img, cv.COLOR_BGR2GRAY)

img_gray = np.float32(img_gray)
dst = cv.cornerHarris(img_gray, 2, 3, 0.04)
print('before dst: ', len(dst))
# 形态学：膨胀操作 dilate()
dst = cv.dilate(dst, None)
print('after dst: ', len(dst))

# 阈值处理
img[dst > 0.01 * dst.max()] = [0, 0 , 255]
cv.imwrite('HarrisCorner.png', img)

尺度不变特征变换

角点具备旋转不变性和光照不变性，角点区域旋转之后还是会被检测到角点，不同程度的光照在角点区域，对其影响不大。

但在不同远近的距离看角点，角点的形状会变化，如图中 A 的区域是角点，B 的区域是近距离观察 A，原先相对尖锐的角点变成了圆滑的边缘，B 区域不会被检测成角点，而是边缘，角点不具备尺度不变性，尺度可以抽象理解为放大和缩小。

后来 David Lowe 提出尺度不变特征变换（Scale Invariant Feature Transform，SIFT），在不同的尺度空间上查找关键点(特征点)，并计算出关键点的方向。

具体步骤如下：

• 构建尺度空间：确保特征与尺度无关
• 关键点定位：识别合适的特征或关键点
• 关键点方向分配：确保关键点是旋转不变的
• 关键点描述符：为每个关键点分配唯一的指纹

构建尺度空间

我们常用高斯模糊技术来减少图像中的噪声，如图所示，左图是原图，右图经过高斯模糊之后，被去除了纹理和次要细节，仅保留形状和边缘等相关信息。

高斯模糊算法成功地去除了图像中的噪声，也突出了图像的重要特征，为确保这些特征不能依赖于尺度，可以通过创建“尺度空间”在多个尺度上搜索这些特征。

尺度空间是从单个图像生成的具有不同尺度的图像的集合。对于尺寸为 $(275,183)$ 的原始图像（original image），将比例缩小一半，尺寸为 $(138,92)$ 的缩放图像（scaled image），再对原始图像和缩放图像进行高斯模型，例如下图。

读者可能在想——需要对图像进行多少次缩放，以及对每个缩放图像进行多少次模糊处理？SIFT 论文提到缩放4次，每个缩放图像进行5次模糊处理。

通过不同的尺度因子的高斯核与图像做卷积获得高斯金字塔，如公式所示
$$L(x,y,\sigma )=G(x,y,\sigma )\ast I(x,y)$$
其中：
$$G(x,y,\sigma )=\frac{1}{2\pi {\sigma }^2}{e}^{-(x^2+y^2)/2{\sigma }^2}$$

$\sigma$ 是尺度空间因子，值大对于大尺度，值小对于小尺度。一个图像的尺度空间 $L(x,y,\sigma )$，定义为原始图像 $I(x,y)$ 与一个可变尺度的 2 维高斯函数 $G(x,y,\sigma )$ 卷积运算。

代码实现如下：

import numpy as np

def gaussian_filter(sigma):
	size = 2*np.ceil(3*sigma)+1
	x, y = np.mgrid[-size//2 + 1:size//2 + 1, -size//2 + 1:size//2 + 1]
	g = np.exp(-((x**2 + y**2)/(2.0*sigma**2))) / (2*np.pi*sigma**2)
	return g/g.sum()

s = 3
sigma = 1.6
k = 2**(1/s)
kernel = gaussian_filter(k * sigma)

def generate_octave(init_level, s, sigma):
	octave = [init_level]
	k = 2**(1/s)   # 1.2599210498948732
	kernel = gaussian_filter(k * sigma)
  for i in range(s+2):
    next_level = convolve(octave[-1], kernel)
    octave.append(next_level)

  return octave
  
def generate_gaussian_pyramid(im, num_octave, s, sigma):
    pyr = []
    for _ in range(num_octave):
        octave = generate_octave(im, s, sigma)

        pyr.append(octave)
        im = octave[-3][::2, ::2]

    return pyr
    
num_octave=4
gaussian_pyr = generate_gaussian_pyramid(im, num_octave, s, sigma)

相邻两层高斯金字塔的比列为 $k$，如第一层的尺度空间因子是 $\sigma$，第二层是 $k\sigma$，第三层是 $k^2\sigma$，以此类推。高斯差分金字塔是不同层的高斯金字塔做减法，公式如下，
$$D(x,y,\sigma )=(G(x,y,k\sigma )-G(x,y,\sigma ))\ast I(x,y)$$
效果图：

代码实现如下：

def generate_DoG_octave(gaussian_octave):
    octave = []
    for i in range(1, len(gaussian_octave)):
        octave.append(gaussian_octave[i] - gaussian_octave[i-1])

    return np.concatenate([o[:,:,np.newaxis] for o in octave], axis=2)

def generate_DoG_pyramid(gaussian_pyramid):
    pyr = []
    for gaussian_octave in gaussian_pyramid:
        pyr.append(generate_DoG_octave(gaussian_octave))

    return pyr

关键点定位

创建高斯金字塔后，下一步是从图像中找到可用于特征匹配的重要关键点。这个想法是找到图像的局部最大值和最小值。这部分分为两个步骤：

1. 找到局部最大值和最小值
2. 删除低对比度关键点（关键点选择）

局部最大值和局部最小值

为了定位局部最大值和最小值，我们遍历图像中的每个像素并将其与其相邻像素进行比较。

“相邻”像素是指该像素在当前图像的周围 8 个像素、以及尺度空间中上下两层中的相邻图像的 18（2x9）个点，总共 26 个像素值。

将每个像素值与其他 26 个像素值进行比较，以确定它是否是局部最大值/最小值。例如，在下图中，将标记为 x 的像素与相邻像素（绿色）进行比较，如果它是相邻像素中最高或最低的，则将其选为关键点：

关键点选择

我们已经成功生成了尺度不变的关键点，但是其中一些关键点可能对噪声不具有鲁棒性。因此，我们将消除对比度低或非常靠近边缘的关键点。
为了处理低对比度关键点，为每个关键点计算二阶泰勒展开。如果结果值小于阈值（一般为 0.03 或 0.04）就会被忽略掉。 在 OpenCV 中这种阈值被称为 contrastThreshold。

DoG 算法对边界非常敏感， 所以我们必须要把边界去除。 Harris 算法除了可以用于角点检测之外还可以用于检测边界。从 Harris 角点检测的算法中，当一个特征值远远大于另外一个特征值时检测到的是边界。那在DoG算法中欠佳的关键点在平行边缘的方向有较大的主曲率，而在垂直于边缘的方向有较小的曲率，两者的比值如果高于某个阈值（在 OpenCV 中叫做边界阈值），就认为该关键点为边界，将被忽略，一般将该阈值设置为10。

将低对比度和边界的关键点去除，得到的就是我们感兴趣的关键点。

关键点方向确定

经过上述两个步骤，图像的关键点就完全找到了，这些关键点具有尺度不变性。为了实现旋转不变性，还需要为每个关键点分配一个方向角度，也就是根据检测到的关键点所在高斯尺度图像的邻域结构中求得一个方向基准。

对于任一关键点，我们采集其所在高斯金字塔图像以 $r$ 为半径的区域内所有像素的梯度特征（幅值和幅角），半径 $r$ 为：
$$r=3\times 1.5\sigma$$

其中 $\sigma$ 是关键点所在 octave 的图像的尺度，可以得到对应的尺度图像。
梯度的幅值和方向的计算公式为：

接下来用一个例子来说明计算过程，如下图所示：

假设我们要找到红色像素值的大小和方向，为此，通过取 55 & 46 和 56 & 42 之间的差来计算 $x$ 和 $y$ 方向的梯度。这分别是 $G_x=9$ 和 $G_y=14$。一旦有了梯度，我们就可以使用公式找到幅值和方向：
$$幅值=\sqrt{G_x^2+G_y^2}=16.64$$
$$方向=arctan(G_y/G_x)=57.17$$

完成关键点梯度计算后，使用直方图统计关键点邻域内像素的梯度幅值和方向。具体做法是，将 360° 分为 36 柱，每 10° 为一柱，然后在以 $r$ 为半径的区域内，将梯度方向在某一个柱内的像素找出来，然后将他们的幅值相加在一起作为柱的高度。因为在 $r$ 为半径的区域内像素的梯度幅值对中心像素的贡献是不同的，因此还需要对幅值进行加权处理，采用高斯加权，方差为 $1.5\sigma$。如下图所示，为简化图中只画了8个方向的直方图。

每个特征点必须分配一个主方向，还需要一个或多个辅方向，增加辅方向的目的是为了增强图像匹配的鲁棒性。当一个柱体的高度大于主方向柱体高度的80%时，则该柱体所代表的的方向就是该特征点的辅方向。

直方图的峰值，即最高的柱代表的方向是特征点邻域范围内图像梯度的主方向，但该柱体代表的角度是一个范围，所以我们还要对离散的直方图进行插值拟合，以得到更精确的方向角度值。利用抛物线对离散的直方图进行拟合，如下图所示：

获得图像关键点主方向后，每个关键点有三个信息(x,y,σ,θ)：位置、尺度、方向。由此我们可以确定一个SIFT特征区域。通常使用一个带箭头的圆或直接使用箭头表示SIFT区域的三个值：中心表示特征点位置，半径表示关键点尺度，箭头表示方向。如下图所示：

关键点描述

通过以上步骤，每个关键点就被分配了位置，尺度和方向信息。接下来我们为每个关键点建立一个描述符，该描述符既具有可区分性，又具有对某些变量的不变性，如光照，视角等。而且描述符不仅仅包含关键点，也包括关键点周围对其有贡献的的像素点。主要思路就是通过将关键点周围图像区域分块，计算块内的梯度直方图，生成具有特征向量，对图像信息进行抽象。
首先在关键点周围取一个 16×16 的邻域。这个 16×16 块进一步分为 4×4 子块，对于每个子块，我们使用幅度和方向生成直方图。在这个阶段，bin 的大小增加了，我们只取了 8 个 bin（不是 36 个）。这些箭头中的每一个都代表 8 个 bin，箭头的长度定义了幅度。因此，我们将为每个关键点提供总共 128 个 bin 值。

OpenCV 中的 SIFT

现在让我们看看 OpenCV 中可用的 SIFT 功能。请注意，这些以前仅在opencv contrib repo中可用，但专利于 2020 年到期。让我们从关键点检测开始并绘制它们。首先，我们必须构造一个 SIFT 对象。我们可以向它传递不同的参数，这些参数是可选的，它们在文档中有很好的解释。

import numpy as np
import cv2 as cv
img = cv.imread('home.jpg')
gray= cv.cvtColor(img,cv.COLOR_BGR2GRAY)
sift = cv.SIFT_create()
kp = sift.detect(gray,None)
img=cv.drawKeypoints(gray,kp,img)
cv.imwrite('sift_keypoints.jpg',img)

OpenCV 还提供了cv.drawKeyPoints()函数，它在关键点的位置上绘制小圆圈。如果将标志cv.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS传递给它，它将绘制一个关键点大小的圆，甚至会显示其方向。请参见下面的示例。

img=cv.drawKeypoints(gray,kp,img,flags=cv.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
cv.imwrite('sift_keypoints.jpg',img)

参考文献：
【1】https://www.analyticsvidhya.com/blog/2019/10/detailed-guide-powerful-sift-technique-image-matching-python/
【2】https://docs.opencv.org/4.x/da/df5/tutorial_py_sift_intro.html
【3】https://medium.com/data-breach/introduction-to-sift-scale-invariant-feature-transform-65d7f3a72d40