Thursday, April 7, 2011

Testing hit or miss with mahotas 0.6.4

Mahotas 0.6.4 is out. I upgrade my distribution with:
sudo easy_install --upgrade mahotas

Let's try to analyse a chromosome skeleton by hit and miss to detect end-points :
Top left:original image. Top right:thresholded image.Bottom left:skeleton. Bottom right:end-points
The skeleton of the binary image is build with pymorph.skelm(), then submitted to hit or miss transformation.
Strangely the skeleton is not a one connected component particle as it was the case here.
The python script is:
# -*- coding: utf-8 -*-
"""
Created on Thu Apr  7 10:08:16 2011

@author: Jean-Patrick Pommier
"""
import mahotas as mah
import pymorph as pm
import numpy as np
import pylab as plb
import os

user=os.path.expanduser("~")
workdir=os.path.join(user,"Applications","ImagesTest","jp","Jpp48","13","DAPI","particules")
file="part14.png"
#structuring elements
N=np.array([[0, 0, 0], [0, 1, 0], [0, 1, 0]])
NW=np.array([[0, 0, 0], [0, 1, 0], [1, 0, 0]])
W=np.array([[0, 0, 0], [1, 1, 0], [0, 0, 0]])
SW=np.array([[1, 0, 0], [0, 1, 0], [0, 0, 0]])
S=np.array([[0, 1, 0], [0, 1, 0], [0, 0, 0]])
SE=np.array([[0, 0, 1], [0, 1, 0], [0, 0, 0]])
E=np.array([[0, 0, 0], [0, 1, 1], [0, 0, 0]])
NE=np.array([[0, 1, 0], [0, 1, 0], [0, 0, 1]])
square=pm.asbinary(np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]]))
if __name__ == "__main__":
    complete_path=os.path.join(workdir,file)
    i1=mah.imread(complete_path)
    bi1=(i1>0)
    sk=pm.skelm(bi1)
    p0=mah.morph.hitmiss(sk,N)
    p1=mah.morph.hitmiss(sk,NW)
    p2=mah.morph.hitmiss(sk,W)
    p3=mah.morph.hitmiss(sk,SW)
    p4=mah.morph.hitmiss(sk,S)
    p5=mah.morph.hitmiss(sk,SE)
    p6=mah.morph.hitmiss(sk,E)
    p7=mah.morph.hitmiss(sk,NE)
    endpoints=p0+p1+p2+p3+p4+p5+p6+p7
    plb.gray()
    plb.subplot(2,2,1)
    plb.imshow(i1,interpolation=None)
    plb.subplot(2,2,2)
    plb.imshow(bi1,interpolation=None)
    plb.subplot(2,2,3)
    plb.imshow(sk,interpolation=None)
    plb.subplot(2,2,4)
    plb.imshow(endpoints,interpolation=None)
    plb.show()

Second trial: let's label the skeleton:
In the second trial, I only try to label the skeleton.


# -*- coding: utf-8 -*-
"""
Created on Thu Apr  7 10:08:16 2011

@author: Jean-Patrick Pommier
"""
from scipy import ndimage as nd
import mahotas as mah
import pymorph as pm
import numpy as np
import pylab as plb
import os

user=os.path.expanduser("~")
workdir=os.path.join(user,"Applications","ImagesTest","jp","Jpp48","13","DAPI","particules")
file="part14.png"
#structuring elements
N=np.array([[0, 0, 0], [0, 1, 0], [0, 1, 0]])
NW=np.array([[0, 0, 0], [0, 1, 0], [1, 0, 0]])
W=np.array([[0, 0, 0], [1, 1, 0], [0, 0, 0]])
SW=np.array([[1, 0, 0], [0, 1, 0], [0, 0, 0]])
S=np.array([[0, 1, 0], [0, 1, 0], [0, 0, 0]])
SE=np.array([[0, 0, 1], [0, 1, 0], [0, 0, 0]])
E=np.array([[0, 0, 0], [0, 1, 1], [0, 0, 0]])
NE=np.array([[0, 1, 0], [0, 1, 0], [0, 0, 1]])
square=pm.asbinary(np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]]))
if __name__ == "__main__":
    complete_path=os.path.join(workdir,file)
    i1=mah.imread(complete_path)
    bi1=(i1>0)
    sk=pm.skelm(bi1)
    sklabel,npart=nd.label(sk)
    print npart
    p0=mah.morph.hitmiss(sk,N)
    p1=mah.morph.hitmiss(sk,NW)
    p2=mah.morph.hitmiss(sk,W)
    p3=mah.morph.hitmiss(sk,SW)
    p4=mah.morph.hitmiss(sk,S)
    p5=mah.morph.hitmiss(sk,SE)
    p6=mah.morph.hitmiss(sk,E)
    p7=mah.morph.hitmiss(sk,NE)
    endpoints=p0+p1+p2+p3+p4+p5+p6+p7
    plb.gray()
    plb.subplot(2,2,1)
    plb.imshow(i1,interpolation=None)
    plb.subplot(2,2,2)
    plb.imshow(bi1,interpolation=None)
    plb.subplot(2,2,3)
    plb.imshow(sk,interpolation=None)
    plb.subplot(2,2,4)
    plb.jet()
    plb.imshow(sklabel,interpolation=None)
    plb.show()

2 comments:

Anonymous said...

How can you get end points coordinates?

dip4fish said...

hmmm, label binary end-points (or branched points) with mahotas.label() then numpy.where(label==k) should do the job (with k=1,2,...).