CV笔记8：角点特征之Harris角点检测算法（基于python-opencv实现）

---

一、背景介绍

1.1 什么是特征

在数据挖掘、计算机视觉等领域经常会提起“特征”一词，那么究竟什么才是特征？数据挖掘中我们可以称被分析数据集的某一属性（一行为一个样本，一列为一个属性）为特征，因为它反映了该数据集的某一种特性，是该数据集的一种固有标识。在计算机视觉中，一张图像的特征其实也有类似的概念，这种特征也能够描述这张图像的全部（或局部）内容或含义。图像特征可以是①有意义的图像区域，该区域具有独特的特征和良好的识别性；也可以是②对图像信息的一种数据抽取和表示。
针对①可以是图像中的角点、边缘、斑点、直线、曲线及高密度区域等。大多数特征检测算法都会涉及图像的角点、边和斑点的识别，也有一些涉及脊向的概念，可以认为脊向是细长物体的对称轴，例如识别图像中的一条路。角点和边都好理解，那什么是斑点呢？斑点通常是指与周围有着颜色和灰度差别的区域。在实际地图中，往往存在着大量这样的斑点，如一颗树是一个斑点，一块草地是一个斑点，一栋房子也可以是一个斑点。由于斑点代表的是一个区域，相比单纯的角点，它的稳定性要好，抗噪声能力要强，所以它在图像配准上扮演了很重要的角色。
针对上述的②，在文章CV笔记7：计算机视觉的通俗理解 第二节中较通俗的解释了对图像特征的理解。

1.2 为什么做特征检测

图像的特征应用很广泛，我们通过图像特征能过够将实现图像的分类、目标检测、图像拼接等。而特征就需要特征检测来提出。我们直接使用下面的例子进行说明。
这里有两张同一山峰不同部位的图像，我们怎么能够将上面的两张图像拼接成一张图像呢？对于我们人来说，将两张图像相同的位置重合折叠在一起就可以了，就像下图，但是对于计算机要怎么实现？

对于计算机来说，并不能像人一样直接对图像进行拼接，而是要根据一些关键点对图像进行处理。我们可以通过一些角点来让计算机认识并拼接图像，就像下面图片中画出的一样，我们找到图像中的一些关键点（角点），并让两张图像中相同的角点处进行拼接。这样计算机也能够像人一样实现图像的拼接了。

上面的例子只是用到了角点特征，其实其他的类似特征，如边缘、线等都能用在不同的场景中。所以，最终达到的目标是：通过特征检测我们实现的就是能够让计算机“读懂”图像，将人类看到的视觉信息转化成计算机能够识别和处理的定量形式，这也就是一种图像特征提取方式。

1.3 角点特征

1.3.1 角点定义

我们在上面图像拼接的说明案例中提到了角点特征，本文也重点介绍角点特征的提取过程。角点一般多为图像轮廓之间的交点，对于同一场景，即使视角发生变化，也通常具备稳定性质的特征，在角点附近区域的像素点无论在梯度方向上还是其梯度幅值上有着较大变化。可以看出，在对图像求导后，极值点处往往就是角点的位置。

再通过图像理解，就是在一幅图像中，如上图，我们考虑图像上的一个很小的窗口。在平坦的区域时，当小窗口进行任意方向的移动，窗口内的像素都没有太大灰度变化；在边缘区域时，当小窗口沿着边缘方向移动，窗口内的像素没有太大灰度变化，当沿着垂直边缘方向移动，窗口内的像素会发生跳变；在角点区域时，当小窗口沿任意方向移动，窗口内的像素都有明显的灰度变化。下图展示了不同类型的角点。

1.3.2 角点优势

角点相对于其他的图像特征具有以下的优势：

• 点特征属于局部特征，对遮挡有一定的鲁棒性
• 通常图像中可以检测到成百上千的点特征，以量取胜
• 点特征辨识度好，不同物体上的点容易区分
• 点特征提取速度很快

