Harris角点检测

特征检测中一般都会涉及角点、边、斑点，而Harris是一种常用的角点检测的方法。

角点是什么呢？下面看一幅图来理解：

上图中标红的点即角点，可以很清楚的看出，物体拐角处的点就是角点（字面上的理解）。

Harris角点检测的思想：使用一个固定滑动窗口在图像上的任意方向进行滑动，比较窗口中的像素灰度变化程度，如果在任意方向上都有较大灰度值变化，那么可以认为该窗口中存在角点。

上面的基本思想可能有点抽象，下面看一组图辅助理解：

根据上图来看：
（1）对于左图，滑动窗口往任意方向滑动，都不会有较大的灰度值变化，因此滑动窗口对应的区域为平坦区域；
（2）对于中图，滑动窗口沿着边线移动，无灰度变化；而其他方向可以有灰度变化，因此属于边缘部分；
（3）对于右图，滑动窗口无论往哪个方向移动，对会带来较大的灰度变化，因此滑窗对应的区域有角点。

下面是Harris角点检测的数学原理：

1.灰度变化描述

当窗口发生[u,v]移动时，那么滑动前与滑动后对应的窗口中的像素点灰度变化描述如下：

$E(u,v)={\sum }_{x,y\in W}w(x,y)[I(x+u,y+v)-I(x,y){]}^2$

其中$[u,v]$是窗口$W$的偏移量，$(x,y)$是窗口所对应的像素点位置，$I(x,y)$是对应位置的像素值，$w(x,y)$这里一般选择的是高斯
加权函数，即离中心点近的梯度贡献更大，距离中心远的点贡献小。

通过上式可以看出，在平坦区域的灰度变化不大，$E(u,v)$很小；而在纹理丰富的区域，$E(u,v)$会比较大。

2. E(u,v)化简

首先对原始的$E(u,v)$进行泰勒展开（省略掉无穷小项），得到：

$E(u,v)={\sum }_{x,y\in W}w(x,y)[I(x,y)+u{I}_x+v{I}_y-I(x,y){]}^2$

$={\sum }_{x,y\in W}w(x,y)(u^2{I}_x^2+2uv{I}_x{I}_y+v^2{I}_y^2)$

$={\sum }_{x,y\in W}w(x,y)\left[\begin{array}{c}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}{cc}u& v\end{array}\right]$

$=\left[\begin{array}{c}u\\ v\end{array}\right]\left({\sum }_{x,y\in W}w(x,y)\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}{cc}u& v\end{array}\right]$

因此，$E(u,v)$可以更新为：

$E(u,v)=\left[\begin{array}{c}u\\ v\end{array}\right]M\left[\begin{array}{cc}u& v\end{array}\right]$

其中$M={\sum }_{x,y\in W}w(x,y)\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]$

$I_x$和$I_y$是滑动窗口内的像素点$(x,y)$在x和y方向上的梯度值。

3. 关键的矩阵M

Harris角点检测并没有直接使用$E(u,v)$来确定角点。而是以x和y方向上的梯度为坐标着，对窗口内每个点的梯度进行统计分析。最终将M的特征值转化为椭圆的两个轴，更详细的原理可以参考这篇博客： https://www.cnblogs.com/zyly/p/9508131.html

得到的结论如下图：

对于矩阵M的特征值${\lambda }_1和{\lambda }_2$：

（1）当${\lambda }_1和{\lambda }_2$都较小时，对应的区域为平坦区域；

（2）当${\lambda }_1或{\lambda }_2$有一个较大，而另一个较小时，对应的区域存在边缘；

（3）当${\lambda }_1和{\lambda }_2$都较大时，对应的区域存在角点。

4. 角点响应的度量

对于每一个滑动窗口计算一个得分R，如果R大于某一个阈值，那么可以认为该滑动窗口对应的区域内存在角点。

