计算机视觉实践之角点检测--Python实现--sobel滤波--Harris矩阵--角点检测器R--非极大值抑制

前言

实现对图像的角点检测，编写Python程序能够对输入图像进行角点检测，并返回角点检测的结果，并且在可视化显示输出结果。
而具体的文件和代码，在如下链接里：计算机视觉实践之角点检测Python实现

使用的是pycharm，提前导入numpy和opencv库

一、方法原理

1. Harris角点检测

算法的核心是利用局部窗口在图像上进行移动，判断灰度是否发生较大的变化。如果窗口内的灰度值（在梯度图上）都有较大的变化，那么这个窗口所在区域就存在角点。
具体的方法是，计算像素点在水平方向和竖直方向的梯度组成的Harris矩阵，其矩阵特征值就反映了该点的像素值变化。
如果该点所对应的两个特征值都比较大，说明像素点的梯度分布比较散，梯度变化程度比较大，符合角点在窗口区域的特点；
如果是平坦区域，那么像素点的梯度所构成的点集比较集中在原点附近，因为窗口区域内的像素点的梯度幅值非常小，此时Harris矩阵的对角化的两个特征值都比较小；
如果是边缘区域，在计算像素点的水平方向和竖直方向上的梯度时，边缘上的像素点的某个方向的梯度幅值变化比较明显，另一个方向上的梯度幅值变化较弱，其余部分的点都还是集中原点附近，这样Harris矩阵的两个特征值理论应该是一个比较大，一个比较小，当然对于边缘，某些情况下致使计算出的特征值并不是都特别的大，但仍跟含有角点的窗口的分布情况具有不同。 Harris角点检测的性质：旋转不变性，光照不变性、比度变化部分不变性；同时不具有尺度不变性

2. Sobel滤波

sobel 算子是一个主要用做边缘检测的离散微分算子，Sobel算子结合了高斯平滑和微分求导，用来计算图像灰度函数的近似梯度。在图像的任何一点使用此算子，将会产生对应的梯度矢量或是其法矢量。这本次实践中，用sobel滤波器计算梯度。

3.角点响应函数R

定义了角点响应函数R ，通过判定R 大小来判断窗口是否有角点，即 R > threshold 时则有角点（threshold是我们自定义的一个阈值。）
角点应该满足基本性质：窗口最小的特征值尽量的大，此时R =λmin ；当然还有更有效的响应函数，计算公式为：R = det(M)−k(trace(M))2 其中的k是超参数
（一般的，增大k的值，降低角点检测的灵敏度，减少被检测角点的数量；减少k值，增加角点检测的灵敏度，增加被检测角点的数量。）

二、代码实现

1. main函数

main函数主要功能对图像进行读取，并进行灰度图转化，同时调用corner_detect函数，对检测结果进行可视化展示

def main():
    # 读取图片
    image = cv2.imread("./test1.jpg")
    height, width, channel = image.shape
    cv2.imshow('image', image)
    cv2.waitKey(1000)
    cv2.destroyAllWindows()

    # 灰度图
    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
	# 调用corner_detect函数计算角点
    dst = corner_detect(gray,block_size=3, nms_size=3, method='harris', k=0.04
)
	
	# 直接用opencv中的角点检测函数  dst = cv2.cornerHarris(gray, blockSize=3, ksize=3,k=0.05) 
    image_dst = image[:, :, :]
    image_dst[dst > 0.01 * dst.max()] = [0, 0, 255]
    #可视化结果
    cv2.imshow('dst', image_dst)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

2. corner_detect函数

主要功能对输入灰度图进行角点检测，并返回角点检测的结果corner