二、Harris角点检测原理

2.1 算法思想

我们已经对角点有了初步的认识，下一步开始进行更深的数理推导，也就是Harris角点检测原理。算法的核心是利用局部窗口在图像上进行移动，判断灰度是否发生较大的变化。如果窗口内的灰度值（在梯度图上）都有较大的变化，那么这个窗口所在区域就存在角点。
这样就可以将 Harris 角点检测算法分为以下三步：

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

2.2 算法推导

我们仍然考虑图像中的一个小窗口，如下图：

2.2.1 灰度变化描述

当窗口发生$[u,v]$移动时，那么滑动前与滑动后对应的窗口中的像素点灰度变化描述如下：
$$E(u,v)=\sum _{(x,y)€W}w(x,y)[I(x+u,y+v)-I(x,y){]}^2$$参数解释：

• $[u,v]$是窗口$W$的偏移量；
• $(x,y)$是窗口$W$所对应的像素坐标位置，窗口有多大，就有多少个位置；
• $I(x,y)$是像素坐标位置$(x,y)$的图像灰度值；
• $I(x+u,y+v)$是像素坐标位置$(x+u,y+v)$的图像灰度值；
• $w(x,y)$是窗口函数，最简单情形就是窗口$W$内的所有像素所对应的$w$权重系数均为1.但有时候，我们会将$w(x,y)$函数设置为以窗口$W$中心为原点的二元正太分布。如果窗口W中心点是角点时，移动前与移动后，该点在灰度变化贡献最大；而离窗口$W$中心(角点)较远的点，这些点的灰度变化几近平缓，这些点的权重系数，可以设定小值，以示该点对灰度变化贡献较小，那么我们自然而然想到使用二元高斯函数来表示窗口函数；
  
  根据上述表达式，当窗口在平坦区域上移动，可以想象得到，灰度不会发生太大变化。$E(u,v)=0$；如果窗口处在纹理比较丰富的区域上滑动，那么灰度变化会很大。算法最终思想就是计算灰度发生较大变化时所对应的位置，当然这个较大是指任意方向上的滑动，并非单指某个方向。
  同时，我们应该注意：$E(u,v)$计算的是窗口$W$向$[u,v]$这一个方向移动所产生的灰度变化。

2.2.2 $E(u,v)$公式化简

为了提高计算效率，利用泰勒级数展开对上述公式进行简化。
对于二维的泰勒展开式公式为：
$$f(x+u,y+v)\approx f(x,y)+u{f}_x(x,y)+v{f}_y(x,y)$$则对于$I(x+u,y+v)$有：
$$I(x+u,y+v)=I(x,y)+u{I}_x+v{I}_y$$其中$I_x$和$I_y$是$I$的微分（偏导），在图像中就是求$x$ 和 $y$ 方向的梯度：
$$I_x=\frac{\partial I(x,y)}{\partial x}$$ $$I_y=\frac{\partial I(x,y)}{\partial y}$$ 那么有：
$$E(u,v)=\sum _{(x,y)€W}w(x,y)[I(x+u,y+v)-I(x,y){]}^2$$ $$\approx \sum _{(x,y)€W}w(x,y)[I(x,y)+u{I}_x+v{I}_y-I(x,y){]}^2$$ $$=\sum _{(x,y)€W}w(x,y)[u^2{I}_x^2+2uv{I}_x{I}_y+v^2{I}_y^2]$$ $$=\sum _{(x,y)€W}w(x,y)[u,v]\left[\begin{array}{cc}{I}_x^2& I_{x}Iy\\ I_y{I}_x& I_y^2\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]$$ $$=[u,v](\sum _{(x,y)€W}w(x,y)\left[\begin{array}{cc}{I}_x^2& I_{x}Iy\\ I_y{I}_x& I_y^2\end{array}\right])\left[\begin{array}{c}u\\ v\end{array}\right]$$从而可以得到：
$$E(u,v)=[u,v]M\left[\begin{array}{c}u\\ v\end{array}\right]$$其中：
$$M=\sum _{(x,y)€W}w(x,y)\left[\begin{array}{cc}{I}_x^2& I_{x}Iy\\ I_y{I}_x& I_y^2\end{array}\right]⟹\left[\begin{array}{cc}A& C\\ C& B\end{array}\right]⟹R^{-1}\left[\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right]R$$最后矩阵形式表达是把实对称矩阵对角化处理后的结果，可以把$R$看成旋转因子，其不影响两个正交方向的变化分量。经对角化处理后，将两个正交方向的变化分量提取出来，就是 ${\lambda }_1$ 和 ${\lambda }_2$（特征值）。

