角点检测：Harris 与 Shi-Tomasi

程序员阿德

图像特征类型
Harris 角点检测器
Shi-Tomasi 角点检测器
OpenCV 代码实现

Harris 角点检测
Shi-Tomasi 角点检测

参考

一、图像特征类型

我们都玩过拼图游戏，将大量混乱的小片段进行正确地排列，形成一个完整的图像。

玩这个游戏的关键点是什么？就是寻找那些区分度比较高的特征片段，它们更容易用于定位。

那这些有效特征到底是什么？

回答这个问题之前，先看一下下面这张图片：

在图像的顶部，给出了它的六个图像片段。我们试试如何在原始图像中找到这些图像片段的确切位置。

如果你亲自试过的话，可以发现：

A和B处于平坦区域，没有什么确切的特征，它们所在的位置有很多种可能；
C和D要相对简单一些，它们是建筑物的边缘，我们可以找到一个大致的位置，但是要定位到精确的位置仍然很难。所以边缘是更好的特征，但还不够好。
E和F是建筑的一些角落，可以很容易地发现它们的位置，因为对于建筑物角落这个图像片段，我们不管朝哪个方向移动，这个片段看起来都会不一样。

简化一下，就可以得到这张图：

蓝色矩形表示一个平坦区域，在各方向移动，窗口内像素值没有变化；
黑色矩形表示一个边缘特征（Edges），如果沿着垂直方向移动(梯度方向)，像素值会发生改变；如果沿着边缘移动(平行于边缘) ，像素值不会发生变化；
对于红色矩形框来说，它是一个角（Corners），不管你把它朝哪个方向移动，像素值都会发生很大变化。

三种情况分别如下图所示：

图像特征提供了图像丰富的信息。角点特征是图像中较好的特征，比边缘特征更好地用于定位。

对于拼图游戏，我们已经知道了什么是好的特征，但是下一个问题出现了，如何找到这些角点特征。

找到这些图像特征的过程被称为特征提取，特征提取决定了最终目标识别效果的好坏。

在图像的所有区域中，那些在所有方向上做微小移动，像素值变化都很大的区域，就是角点特征所在区域。

检测角点的应用主要有：图像对齐、图像拼接(拍摄全景图)、目标识别、3D重建、运动跟踪等等。

二、Harris 角点检测器

在图像处理中，检测角点特征的算法有很多，这里先介绍最常用的，也是最基础的 Harris 角点检测器(Harris Corner Detection)。

角点是两条边缘的交点，它表示两条边方向改变的地方，所以角点在任意一个方向上做微小移动，都会引起该区域的梯度图的方向和幅值发生很大变化。

也就是一阶导数(即灰度图的梯度)中的局部最大所对应的像素点就是角点。

于是我们可以利用这一点来检测角点。

使一个固定尺寸的窗口在图像上某个位置以任意方向做微小滑动，如果窗口内的灰度值（在梯度图上）都有较大的变化，那么这个窗口所在区域就存在角点。

这样就可以将 Harris 角点检测算法分为以下三步：

当窗口（小的图像片段）同时向 x 和 y 两个方向移动时，计算窗口内部的像素值变化量 E(u, v) ；
对于每个窗口，都计算其对应的一个角点响应函数 R；
然后对该函数进行阈值处理，如果 R > threshold，表示该窗口对应一个角点特征。

接下来对每一步进行详细介绍。

2.1、第一步

如何确定哪些窗口会引起较大的灰度值变化？

让一个窗口的中心位于灰度图像的一个位置 (x, y) ，这个位置的像素灰度值为 I(x, y) ，如果这个窗口分别向 x 和 y 方向移动一个小的位移u和v，到一个新的位置 (x+u, y+v) ，这个位置的像素灰度值就是 I(x+u, y+v) 。

[I(x+u, y+v) - I(x, y)] 就是窗口移动引起的灰度值的变化值。

设 w(x,y) 为位置 (x,y) 处的窗口函数，表示窗口内各像素的权重，最简单的就是把窗口内所有像素的权重都设为1。

有时也会把 w(x,y) 设定为以窗口中心为原点的高斯分布（二元正态分布）。如果窗口中心点像素是角点，那么窗口移动前后，中心点的灰度值变化非常强烈，所以该点权重系数应该设大一点，表示该点对灰度变化的贡献较大；而离窗口中心（角点）较远的点，这些点的灰度变化比较小，于是将权重系数设小一点，表示该点对灰度变化的贡献较小。

则窗口在各个方向上移动 (u,v) 所造成的像素灰度值的变化量公式如下：

