一、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