<span style="font-size:18px;">function [ptx,pty] = HarrisPoints(ImgIn,threshold)
%   Harris角點提取演算法  
%計算影象亮度f(x,y)在點(x,y)處的梯度-----------------------------------------
fx = [5 0 -5;8 0 -8;5 0 -5];  % 高斯函式一階微分，x方向(用於改進的Harris)                   
%fx = [-2 -1 0 1 2];            % x方向梯度運算元(用於Harris角點提取演算法)  
Ix = filter2(fx, ImgIn);      % x方向濾波  
fy = [5 8 5;0 0 0;-5 -8 -5];  % 高斯函式一階微分，y方向(用於改進的Harris)  
%fy = [-2; -1; 0; 1; 2];        % y方向梯度運算元(用於Harris角點提取演算法)  
Iy = filter2(fy, ImgIn);      % y方向濾波  
%構造自相關矩陣-------------------------------------------------------------
Ix2 = Ix .^ 2;  
Iy2 = Iy .^ 2;  
Ixy = Ix .* Iy;  
clear Ix;  
clear Iy;  
h= fspecial('gaussian', [7 7], 2);% 產生7*7的高斯窗函式，sigma=2  
Ix2 = filter2(h,Ix2);  
Iy2 = filter2(h,Iy2);  
Ixy = filter2(h,Ixy);  
%提取特徵點----------------------------------------------------------------- 
height = size(ImgIn, 1);  
width = size(ImgIn, 2);  
result = zeros(height, width);% 紀錄角點位置，角點處值為1  
R = zeros(height, width);  
Rmax = 0;  % 影象中最大的R值  
k = 0.05; %k為常係數，經驗取值範圍為0.04~0.06  
for i = 1 : height  
    for j = 1 : width  
        M = [Ix2(i, j) Ixy(i, j); Ixy(i, j) Iy2(i, j)];
        R(i,j) = det(M) - k * (trace(M)) ^ 2; % 計算R  
        if R(i,j) > Rmax  
            Rmax = R(i, j);  
        end;  
    end;  
end;  
T = threshold* Rmax;%固定閾值，當R(i, j)>T時，則被判定為候選角點 
%在計算完各點的值後，進行區域性非極大值抑制-------------------------------------  
cnt = 0;  
for i = 2 : height-1  
    for j = 2 : width-1  
        % 進行非極大抑制，視窗大小3*3  
        if (R(i,j)>T && 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;  
[pty, ptx] = find(result == 1);  %row 行；column 列；
end</span>
 

6.實驗結果展示