# 对输入图像进行Sobel滤波并计算梯度、计算每个像素的H矩阵、计算角点检测响应图R。
def corner_detect(img, block_size, nms_size, method, k):
    """
    img:待检测的灰度图
    block_size:领域范围
    nms_size：非最大值抑制的邻域范围
    method:检测角点的方法
    k：超参数
    """
    H, W = img.shape[:2]

    # 1.sobel滤波
    pad = block_size // 2  # pad=1
    out = np.zeros((H + pad * 2, W + pad * 2), dtype=float)
    out[pad:pad + H, pad:pad + W] = img.copy().astype(float)

    tmp = out.copy()
    out_v = out.copy()
    out_h = out.copy()

    K_v = [[1, 2, 1], [0, 0, 0], [-1, -2, -1]]
    K_h = [[1, 0, -1], [2, 0, -2], [1, 0, -1]]
    for y in range(H):
        for x in range(W):
            out_v[pad + y, pad + x] = sum(K_v * (tmp[y:y + block_size, x:x + block_size]))
            out_h[pad + y, pad + x] = sum(K_h * (tmp[y:y + block_size, x:x + block_size]))
    out_v = np.clip(out_v, 0, 255)
    out_h = np.clip(out_h, 0, 255)
    out_v = out_v[pad:pad + H, pad:pad + W].astype(np.uint8)
    out_h = out_h[pad:pad + H, pad:pad + W].astype(np.uint8)

    # 2.计算梯度
    grad = np.zeros((H, W, 2), dtype=float)
    grad[:, :, 0] = out_v
    grad[:, :, 1] = out_h
    # 3.计算H矩阵
    H_m = np.zeros((H, W, 3), dtype=float)
    H_m[:, :, 0] = grad[:, :, 0] ** 2
    H_m[:, :, 1] = grad[:, :, 1] ** 2
    H_m[:, :, 2] = grad[:, :, 0] * grad[:, :, 1]
    # 取最小特征值作为检测值，记在H_e矩阵
    H_e = np.zeros((H, W), dtype=float)
    for i in range(H):
        for j in range(W):
            if H_m[i, j, 0] > H_m[i, j, 1]:
                H_e[i, j] = H_m[i, j, 1]
            else:
                H_e[i, j] = H_m[i, j, 0]

    # 4.计算角点检测响应图R
    H_m = [np.array([[H_m[i, j, 0], H_m[i, j, 2]], [H_m[i, j, 2], H_m[i, j, 1]]]) for i in range(H) for j in range(W)]
    D, T = list(map(np.linalg.det, H_m)), list(map(np.trace, H_m))
    R = np.array([abs(d - k * t ** 2) for d, t in zip(D, T)])

    # 5.将计算出响应函数的值R进行非极大值抑制，滤除一些不是角点的点，同时要满足大于设定的阈值
    R_max = np.max(R)
    R = R.reshape(H, W)
    H_max = np.max(H_e)
    corner = np.zeros_like(R, dtype=uint8)

    for i in range(H):
        for j in range(W):
            if method == 'harris':
                # 角点检测器，以及对3x3邻域内非极大值进行抑制
                if R[i, j] > R_max * 0.01 and R[i, j] == np.max(
                        R[max(0, i - 1):min(i + nms_size - 1, H - 1), max(0, j - 1):min(j + nms_size - 1, W - 1)]):
                    corner[i, j] = 255
            else:
                # 依据H矩阵最小特征值进行检测
                if H_e[i, j] > H_max * 0.01:
                    corner[i, j] = 255

    cv2.imshow('corner', corner)
    # cv2.waitKey(0)
    return corner

三、结果对比

1.不同检测方法（对比一下使用角点检测器R与H矩阵最小特征值）

图1为角点检测器R，图2为只使用最小特征值，显然使用检测器R的效果不如最小特征值，检测出的角点更多倾向于边缘线条而不是单独的点

2.不同超参数k（使用角点检测器R，改变k值）

图3为k=0.05，图4为k=0.03，可知增大k的值，降低角点检测的灵敏度，减少被检测角点的数量；减少k值，增加角点检测的灵敏度，增加被检测角点的数量。

3.非最大值抑制（对3x3的邻域范围内非最大值抑制）

加入之后能解决角点检测器R的线条问题，对于检测连续的点进行抑制，角点效果更好
（几乎看不出原图的美女有角点）

四、参考资料

1. CSDN博主「果壳乄」的文章《角点特征之Harris角点检测算法》原文链接：
   https://blog.csdn.net/my_kun/article/details/106918857

2. CSDN博主「weixin_39745345」的文章《python sobel滤波_Sobel滤波器》原文链接：https://blog.csdn.net/weixin_39745345/article/details/111515227

3. 博客园博主wyu123《学习OpenCV中边缘检测的各种算子和滤波器——Canny算子,Sobel算子,Laplace算子以及Scharr滤波器》原文链接：https://www.cnblogs.com/wyuzl/p/6286083.html