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```