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)
pylab.imshow(img, interpolation='nearest')
show()

lab=np.array([\
[0,1,1,1,0,0,0,0],\
[0,1,1,1,1,0,0,0],\
[0,1,1,1,1,1,1,0],\
[0,1,2,2,0,0,1,0],\
[0,0,2,0,0,2,2,2],\
[0,0,2,2,2,2,2,2],\
[0,0,2,2,2,2,2,2],\
[0,2,2,2,2,2,2,2]],np.int16)
pylab.imshow(lab, interpolation='nearest')
show()

#
#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)


Particles=extractParticles(img,lab)
#
from scipy import ndimage as nd

seg2 = nd.find_objects(lab==2)
seg1 = nd.find_objects(lab==1)

subplot(242, xticks=[],yticks=[])
imshow(lab, interpolation='nearest')

subplot(243, xticks=[],yticks=[])
imshow(Particles[0]>0, interpolation='nearest')

subplot(244, xticks=[],yticks=[])
imshow(Particles[1]>0, interpolation='nearest')

subplot(247, xticks=[],yticks=[])
imshow(lab[seg1[0]], interpolation='nearest')

subplot(248, xticks=[],yticks=[])
imshow(lab[seg2[0]], interpolation='nearest')

show()