图像局部特征描述子——Harris角点及其Python代码实现

前言：参考了《Python计算机视觉编程》一书，部分知识点理解有误请多谅解( • ̀ω•́ )✧

原理部分叙述

我们需要利用一定的准则来获取图片中的某些特征点，并获取特征点周围的图像块，最后进行图片的匹配，特征应具有几何、光度等的不变性，本文章介绍的是利用Harris算子进行特征获取。
角点（corner points）：
• 局部窗口沿各方向移动，均产生明显变化的点
• 图像局部曲率突变的点
以下是不同类型的角点

如何探查出角点？
假设有一块矩形区域在平面上移动，那么矩形内的图案有以下三种情况：

此时我们可以假设：
将图像窗口平移[u,v]产生灰度变化E(u,v)，那么：

函数I(x,y)表示图像的灰度信息，将I(x+u, y+v)函数在（x, y）处泰勒展开可以得到

最后，我们将E(u,v)用下面的式子表示：

并且，M等于：

记M的特征值为K1和K2，当K1和K2都比较小的时候可以认为是在平坦区域，K1或K2其中有一个远大于另一个值，则可能是在边缘区域，当K1和K2都比较大时，表示位于角点。
怎样确定K1和K2的大小关系？使用角点响应函数R

• R 只与M的特征值有关
• 角点：R 为大数值正数
• 边缘：R为大数值负数
• 平坦区：R为小数值

Python代码实现

我们需要建立一个文件用于写入角点响应函数，我这里设为Harris.py，另外需要另一个文件用来运行我们的命令，我设为test.py。
下面是Harris.py文件内容：

#<pre name="code" class="python">
#coding:utf-8
from scipy.ndimage import filters
from numpy import *
from pylab import *
def compute_harris_response(im,sigma=3):
    imx=zeros(im.shape)#计算导数
    filters.gaussian_filter(im,(sigma,sigma),(0,1),imx)
    imy=zeros(im.shape)
    filters.gaussian_filter(im,(sigma,sigma),(1,0),imy)
    Wxx=filters.gaussian_filter(imx*imx,sigma)
   #计算harris矩阵分量   
    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响应图像中返回角点"""
	conner_threshold=harrisim.max()*threshold
    harrisim_t=(harrisim>conner_threshold)*1
    
    coords=array(harrisim_t.nonzero()).T
    candidate_values=[harrisim[c[0],c[1]] for c in coords]
    index=argsort(candidate_values)
    allowed_locations=zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist]=1
    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#此处保证min_dist*min_dist只有一个harris特征点
    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()

接下来是测试文件test.py

from PIL import Image
 
from numpy import *
import harris
from pylab import *
from scipy.ndimage import filters
#改变图片路径，测试的时候路径里面有中文好像不可以
im=array(Image.open('C:/Users/lenovo/Desktop/house.jpg').convert('L')) 
harrisim=harris.compute_harris_response(im)
filtered_coords=harris.get_harris_points(harrisim)
harris.plot_harris_points(im,filtered_coords)

我的测试结果如下：