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
