【论文阅读】Harris角点算法

#coding=utf-8

from PIL import Image

import numpy as np

from scipy.ndimage import filters

import matplotlib.pyplot as plt

import scipy.signal

def mean_filter2d(arr):

    n = 3

    # 3*3 滤波器, 每个系数都是 1/9

    w = np.ones((n, n)) / n ** 2

    # 使用滤波器卷积图像

    # mode = same 表示输出尺寸等于输入尺寸

    # boundary 表示采用对称边界条件处理图像边缘

    s = scipy.signal.convolve2d(arr, w, mode='same', boundary='symm')

    return s

def harris_response1(im,sigma=1.1):

    """计算图像的harris响应函数"""

    #计算导数

    imx = np.zeros(im.shape)

    imx = scipy.ndimage.sobel(im,axis=0,mode='reflect')

    # imx = filters.gaussian_filter(im,(sigma,sigma),(0,1),imx)

    imy = np.zeros(im.shape)

    imy = scipy.ndimage.sobel(im,axis=1,mode='reflect')

    # imy = filters.gaussian_filter(im, (sigma, sigma), (1, 0), imy)

    fig, ax = plt.subplots(1,2)

    ax[0].imshow(imx, cmap='gray')

    ax[1].imshow(imy, cmap='gray')

    plt.show()

    #计算Harris的各个分量

    wxx = filters.gaussian_filter(imx*imx,sigma)

    wxy = filters.gaussian_filter(imx*imy,sigma)

    wyy = filters.gaussian_filter(imy*imy,sigma)

    #计算像素的角点响应函数

    # return (wxx*wyy - 2*wxy)/(wxx + wyy)

    return wxx*wyy - wxy*wxy - 0.04*((wxx + wyy)**2)

def harris_response2(im,sigma=1.1):

    """计算图像的harris响应函数"""

    #计算导数

    imx = np.zeros(im.shape)

    imx = scipy.ndimage.sobel(im,axis=0,mode='reflect')

    # imx = filters.gaussian_filter(im,(sigma,sigma),(0,1),imx)

    imy = np.zeros(im.shape)

    imy = scipy.ndimage.sobel(im,axis=1,mode='reflect')

    # imy = filters.gaussian_filter(im, (sigma, sigma), (1, 0), imy)

    fig, ax = plt.subplots(1,2)

    ax[0].imshow(imx, cmap='gray')

    ax[1].imshow(imy, cmap='gray')

    plt.show()

    #计算Harris的各个分量

    wxx = mean_filter2d(imx*imx)

    wxy = mean_filter2d(imx*imy)

    wyy = mean_filter2d(imy*imy)

    #计算像素的角点响应函数

    return wxx*wyy - wxy*wxy - 0.04*((wxx + wyy)**2)

def get_harris_points(harrism,min_dist = 10,thresold = 0.1):

    """从一幅Harrisim响应中返回角点，min_dist为分割角点和图像边界的最少像素数目"""

    corner_thsold = harrism.max()*thresold

    harrism_t = (harrism > corner_thsold) * 1

    #得到候选点的坐标

    coords = np.array(harrism_t.nonzero()).T#返回非零值的坐标的矩阵

    #他们的Harris响应值

    candidate_values = [harrism[c[0],c[1]] for c in coords]

    #对候选点进行harris响应值进行排序

    index = np.argsort(candidate_values)[::-1]#将x中的元素从小到大排列，提取其对应的index(索引)，然后输出到y

    #将可行点的位置保存在数组里

    allowed_locations = np.zeros(harrism.shape)

    allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1

    #按照min_distance原则，选择最佳harris点

    filters_coords = []

    for i in index:

        if allowed_locations[coords[i,0],coords[i,1]] == 1:

            filters_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 filters_coords

im = np.array(Image.open(r'D:\cvImageSamples\lena.png'),dtype=np.float32)

hr1 = harris_response1(im[:,:,0])

hr2 = harris_response2(im[:,:,0])

fig, ax = plt.subplots(1,3)

ax[0].imshow(im[:,:,0],cmap='gray')

ax[1].imshow(hr1, cmap='gray')

ax[2].imshow(hr2, cmap='gray')

plt.show()

临时起意写的文章，没有写成md格式。只好截图了。文字版本访问https://zhuanlan.zhihu.com/p/148127081