Harris角点提取并精确至亚像素级,采用标准相机标定图片进行测试-使用matlab2021a或者以上版本测试-源码

% Harris角点提取算法并精确至亚像素级 t='z1.jpg'; % 提示输入要处理的图像 tic; %计时开始 Image = imread(t); % 读取图像 Image = rgb2gray(Image); % 转化为灰度图像 HdImage=Image; % 转化后的灰度图像 fx = [-1 0 1;-1 0 1;-1 0 1]; % x方向的Prewitt算子，用于对x方向滤波 Ix = filter2(fx,HdImage); % 对x方向滤波 fy = [-1 -1 -1;0 0 0;1 1 1]; % y方向的Prewitt算子，用于对y方向滤波 Iy = filter2(fy,HdImage); % 对y方向滤波 Ix2 = Ix.^2; % .^2用来求数组的平方 Iy2 = Iy.^2; % .^2用来求数组的平方 Ixy = Ix.*Iy; % 数组相乘 h= fspecial('gaussian',[9 9],2); % 产生9*9的高斯窗口，sigma=2，产生的窗口越大，得到的角点越少（7-11） A = filter2(h,Ix2); % 用产生的高斯窗口处理Ix2得到A B = filter2(h,Iy2); % 用产生的高斯窗口处理Iy2得到B C = filter2(h,Ixy); % 用产生的高斯窗口处理Ixy得到C height = size(HdImage,1); % 计算图像的第一维的元素数，即行数 width = size(HdImage,2); % 计算图像的第二维的元素数，即列数 CRF = zeros(height,width); % 生成一个和图像大小一致的全0的double型数组，用来保存角点响应函数值 CRFmax = 0; % 保存图像中最大的角点响应函数值 %M = [A(i,j) C(i,j);C(i,j) B(i,j)]; %CRF(i,j) = det(M)-0.05*(trace(M))^2; % 计算角点响应函数值，k的取值一般在0.04--0.06 CRF=(A.*B - C.^2) - 0.05*(A + B).^2; %代码的优化把for循环改为向量循环，k=0.05 CRFmax=max(max(CRF)); %找到最大的角点响应函数（用来设置阈值时用） l=ordfilt2(CRF,7^2,ones(7)); %生成在7*7的窗口进行非最大值抑制（排序滤波器） k=(l==CRF)&(CRF>0.001*CRFmax); %设定阈值为0.01*CRFmax,只有是局部最大值并且角点响应函数值大于阈值才是角点 [v, u] = find(k); % 找到角点是的位置，并保存 count = size(u,1); % 用来记录角点的个数 disp('检测到的角点个数为：') disp(count) % 输出角点个数 figure,imshow(HdImage); % 显示灰度图像 hold on; plot(u,v,'r.'); % 在灰度图像上用红色‘.‘标出角点的位置 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 以下程序把角点精确到亚像素级，所用窗口为3*3+5*5 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% u1=ones(count,1); %保存改进后的角点坐标 v1=ones(count,1); %保存改进后的角点坐标 HdImage=double(HdImage); for i=1:count H=zeros(1,2); %保存梯度 m=zeros(1,2); n=0; if((u(i)+3<=width)&(v(i)+3<=height)&(u(i)>3)&(v(i)>3)) for p=(v(i)-1):(v(i)+1) for q=(u(i)-1):(u(i)+1) H=[HdImage(p,q+1)-HdImage(p,q),HdImage(p+1,q)-HdImage(p,q)]; m=m+H*H'*[q,p]; n=n+H*H'; end; end; for p=(v(i)-2):1:(v(i)+2) for q=(u(i)-2):1:(u(i)+2) H=[HdImage(p,q+1)-HdImage(p,q),HdImage(p+1,q)-HdImage(p,q)]; m=m+H*H'*[q,p]; n=n+H*H'; end; end; g=m/n; u1(i,1)=g(1); v1(i,1)=g(2); else u1(i,1)=u(i,1); v1(i,1)=v(i,1); end end; hold on; plot(u1,v1,'g.'); % 在灰度图像上用绿色‘.‘标出角点优化后的位置 disp('角点的坐标为： u--v 改进后 ') disp([u,v,u1,v1]) %显示原始的和改进后的角点坐标 toc %计时结束