Dear all, I try to use watershed implemenation in skimage, compared with scipy as follow: # -*- coding: utf-8 -*- """ Created on Thu Dec 22 13:45:46 2011 @author: Jean-Patrick Pommier Telomeres quantification at low resolution with skimage """ import scipy.ndimage as nd import pylab as plb import numpy as np import skimage as sk import skimage.morphology as mo import skimage.io as io import skimage.filter.edges as edges #import mahotas as mh #import pymorph def gray12_to_8(im): '''Converts a 12 bits image to 8 bits ''' i=0.062271062*im return np.uint8(i) def Hipass(im,size=9): '''returns a hipass filtered image with background >=0 with a gaussian filter of size size=9 by default''' blur=blur=nd.gaussian_filter(im,size) hi=im-1.0*blur-5 mask=np.copy(hi) hi[mask<0]=0 return hi def gray_dist(im): '''Build a pseudo distance map from a high pass filtered image suitable for watershed''' f=Hipass(im,9) dist=f.max()-f dist -= dist.min() dist = dist/float(dist.ptp()) * 255 dist = dist.astype(np.uint8) return dist #============================================================================== # images #============================================================================== path="/home/claire/Applications/ImagesTest/JPP50/16/" CC="DAPI/" PR="CY3/" im="1.TIF" #sk.io.use_plugin('gdal', 'imread') sk.io.use_plugin('freeimage', 'imread') sk.io.use_plugin('matplotlib', 'imshow') dapi=io.imread(path+CC+im) cy3=io.imread(path+PR+im) #============================================================================= # Extracting known nuclei #============================================================================== nuc=dapi[215:264,1151:1198] telo=cy3[215:264,1151:1198] telo=np.uint16(telo) #============================================================================== # Telomeres Segmentation with skimage & ndimage print telo.dtype se=mo.disk(7) top=mo.greyscale_white_top_hat(gray12_to_8(telo),se) #print top.min(),top.max(),top.dtype hi=np.uint8(Hipass(top,size=9)) #print hi.min(), hi.max(),hi.dtype dist=nd.distance_transform_bf(top>10) bassin_dist=np.uint16(dist.max()-dist) bassin_hi=hi.max()-hi remax=mo.is_local_maximum(bassin_hi,hi>6) markers,n=nd.label(remax) print n, 'spots detected' edgesmap=edges.sobel(hi) seg=sk.morphology.watershed(edgesmap,markers,top>30) segn=nd.watershed_ift(bassin_hi,markers) segn_mark_dist=nd.watershed_ift(bassin_dist,markers) #============================================================================== # Quantification #============================================================================== print top.shape,top.dtype #print np.shape(seg),np.dtype(seg) #spots_I=nd.measurements.sum(top,seg) #print spots_I #============================================================================== # #============================================================================== plb.figure(1) plb.subplot(341,frameon=False, xticks=[], yticks=[]) plb.title('raw 12bits') plb.imshow(telo) plb.subplot(342,frameon=False, xticks=[], yticks=[]) plb.title('top Hat') plb.imshow(top) plb.subplot(343,frameon=False, xticks=[], yticks=[]) plb.title('hiP') plb.imshow(hi) plb.subplot(344,frameon=False, xticks=[], yticks=[]) plb.title('regional max') plb.imshow(remax,interpolation=None) plb.subplot(345,frameon=False, xticks=[], yticks=[]) plb.title('markers') plb.imshow(markers,interpolation=None) plb.subplot(346,frameon=False, xticks=[], yticks=[]) plb.title('bassin=hi') plb.imshow(bassin_hi,interpolation=None) plb.subplot(347,frameon=False, xticks=[], yticks=[]) plb.title('bassin dist') plb.imshow(bassin_dist,interpolation=None) plb.subplot(348,frameon=False, xticks=[], yticks=[]) plb.title('hi+mark (nd)') plb.imshow(segn,interpolation=None) plb.subplot(349,frameon=False, xticks=[], yticks=[]) plb.title('sobel (sk)') plb.imshow(edgesmap,interpolation=None) plb.subplot(3,4,10,frameon=False, xticks=[], yticks=[]) plb.title('sob+mark (sk)') plb.imshow(seg,interpolation=None) plb.subplot(3,4,12,frameon=False, xticks=[], yticks=[]) plb.title('dist+mark (nd)') plb.imshow(segn_mark_dist,interpolation=None) plb.show() <https://lh5.googleusercontent.com/-OBPfUOGJNok/TvtKO_6q20I/AAAAAAAAAqw/6VX7eIrySXA/s1600/watershedissue.png> I get a strange result, see the "sob+mark (sk)" image, compared to the segmentations obtained with the scipy implementation. Is it a skimage bug ? Best regards. Jean-Patrick
Hello Jean-Patrick, could you please send a minimal self-contained example (including data) that reproduces the behaviour you observe? Right now it's quite difficult to understand what's going on, because there are a lot of pre-processing steps before the watershed, and we don't have the data to run the script. In the code below you don't use the same elevation function (1st argument of the watershed) for scipy's and skimage's watersheds, do you also get different results when you use the same array for the two watersheds? Cheers, Emmanuelle On Wed, Dec 28, 2011 at 11:05:50AM -0800, jip wrote:
Dear all, I try to use watershed implemenation in skimage, compared with scipy as follow:
# -*- coding: utf-8 -*- """ Created on Thu Dec 22 13:45:46 2011
@author: Jean-Patrick Pommier Telomeres quantification at low resolution with skimage """
import scipy.ndimage as nd import pylab as plb import numpy as np
import skimage as sk import skimage.morphology as mo import skimage.io as io import skimage.filter.edges as edges
#import mahotas as mh #import pymorph
def gray12_to_8(im): ��� '''Converts a 12 bits image to 8 bits ''' ��� i=0.062271062*im ��� return np.uint8(i) ��� def Hipass(im,size=9): ��� '''returns a hipass filtered image with background >=0 with a gaussian filter ��� of size size=9 by default''' ��� blur=blur=nd.gaussian_filter(im,size) ��� hi=im-1.0*blur-5 ��� mask=np.copy(hi) ��� hi[mask<0]=0 ��� return hi ��� def gray_dist(im): ��� '''Build a pseudo distance map from a high pass filtered image suitable ��� for watershed''' ��� f=Hipass(im,9) ��� dist=f.max()-f ��� dist -= dist.min() ��� dist = dist/float(dist.ptp()) * 255 ��� dist = dist.astype(np.uint8) ��� return dist #============================================================================== # images #============================================================================== path="/home/claire/Applications/ImagesTest/JPP50/16/" CC="DAPI/" PR="CY3/" im="1.TIF"
#sk.io.use_plugin('gdal', 'imread') sk.io.use_plugin('freeimage', 'imread') sk.io.use_plugin('matplotlib', 'imshow') dapi=io.imread(path+CC+im) cy3=io.imread(path+PR+im)
#============================================================================= # Extracting known nuclei #============================================================================== nuc=dapi[215:264,1151:1198] telo=cy3[215:264,1151:1198] telo=np.uint16(telo) #============================================================================== # Telomeres Segmentation with skimage & ndimage print telo.dtype
se=mo.disk(7) top=mo.greyscale_white_top_hat(gray12_to_8(telo),se) #print top.min(),top.max(),top.dtype hi=np.uint8(Hipass(top,size=9)) #print hi.min(), hi.max(),hi.dtype dist=nd.distance_transform_bf(top>10) bassin_dist=np.uint16(dist.max()-dist) bassin_hi=hi.max()-hi remax=mo.is_local_maximum(bassin_hi,hi>6) markers,n=nd.label(remax) print n, 'spots detected' edgesmap=edges.sobel(hi) seg=sk.morphology.watershed(edgesmap,markers,top>30) segn=nd.watershed_ift(bassin_hi,markers) segn_mark_dist=nd.watershed_ift(bassin_dist,markers) #============================================================================== # Quantification #============================================================================== print top.shape,top.dtype #print np.shape(seg),np.dtype(seg) #spots_I=nd.measurements.sum(top,seg) #print spots_I #==============================================================================
#============================================================================== plb.figure(1) plb.subplot(341,frameon=False, xticks=[], yticks=[]) plb.title('raw 12bits') plb.imshow(telo) plb.subplot(342,frameon=False, xticks=[], yticks=[]) plb.title('top Hat') plb.imshow(top) plb.subplot(343,frameon=False, xticks=[], yticks=[]) plb.title('hiP') plb.imshow(hi) plb.subplot(344,frameon=False, xticks=[], yticks=[]) plb.title('regional max') plb.imshow(remax,interpolation=None) plb.subplot(345,frameon=False, xticks=[], yticks=[]) plb.title('markers') plb.imshow(markers,interpolation=None) plb.subplot(346,frameon=False, xticks=[], yticks=[]) plb.title('bassin=hi') plb.imshow(bassin_hi,interpolation=None) plb.subplot(347,frameon=False, xticks=[], yticks=[]) plb.title('bassin dist') plb.imshow(bassin_dist,interpolation=None) plb.subplot(348,frameon=False, xticks=[], yticks=[]) plb.title('hi+mark (nd)') plb.imshow(segn,interpolation=None) plb.subplot(349,frameon=False, xticks=[], yticks=[]) plb.title('sobel (sk)') plb.imshow(edgesmap,interpolation=None) plb.subplot(3,4,10,frameon=False, xticks=[], yticks=[]) plb.title('sob+mark (sk)') plb.imshow(seg,interpolation=None) plb.subplot(3,4,12,frameon=False, xticks=[], yticks=[]) plb.title('dist+mark (nd)') plb.imshow(segn_mark_dist,interpolation=None)
plb.show()
[1][IMG]
I get a strange result, see the "sob+mark (sk)" image, compared to the segmentations obtained with the scipy implementation.
Is it a skimage bug ? Best regards.
Jean-Patrick
References
Visible links 1. https://lh5.googleusercontent.com/-OBPfUOGJNok/TvtKO_6q20I/AAAAAAAAAqw/6VX7e...
Dear Emmanuelle, Sorry for the delay. I have rewritten a clearer script <https://docs.google.com/open?id=0B-YDFMbEy1grODU1Y2U0OTgtMzk1MS00MzkxLTliOWItYmMyOTNlOGI3ODRj>I hope. First, it takes a 12bits image<https://docs.google.com/open?id=0B-YDFMbEy1grNWJlYTkwYzYtMjk0ZS00NGVlLWFjMzktMmU3OWU1ODgwMDhm>(cy3 in the archive), isolates a region andremoves the background<https://docs.google.com/open?id=0B-YDFMbEy1grNDZjOWZjZjAtMGVmMS00ZTM2LTgwYzctMjZkM2JiMjFmYjBh> . Regional maxima are found to get markers for watershed. Three kind of maps are also prepared. Each map, with the markers image, is submitted to the watershed algorithm.Segmentation is performed with watershed as indicated here<http://scikits-image.org/docs/dev/auto_examples/plot_watershed.html#example-plot-watershed-py>, as follow: segmentation=watershed(map,markers,mask) map:grayscale image markers:label image mask:binary image The skimage and the scipy implementation are tested with the following result (see image<https://docs.google.com/open?id=0B-YDFMbEy1grZmExNjU5NTAtMWIyZi00ZTAxLWJjZGUtYmNiY2ZkODdmN2Fh> ). Thank you for your help. Best regards Jean-Patrick Pommier http://dip4fish.blogspot.com/
On Fri, Dec 30, 2011 at 5:39 AM, jip <jeanpatrick.pommier@gmail.com> wrote:
Dear Emmanuelle,
Sorry for the delay. I have rewritten a clearer script <https://docs.google.com/open?id=0B-YDFMbEy1grODU1Y2U0OTgtMzk1MS00MzkxLTliOWItYmMyOTNlOGI3ODRj>I hope. First, it takes a 12bits image<https://docs.google.com/open?id=0B-YDFMbEy1grNWJlYTkwYzYtMjk0ZS00NGVlLWFjMzktMmU3OWU1ODgwMDhm>(cy3 in the archive), isolates a region andremoves the background<https://docs.google.com/open?id=0B-YDFMbEy1grNDZjOWZjZjAtMGVmMS00ZTM2LTgwYzctMjZkM2JiMjFmYjBh> . Regional maxima are found to get markers for watershed. Three kind of maps are also prepared. Each map, with the markers image, is submitted to the watershed algorithm.Segmentation is performed with watershed as indicated here<http://scikits-image.org/docs/dev/auto_examples/plot_watershed.html#example-plot-watershed-py>, as follow: segmentation=watershed(map,markers,mask) map:grayscale image markers:label image mask:binary image
There's a small typo here: watershed should be called as: watershed(map, makers, mask=mask) More specifically, in your code, you have: sk.morphology.watershed(bassin_hi,markers,top>30) but it should be called as: sk.morphology.watershed(bassin_hi,markers,mask=top>30) Without specifying `mask=`, you're setting the connectivity parameter in watershed. Even with this correction, something seems to be off since watershed seems to return blank images. Without the masks, it gives what I would expect. I'm not sure what's going on the masked output. -Tony
participants (3)
-
Emmanuelle Gouillart -
jip -
Tony Yu