How can I find local maxima in an image in MATLAB?

Question

I have an image in MATLAB:

y = rgb2gray(imread('some_image_file.jpg'));

and I want to do some processing on it:

pic = some_processing(y);

and find the local maxima of the output. That is, all the points in y that are greater than all of their neighbors.

I can't seem to find a MATLAB function to do that nicely. The best I can come up with is:

[dim_y,dim_x]=size(pic);
enlarged_pic=[zeros(1,dim_x+2);
              zeros(dim_y,1),pic,zeros(dim_y,1);
              zeros(1,dim_x+2)];

% now build a 3D array
% each plane will be the enlarged picture
% moved up,down,left or right,
% to all the diagonals, or not at all

[en_dim_y,en_dim_x]=size(enlarged_pic);

three_d(:,:,1)=enlarged_pic;
three_d(:,:,2)=[enlarged_pic(2:end,:);zeros(1,en_dim_x)];
three_d(:,:,3)=[zeros(1,en_dim_x);enlarged_pic(1:end-1,:)];
three_d(:,:,4)=[zeros(en_dim_y,1),enlarged_pic(:,1:end-1)];
three_d(:,:,5)=[enlarged_pic(:,2:end),zeros(en_dim_y,1)];
three_d(:,:,6)=[pic,zeros(dim_y,2);zeros(2,en_dim_x)];
three_d(:,:,7)=[zeros(2,en_dim_x);pic,zeros(dim_y,2)];
three_d(:,:,8)=[zeros(dim_y,2),pic;zeros(2,en_dim_x)];
three_d(:,:,9)=[zeros(2,en_dim_x);zeros(dim_y,2),pic];

And then see if the maximum along the 3rd dimension appears in the 1st layer (that is: three_d(:,:,1)):

(max_val, max_i) = max(three_d, 3);
result = find(max_i == 1);

Is there any more elegant way to do this? This seems like a bit of a kludge.

Solution 1

bw = pic > imdilate(pic, [1 1 1; 1 0 1; 1 1 1]);

Solution 2

If you have the Image Processing Toolbox, you could use the IMREGIONALMAX function:

BW = imregionalmax(y);

The variable BW will be a logical matrix the same size as y with ones indicating the local maxima and zeroes otherwise.

NOTE: As you point out, IMREGIONALMAX will find maxima that are greater than or equal to their neighbors. If you want to exclude neighboring maxima with the same value (i.e. find maxima that are single pixels), you could use the BWCONNCOMP function. The following should remove points in BW that have any neighbors, leaving only single pixels:

CC = bwconncomp(BW);
for i = 1:CC.NumObjects,
  index = CC.PixelIdxList{i};
  if (numel(index) > 1),
    BW(index) = false;
  end
end

Solution 3

Alternatively, you can use nlfilter and supply your own function to be applied to each neighborhood.

This "find strict max" function would simply check if the center of the neighborhood is strictly greater than all the other elements in that neighborhood, which is always 3x3 for this purpose. Therefore:

I = imread('tire.tif');
BW = nlfilter(I, [3 3], @(x) all(x(5) > x([1:4 6:9])) );
imshow(BW)

Solution 4

In addition to imdilate, which is in the Image Processing Toolbox, you can also use ordfilt2.

ordfilt2 sorts values in local neighborhoods and picks the n-th value. (The MathWorks example demonstrates how to implemented a max filter.) You can also implement a 3x3 peak finder with ordfilt2 with the following logic:

1. Define a 3x3 domain that does not include the center pixel (8 pixels).

   >> mask = ones(3); mask(5) = 0 % 3x3 max
   mask =
        1     1     1
        1     0     1
        1     1     1
   

2. Select the largest (8th) value with ordfilt2.

   >> B = ordfilt2(A,8,mask)
   B =
        3     3     3     3     3     4     4     4
        3     5     5     5     4     4     4     4
        3     5     3     5     4     4     4     4
        3     5     5     5     4     6     6     6
        3     3     3     3     4     6     4     6
        1     1     1     1     4     6     6     6
   

3. Compare this output to the center value of each neighborhood (just A):

   >> peaks = A > B
   peaks =
        0     0     0     0     0     0     0     0
        0     0     0     0     0     0     0     0
        0     0     1     0     0     0     0     0
        0     0     0     0     0     0     0     0
        0     0     0     0     0     0     1     0
        0     0     0     0     0     0     0     0
   

Solution 5

or, just use the excellent: extrema2.m

View more solutions

Comments

[Nathan Fellman]

Thanks! I see that imregionalmax finds maxima that are greater than or equal to their neighbors. Do you know how I can find only those that are greater and not equal to their neighbors?

[gnovice]

@Nathan: So, if you were to find a set of neighboring maxima that are equal, would you want to just pick one of them, or exclude all of them?

[Nathan Fellman]

oh... and I fixed the question to show that I'm working with grayscale.

[Amro]

yep, this one is even faster :)

[gnovice]

+1 I had forgotten how IMDILATE would work with grayscale images ( I usually only use it with logical masks).

[gnovice]

@Nathan: I modified my answer to show you how to remove maxima that have neighbors, but you should also check out Steve Eddins answer using IMDILATE.

[gnovice]

@Nathan: IMDILATE operates on each pixel of the grayscale image. The center of the 3-by-3 matrix is positioned at each pixel, and the pixel value is replaced by the maximum value found at the neighboring pixels where there is a value of 1 in the 3-by-3 matrix. The call to IMDILATE therefore returns a new matrix where each point is replaced by the maximum value of its 8 neighbors (zero padded at the edges as needed), and the points where the original matrix is larger indicates a local maxima.

[Amro]

imdilate goes over each pixel and computes the max of the neighboring pixels centered around it and specified by the mask given (notice the zero in the middle to exclude the pixel itself). Then we compare the resulting image with the original to check whether each pixel is strictly greater than the max of its neighborhood. Make sure to read the documentation page on morphological operations: mathworks.com/access/helpdesk/help/toolbox/images/…

[Nathan Fellman]

Steve's answer is indeed more elegant.

[Steve Eddins]

Thanks for explaining how this works, gnovice and Amro. I guess my suggestion was a bit cryptic. ;-) One clarification ... imdilate pads by 0 at the edges only for the unsigned integer data types. Generally speaking, imdilate pads by the lowest value for the data type of the input image. That means the technique will work at the edges for all data types - signed and unsigned ints, single-precision and double-precision floating point - even if the array pic has negative values.

[shabbychef]

seems like imdilate is in the image processing toolbox. is there a native matlab solution?

[Dynamite]

Is there any similar way to find local minima ??

[Dynamite]

Is there any similar way to find local miinima ?

[gnovice]

@Dynamite: you can invert the image first so the minima become maxima, then use the same approaches as above. For example, if you have an unsigned 8-bit integer image, you can invert it with "255 - y".

[bla]

If speed is of essence, you may want to try this 2d peak finder It is faster than IMREGIONALMAX by a factor of 6... mathworks.com/matlabcentral/fileexchange/…

[sparkonhdfs]

This is the most correct solution here. It is natively in Matlab and takes alot less time to compute than nfilter does.

[chappjc]

@Franzd'Anconia But I answered 5 years late, so here it is at the bottom. :)

[rayryeng]

Great answer. Is it possible to include the original matrix A? It seems to be missing from your processing chain. I can easily reverse engineer it but it would be nice to include what it was for self-containment :). Thanks!