Harris角点检测

一、harris角点检测的原理和过程
参考课件：http://download.csdn.net/detail/zhangchen1003/9071467

---

二、代码实现

clear;
clc;
I=imread('10.tif');
I=double(I);
%计算差分Ix和Iy
sobelx=[-1 0 1;-2 0 2; -1 0 1];
sobely=[-1 -2 -1;0 0 0;1 2 1];
h=fspecial('gaussian',[3,3],2);
Ix=imfilter(I,sobelx);
Iy=imfilter(I,sobely);
Ix2=Ix.^2;
Iy2=Iy.^2;
Ixy=Ix.*Iy;
%乘以窗口
Ix2=imfilter(Ix2,h);
Iy2=imfilter(Iy2,h);
Ixy=imfilter(Ixy,h);
%
[m,n]=size(I);
R=zeros(m,n);
Rmax=0;
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.04*(trace(M))^2;%0.04代表一个参数
        if R(i,j)>Rmax
            Rmax=R(i,j);
        end
    end
end
%局部非极大值抑制处理
cnt=0;
mark=zeros(m,n);
for i=2:m-1
    for j=2:n-1

        if R(i,j)>0.01*Rmax &&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)
               mark(i,j)=1;
               cnt=cnt+1;
        end

    end
end

[x,y]=find(mark==1);
figure,imshow(uint8(I));
hold on
plot(y,x,'w*');
hold off
cnt