Harris角点特征提取MATLAB代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Harris角点提取算法                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;

%filename ='Lena.jpg';
filename='object0024.view01.png';
X =imread(filename);    % 读取图像
% imshow(X);
Info = imfinfo(filename); %获取图像相关信息
if (Info.BitDepth > 8)
    f =rgb2gray(X);
end

%《基于特征点的图像配准与拼接技术研究》
%计算图像亮度f(x,y)在点(x,y)处的梯度-----------------------------------------------
% fx = [5 0 -5;8 0 -8;5 0-5];         %高斯函数一阶微分，x方向(用于改进的Harris角点提取算法)
ori_im = double(f) /255;                  %unit8转化为64为双精度double64
fx = [-2-1 0 12];                    % x方向梯度算子(用于Harris角点提取算法)
Ix =filter2(fx,ori_im);               % x方向滤波

% fy = [5 8 5;0 0 0;-5-8-5];         %高斯函数一阶微分，y方向(用于改进的Harris角点提取算法)
fy = [-2; -1; 0; 1;2];                    % y方向梯度算子(用于Harris角点提取算法)
Iy =filter2(fy,ori_im);               % y方向滤波

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

clear Ix;
clear Iy;

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

Ix2 =filter2(h,Ix2);
Iy2 = filter2(h,Iy2);
Ixy = filter2(h,Ixy);

%提取特征点---------------------------------------------------------------
height = size(ori_im, 1);
width = size(ori_im, 2);
result = zeros(height,width);          % 纪录角点位置，角点处值为1

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

for i = 1 : height
    for j = 1 :width
       M = [Ix2(i, j) Ixy(i, j); Ixy(i, j) Iy2(i,j)];            % auto correlation matrix
       R(i,j) = det(M) - k * (trace(M)) ^2;                    %计算R
       if R(i,j) > Rmax
           Rmax = R(i, j);
       end;
    end;
end;

%T= 0.01 * Rmax;%固定阈值，当R(i, j) >T时，则被判定为候选角点
T = 0.1 * Rmax;%固定阈值，当R(i, j)> T时，则被判定为候选角点

%在计算完各点的值后，进行局部非极大值抑制-------------------------------------
cnt = 0;
for i = 2 : height-1
    for j = 2 :width-1
       %进行非极大抑制，窗口大小3*3
       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))
           result(i, j) = 1;
           cnt = cnt+1;
       end;
    end;
end;

i = 1;
    for j = 1 :height
       for k = 1 : width
           if result(j, k) == 1;
               corners1(i, 1) = j;
               corners1(i, 2) = k;
               i = i + 1;
           end;
       end;
    end;

[posc, posr] = find(result == 1);

figure,imshow(ori_im);
hold on;
plot(posr, posc, 'r+');

实验结果如图所示。