Harris角点检测

角点的定义：

角点检测算法基本思想是使用一个固定窗口（取某个像素的一个邻域窗口）在图像上进行任意方向上的滑动，比较滑动前与滑动后两种情况，窗口中的像素灰度变化程度，如果存在任意方向上的滑动，都有着较大灰度变化，那么我们可以认为该窗口中存在角点。

角点：最直观的印象就是在水平和竖直两个方向变化均较大的两个点，即X,Y方向梯度都较大

边缘：仅在水平或竖直方向有较大的变化量，即X、Y方向梯度只有一个较大

平坦区域：在水平和竖直方向的变化量均较小，即X、Y方向梯度都较小

在特征点检测中经常提出尺度不变、旋转不变、抗噪声影响等，这些是判断特征点是否稳定的指标。

性能较好的角点：

1. 检测出图像中“真实”的角点
2. 准确的定位性能
3. 很高的重复检测率
4. 噪声的鲁棒性
5. 较高的计算效率

今天看到角点检测的时候真的是一脸蒙蔽，服了，书上直接给出的公式。整的我直接懵了。参考晚上其他大佬写的博客后逐渐明白了一点

书上给的公式如下：
$E(u,v)=\sum_{x,y}w(x,y) \left[] \begin{matrix} Ix^2 & IxIy\\ IxIy & Iy^2 \end{matrix} \right] \tag{B}$E(u,v)=x,y∑​w(x,y)[]Ix2IxIy​IxIyIy2​](B)

w(x,y)是窗口函数，最简单情形就是窗口WW内的所有像素所对应的w权重系数均为1.但有时候，我们会将w(x,y)函数设置为以窗口WW中心为原点的二元正态分布。如果窗口WW中心点是角点时，移动前与移动后，该点在灰度变化贡献最大；而离窗口WW中心(角点)较远的点，这些点的灰度变化几*缓，这些点的权重系数，可以设定小值，以示该点对灰度变化贡献较小，那么我们自然而然想到使用二元高斯函数来表示窗口函数；

Harris矩阵推导过程:

使用一个固定窗口在图像上进行任意方向上的滑动。滑动分别为u和v.

比较滑动前与滑动后俩种情况。所以可以得到以下公式：
$E(u,v)=\sum_{x,y}w(x,y)[I(x+u,y+v)-I(x,y)]^2$E(u,v)=x,y∑​w(x,y)[I(x+u,y+v)−I(x,y)]2

高数中将这个式子进行泰勒展开：
$f(x+u,y+v)\approx f(x,y)+ufx(x,y)+vfy(x,y)$f(x+u,y+v)≈f(x,y)+ufx(x,y)+vfy(x,y)

$所以\sum[I(x+u,y+v)-I(x,y)]^2=\sum[I(x,y)+uIx+vIy-I(x,y)]^2$所以∑[I(x+u,y+v)−I(x,y)]2=∑[I(x,y)+uIx+vIy−I(x,y)]2
$=\\\sum(u^2Ix^2)+2uvIxIy+v^2Iy^2=\\\sum[u,v] \left[ \begin{matrix} Ix^2 & IxIy\\ IxIy & Iy^2 \end{matrix} \right] \tag{B} \left[ \begin{matrix} u\\ v\\ \end{matrix} \right]$=∑(u2Ix2)+2uvIxIy+v2Iy2=∑[u,v][Ix2IxIy​IxIyIy2​][uv​](B)

• [u,v][u,v]是窗口WW的偏移量；
• (x,y)(x,y)是窗口WW所对应的像素坐标位置，窗口有多大，就有多少个位置；
• I(x,y)是像素坐标位置(x,y)的图像灰度值；
• I(x+u,y+v)是像素坐标位置(x+u,y+v)的图像灰度值；

Ix,Iy分别为窗口内像素点(x,y)在x方向上和y方向上的梯度值。 其中**窗口函数（权重矩阵）**可以是平坦的，通常为高斯滤波器Gσ：

