基于MATLAB图像harris角点提取算法仿真

figure;
ori_im2 = imread('ucsb1.jpg');     % 读取图像
ori_im2=rgb2gray(ori_im2);
 fx = [6 0 -6;8 0 -8;6 0 -6];          % 高斯函数一阶微分，x方向(用于改进的Harris角点提取算法)
%fx = [-2 -1 0 1 2];                 % x方向梯度算子(用于Harris角点提取算法)
Ix = filter2(fx,ori_im2);              % x方向滤波
 fy = [6 8 6;0 0 0;-6 -8 -6];          % 高斯函数一阶微分，y方向(用于改进的Harris角点提取算法)
%fy = [-2;-1;0;1;2];                 % y方向梯度算子(用于Harris角点提取算法)
Iy = filter2(fy,ori_im2);              % 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_im2,1);
width = size(ori_im2,2);
result = zeros(height,width);         % 纪录角点位置，角点处值为1

R = zeros(height,width);
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)-0.06*(trace(M))^2;   
         end;
end;
cnt = 0;
for i = 2:height-1
    for j = 2:width-1
        % 进行非极大抑制，窗口大小3*3
        if  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;
Rsort=zeros(cnt,1);
[posr, posc] = find(result == 1);
for i=1:cnt
    Rsort(i)=R(posr(i),posc(i));
end;
[Rsort,ix]=sort(Rsort,1);
Rsort=flipud(Rsort);
ix=flipud(ix);
ps=200;
posr2_pic2=zeros(ps,1);
posc2_pic2=zeros(ps,1);
for i=1:ps
    posr2_pic2(i)=posr(ix(i));
    posc2_pic2(i)=posc(ix(i));
end;
    
imshow(ori_im2);
hold on;
plot(posc2_pic2,posr2_pic2,'r+');

 D166