2.2.3 矩阵$M$的理解

虽然我们已经得到了$E(u,v)$，但我们并不直接使用它来进行判断当前窗口是否含有角点。Harris角点检测而是通过对窗口内的每个像素的$x$方向上的梯度与$y$方向上的梯度进行统计分析，并结合了矩阵$M$的性质进行角点判定的。
我们以$I_x$和$I_y$为坐标轴，每个像素的梯度坐标可以表示成$(I_x,I_y)$，在此基础上针对平坦区域，边缘区域以及角点区域和斜边缘区域四种情形进行分析：

针对上图中四种窗口区域，统计对应像素的梯度分布情况，得到下图，其中$I_x$为横轴$I_y$为纵轴：

我们从上面能够观察到这几种区域的特点：

• 平坦区域上的每个像素点所对应的$(I_x,I_y)$坐标分布在原点附近，这是因为平坦区域像素梯度方向虽然各异，但是其幅值都不是很大，所以均聚集在原点附近
• 边缘区域沿一个方向分布较散，至于是哪一个方向不能一概而论，这要视边缘在图像上的具体位置而定，如果边缘是水平或者垂直方向，那么$I_y$轴方向或者$I_x$方向上的数据分布就比较散
• 角点区域的在$x$、$y$方向上的梯度分布都比较散
  我们再回头看$E(u,v)$，根据上面的推导，$M$可以表示为如下的形式:$$\left[\begin{array}{cc}A& C\\ C& B\end{array}\right]$$所以，我们可以将$E(u,v)$近似为二项函数：
  $$E(u,v)=A{u}^2+2Cuv+B{v}^2$$ 其中：
  $$A=\sum _{(x,y)€W}w(x,y)\ast I_x^2$$$$B=\sum _{(x,y)€W}w(x,y)\ast I_y^2$$$$C=\sum _{(x,y)€W}w(x,y)\ast I_x{I}_y$$二次项函数本质上就是一个椭圆函数。椭圆的扁率和尺寸是由其特征值λ1、λ2决定的,方向是由其特征向量决定的。以下图为例，设椭圆方程为：
  $$[u,v]M\left[\begin{array}{c}u\\ v\end{array}\right]=1$$
  可以看出$M$的特征根决定了椭圆的长短轴的长度，对应的特征向量决定了椭圆的方向（因为椭圆的两个轴指向特征向量的方向）。所以，我们知道了求得的$E(u,v)$是可以通过椭圆的形式来表达的，我们对之前的数据集用椭圆形式表示，绘制的图像如下图所示：
  
  之所以能够使用椭圆描述上面的数据集，是因为矩阵$M$本身的形式和协方差矩阵就有着千丝万缕的联系。我们将$I_x,I_y$看成两个字段，假设窗口内有$m$个像素点，也就是等价于有$m$个样本，我们先计算每个字段的均值：
  $$\overline{I_x}=\sum _{i=1}^m{I}_{xi}$$ $$\overline{I_y}=\sum _{i=1}^m{I}_{yi}$$ 我们仍然使用$(I_{xi},I_{yi})$表示样本$(I_{xi},I_{yi})$去均值后的值，则由这$m$个样本组成的矩阵为：
  $$X=\left[\begin{array}{cccc}{I}_{x1}& I_{x2}& ...& I_{xm}\\ I_{y1}& I_{y2}& ...& I_{ym}\end{array}\right]$$则对应的协方差矩阵为：
  $$C=\frac{1}{m}X{X}^T=\frac{1}{m}\sum _{i=1}^m\left[\begin{array}{cc}{I}_x^2& I_{x}Iy\\ I_y{I}_x& I_y^2\end{array}\right]$$我们可以看到，上式中的$C$中，我们先进行了各维的零均值化，这样各维所对应的随机变量的均值为0，协方差矩阵就大大简化，简化的最终结果就是矩阵$M$(注意：在这里为了简化运算，我们先假设了$M$矩阵中的权重系数$w(x,y)=1$,并且忽略了最后除$m$样本数的操作)。到这里是否已经明白了我们的目的：我们是来分析图像导数数据的主要成分。
  先前我们已经对$M$进行了对角化操作：
  $$M=\sum _{(x,y)€W}w(x,y)\left[\begin{array}{cc}{I}_x^2& I_{x}Iy\\ I_y{I}_x& I_y^2\end{array}\right]⟹R^{-1}\left[\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right]R$$是否也让你想起了PCA中的操作？不明白可以复习一下PCA（注：协方差矩阵需要大家深刻理解一下）。所以再结合上面数据集分布图像，可以知道：
