【懒懒的计算机视觉笔记之Harris角点检测一】

Harris角点检测算法，又称（Harris&Stephens角点检测器）是较为简单的一种角点检测算法，主要思想是：如果某点存在多于一个方向的边，则认为该点是兴趣点，称为角点。

我们将图像上某点x上对称的半定矩阵MI定义为：

其中为包含导师Ix和Iy的图像梯度，根据这个公式MI的秩为1，特征值为对于图像上的每个点，我们都可以计算出该矩阵。
选择权重矩阵W（通常选择高斯滤波器G），我们可以得到卷积：

卷积的目的是得到矩阵MI在周围像素上的局部平均，计算出来的称为Harris矩阵。权重矩阵的宽度决定像素x周围的感兴趣区域的范围，之所以取平均的原因是因为Harris矩阵的特征值是受局部像素点的特性变化而变化的。如果图像的梯度在该区域内变化，则Harris矩阵的第二个特征值不再为0，如果图像梯度没有变化，则Harris矩阵特征值也不会与变化。

取决于该区域的，Harris矩阵的特征值有三种情况：
①如果和都是很大的正数，则该点位角点；
②如果很大，而为0，则该区域内存在一条边，且该区域内特征值变化不会很大
③如果和都为0，则该区域为空。
在实际过程中，我们使用指示函数作为指示器：

下面我们使用Python编写Harris角点检测程序，对于这个函数我们要使用scipy.ndimage.filters模块中的高斯滤波器来计算导数，也可以使用Sobel()算子来计算。

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()