一、代码

close all;
clc;
% 读取图像信息（原图为灰度图）
img = imread('lena.bmp');
[m,n] = size(img);
% 先在原图外围扩展一圈
tmp = zeros(m+2,n+2);
tmp(2:m+1,2:n+1) = img;
% 初始化各一阶偏导矩阵
Ix = zeros(m+2,n+2);
Iy = zeros(m+2,n+2);
E = zeros(m+2,n+2);
% 求偏导
Ix(:,2:n) = tmp(:,3:n+1) - tmp(:,1:n-1);
Iy(2:m,:) = tmp(3:m+1,:) - tmp(1:m-1,:);
Ix2 = Ix(2:m+1,2:n+1) .^ 2;
Iy2 = Iy(2:m+1,2:n+1) .^ 2;
Ixy = Ix(2:m+1,2:n+1) .* Iy(2:m+1,2:n+1);
%生成高斯卷积核，便于对Ix2、Iy2、Ixy进行平滑
% sigma = 2
h = fspecial('gaussian',[3 3],2);
Ix2 = filter2(h,Ix2);
Iy2 = filter2(h,Iy2);
Ixy = filter2(h,Ixy);% 初始化Rmax
Rmax = 0;
R = zeros(m,n);
for i = 1 : mfor j = 1 : nM = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];R(i,j) = det(M) - 0.06 * (trace(M))^2;if R(i,j) > RmaxRmax = R(i,j);endend
end
% 显示图像
imshow(img);
title('角点检测');
hold on;% 求角点
tmp(2:m+1,2:n+1) = R;
result = zeros(m+2,n+2);
result(2:m+1,2:n+1) = img;
for i = 2 : m + 1for j = 2 : n + 1% 阈值为0.02*Rmax% 认为R值大于阈值的点为角点% 求当前像素点的邻域current = [tmp(i-1,j-1) tmp(i-1,j) tmp(i-1,j+1);tmp(i,j-1)   tmp(i,j)   tmp(i,j+1);tmp(i+1,j-1) tmp(i+1,j) tmp(i+1,j+1)];% 若当前像素点的R值大于阈值且它是其八邻域中R值最大的点，则它为角点if tmp(i,j) >= 0.02 * Rmax && tmp(i,j) >= max(max(current))result(i,j) = 255;% plot绘制点的时候是以左上角为原点，水平向右为x正半轴轴，竖直向下为y正半轴% 这和我们对于图像矩阵坐标的直观印象恰好相反plot(j,i,'go')endend
end% 这是测试plot绘制点的代码
% for i = 1 : m
%     for j = 1 : n
%         plot(i,j,'b+');
%         pause;
%     end
% end

二、结果