前言
在传统目标识别中，特征提取是最终目标识别效果好坏的一个重要决定因素，因此，在这项工作里，有很多研究者把主要精力都放在特征提取方向。在传统目标识别中，主要使用的特征主要有如下几类：
- 边缘特征
- 纹理特征
- 区域特征
- 角点特征
本文要讲述的Harris角点检测就是焦点特征的一种。
目前角点检测算法主要可归纳为3类：
- 基于灰度图像的角点检测
- 基于二值图像的角点检测
- 基于轮廓的角点检测
因为角点在现实生活场景中非常常见，因此，角点检测算法也是一种非常受欢迎的检测算法，尤其本文要讲的Harris角点检测，可以说传统检测算法中的经典之作。
Harris角点检测
什么是角点？
要想弄明白角点检测，首先要明确一个问题，什么是角点？
这个在现实中非常常见，例如图中标记出的飞机的角点，除此之外例如桌角、房角等。这样很容易理解，但是该怎么用书面的语言阐述角点？
角点就是轮廓之间的交点。
如果从数字图像处理的角度来描述就是：像素点附近区域像素无论是在梯度方向、还是在梯度幅值上都发生较大的变化。
这句话是焦点检测的关键，也是精髓，角点检测算法的思想就是由此而延伸出来的。
角点检测的算法思想是：选取一个固定的窗口在图像上以任意方向的滑动，如果灰度都有较大的变化，那么久认为这个窗口内部存在角点。
要想实现角点检测，需要用数学语言对其进行描述，下面就着重用数学语言描述一下角点检测算法的流程和原理。
用w(x,y)表示窗口函数，[u,v]为窗口平移量，像素在窗口内的变化量为，
其I(x,y)为平移前的像素灰度值，I(x+u,y+v)为平移后的像素灰度值，
通过对灰度变化部分进行泰勒展开得到，
因此得到，
现在在回过头来看一下角点与其他类型区域的不同之处：
- 平坦区域：梯度方向各异，但是梯度幅值变化不大
- 线性边缘：梯度幅值改变较大，梯度方向改变不大
- 角点：梯度方向和梯度幅值变化都较大
明白上述3点之后看一下怎么利用其矩M进行角点检测。
根据主成分分析(PCA)的原理可知，如果对矩M对角化，那么，特征值就是主分量上的方差，M是二维的方阵，有两个主分量，如果在窗口区域内是角点，那么梯度变化会较大，像素点的梯度分布比较离散，这体现在特征值上就是特征值比较大。
换句话说，
- 如果矩阵对应的两个特征值都较大，那么窗口内含有角点
- 如果特征值一个大一个小，那么窗口内含有线性边缘
- 如果两个特征值都很小，那么窗口内为平坦区域
读到这里就应该明白了，角点的检测转化为数学模型，就是求解窗口内矩阵M的特征值并且判断特征值的大小。
如果要评价角点的强度，可以用下方公式，
其中，
编程实践
1requirement: OpenCV/scipy/matplotlib
因为Harris角点检测算法非常经典，因此，一些成熟的图像处理或视觉库都会直接提供Harris角点检测的算法，以OpenCV为例，
1import cv2 2import numpy as np 3 4filename = '2007_000480.jpg' 5 6img = cv2.imread(filename) 7gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY) 8gray = np.float32(gray) 9dst = cv2.cornerHarris(gray,blockSize,ksize,k) 10 11""" 12其中， 13gray为灰度图像， 14blockSize为邻域窗口的大小， 15ksize是用于Soble算子的参数， 16k是一个常量，取值为0.04~0.06 17"""
因为本文要讲一步一步实现Harris角点检测算法，因此，对OpenCV提供的函数不多阐述，下面开始一步一步实现Harris角点检测算法。
检点检测算法的流程如下：
- 利用公式(1)求出输入图像每个位置的角点强度相应
- 给定阈值，当一个位置的强度大于阈值则认为是角点
- 画出角点
首先是第一步，根据上述提到的公式求矩阵的特征值和矩阵的迹，然后计算图像的角点强度，这里选取常数k=0.04，
1def calculate_corner_strength(img, scale=3, k=0.06): 2 # 计算图像在x、y方向的梯度 3 # 用滤波器采用差分求梯度的方式 4 gradient_imx, gradient_imy = zeros(img.shape), zeros(img.shape) 5 filters.gaussian_filter(img, (scale, scale), (0, 1), gradient_imx) 6 filters.gaussian_filter(img, (scale, scale), (1, 0), gradient_imy) 7 8 # 计算矩阵M的每个分量 9 I_xx = filters.gaussian_filter(gradient_imx*gradient_imx, scale) 10 I_xy = filters.gaussian_filter(gradient_imx*gradient_imy, scale) 11 I_yy = filters.gaussian_filter(gradient_imy*gradient_imy, scale) 12 13 # 计算矩阵的迹、特征值和响应强度 14 det_M = I_xx * I_yy - I_xy ** 2 15 trace_M = I_xx + I_yy 16 return det_M + k * trace_M ** 2
接下来完成第2步，根据给定阈值，获取角点，
1def corner_detect(img, min=15, threshold=0.04): 2 # 首先对图像进行阈值处理 3 _threshold = img.max() * threshold 4 threshold_img = (img > _threshold) * 1 5 corners = array(threshold_img.nonzero()).T 6 candidate_values = [img[c[0], c[1]] for c in corners] 7 index = argsort(candidate_values) 8 9 # 选取领域空间，如果邻域空间距离小于min的则只选取一个角点 10 # 防止角点过于密集 11 neighbor = zeros(img.shape) 12 neighbor[min:-min, min:-min] = 1 13 detect_corner= [] 14 for i in index: 15 if neighbor[corners[i, 0], corners[i, 1]] == 1: 16 detect_corner.append(corners[i]) 17 neighbor[(corners[i, 0] - min):(corners[i, 0] + min), 18 (corners[i, 1] - min):(corners[i, 1] + min)] = 0 19 return detect_corner
然后是画出角点，
1def corner_plot(image, corner_img): 2 figure() 3 gray() 4 imshow(image) 5 plot([p[1] for p in corner_img], [p[0] for p in corner_img], 'ro') 6 axis('off') 7 show()
检测结果，
完整代码如下，
1from scipy.ndimage import filters 2import cv2 3from matplotlib.pylab import * 4 5 6class Harris(object): 7 def __init__(self, img_path): 8 self.img = cv2.imread(img_path, 0) 9 10 def calculate_corner_strength(self): 11 # 计算图像在x、y方向的梯度 12 # 用滤波器采用差分求梯度的方式 13 scale = self.scale 14 k = self.k 15 img = self.img 16 gradient_imx, gradient_imy = zeros(img.shape), zeros(img.shape) 17 filters.gaussian_filter(img, (scale, scale), (0, 1), gradient_imx) 18 filters.gaussian_filter(img, (scale, scale), (1, 0), gradient_imy) 19 20 # 计算矩阵M的每个分量 21 I_xx = filters.gaussian_filter(gradient_imx*gradient_imx, scale) 22 I_xy = filters.gaussian_filter(gradient_imx*gradient_imy, scale) 23 I_yy = filters.gaussian_filter(gradient_imy*gradient_imy, scale) 24 25 # 计算矩阵的迹、特征值和响应强度 26 det_M = I_xx * I_yy - I_xy ** 2 27 trace_M = I_xx + I_yy 28 return det_M + k * trace_M ** 2 29 30 def corner_detect(self, img): 31 # 首先对图像进行阈值处理 32 _threshold = img.max() * self.threshold 33 threshold_img = (img > _threshold) * 1 34 corners= array(threshold_img.nonzero()).T 35 candidate_values = [img[c[0], c[1]] for c in corners] 36 index = argsort(candidate_values) 37 38 # 选取领域空间，如果邻域空间距离小于min的则只选取一个角点 39 # 防止角点过于密集 40 neighbor = zeros(img.shape) 41 neighbor[self.min:-self.min, self.min:-self.min] = 1 42 detect_corner= [] 43 for i in index: 44 if neighbor[corners[i, 0], corners[i, 1]] == 1: 45 detect_corner.append(corners[i]) 46 neighbor[(corners[i, 0] - self.min):(corners[i, 0] + self.min), 47 (corners[i, 1] - self.min):(corners[i, 1] + self.min)] = 0 48 return detect_corner 49 50 def corner_plot(self, img, corner_img): 51 figure() 52 gray() 53 imshow(img) 54 plot([p[1] for p in corner_img], [p[0] for p in corner_img], 'ro') 55 axis('off') 56 show() 57 58 def __call__(self, k=0.04, scale=3, min=15, threshold=0.03): 59 self.k = k 60 self.scale = scale 61 self.min = min 62 self.threshold = threshold 63 strength_img = self.calculate_corner_strength() 64 corner_img = self.corner_detect(strength_img) 65 self.corner_plot(self.img, corner_img) 66 67 68if __name__ == '__main__': 69 harris = Harris("2007_002619.jpg") 70 harris()
