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")'''