Python 实现Harris 角点检测

具体原理见Feature extraction & image processing for computer vision (3rd edition) 4.4.1.4小节

# Harris Corner Detector
# Code - 4.18
# -*- coding:gb2312 -*-
import cv2
import numpy as np
import matplotlib.pyplot as plt
from skimage.feature import canny
# from scipy.ndimage import *
def Harris(img,op='H',w=4,k=0.1):
    # below is wrong for trying to create a
    # ndarray with demision number greater than 32
    # img = np.ndarray(img, dtype=np.float32)
    # if use np.mat, be ware of the difference in matrix multiply
    img = np.array(img,dtype=np.float32)
    # row,col = img.shape
    # out = np.zeros(img.shape)
    difx,dify=np.gradient(img)

    # detect edge with cv2.Sobel
    # EdgeX = cv2.Sobel(img, cv2.CV_32FC1, 1, 0)
    # EdgeY = cv2.Sobel(img, cv2.CV_32FC1, 0, 1)
    # Mag = np.sqrt(EdgeX ** 2 + EdgeY ** 2)

    # or use canny edge detector from skimage.feature
    Mag = canny(img, 1, 0.4, 0.2)

    # compute autocorralation
    difx2=difx**2
    dify2=dify**2
    difxy=difx*dify

    # or use cv2.multiply
    # difx2 = cv2.multiply(difx,difx)
    # dify2 = cv2.multiply(dify,dify)
    # difxy = cv2.multiply(difx, dify)

    # mean filter in scipy.ndimage
    # A = uniform_filter(difx2,size=w)
    # B = uniform_filter(dify2,size=w)
    # C = uniform_filter(difxy,size=w)

    # or use mean filter in cv2
    A = cv2.blur(difx2,(w,w))
    B = cv2.blur(dify2,(w,w))
    C = cv2.blur(difxy,(w,w))

    if op =='H':
        out = A*B - C**2 -k*(A+B)**2
        out[Mag == 0] = 0
    else:
        out = (A*dify2-2*C*difxy+B*difx2)/(difx2+dify2+1)
        out[difx2+dify2==0]=0

    # next section for debug
    # plt.subplot(221);plt.imshow(img)
    # plt.subplot(222);plt.imshow(difx)
    # plt.subplot(223);plt.imshow(difx2);plt.show()
    # # plt.subplot(224);plt.imshow(C);plt.show()
    return out

if __name__ == '__main__':
    image = cv2.imread('u.tif', flags=0)
    image = cv2.resize(image, (100, 100))
    corimg = Harris(image, op='TI')
    plt.subplot(121);plt.imshow(image);plt.title('Original')
    plt.subplot(122);plt.imshow(corimg);plt.title('Harris Corner')
    plt.show()

图像处理结果：

有几个问题：
1. numpy.array 和numpy.matrix的矩阵运算符定义不同。前者*表示矩阵点乘，后者表示矩阵相乘。
2. 计算梯度时，需要将图像数据转换成有符号数据类型。因为梯度有可能是负数。
3. Harris 算法的二次方程的旋转蕴含了主元分析的思想。