Harris角点检测算法

Harris角点检测算法
软工1303陈伟峰
1.算法介绍
1988年Harris在Moravec算法的基础上提出了Harris算法。

Harris算法是对moravec算法的改进和提高，harris算法使用高斯函数替代moravec算法的二值窗口函数，另外在moravec中只考虑每个45度的方向的灰度变化，二harris利用泰勒展开式，去近似计算每个方向的灰度变化情况。

2.Harris算法
(1)算法思想
Harris角点检测算法思想就是拿一个小窗在图像中移动，通过考察这个小窗口内图像灰度的平均变换值来确定角点。

相应的会有三种情况发生。

（1）如果窗口内区域图像的灰度值恒定，那么所有不同方向的偏移几乎不发生变化，这个区域属于平坦区域；
（2）如果窗口跨越一条边，那么沿着这条边的偏移几乎不发生变化，但是与边垂直的偏移会发生很大的变化，这个区域属于边缘区域；
（3）如果窗口包含一个孤立的点或者角点，那么所有不同方向的偏移会发生很大的变化。

图1 Harris角点检测基本思想
(2)算法推导
假设窗口W发生位置偏移(u,v)；比较偏移前后窗口中每一个像素点的灰度变化值；使用灰度误差平方和来构造一个误差函数E(u,v)，其中的窗口函数是用来滤波的。

平坦区域：
任意方向移动，无灰度变化边缘：
沿着边缘方向移
动，无灰度变化
角点：
沿任意方向移动，
明显灰度变化
其中w(x,y)为窗口函数，I(x+u,y+v)为平移后的灰度值，I(x,y)为平移前的灰度值。

由Taylor展开式可以得到：
我们定义：
H称为自相关矩阵, λmax和λmin是自相关矩阵的特征值。

如图2所示，其中E（u,v）是一个二次型函数，二次型函数的本质就是一个椭圆，椭圆的扁率和尺寸是由H的特征值λmax和λmin决定的，椭圆的方向由H的特征向量决定。

图2 E(u,v)的椭圆形式
图3 椭圆与点线面的关系
根据图3我们可以看到，当λmax和λmin两者都比较大，并且大小相当时对应点为角点，两者都非常小时为平坦区域；一大一小时为边界区域。

图4 两个特征值的大小对图像点进行分类
因为λmax和λmin在编程时比较难求解，因此1988年，哈里斯在其论文《A combined corner and edge detector》里给出了更有效的角点响应函数：
R为正值时，检测到的是角点；R为负时检测到的是边；R很小时检测到的是平坦区域，由此也就有了更便于计算的数学公式。

3.算法代码实现
clear;
filenema = '图1.jpg';
img = imread(filenema); % 读取图像
Info = imfinfo(filenema); %获取图像相关信息
%%灰度转换
if (Info.BitDepth > 8)
img = rgb2gray(img);
end
%M=[Ix2,Ixy;Ixy,Iy2]
ori_im = double(img) / 255; %unit8转化为64为双精度double64
fx = [-2 -1 0 1 2]; % x方向梯度算子
%filter2表示二维滤波器
Ix = filter2(fx, ori_im); % x方向滤波
fy = [-2; -1; 0; 1; 2]; % y方向梯度算子
Iy = filter2(fy, ori_im); % y方向滤波
Ix2 = Ix .^ 2;
Iy2 = Iy .^ 2;
Ixy = Ix .* Iy;
%%创建一个高斯平滑滤波器
%fspecial创建一个滤波器，参数说明，‘gaussian’表示滤波器类型，Gaussian lowpass %filter，[7,7]表示窗口大小，2表示标准误差
h= fspecial('gaussian', [7,7], 2); % 产生7*7的高斯窗函数，sigma=2
Ix2 = filter2(h,Ix2);
Iy2 = filter2(h,Iy2);
Ixy = filter2(h,Ixy);
[height,width] = size(ori_im);
result = zeros(height, width); %%创建一个矩阵记录角点位置，角点处值为1
R = zeros(height, width);
Rmax = 0; % %设置最大的R值，初始为0
k = 0.04; %k为常系数，经验取值范围为0.04~0.06
%%利用M计算对应于每个像素的角点响应函数R=det(M)-k*(trace(M))^2
for i = 1 : height
for j = 1 : width
M = [Ix2(i, j) Ixy(i, j); Ixy(i, j) Iy2(i, j)]; % 自相关矩阵
% 计算R值,det()求一个方阵的行列式;trace()求方阵的迹，即该方阵对角线上元素之和
R(i,j) = det(M) - k * (trace(M)) ^ 2;
if R(i,j) > Rmax
Rmax = R(i, j);
end;
end;
end;
T = 0.1 * Rmax;%固定阈值，当R(i, j) > T时，则被判定为候选角点
%%下面是进行角点检测记录和局部非极大的检验
winr = floor(7 / 2);
for i = (winr + 1) : (height - winr)
for j = (winr + 1) : (width - winr)
subr = R(((i-winr) : (i + winr)),((j-winr) : (j+ winr)));
submax = max(max(subr));
if((R(i,j) > T) & (R(i,j) == submax))
result(i,j) = 1; %result矩阵中为1表示是角点
end
end
end
[posc, posr] = find(result == 1);%记录角点的位置
figure,imshow(ori_im);
hold on;
plot(posr, posc, 'r+');
实验结果：
第一组（尺寸694*502）
处理前处理后
角点个数：372
第二组（尺寸：556*402）
处理前处理后
角点个数：285
4.总结
Harris角点检测算法是对Moravec算法的改进和提高，在两次实验中我也深有体会，相对来说，Moravec算法比较简单，容易编程，但是结果相对来说不是很准确，而Harris算法，利用高斯窗函数和泰勒展开式，增加了角点检测算法的准确性，相对编程的难度也有所增加。

在没有图像尺度变化的条件下, Harris 角点检测子对图像的旋转变化、视角变化以及图像噪声具有比同类检测子更好的稳定性.。

然而Harris 角点检测子对图像的尺度变化非常敏感，小的尺度变化就会造成Harris 角点的复现率快速下降。

这个特点在上面的实验可以看出来。

图像的尺度伸缩将使得由Harris 检测子提取的特征不具有可重复性。