matlab编程实现Harris角点检测算法

img=imread('Lena.bmp');[m,n]=size(img);tmp=zeros(m+2,n+2);tmp(2:m+1,2:n+1)=img;%扩展图像边缘1个像素Ix=zeros(m+2,n+2);Iy=zeros(m+2,n+2);Ix(:,2:n+1)=tmp(:,3:n+2)-tmp(:,1:n);%x方向上的差分Iy(2:m+1,:)=tmp(3:m+2,:)-tmp(1:m,:);%y方向上的差分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',[77],2);%生成高斯滤波器，消弱噪声的影响Ix2=filter2(h,Ix2);%滤波Iy2=filter2(h,Iy2);Ixy=filter2(h,Ixy);R=zeros(m,n);fori=1:mforj=1:nM=[Ix2(i,j)Ixy(i,j);Ixy(i,j)Iy2(i,j)];%构造Hessian矩阵R(i,j)=det(M)-0.06*(trace(M))^2;%角点判定公式endendRmax=max(max(R));loc=[];%记录角点位置tmp(2:m+1,2:n+1)=R;%扩展图像边缘1个像素fori=2:m+1forj=2:n+1iftmp(i,j)0.01*Rmax%要求矩阵R中的每个元素值大于0.01倍的最大值sq=tmp(i-1:i+1,j-1:j+1);sq=reshape(sq,1,9);sq=[sq(1:4),sq(6:9)];iftmp(i,j)sq%局部非极大值抑制loc=[loc;[j-1,i-1]];endendendend%以下代码显示获取的Harris角点X=loc(:,1);Y=loc(:,2);subplot(1,2,1);imshow(img);title('原图')subplot(1,2,2);imshow(img);title('角点检测')holdonplot(X,Y,'*','color','g');holdoff