设:
$Mi =\left[ \begin{matrix} Ix^2 & IxIy\\ IxIy & Iy^2 \end{matrix} \right]\\ M = \sum_{x,y}w(x,y) * Mi$Mi=[Ix2IxIy​IxIyIy2​]M=x,y∑​w(x,y)∗Mi

所以可以得到卷积：M

该卷积的目的是得到MI在周围像素上的局部平均。**矩阵M又称为Harris矩阵。**W 的宽度决定了在像素x 周围的感兴趣区域。

根据上述表达式，当窗口在平坦区域上移动，可以想象得到，灰度不会发生什么变换。E(u,v)=0；如果窗口处在纹理比较丰富的区域上滑动，那么灰度变化会很大。算法最终思想就是计算灰度发生较大变化时所对应的位置，当然这个较大是指任意方向上的滑动，并非单指某个方向。

接下来我们对矩阵M的特征值进行分析。设矩阵M的特征值为λ1,λ2.

则有以下三种情况：

具体为什么有这三种情况以下就是推导过程：

我们可以将E(u,v)近似为二项函数：

​
$E(u,v)=Au^2+2Cuv+Bv^2=\left[ \begin{matrix} A & C\\ C & A\\ \end{matrix} \right]\\ A=\sum_{x,y}w(x,y)*Ix^2\\ B=\sum_{x,y}w(x,y)*Iy^2\\ C=\sum_{x,y}w(x,y)*IxIy$E(u,v)=Au2+2Cuv+Bv2=[AC​CA​]A=x,y∑​w(x,y)∗Ix2B=x,y∑​w(x,y)∗Iy2C=x,y∑​w(x,y)∗IxIy
二次项函数本质上就是一个椭圆函数。椭圆的长和宽是由MM的特征值λ1,λ2决定的(椭圆的长短轴正是矩阵MM特征值平方根的倒数)，椭圆的方向是由M的特征向量决定的，椭圆方程为：
$\left[ \begin{matrix} u & v\\ \end{matrix} \right]*M*\left[ \begin{matrix} u\\ v\\ \end{matrix} \right]=1$[u​v​]∗M∗[uv​]=1
虽然我们利用E(u,v)来描述角点的基本思想，然而最终我们仅仅使用的是矩阵M。让我们看看矩阵M形式，是不是跟协方差矩阵形式很像，像归像，但是还是有些不同，哪儿不同？一般协方差矩阵对应维的随机变量需要减去该维随机变量的均值：

