原理

Harris 角点检测的基本思路如下：考虑一个局部的区域，将其作为一个窗口四处移动，若窗口灰度发生了较大的变化，那么，就认为窗口内存在角点，否则窗口内就不存在角点。

对于图像$I(\vec{p})$，当在点$\vec{p}=(x,y)$处平移$\Delta \vec{p}=(\Delta x,\Delta y)$后的自相似性，可以通过下面的自相关函数给出

$$S_W(\Delta x,\Delta y)=\sum _{x_i\in W}\sum _{y_i\in W}w({\vec{p}}_i)(I({\vec{p}}_i+\Delta \vec{p})-I({\vec{p}}_i){)}^2$$

其中，$W$是以点$(x,y)$为中心的窗口，$w({\vec{p}}_i)$是点${\vec{p}}_i$处的加权函数，一般取高斯函数$e^{-\frac{(x_i-x{)}^2+(y_i-y{)}^2}{2{\sigma }^2}}$。当$\Delta \vec{p}$取不同值时，若$S_W$都有较大变化，则可认为$\vec{p}$是一个角点，此即Harris角点检测的基本思路。

对平移后的图像进行泰勒展开，可得

$$I({\vec{p}}_i+\Delta \vec{p})\approx I({\vec{p}}_i)+I_x({\vec{p}}_i)\Delta x+I_y({\vec{p}}_i)\Delta y$$

其中$I_x,I_y$为偏导数，将其带入$S_W$，有

$$\begin{array}{rl}{S}_W(\Delta x,\Delta y)& \approx \sum _{x_i\in W}\sum _{y_i\in W}w({\vec{p}}_i)(I_x({\vec{p}}_i)\Delta x+I_y({\vec{p}}_i)\Delta y{)}^2\\ & =\sum _{x_i\in W}\sum _{y_i\in W}w({\vec{p}}_i)\left[I_x^2({\vec{p}}_i)\Delta x^2+I_y^2({\vec{p}}_i)\Delta y^2+2I_x({\vec{p}}_i)I_y({\vec{p}}_i)\Delta x\Delta y\right]\\ & =\sum _{x_i\in W}\sum _{y_i\in W}w({\vec{p}}_i)\left[\begin{array}{cc}\Delta x& \Delta y\end{array}\right]\left[\begin{array}{cc}{I}_x^2({\vec{p}}_i)& I_x({\vec{p}}_i)I_y({\vec{p}}_i)\\ I_x({\vec{p}}_i)I_y({\vec{p}}_i)& I_y^2({\vec{p}}_i)\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta y\end{array}\right]\end{array}$$

将其记作

$$S_W(\Delta x,\Delta y)=\left[\begin{array}{cc}\Delta x& \Delta y\end{array}\right]M(x,y)\left[\begin{array}{c}\Delta x\\ \Delta y\end{array}\right]$$

其中$M(\vec{p})$只与当前点$\vec{p}$有关，可记作

$$\begin{array}{rl}M(\vec{p})& =\sum _{x_i\in W}\sum _{y_i\in W}w({\vec{p}}_i)\left[\begin{array}{cc}{I}_x^2({\vec{p}}_i)& I_x({\vec{p}}_i)I_y({\vec{p}}_i)\\ I_x({\vec{p}}_i)I_y({\vec{p}}_i)& I_y^2({\vec{p}}_i)\end{array}\right]\\ & =\left[\begin{array}{cc}{\sum }_{x_i\in W}{\sum }_{y_i\in W}w({\vec{p}}_i)I_x^2({\vec{p}}_i)& {\sum }_{x_i\in W}{\sum }_{y_i\in W}w({\vec{p}}_i)I_x({\vec{p}}_i)I_y({\vec{p}}_i)\\ {\sum }_{x_i\in W}{\sum }_{y_i\in W}w({\vec{p}}_i)I_x({\vec{p}}_i)I_y({\vec{p}}_i)& {\sum }_{x_i\in W}{\sum }_{y_i\in W}w({\vec{p}}_i)I_y^2({\vec{p}}_i)\end{array}\right]\end{array}$$

$M(\vec{p})$有四项，可写作

$$M(\vec{p})=\left[\begin{array}{cc}A& C\\ C& B\end{array}\right]$$

从而自相关函数可近似为

$$\begin{array}{rl}{S}_W(\Delta x,\Delta y)& \approx \left[\begin{array}{cc}\Delta x& \Delta y\end{array}\right]\left[\begin{array}{cc}A& C\\ C& B\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta y\end{array}\right]\\ & =A\Delta x^2+2C\Delta x\Delta y+B\Delta y^2\end{array}$$

以$\Delta x,\Delta y$分别作为横轴与纵轴，则$A,B,C,S_W$固定时，上式构成一个椭圆方程，其短轴和长轴分别为$\sqrt{\frac{S_W}{{\lambda }_M}}，\sqrt{\frac{S_W}{{\lambda }_m}}$，${\lambda }_M,{\lambda }_m$为$M(\vec{p})$的大小两个特征值。

参数一经固定，那么$\Delta x,\Delta y$就只能在椭圆上取值，才能使等式成立。而若换一种思路，点$\Delta \vec{p}=(\Delta x,\Delta y)$，将存在三种位置，即椭圆内部、椭圆上和椭圆外部。对于角点来说，应该让$\Delta \vec{p}$挪动尽量小的情况下，从椭圆内部跳到椭圆外部，换言之，小窗移动因其较大变化。因此，椭圆的轴越短，则这种变化越明显。从而，当矩阵$M(x)$的特征值较大时，可以认为$\vec{p}=(x,y)$是角点。

由于求特征值比较耗时，Harris定义了响应值$R$

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

其中trace表示迹，det为行列式，$k$为常数，一般取$0.04\to 0.06$。$R$值与点的关系如下

• $|R|$很小时，两个特征值都很小，对应区域为平面。
• $R<0$时，两个特征值相差较大，对应区域为直线
• $R$较大时，两个特征值都很大，且近似相等，对应区域为角点。

测试

python-opencv中提供了Harris角点检测函数，其$R$值计算结果以及角点检测结果如下

代码为

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

path = 'lena.jpg'
img = plt.imread(path)
gray = img[:,:,0]

harris = cv2.cornerHarris(gray, 2, 3, 0.04)

fig = plt.figure()
ax = fig.add_subplot(121)
ax.imshow(harris)

ax = fig.add_subplot(122)
dst = cv2.dilate(harris, None) # 腐蚀harris结果
im1 = img + 0

im1[dst > 0.1 * dst.max()] = [255, 0, 0]
ax.imshow(im1)
plt.tight_layout()
plt.show()

【cornerHarris】即为harris角点计算函数，其输入的四个参数分别是待检测图像；角点检测时的移动范围；Sobel导数的尺寸（为奇数）以及$k$值。