21

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.

Assad Ebrahim
  • 6,234
  • 8
  • 42
  • 68
Nathan Fellman
  • 122,701
  • 101
  • 260
  • 319
  • Related question: [How can I find many local maxima in a noisy image?](http://stackoverflow.com/questions/2706528/finding-many-local-max-in-an-image-using-matlab) – Jonas Heidelberg Dec 29 '11 at 13:35
  • Does this answer your question? [how to find local maxima in image](https://stackoverflow.com/questions/22218037/how-to-find-local-maxima-in-image) – NoDataDumpNoContribution Jun 02 '22 at 07:24

5 Answers5

37
bw = pic > imdilate(pic, [1 1 1; 1 0 1; 1 1 1]);
Steve Eddins
  • 1,311
  • 7
  • 7
  • yep, this one is even faster :) – Amro Dec 06 '09 at 21:45
  • +1 I had forgotten how IMDILATE would work with grayscale images ( I usually only use it with logical masks). – gnovice Dec 06 '09 at 22:34
  • 4
    @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. – gnovice Dec 07 '09 at 05:07
  • 2
    `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: http://www.mathworks.com/access/helpdesk/help/toolbox/images/f18-12508.html – Amro Dec 07 '09 at 05:10
  • 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. – Steve Eddins Dec 07 '09 at 12:47
  • 3
    seems like imdilate is in the image processing toolbox. is there a native matlab solution? – shabbychef Dec 08 '09 at 17:17
  • Is there any similar way to find local miinima ? – Dynamite Feb 08 '14 at 21:43
18

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
gnovice
  • 125,304
  • 15
  • 256
  • 359
  • 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? – Nathan Fellman Dec 06 '09 at 19:16
  • @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? – gnovice Dec 06 '09 at 19:19
  • oh... and I fixed the question to show that I'm working with grayscale. – Nathan Fellman Dec 06 '09 at 19:47
  • @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 Dec 06 '09 at 22:36
  • 1
    Steve's answer is indeed more elegant. – Nathan Fellman Dec 07 '09 at 06:18
  • Is there any similar way to find local minima ?? – Dynamite Feb 08 '14 at 21:38
  • 1
    @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". – gnovice Feb 08 '14 at 21:45
  • If speed is of essence, you may want to try this 2d peak finder It is faster than IMREGIONALMAX by a factor of 6... http://www.mathworks.com/matlabcentral/fileexchange/37388-fast-2d-peak-finder – bla Oct 03 '14 at 23:08
11

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)
Amro
  • 123,847
  • 25
  • 243
  • 454
3

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
chappjc
  • 30,359
  • 6
  • 75
  • 132
  • 3
    This is the most correct solution here. It is natively in Matlab and takes alot less time to compute than nfilter does. – sparkonhdfs Sep 13 '15 at 22:58
  • @Franzd'Anconia But I answered 5 years late, so here it is at the bottom. :) – chappjc Sep 13 '15 at 23:13
  • 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! – rayryeng Jan 02 '18 at 16:41
2

or, just use the excellent: extrema2.m

Alex Cohen
  • 486
  • 4
  • 3