Harris角点检测及matlab程序

close all;
clear all;
clc;

img=imread('lena.png');
imshow(img);
img = rgb2gray(img);
img =double(img);
[m n]=size(img);

tmp=zeros(m+2,n+2);
tmp(2:m+1,2:n+1)=img;
Ix=zeros(m+2,n+2);
Iy=zeros(m+2,n+2);

E=zeros(m+2,n+2);

Ix(:,2:n)=tmp(:,3:n+1)-tmp(:,1:n-1);
Iy(2:m,:)=tmp(3:m+1,:)-tmp(1:m-1,:);

Ix2=Ix(2:m+1,2:n+1).^2;
Iy2=Iy(2:m+1,2:n+1).^2;
Ixy=Ix(2:m+1,2:n+1).*Iy(2:m+1,2:n+1);

h=fspecial('gaussian',[7 7],2);
Ix2=filter2(h,Ix2);
Iy2=filter2(h,Iy2);
Ixy=filter2(h,Ixy);

Rmax=0;
R=zeros(m,n);
for i=1:m
for j=1:n
M=[Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];
R(i,j)=det(M)-0.06*(trace(M))^2;

if R(i,j)>Rmax
Rmax=R(i,j);
end
end
end
re=zeros(m+2,n+2);

tmp(2:m+1,2:n+1)=R;
img_re=zeros(m+2,n+2);
img_re(2:m+1,2:n+1)=img;
for i=2:m+1
for j=2:n+1

if tmp(i,j)>0.02*Rmax &&...
tmp(i,j)>tmp(i-1,j-1) && tmp(i,j)>tmp(i-1,j) && tmp(i,j)>tmp(i-1,j+1) &&...
tmp(i,j)>tmp(i,j-1) && tmp(i,j)>tmp(i,j+1) &&...
tmp(i,j)>tmp(i+1,j-1) && tmp(i,j)>tmp(i+1,j) && tmp(i,j)>tmp(i+1,j+1)
img_re(i,j)=255;
end
end
end

img_re=mat2gray(img_re(2:m+1,2:n+1));
figure,imshow(img_re);