$R=det(M)-k(trace(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来说：

（1）当滑窗内存在角点是，R为正数且不会太小；

（2）当滑窗内为平坦区域，R是一个很小的值；

（3）当滑窗内存在边缘，R值为负。

因此可以设定一个合理的阈值，当R大于该阈值时，认为滑窗内存在角点。

Harris角点检测实例

下面的一个例子是检测一副RGB图像的角点。

（1）将RGB三通道图像转化为灰度图；

（2）利用Sobel算子计算灰度图在x和y方向上的梯度矩阵；

（3）对得到的两个梯度矩阵进行高斯滤波，得到新的梯度矩阵；

（4）计算梯度矩阵上每个点的R值（每个点对应之前的一个滑窗），如果R值大于一个设定的阈值，则该点是角点。

该例子的实现源代码如下：

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


# Harris corner detection
def Harris_corner(img):

	## Grayscale
	def BGR2GRAY(img):
		gray = 0.2126 * img[..., 2] + 0.7152 * img[..., 1] + 0.0722 * img[..., 0]
		gray = gray.astype(np.uint8)
		return gray

	## Sobel
	def Sobel_filtering(gray):
		# get shape
		H, W = gray.shape

		# sobel kernel
		sobely = np.array(((1, 2, 1),
						(0, 0, 0),
						(-1, -2, -1)), dtype=np.float32)

		sobelx = np.array(((1, 0, -1),
						(2, 0, -2),
						(1, 0, -1)), dtype=np.float32)

		# padding
		tmp = np.pad(gray, (1, 1), 'edge')

		# prepare
		Ix = np.zeros_like(gray, dtype=np.float32)
		Iy = np.zeros_like(gray, dtype=np.float32)

		# get differential
		for y in range(H):
			for x in range(W):
				Ix[y, x] = np.mean(tmp[y : y  + 3, x : x + 3] * sobelx)
				Iy[y, x] = np.mean(tmp[y : y + 3, x : x + 3] * sobely)
			
		Ix2 = Ix ** 2
		Iy2 = Iy ** 2
		Ixy = Ix * Iy

		return Ix2, Iy2, Ixy


	# gaussian filtering
	def gaussian_filtering(I, K_size=3, sigma=3):
		# get shape
		H, W = I.shape

		## gaussian
		I_t = np.pad(I, (K_size // 2, K_size // 2), 'edge')

		# gaussian kernel
		K = np.zeros((K_size, K_size), dtype=np.float)
		for x in range(K_size):
			for y in range(K_size):
				_x = x - K_size // 2
				_y = y - K_size // 2
				K[y, x] = np.exp( -(_x ** 2 + _y ** 2) / (2 * (sigma ** 2)))
		K /= (sigma * np.sqrt(2 * np.pi))
		K /= K.sum()

		# filtering
		for y in range(H):
			for x in range(W):
				I[y,x] = np.sum(I_t[y : y + K_size, x : x + K_size] * K)
				
		return I

	# corner detect
	def corner_detect(gray, Ix2, Iy2, Ixy, k=0.04, th=0.1):
		# prepare output image
		out = np.array((gray, gray, gray))
		out = np.transpose(out, (1,2,0))

		# get R
		R = (Ix2 * Iy2 - Ixy ** 2) - k * ((Ix2 + Iy2) ** 2)

		# detect corner
		out[R >= np.max(R) * th] = [0, 0, 255]

		out = out.astype(np.uint8)

		return out

	
	# 1. grayscale
	gray = BGR2GRAY(img)

	# 2. get difference image
	Ix2, Iy2, Ixy = Sobel_filtering(gray)

	# 3. gaussian filtering
	Ix2 = gaussian_filtering(Ix2, K_size=3, sigma=3)
	Iy2 = gaussian_filtering(Iy2, K_size=3, sigma=3)
	Ixy = gaussian_filtering(Ixy, K_size=3, sigma=3)

	# 4. corner detect
	out = corner_detect(gray, Ix2, Iy2, Ixy)

	return out


# Read image
img = cv2.imread("thorino.jpg").astype(np.float32)

# Harris corner detection
out = Harris_corner(img)

cv2.imwrite("out.jpg", out)
cv2.imshow("result", out)
cv2.waitKey(0)

代码来自： https://github.com/gzr2017/ImageProcessing100Wen/blob/master/Question_81_90/README.md

里面的图片资源可在上面的链接中下载。