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.h... 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
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.h...
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
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 <gavet@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.h... 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
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. Sorry about this. Maybe this function should be marked as deprecated in favor of local_maxima. 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 <gavet@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.h...
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
Hi, I haven't tested the difference between both functions, but I can say that peak_local_max and _prominent_peaks are used in several places internally. Another point I would like to make: recently, a scipy contributor added a beautiful function to detect peaks. I'm already using this feature in my code. The first PR is here: https://github.com/scipy/scipy/pull/8264 Others are for performance enhancements. Maybe it would be useful to open an issue and kindly ask this contributor if he is willing to bring his experience on our code? Just a thought I had recently, independently of this topic. -- François Boulogne. http://www.sciunto.org GPG: 32D5F22F
On Wed, 11 Apr 2018 16:07:55 +0200, François Boulogne wrote:
Another point I would like to make: recently, a scipy contributor added a beautiful function to detect peaks. I'm already using this feature in my code. The first PR is here: https://github.com/scipy/scipy/pull/8264 Others are for performance enhancements. Maybe it would be useful to open an issue and kindly ask this contributor if he is willing to bring his experience on our code?
Please do that! The plateau handling in our peak finding seems sub-optimal and confusing. Best regards, Stéfan
Please do that! The plateau handling in our peak finding seems sub-optimal and confusing.
Stephan: it's done here :) https://github.com/scikit-image/scikit-image/issues/3016 Best -- François Boulogne. http://www.sciunto.org GPG: 32D5F22F
On Wed, 11 Apr 2018 12:44:45 +1000, Juan Nunez-Iglesias wrote:
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 [15]: feature.peak_local_max(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.]])
That output in column 1 looks highly suspect! This is a great example for a regression test, thanks Yann. Stéfan
participants (4)
-
François Boulogne -
Juan Nunez-Iglesias -
Stefan van der Walt -
Yann GAVET