• 如果两个字段$(I_x,I_y)$所对应的特征值都比较大，说明像素点的梯度分布比较散，梯度变化程度比较大，符合角点在窗口区域的特点；
• 如果是平坦区域，那么像素点的梯度所构成的点集比较集中在原点附近，因为窗口区域内的像素点的梯度幅值非常小，此时矩阵M的对角化的两个特征值比较小；
• 如果是边缘区域，在计算像素点的$x、y$方向上的梯度时，边缘上的像素点的某个方向的梯度幅值变化比较明显，另一个方向上的梯度幅值变化较弱，其余部分的点都还是集中原点附近，这样$M$对角化后的两个特征值理论应该是一个比较大，一个比较小，当然对于边缘，某些情况下致使计算出的特征值并不是都特别的大，但仍跟含有角点的窗口的分布情况具有不同。
  
  所以我们可以直接根据$M$两个特征值的大小对图像点进行分类，如上图所示：
• 特征值都比较大时，即窗口中含有角点；
• 特征值一个较大，一个较小，窗口中含有边缘；
• 特征值都比较小，窗口处在平坦区域；

2.2.4 角点响应函数$R$

我们已经知道了什么样的窗口含有角点又或者是边缘等。在实际应用中，为了更好的应用到编程中，我们有定义了角点响应函数$R$，通过判定$R$大小来判断窗口是否有角点，即$R>threshold$时则有角点，$threshold$是我们自定义的一个阈值。
最简单的一种是，角点应该满足基本性质：窗口最小的特征值尽量的大。此时定义$R={\lambda }_{min}$。
还有比上式更有效的角点相应函数：
$$R=det(M)-k(trace(M){)}^{2}k\in [0.04,0.06](1)$$ $$R={\lambda }_1-k{\lambda }_{2}k=0.05(2)$$ $$R=\frac{det(M)}{trace(M)}=\frac{{\lambda }_1{\lambda }_2}{{\lambda }_1+{\lambda }_2}(3)$$上面新给出了三种角点响应函数，我们以(1)进行讲解：
$$R=det(M)-k(trace(M){)}^{2}k\in [0.04,0.06]$$ $$det(M)={\lambda }_1{\lambda }_2$$ $$trace(M)={\lambda }_1+{\lambda }_2$$ 这里${\lambda }_1,{\lambda }_2$是矩阵$M$的2个特征值，$k$是一个指定值，这是一个经验参数，需要实验确定它的合适大小，通常它的值在0.04和0.06之间,它的存在只是调节函数的形状而已。注：一般的，增大k的值，降低角点检测的灵敏度，减少被检测角点的数量；减少k值，增加角点检测的灵敏度，增加被检测角点的数量。
为什么可以使用这个函数进行判定呢？我们将其图像画出进行理解，图像如下，可以看出这个函数图形正好满足角点、边缘和平坦区域的特征。

