how to find local maxima in image

Question

The question is about feature detection concept. I'm stuck after I finding the corner of image and I want to know how to finding the feature point within the computed corners.

Suppose I have grayscale image that have data like this

A = [ 1 1 1 1 1 1 1 1;
      1 3 3 3 1 1 4 1;
      1 3 5 3 1 4 4 4;
      1 3 3 3 1 4 4 4; 
      1 1 1 1 1 4 6 4;
      1 1 1 1 1 4 4 4]

if I use

B = imregionalmax(A);

the result would be like this

B = [ 0 0 0 0 0 0 0 0;
      0 1 1 1 0 0 1 0;
      0 1 1 1 0 1 1 1;
      0 1 1 1 0 1 1 1;
      0 0 0 0 0 1 1 1;
      0 0 0 0 0 1 1 1]

The question is how do I pick the highest peak inside max local region (in sample how did I chose 5 from 3 and 6 from 4)?

My idea was using B to detect each region and use imregionalmax() again but I'm not good at coding and I need some advice or other ideas.

Answer 1

There are a couple of other easy ways to implement a 2D peak finder: ordfilt2 or imdilate.

ordfilt2

The most direct method is to use ordfilt2, which 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 by, (1) using a 3x3 domain that does not include the center pixel, (2) selecting the largest (8th) value and (3) comparing to the center value:

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

There are 8 values considered in this mask, so the 8-th value is the max. The filter output:

>> 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

The trick is compare this to A, the center value of each neighborhood:

>> 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

imdilate

Image dilation is usually done on binary images, but grayscale image dilation is simply a max filter (see Definitions section of imdilate docs). The same trick used with ordfilt2 applies here: define a neighborhood that does not include the center neighborhood pixel, apply the filter and compare to the unfiltered image:

B = imdilate(A, mask);
peaks = A > B;

NOTE: These methods only find a single pixel peak. If any neighbors have the same value, it will not be a peak.

Answer 2

The function imregionalmax gives you the 8-connected region containing the maximum and its 8 neighbours (i.e. the 3x3-regions you are seeing). You could then use morphological operations with the same 3x3 structural element to thin out those regions to their centers. E.g.

 B = imregionalmax(A);
 C = imerode(B, ones(3));

or equivalently

 B = imregionalmax(A);
 D =  bwmorph(B, 'erode');

Alternatively you could write your own maximum finding function using block-processing:

 fun = @(block) % your code working on 'block' goes here ...
 B = blockproc(A, ones(3), fun)

But most likely this will be slower than the built-in functions. (I don't have the toolbox available right now, so I can't try that out.)

Also have a look here and here.