Code covered by the BSD License  

Highlights from
Find Image Extrema

image thumbnail

Find Image Extrema

by

Tristan Ursell (view profile)

 

Find local maxima, local minima, saddle points, and flat regions in an image.

[x,y,z,c]=imextrema(im1,varargin)
%Tristan Ursell
%Image Extrema Finder
%May 2013
%
% [x,y,z,c]=imextrema(im1);
% [x,y,z,c]=imextrema(im1,hood);
%
%Estimate extrema in an image at pixel resolution.  The input image 'im1'
%is a grayscale image of any class.  The outputs 'x' and 'y' specify the
%pixel positions of the extrema, 'z' specifies the value of the image at
%the extrema positions, and 'c' classifies the extrema, with:
%
% c = -1 --> local minimum
% c =  0 --> saddle point
% c = +1 --> local maximum
% c = +2 --> locally flat / extrema undefined
%
%The optional variable 'hood' specifies the topology of the pixel
%neighborhood used to find extrema, the options are hood = 4 or hood = 8,
%8 is the default.
%
%Some extrema may not be mathematically defined at a single point, i.e. in
%the example below the saddle points come in pixel doublets -- this is a
%property of image discretization, not an issue with the algorithm.
%
%Image extrema can not occur on the image edge.
%
%If noise is an issue you might consider smoothing the input image with a
%Gaussian blur filter or a hmin (matlab: imhmin) transform.
%
%The algorithm will not find ridge lines.
%
%EXAMPLE 1:
%dx=500;
%Xpos=ones(dx,1)*(1:dx);
%Ypos=Xpos';
%Im1=cos(4*pi*Xpos/dx).*cos(6*pi*Ypos/dx)+5*((Xpos-dx/2).^2/dx^2+(Ypos-dx/2).^2/dx^2);
%Im1(Im1<-0.9)=-0.9;
%
%[x,y,z,c]=imextrema(Im1);
%
%figure
%hold on
%imagesc(Im1)
%plot(x(c==1),y(c==1),'w.')
%plot(x(c==-1),y(c==-1),'k.')
%plot(x(c==0),y(c==0),'kx')
%plot(x(c==2),y(c==2),'wx')
%axis equal tight
%title(['white dots = maxima, black dots = minima, black x = saddles, white x = flat'])
%
%
%EXAMPLE 2:
%Im1=imread('rice.png');
%Im1(Im1<100)=0;
%
%[x,y,z,c]=imextrema(Im1);
%
%figure
%hold on
%imagesc(Im1)
%plot(x(c==1),y(c==1),'rx')
%axis equal tight
%title(['red x = maxima'])
%colormap(gray)

function [x,y,z,c]=imextrema(im1,varargin)

%convert class
im1=double(im1);

%set pixel neighborhood topology
if ~isempty(varargin)
    if varargin{1}==8
        hood8=1;
    else
        hood8=0;
    end
else
    hood8=1;
end

%create edge filter
edge1=ones(size(im1));
edge1(1,:)=0;
edge1(end,:)=0;
edge1(:,1)=0;
edge1(:,end)=0;

