前言

提到图像匹配，就不得不提到SIFT(尺度不变特征变换) 算法，作为图像特征提取的典型算法，SIFT算法单单在论文上的引用就已经到了6万多次，更不用提工程项目上的使用次数了。不过在介绍经典算法之前，我先介绍下啥是特征点和一个简单的Harris角点检测算法。
本篇博客代码参考，来源于书本Python计算机视觉编程：https://github.com/jesolem/PCV

1.Harris角点检测算法

1.1.特征点

特征点就像我们平时辨认两个物品分别是什么的判断依据，一个物品是正立方体，一个物品是球，那么我们会说一个物品有棱角，一个物品很圆滑。但对于计算机来说，棱角和圆滑都是无法解析的词，于是计算机需要一个计算机来识别不同物品的方法，而起到让计算机分辨作用的点就是特征点。

1.2.Harris角点

很显然，不是每个像素点都能作为特征点的,特征点也不是仅用一个像素点就能概括的，为了便于理解，现在说的特征点仅指一个像素点，便于大家理解。
而角点顾名思义就是转角的点，比如三角形的三个顶点，正方体的八个顶点等。

1.3.Harris角点检测基本思想

在这里插入图片描述
如图所示，我们用一个四边形方框不断移动来判断这个正方形方框是否有巨大变化；第一幅图中，不管我们怎么移动，图像的的内容都不会有太大变换；第二幅图中，如果我们左右移动，图像的内容就会有比较明显的改变，但如果上下移动电话，图像又没有多大的变化；第三幅图中，不管我们往哪个方向移动，图像都会有明显变化。
我们要明确一点，图像的明显变化并不是我们看出来的，而是让计算机对比图像对应位置的像素值之差的累计之和得到的。

1.4.Harris角点检测数学原理

假设方框的大小为x*y的矩形，方框变化的距离为u，v,产生的灰度变化为E(u,v)，那么会有以下公式：
在这里插入图片描述
式子中的I(x,y)为坐标x,y的像素值，w(x,y)为每个像素值的权重，比如我们认为中间部分的像素值更重要，那么我们就把中间部分的w(x,y)调大。
使用泰勒公式展开化简后可以得到这个式子：

于是对于局部微小的移动量，我们可以得到以下式子
在这里插入图片描述

我们记M的特征值为λ1（y向）和λ2（x向），那么就有

如果λ1和λ2都很小，说明图像窗口在所有方向上移动都无明显灰度变化。

1.5.Harris焦点检测算法代码解析

还是以之前拍的照片为例
在这里插入图片描述

#在一幅灰度图像中，计算每一个像素点的Harris角点
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)

        # 可视化原图
    plt.subplot(2, 3, 1)
    plt.title('gray')
    plt.axis('off')
    plt.imshow(img, plt.cm.gray)

    # 可视化x方向导数图像
    plt.subplot(2, 3, 2)
    plt.title('Derivative in X direction')
    plt.axis('off')
    plt.imshow(imx, plt.cm.gray)

    # 可视化y方向导数图像
    plt.subplot(2, 3, 3)
    plt.title("Derivative in y direction")
    plt.axis('off')
    plt.imshow(imy, plt.cm.gray)

    plt.subplot(2, 3, 4)
    plt.title("Second derivative in X direction")
    plt.axis('off')
    plt.imshow(imx * imx, plt.cm.gray)

    plt.subplot(2, 3, 5)
    plt.title("Second derivative in X/Y direction")
    plt.axis('off')
    plt.imshow(imx * imy, plt.cm.gray)

    plt.subplot(2, 3, 6)
    plt.title("Second derivative in y direction")
    plt.axis('off')
    plt.imshow(imy * imy, plt.cm.gray)
    plt.show()

    # 计算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*Wtr)

在这里插入图片描述

#从一幅Harris响应图像中返回角点。min_dist为分隔角点和图像边界的最少像素数目
def get_harris_points(harrisim,min_dist=10,threshold=0.1):
    
    # 寻找高于阈值的候选角点
    corner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim > corner_threshold) * 1
    
    # 得到候选点的坐标
    coords = array(harrisim_t.nonzero()).T
    
    # 得到候选点的Harris响应值
    candidate_values = [harrisim[c[0],c[1]] for c in coords]
    
    # 对候选点按照Harris响应值进行排序
    index = argsort(candidate_values)
    
    # 将可行点的位置保存到数组中
    allowed_locations = 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()

这是get_harris_points(harrisim,6,0.1)的效果
在这里插入图片描述
检测的角点数目有点过多了，我们换成get_harris_points(harrisim,6,0.5)再试一次

还是骗多了，可能是照片本身偏暗了，换张亮的试试

threshold=0.1

threshold=0.5

仔细观察了一下，发现可能是因为最开始的照片像素值过大了，将其resize后
在这里插入图片描述
虽然角点还是过多，但至少能看出原来的样子了