CV之Harris特征点检测器-兴趣点检测（详解）

1.1 简介

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

本次博客较为常用而且较为基础的Harris角点检测算法，它的思想以及数学理论能够很好地帮助我们了解兴趣点检测的相关原理。

1.2内容介绍

1.2.1基础知识

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

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

所以，右图是一个角点。

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

3.图像梯度

图像的梯度与数学中的梯度在形式上是由差异的，原因在于图像的特点，图像是一个离散的二维函数，接下来会尝试说明这一点，但从意义上来说是相同的，是为了表现图像灰度的变化率。
在数学微积分中，一维函数的一阶微分的基本定义是：
dfdx=lim⁡ϵ→0f(x+ϵ)−f(x)ϵ\frac{d f}{d x}=\lim _{\epsilon \rightarrow 0} \frac{f(x+\epsilon)-f(x)}{\epsilon}dxdf=ϵ→0limϵf(x+ϵ)−f(x)
二维函数的一阶微分的基本定义是：
∂f(x,y)∂x=lim⁡ϵ→0f(x+ϵ,y)−f(x,y)ϵ∂f(x,y)∂y=lim⁡ϵ→0f(x,y+ϵ)−f(x,y)ϵ\begin{aligned} \frac{\partial f(x, y)}{\partial x} &=\lim _{\epsilon \rightarrow 0} \frac{f(x+\epsilon, y)-f(x, y)}{\epsilon} \\ \frac{\partial f(x, y)}{\partial y} &=\lim _{\epsilon \rightarrow 0} \frac{f(x, y+\epsilon)-f(x, y)}{\epsilon} \end{aligned}∂x∂f(x,y)∂y∂f(x,y)=ϵ→0limϵf(x+ϵ,y)−f(x,y)=ϵ→0limϵf(x,y+ϵ)−f(x,y)
对于灰度图这样二维数组的图像，它其实就是一个离散的二位函数，说它离散，是因为每个灰度值取值为整数而不是小数，因此ϵ不能无限小，而ϵ的最小单位即是1像素。因此当ϵ取最小值1时，灰度图（离散的二维函数）的一阶微分基本定义是：
∂f(x,y)∂x=f(x+1,y)−f(x,y)=gx∂f(x,y)∂y=f(x,y+1)−f(x,y)=gy\begin{array}{l} \frac{\partial f(x, y)}{\partial x}=f(x+1, y)-f(x, y)=g x \\ \frac{\partial f(x, y)}{\partial y}=f(x, y+1)-f(x, y)=g y \end{array}∂x∂f(x,y)=f(x+1,y)−f(x,y)=gx∂y∂f(x,y)=f(x,y+1)−f(x,y)=gy
我们可以观察到，图像在(x, y)点处x方向和y方向上的梯度相当于2个相邻像素之间的差值。而且一般地，我们会使用两个方向上梯度的绝对值，因为我们只关心图像的灰度值变化大小，而梯度方向只需要知道是水平还是垂直即可，无须知晓到底是具体是向左还是向右，向下还是向上。
计算图像某点的梯度。
上面，我们介绍了图像某个点的水平梯度gx和垂直梯度gy的计算方式，而将两者结合起来即可得到图像中该点的梯度的大小与方向。
a. 计算梯度的大小，在数学上，这种结合的公式是：
M(x,y)=(gx)2+(gy)2M(x, y)=\sqrt{(g x)^{2}+(g y)^{2}}M(x,y)=(gx)2+(gy)2
但是图像处理中经常使用以下方式（用绝对值加和来近似平方和平方根）来减小计算开支。
M(x,y)=∣gx∣+∣gy∣M(x, y)=|g x|+|g y|M(x,y)=∣gx∣+∣gy∣

再进一步扩展，图像的梯度也常用来判断某个窗口区域是平坦、边缘还是角点。而在边缘检测中，图像某一点的梯度方向是与边缘方向相垂直的，所以知道了梯度的方向，就晓得了边缘的方向，在以后的边缘检测、轮廓检测等知识中会用到。
还有，希望大家记住梯度的计算是以每一个像素点为中心，而非一个区域（比如 3×3）。

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

1.2.2 Sobel算子简介

