【局部图像描述子】Harris角点检测器

声明：笔记来源于【python计算机视觉】一书。

什么是角点？
1. 角点是一阶导数（即灰度的梯度）的局部最大所对应的像素点；
2. 角点是两条及两条以上边缘的交点;
3. 角点指示了物体边缘变化不连续的方向;
4. 角点处的一阶导数最大，二阶导数为0；
5. 角点是指图像中梯度值和梯度方向的变化速率都很高的点。
  
  来源：https://wenku.baidu.com/view/f61bc369561252d380eb6ef0.html
Harris角点检测算法

数学推导

设图像窗口平移量为（u,v）,产生的灰度变化为E(u,v),
有

E(u,v)=sum[w(x,y)[I(x+u,y+v)−I(x,y)]2](1)

E(u,v)=sum[w(x,y)[I(x+u,y+v)-I(x,y)]^2]\tag{1}

w(x,y)为窗口函数
I(x+u,y+v)为平移后的灰度值
I(x,y)为平移前的灰度值

有泰勒公式展开可得：

I(x+u,y+v)=I(x,y)+Ix∗u+Iy∗v+O(u2,v2)(2)

I(x+u,y+v)=I(x,y)+Ix*u+Iy*v+O(u^2,v^2) \tag{2}

Ix,Iy 分别为偏导数，在图像中为图像的方向导数。

E(u,v)=sum[w(x,y)[Ix∗u+Iy∗v+O(u2,v2)]2](3)

E(u,v)=sum[w(x,y) [Ix*u+Iy*v+O(u^2,v^2)]^2]\tag{3}

可以近似得到

E(u,v)=sum[w(x,y)[Ix∗u+Iy∗v]2](4)

E(u,v)=sum[w(x,y) [Ix*u+Iy*v]^2]\tag{4}

即

E(u,v)=[u,v][Ix2Ix∗IyIx∗IyIy2][uv](5)

E(u,v)=[u,v] \left[\begin{matrix}Ix^2 & Ix*Iy \\Ix*Iy & Iy^2\end{matrix}\right]\left[\begin{matrix}u \\v \end{matrix}\right]\tag{5}
令

M=[Ix2Ix∗IyIx∗IyIy2](6)

M = \left[\begin{matrix}Ix^2 & Ix*Iy \\Ix*Iy & Iy^2\end{matrix}\right]\tag{6}

因此最后对角点的检测成了对矩阵M的特征值的分析了，令M其特征值为x1,x2;

当x1>>x2或者x2>>x1，则检测到的是边缘部分；
当x1,x2都很小，图像窗口在所有移动的方向上移动灰度级都无明显变化.
当X1,X2都很大时且相当，检测到的是角点。

编程时用x1,x2不方便，因此定义角点响应函数;

R=det(M)−k(trace(M))2(7)

R=det(M)-k(trace(M))^2\tag{7}

其中det(M)为矩阵M的行列式，trace(M)为矩阵M的迹。

更具体数学公式实际编程的步骤

利用水平，竖直差分算子对图像的每个像素进行滤波以求得Ix,Iy,进而求得M中的四个元素的值。

M=[Ix2Ix∗IyIx∗IyIy2](6)

M = \left[\begin{matrix}Ix^2 & Ix*Iy \\Ix*Iy & Iy^2\end{matrix}\right]\tag{6}
对M的四个元素进行高斯平滑滤波，为的是消除一些不必要的孤立点和凸起，得到新的矩阵M。
接下来利用M计算对应每个像素的角点响应函数R，即：

R=det(M)−k(trace(M))2(7)

R=det(M)-k(trace(M))^2\tag{7}
也可以使用改进的R：

R=（Ix2∗Iy2−(Ix∗Iy)2])/(Ix2+Iy2)(8)

R=（Ix^2*Iy^2-(Ix*Iy)^2])/(Ix^2+Iy^2)\tag{8}
【Python计算机视觉】书中选用的是后一种方法计算角点的响应函数。

里面没有随意给定的参数k，取值应当比第一个令人满意。
4.在矩阵R中，同时满足R(i,j)大于一定阈值threshold和R(i,j)是某领域内的局部极大值，则被认为是角点。

完整代码

# -*- coding:utf-8 -*-
from PIL import Image
from numpy import *
from pylab import *
from scipy.ndimage import filtersdef compute_harris_response(im, sigma=3):""" 在一幅灰度图像中，对每一个像素计算Harris角点检测器响应函数im:(数组图像）  sigma=3:标准差为3"""#计算x方向的导数imx = zeros(im.shape) #zeros()生成全零数组filters.gaussian_filter(im, (sigma,sigma),(0,1),imx)#计算x方向的导数imy = zeros(im.shape) #zeros()生成全零数组filters.gaussian_filter(im, (sigma,sigma),(0,1),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*2Wtr = Wxx + Wyyreturn Wdet / Wtr#获取角点
def get_harris_points(harrism, min_dist=10, threshold=0.1):"""从一幅Harris响应图像中返回角点。min_dist为分割角点和图像边界的最少像素数目"""#寻找高于阈值的候选角点corner_threshold = harrism.max()*thresholdharrisim_t = (harrisim > corner_threshold) * 1#得到候选点的坐标coords = array(harrisim_t.nonzero()).T#以及它们的Harris响应值candidate_values = [harrism[c[0],c[1]] for c in coords]#对侯选点按照Harris响应值进行排序index = argsort(candidate_values)#将可行点的位置保存到数组中allowed_locations = zeros(harrism.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)] = 0return filtered_coordsdef 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()if __name__ == "__main__":im = array(Image.open("/home/xuna/桌面/image/2.jpg").convert("L"))harrisim = compute_harris_response(im)filtered_coords = get_harris_points(harrisim,6)plot_harris_points(im, filtered_coords)

运行效果：

总结

数学，数学是最重要的，第一章还容易理解，到了第二章就有点不太懂了，harris算法学了一将近一周，书上说的知识感觉很朦胧，期间查了不少资料，打算先总结记录，这么经典的算法也不是一时可以理解透彻的，先慢慢消化。
参考资料：
https://wenku.baidu.com/view/89a54282a6c30c2258019e39.html?from=search
https://wenku.baidu.com/view/fafa9cf4ec3a87c24128c46c.html?from=search