把Ix,Iy看成两个字段，假设窗口内有m个像素点，也就是等价于有m个样本，我们先计算每个字段的均值：
$\overline{Ix}=\sum_{i}^mIxi\\ \overline{Iy}=\sum_{i}^mIyi$Ix=i∑m​IxiIy​=i∑m​Iyi
我们仍然使用(Ixi,Iyi)表示样本(Ixi,Iyi)去均值后的值，则由这m个样本组成的矩阵：
$x=\left[ \begin{matrix} Ix1 & Ix2 & ... & Ixm \\ Iy1 & Iy2 & ... & Iym\\ \end{matrix} \right]$x=[Ix1Iy1​Ix2Iy2​......​IxmIym​]
则对应协方差矩阵可以写成：
$C=\frac{1}{m}*X*X^T=\frac{1}{m}*\sum_{i=1}\left[ \begin{matrix} Ix^2 & IxIy\\ IxIy & Iy^2 \end{matrix} \right]$C=m1​∗X∗XT=m1​∗i=1∑​[Ix2IxIy​IxIyIy2​]
但矩阵M中并没有这样做，所以在矩阵M里，我们先进行各维的均值化处理，那么各维所对应的随机变量的均值为0，协方差矩阵就大大简化了，简化的最终结果就是矩阵M，是否明白了(注意为了简化运算，我们先假设M矩阵中的权重系数w(x,y)=1，并且省略掉了求均值.我们的目的是分析数据的主要成分

如果我们对协方差矩阵M进行对角化，很明显，其对角线就是各个字段的方差，这点大家应该明白吧？不明白的话可以复习下PCA原理。

• 如果两个字段(Ix,Iy)所对应的特征值都比较大，说明什么？ 像素点的梯度分布比较散，梯度变化程度比较大，符合角点在窗口区域的特点；
• 如果是平坦区域，那么像素点的梯度所构成的点集比较集中在原点附近，因为窗口区域内的像素点的梯度幅值非常小，此时矩阵M的对角化的两个特征值比较小；
• 如果是边缘区域，在计算像素点的x、y方向上的梯度时，边缘上的像素点的某个方向的梯度幅值变化比较明显，另一个方向上的梯度幅值变化较弱，其余部分的点都还是集中原点附近，这样MM对角化后的两个特征值理论应该是一个比较大，一个比较小，当然对于边缘这种情况，可能是呈45°的边缘，致使计算出的特征值并不是都特别的大，总之跟含有角点的窗口的分布情况还是不同的。

因此可以得出下列结论：

• 特征值都比较大时，即窗口中含有角点；

• 特征值一个较大，一个较小，窗口中含有边缘；

• 特征值都比较小，窗口处在平坦区域；

  那么根据书上又提出了一个不用实际计算特征值的方法。也就是度量角点响应：

  对每一个窗口计算得到一个分数R，根据R的大小来判定窗口内是否存在harris特征角。分数R根据下面公式计算得到：
  $R=det(M)-k(trace(M))^2\\ 其中：det(M)=\lambda1*\lambda2\\ trace(M)=\lambda1+\lambda2\\ 另外我在别的地方看到简化后的R为： R=[(Ix^2+Iy^2)-(IxIy)^2]/(Ix^2+Iy^2)$R=det(M)−k(trace(M))2其中：det(M)=λ1∗λ2trace(M)=λ1+λ2另外我在别的地方看到简化后的R为：R=[(Ix2+Iy2)−(IxIy)2]/(Ix2+Iy2)

得到Harris角点的步骤

所以得到Harris角点的步骤为：

1、计算图像I(x,y)I(x,y)在x和y两个方向的梯度Ix,Iy;

2、计算图像两个方向梯度的乘积;

3、使用高斯函数对第二步算出的数进行高斯加权，计算中心点为(x,y)的窗口W对应的矩阵M

4、计算每个像素点(x,y)处的(x,y)处的Harris响应值R；

5：过滤大于某一阈值t的R值；

代码实现

python代码如下：

from scipy.ndimage import filters
from PIL import Image
from pylab import *
from imtools import imresize
np.set_printoptions(threshold=np.inf)#解决printf打印矩阵不完全的问题
def harris_response(im,sigma=3):
    """计算图像的harris响应函数"""
    #计算导数
    imx = zeros(im.shape)
    imx = filters.gaussian_filter(im,(sigma,sigma),(0,1),imx)
    imy = zeros(im.shape)
    imy = 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)
    #计算像素的角点响应函数
    return (wxx*wyy - 2*wxy)/(wxx + wyy)
def get_harris_points(harrism,min_dist = 10,thresold = 0.1):
    """从一幅Harrisim响应中返回角点，min_dist为分割角点和图像边界的最少像素数目"""
    corner_thsold = harrism.max()*thresold
    harrism_t = (harrism > corner_thsold) * 1
    #得到候选点的坐标
    coords = array(harrism_t.nonzero()).T#返回非零值的坐标的矩阵
    #他们的Harris响应值
    candidate_values = [harrism[c[0],c[1]] for c in coords]
    #对候选点进行harris响应值进行排序
    index = argsort(candidate_values)[::-1]#将x中的元素从小到大排列，提取其对应的index(索引)，然后输出到y
    #将可行点的位置保存在数组里
    allowed_locations = zeros(harrism.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
    #按照min_distance原则，选择最佳harris点
    filters_coords = []
    for i in index:
            if allowed_locations[coords[i,0],coords[i,1]] == 1:
                filters_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 filters_coords