watershed issue

[jip]

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 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scikit-image/attachments/20111228/30dcec97/attachment.html>