MatLab角点检测(harris经典程序)

%MatLab角点检测程序harris。 ori_im2=rgb2gray(imread('2.bmp'));     %ori_im2=imresize(ori_im2',0.50,'bicubic');  %加上这句图就变成竖着的了  fx = [5 0 -5;8 0 -8;5 0 -5];          % % la gaucienne，ver axe xIx = filter2(fx,ori_im2);              % la convolution vers axe xfy = [5 8 5;0 0 0;-5 -8 -5];          % la gaucienne，ver axe yIy = filter2(fy,ori_im2);              % la convolution vers axe yIx2 = Ix.^2;Iy2 = Iy.^2;Ixy = Ix.*Iy;clear Ix;clear Iy;h= fspecial('gaussian',[3 3],2);      % générer une fonction gaussienne，sigma=2Ix2 = filter2(h,Ix2);Iy2 = filter2(h,Iy2);Ixy = filter2(h,Ixy);height = size(ori_im2,1);width = size(ori_im2,2);result = zeros(height,width);         % enregistrer la position du coinR = zeros(height,width);K=0.04;Rmax = 0;                              % chercher la valeur maximale de Rfor 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;                     % % calcule R        if R(i,j) > Rmax           Rmax = R(i,j);        end;    end;end;cnt = 0;for i = 2:height-1    for j = 2:width-1        % réduire des valuers minimales ，la taille de fenetre 3*3        if R(i,j) > 0.01*Rmax && 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;[posr2, posc2] = find(result == 1);cnt                                      % compter des coinsfigureimshow(ori_im2);hold on;plot(posc2,posr2,'w*');

 

harris优化的角点检测

 

%%%Prewitt Operator Corner Detection.m %%%时间优化--相邻像素用取差的方法 %% clear; Image = imread('15.bmp');                 % 读取图像 Image = im2uint8(rgb2gray(Image));     dx = [-1 0 1;-1 0 1;-1 0 1];  %dx：横向Prewitt差分模版 Ix2 = filter2(dx,Image).^2;   Iy2 = filter2(dx',Image).^2;                                         Ixy = filter2(dx,Image).*filter2(dx',Image); %生成 9*9高斯窗口。窗口越大，探测到的角点越少。 h= fspecial('gaussian',9,2);     A = filter2(h,Ix2);       % 用高斯窗口差分Ix2得到A B = filter2(h,Iy2);                                 C = filter2(h,Ixy);                                  nrow = size(Image,1);                            ncol = size(Image,2);                             Corner = zeros(nrow,ncol); %矩阵Corner用来保存候选角点位置,初值全零，值为1的点是角点                            %真正的角点在137和138行由(row_ave,column_ave)得到 %参数t:点(i,j)八邻域的“相似度”参数，只有中心点与邻域其他八个点的像素值之差在 %（-t,+t）之间，才确认它们为相似点，相似点不在候选角点之列 t=20; %我并没有全部检测图像每个点，而是除去了边界上boundary个像素， %因为我们感兴趣的角点并不出现在边界上 boundary=8; for i=boundary:nrow-boundary+1     for j=boundary:ncol-boundary+1         nlike=0; %相似点个数         if Image(i-1,j-1)>Image(i,j)-t && Image(i-1,j-1)<Image(i,j)+t             nlike=nlike+1;         end         if Image(i-1,j)>Image(i,j)-t && Image(i-1,j)<Image(i,j)+t              nlike=nlike+1;         end         if Image(i-1,j+1)>Image(i,j)-t && Image(i-1,j+1)<Image(i,j)+t              nlike=nlike+1;         end          if Image(i,j-1)>Image(i,j)-t && Image(i,j-1)<Image(i,j)+t              nlike=nlike+1;         end         if Image(i,j+1)>Image(i,j)-t && Image(i,j+1)<Image(i,j)+t              nlike=nlike+1;         end         if Image(i+1,j-1)>Image(i,j)-t && Image(i+1,j-1)<Image(i,j)+t              nlike=nlike+1;         end         if Image(i+1,j)>Image(i,j)-t && Image(i+1,j)<Image(i,j)+t              nlike=nlike+1;         end         if Image(i+1,j+1)>Image(i,j)-t && Image(i+1,j+1)<Image(i,j)+t              nlike=nlike+1;         end         if nlike>=2 && nlike<=6             Corner(i,j)=1;%如果周围有0，1，7，8个相似与中心的（i,j）                           %那(i,j)就不是角点，所以，直接忽略         end;     end; end;CRF = zeros(nrow,ncol);    % CRF用来保存角点响应函数值,初值全零 CRFmax = 0;                % 图像中角点响应函数的最大值，作阈值之用 t=0.05;   % 计算CRF %工程上常用CRF(i,j) =det(M)/trace(M)计算CRF，那么此时应该将下面第105行的 %比例系数t设置大一些，t=0.1对采集的这几幅图像来说是一个比较合理的经验值 for i = boundary:nrow-boundary+1 for j = boundary:ncol-boundary+1     if Corner(i,j)==1  %只关注候选点         M = [A(i,j) C(i,j);              C(i,j) B(i,j)];               CRF(i,j) = det(M)-t*(trace(M))^2;         if CRF(i,j) > CRFmax             CRFmax = CRF(i,j);            end;                end end;             end;  %CRFmax count = 0;       % 用来记录角点的个数 t=0.01;         % 下面通过一个3*3的窗口来判断当前位置是否为角点 for i = boundary:nrow-boundary+1 for j = boundary:ncol-boundary+1         if Corner(i,j)==1  %只关注候选点的八邻域             if CRF(i,j) > t*CRFmax && CRF(i,j) >CRF(i-1,j-1) ......                && CRF(i,j) > CRF(i-1,j) && CRF(i,j) > CRF(i-1,j+1) ......                && CRF(i,j) > CRF(i,j-1) && CRF(i,j) > CRF(i,j+1) ......                && CRF(i,j) > CRF(i+1,j-1) && CRF(i,j) > CRF(i+1,j)......                && CRF(i,j) > CRF(i+1,j+1)             count=count+1;%这个是角点，count加1             else % 如果当前位置（i,j）不是角点，则在Corner(i,j)中删除对该候选角点的记录                 Corner(i,j) = 0;                 end;         end; end; end; % disp('角点个数'); % disp(count) figure,imshow(Image);      % display Intensity Image hold on; % toc(t1) for i=boundary:nrow-boundary+1 for j=boundary:ncol-boundary+1         column_ave=0;         row_ave=0;         k=0;        if Corner(i,j)==1             for x=i-3:i+3  %7*7邻域                 for y=j-3:j+3                     if Corner(x,y)==1 % 用算数平均数作为角点坐标，如果改用几何平均数求点的平均坐标，对角点的提取意义不大                         row_ave=row_ave+x;                         column_ave=column_ave+y;                         k=k+1;                     end                 end             end         end         if k>0 %周围不止一个角点                        plot( column_ave/k,row_ave/k ,'g.');         end end; end; %end