# Harris 兴趣点检测

# 什么是 “兴趣点”

在 Harris 兴趣点检测中，“兴趣点” 通常指图像中那些在水平和垂直两个方向上都具有显著灰度变化的像素点，也就是说在局部窗口在任意方向做小幅平移都会引起较大像素差异的点。这类点往往对应于图像中的角点或局部纹理显著的区域，因此被认为是具有判别性和稳定性的特征点。

一大片平坦区域，所有方向都变化小，不是兴趣点
高楼顶部的转角处、十字路口的交叉点，在所有方向都变化大这是兴趣点

角点会在两个方向同时产生剧烈变化，边界只会在一个方向上产生剧烈变化，平坦区则都不会。

# 数学原理

# 从上帝视角理解

如何检测本质在于计算窗口滑动前后，其像素的变化情况。

假设一个窗口包含四个像素，那么这里会采集到这四个值：

$x_1,x_2,x_3,x_4$

再次假设，白色部分数值为 0，蓝色部分数值为 1，那么有：

$x_1=0,x_2=1,x_3=0,x_4=1$

水平方向移动：

此时有：

$x_1=0,x_2=0,x_3=0,x_4=0$

数值产生剧烈变化。

垂直方向移动：

$$ x_1=0,x_2=1,x_3=0,x_4=1 $$ 数值不变。

噢～这是边界！

如果把这个放在蓝色矩形的左上角呢？可以很轻松的想象到，无论是在水平方向还是垂直方向，数值都会产生剧烈变化。噢～这是角点！

如果把这个放在蓝色矩形的中间呢？可以更很轻松的想象到，无论是在水平方向还是垂直方向，数值都不会产生剧烈变化。噢～这是平坦区！

我们一看就知道，噢，数值产生剧烈变化了，再思考一下我们就知道这是角点还是边界了，那有没有一种数学方法让计算机也知道这是角点还是边界亦或是平坦区呢？

有的兄弟有的，这就是 Harris 兴趣点检测法。

# Harris 检测方法

对于图像 $I(x,y)$ ，我们可以判断在点 $(x,y)$ 处平移 $(\Delta x,\Delta y)$ 后的自相似性。此处 $I(x,y)$ 是一个灰度图像。

有自相似性公式：

$c(x,y,\Delta x,\Delta y)=\sum_{(u,v)\in W(x,y)}W(u,v)(I(u,v)-I(u+\Delta x,v+\Delta y))^2$

其中 $I(u,v)$ 是原来的灰度值， $I(u+\Delta x,v+\Delta y)$ 是经过平移变换后的灰度值，做减法是为了看看平移操作后其与之前的差异是多大的，有时候我们的灰度变化是上升的，而有时候我们的灰度变化是下降的，我们只关心他变化了多少，因此这里做了平方操作，经过平方后不仅化负为正，还进一步强化了前后的差异。

而 $W(u,v)$ 是一个权重设计，在整个窗口中，我们没有必要平等的看待每一个像素（ps：当然也可以平等的看待，此时 $W(u,v)$ 每一项均为同一常数），通常我们使用高斯加权函数，在图像上，最常用的是二维版本：

$w(x, y) = \frac{1}{2\pi \sigma^2} \exp\left(-\frac{x^2 + y^2}{2\sigma^2}\right)$

其中 $(x, y)$ ：相对于窗口中心的坐标， $\sigma$ 控制模糊程度或者说是影响范围，二维权重的三维表示如图所示：

像图中表示的一样，他将会更加重视窗口中心的元素，一个 $3\times 3$ 的高斯矩阵形如：

$\begin{bmatrix} 1 & 2 & 1\\ 2 & 4 & 2\\ 1 & 2 & 1 \end{bmatrix}$

对 $I(u+\Delta x,v+\Delta y)$ 进行泰勒展开：

$I(u+\Delta x,v+\Delta y)=I(u,v)+I_x(u,v)\Delta x+I_y(u,v)\Delta y+O(\Delta x^2,\Delta y^2)$

如此展开到一阶就够用了，丢弃高阶无穷小后，原式近似等于：

$I(u+\Delta x,v+\Delta y)\approx I(u,v)+I_x(u,v)\Delta x+I_y(u,v)\Delta y$

其中 $I_x$ 表示其对 $x$ 的偏导，其中 $I_y$ 表示其对 $y$ 的偏导。

