一、代码
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 = 2h = fspecial('gaussian',[3 3],2);Ix2 = filter2(h,Ix2);Iy2 = filter2(h,Iy2);Ixy = filter2(h,Ixy);% 初始化RmaxRmax = 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);endendend% 显示图像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')endendend% 这是测试plot绘制点的代码% for i = 1 : m%for j = 1 : n% plot(i,j,'b+');% pause;%end% end
