欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

数字摄影测量之——Matlab实现Harris算子

程序员文章站 2024-01-19 19:15:16
...

Harris算子理论上来说是Moravec算子的改进,不过实际在求解过程中与Moravec算子的算法有很大不同。首先是使用高斯滤波,对整幅图像进行处处理,以达到降噪的目的。由于Matlab中可以很容易地实现滤波操作,因此在对一阶梯度进行求解时,可以使用[-1 0 1][-1 0 1]T的滤波核对整幅图像进行滤波处理,得到的即是一阶梯度的数据。之后根据得到的一阶梯度,求解图像的自相关矩阵M,通过自相关矩阵M的迹和行列式的值,计算响应因子R,从而根据响应因子R来判定其是否为角点。

具体步骤:

1、对整幅图像做高斯滤波进行降噪

2、使用-1 0 1][-1 0 1]T作为滤波核,分别求出图像的水平和竖直方向的一阶梯度。

3、基于上一步的梯度结果,计算自相关矩阵M,进而计算响应因子R

4、设定阈值,此处选择阈值T=40000

5、进行抑制局部非最大操作,将符合条件的点显示在图像上。


close all;
clear all;
clc;
img=(imread('lena512color.tif'));
img =double( rgb2gray(img));
[X,Y]=size(img);
img2 = imread ('lena512color.tif');
img2 = rgb2gray(img2);
figure;
subplot (1,2,1);
imshow(img2);

imgn=zeros(X,Y);
% 定义高斯滤波
sigma = 1.6;
gausFilter = fspecial('gaussian',5,sigma);
% blur = imfilter(img,gausFilter,'conv');
% 对图像求一阶梯度
fx=[-1,0,1];
Ix=filter2(fx,img);
fy=[-1;0;1];
Iy=filter2(fy,img);
Ix2=Ix.^2;
Iy2=Iy.^2;
Ixy=Ix.*Iy;
gx=filter2(gausFilter,Ix2);
gy=filter2(gausFilter,Iy2);
gxy=filter2(gausFilter,Ixy);
% 遍历所有像元,计算响应因子R 
count1 = 0 ; 
for x=1:X
   for y=1:Y
      %       构建自相关矩阵M
      M=[gx(x,y),gxy(x,y);gxy(x,y),gy(x,y)];
      %      计算响应算子 
      k = 0.04;
      R = det( M )-k*trace(M)^2;    
      imgn(x,y)=R ;
      if(R>0)
      count1 = count1 + 1;
      end
   end
end
% %设阈值,小于均置零
% % % % % %  多种阈值的显示%%%%%%
sortedarray= sort(imgn(:),'descend');%降序排序
% T = mean(sortarray(1:(count1*0.9)));
% T = sum(sortarray(1:(count1*0.2)))/count1;%求前百分之20的均值
% T=0.01*max(max(imgn)); %最大值的0.01为阈值
T=40000;%手动经验阈值
% T  = mean(imgn(:)); %整幅图像均值
% T = mean(sortedarray(1:count1));
% T = 0;
ind=find(imgn<T);
imgn(ind)=0;

subplot(1,2,2);
imshow(img2);hold on;
m=3;
count = 0;
for x=1+m:X-m           %选局部最大(窗口为11*11)且非零值作为特征点
    for y=1+m:Y-m
        window2=imgn(x-m:x+m,y-m:y+m);
        if max(window2(:))==imgn(x,y) && imgn(x,y)>0
            plot(y,x,'+','color','red');
            count=count+1; 
        end
    end
end



程序结果如下:

数字摄影测量之——Matlab实现Harris算子数字摄影测量之——Matlab实现Harris算子