计算机视觉基础（七）—— 一文解析Harris角点检测

在图像处理领域中，特征点又被称为兴趣点或者角点，它通常具有旋转不变性和光照不变性和视角不变性等优点，是图像的重要特征之一，常被应用到目标匹配、目标跟踪、三维重建等应用中。点特征主要指图像中的明显点，如突出的角点、边缘端点、极值点等等，用于点特征提取的算子称为兴趣点提取（检测）算子，常用的有Harris角点检测、FAST特征检测、SIFT特征检测及SURF特征检测。

1 什么是角点检测

1.1 角点

定义：角点就是轮廓之间的交点。
如果从数字图像处理的角度来描述就是：像素点附近区域像素无论是在梯度方向、还是在梯度幅值上都发生较大的变化。
基本思想：选取一个固定的窗口在图像上以任意方向的滑动，如果灰度都有较大的变化，那么久认为这个窗口内部存在角点。

使用一个滑动窗口在下面三幅图中滑动，可以得出以下结论：

• 左图表示一个平坦区域，在各方向移动，窗口内像素值均没有太大变化；
• 中图表示一个边缘特征（Edges），如果沿着水平方向移动(梯度方向)，像素值会发生跳变；如果沿着边缘移动(平行于边缘) ，像素值不会发生变化；
• 右图表示一个角（Corners），不管你把它朝哪个方向移动，像素值都会发生很大变化。
  
  所以，右图是一个角点。

1.2 角点类型

下图展示了不同角点的类型，可以发现：如果使用一个滑动窗口以角点为中心在图像上滑动，存在朝多个方向上的移动会引起该区域的像素值发生很大变化的现象。

2 图像梯度

像素值发生很大变化这一现象可以用图像梯度进行描述。在图像局部内，图像梯度越大表示该局部内像素值变化越大（灰度的变化率越大）。
而图像的梯度在数学上可用微分或者导数来表示。对于数字图像来说，相当于是二维离散函数求梯度，并使用差分来近似导数：
$$G_x(x,y)=H(x+1,y)-H(x-1,y)$$ $$G_y(x,y)=H(x,y+1)-H(x,y-1)$$
在实际操作中，对图像求梯度通常是考虑图像的每个像素的某个邻域内的灰度变化，因此通常对原始图像中像素某个邻域设置梯度算子，然后采用小区域模板进行卷积来计算，常用的有Prewitt算子、Sobel算子、Robinson算子、Laplace算子等。

Sobel算子

我们可以使用$3\times 3$ 的卷积核来进行图像求导：
$$\begin{array}{rlr}& G^{\prime }y=\left[\begin{array}{ccc}+1& +2& +1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right]\ast I& \text{ 和 }{G}^{\prime }x=\left[\begin{array}{ccc}+1& +2& +1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right]\ast I\end{array}$$
中$\mathbf{I}$表示原图片，${\mathbf{G}}_x^{\prime }$和${\mathbf{G}}_y^{\prime }$分别表示沿图片水平和竖直方向上的变化，$\ast$表示卷积操作。

下面以Sobel算子为例讲述如何计算梯度
x和y方向的Sobel算子分别为：
$$G_x=\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]G_y=\left[\begin{array}{ccc}1& 2& 1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right]$$

$\mathbf{G}$中一个3x3的窗口为A，要计算梯度的像素点为e，则和Sobel算子进行卷积之后，像素点e在x和y方向的梯度值分别为：
$$G^{\prime }x=Gx\ast A=\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]\ast \left[\begin{array}{ccc}a& b& c\\ d& e& f\\ g& h& i\end{array}\right]=\mathrm{sum}\left(\left[\begin{array}{ccc}-a& 0& c\\ -2d& 0& 2f\\ -g& 0& i\end{array}\right]\right)$$ $$G_y^{\prime }=G_y\ast A=\left[\begin{array}{ccc}1& 2& 1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right]\ast \left[\begin{array}{ccc}a& b& c\\ d& e& f\\ g& h& i\end{array}\right]=\mathrm{sum}\left(\left[\begin{array}{ccc}a& 2b& c\\ 0& 0& 0\\ -g& -2h& -i\end{array}\right]\right)$$
中 * 为卷积符号，sum表示矩阵中所有元素相加求和。

3 Harris角点检测算法原理

算法思想
算法的核心是利用局部窗口在图像上进行移动，判断灰度是否发生较大的变化。如果窗口内的灰度值（在梯度图上）都有较大的变化，那么这个窗口所在区域就存在角点。
这样就可以将 Harris 角点检测算法分为以下三步：

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

3.1 Step 01 — 建立数学模型

