最佳化的Harris角點檢測演算法的MATLAB實現

%%%Prewitt Operator Corner Detection.m
%%%時間最佳化--相鄰像素用取差的方法
%%
clear;
for nfigure=1:6
   
t=input('input your figure’s name(including its extern name):','s');

% t1 = tic;                        %測算時間
FileInfo = imfinfo(t);             % 儲存映像的所有資訊
Image = imread(t);                 % 讀取映像
% 轉為灰階值映像(Intensity Image)
if(strcmp('truecolor',FileInfo.ColorType) == 1) %轉為灰階值映像
Image = im2uint8(rgb2gray(Image));  
end   

dx = [-1 0 1;-1 0 1;-1 0 1];  %dx：橫向Prewitt差分模版
Ix2 = filter2(dx,Image).^2;  
Iy2 = filter2(dx',Image).^2;                                        
Ixy = filter2(dx,Image).*filter2(dx',Image);

%產生 9*9高斯視窗。視窗越大，探測到的角點越少。
h= fspecial('gaussian',9,2);    
A = filter2(h,Ix2);       % 用高斯視窗差分Ix2得到A
B = filter2(h,Iy2);                                
C = filter2(h,Ixy);                                 
nrow = size(Image,1);                           
ncol = size(Image,2);                            
Corner = zeros(nrow,ncol); %矩陣Corner用來儲存候選角點位置,初值全零，值為1的點是角點
                           %真正的角點在137和138行由(row_ave,column_ave)得到
%參數t:點(i,j)八鄰域的“相似性”參數，只有中心點與鄰域其他八個點的像素值之差在
%（-t,+t）之間，才確認它們為相似點，相似點不在候選角點之列
t=20;
%我並沒有全部檢測映像每個點，而是除去了邊界上boundary個像素，
%因為我們感興趣的角點並不出現在邊界上
boundary=8;
for i=boundary:nrow-boundary+1
    for j=boundary:ncol-boundary+1
        nlike=0; %相似點個數
        if Image(i-1,j-1)>Image(i,j)-t && Image(i-1,j-1)<Image(i,j)+t
            nlike=nlike+1;
        end
        if Image(i-1,j)>Image(i,j)-t && Image(i-1,j)<Image(i,j)+t 
            nlike=nlike+1;
        end
        if Image(i-1,j+1)>Image(i,j)-t && Image(i-1,j+1)<Image(i,j)+t 
            nlike=nlike+1;
        end 
        if Image(i,j-1)>Image(i,j)-t && Image(i,j-1)<Image(i,j)+t 
            nlike=nlike+1;
        end
        if Image(i,j+1)>Image(i,j)-t && Image(i,j+1)<Image(i,j)+t 
            nlike=nlike+1;
        end
        if Image(i+1,j-1)>Image(i,j)-t && Image(i+1,j-1)<Image(i,j)+t 
            nlike=nlike+1;
        end
        if Image(i+1,j)>Image(i,j)-t && Image(i+1,j)<Image(i,j)+t 
            nlike=nlike+1;
        end
        if Image(i+1,j+1)>Image(i,j)-t && Image(i+1,j+1)<Image(i,j)+t 
            nlike=nlike+1;
        end
        if nlike>=2 && nlike<=6
            Corner(i,j)=1;%如果周圍有0，1，7，8個相似與中心的（i,j）
                          %那(i,j)就不是角點，所以，直接忽略
        end;
    end;
end;

CRF = zeros(nrow,ncol);    % CRF用來儲存角點響應函數值,初值全零
CRFmax = 0;                % 映像中角點響應函數的最大值，作閾值之用
t=0.05;  
% 計算CRF
%工程上常用CRF(i,j) =det(M)/trace(M)計算CRF，那麼此時應該將下面第105行的
%比例係數t設定大一些，t=0.1對採集的這幾幅映像來說是一個比較合理的經驗值
for i = boundary:nrow-boundary+1
for j = boundary:ncol-boundary+1
    if Corner(i,j)==1  %只關注候選點
        M = [A(i,j) C(i,j);
             C(i,j) B(i,j)];     
         CRF(i,j) = det(M)-t*(trace(M))^2;
        if CRF(i,j) > CRFmax
            CRFmax = CRF(i,j);   
        end;           
    end
end;            
end; 
%CRFmax
count = 0;       % 用來記錄角點的個數
t=0.01;        
% 下面通過一個3*3的視窗來判斷當前位置是否為角點
for i = boundary:nrow-boundary+1
for j = boundary:ncol-boundary+1
        if Corner(i,j)==1  %只關注候選點的八鄰域
            if CRF(i,j) > t*CRFmax && CRF(i,j) >CRF(i-1,j-1) ......
               && CRF(i,j) > CRF(i-1,j) && CRF(i,j) > CRF(i-1,j+1) ......
               && CRF(i,j) > CRF(i,j-1) && CRF(i,j) > CRF(i,j+1) ......
               && CRF(i,j) > CRF(i+1,j-1) && CRF(i,j) > CRF(i+1,j)......
               && CRF(i,j) > CRF(i+1,j+1)
            count=count+1;%這個是角點，count加1
            else % 如果當前位置（i,j）不是角點，則在Corner(i,j)中刪除對該候選角點的記錄
                Corner(i,j) = 0;    
            end;
        end;
end;
end;
% disp('角點個數');
% disp(count)
figure,imshow(Image);      % display Intensity Image
hold on;
% toc(t1)
for i=boundary:nrow-boundary+1
for j=boundary:ncol-boundary+1
        column_ave=0;
        row_ave=0;
        k=0;
       if Corner(i,j)==1
            for x=i-3:i+3  %7*7鄰域
                for y=j-3:j+3
                    if Corner(x,y)==1
% 用算數平均數作為角點座標，如果改用幾何平均數求點的平均座標，對角點的提取意義不大
                        row_ave=row_ave+x;
                        column_ave=column_ave+y;
                        k=k+1;
                    end
                end
            end
        end
        if k>0 %周圍不止一個角點          
             plot( column_ave/k,row_ave/k ,'g.');
        end
end;
end;

end