写在前面：

七月，骄阳似火，好想dying in the sun。黄宁然，为何还在这里。

参考文献镇楼：

[1]袁志聪，基于harris特征的点云配准方法研究
[2]高亭，基于改进Harris角点检测的印刷体文档图像检索技术
[3]景庆阳，基于harris角点检测的KCF算法研究.
[4]zoukankan，OpenCV计算机视觉学习（13）——图像特征点检测（Harris角点检测，sift算法），
[5] Harris角点检测原理【OpenCV教程】
[6] 勿在浮砂筑高台，【特征匹配】Harris及Shi-Tomasi原理及源码解析
注：csdn发文助手提示“外链过多”，故文献链接网址请见评论区。

问题来源：

笔者的执念。

1、原理简介

Harris算法是chris harris等提出的一种基于信号的点特征提取算法，是对Moravec算子的改良和优化[1]。Moravec算子主要是研究图像的固定窗口在四个不同方向上的移动，比较窗口中的像素灰度变化情况的大小，而harris角点检测则是比较窗口在任意方向滑动前后的灰度变化大小的情况，并且用数学表达式确定特征点的位置；同时引入了平滑因子，进一步增强了抗干扰能力[2]。若窗口移动前后，灰度发生了较大变化，则判断窗口中存在角点，否则窗口内不存在角点[1]。根据文献[3]，具体公式推导如下。
当窗口的偏移量为[j,q]时，滑动前后窗口的灰度变化为：
在这里插入图片描述

式中，Window(x,y)是用来滤波的高斯函数，I(x+j,y+q)为移动后的图像灰度，I(x,y)为原始图像灰度。
根据泰勒公式，二次展开为：
在这里插入图片描述
则式（1）变换为：

其中，Ix、Iy分别为x、y方向上的偏导数
可将式（4）简写成：
在这里插入图片描述
其中

矩阵M也可以进一步写成：

其中，A、B、C分别是 $I_x^2$ 、 $I_y^2$ 、 $I_xy$ 在窗口w内的求和，即：

上述就得到自相关函数的表达式。据文献4，这是一椭圆的矩阵表达形式（非标准椭圆），所以系数矩阵M的特征值与椭圆的半轴长短有关。而椭圆的半轴长短也影响到曲率，也即M的特征值会影响到自相关函数的变化快慢。如何影响？可看文献4的讲解（笔者数学不好，这里也是一知半解）。
根据文献4，对于一标准椭圆，其表达式为：
在这里插入图片描述
矩阵形式为：

则系数矩阵：

则特征值与半轴的关系为：

结论是：特征值大，半轴短。

便得到边界、平面、角点几种情况下，特征值的大小情况：
（1）边界。一个特征值大，另一个特征值小，λ1 >> λ2 或者 λ2 >> λ1。自相关函数值在某一个方向上大，在其他方向上小。
（2）平面。两个特征值都小，且近似相等；自相关函数数值在各个方向上都小。
（3）角点。两个特征值都很大且近似相等，自相关函数在所有方向都增大。
在这里插入图片描述
实际中，Harris算法并不需要计算具体的特征值，而是根据角点响应值R来判断角点[2]，R的定义为：

其中，detM为矩阵M的行列式，traceM为矩阵M的迹，k为Harris算子参数，据文献3，其一般取值为[0.04，0.06]。据文献5，角点响应值判断角点的依据如下：
在这里插入图片描述
即：
（1）当R为大数值的正数时是角点
（2）当R为大数值的负数时是边界
（3）当R为小数是认为是平坦区域
综上，Harris 算法分以下几个步骤[2]：
（1）计算图像在水平方向和垂直方向的导数 $I x$ 和 $I y$ 以及 $I_x^2$ 、 $I_y^2$ 、 $I_xy$ 。
（2）根据式（8）对上述元素进行平滑滤波，以得到系数 A、B、C，得到M。
（3）将求得的系数带入式（13）~式（15），来计算角点响应值 R。
（4）对于各响应值R，若大于设定好的阈值，则该R对应的坐标即为角点。

