%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Harris角點提取算法                                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
filename='22.bmp';
X= imread(filename);     % 讀取圖像
Info=imfinfo(filename);
% if Info.BitDepth>8
%     f=rgb2gray(X);
f=X;
%end
%
% fx = [5 0 -5;8 0 -8;5 0 -5];          % 高斯函數(shù)一階微分，x方向(用于改進的Harris角點提取算法)
ori_im=double(f)/255;                   %unit8轉(zhuǎn)化為64為雙精度double64
fx = [-2 -1 0 1 2];                     % x方向梯度算子(用于Harris角點提取算法)
Ix = filter2(fx,ori_im);                % x方向濾波
% fy = [5 8 5;0 0 0;-5 -8 -5];          % 高斯函數(shù)一階微分，y方向(用于改進的Harris角點提取算法)
fy = [-2;-1;0;1;2];                     % y方向梯度算子(用于Harris角點提取算法)
Iy = filter2(fy,ori_im);                % y方向濾波
Ix2 = Ix.^2;
Iy2 = Iy.^2;
Ixy = Ix.*Iy;
clear Ix;
clear Iy;

h= fspecial('gaussian',7,2);        % 產(chǎn)生7*7的高斯窗函數(shù)，sigma=2

Ix2 = filter2(h,Ix2);
Iy2 = filter2(h,Iy2);
Ixy = filter2(h,Ixy);

height = size(ori_im,1);
width = size(ori_im,2);
result = zeros(height,width);           % 紀(jì)錄角點位置，角點處值為1

R = zeros(height,width);

Rmax = 0;                              % 圖像中最大的R值
for i = 1:height
    for j = 1:width
        M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];             % auto correlation matrix
        R(i,j) = det(M)-0.06*(trace(M))^2;                     % 計算R
        if R(i,j) > Rmax
            Rmax = R(i,j);
        end;
    end;
end;

cnt = 0;
for i = 3:height-2
    for j = 3:width-2
        % 進行非極大抑制，窗口大小3*3
        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)
            result(i,j) = 1;
            cnt = cnt+1;
        end;
    end;
end;

i=1;
    for j=1:height
        for k=1:width
            if result(j,k)==1;
                corners1(i,1)=j;
                corners1(i,2)=k;
                i=i+1;
            end;
        end;
    end;

[posc, posr] = find(result == 1);
cnt                                      % 角點個數(shù)
imshow(ori_im)
hold on;
plot(posr,posc,'r+');