Matlab simulation code for Harris Corner Point Detection

1. % Harris detector
2. % The code calculates
3. % the Harris Feature Points(FP) 
4. % 
5. % When u execute the code, the test image file opened
6. % and u have to select by the mouse the region where u
7. % want to find the Harris points, 
8. % then the code will print out and display the feature
9. % points in the selected region.
10. % You can select the number of FPs by changing the variables 
11. % max_N & min_N
12. 
    
13. load Imag;
14. 
    
15. I =double(frame);
16. 
    
17. %****************************
18. imshow(frame);
19. 
    
20. k = waitforbuttonpress;
21. point1 = get(gca,'CurrentPoint');  %button down detected
22. rectregion = rbbox;  %%%return figure units
23. point2 = get(gca,'CurrentPoint');%%%%button up detected
24. point1 = point1(1,1:2); %%% extract col/row min and maxs
25. point2 = point2(1,1:2);
26. lowerleft = min(point1, point2);
27. upperright = max(point1, point2); 
28. ymin = round(lowerleft(1)); 
29. 
    
30. %%% Length and height of the bounded rectangle region
31. ymax = round(upperright(1));
32. xmin = round(lowerleft(2));
33. xmax = round(upperright(2));
34. %***********************************
35. 
    
36. Aj=6;
37. cmin=xmin-Aj; cmax=xmax+Aj; rmin=ymin-Aj; rmax=ymax+Aj;
38. min_N=12;max_N=16;
39. 
    
40. %%%%%%%%%%%%%%Interest Points%%%%%%%%%%%%%%%%%%
41. sigma=2; Thrshold=20; r=6; disp=1;
42. dx = [-1 0 1; -1 0 1; -1 0 1]; % The Mask 
43.     dy = dx';
44.    
45.     %%%%%% Derivatives of Ix and Iy
46.     Ix = conv2(I(cmin:cmax,rmin:rmax), dx, 'same');   
47.     Iy = conv2(I(cmin:cmax,rmin:rmax), dy, 'same');
48.     
49.      %%%%%% Gaussien Filter
50.     g = fspecial('gaussian',max(1,fix(6*sigma)), sigma);
51.     
52.     %%%%% Convulation of Ix2, Iy2,Ixy using Gaussian window
53.     Ix2 = conv2(Ix.^2, g, 'same');  
54.     Iy2 = conv2(Iy.^2, g, 'same');
55.     Ixy = conv2(Ix.*Iy, g,'same');
56.    
57.     %%%%%%%%%%%%%% Compute the cornerness values
58.     k = 0.04;
59.     R11 = (Ix2.*Iy2 - Ixy.^2) - k*(Ix2 + Iy2).^2;
60.     R11=(1000/max(max(R11)))*R11;
61.     R=R11;
62.     ma=max(max(R));
63.     sze = 2*r+1; 
64.     MX = ordfilt2(R,sze^2,ones(sze));
65.     R11 = (R==MX)&(R>Thrshold); 
66.     count=sum(sum(R11(5:size(R11,1)-5,5:size(R11,2)-5)));
67.     
68.     %%%%% Find local maxima above some threshold to detect interest points 
69.     loop=0;
70.     while (((count<min_N)|(count>max_N))&(loop<30))
71.         if count>max_N
72.             Thrshold=Thrshold*1.5;
73.         elseif count < min_N
74.             Thrshold=Thrshold*0.5;
75.         end
76.         
77.         R11 = (R==MX)&(R>Thrshold); 
78.         count=sum(sum(R11(5:size(R11,1)-5,5:size(R11,2)-5)));
79.         loop=loop+1;
80.     end
81.     
82.     
83. R=R*0;
84.     R(5:size(R11,1)-5,5:size(R11,2)-5)=R11(5:size(R11,1)-5,5:size(R11,2)-5);
85. [r1,c1] = find(R);
86.     PIP=[r1+cmin,c1+rmin]
87.    
88. 
    
89.    %%%%%%%%%%%%%%%%%%%% Display
90.    
91.    Size_PI=size(PIP,1);
92.    for r=1: Size_PI
93.    I(PIP(r,1)-2:PIP(r,1)+2,PIP(r,2)-2)=255;
94.    I(PIP(r,1)-2:PIP(r,1)+2,PIP(r,2)+2)=255;
95.    I(PIP(r,1)-2,PIP(r,2)-2:PIP(r,2)+2)=255;
96.    I(PIP(r,1)+2,PIP(r,2)-2:PIP(r,2)+2)=255;
97.    
98.    end
99.    
100.    imshow(uint8(I))