E(u,v) = \sum\limits_{(x,y)} w(x,y) \times [I(x+u, y+v) - I(x,y)]^2 \\

对于一个角点来说， E(u,v) 会非常大。因此，我们可以最大化上面这个函数来得到图像中的角点。

用上面的函数计算 E(u,v) 会非常慢。因此，我们使用泰勒展开式（只有一阶）来得到这个公式的近似形式。

对于二维的泰勒展开式公式为：

T(x,y) \approx f(u,v) + (x-u)f_x(u,v) + (y-v)f_y(u,v) + ... \\

将 \displaystyle I(u+x,v+y) 套用上面的公式，可以得到：

I(x+u,y+v) \approx I(x,y) + uI_x + vI_y \\

其中 I_x 和 I_y 是 I 的偏微分，在图像中就是在 x 和 y 方向的梯度图（可以通过cv2.Sobel()来得到）：

I_x = \frac {\partial I(x,y)}{\partial x},\ \ I_y = \frac {\partial I(x,y)}{\partial y} \\

接下来继续推导：

\begin{aligned} E(u,v) &= \sum\limits_{(x,y)} w(x,y) \times [I(x,y) + uI_x + vI_y - I(x,y)]^2 \\ &= \sum\limits_{(x,y)} w(x,y) \times (uI_x + vI_y)^2 \\ &= \sum\limits_{(x,y)} w(x,y) \times (u^2I_x^2 + v^2I_y^2 + 2uvI_xI_y) \end{aligned} \\

把 u 和 v 拿出来，得到最终的形式：

E(u,v) \approx \begin{bmatrix} u, v \end{bmatrix} M \begin{bmatrix} u \\ v \end{bmatrix} \\

其中矩阵M为：

M = \sum\limits_{(x,y)} w(x,y) \begin{bmatrix} I_x^2 & I_xI_y \\ I_xI_y & I_y^2 \end{bmatrix} \rightarrow R^{-1} \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix} R \\

最后是把实对称矩阵对角化处理后的结果，可以把R看成旋转因子，其不影响两个正交方向的变化分量。

经对角化处理后，将两个正交方向的变化分量提取出来，就是 λ1 和 λ2（特征值）。

公式推导完了，现在回顾一下：

对于图像的每一个像素点 (x,y) ，对应一个以该像素为中心的窗口 w(x,y) ，然后该像素平移 (u,y) 得到新的像素点 (x+u, y+v) ，而 E(u,v) 就是窗口中所有像素的加权和乘以不同位置像素的灰度差值。

2.2、第二步

现在我们已经得到 E(u,v) 的最终形式，别忘了我们的目的是要找到会引起较大的灰度值变化的那些窗口。

灰度值变化的大小则取决于矩阵M，那么如何找到这些窗口，我们可以使用矩阵的特征值来实现。

计算每个窗口对应的得分（角点响应函数R）：

R = \text {det}(M) - k (\text {trace(M)})^2 \\

其中 \text {det}(M) = \lambda_1 \lambda_2 是矩阵的行列式， \text {trace(M)} = \lambda_1 + \lambda_2 是矩阵的迹。

λ1 和 λ2 是矩阵M的特征值， k 是一个经验常数，在范围 (0.04, 0.06) 之间。

2.3、第三步

根据 R 的值，将这个窗口所在的区域划分为平面、边缘或角点。为了得到最优的角点，我们还可以使用非极大值抑制。

注意：Harris 检测器具有旋转不变性，但不具有尺度不变性，也就是说尺度变化可能会导致角点变为边缘，如下图所示：

（插播一句，想要尺度不变特性的话，可以关注SIFT特征）

因为特征值 λ1 和 λ2 决定了 R 的值，所以我们可以用特征值来决定一个窗口是平面、边缘还是角点：

平面: 该窗口在平坦区域上滑动，窗口内的灰度值基本不会发生变化，所以 |R| 值非常小，在水平和竖直方向的变化量均较小，即 I_x 和 I_y 都较小，那么 λ1 和 λ2 都较小；
边缘: R 值为负数，仅在水平或竖直方向有较大的变化量，即 I_x 和 I_y 只有一个较大，也就是 λ1>>λ2 或 λ2>>λ1；
角点: R 值很大，在水平、竖直两个方向上变化均较大的点，即 I_x 和 I_y 都较大，也就是 λ1 和 λ2 都很大

用图片表示如下：

Harris 角点检测的结果是带有这些分数 R 的灰度图像，设定一个阈值，分数大于这个阈值的像素就对应角点。

三、Shi-Tomasi 角点检测器