if hood8
    %create shift images (3 x 3 neighborhood)
    shift_1=circshift(im1,[-1 0]);
    shift_2=circshift(im1,[1 0]);
    shift_3=circshift(im1,[0 1]);
    shift_4=circshift(im1,[0 -1]);
    shift_5=circshift(im1,[-1 -1]);
    shift_6=circshift(im1,[1 1]);
    shift_7=circshift(im1,[-1 1]);
    shift_8=circshift(im1,[1 -1]);
    
    %get inequality conditions
    cond1_gt=im1>shift_1;
    cond1_lt=im1<shift_1;
    cond2_gt=im1>shift_2;
    cond2_lt=im1<shift_2;
    cond3_gt=im1>shift_3;
    cond3_lt=im1<shift_3;
    cond4_gt=im1>shift_4;
    cond4_lt=im1<shift_4;
    cond5_gt=im1>shift_5;
    cond5_lt=im1<shift_5;
    cond6_gt=im1>shift_6;
    cond6_lt=im1<shift_6;
    cond7_gt=im1>shift_7;
    cond7_lt=im1<shift_7;
    cond8_gt=im1>shift_8;
    cond8_lt=im1<shift_8;
    
    cond1_e=im1==shift_1;
    cond2_e=im1==shift_2;
    cond3_e=im1==shift_3;
    cond4_e=im1==shift_4;
    cond5_e=im1==shift_5;
    cond6_e=im1==shift_6;
    cond7_e=im1==shift_7;
    cond8_e=im1==shift_8;
    
    %find peaks
    highs=cond1_gt.*cond2_gt.*cond3_gt.*cond4_gt...
        .*cond5_gt.*cond6_gt.*cond7_gt.*cond8_gt.*edge1;
    
    %find troughs
    lows=cond1_lt.*cond2_lt.*cond3_lt.*cond4_lt...
        .*cond5_lt.*cond6_lt.*cond7_lt.*cond8_lt.*edge1;
    
    %find saddles
    sads=or(...
        or(cond1_gt.*cond2_gt.*cond3_lt.*cond4_lt,...
        cond1_lt.*cond2_lt.*cond3_gt.*cond4_gt),...
        or(cond5_gt.*cond6_gt.*cond7_lt.*cond8_lt,...
        cond5_lt.*cond6_lt.*cond7_gt.*cond8_gt)).*edge1;
    
    %find locally flat
    flats=cond1_e.*cond2_e.*cond3_e.*cond4_e.*...
        cond5_e.*cond6_e.*cond7_e.*cond8_e.*edge1;
else
    %create shift images (3 x 3 neighborhood)
    shift_1=circshift(im1,[-1 0]);
    shift_2=circshift(im1,[1 0]);
    shift_3=circshift(im1,[0 1]);
    shift_4=circshift(im1,[0 -1]);
    
    %get conditions binary
    cond1_gt=im1>shift_1;
    cond1_lt=im1<shift_1;
    cond2_gt=im1>shift_2;
    cond2_lt=im1<shift_2;
    cond3_gt=im1>shift_3;
    cond3_lt=im1<shift_3;
    cond4_gt=im1>shift_4;
    cond4_lt=im1<shift_4;
    
    cond1_e=im1==shift_1;
    cond2_e=im1==shift_2;
    cond3_e=im1==shift_3;
    cond4_e=im1==shift_4;
    
    %find highs
    highs=cond1_gt.*cond2_gt.*cond3_gt.*cond4_gt.*edge1;    
    
    %find lows
    lows=cond1_lt.*cond2_lt.*cond3_lt.*cond4_lt.*edge1;
    
    %find saddles
    sads=or(cond1_gt.*cond2_gt.*cond3_lt.*cond4_lt,...
        cond1_lt.*cond2_lt.*cond3_gt.*cond4_gt).*edge1;

    %find locally flat
    flats=cond1_e.*cond2_e.*cond3_e.*cond4_e.*edge1;
end

%find locations of extrema
[y_h,x_h]=find(highs);
[y_l,x_l]=find(lows);
[y_s,x_s]=find(sads);
[y_f,x_f]=find(flats);

%find values
z_h=im1(highs==1);
z_l=im1(lows==1);
z_s=im1(sads==1);
z_f=im1(flats==1);

%create classification vector
c_h=ones(length(y_h),1);
c_l=-ones(length(y_l),1);
c_s=zeros(length(y_s),1);
c_f=2*ones(length(y_f),1);

%create outputs
x=[x_h;x_s;x_l;x_f];
y=[y_h;y_s;y_l;y_f];
z=[z_h;z_s;z_l;z_f];
c=[c_h;c_s;c_l;c_f];






Contact us