文章目录

1、Harris角点检测

1.1 基本概念
1.2 数学表达

2、实现与测试

2.1 代码
2.2 运行结果
2.3 参数变化

2.3.1 其他不变，改变k（此时sigma=3）
2.3.2 其他不变，改变sigma（k=0.04）

1、Harris角点检测

1.1 基本概念

角点：窗口向任意方向的移动都导致图像灰度的明显变化。

在这里插入图片描述

1.2 数学表达

将图像窗口平移[u,v]产生灰度变化E(u,v)：

$E(u,v)=\sum_{x,y}^{ }w(x,y)[I(x+u,y+v)-I(x,y)]^{2}$
那么，如何求解平移后的图像灰度I(x+u,y+v)，以及灰度变化E(u,v)呢？

将I(x+u, y+v)函数在(x, y)处泰勒展开，得：
$I(x+u,y+v)=I(x,y)+I_{x}u+I_{y}v+O(u^{2},v^{2})$
则可求得灰度变化:
$E(u,v)=\sum_{x,y}^{ }w(x,y)[I_{x}u+I_{y}v+O(u^{2},v^{2})]^{2}\approx \sum_{x,y}^{ }w(x,y)[I_{x}u+I_{y}v]^{2}$
$[I_{x}u+I_{y}v]^{2}=[u,v]\begin{bmatrix} I_{x}^{2} & I_{x}I_{y}\\ I_{x}I_{y} & I_{y}^{2} \end{bmatrix}\begin{bmatrix} u\\ v \end{bmatrix}$
于是对于局部微小的移动量 [u,v]，可以近似得到下面的表达：
$E(u,v)\cong [u,v]M\begin{bmatrix} u\\v \end{bmatrix}$
其中M是 2*2 矩阵，可由图像的导数求得：
$M=\sum_{x,y}^{ }w(x,y)\begin{bmatrix} I_{x}^{2} &I_{x}I_{y} \\ I_{x}I_{y}& I_{y}^{2} \end{bmatrix}$
窗口移动导致的图像变化量：实对称矩阵M的特征值分析，记M的特征值为 $\lambda_{1}$ 、 $\lambda_{2}$ 。
且定义角点响应函数R:
$R=detM-k(traceM)^{2}$
$detM=\lambda _{1}\lambda _{2}$
$traceM=\lambda _{1}+\lambda _{2}$
$(k=0.04\sim 0.06)$

2、实现与测试

2.1 代码

# -*- coding: utf-8 -*-
from pylab import *
from PIL import Image
from PCV.localdescriptors import harris
from scipy.ndimage import filters

def compute_harris_response(im, sigma=3):
    # 在一幅灰度图像中，对每个像素计算Harris角点检测器响应函数
    
    # 计算导数
    k = 0.04
    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)
    
    # 计算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
    
    # 返回像素值为 Harris 响应函数值的一幅图像
    print(Wdet-(k*Wtr)) # 输出角点响应函数R
    return Wdet / Wtr  # 此处可消除参数k的影响
    
# 读入图像(读取图像到数组中)
im = array(Image.open('p3.jpg').convert('L'))  # 括号内：读取一幅图像，并将其转换成灰度图像

# 检测harris角点
harrisim = compute_harris_response(im)
                # compute_harris_response(im, sigma):
                #     在一幅灰度图像中，对每一个像素计算Harris角点检测器响应函数
                #     im:(数组图像）  sigma:标准差默认为3
                
# Harris响应函数
# print (harrisim)
harrisim1 = 255 - harrisim
figure()
gray()

#画出Harris响应图
subplot(141)
imshow(harrisim1)
print harrisim1.shape
axis('off')
axis('equal')

threshold = [0.01, 0.05, 0.1]
for i, thres in enumerate(threshold):
    # enumerate()函数用于将一个可遍历的数据对象(如列表、元组或字符串)组合为一个索引序列，同时列出数据和数据下标，一般用在for 循环当中
        filtered_coords = harris.get_harris_points(harrisim, 6, thres)
                                # get_harris_points(harrisim, min_dist=10, threshold=0.1):
                                #       从一幅Harris 响应图像harrisim中返回角点。min_dist 为分割角点和图像边界的最少像素数目（默认为10）。
                                #       threshold 为阀值，大于阀值的harris响应函数值被认为是可能的角点，并在harrism_t矩阵中相应的位置1,其余地方置0
        subplot(1, 4, i+2)
        imshow(im)
        print im.shape
        plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*') # 标出角点
        axis('off')
        
show()

上面代码中重写compute_harris_response()函数，以便将角点响应函数R输出。

2.2 运行结果

（1）输入图片特征为纹理角点丰富：
在这里插入图片描述

在这里插入图片描述
（2）输入图片特征为纹理平坦丰富：

在这里插入图片描述

（3）输入图片特征为纹理边缘丰富：
在这里插入图片描述

在这里插入图片描述

分析：由上面三种不同场景的输入图片所得到的不同运行结果，可得出：

当输入图片特征为纹理角点丰富时，角点响应函数R为大数值正数；
当输入图片特征为纹理平坦丰富时，角点响应函数R为小数值；
当输入图片特征为纹理边缘丰富时，角点响应函数R为大数值负数

2.3 参数变化

2.3.1 其他不变，改变k（此时sigma=3）

(1)k=0.02
在这里插入图片描述
(2)k=0.04

(3)k=0.06

(4)k=0.08

(5)k=0.1

分析：当k取值在[0,0.04]U[0.06-0.1]之间，角点响应函数R中各数值很接近，不便于角点的选取。因此，k取值范围为[0.04,0.06]时，角点响应函数R效果最好。
注意：此处不研究Harris响应图的原因——函数compute_harris_response()的返回值为Wdet / Wtr，与k无关。

2.3.2 其他不变，改变sigma（k=0.04）

（1）sigma=1
在这里插入图片描述

（2）sigma=3
在这里插入图片描述

（3）sigma=4
在这里插入图片描述

（4）sigma=5

在这里插入图片描述

（5）sigma=7

在这里插入图片描述

（6）sigma=9
在这里插入图片描述

分析：当sigma取值为1时，角点响应函数R过低，找不到任何角点；当sigma取值为7及其以上，角点数量过多；且由Harris响应图可看出，当sigma越大时，角点数量越多；当sigma越小时，角点数量越少。但Harris响应图并非是角点越多越好，故sigma的取值范围落在[3,5]时，Harris响应图效果最好。