# Re: [scikit-image] local maxima improvements

```Thank you Juan,

```
all my researches from a well-known web search engine gave me peak_local_max as a solution, in particular, the example attached to the watershed function.

```
```
```
```
About peak_local_max, there is an option indices=False that results in an array with the same shape as the image in argument.
```
yann

Le 11/04/2018 à 04:44, Juan Nunez-Iglesias a écrit :
```
```Hi Yann, and thanks for the interest!

```
We actually already have this algorithm implemented in skimage.morphology.local_maxima.
```
```
peak_local_max is a bit different and I must admit I don’t understand the logic in it. I *particularly* don’t understand the following result:
```
In [1]: def rmax(I):
...:     """
...:     Own version of regional maximum
...:     This avoids plateaus problems of peak_local_max
...:     I: original image, int values
...:     returns: binary array, with 1 for the maxima
...:     """
...:     I = I.astype('float');
...:     I = I / np.max(I) * (2**31 -2);
...:     I = I.astype('int32');
...:     h = 1;
...:     rec = morphology.reconstruction(I, I+h);
...:     maxima = I + h - rec;
...:     return maxima
...:
...:

In [2]: image = np.zeros((6, 6))
In [3]: image[1, 1] = 1
In [4]: image[3:] = 2
In [6]: image[3:-1, 3:-1] = 4

In [7]: image
Out[7]:
array([[ 0.,  0.,  0.,  0.,  0.,  0.],
[ 0.,  1.,  0.,  0.,  0.,  0.],
[ 0.,  0.,  0.,  0.,  0.,  0.],
[ 2.,  2.,  2.,  4.,  4.,  2.],
[ 2.,  2.,  2.,  4.,  4.,  2.],
[ 2.,  2.,  2.,  2.,  2.,  2.]])

In [8]: rmax(image)
Out[8]:
array([[ 0.,  0.,  0.,  0.,  0.,  0.],
[ 0.,  1.,  0.,  0.,  0.,  0.],
[ 0.,  0.,  0.,  0.,  0.,  0.],
[ 0.,  0.,  0.,  1.,  1.,  0.],
[ 0.,  0.,  0.,  1.,  1.,  0.],
[ 0.,  0.,  0.,  0.,  0.,  0.]])

In [12]: morphology.local_maxima(image)
Out[12]:
array([[0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 1, 0],
[0, 0, 0, 1, 1, 0],
[0, 0, 0, 0, 0, 0]], dtype=uint8)

In [15]: feature.peak_local_max(image)
Out[15]:
array([[4, 4],
[4, 3],
[4, 1],
[3, 4],
[3, 3],
[3, 1],
[1, 1]])

In [16]: image_peak = np.zeros_like(image)

In [17]: image_peak[tuple(feature.peak_local_max(image).T)] = 1

In [18]: image_peak
Out[18]:
array([[ 0.,  0.,  0.,  0.,  0.,  0.],
[ 0.,  1.,  0.,  0.,  0.,  0.],
[ 0.,  0.,  0.,  0.,  0.,  0.],
[ 0.,  1.,  0.,  1.,  1.,  0.],
[ 0.,  1.,  0.,  1.,  1.,  0.],
[ 0.,  0.,  0.,  0.,  0.,  0.]])

Anyone else on the list care to comment???

Juan.

On 10 Apr 2018, 11:49 PM +1000, Yann GAVET <ga...@emse.fr>, wrote:
```
First mistake, this should work, but the discretization of the 'continuous' values should be handled with care.
```
def rmax(I):

"""

Own version of regional maximum

This avoids plateaus problems of peak_local_max

I: original image, int values

returns: binary array, with 1 for the maxima

"""

I = I.astype('float');

I = I / np.max(I) * (2**31 -2);

I = I.astype('int32');

h = 1;

rec = morphology.reconstruction(I, I+h);

maxima = I + h - rec;

return maxima

Le 10/04/2018 à 15:35, Yann GAVET a écrit :
```
```
Dear all,

```
I have been playing around with the watershed segmentation by markers with the code proposed as example:
```
http://scikit-image.org/docs/dev/auto_examples/segmentation/plot_watershed.html

```
Unfortunately, if we use for example floating values for the radii of the circles (like r1, r2 = 20.7, 24.7), the separation is not perfect, as it gives 4 labels.
```
```
If we use the chamfer distance transform instead of the Euclidean distance transform, it is even worse.
```
```
It appears that the markers detection by regional maximum (peak_local_max) fails in the presence of plateaus. Its algorithm is basically D(I)==I, where D is the morphological dilation.
```
```
A better algorithm would be to use morphological reconstruction (see SOILLE, Pierre. /Morphological image analysis: principles and applications/. Springer Science & Business Media, 2003, p202, Eq 6.13). A proposition of the code can be the following (it should deal with float values):
```

import numpy as np

from skimage import morphology

def rmax(I):

"""

This avoids plateaus problems of peak_local_max

I: original image, float values

returns: binary array, with True for the maxima

"""

I = I.astype('float');

I = I / np.max(I) * 2**31;

I = I.astype('int32');

rec = morphology.reconstruction(I-1, I);

maxima = I - rec;

return maxima>0

```
This code is relatively fast. Notice that the matlab function imregionalmax seem to work the same way (the help is not explicit, but the results on a few tests seem to be similar).
```
```
I am afraid I do not have time to integrate it on gitlab, but this should be a good start if someone wants to work on it. If you see any problem with this code, please correct it.
```
thank you

best regards

--
Yann

_______________________________________________
scikit-image mailing list
scikit-image@python.org
https://mail.python.org/mailman/listinfo/scikit-image
```
```
--
Yann GAVET
Assistant Professor - Ecole Nationale Supérieure des Mines de Saint-Etienne
158 Cours Fauriel, CS 62362, 42023 SAINT-ETIENNE cedex 2
Tel: (33) - 4 7742 0170
_______________________________________________
scikit-image mailing list
scikit-image@python.org
https://mail.python.org/mailman/listinfo/scikit-image
```
```

_______________________________________________
scikit-image mailing list
scikit-image@python.org
https://mail.python.org/mailman/listinfo/scikit-image
```
```
--
Yann GAVET
Assistant Professor - Ecole Nationale Supérieure des Mines de Saint-Etienne
158 Cours Fauriel, CS 62362, 42023 SAINT-ETIENNE cedex 2
Tel: (33) - 4 7742 0170

```
```_______________________________________________
scikit-image mailing list
scikit-image@python.org
https://mail.python.org/mailman/listinfo/scikit-image
```