CV1-Harris Corner Detector 的 Python3 实现

算法

Harris Corner Detector 的原理不讲，详见谭平老师的计算机视觉P9，觉得讲得不错，也可以从百度网盘下载该算法的文档 提取码: dixx。算法如下：

结果如下

第一行从左到右分别是：原图，水平方向梯度图（Gaussian梯度算子），竖直方向梯度图（Gaussian梯度算子）
第二行从左到右分别是：Harris算法检测到的边缘，Harris检测到的角点，nonmax后检测到的角点（有的地方只有一个像素看不清）

代码

import numpy as np 
import matplotlib.pyplot as plt 
import cv2

def scale01(img):
    """ 把图像缩放到0-1之间，便于plt显示 """
    return (img - np.min(img))/(np.max(img) - np.min(img))

def nonmax_suppression(I):
    # 对I的每一列，如果某个像素左右与之相等，则只取一个像素点
    for x in range(I.shape[0]):
        for y in range(1,I.shape[1]-1):
            if I[x,y-1] >= I[x,y] or I[x,y+1] >= I[x,y]:
                I[x,y] = 0.0
    return I
            
# 读取图片并放缩到0-1之间，imread 的参数 0 表示转换成灰度图
I = cv2.imread("D:/picture/test/qipan2.jpg", 0)
I = np.array(I, dtype=np.float32) / 255 

# 对应算法第 1 步中的 Gx Gy 两个梯度滤波器
GxFilter = np.array([[-1,-3,-1],[0,0,0],[1,3,1]], dtype=np.float32)
GyFilter = np.array([[-1,0,1],[-3,0,3],[-1,0,1]], dtype=np.float32)
# filter2D 参数 -1 表示保持图像 dtype 不变
Ix = cv2.filter2D(I, -1, GxFilter)
Iy = cv2.filter2D(I, -1, GyFilter)

Ix2 = Ix * Ix 
Iy2 = Iy * Iy 
Ixy = Ix * Iy 

# Gaussian 7x7
ksize = (7,7)
Sx2 = cv2.GaussianBlur(Ix2, ksize, 0)
Sy2 = cv2.GaussianBlur(Iy2, ksize, 0)
Sxy = cv2.GaussianBlur(Ixy, ksize, 0)

# 建议 k=0.04~0.06
k = 0.05

# IR 图像对应的 R 值，通过逐像素计算
IR = np.zeros_like(I)
for x in range(I.shape[0]):
    for y in range(I.shape[1]):
        H = np.array([
            [Sx2[x,y], Sxy[x,y]],
            [Sxy[x,y], Sy2[x,y]]
        ])
        IR[x,y] = np.linalg.det(H) - k* H.trace()

# IR 的最小值，平均值，中位数，最大值
print(np.min(IR), np.mean(IR), np.percentile(IR, 50), np.max(IR))

# IR 较小的是边，较大的是角点，不大不小的什么都不是
# 用分位数来确定边（小于10%）和角点（大于90%）, 根据整体 IR 的大小决定的
edges = np.array(IR<np.percentile(IR, 10), dtype=np.float32)
points = np.array(IR>np.percentile(IR, 90), dtype=np.float32)

# 画图
fig, axes = plt.subplots(2,3,figsize=(10,8))
plt.axis('off')
axes[0,0].imshow(I, cmap='gray')
axes[0,0].set_title("Source Image")
axes[0,1].imshow(scale01(Ix), cmap='gray')
axes[0,1].set_title("X derivative")
axes[0,2].imshow(scale01(Iy), cmap='gray')
axes[0,2].set_title("Y derivative")
axes[1,0].imshow(edges, cmap='gray')
axes[1,0].set_title("harris edges")
axes[1,1].imshow(points, cmap='gray')
axes[1,1].set_title("harris corner")
axes[1,2].imshow(nonmax_suppression(points), cmap='gray')
axes[1,2].set_title("harris corner nonmax")
plt.show()