How to make ndimage.filters.maximum_filter work like MATLAB's imregionalmax function?

Question

After reading this post, and play also with SciKit-image I found a difference in Python compared to MATLAB's function imregionalmax.

I have these lines of code:

from skimage.feature import peak_local_max

manos = np.ones([5,5])
manos[2,2] = 0.
manos[2,4] = 2.

giannis = peak_local_max(manos,min_distance=1, indices=False, exclude_border=False)

giorgos = ndimage.filters.maximum_filter(manos, footprint=np.ones([3,3]))
giorgos = (giorgos == manos)

I would expect a 2D array with only one True value ([2,4]) for variables giannis or giorgos as I get in MATLAB. Instead I take more than one maximum.

Any idea why this works this way and how to make it work like in MATLAB?

Answer 1

Both giannis and giorgos are similar in that they find pixels that are equal or larger than the other pixels in the 3x3 neighborhood. I believe giannis would have some additional thresholding.

Neither of these methods guarantee that the pixels found are actually local maxima. Note where I said "larger or equal" above. Any plateau in your image (a region where all pixels have the same value) that is large enough will be marked by the algorithm, no matter if they are local maxima, local minima or somewhere in between.

For example:

import numpy as np
import matplotlib.pyplot as pp
import scipy.ndimage as ndimage

manos = np.sin(np.arange(100)/10)
manos = np.round(30*manos)/30     # Rounding to create plateaus

giorgos = ndimage.filters.maximum_filter(manos, footprint=np.ones([3]))
giorgos = (giorgos == manos)

pp.plot(manos);
pp.plot(giorgos);
pp.show()

Note how the filter identified three points near the local minimum of the sinusoid. The middle one of these is the actual local minimum, the other two are plateaus that are neither local maxima nor minima.

In contrast, the MATLAB function imregionalmax identifies all plateaus that are surrounded by pixels with a lower value. The algorithm required to do this is very different from the one above. It can be efficiently accomplished using a Union-Find algorithm, or less efficiently using a flood-fill-type algorithm. The main idea is to find a pixel that is not lower than any neighbor, then expand from it to its equal-valued neighbors until the whole plateau has been explored or until you find one of the pixels in the plateau with a higher-valued neighbor.

One implementation available from Python is in DIPlib (note: I'm an author):

import diplib as dip

nikos = dip.Maxima(manos)

pp.plot(manos);
pp.plot(nikos);
pp.show()

Another implementation is in SciKit-Image (Thanks to Juan for pointing this out):

nikos = skimage.morphology.local_maxima(manos)