2、python源码实现

参阅文献6，基于python自行编写代码实现harris角点检测。几点需要注意的地方：
（1）通过sobel算子，快速实现图像x、y方向上导数的求取；
（2）opencv在实现harris算子时，窗口函数并非取的高斯函数，而是普通的boxFilter函数（即权值均为1，且求和后，不归一化）。
引用库：

#coding=utf8
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import math
import cv2
import os,sys
import scipy.ndimage
import time
import scipy

harris角点检测代码：

```bash
def harris_corner_detect(img_src,block_size=2,aperture_size=3,k=0.04,borderType=cv2.BORDER_DEFAULT):
	if img_src.dtype!=np.uint8:
		raise ("input image shoud be uint8 type")
	R_arr = np.zeros(img_src.shape,dtype=np.float32)#用来存储角点响应值
	img = img_src.astype(np.float32)
	scale = 1.0/( (aperture_size-1)*2*block_size*255 )#参考opencv实现源码，在sobel算子时乘以这个系数
	#借助sobel算子，计算x、y方向上的偏导数
	Ix = cv2.Sobel(img,-1,dx=1,dy=0,ksize=aperture_size,scale=scale,borderType=borderType)
	Iy = cv2.Sobel(img,-1,dx=0,dy=1,ksize=aperture_size,scale=scale,borderType=borderType)
	Ixx = Ix**2
	Iyy = Iy**2
	Ixy = Ix*Iy
	#借助boxFilter函数，以block_size为窗口大小，对窗口内的数值求和，且不归一化
	f_xx = cv2.boxFilter(Ixx,ddepth=-1,ksize=(block_size,block_size) ,anchor =(-1,-1),normalize=False, borderType=borderType)
	f_yy = cv2.boxFilter(Iyy,ddepth=-1,ksize=(block_size,block_size),anchor =(-1,-1),normalize=False,borderType=borderType)
	f_xy = cv2.boxFilter(Ixy, ddepth=-1,ksize=(block_size,block_size),anchor =(-1,-1),normalize=False,borderType=borderType)
	# 也可以尝试手动求和
	radius = int((block_size - 1) / 2)  # 考虑blocksize为偶数的情况，奇数时，前、后的数量一样；为偶数时，后比前多一个
	N_pre = radius
	N_post = block_size - N_pre - 1
	row_s, col_s = N_pre, N_pre
	row_e, col_e = img.shape[0] - N_post, img.shape[1] - N_post
	#开始计算每一个坐标下的响应值
	for r in range(row_s,row_e):
		for c in range(col_s,col_e):
			#手动对窗口内的数值求和
			#sum_xx = Ixx[r-N_pre:r+N_post+1,c-N_pre:c+N_post+1].sum()
			#sum_yy = Iyy[r-N_pre:r+N_post+1,c-N_pre:c+N_post+1].sum()
			#sum_xy = Ixy[r-N_pre:r+N_post+1,c-N_pre:c+N_post+1].sum()
			#或者直接使用boxFilter的结果
			sum_xx = f_xx[r,c]
			sum_yy = f_yy[r, c]
			sum_xy = f_xy[r, c]
			#根据行列式和迹求响应值
			det = (sum_xx * sum_yy) - (sum_xy ** 2)
			trace = sum_xx + sum_yy
			res = det - k * (trace ** 2)
			# 或者用如下方式求行列式和迹
			#M = np.array([[sum_xx,sum_xy],[sum_xy,sum_yy]])
			#res = np.linalg.det(M) - (k *  (np.trace(M))**2 )
			R_arr[r,c] = res
	return R_arr

主程序调用：
这里，响应值的门限选择为最大响应值得0.01倍。

if __name__ == '__main__':
	img_src = cv2.imread('susan_input1.png',cv2.IMREAD_GRAYSCALE)
	#source code
	res_arr = harris_corner_detect(img_src,block_size=2,aperture_size=3,k=0.04)
	max_v = np.max(res_arr)
	res_arr[res_arr<0.01*max_v]=0
	img_show = img_src.copy()
	if(len(img_show.shape)==2):
		img_show = cv2.cvtColor(img_show,cv2.COLOR_GRAY2BGR)
	img_show[res_arr!=0] = (255,0,0)
	print(len(np.where(res_arr!=0)[0]))
	print(np.max(res_arr))
	plt.figure()
	plt.title("corners-raw")
	plt.imshow(img_show, cmap=cm.gray)
	plt.show()

检测结果如下：
在这里插入图片描述

3、基于opencv实现

Opencv自带harris算子，可以直接求取响应矩阵。这里使用的opencv版本为3.4.2.16
调用方式为：
dst = cv2.cornerHarris(img_src,blockSize= 2,ksize= 3,k= 0.04)
主程序调用：

if __name__ == '__main__':
	img_src = cv2.imread('susan_input1.png',cv2.IMREAD_GRAYSCALE)
	#source code
	res_arr = harris_corner_detect(img_src,block_size=2,aperture_size=3,k=0.04)
	max_v = np.max(res_arr)
	res_arr[res_arr<0.01*max_v]=0
	img_show = img_src.copy()
	if(len(img_show.shape)==2):
		img_show = cv2.cvtColor(img_show,cv2.COLOR_GRAY2BGR)
	img_show[res_arr!=0] = (255,0,0)
	print(len(np.where(res_arr!=0)[0]))
	print(np.max(res_arr))
	plt.figure()
	plt.title("corners-raw")
	plt.imshow(img_show, cmap=cm.gray)
	#opencv库函数调用
	dst = cv2.cornerHarris(img_src,blockSize= 2,ksize= 3,k= 0.04) #blockSize为窗口大小，ksize为sobel算子的核大小，k为harris算子参数
	img_show2 = img_src.copy()
	if (len(img_show2.shape) == 2):
		img_show2 = cv2.cvtColor(img_show2, cv2.COLOR_GRAY2BGR)
	dst2 = dst.copy()
	dst2[dst<= 0.01* dst.max()]=0
	img_show2[dst2!=0] = (255, 0, 0)
	print(len(np.where(dst2 != 0)[0]))
	print(np.max(dst2))
	plt.figure()
	plt.title("opencv")
	plt.imshow(img_show2, cmap=cm.gray)
	plt.show()

opencv库函数检测结果如下：
在这里插入图片描述
将自行编写的harris角点检测代码运行结果与opencv运行结果相对比，二者一致（仅仅比较了数量、最大值）。

4、非极大值抑制

可根据需要，执行非极大值抑制。

非极大值抑制代码：

def corner_nms(corner,kernal=3):
	out = corner.copy()
	row_s = int(kernal/2)
	row_e = out.shape[0] - int(kernal/2)
	col_s,col_e = int(kernal/2),out.shape[1] - int(kernal/2)
	for r in range(row_s,row_e):
		for c in range(col_s,col_e):
			if corner[r,c]==0: #不是可能的角点
				continue
			zone = corner[r-int(kernal/2):r+int(kernal/2)+1,c-int(kernal/2):c+int(kernal/2)+1]
			index = corner[r,c]<zone
			(x,y) = np.where(index==True)
			if len(x)>0 : #说明corner[r,c]不是最大，直接归零将其抑制
				out[r,c] = 0
	return out

在上文基础上，调用nms：

	#nms
	score_nms = corner_nms(res_arr)
	img_show3 = cv2.cvtColor(img_src,cv2.COLOR_GRAY2BGR)
	img_show3[score_nms != 0] = (255, 0, 0)
	plt.figure()
	plt.title("corners-nms")
	plt.imshow(img_show3, cmap=cm.gray)

	plt.show()