my_harris.m

自已编的一个matlab的harris角点检测的代码

clear
clc
load J
I=double(J);
imshow(J)%读取图像
k = waitforbuttonpress;%鼠标响应
%%%%%%截取图像
point1=get(gca,'CurrentPoint');
rectregion = rbbox;
point2=get(gca,'CurrentPoint');
%%%%%%%找到截得的图像的xmax,xmin,ymin,ymax
point1=point1(1,1:2);
point2=point2(1,1:2);
lowerup=min(point1,point2);
rightlower=max(point1,point2);
xmin=round(lowerup(2));
xmax=round(rightlower(2));
ymin=round(lowerup(1));
ymax=round(rightlower(1));
%%%%%%%%将截得的图像扩大Aj倍
Aj=6;
levelmin=xmin-Aj;
levelmax=xmax+Aj;
highmin=ymin-Aj;
highmax=ymax+Aj;
%%%%%%%%%求Ix,Iy,Ix2,Iy2,IxIy
sigma=2;
h=[-1 0 1;-1 0 1;-1 0 1];
Ix=conv2(I(levelmin:levelmax,highmin:highmax),h, 'same');
Iy=conv2(I(levelmin:levelmax,highmin:highmax),h',  'same');
g=fspecial('gaussian',max(1,fix(sigma*6)),sigma);
Ix2=conv2(Ix.^2,g,'same');
Iy2=conv2(Iy.^2,g,'same');
Ixy=conv2(Ix.*Iy,g,'same');
%%%%%%%%%%%
k=0.04;r=6;sze=2*r+1;Threshold=20;
R=(Ix2.*Iy2-Ixy.^2)-k*(Ix2+Iy2).^2;%%%%%计算角点量
%R=(Ix2.*Iy2-(Ixy).^2)/(Ix2+Iy2);
R=(1000/max(max(R)))*R;
R11=R;
MX=ordfilt2(R,sze^2,ones(sze));
R=(R11==MX)&(R11>Threshold);
cout=sum(sum(R(5:size(R,1)-5,5:size(R,2)-5)));
c_min=50;
c_max=100;
loop=0;
while ((cout<c_min)|(cout>c_max))&(loop<=30)
    if cout<c_min
        Threshold=Threshold*0.5;
    elseif cout>c_max
            Threshold=Threshold*1.5;
        end

    loop=loop+1;
 R=(R11==MX)&(R11>Threshold);
 cout=sum(sum(R(5:size(R,1)-5,5:size(R,2)-5)));
end
R11=R11*0;
R11(5:size(R11,1)-5,5:size(R11,2)-5)=R(5:size(R,1)-5,5:size(R,2)-5);
[iv,ju]=find(R11);
IP=[iv+levelmin,ju+highmin]
%%%%%将角点标出
Size_IP=size(IP,1);
for i=1:Size_IP
   I(IP(i,1)-2:IP(i,1)+2,IP(i,2)-2)=255;
   I(IP(i,1)-2:IP(i,1)+2,IP(i,2)+2)=255;
   I(IP(i,1)-2,IP(i,2)-2:IP(i,2)+2)=255;
   I(IP(i,1)+2,IP(i,2)-2:IP(i,2)+2)=255;
end
imshow(uint8(I))