sobel算子用来检测主体与主体之间的边缘特征，其效果有点类似高斯滤波，能够使得离中心点越近的像素权重越大，越重要。
需要解释的是，我们刚才介绍了，图像在某一点A的梯度，只需要计算其相邻两个像素点的差值即可得到X方向Y方向上的梯度；
但为了更加合理的判断某一点的梯度大小，需要借助算子为这个以点A为中心的ksize×ksize区域（一般是ksize取3或5或者其他奇数）添加权重，来扩大差异，增强边缘检测效果。
而且在运用中，因为计算一点的梯度需要借助周边区域的多个点的灰度值，因而对图像中的噪点比较敏感，因此时常在调用sobel方法进行计算之前，需要先借助高斯滤波器或者其他滤波器对原图像进行平滑\模糊处理以降噪。
所以完整的Sobel梯度计算过程是：1.高斯滤波器对原图像平滑降噪GaussianBlur 2.转灰度图cvtColor3.求X和Y方向的梯度Sobel。
sobel算子的形式
sx=[−101−202−101]sy=[−12−1000121]\begin{array}{c} s_{x}=\left[\begin{array}{ccc} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{array}\right] \\ s_{y}=\left[\begin{array}{ccc} -1 & 2 & -1 \\ 0 & 0 & 0 \\ 1 & 2 & 1 \end{array}\right] \end{array}sx=⎣⎡−1−2−1000121⎦⎤sy=⎣⎡−101202−101⎦⎤
在OpenCV - Python 中，调用cv2.Sobel方法，一般的sobel算子是，右减左，下减上，进而得到图像中某一点的x方向上的梯度或者y方向上的梯度。另外，因为进行减法操作后不一定得到的是正值，也可能是负值，但是uint8不允许负值的存在，因此要指定ddepth = cv2.CV_64F 以包容负值的存在（eg:cv2.Sobel(img,cv2.CV_64F,1,0,ksize=3)），此时进行图像的展示，会先将数组中的负值视为0。但是随后会调用其他方法(cv2.convertScaleAbs)以进行绝对值化的操作，此时进行展示时，会使得负值转化为正值，得到我们完整的特征图。
sobel算子的计算梯度的数学计算类似于卷积神经网络的卷积计算，而具体的计算过程参考sobel算子原理与实现.在此需要注意，sobel算子是算梯度的，虽然计算涉及很多的像素点，但是它计算的结果是一个值，是计算区域中心坐标的梯度值！
sobel算子的计算梯度时，对于边缘点的梯度，一般会做填充处理，即如果算子核大小为3的话，会在图像四边界填充一层像素值为0的点。
不明之处：我们这篇文章其实是打算讲CornerHarris的参数的。而sobel算子是Harris角点检测过程中使用的算子，但是CornerHarris的参数中似乎只能指定sobel算子的大小即ksize，而不能指定具体使用的核是啥，而这一点有待深挖，因为我也看到过自定义算子，比如，高斯滤波器曾作为图像降噪的方式，也可用来计算梯度并赋予到某些方法（eg.erode）中的情况。

1.2.3滤波器

主要想推荐一篇文章图像平滑去噪之高斯滤波器。具体的滤波器内容内容比较广泛，读者自行学习。
高斯滤波器的形式。由于高斯滤波实质是一种加权平均滤波，为了实现平均，核还带有一个系数，例如上图中的十六分之一、八十四分之一。这些系数等于矩阵中所有数值之和的倒数。相比之下，sobel算子的核在计算时，只是简单的加权求和。
116[121242121]184[1232125652368632565212321]\frac{1}{16}\left[\begin{array}{lll} 1 & 2 & 1 \\ 2 & 4 & 2 \\ 1 & 2 & 1 \end{array}\right] \quad \frac{1}{84}\left[\begin{array}{lllll} 1 & 2 & 3 & 2 & 1 \\ 2 & 5 & 6 & 5 & 2 \\ 3 & 6 & 8 & 6 & 3 \\ 2 & 5 & 6 & 5 & 2 \\ 1 & 2 & 3 & 2 & 1 \end{array}\right]161⎣⎡121242121⎦⎤841⎣⎢⎢⎢⎢⎡1232125652368632565212321⎦⎥⎥⎥⎥⎤
高斯滤波的特点。因为高斯滤波是以空间距离来确定权重大小的，但并没有考虑用颜色距离来确定权重，这样就导致了高斯滤波在去除噪声的同时，也一定程度上模糊了边界。

1.2.4 Harris角点检测算法原理

1. 算法思想
算法的核心是利用局部窗口在图像上进行移动，判断灰度是否发生较大的变化。如果窗口内的灰度值（在梯度图上）都有较大的变化，那么这个窗口所在区域就存在角点。

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

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

2. 第一步 — 建立数学模型

