Extracting particles from an image: a toy model with numpy

Image segmentation can be performed by an image labelling process. Once segmented, it could necessary to isolate in independant images the labelled regions. Think of chromosomes in a metaphase image, if resolving chromosomes overlapps or karyotyping must be performed.

The following script takes "two images", (in fact two small numpy arrays : img is supposed to be a grey level image and lab, the label image which could be obtained after segmentation) and returns two particles, that is two images corresponding to the labelled regions (1 or 2. 0 corresponds to the background):

 

The left image corresponds to a labelled image. The two upper images were obtained with extractParticles() and the two lower images were extracted with:

scipy.ndimage.find_objects()

find_objects() doesn't extract the particles properly.

Update 2014-29-01
The code is on gist 
If ipython is not installed on your machine, play with it on cloudsagemath.

 #pseudo 16 bits gray level image
img=np.array([\
[75,167,181,105,33,5,2,2],\
[72,163,193,136,54,10,1,2],\
[54,123,154,118,51,11,2,2],\
[23,54,69,58,33,18,11,4],\
[5,15,27,43,62,69,53,23],\
[5,24,61,104,147,162,124,54],\
[10,54,128,186,220,222,168,72],\
[14,72,169,230,247,238,177,76]],np.int16)
 

#pseudo 16 bits label image
lab=np.array([\
[0,1,1,1,0,0,0,0],\
[0,1,1,1,1,0,0,0],\
[0,1,1,1,1,0,0,0],\
[0,1,1,1,0,0,0,0],\
[0,0,0,0,0,2,2,2],\
[0,0,0,2,2,2,2,2],\
[0,0,2,2,2,2,2,2],\
[0,2,2,2,2,2,2,2]],np.int16)
 

#
#March14, 2011, A script to extract regions (particles) from an array (image)
#Jean-Patrick Pommier
#
import numpy as np
import os
def extractParticles(grayIm,labIm):
    ''' give a grayscaled and a labelled image, extract the segmented particles
    ,returns a list of flattened particles'''
    #grayIm and labIm should have the same size
    def unflattenParticles(flatParticleList):
        '''take a list of flat particles and unflat them to yield an image'''
        unflatList=[]
        lenFlatList=len(flatParticleList)
        for i in range(0,lenFlatList):
            #get the i particle:current Particle
            curPart=flatParticleList[i]#current particle
            #x values(col) are stored in the third col (3-1)
            colmax=curPart[:,2].max()
            colmin=curPart[:,2].min()
            #y values(li) are stored in the fourth col (4-1)
            limax=curPart[:,3].max()
            limin=curPart[:,3].min()
            unflatIm=np.zeros((limax-limin+1,colmax-colmin+1),np.int16)
            #number of pixels in the particle
            nbPixel=len(curPart[:,1])#count how many lines at col=1
            for line in range(0,nbPixel):
                col=curPart[line,2]
                li=curPart[line,3]
                pixVal=curPart[line,1]
                unflatIm[li-limin,col-colmin]=pixVal
            unflatList.append(unflatIm)
        return unflatList
                   
    sx=grayIm.shape[0]
    sy=grayIm.shape[1]
    #flatten grayIm
    fg=grayIm.flatten()
    fl=labIm.flatten()
    labmax=fl.max()
    #print fg
    #print fl
    #build two 2D array containing x and y
    #of each pixel of the grayIm
    ax=np.zeros((sx,sy),np.int16)
    ay=np.zeros((sx,sy),np.int16)
    #vectorization with numpy may be 
    #more efficient than two loops
    for j in range(0,sy):
        for i in range(0,sx):
            ax[i,j]=j#filling ax with x=col
            ay[i,j]=i#filling ay with y values y=li
    #flat arrays of coordinates
    fax=ax.flatten()
    fay=ay.flatten()
    #1D merge graylevel, label and coordinates 
    #in one 1D array of 4-uplet
    extract=np.vstack((fl,fg,fax,fay))
    #transpose to watch it easily
    eT=extract.T
    #create a list of flatten particles
    #labIndex takes the value from 1 (the first particle to labmax the\
    #label of the last particle
    flatParticleList=[]#from Matthieu Brucher
    for labIndex in range(1,labmax+1):
        flatParticleList.append(eT[eT[:,0]==labIndex])#from Matthieu Brucher
    return unflattenParticles(flatParticleList)

if __name__ == "__main__":
    Particles=extractParticles(img,lab)
    for i in range(0,len(Particles)):
        print Particles[i]
        print    
 

Since that script is supposed to work on true images with a least 1536x1024 pixels, some part of that python script may be optimized with numpy, as:

#vectorization with numpy may be 
    #more efficient than two loops
    for j in range(0,sy):
        for i in range(0,sx):
            ax[i,j]=j#filling ax with x=col
            ay[i,j]=i#filling ay with y values y=li