<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.实验结果展示