Harris角点检测算法,又称(Harris&Stephens角点检测器)是较为简单的一种角点检测算法,主要思想是:如果某点存在多于一个方向的边,则认为该点是兴趣点,称为角点。
我们将图像上某点x上对称的半定矩阵MI定义为:
其中为包含导师Ix和Iy的图像梯度,根据这个公式MI的秩为1,特征值为对于图像上的每个点,我们都可以计算出该矩阵。
选择权重矩阵W(通常选择高斯滤波器G),我们可以得到卷积:
卷积的目的是得到矩阵MI在周围像素上的局部平均,计算出来的称为Harris矩阵。权重矩阵的宽度决定像素x周围的感兴趣区域的范围,之所以取平均的原因是因为Harris矩阵的特征值是受局部像素点的特性变化而变化的。如果图像的梯度在该区域内变化,则Harris矩阵的第二个特征值不再为0,如果图像梯度没有变化,则Harris矩阵特征值也不会与变化。
取决于该区域的,Harris矩阵的特征值有三种情况:
①如果和都是很大的正数,则该点位角点;
②如果很大,而为0,则该区域内存在一条边,且该区域内特征值变化不会很大
③如果和都为0,则该区域为空。
在实际过程中,我们使用指示函数作为指示器:
def compute_harris_response(im, sigma=3):
'''
对每个像素值计算Harris角点检测器响应函数
:param im:
:param sigma:
:return:
'''
# 计算导数
im_x = np.zeros(im.shape)
im_y = np.zeros(im.shape)
filters.gaussian_filter(im, (sigma, sigma), (0, 1), im_x)
filters.gaussian_filter(im, (sigma, sigma), (1, 0), im_y)
# 计算Harris矩阵分量
Ixx = filters.gaussian_filter(im_x * im_x, sigma)
Ixy = filters.gaussian_filter(im_x * im_y, sigma)
Iyy = filters.gaussian_filter(im_y * im_y, sigma)
# 计算特征值和迹
Idet = Ixx * Iyy - Ixy ** 2
Itrace = Ixx + Iyy
return Idet / Itrace
上面的程序会返回像素值为Harris响应函数值得一副图像,但是实际上我们需要从这幅图像上选择需要的信息,首先需要选择像素值高于阈值的的所有像素点,然后再加上额外的限制,即角点之间的间隔必须大于设定的最小值,为此我们需要将像素点排序筛选,并存入数组中,实现代码如下:
def get_harris_points(harrisim, min_dist=10, threshold=0.1):
'''
从Harrris响应图像中筛选角点
:param harrisim:
:param min_dist:
:param threshold:
:return:
'''
corner_threshold = harrisim.max() * threshold
harrisim_t = (harrisim > corner_threshold) * 1
# 获得候选点的坐标和对应的值
coords = np.array(harrisim_t.nonzero()).T
candidate_values = [harrisim[c[0], c[1]] for c in coords]
# 对候选点排序
index = np.argsort(candidate_values)
# 将可行点存储在数组中
allowed_locations = np.zeros(harrisim.shape)
allowed_locations[min_dist:-min_dist, min_dist:-min_dist] = 1
# 按照最小间距原则选择最佳Harris角点
filtered_coords = []
for i in index:
if allowed_locations[coords[i,0], coords[i,1]] == 1:
filtered_coords.append(coords[i])
allowed_locations[(coords[i, 0]-min_dist):(coords[i, 0]+min_dist), (coords[i, 1]-min_dist):(coords[i, 1]+min_dist)] = 0
return filtered_coords
为了在图像上显示角点,我们需要利用pylab模块来显示图像和角点,代码如下:
def plot_corner_points(image, filtered_coords ):
'''
绘制图像中的角点
:param image:
:param filtered_coords:
:return:
'''
imshow(image)
plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*')
axis('off')
最后附上全部的代码和运行结果图:
# coding=UTF-8
from PIL import Image
import numpy as np
from pylab import *
from scipy.ndimage import filters
def compute_harris_response(im, sigma=3):
'''
对每个像素值计算Harris角点检测器响应函数
:param im:
:param sigma:
:return:
'''
# 计算导数
im_x = np.zeros(im.shape)
im_y = np.zeros(im.shape)
filters.gaussian_filter(im, (sigma, sigma), (0, 1), im_x)
filters.gaussian_filter(im, (sigma, sigma), (1, 0), im_y)
# 计算Harris矩阵分量
Ixx = filters.gaussian_filter(im_x * im_x, sigma)
Ixy = filters.gaussian_filter(im_x * im_y, sigma)
Iyy = filters.gaussian_filter(im_y * im_y, sigma)
# 计算特征值和迹
Idet = Ixx * Iyy - Ixy ** 2
Itrace = Ixx + Iyy
return Idet / Itrace
def get_harris_points(harrisim, min_dist=10, threshold=0.1):
'''
从Harrris响应图像中筛选角点
:param harrisim:
:param min_dist:
:param threshold:
:return:
'''
corner_threshold = harrisim.max() * threshold
harrisim_t = (harrisim > corner_threshold) * 1
# 获得候选点的坐标和对应的值
coords = np.array(harrisim_t.nonzero()).T
candidate_values = [harrisim[c[0], c[1]] for c in coords]
# 对候选点排序
index = np.argsort(candidate_values)
# 将可行点存储在数组中
allowed_locations = np.zeros(harrisim.shape)
allowed_locations[min_dist:-min_dist, min_dist:-min_dist] = 1
# 按照最小间距原则选择最佳Harris角点
filtered_coords = []
for i in index:
if allowed_locations[coords[i,0], coords[i,1]] == 1:
filtered_coords.append(coords[i])
allowed_locations[(coords[i, 0]-min_dist):(coords[i, 0]+min_dist), (coords[i, 1]-min_dist):(coords[i, 1]+min_dist)] = 0
return filtered_coords
def plot_corner_points(image, filtered_coords ):
'''
绘制图像中的角点
:param image:
:param filtered_coords:
:return:
'''
imshow(image)
plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*')
axis('off')
if __name__ =='__main__':
im = array(Image.open('1.jpg').convert('L'))
harrsim = compute_harris_response(im)
filter_coords1 = get_harris_points(harrsim, 10, 0.1)
filter_coords2 = get_harris_points(harrsim, 10, 0.4)
filter_coords3 = get_harris_points(harrsim, 10, 0.7)
figure()
gray()
subplot(131)
plot_corner_points(im, filter_coords1)
title('threshold = 0.1')
subplot(132)
plot_corner_points(im, filter_coords2)
title('threshold = 0.4')
subplot(133)
plot_corner_points(im, filter_coords3)
title('threshold = 0.7')
show()