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

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.

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

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

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)

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
   

or, just use the excellent: extrema2.m