知道了什么是 Harris 角点检测，后来有大佬在论文《Good_Features_to_Track》中提出了它的改进版——Shi-Tomasi 角点检测，Shi-Tomasi 方法在很多情况下可以得到比 Harris 算法更好的结果。

Harris 角点检测中每个窗口的分数公式是将矩阵 M 的行列式与 M 的迹相减：

R = λ_1λ_2 − k (λ_1 + λ_2)^2 \\

由于 Harris 角点检测算法的稳定性和 k 值有关，而 k 是个经验值，不好设定最佳值。

Shi-Tomasi 发现，角点的稳定性其实和矩阵 M 的较小特征值有关，于是直接用较小的那个特征值作为分数。这样就不用调整k值了。

所以 Shi-Tomasi 将分数公式改为如下形式：

R = min (λ_1, λ_2) \\

和 Harris 一样，如果该分数大于设定的阈值，我们就认为它是一个角点。

我们可以把它绘制到 λ1 ～ λ2 空间中，就会得到下图：

四、OpenCV 代码实现

4.1 Harris 角点检测

在opencv中有提供实现 Harris 角点检测的函数 cv2.cornerHarris，我们直接调用的就可以，非常方便。

函数原型：cv2.cornerHarris(src, blockSize, ksize, k[, dst[, borderType]])

对于每一个像素 (x,y)，在 (blockSize x blockSize) 邻域内，计算梯度图的协方差矩阵 M(x,y) ，然后通过上面第二步中的角点响应函数得到结果图。图像中的角点可以为该结果图的局部最大值。

即可以得到输出图中的局部最大值，这些值就对应图像中的角点。

涉及到几个参数：

src - 输入灰度图像，float32类型
blockSize - 用于角点检测的邻域大小，就是上面提到的窗口的尺寸
ksize - 用于计算梯度图的Sobel算子的尺寸
k - 用于计算角点响应函数的参数k，取值范围常在0.04~0.06之间

import cv2 as cv
import numpy as np

# detector parameters
block_size = 5
sobel_size = 3
k = 0.04

image = cv2.imread('bird.jpg')
gray_img = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

# modify the data type setting to 32-bit floating point 
gray_img = np.float32(gray_img)

# detect the corners with appropriate values as input parameters
corners_img = cv2.cornerHarris(gray_img, block_size, sobel_size, k)

# result is dilated for marking the corners, not necessary
dst = cv.dilate(corners_img, None)

# Threshold for an optimal value, marking the corners in Green
image[corners_img>0.01*corners_img.max()] = [0,0,255]

cv2.imwrite('new_bird.jpg', image)

4.2、Shi-Tomasi 角点检测

OpenCV 提供了 Shi-Tomasi 的函数： cv2.goodFeaturesToTrack()，来获取图像中前 N 个最好的角点。

函数原型如下：

goodFeaturesToTrack(image, maxCorners, qualityLevel, minDistance[, corners[, mask[, blockSize[, useHarrisDetector[, k]]]]])

其中的参数如下：

image：输入灰度图像，float32类型
maxCorners：返回角点的最大数目，值为0表表示没有设置最大值限制，返回所有检测到的角点。
qualityLevel：质量系数（小于1.0的正数，一般在0.01-0.1之间），表示可接受角点的最低质量水平。该系数乘以最好的角点分数（也就是上面较小的那个特征值），作为可接受的最小分数；例如，如果最好的角点分数值为1500且质量系数为0.01，那么所有质量分数小于15的角都将被忽略。
minDistance：角之间最小欧式距离，忽略小于此距离的点。
corners：输出角点坐标
mask：可选的感兴趣区域，指定想要检测角点的区域。
blockSize：默认为3，角点检测的邻域大小（窗口尺寸）
useHarrisDetector：用于指定角点检测的方法，如果是true则使用Harris角点检测，false则使用Shi Tomasi算法。默认为False。
k：默认为0.04，Harris角点检测时使用。

设定好这些参数，函数就能在图像上找到角点。所有低于质量水平的角点都会被忽略，然后再把合格角点按角点质量进行降序排列。

然后保留质量最高的一个角点，将它附近（最小距离之内）的角点都删掉（类似于非极大值抑制），按这样的方式最后得到 N 个最佳角点。

import numpy as np
import cv2

maxCorners = 100
qualityLevel = 0.01
minDistance = 10

img = cv2.imread('bird.jpg')
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
corners = cv2.goodFeaturesToTrack(gray, maxCorners, qualityLevel, minDistance)

corners = np.int0(corners)
for i in corners:
    x,y = i.ravel()
    cv2.circle(img,(x,y),2,(0,0,255),-1)
    
cv2.imwrite('new_bird.jpg', img)