之后，我们根据需求设定$R$的阈值，就可以进行角点的判断了。后续还可以增加很多处理，比如如果需要可以在3×3或者5×5的邻域进行非最大值抑制，则局部最大值点即为图像中的角点。

2.3 Harris角点的算法流程及性质

算法流程：

• 将原图像$I$使用$w(x,y)$进行卷积，并计算图像的梯度$I_x$和$I_y$;
• 计算每一个图像像素点的自相关矩阵$M$;
• 计算角点相应$R$;
• 选择$R$大于某一阈值的点作为角点；
• 根据需要在图像区域内进行角点的非极大值抑制；

Harris角点缺点：
Harris角点检测获取的角点在图像中分布不均匀（对比度高的区域角点多）

Harris角点性质：

• 旋转不变性
  Harris角点检测算子使用的是角点附近的区域灰度二阶矩矩阵。而二阶矩矩阵可以表示成一个椭圆，椭圆的长短轴正是二阶矩矩阵特征值平方根的倒数。当特征椭圆转动时，特征值并不发生变化，所以判断角点响应值R也不发生变化，由此说明Harris角点检测算子具有旋转不变性。
• 光照不变性、比度变化部分不变性
  这是因为在进行Harris角点检测时，使用了微分算子对图像进行微分运算，而微分运算对图像密度的拉升或收缩和对亮度的抬高或下降不敏感。换言之，对亮度和对比度的仿射变换并不改变Harris响应的极值点出现的位置，但是，由于阈值的选择，可能会影响角点检测的数量。
  
• 不具有尺度不变性
  如下图所示，当图像被缩小时，在检测窗口尺寸不变的前提下，在窗口内所包含图像的内容是完全不同的。左侧的图像可能被检测为边缘或曲线，而右侧的图像则可能被检测为一个角点。
  

三、基于python-opencv实现Harris角点检测

3.1 opencv API实现

python-opencv提供了Harris角点检测的函数:

