Harris特征点检测算法在matlab中的实现（xk)_xk的博客

http://blog.sina.com.cn/u/2971621807

首页博文目录关于我

个人资料

微博

加好友发纸条

写留言加关注

博客等级：
博客积分：

博客访问：
关注人气：
获赠金笔：0支
赠出金笔：0支
荣誉徽章：

正文字体大小：大中小

Harris特征点检测算法在matlab中的实现（xk)

(2013-02-01 20:05:08)

标签：

harris

算法

matlab

特征点

代码

分类：技术博文

Harris特征点检测算法在matlab中的实现（xk）

一、Harris特征点检测算法是一种经典的特征提取算法。其思想是：

1、指定一图像，设该图像函数为f(x, y)，图像像素点为(x, y)，分别取像素点在x,y方向上的一阶微分导数I_x，I_y，由这两个导数及滤波函数相卷积构造Harris相关矩阵M。

2、将图像中的像素点划分为三种，普通点，边缘点，特征点（也叫角点，以下统称特征点）。

3、矩阵M有两个特征值。并且这这个矩阵为半正定矩阵，所以两个特征值大于等于0。根据这两个特征值的大小可以区分这三种点，找出特征点。

如果两个特征值都比较小，则该点是普通点；

如果两个特征值一个很大，一个很小，则该点为边缘点；

如果两个特征值都较大，则该点是特征点。

4、构造Harris算子，求出一个R值，设定一个阈值T，当R值在取定的局部区域为极值且大于设定的阈值T时，该点为特征点。

说明：Harris算法不受图像的平移旋转、光照变化等的影响，能够稳定的检测特征点。

二、Harris算法的数学表达式及说明如下：

（1）像素点分别在x，y方向上的一阶微分导数：

http://s3/bmiddle/b11f59aftd4a8aa3bca02&690

说明：① w为滤波函数，是均值为零、方差为σ的高斯模板，用于进行高斯平滑，消除孤立噪声点的影响。

② _*为卷积。

（3）构造Harris算子：

http://s2/mw690/b11f59aftd4a8b1580af1&690

说明：① det(M)是矩阵的行列式，trc(M)为迹。

② k为常数，是一个经验值，通常在0.04-0.06之间。

③ k值越小，R值越大，检测到的特征点也就越多

三、matlab代码及注释说明：

clear; %清除以前使用过的函数，防止名字冲突

filename='1.jpg';                          %将要导入的图片定义为filename
X = imread(filename);                 %读取图像到图像X
imshow(X);                                %显示图像X，这一步可有可无
Info = imfinfo(filename);            %获取图像相关信息
if (Info.BitDepth > 8)                  %判断该图像的最大深度是否大于8
    f = rgb2gray(X);                            %将真彩色的图像X转换成灰度图像，并命名为f
end                                            %结束if语句

newim = double(f) / 255;             %unit8转化为64为双精度double64
fx = [-2 -1 0 1 2];                              % x方向梯度算子(用于Harris特征点提取算法)
Ix = filter2(fx, newim);               % x方向滤波，返回newim经过fx算子滤波后的结果Ix

fy = [-2; -1; 0; 1; 2]; % y方向梯度算子(用于Harris特征点提取算法)
Iy = filter2(fy, newim); % y方向滤波

%-----------------------------------------------构造自相关矩阵------------------------------------------------
Ix2 = Ix .^ 2;                              %定义变量Ix2
Iy2 = Iy .^ 2;                              %定义变量Iy2
Ixy = Ix .* Iy;                             %定义变量Ixy

clear Ix; %清除变量Ix
clear Iy; %清除变量Iy

h= fspecial('gaussian', [7 7], 2); %建立滤波算子h，产生7*7的高斯窗函数，sigma=2

Ix2 = filter2(h,Ix2);                            %返回Ix2经过h算子滤波后的结果Ix2
Iy2 = filter2(h,Iy2);                            %返回Iy2经过h算子滤波后的结果Iy2
Ixy = filter2(h,Ixy);                            %返回Ixy经过h算子滤波后的结果Ixy

%--------------------------------------------------提取特征点---------------------------------------------------
height = size(newim, 1);                     %将height定为newim矩阵的行数
width = size(newim, 2);               %将width定位newim矩阵的列数
result = zeros(height, width);        %纪录特征点位置，特征点处值为1

R = zeros(height, width);             %得出R值
Rmax = 0;                                         %让图像中最大的R值为0
k = 0.05;                                   %k为常系数，经验取值范围为0.04~0.06，这里取0.05

for i = 1 : height                                              %在1到height的值之间对以下语句求值
    for j = 1 : width                                           %在1到width的值之间对以下语句求职
        M = [Ix2(i, j) Ixy(i, j); Ixy(i, j) Iy2(i, j)];   %获得矩阵M
        R(i,j) = det(M) - k * (trace(M)) ^ 2;           %计算个点的R值
        if R(i,j) > Rmax                                       %判断个点的R值是否大于0

            Rmax = R(i, j);                                    %如果大于0，则是其成为Rmax
        end;                                                               %结束if
    end;                                                            %结束第二个for
end;                                                                %结束第一个for
T = 0.1 * Rmax;                                               %固定阈值T，当R>T，则定为候选特征点

%----------------------------在计算完各点的值后，进行局部非极大值抑制-----------------------------

for i = 2 : height-1 %在2到height-1的值之间对以下语句求值

for j = 2 : width-1 %在2到width-1的值之间对以下语句求值

if (R(i, j) > T && R(i, j) > R(i-1, j-1) && R(i, j) > R(i-1, j) && R(i, j) > R(i-1, j+1) &&R(i, j) > R(i, j-1) &&R(i, j) > R(i, j+1) && R(i, j) > R(i+1, j-1) && R(i, j) > R(i+1, j) && R(i, j) > R(i+1, j+1)) %进行非极大抑制