第一步是通过建立数学模型，确定哪些窗口会引起较大的灰度值变化。让一个窗口的中心位于灰度图像的一个位置(x,y)(x,y)(x,y)，这个位置的像素灰度值为I(x,y)I(x,y)I(x,y) ，如果这个窗口分别向 xxx 和 yyy 方向移动一个小的位移uuu和vvv，到一个新的位置 (x+u,y+v)(x+u,y+v)(x+u,y+v) ，这个位置的像素灰度值就是I(x+u,y+v)I(x+u,y+v)I(x+u,y+v)

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

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

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

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

若窗口内是一个角点，则E(u,v)E(u,v)E(u,v)的计算结果将会很大。

为了提高计算效率，对上述公式进行简化，利用泰勒级数展开来得到这个公式的近似形式：

对于二维的泰勒展开式公式为：
T(x,y)=f(u,v)+(x−u)fx(u,v)+(y−v)fy(u,v)+....T(x,y)=f(u,v)+(x-u)f_x(u,v)+(y-v)f_y(u,v)+....T(x,y)=f(u,v)+(x−u)fx(u,v)+(y−v)fy(u,v)+....

则I(x+u,y+v)I(x+u,y+v)I(x+u,y+v) 为：
I(x+u,y+v)=I(x,y)+uIx+vIyI(x+u,y+v)=I(x,y)+uI_x+vI_yI(x+u,y+v)=I(x,y)+uIx+vIy

其中IxI_xIx和IyI_yIy是III的微分（偏导），在图像中就是求xxx 和 yyy 方向的梯度图：

Ix=∂I(x,y)∂xI_x=\frac{\partial I(x,y)}{\partial x}Ix=∂x∂I(x,y)

Iy=∂I(x,y)∂yI_y=\frac{\partial I(x,y)}{\partial y}Iy=∂y∂I(x,y)

将I(x+u,y+v)=I(x,y)+uIx+vIyI(x+u,y+v)=I(x,y)+uI_x+vI_yI(x+u,y+v)=I(x,y)+uIx+vIy代入E(u，v)E(u，v)E(u，v)可得：

E(u,v)=∑(x,y)w(x,y)×[I(x,y)+uIx+vIy−I(x,y)]2=∑(x,y)w(x,y)×(uIx+vIy)2=∑(x,y)w(x,y)×(u2Ix2+v2Iy2+2uvIxIy)\begin{aligned} E(u, v) &=\sum_{(x, y)} w(x, y) \times\left[I(x, y)+u I_{x}+v I_{y}-I(x, y)\right]^{2} \\ &=\sum_{(x, y)} w(x, y) \times\left(u I_{x}+v I_{y}\right)^{2} \\ &=\sum_{(x, y)} w(x, y) \times\left(u^{2} I_{x}^{2}+v^{2} I_{y}^{2}+2 u v I_{x} I_{y}\right) \end{aligned}E(u,v)=(x,y)∑w(x,y)×[I(x,y)+uIx+vIy−I(x,y)]2=(x,y)∑w(x,y)×(uIx+vIy)2=(x,y)∑w(x,y)×(u2Ix2+v2Iy2+2uvIxIy)

提出 u 和 v ，得到最终二次型的近似形式：

E(u,v)≈[u,v]M[uv]E(u, v) \approx[u, v] M\left[\begin{array}{l}u \\ v\end{array}\right]E(u,v)≈[u,v]M[uv]

其中矩阵M为：M=∑(x,y)w(x,y)[Ix2IxIyIxIyIy2]→R−1[λ100λ2]RM=\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] \rightarrow R^{-1}\left[\begin{array}{cc}\lambda_{1} & 0 \\ 0 & \lambda_{2}\end{array}\right] RM=∑(x,y)w(x,y)[Ix2IxIyIxIyIy2]→R−1[λ100λ2]R

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

（接下来的一小段内容是说：该自相似函数可以视为一个椭圆，而椭圆的形状与特征值的大小相关）
而对于
c(s,y;Δx,Δy)≈AΔx2+2CΔxΔy+BΔy2,c(s, y ; \Delta x, \Delta y) \approx A \Delta x^{2}+2 C \Delta x \Delta y+B \Delta y^{2},c(s,y;Δx,Δy)≈AΔx2+2CΔxΔy+BΔy2, 如果把c视为一个常数, 那么这就是一个二次
型方程，而二次型方程在二维图像中就是一个杨圆。而常数的具体取值并不影响梅圆的形状。而梅圆的形状又恰好是我们要重点研究的，**准确的说是它的长短轴与我们想要求得特征值密切相关。**因此在此我们姑且将该二次型方程视为：