cornerHarris(src, blockSize, ksize, k[, dst[, borderType]]) -> dst
.
The function runs the Harris corner detector on the image. Similarly to cornerMinEigenVal and ornerEigenValsAndVecs , for each pixel $(x,y)$ it calculates a $2\times 2$ gradient covariance matrix $M^{(x,y)}$ over a $\text{blockSize}\times \text{blockSize}$ neighborhood. Then, it computes the following characteristic:
. $\text{dst}(x,y)=\mathrm{d}\mathrm{e}\mathrm{t}{M}^{(x,y)}-k\cdot {\left(\mathrm{t}\mathrm{r}{M}^{(x,y)}\right)}^2$
. Corners in the image can be found as the local maxima of this response map.
.
. @param src Input single-channel 8-bit or floating-point image. 输入单通道8位或者浮点型图像
. @param dst Image to store the Harris detector responses. It has the type CV_32FC1 and the same size as src .输出为角点响应图
. @param blockSize Neighborhood size (see the details on #cornerEigenValsAndVecs ). 扫描时窗口大小
. @param ksize Aperture parameter for the Sobel operator. 使用Sobel算子，该参数定义了Sobel算子的中孔。简单来说，该函数定义了角点检测的敏感度，其值必须介于3~31之间的奇数。
. @param k Harris detector free parameter. See the formula below. 响应函数中的k值，一般取0.04~0.06
. @param borderType Pixel extrapolation method. See #BorderTypes. 像素插值方法

import numpy as np
import cv2

# original image
image = cv2.imread('./timg.jpg')
h, w, c = image.shape
print('image shape --> h:%d  w:%d  c:%d' % (h, w, c))
cv2.imshow('image', image)
cv2.waitKey(2000)
cv2.destroyAllWindows()

# harris dst
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
# gray = np.float32(gray)
dst = cv2.cornerHarris(gray, blockSize=3, ksize=5, k=0.05)
image_dst = image[:, :, :]
image_dst[dst > 0.01 * dst.max()] = [0, 0, 255]
cv2.imwrite('./dst.jpg', image_dst)
cv2.imshow('dst', image_dst)
cv2.waitKey(0)
cv2.destroyAllWindows()

处理图像对比：

3.2 python实现

def harris_det(img, block_size=3, ksize=3, k=0.04, threshold = 0.01, WITH_NMS = False):
    '''
    params:
        img:单通道灰度图片
        block_size:权重滑动窗口
        ksize：Sobel算子窗口大小
        k:响应函数参数k
        threshold:设定阈值
        WITH_NMS:非极大值抑制
    return：
        corner：角点位置图，与源图像一样大小，角点处像素值设置为255
    '''
    h, w = img.shape[:2]
    # 1.高斯权重
    gray = cv2.GaussianBlur(img, ksize=(ksize, ksize), sigmaX=2)

    # 2.计算梯度
    grad = np.zeros((h,w,2),dtype=np.float32)
    grad[:,:,0] = cv2.Sobel(gray,cv2.CV_16S,1,0,ksize=3)
    grad[:,:,1] = cv2.Sobel(gray,cv2.CV_16S,0,1,ksize=3)

    # 3.计算协方差矩阵
    m = np.zeros((h,w,3),dtype=np.float32)
    m[:,:,0] = grad[:,:,0]**2
    m[:,:,1] = grad[:,:,1]**2
    m[:,:,2] = grad[:,:,0]*grad[:,:,1]
    m = [np.array([[m[i,j,0],m[i,j,2]],[m[i,j,2],m[i,j,1]]]) for i in range(h) for j in range(w)]

    # 4.计算局部特征结果矩阵M的特征值和响应函数R(i,j)=det(M)-k(trace(M))^2  0.04<=k<=0.06
    D,T = list(map(np.linalg.det,m)),list(map(np.trace,m))
    R = np.array([d-k*t**2 for d,t in zip(D,T)])

    # 5.将计算出响应函数的值R进行非极大值抑制，滤除一些不是角点的点，同时要满足大于设定的阈值
    #获取最大的R值
    R_max = np.max(R)
    #print(R_max)
    #print(np.min(R))
    R = R.reshape(h,w)
    corner = np.zeros_like(R,dtype=np.uint8)
    for i in range(h):
        for j in range(w):
            if WITH_NMS:
                #除了进行进行阈值检测 还对3x3邻域内非极大值进行抑制(导致角点很小，会看不清)
                if R[i,j] > R_max*threshold and R[i,j] == np.max(R[max(0,i-1):min(i+2,h-1),max(0,j-1):min(j+2,w-1)]):
                    corner[i,j] = 255
            else:
                #只进行阈值检测
                if R[i,j] > R_max*threshold :
                    corner[i,j] = 255
    return corner

if __name__ == "__main__":
    image = cv2.imread('./timg.jpg')
    height, width, channel = image.shape
    print('image shape --> h:%d  w:%d  c:%d' % (height, width, channel))
    cv2.imshow('mount', image)
    cv2.waitKey(2000)
    cv2.destroyAllWindows()

    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
    # gray = np.float32(gray)
    dst = harris_detect(gray)

    image_dst = image[:, :, :]
    image_dst[dst > 0.01 * dst.max()] = [0, 0, 255]
    cv2.imwrite('./dsti8_1.jpg', image_dst)
    cv2.imshow('dst', dst)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

测试结果：

参考

第十一节、Harris角点检测原理(附源码)
PCA的数学原理
Harris角点