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+');