Harris角点检测(Python-OpenCV)
author@jason
1 角点
- 角点：是一类含有足够信息且能从当前帧和下一帧中都能提取出来的点。
- 最普遍使用的角点的定义是由
Harris提出的。 - 典型的角点检测算法：
Harris角点检测、CSS角点检测 - 好的角点检测算法的特点：1、检测出图像中“
真实的”角点；2、准确的定位性能；3、很高的重复检测率（稳定性好）；4、具有对噪声的鲁棒性；5、具有较高的计算效率
2 Harris角点检测原理
图像：I
图像的灰度值(指单个像素值)：I(x,y)
图像的平移向量：[u,v]
平移产生的灰度变化：E(u,v)
高斯窗口系数：w(x,y)
在我们解决问题时，往往希望找到特征点，“特征”顾名思义，指能描述物体本质的东西，还有一种解释就是这个特征微小的变化都会对物体的某一属性产生重大的影响。而角点就是这样的特征。
观察日常生活中的“角落”就会发现，“角落”可以视为所有平面的交汇处，或者说是所有表面的发起处。假设我们要改变一个墙角的位置，那么由它而出发的平面势必都要有很大的变化。所以，这就引出了图像角点的定义
“如果某一点在任意方向的一个微小变动都会引起灰度很大的变化，那么我们就把它称之为角点”.
由上面的定义，检测图像像素的灰度变化情况，即求解E(u,v)
E(u,v)=∑x,yw(x,y)[I(x+u,y+v)−I(x,y)]2(1)
要求解E(u,v)的最大值的点，首先对I(x+u,y+v)进行泰勒展开：
I(x+u,y+v)=I(x,y)+Ixu+Iyv+O(u2,v2)(2)
将（2）式代入（1）式：
E(u,v)=∑x,yw(x,y)[I(x,y)+Ixu+Iyv−I(x,y)]2=∑x,yw(x,y)[Ixu+Iyv]2=∑x,yw(x,y)[u,v][I2xIxIyIxIyI2y][uv]=[u,v]{∑x,yw(x,y)[I2xIxIyIxIyI2y]}[uv]=[u,v]M[uv](3)
其中
M=∑x,yw(x,y)[I2xIxIy IxIyI2y]
Harris算子：
R=det(M)−k∗trace2(M)
也可表示为：
R=|M|−k∗tr2(M)
即
R=λ1λ2−k(λ1+λ2)2
其中trace(M)表示矩阵M的迹，即对角线元素的和。k为经验值，一般在0.04-0.06之间，只要在某一点(x,y)处，λ1,λ2为M的特征值。
R超过某一设定的阈值，则认为该点为角点。阈值的设置是Harris角点检测算法的关键，阈值设置过大会产生伪角点，阈值过小则会漏掉角点。
像上面这样去求解，原则上是没有问题的，可是，众所周知，原始的Harris角点检测算法不具有尺度不变性，也就是说如果图像的尺度发生变化，那么可能原来是角点的点在新的尺度就不是角点了。所以，我们在进行运算的开始，先将图像转换到尺度空间表示，即将原图像进行尺度变换，而尺度变换的方式就是问题的输入信号与尺度核函数做卷积运算：
fout=K∗fin，这里的运算为卷积运算，不是乘法运算，即L(x,y,σ)=G(x,y,σ)∗I(x,y),其中σ表示尺度，然后我们用L代替原图像去进行运算，而尺度成了我们运算的参数了。
- Harris角点检测实现步骤：
1、求出I(x,y)在x、y方向上的梯度Ix、Iy
2、分别求出在x、y方向上的梯度乘积，I2x=Ix∗Ix,I2y=Iy∗Iy,Ixy=Ix∗Iy
3、对I2x，I2y，Ixy进行高斯加权，从而产生A，B，C三个元素，如下：
A=g(I2x)=I2x∗w，B=g(I2y)=I2y∗w，C=g(I2xy)=I2xy∗w
4、求出每个像素的Harris响应值R，令小于阈值（阈值一般为0.01*np.max(R)）的响应值R为0
5、进行3*3邻域非极大值抑制，即如果该点的值比3*3邻域的其它角点小，则删除该角点。
6、记录下角点在原图像的位置，即图像角点所在的位置3 Harris角点检测算法代码
import cv2
import numpy as np
import matplotlib
import math
from matplotlib import pyplot as plt
#根据一阶锐化算子，求x，y的梯度，显示锐化图像
#读取图片
filename = 'girl.jpg'
tu = cv2.imread(filename)
#转换为灰度图
gray = cv2.cvtColor(tu, cv2.COLOR_RGB2GRAY)
#获取图像属性
print '获取图像大小: '
print gray.shape
print '\n'
#打印数组gray
print '灰度图像数组：\n %s \n \n' % (gray)
#输出n*n的数组
#print gray[:2,:2]
#转换为矩阵
m = np.matrix(gray)
#计算x方向的梯度的函数（水平方向锐化算子）
delta_h = m
def grad_x(h):
a = int(h.shape[0])
b = int(h.shape[1])
for i in range(a):
for j in range(b):
if i-1>=0 and i+1<a and j-1>=0 and j+1<b:
c = abs(int(h[i-1,j-1]) - int(h[i+1,j-1]) + 2*(int(h[i-1,j]) - int(h[i+1,j])) + int(h[i-1,j+1]) - int(h[i+1,j+1]))
# print c
if c>255:
# print c
c = 255
delta_h[i,j] = c
else:
delta_h[i,j] = 0
print 'x方向的梯度：\n %s \n' %delta_h
return delta_h
##计算y方向的梯度的函数（水平方向锐化算子）
def grad_y(h):
a = int(h.shape[0])
b = int(h.shape[1])
for i in range(a):
for j in range(b):
if i-1>=0 and i+1<a and j-1>=0 and j+1<b:
c = abs(int(h[i-1,j-1]) - int(h[i-1,j+1]) + 2*(int(h[i,j-1]) - int(h[i,j+1])) + (int(h[i+1,j-1]) - int(h[i+1,j+1]))) #注意像素不能直接计算，需要转化为整型
# print c
if c > 255:
c = 255
delta_h[i,j] = c
else:
delta_h[i,j] = 0
print 'y方向的梯度：\n %s \n' %delta_h
return delta_h
# Laplace 算子
img_laplace = cv2.Laplacian(gray, cv2.CV_64F, ksize=3)
dx = np.array(grad_x(gray))
dy = np.array(grad_y(gray))
#dxy = dx + dy
#print 'dxy1:'
#print dxy
A = dx * dx
B = dy * dy
C = dx * dy
print A
print B
print C
A1 = A
B1 = B
C1 = C
A1 = cv2.GaussianBlur(A1,(3,3),1.5)
B1 = cv2.GaussianBlur(B1,(3,3),1.5)
C1 = cv2.GaussianBlur(C1,(3,3),1.5)
print A1
print B1
print C1
a = int(gray.shape[0])
b = int(gray.shape[1])
R = np.zeros(gray.shape)
for i in range(a):
for j in range(b):
M = [[A1[i,j],C1[i,j]],[C1[i,j],B1[i,j]]]
R[i,j] = np.linalg.det(M) - 0.06 * (np.trace(M)) * (np.trace(M))
print R
cv2.namedWindow('R',cv2.WINDOW_NORMAL)
cv2.imshow('R',R)
cv2.waitKey(0)
cv2.destroyAllWindows()参考：
(1)、Harris角点检测原理与流程(2)、Harris角点检测算法原理及其MATLAB编程实现
http://wenku.baidu.com/view/3c0dfdf7680203d8cf2f246d.html?from=search
4 FAST角点检测算法
链接：article/details/65662648
