find_extrema.m
上传用户:trade789
上传日期:2018-05-10
资源大小:603k
文件大小:3k
源码类别:

2D图形编程

开发平台:

Matlab

  1. %/////////////////////////////////////////////////////////////////////////////////////////////
  2. %
  3. % find_extrema - finds local maxima within a grayscale image.  Each point is 
  4. %                checked against all of the pixels within a given radius to be a local max/min.
  5. %                The magnitude of pixel values must be above the given threshold to be picked 
  6. %                as a valid maxima or minima. 
  7. %
  8. % Usage:  m = find_extrema(img,thresh,radius,min_separation)
  9. %
  10. % Parameters:   
  11. %            img :      image matrix
  12. %            thresh :   threshold value
  13. %            radius :   pixel radius          
  14. %
  15. % Returns:
  16. %
  17. %            m       - an nX2 matrix of row,column coordinates of selected points
  18. %
  19. % Author: 
  20. % Scott Ettinger
  21. % scott.m.ettinger@intel.com
  22. %
  23. % May 2002
  24. %/////////////////////////////////////////////////////////////////////////////////////////////
  25. function [mx] = find_extrema(img,thresh,radius)
  26. %img = abs(img);
  27. [h,w] = size(img);
  28. % get interior image subtracting radius pixels from border
  29. p = img(radius+1:h-radius, radius+1:w-radius);  
  30. %get pixels above threshold
  31. [yp,xp] = find(p>thresh);     
  32. yp=yp+radius; xp=xp+radius;
  33. pts = yp+(xp-1)*h;
  34. %create offset list for immediate neighborhood
  35. z=ones(3,3);                    
  36. z(2,2)=0;
  37. [y,x]=find(z);
  38. y=y-2; x=x-2;
  39. if size(pts,2)>size(pts,1)
  40.     pts = pts';
  41. end
  42. %test for max within immediate neighborhood
  43. if size(pts,1)>0
  44.     maxima=ones(length(pts),1);
  45.     for i=1:length(x)
  46.         pts2 = pts + y(i) + x(i)*h;
  47.         maxima = maxima & img(pts)>img(pts2);
  48.     end
  49.     
  50.     xp = xp(find(maxima));                          %save maxima 
  51.     yp = yp(find(maxima));
  52.     pts = yp+(xp-1)*h;                              %create new index list of good points
  53. end
  54.     
  55. %create offset list for radius of pixels
  56. [y,x]=meshgrid(-20:20,-20:20);  
  57. z = (x.^2+y.^2)<=radius^2 & (x.^2+y.^2)>(1.5)^2;    %include points within radius without immediate neighborhood
  58. [y,x]=find(z);
  59. x=x-21; y=y-21;
  60. %create offset list for ring of pixels
  61. [y2,x2]=meshgrid(-20:20,-20:20);    
  62. z = (x2.^2+y2.^2)<=(radius)^2 & (x2.^2+y2.^2)>(radius-1)^2;
  63. [y2,x2]=find(z);
  64. x2=x2-21; y2=y2-21;
  65. maxima = ones(length(pts),1);
  66. %test within radius of pixels (done after first test for slight speed increase)
  67. if size(pts,1)>0
  68.     for i = 1:length(x)
  69.         pts2 = pts + y(i) + x(i)*h;
  70.         maxima = maxima & img(pts)>img(pts2);       %test points   
  71.     end
  72.     xp = xp(find(maxima));                          %save maxima from immediate neighborhood
  73.     yp = yp(find(maxima));
  74.     pts = yp+(xp-1)*h;                              %create new index list
  75.     mx = [yp xp];
  76.     
  77. else
  78.     mx = [];
  79. end
  80.