对于原式 $I(u,v)-I(u+\Delta x,v+\Delta y)$ 我们可以观察到，其包含了 $I(u,v)$ ，因此可以直接被消掉，于是我们重写 $c(x,y,\Delta x,\Delta y)$ ，近似得到：

$c(x,y,\Delta x,\Delta y)\approx \sum_{(u,v)\in W(x,y)}W(u,v)(I_x(u,v)\Delta x+I_y(u,v)\Delta y)^2$

简单起见，我们将 $\sum_{(u,v)\in W(x,y)}W(u,v)$ ，收缩至求和符号中记作：\sum_

进一步，我们可以把上式改写成矩阵的形式，设矩阵 $M(x,y)$ 为：

$M(x,y)=\sum_{w}\begin{bmatrix} I_x(x,y)^2 & I_x(x,y)I_y(x,y)\\ I_x(x,y)I_y(x,y) & I_y(x,y)^2 \end{bmatrix}=\begin{bmatrix} \sum_{w} I_x(x,y)^2 & \sum_{w} I_x(x,y)I_y(x,y)\\ \sum_{w} I_x(x,y)I_y(x,y) &\sum_{w} I_y(x,y)^2 \end{bmatrix}$

观察到矩阵的副对角线完全一样，为进一步简化 $M(x,y)$ ，设 $A=I_x(x,y)^2$ 、 $B=I_x(x,y)^2$ 、 $C=I_x(x,y)I_y(x,y)$ ，于是得到：

$M(x,y)=\begin{bmatrix} A & C\\ C & B \end{bmatrix}$

显然， $M(x,y)$ 是一个对称阵，因此得到：

$c(x,y,\Delta x,\Delta y)\approx \begin{bmatrix} \Delta x, \Delta y \end{bmatrix}M(x,y)\begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix}=\begin{bmatrix} \Delta x, \Delta y \end{bmatrix}\begin{bmatrix} A & C\\ C & B \end{bmatrix}\begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix}$

化简可得：

$c(x,y,\Delta x,\Delta y)\approx A\Delta x^2+2C\Delta x\Delta y+B\Delta y^2$

其中 $A=\sum_{w}I_x^2$ ， $B=\sum_{w}I_y^2$ ， $C=\sum_{w}I_x I_y$ 。

这里的 $A\Delta x^2+2C\Delta x\Delta y+B\Delta y^2$ ，非常像椭圆方程，回忆椭圆一般二次型： $Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0$ ，出现交叉项时，其变为非标准形，形如：

我们需要把整个非标准椭圆做标准化变换。

为什么 $A\Delta x^2+2C\Delta x\Delta y+B\Delta y^2$ 并不标准，因为这里存在交叉项 $2C\Delta x\Delta y$ ，如果这一项能被消掉那就变成标准的了，最直接的方法就是让 $C$ 等于 $0$ ，此处的 $C$ 是由 $M(x,y)$ 产生的，观察到 $M$ 矩阵是一个实对称矩阵，那么我们就可以对 $M$ 进行相似对角化，消掉 $C$ ，变标准椭圆。

下面假设已经完成了对矩阵 $M$ 的形似对角化：

$P^{-1}MP\sim \begin{bmatrix} \lambda _1& \\ & \lambda _2 \end{bmatrix}$

其中， $\lambda _1$ ， $\lambda _2$ ，是矩阵 $M$ 的特征值，因此可进一步写做：

$c(x,y,\Delta x,\Delta y)\approx\begin{bmatrix} \Delta x, \Delta y \end{bmatrix}\begin{bmatrix} A & C\\ C & B \end{bmatrix}\begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix}\approx \begin{bmatrix} \Delta x' , \Delta y' \end{bmatrix} \begin{bmatrix} \lambda _1& \\ & \lambda _2 \end{bmatrix} \begin{bmatrix} \Delta x' \\ \Delta y' \end{bmatrix}=\lambda_1 (\Delta x')^2 + \lambda_2 (\Delta y')^2$

于是我们可以用 $\lambda_1,\lambda_2$ 的大小关系区分平坦、边缘亦或是角点。相似对角化不是为了方便计算，而是为了方便理解与分类，实际实现 Harris 的时候，基本不会显式去算特征值、特征向量，对角化是为了看清 $\lambda_1,\lambda_2$ 的几何意义。

