Harris角点检测器

1. Moravec角点检测算子
Moravec角点检测的思想很简单，在图像上取$W*W$的滑动窗口，不断移动窗口并检测窗口中像素变化情况E。像素变化情况可分三种：(1)如果在窗口中图像是平坦的，那么E的变化不大；(2)如果在窗口中的图像是一条边，那么在沿这条边滑动时E变化不大，而在沿垂直于这条边滑动时E变化很大；(3)如果在窗口中图像是一个角点，窗口沿任何方向移动E都会发生很大变化，即

$$\begin{equation*} E_{x,y}=\sum_{u,v}{w_{u,v}{\vert{I_{x+u,y+v}-I_{u,v}}\vert}^2} \end{equation*}$$
其中，(x,y)就表示八个移动方向( $\pm1,0$),( $0,\pm1$),( $\pm1,\pm1$),E是像素的变化值。

2. Harris角点检测算子
Harris角点检测算子实际上是对Moravec算子的改良，在原文中，作者提出了三点Moravec算子的缺陷并且给出了改良方法：
(1)Moravec算子对方向依赖性太强，Moravec算子实际上只移动了水平，垂直，和对角方向，为此，作者采用微分的思想：

$$\begin{align} E_{x,y}=\sum_{u,v}{{\vert{I_{x+u,y+v}-I_{u,v}}\vert}^2} =\sum_{u,v}{w_{u,v}[xX+yY+o(x^2,y^2)]^2} \end{align}$$
其中，
$$\begin{equation} X = I \ast(-1,0,1) \approx \frac{\partial{I}}{\partial{x}} Y = I \ast(-1,0,1)^T \approx \frac{\partial{I}}{\partial{y}} \end{equation}$$
所以，E就可表示为：
$$\begin{equation} E(x,y)=Ax^2+2Cxy+By^2 \end{equation}$$
其中，
$$\begin{aligned} A&=X^2\ast w\\ B&=Y^2\ast w\\ C&=XY\ast w \end{aligned}$$

(2)Moravec算子采用的是方形的窗口，E的响应比较容易受到干扰，Harris采用高斯核函作为平滑的窗口

$$\begin{aligned} w_{u,v}=e^{{-\frac{u^2+v^2}{2\sigma^2}}} \end{aligned}$$

(3)Moravec算子对于边缘响应过于灵敏(这个不是很明白作者的意思)

$$\begin{equation} E(x,y)=(x,y)M(x,y)^T \end{equation}$$
其中，
$$\begin{equation} M=\left(\begin{array}{cc}A&C\\C&B\end{array}\right) \end{equation}$$
用 $\alpha,\beta$表示矩阵的特征值，A.如果 $\alpha,\beta$都很小，说明高斯windows中图像接近平坦；B.如果一个大一个小，则表示检测到边；C. 如果都很大，那么表示检测到了角点。

Harris定义了一个响应函数来表示$\alpha,\beta$的大小关系：

$$\begin{equation} R=\alpha\beta-k(\alpha+\beta)^2=det(M)-kTr(M) \end{equation}$$
k一般是一个远小于1的数，opencv的默认推荐值是0.04
R为正值时，检测到的为角点，R为负值时检测到的是边，R很小时检测到的是平坦区域

Harris角点检测算法具有旋转不变形，不具备尺度不变性

# congding:utf-8

from scipy.ndimage import filters
import numpy as np
from matplotlib.pyplot import *

def compute_harris_response(im,sigma=3):
    '''在一幅灰度图像中，对每个像素计算Harris 角点检测器响应函数'''
    # 计算导数
    imx = np.zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (0,1), imx)
    imy = np.zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (1,0), imy)
    # 计算Harris 矩阵的分量
    Wxx = filters.gaussian_filter(imx*imx,sigma)
    Wxy = filters.gaussian_filter(imx*imy,sigma)
    Wyy = filters.gaussian_filter(imy*imy,sigma)
    # 计算特征值和迹
    Wdet = Wxx*Wyy - Wxy**2
    Wtr = Wxx + Wyy
    return Wdet / Wtr

def get_harris_points(harrisim,min_dist=10,threshold=0.1):
    """ 从一幅Harris 响应图像中返回角点。min_dist 为分割角点和图像边界的最少像素数目"""
    # 寻找高于阈值的候选角点
    corner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim > corner_threshold) * 1
    # 得到候选点的坐标
    coords = np.array(harrisim_t.nonzero()).T
    # 以及它们的Harris 响应值
    candidate_values = [harrisim[c[0],c[1]] for c in coords]
    # 对候选点按照Harris 响应值进行排序
    index = np.argsort(candidate_values)
    # 将可行点的位置保存到数组中
    allowed_locations = np.zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
    # 按照min_distance 原则，选择最佳Harris 点
    filtered_coords = []
    for i in index:
        if allowed_locations[coords[i,0],coords[i,1]] == 1:
            filtered_coords.append(coords[i])
            allowed_locations[(coords[i,0]-min_dist):(coords[i,0]+min_dist),(coords[i,1]-min_dist):(coords[i,1]+min_dist)] = 0
    return filtered_coords

def plot_harris_points(image,filtered_coords):
    """ 绘制图像中检测到的角点"""
    figure()
    gray()
    imshow(image)
    plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords],'*')
    axis('off')
    show()