Wednesday, April 20, 2011

Hit or miss: skeleton and convex hull with Python

The following script is used to detect both branching points (green) and endpoints (blue)in skeletons (red) of touching and overlapping chromosomes, before and after high pass filtering of the image of the chromosomes.

Hit and miss morphological operator is also used to compute the convex hull of a particle. I try to implement the algorithm without success, the convex hull calculated (red) is not different from the initial binary particle (white):
Upper and middle image: Skeleton analysis. Lower image: fail to produce a convex hull (red)
 The set of structuring elements used in the convex hull function, corresponds to all the possible corners in 3x3 pixel array.


Convex hull algorithm (II) :
The first trial was not successful due to a boolean type issue. I tried to implement the algorithm again according to the definition found in Gonzalez & Woods, "Digital Image Processing", using the hit or miss operator provided by the mahotas library.
Lower image:result in red of the convex hull found
Clearly the function ConvexHull_2, I implemented, gives a wrong result. 
For comparison, the convex hull obtained with ImageJ is:
Convex hull (blue)
De DIP4FISH
However, the result seems to be correct on the "west" side or on the "east" side of the image whereas "North" and "South" are wrong. The structuring elements may be implicated, or not .... (to be continued).

def ConvexHull_2(bim):
    ''' 
    From Gonzalez & Woods p545 (Addison Wesley 1992), the algorithm is:
        X(i,0)=A        
        X(i,k)=[X(i,k-1)*B(i)] u A
        D(i)=X(i,conv) conv such X(i,k)=X(i,k-1)
        k=1,2,3,....
        i=1,2,3,4 :one of the four structuring element B of type III
        convexhull of A: C(A)=u(from i=1 to 4) D(i)
        *: hit or miss operator
        u:union
    '''
    #
    #           3x3 structuring elements: 
    #   0:no pixels; 1: one pixel; 2:don't care
    #
    cornerS=np.array([[2,2,2],
                      [2,0,2],
                      [1,1,1]])
    
    cornerE=np.array([[2,2,1],
                      [2,0,1],
                      [2,2,1]])
    
    cornerN=np.array([[1,1,1],
                      [2,0,2],
                      [2,2,2]])
    
    cornerW=np.array([[1,2,2],
                      [1,0,2],
                      [1,2,2]])
    def countTruePix(bim):
        '''bim must be a binary array 
      Counts the occurence of True (pixel)'''
        return np.sum(bim[:,:]==True)
    
    def HitMissConfigUpToConv(bim,config):
        '''
        config is the "i" of the algo
        '''
        if config=="N":B=cornerN
        if config=="W":B=cornerW
        if config=="S":B=cornerS
        if config=="E":B=cornerE
        X=bim
        idempotence=False
        k=1
        while not(idempotence):
            size=countTruePix(X)
            hitpoints=mah.morph.hitmiss(X.astype(bool),B)
            X=X+hitpoints
            newsize=countTruePix(X)
            #size==newsize supposes that idempotence (convergence) is reached
            #this may not be true:same size but different shape
            idempotence=(newsize-size==0)
            k=k+1
        return X
    #starts here
    A=bim
    D1=HitMissConfigUpToConv(A,"N")
    D2=HitMissConfigUpToConv(A,"W")
    D3=HitMissConfigUpToConv(A,"S")
    D4=HitMissConfigUpToConv(A,"S")
    return D1+D2+D3+D4