Segmentation of a complex metaphase with pymorph

Part of the original 16bits TIF image

Labelled image displayed with pylab

The mean gray level value is also calculated for each labelled region with pymorph.grain:

[ 1301.88575849    23.00804348   704.34775673   629.5625966   1410.2556391
   871.05686546   946.11368209   841.65226782   860.98429319  1526.48712121
   705.22321429   608.3359841    929.86325301   858.11372549  1000.02900378
   831.46130952   717.89433962   790.23236741   919.48754448   603.34122288
   830.35550936  1325.0239521    843.45842697   880.9575       864.25319693
   970.12876254   901.85921959   914.75644699   902.17860963   864.04659498
   960.92638037   892.95342466   946.98757358  1119.16554163   835.44263629
   872.87977369   838.04045734   734.5129683   1147.49576784   995.61397059
   839.4835443    823.52125984   795.88679245   898.91109422   748.07751938
   635.60344828]

the python code used for that result:

from scipy import ndimage as ND
import scipy
import numpy as np
import readmagick
import mahotas
import pymorph
import pylab
import os
def SegmentMetaphase(image):
    #print "def image:",image.mean()
    #print "def background:",back.mean()
    #print "def noback:",im.mean()
    noBack=RemoveModalBackground(image)
    print "from Segmeta noBack:",noBack.min(),noBack.mean()
    blurLowRes=ND.filters.gaussian_filter(noBack,13)
    blurHiRes=ND.filters.gaussian_filter(noBack,1)
    midPass=pymorph.subm(blurHiRes,0.70*blurLowRes)
    bin=(midPass>1.5*midPass.mean())
    binLowRes=pymorph.open(bin,pymorph.sedisk(4))
    binLR=pymorph.edgeoff(binLowRes)
    blur=ND.filters.gaussian_filter(noBack,9)
    hiPass=pymorph.subm(noBack,1.0*blur)
    gradLowRes=pymorph.gradm(noBack)
    #build seeds for watershed
    bSeeds1=(hiPass>hiPass.mean())
    bSeeds2=scipy.logical_and(bSeeds1,binLR)
    #bSeeds3=ND.binary_fill_holes(bSeeds2)
    bSeeds4=ND.binary_opening(bSeeds2,iterations=4)
    outline=mahotas.bwperim(bin, n=4)
    #display=pymorph.overlay(image,outline)#Hugly
    markers,nr_obj=ND.label(bSeeds4)
    labeled=mahotas.cwatershed(gradLowRes, markers, Bc=None, return_lines=False)
    return labeled
def RemoveModalBackground(image):
    mode=ModalValue(image)
    back = np.zeros(image.shape, image.dtype)
    back.fill(mode)
    #print "def background:",back.mean()
    im=pymorph.subm(image,back)
    return im
def ModalValue(image):
    '''look for the modal value of an image'''
    #print image.dtype
    if image.dtype=="uint8":
        depthmax=255
        print "8bits"
    if image.dtype=="uint16":
        depthmax=65535
        print "16bits"
    histo=mahotas.fullhistogram(image)
    countmax=histo.max()
    print "countmax:",countmax
    print "image max",image.max()
    mig=image.min()#image min graylevel
    mag=image.max()#image max gray level
    mode=0
    countmax=0#occurence of a given grayscale
    print "mig=",mig,"  mag=",mag
    for i in range(mig,mag-1,1):
        test=histo[i]>countmax
        #print "test:",test,"histo(",i,")=", histo[i],"max",countmax
        if  test:
            countmax=histo[i]
            mode=i
            #print "mode",mode
    return mode
def BorderKill(imlabel):
    '''remove labelled objects touching the image border'''
    #from pythonvision.org
    #pymorph.edgeoff should do the same
    whole = mahotas.segmentation.gvoronoi(imlabel)
    borders = np.zeros(imlabel.shape, np.bool)
    borders[ 0,:] = 1
    borders[-1,:] = 1
    borders[:, 0] = 1
    borders[:,-1] = 1
    at_border = np.unique(imlabel[borders])
    for obj in at_border:
        whole[whole == obj] = 0
    return whole
