OpenCV --Harris corner detection

为什么需要角点

在现实世界中,要想提取目标物体的主要特征(如位置特征),最简单的可以使用模板匹配的方法进行像素级别的比较。但是模板匹配对于模板的大小要求是非常严苛的,因此又可以考虑通过边缘检测来直接提取方位。对于直线或者圆形的提取可以使用Hough Transform的投票算法,甚至对应任意多边形还可以使用Generalized Hough Transform方法。更进一步,如果物体的角度发生变化或者同一个物体的不同部分进行匹配又当如何呢?我们就需要提取特征点(feature)。特征点通常具备独立性(一个角一个点),局部性(粒度足够细)等特点,而角点也是特征点的一种,主要体现的也是一种位置信息。(密集恐惧症患者请撤离现场!算了已经晚了)
微信截图_20200528162231.png


角点的定义

  1. 角点是具有两个边缘的特征点
  2. 在角点处有多于一个主方向

角点的研究方法

  1. 通过图像边缘编码,对图像分割,提取边缘。但是如果局部的目标发生变化就会出现较大的问题。这种研究方法早期有Rosenfeld和Freeman等人的研究,后期主要是CSS等方法。
  2. 通过梯度计算提取边缘来避免第一类缺陷。那么就需要设立响应的算子。主要有Moravec算子、Forstner算子、Harris算子和SUSAN算子等等。本文主要记录Harris算子的研究心得。

Harris角点

基本原理

Harris角点在二维图像的基本特征在于两个方向上都存在梯度。因此位于边缘和角点的梯度特征应该是不同的。
081646497655888.png

假设在某像素位置的窗口(x, y)移动(u, v)
Harris据此提出了一个Error结构:
1.jpg
由:
2.jpg
得到:
3.jpg
从上面的表述很容易看出中间的(2*2)的矩阵是Hessian矩阵。其中如果$\lambda_1$、$\lambda_2$是矩阵的两个特征值,那么就代表着图像的$x$、$y$两个方向上的梯度变化。因此,如果:

  1. $\lambda_1$、$\lambda_2$都比较小:说明是非边缘的中间的区域。
  2. $\lambda_1$ >> $\lambda_2$或者$\lambda_1$ << $\lambda_2$:边缘区域。
  3. $\lambda_1$、$\lambda_2$都比较大:图像的角点。

因此Harris定义了一个响应值来定义角点的特征:
4.jpg

Harris算法过程

  1. 对图像使用梯度算子,并计算$I_x2$、$I_y2$、$I_xy$,在过程中还发现了对梯度算子解释的非常好的一篇博文。文章里面将梯度的计算归类为距离的定义,有连续距离定义方式(欧氏距离)和离散距离定义方式(城市距离)关于算子的解释
  2. 对梯度图像进行模糊,去除噪声干扰。
  3. 计算每个点的响应值$R$。
  4. 调试阈值参数$T$,选择高过$T$的响应值点就是角点。

opencv-python实现方法

import cv2
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline

img_chess = cv2.imread('./chessboard.jpg')		## chessboard image
img_chess_gray = cv2.cvtColor(img_chess, cv2.COLOR_BGR2GRAY)
img_chess_gray = np.float32(img_chess_gray)
dst = cv2.cornerHarris(img_chess_gray, 2, 5, 0.04)		## 参数表(src, block_size, k_size, k)分别是float32类型图像,在求E过程中的邻域大小,Sobel算子kernel大小和Harris equation的自由参数(alpha)
dst = cv2.dilate(dst,None)		## 对点做膨胀,不是必要步骤

img = img_chess.copy()
img[dst > 0.01 * dst.max()] = [255, 0, 0]
plt.imshow(img)

调包可能也就是一行的问题...下面是自己对Harris Corner Detection的检测算法实现,非极大值抑制部分写的不是很高效。

def cornerHarris(img_chess_gray, alpha=0.06):
#     fx = np.array([[5, 0, -5],
#           [8, 0, -8], 
#           [5, 0, -5]])
    fx = np.array([[1, 0, -1],		## Sobel Gx
                   [2, 0, -2],
                  [1, 0, -1]])
#     fy = np.array([[5, 8, 5],
#           [0, 0, 0],
#           [-5, -8, -5]])
    fy = np.array([[1, 2, 1], 		## Sobel Gy
                   [0, 0, 0], 
                   [-1, -2, -1]])
    if img_chess_gray is None:
        return -1
    else:
        img_gray = img_chess_gray.copy()
    Ix = np.zeros_like(img_gray)
    Iy = np.zeros_like(img_gray)
    Ix = cv2.filter2D(img_gray, -1, fx)
    Iy = cv2.filter2D(img_gray, -1, fy)
    Ix2 = np.multiply(Ix, Ix)
    Iy2 = np.multiply(Iy, Iy)
    Ixy = np.multiply(Ix, Iy)     # cross derivate
    
    ## Blur
    Ix2 = cv2.filter2D(Ix2, -1, cv2.getGaussianKernel(5, 2))		## 图像噪声比较大,提前处理一下,对实验结果影响比较大
    Iy2 = cv2.filter2D(Iy2, -1, cv2.getGaussianKernel(5, 2))
    Ixy = cv2.filter2D(Ixy, -1, cv2.getGaussianKernel(5, 2))
    result = np.zeros_like(img_gray)     # record the corner point
    lamb = np.zeros((img_gray.shape[0]*img_gray.shape[1], 2))
    R = np.zeros_like(img_gray)
    
    ## calculate R
    k = 0
    for i in range(img_gray.shape[0]):
        for j in range(img_gray.shape[1]):
            M = np.array([[Ix2[i, j], Ixy[i, j]], [Ixy[i, j], Iy2[i, j]]])
            K = np.linalg.det(M)
            H = np.trace(M)
            R[i, j] = K - alpha * H**2
            lamb[k, :] = [K, H]
            k += 1
#     plt.scatter(lamb[:, 0], lamb[:, 1])
#     plt.xlabel('det')
#     plt.ylabel('trace')

    ## NMS
#     R_max = np.max(R)
    for i in range(2, img_gray.shape[0]-1):
        for j in range(2, img_gray.shape[1]-1):
            if (R[i, j] > 0) and (R[i, j] >= np.max(R[i-1:i+2, j-1:j+2])):
                result[i, j] = 1
    img2 = img_chess.copy()
    img2[result==1] = [255, 0, 0]
    plt.imshow(img2)

参考文献