【论文阅读】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