第一步是通过建立数学模型，确定哪些窗口会引起较大的灰度值变化。
让一个窗口的中心位于灰度图像的一个位置$(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 _{x,y}w(x,y)[I(x+u,y+v-I(x,y){]}^2$$
其中灰度变化部分：
$$\begin{array}{rl}\sum [I(x+u,y+v-I(x,y){]}^2& \approx \sum {\left[I(x,y)+u{I}_x+v{I}_y-I(x,y)\right]}^2\\ & =\sum u^2{I}_x^2+2uv{I}_x{I}_y+v^2{I}_y^2\\ & =\sum \left[\begin{array}{ll}u& v\end{array}\right]\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]\left[\begin{array}{l}u\\ v\end{array}\right]\\ & =\left[\begin{array}{ll}u& v\end{array}\right]\left(\sum \left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]\right)\left[\begin{array}{l}u\\ v\end{array}\right]\end{array}$$

因此得到：
$$\begin{array}{c}E(u,v)=\left[\begin{array}{ll}u& v\end{array}\right]M\left[\begin{array}{l}u\\ \\ v\end{array}\right]\end{array}$$

$$M=\sum _{(x,y)}w(x,y)\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]\to R^{-1}\left[\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right]R$$

函数 $w(x,y)$，最简单的取值为1，一般为以中心为原点的二元正态分布。
最后是把实对称矩阵对角化处理后的结果，可以把R看成旋转因子，其不影响两个正交方向的变化分量。
经对角化处理后，将两个正交方向的变化分量提取出来，就是 λ1 和 λ2（特征值）。这里利用了线性代数中的实对称矩阵对角化的相关知识。

3.2 Step 02 — 角点响应函数R

现在我们已经得到 $E(u,v)$的最终形式，我们的目的是要找到会引起较大的灰度值变化的那些窗口。
灰度值变化的大小则取决于矩阵M，M为梯度的协方差矩阵。在实际应用中为了能够应用更好的编程，所以定义了角点响应函数R，通过判定R大小来判断像素是否为角点。
计算每个窗口对应的得分（角点响应函数R定义）：
$$R=\det (M)-k(\mathrm{trace}(\mathrm{M}){)}^2$$ 其中 $det(M)={\lambda }_1{\lambda }_2$是矩阵的行列式， $trace(M)={\lambda }_1+{\lambda }_2$ 是矩阵的迹。

$\lambda 1$ 和 $\lambda 2$ 是矩阵$M$的特征值， $k$是一个经验常数，在范围 (0.04, 0.06) 之间。
$R$的值取决于$M$的特征值，对于角点$|R|$很大，平坦的区域$|R|$很小，边缘的$R$为负值。

3.3 Step 03 — 角点判定

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

注意：Harris 检测器具有旋转不变性，但不具有尺度不变性，也就是说尺度变化可能会导致角点变为边缘。想要尺度不变特性的话，可以关注SIFT特征。
因为特征值 λ1 和 λ2 决定了 R 的值，所以我们可以用特征值来决定一个窗口是平面、边缘还是角点：

• 平面:：该窗口在平坦区域上滑动，窗口内的灰度值基本不会发生变化，所以 $|R|$ 值非常小，在水平和竖直方向的变化量均较小，即 $I_x$和 $I_y$都较小，那么 λ1 和 λ2 都较小；
• 边缘：$|R|$值为负数，仅在水平或竖直方向有较大的变化量，即 $I_x$和 $I_y$只有一个较大，也就是 λ1>>λ2 或 λ2>>λ1；
• 角点：[公式] 值很大，在水平、竖直两个方向上变化均较大的点，即 $I_x$和 $I_y$ 都较大，也就是 λ1 和 λ2 都很大。
  如下图所示：
  
  Harris 角点检测的结果是带有这些分数 R 的灰度图像，设定一个阈值，分数大于这个阈值的像素就对应角点。

4 基于OpenCV的实现

4.1 API

在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之间
• dst - 函数调用后的运算结果存在这里，即这个参数用于存放Harris角点检测的输出结果，和原图片有一样的尺寸和类型
• borderType - int类型，图像像素的边界模式，注意它有默认值BORDER_DEFAULT

4.2 C++代码示例

#include <opencv2/opencv.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
using namespace cv;

#define WINDOW_NAME1 "窗口1"
#define WINDOW_NAME2 "窗口2"

Mat g_srcImage, g_srcImage1, g_grayImage;
int g_nThresh = 30; //当前阀值
int g_nMaxThresh = 175;