图来自二次型的意义是什么？有什么应用？ - 马同学的回答 - 知乎非常赞的回答！推荐阅读。

（接下里这段是说，求特征值为什么借助矩阵，且特征值的意义是什么）

首先，求特征值为什么借助矩阵，这个问题我不该问的，因为求特征值本身就需要借助矩阵运算。
其次，特征值或者确切说特征矩阵意味这什么？首先，矩阵代表着运动，且当矩阵维度没有发生改变时，一个矩阵就包含了两个运动：旋转和拉伸。而如果想把这两个运动分解，可以对矩阵进行正交化操作和特征值分解，由此得到两个矩阵正交矩阵和对角矩阵。

而正交矩阵就单纯的意味着旋转运动，而对角矩阵就单纯地意味着左右上下拉伸运动。
而对于一个二次方程的二次型矩阵即M(x,y)更换成它对应的对角矩阵，则一个倾斜的椭圆就扶正了。
结合下图进行理解上一段话，如果没能理解，那一定不是你的原因，一定是我没说清楚，具体可以参考文章：二次型的意义是什么？有什么应用？ - 马同学的回答 - 知乎。

经对角化处理后，将两个正交方向的变化分量提取出来，就是 λ1 和 λ2（特征值）。这里利用了线性代数中的实对称矩阵对角化的相关知识，有兴趣的同学可以进一步查阅相关资料。

3. 第二步—角点响应函数R
现在我们已经得到 E(u,v)E(u,v)E(u,v)的最终形式，别忘了我们的目的是要找到会引起较大的灰度值变化的那些窗口。

灰度值变化的大小则取决于矩阵M，M为梯度的协方差矩阵。在实际应用中为了能够应用更好的编程，所以定义了角点响应函数R，通过判定R大小来判断像素是否为角点。

计算每个窗口对应的得分（角点响应函数R定义）：R=det⁡(M)−k(trace⁡(M))2R=\operatorname{det}(M)-k(\operatorname{trace}(M))^{2}R=det(M)−k(trace(M))2

其中 det(M)=λ1λ2det(M)=\lambda_1\lambda_2det(M)=λ1λ2是矩阵的行列式， trace(M)=λ1+λ2trace(M)=\lambda_1+\lambda_2trace(M)=λ1+λ2 是矩阵的迹。

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

RRR的值取决于MMM的特征值，对于角点∣R∣|R|∣R∣很大，平坦的区域∣R∣|R|∣R∣很小，边缘的RRR为负值。

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

注意：Harris 检测器具有旋转不变性，但不具有尺度不变性，也就是说尺度变化可能会导致角点变为边缘。想要尺度不变特性的话，可以关注SIFT特征。

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

平面:：该窗口在平坦区域上滑动，窗口内的灰度值基本不会发生变化，所以 ∣R∣|R|∣R∣ 值非常小，在水平和竖直方向的变化量均较小，即 IxI_xIx和 IyI_yIy都较小，那么 λ1 和 λ2 都较小；

边缘：∣R∣|R|∣R∣值为负数，仅在水平或竖直方向有较大的变化量，即 IxI_xIx和 IyI_yIy只有一个较大，也就是 λ1>>λ2 或 λ2>>λ1；

角点：∣R∣|R|∣R∣ 值很大，在水平、竖直两个方向上变化均较大的点，即 IxI_xIx和 IyI_yIy 都较大，也就是 λ1 和 λ2 都很大。

如下图所示：

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

1.3 基于OpenCV的实现

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

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

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

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

参数解释：

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

import cv2 as cv
from matplotlib import pyplot as plt
import numpy as np# detector parameters
block_size = 3
sobel_size = 3
k = 0.06image = cv.imread('E:/python-project/deep-learning/picture/test1.jpg')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 = cv.cvtColor(image, cv.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 = cv.cornerHarris(gray_img, block_size, sobel_size, k)# result is dilated for marking the corners, not necessary
kernel = cv.getStructuringElement(cv.MORPH_RECT,(3, 3))
dst = cv.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 = cv.cvtColor(image, cv.COLOR_BGR2RGB)
plt.imshow(image)
plt.show()

结果：

1.3 总结

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

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

论文：《C.Harris, M.Stephens. “A Combined Corner and Edge Detector”. Proc. of 4th Alvey Vision Conference》