user=os.path.expanduser("~")
workdir=os.path.join(user,"Applications","ImagesTest","JPP60-3","DAPI")
file="1.tif"
complete_path=os.path.join(workdir,file)
if __name__ == "__main__":
    #print complete_path
    dapi=readmagick.readimg(complete_path)
    readmagick.writeimg(dapi,os.path.join(workdir,"dapi.png"))
    NoBackIm=RemoveModalBackground(dapi)
    dapilab1=SegmentMetaphase(dapi)
    readmagick.writeimg(dapilab1,os.path.join(workdir,"dapilabel.png"))
    print pymorph.grain(NoBackIm,dapilab1, 'mean', option='data')
    pylab.imshow(dapilab1)
    pylab.show()

Metaphasic chromosomes segmentation under python

The following python script implements three ways to segment metaphasic chromosomes through functions relying on scipy.ndimage and pymorph/mahotas. Starting from initial image:

 

Human DAPI stained metaphasic chromosomes recorded with KAF1600 Peltier cooled CCD camera.

• The first idea consists in performing a hipass filtering followed by a simple thresholding step and finally cleaning the binary image with a morphological opening. This is achieved with the function def HighPassBinary(image):

Labelled image displayed with pylab. Hipass filtering results in fragmented nuclei

•  The second idea, implemented in def metaphaseSeg(image):, tries to preserve both nuclei and chromosomes shape using a region growing algorithm (watershed + markers). The markers are obtained from a regional maxima on a blured DAPI image. However, according to the size of the gaussian kernel the segmented image can be over or under segmented with undetected chromosomes as shown in the following animation:

gif animation of a region based segmentation with increasing gaussian kernel size (3, 5, 7, 9,11, 13, 19). A small kernel (3) yields an oversegmented image.

• A third idea tries to combine both methods using the result of the first method as markers to initiate a region growing, implemented in def HighPassMetaphaseSeg(image):

Artifacts close to the nuclei yield background artifacts near the nuclei.

The result of this third method wouldn't be too bad if additionnal image cleaning was performed on the first image.

 The python script runs under ubuntu linux , I didn't take care of path to access to the image file to make it portable for windows or MacOS. To use one of the function comment/uncomment the code.

from scipy import ndimage as ND
import readmagick
import mahotas
import pymorph
import pylab
def HighPassBinary(image):
    '''depuis une image DAPI,  genere une image passe haut avec fond a zero'''
    '''necessite scipy.ndimage'''
    blur=ND.filters.gaussian_filter(image, 9)
    sharp=image-1.0*blur
    '''float 64bits'''
    print sharp.dtype, sharp.shape
    smax=sharp.max()
    '''threshold above negative values'''
    imBin=(sharp>0)
    imBin2=ND.binary_fill_holes(imBin)
    imBin3=ND.binary_opening(imBin2, iterations=4)
    labeled,nr_obj=ND.label(imBin3)
    print nr_obj
    return labeled

def metaphaseSeg(image):
    blur=ND.filters.gaussian_filter(image, 15)
    remax=pymorph.regmax(blur)
    binRemax=(remax>0)
    seeds,nr_obj=ND.label(binRemax)
    print nr_obj
    grad=pymorph.gradm(image)
    label=mahotas.cwatershed(grad, seeds, Bc=None, return_lines=False)
    return label
def HighPassMetaphaseSeg(image):
    '''depuis une image DAPI,  genere une image passe haut avec fond a zero'''
    '''necessite scipy.ndimage'''
    blur=ND.filters.gaussian_filter(image, 9)
    sharp=image-1.0*blur
    '''float 64bits'''
    print sharp.dtype, sharp.shape
    smax=sharp.max()
    '''threshold negative values'''
    imBin=(sharp>0)
    imBin2=ND.binary_fill_holes(imBin)
    imBin3=ND.binary_opening(imBin2, iterations=4)
    seeds,nr_obj=ND.label(imBin3)
    print nr_obj
    grad=pymorph.gradm(ND.filters.gaussian_filter(image, 1))
    labeled=mahotas.cwatershed(grad, seeds, Bc=None, return_lines=False)
    return labeled
   
if __name__ == "__main__":
    dapi=readmagick.readimg("/home/jeanpat/Applications/Images_test/png/essai09/humain.png")
    '''dapiLabeled=metaphaseSeg(dapi)'''
    dapiLabeled1=HighPassBinary(dapi)
    '''dapiLabeled2=HighPassMetaphaseSeg(dapi)'''
    pylab.imshow(dapiLabeled1)
    pylab.show()
    '''pylab.imshow(dapiLabeled1)
    pylab.show()
    pylab.imshow(dapiLabeled2)
    pylab.show()'''
    '''readmagick.writeimg(dapiBin,"/home/jeanpat/Applications/Images_test/png/essai09/humainBin.png")'''