void on_fCornerHarris(int, void*)
{
 Mat dstImage; //目标图
 Mat normImage; //归一化后的图
 Mat scaledImage; //线性变换后的八位无符号整形的图
 //置零当前需要显示的两幅图，即清除上一次调用此函数时他们的值
 dstImage = Mat::zeros(g_srcImage.size(), CV_32FC1);
 g_srcImage1 = g_srcImage.clone();
 //角点检测
 cornerHarris(g_grayImage, dstImage, 2, 3, 0.04, BORDER_DEFAULT);
 normalize(dstImage, normImage, 0, 255, NORM_MINMAX, CV_32FC1, Mat());
 //将归一化后的图线性变换成8位无符号整形
 convertScaleAbs(normImage, scaledImage);
 //将检测到的，且符合阀值条件的角点绘制出来
 for (int i = 0; i < normImage.rows; i++)
 {
  for (int j = 0; j < normImage.cols; j++)
  {
   if ((int)normImage.at<float>(i, j) > g_nThresh + 80)
   {
    circle(g_srcImage1, Point(j, i), 5, Scalar(10, 10, 255), 2, 8, 0);
    circle(scaledImage, Point(j, i), 5, Scalar(0, 10, 255), 2, 8, 0);
   }
  }
 }
 imshow(WINDOW_NAME1, g_srcImage1);
 imshow(WINDOW_NAME2, scaledImage);
}

int main()
{
 g_srcImage = imread("sample.jpg");
 g_srcImage.copyTo(g_srcImage1);
 cvtColor(g_srcImage1, g_grayImage, COLOR_BGR2GRAY);
 namedWindow(WINDOW_NAME1, WINDOW_AUTOSIZE);
 createTrackbar("阀值：", WINDOW_NAME1, &g_nThresh, g_nMaxThresh, on_fCornerHarris);
 on_fCornerHarris(0, NULL);
 waitKey(0);
 return 0;
}

4.3 Python代码示例

import cv2
from matplotlib import pyplot as plt
import numpy as np

# detector parameters
block_size = 3
sobel_size = 3
k = 0.06

image = cv2.imread('image/sample.jpg')
img = cv2.cvtColor(image, cv2.COLOR_BGR2RGB) 
plt.rcParams['savefig.dpi'] = 600 #图片像素
plt.rcParams['figure.dpi'] = 600 #分辨率
plt.subplot(121)
plt.imshow(img)

print(image.shape)
height = image.shape[0]
width = image.shape[1]
channels = image.shape[2]
print("width: %s  height: %s  channels: %s"%(width, height, channels))
   
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
kernel = cv2.getStructuringElement(cv2.MORPH_RECT,(3, 3))
dst = cv2.dilate(corners_img, kernel)

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

for r in range(height):
        for c in range(width):
            pix=dst[r,c]
            if pix>0.05*dst.max():
               cv2.circle(image,(c,r),5,(0,0,255),0)

image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)

plt.subplot(122)
plt.imshow(image)
plt.show()

5 Harris角点检测的性质

Harris角点检测的性质可总结如下：

• 阈值决定角点的数量
• Harris角点检测对亮度和对比度的变化不敏感（光照不变性）
  在进行Harris角点检测时，使用了微分算子对图像进行微分运算，而微分运算对图像密度的拉升或收缩和对亮度的抬高或下降不敏感。换言之，对亮度和对比度的仿射变换并不改变Harris响应的极值点出现的位置，但是，由于阈值的选择，可能会影响角点检测的数量。
  
• Harris角点检测算子具有旋转不变性
  Harris角点检测算子使用的是角点附近的区域灰度二阶矩矩阵。而二阶矩矩阵可以表示成一个椭圆，椭圆的长短轴正是二阶矩矩阵特征值平方根的倒数。当特征椭圆转动时，特征值并不发生变化，所以判断角点响应值也不发生变化，由此说明Harris角点检测算子具有旋转不变性。
• Harris角点检测算子不具有尺度不变性
  尺度的变化会将角点变为边缘，或者边缘变为角点，Harris的理论基础并不具有尺度不变性。
  

参考

• https://blog.csdn.net/weixin_40647819
• 论文：《C.Harris, M.Stephens. “A Combined Corner and Edge Detector”. Proc. of 4th Alvey Vision Conference》
• Harris角点算法
• 角点检测：Harris 与 Shi-Tomasi
• https://www.cnblogs.com/ronny/p/4009425.html
• https://www.jianshu.com/p/f0cd0f6e63a9
• https://blog.csdn.net/qq_27396861/article/details/87898892

关于Datawhale：

Datawhale是一个专注于数据科学与AI领域的开源组织，汇集了众多领域院校和知名企业的优秀学习者，聚合了一群有开源精神和探索精神的团队成员。Datawhale以“for the learner，和学习者一起成长”为愿景，鼓励真实地展现自我、开放包容、互信互助、敢于试错和勇于担当。同时Datawhale 用开源的理念去探索开源内容、开源学习和开源方案，赋能人才培养，助力人才成长，建立起人与人，人与知识，人与企业和人与未来的联结。