像这样的 $\lambda_1 (\Delta x')^2 + \lambda_2 (\Delta y')^2$ ，我们可以把它改变一下，让其与标准式对齐，得到：

$\frac{\Delta x'^2}{(\sqrt{\frac{1}{\lambda _1} } )^2} +\frac{\Delta y'^2}{(\sqrt{\frac{1}{\lambda _2} } )^2}$

太复杂了，简单点写做：

$\frac{\Delta x'^2}{(\lambda_1^{-\frac{1}{2} })^2} +\frac{\Delta y'^2}{(\lambda_2^{-\frac{1}{2} })^2}$

现在有三种情况：

平坦区域： $\lambda_1 \approx 0,\ \lambda_2 \approx 0$ （两个都很小）
边缘区域： $\lambda_1 \gg 0,\ \lambda_2 \approx 0$ 或者 $\lambda_2 \gg 0,\ \lambda_1 \approx 0$ （一大一小）
角点： $\lambda_1 \gg 0,\ \lambda_2 \gg 0$ （两个都很大）

然后，使用 $\det(M)$ 和 $\text{trace}(M)$ 构造角点响应 $R$ 值：

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

其中 $\det(M)$ 是矩阵 $M$ 的行列式， $\text{trace}(M)$ 是矩阵 $M$ 的迹， $\det(M) = \lambda_1\lambda_2$ ， $\text{trace}(M) = \lambda_1 + \lambda_2$ ，这就是 Harris 设计出来的一个角点评分公式。

进一步重述上述三种情况：

平坦区域： $\lambda_1 \approx 0,\ \lambda_2 \approx 0$ （两个都很小）， $\det = \lambda_1\lambda_2 \approx 0$ ， $\text{trace}^2 \approx 0 ⇒ (R \approx 0)$
边缘区域： $\lambda_1 \gg 0,\ \lambda_2 \approx 0$ 或者 $\lambda_2 \gg 0,\ \lambda_1 \approx 0$ （一大一小）， $\det = \lambda_1\lambda_2$ 不算大， $\text{trace}^2 = (\lambda_1+\lambda_2)^2 \approx \lambda_1^2$ 很大，$⇒ R $ 往往是负的（这里 Harris 故意用 $-k \cdot \text {trace}^2 $ 把边缘压下去）
角点： $\lambda_1 \gg 0,\ \lambda_2 \gg 0$ （两个都很大）， $\det = \lambda_1\lambda_2$ 很大，即便减去了 $k(\lambda_1+\lambda_2)^2$ ，但整体还是会是个比较大的正数

综上所述， $R$ 有如下特性：

$R \approx 0$ ⇒ 平坦区
$R < 0$ ⇒ 更像边缘
$R \gg 0$ ⇒ 很像角点（兴趣点）

$\det(M)$ 大，就说明两个方向上变化都大， $\text{trace}(M)^2$ 大，可能只是有一个方向特别大，比如边缘， $R = \det - k \cdot \text{trace}^2$ ，用 $\det$ 奖励两个方向都大的情况，用 $-k\text{trace}^2$ 惩罚只有一个方向大、另一个方向小的情况。

# 在 OpenCV 中的使用

	import cv2
	img = cv2.imread('1.jpg')
	gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
	dst = cv2.cornerHarris(gray, 2, 3, 0.04)
	img[dst > 0.1 * dst.max()]=[0, 0, 255]
	cv2.imshow('dst', img)
	cv2.imwrite('1_corner.png', img)
	cv2.waitKey(0)
	cv2.destroyAllWindows()

blockSize=2 是窗口大小
ksize =3 是 Sobel 求导算子核大小，Harris 需要先算图像在 x，y 方向的梯度 $I_x, I_y$ ，OpenCV 用的是 Sobel 算子。
k = 0.04 是 Harris 公式里的经验参数 k， $R = \det(M) - k \cdot \text{trace}^2(M)$

# 对比

# 简单场景

# 复杂场景

# 参考

[1] Harris C G, Stephens M. A combined corner and edge detector[C] Alvey vision conference. 1988, 15(50): 10-5244.

[2] Harris 角点检测与 SIFT 特征匹配全解析 https://www.bilibili.com/video/BV1Zi1eYLEiR

[3] Harris Corner Detection https://docs.opencv.org/4.x/dc/d0d/tutorial_py_features_harris.html

[4] Harris corner detector https://en.wikipedia.org/wiki/Harris_corner_detector?utm_source=chatgpt.com