python实现harris角点检测以及边缘提取（没有使用现成的cv函数）

原理：特征提取算法（4）

#%%导入库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import cv2 as cv
import math
# import scipy.io as sc

#%%数据初始处理
im_1_path=r'C:\Users\cz\Desktop\ceshi.jpg'
im_2_path=r'C:\Users\cz\Desktop\WIN_20200813_14_31_15_Pro.jpg'
im_3_path=r'C:\Users\cz\Desktop\WIN_20200815_15_22_34_Pro.jpg'
im_4_path=r'C:\Users\cz\Desktop\20200815152820.jpg'
im_5_path=r'C:\Users\cz\Desktop\nnn.jpg'
src = cv.imread(im_2_path)                         
plt.imshow(src[:,:,[2,1,0]])                          

#%%准备图片
im1=cv.cvtColor(src,cv.COLOR_BGR2GRAY).astype(float)
plt.imshow(im1,cmap='gray')

#%% 定义函数
def gauss_5x5(a):#高斯滤波算子计算函数
    computer=np.ones(dtype=float,shape=(5,5))
    for i in range(-2,3):
        for j in  range(-2,3):
            computer[i+2,j+2]=pow(math.e,(-i*i-j*j)/(2*a*a))
    computer=computer/np.min(computer)
    computer=computer.astype(int)+1
    computer=np.where(computer==2,1,computer)
    computer[2,2]=56
    computer=computer.astype(float)
    return computer

def convolution(computer,image_1):#泛用卷积函数
    image_2=image_1.copy()
    size_i=int(computer.shape[0]/2)
    a=[0,0]
    for b in range(0,2):
        if computer.shape[b]%2==0:
            a[b]=0
        else:
            a[b]=1
    size_j=int(computer.shape[1]/2)
    sum_1=np.sum(computer)
    if sum_1==0:
        sum_1=1
    for i_1 in range(size_i,image_1.shape[0]-size_i):
        for j_1 in range(size_j,image_1.shape[1]-size_j):
            image_2[i_1,j_1]=np.sum(image_1[i_1-size_i:i_1+size_i+a[1],j_1-size_j:j_1+size_j+a[0]]*computer)/sum_1
    return image_2

def edg(edg_cp_x,edg_cp_y,im):
    im_1=im.copy()
    im_a=convolution(edg_cp_x,im_1)
    im_b=convolution(edg_cp_y,im_1)
    im_2=np.sqrt(im_a*im_a+im_b*im_b)
    return im_2

def point_extract(image,a,b=255):#返回角点检测地址
    image=image*255/np.max(image)
    # image=np.absolute(image)
    arry=np.where(image>=a)
    arry=np.array(arry)
    return arry

#%%进行滤波
detect=gauss_5x5(1)
im2=convolution(detect,im1)
plt.imshow(im2,cmap='gray')

#%%roberts边缘算子计算2x2
edg1_cp_x=np.array([1,0,0,-1]).reshape(2,2)
edg1_cp_y=np.array([0,-1,1,0]).reshape(2,2)
im5=edg(edg1_cp_x,edg1_cp_y,im2)

#%%sobel边缘算子计算3x3
edg2_cp_x=np.array([-1,0,1,-1,0,1,-1,0,1]).reshape(3,3)
edg2_cp_y=-edg2_cp_x.T
im6=edg(edg2_cp_x,edg2_cp_y,im2)

#%% Harris角点检测
H_x=np.array([-1,0,1,-1,0,1,-1,0,1],dtype=float).reshape(3,3)
H_y=-H_x.T

H_i_x=convolution(H_x,im2)
H_i_y=convolution(H_y,im2)

H_A=convolution(detect,H_i_x*H_i_x)
H_B=convolution(detect,H_i_y*H_i_y)
H_C=convolution(detect,H_i_x*H_i_y)

M=np.array([[H_A,H_C],
            [H_C,H_B]])
M_T=M.transpose((2,3,0,1))
M_det=np.linalg.det(M_T)
M_trace=np.trace(M)
M_trace=M_trace.astype(float)
im3=M_det-0.04*M_trace*M_trace

#%%归一化角点提取
l1=point_extract(im3,110)

im7=np.float32(im1)
im4=cv.cornerHarris(im7,2,3,0.04)

l2=point_extract(im4,160)

#%% 显示图像
fig, axs = plt.subplots(3,2, figsize=(40,40))
axs[0,0].imshow(im1,cmap='gray')
axs[0,0].scatter(l1[1,:],l1[0,:],marker='o',color='red',s=200)
axs[0,1].imshow(im2,cmap='gray')
axs[0,1].scatter(l2[1,:],l2[0,:],marker='o',color='red',s=200)
axs[1,0].imshow(im3,cmap='gray')
axs[1,1].imshow(im4,cmap='gray')
axs[2,0].imshow(im5,cmap='gray')
axs[2,1].imshow(im6,cmap='gray')
plt.show()

不完善的地方：代码实现后和cv库的harris检测函数检测结果基本相同，但是同一个点受到分辨率影响可能出现多个检测点，不知道如何解决，怎么才能计算正确的对应点。（好像可以用局部非最大值抑制）

编译过程中出现的问题：最开始没有注意矩阵平方计算的时候使用的整数数据类型，推测是整数平方后数值过大溢出导致计算错误而且没有报错，将图片所有的数据强制转为浮点数之后计算正确。

检测结果如图：

欢迎大佬提出修改建议！