写在前面：
黄宁然， 七月，骄阳似火。
参考文献镇楼：
[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_x^2 I x 2 、I y 2 I_y^2 I y 2 、I x y I_xy I x y 在窗口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 Ix I x 和I y Iy I y 以及I x 2 I_x^2 I x 2 、I y 2 I_y^2 I y 2 、I x y I_xy I x y 。
（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，且求和后，不归一化）。
引用库：
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角点检测代码：
主程序调用：
这里，响应值的门限选择为最大响应值得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()
经过nms抑制后的结果：
5、python源码下载
Python程序源码下载地址
https://download.csdn.net/download/xiaohuolong1827/85847444
Original: https://blog.csdn.net/xiaohuolong1827/article/details/125559091
Author: 小火龙的马甲
Title: Harris角点检测python实现及基于opencv实现
原创文章受到原创版权保护。转载请注明出处：https://www.johngo689.com/672225/
转载文章受原作者版权保护。转载请注明原作者出处！