Harris corner detector

 

The Harris corner detector is a popular interest point detector due to its strong invariance to rotation, scale, illumination variation and image noise. The Harris corner detector is based on the local auto-correlation fuction of a signal; where the local auto-correlation fuction measures the local changes of the signal with patches shifted by a small amount in different directions.

 

 

 


function PIP = harris(I)

 

% Harris detector

% The code calculates

% the Harris Feature Points(FP)

%

% When u execute the code, the test image file opened

% and u have to select by the mouse the region where u

% want to find the Harris points,

% then the code will print out and display the feature

% points in the selected region.

% You can select the number of FPs by changing the variables

% max_N & min_N

% A. Ganoun

 

 

%Aj=6;

%cmin=xmin-Aj; cmax=xmax+Aj; rmin=ymin-Aj; rmax=ymax+Aj;

 

 


 

 

cmin=1; cmax=size(c,1); rmin=1; rmax=size(I,2);

min_N=12;max_N=16;

sigma=2; Thrshold=0; r=6; disp=1;

dx = [-1 0 1; -1 0 1; -1 0 1]; % The Mask

dy = dx';

 


Ix = conv2(I(cmin:cmax,rmin:rmax), dx, 'same');   

Iy = conv2(I(cmin:cmax,rmin:rmax), dy, 'same');

g = fspecial('gaussian',max(1,fix(6*sigma)), sigma); %%%%%% Gaussien Filter

 


A = conv2(Ix.^2, g, 'same');  

B = conv2(Iy.^2, g, 'same');

C = conv2(Ix.*Iy, g,'same');

 


 

k = 0.04;

%k = 0.06;

R11 = (A.*B - C.^2) - k*(A + B).^2;


 

 

R11=(1000/max(max(R11)))*R11;

 

 

R=R11;

ma=max(max(R));

sze = 2*r+1;

MX = ordfilt2(R,sze^2,ones(sze));

R11 = (R==MX)&(R>Thrshold);

  

 

 

count=sum(sum(R11(5:size(R11,1)-5,5:size(R11,2)-5)));

loop=0;

while (((count<min_N)|(count>max_N))&(loop<30))

     if count>max_N

            Thrshold=Thrshold*1.5;

        elseif count < min_N

            Thrshold=Thrshold*0.5;

        end

        

        R11 = (R==MX)&(R>Thrshold);

        

       

        count=sum(sum(R11(5:size(R11,1)-5,5:size(R11,2)-5)));

        loop=loop+1;

end

 

 

R=R*0;

R(5:size(R11,1)-5,5:size(R11,2)-5)=R11(5:size(R11,1)-5,5:size(R11,2)-5);

       

[r1,c1] = find(R);

PIP=[r1+cmin,c1+rmin]%% IP     

end