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/6VX7e... 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