声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 4740|回复: 4

[共享资源] Harris角点提取算法和改进的Harris角点提取算法

[复制链接]
发表于 2006-9-24 07:37 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
  1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  2. %   Harris角点提取算法                                          %
  3. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  4. clear;

  5. ori_im = imread('arroyo-r.tiff');     % 读取图像

  6. % fx = [5 0 -5;8 0 -8;5 0 -5];          % 高斯函数一阶微分,x方向(用于改进的Harris角点提取算法)
  7. fx = [-2 -1 0 1 2];                 % x方向梯度算子(用于Harris角点提取算法)
  8. Ix = filter2(fx,ori_im);              % x方向滤波
  9. % fy = [5 8 5;0 0 0;-5 -8 -5];          % 高斯函数一阶微分,y方向(用于改进的Harris角点提取算法)
  10. fy = [-2;-1;0;1;2];                 % y方向梯度算子(用于Harris角点提取算法)
  11. Iy = filter2(fy,ori_im);              % y方向滤波
  12. Ix2 = Ix.^2;
  13. Iy2 = Iy.^2;
  14. Ixy = Ix.*Iy;
  15. clear Ix;
  16. clear Iy;

  17. h= fspecial('gaussian',[7 7],2);      % 产生7*7的高斯窗函数,sigma=2

  18. Ix2 = filter2(h,Ix2);
  19. Iy2 = filter2(h,Iy2);
  20. Ixy = filter2(h,Ixy);

  21. height = size(ori_im,1);
  22. width = size(ori_im,2);
  23. result = zeros(height,width);         % 纪录角点位置,角点处值为1

  24. R = zeros(height,width);

  25. Rmax = 0;                              % 图像中最大的R值
  26. for i = 1:height
  27.     for j = 1:width
  28.         M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];             % auto correlation matrix
  29.         R(i,j) = det(M)-0.06*(trace(M))^2;                     % 计算R
  30.         if R(i,j) > Rmax
  31.             Rmax = R(i,j);
  32.         end;
  33.     end;
  34. end;

  35. cnt = 0;
  36. for i = 2:height-1
  37.     for j = 2:width-1
  38.         % 进行非极大抑制,窗口大小3*3
  39.         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)
  40.             result(i,j) = 1;
  41.             cnt = cnt+1;
  42.         end;
  43.     end;
  44. end;

  45. [posc, posr] = find(result == 1);
  46. cnt                                      % 角点个数
  47. imshow(ori_im);
  48. hold on;
  49. plot(posr,posc,'r+');
复制代码
回复
分享到:

使用道具 举报

 楼主| 发表于 2006-9-24 07:38 | 显示全部楼层
  1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  2. %   改进的Harris角点提取算法                              %
  3. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  4. clear;

  5. ori_im = imread('arroyo-r.tiff');     % 读取图像

  6. % fx = [5 0 -5;8 0 -8;5 0 -5];          % 高斯函数一阶微分,x方向(用于改进的Harris角点提取算法)
  7. fx = [-2 -1 0 1 2];                 % x方向梯度算子(用于Harris角点提取算法)
  8. Ix = filter2(fx,ori_im);              % x方向滤波
  9. % fy = [5 8 5;0 0 0;-5 -8 -5];          % 高斯函数一阶微分,y方向(用于改进的Harris角点提取算法)
  10. fy = [-2;-1;0;1;2];                 % y方向梯度算子(用于Harris角点提取算法)
  11. Iy = filter2(fy,ori_im);              % y方向滤波
  12. Ix2 = Ix.^2;
  13. Iy2 = Iy.^2;
  14. Ixy = Ix.*Iy;
  15. clear Ix;
  16. clear Iy;

  17. h= fspecial('gaussian',[7 7],2);      % 产生7*7的高斯窗函数,sigma=2

  18. Ix2 = filter2(h,Ix2);
  19. Iy2 = filter2(h,Iy2);
  20. Ixy = filter2(h,Ixy);

  21. height = size(ori_im,1);
  22. width = size(ori_im,2);
  23. result = zeros(height,width);         % 纪录角点位置,角点处值为1

  24. R = zeros(height,width);
  25. for i = 1:height
  26.     for j = 1:width
  27.         M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];             % auto correlation matrix
  28.         R(i,j) = det(M)-0.06*(trace(M))^2;   
  29.          end;
  30. end;
  31. cnt = 0;
  32. for i = 2:height-1
  33.     for j = 2:width-1
  34.         % 进行非极大抑制,窗口大小3*3
  35.         if  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)
  36.             result(i,j) = 1;
  37.             cnt = cnt+1;
  38.         end;
  39.     end;
  40. end;
  41. Rsort=zeros(cnt,1);
  42. [posr, posc] = find(result == 1);
  43. for i=1:cnt
  44.     Rsort(i)=R(posr(i),posc(i));
  45. end;
  46. [Rsort,ix]=sort(Rsort,1);
  47. Rsort=flipud(Rsort);
  48. ix=flipud(ix);
  49. ps=100;
  50. posr2=zeros(ps,1);
  51. posc2=zeros(ps,1);
  52. for i=1:ps
  53.     posr2(i)=posr(ix(i));
  54.     posc2(i)=posc(ix(i));
  55. end;
  56.    
  57. imshow(ori_im);
  58. hold on;
  59. plot(posc2,posr2,'r+');
复制代码

评分

1

查看全部评分

发表于 2009-3-16 21:01 | 显示全部楼层
谢谢,大有收获
发表于 2011-3-6 13:08 | 显示全部楼层
发表于 2013-4-6 11:55 | 显示全部楼层
楼主太厉害了
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-15 21:55 , Processed in 0.064524 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表