Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
This python code implements the morphological "rose des directions" operator, based on morphological erosion. readmagick is used to load an image, another lib could be used.
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 7 15:19:19 2011
@author: Jean-Patrick Pommier
http://www.dip4fish.blogspot.com
"""
from __future__ import division
import numpy as np
from scipy import ndimage as nd
import readmagick
import pymorph
import mahotas
import mahotas.polygon
import pylab
import os
def HighPass(im,size):
blur=nd.gaussian_filter(im, size)
hi=pymorph.subm(im,blur)
return hi
def PeakByModalValue(array):
'''look for the modal value of a 2D array'''
x=array[0,:]
y=array[1,:]
print "y.shape",y.shape[0]
print "y[0]",y[0]
print "y[y.shape[0]-1]",y[y.shape[0]-1]
print "Searching modal value"
xmin=x.min()#image min graylevel
xmax=x.max()#image max gray level
mode=xmin
countmax=0#occurence of a given grayscale
#print "mig=",xmin," mag=",xmax
for i in range(0,y.shape[0]-1):
test=y[i]>countmax
#print "test:",test,"histo(",i,")=", y[i],"max",countmax
if test:
countmax=y[i]
mode=x[i]
#print "mode",mode
return mode
class particle(object):
rotatedIm=np.array([[0,0],[0,0]])
rotatedFlag=False
#ratio of particle area by convexhull area
#close to 1 for a convex particle
#lower for other as touching chromosomes
CvxhParticleArea_ratio=0
def __init__(self,arrayIm):
'''
Initially the rotated image is empty and the flag
indicates that no rotation is performed
'''
self.particuleImage=arrayIm
self.rotatedIm=np.array([[0,0],[0,0]])
self.compassTable=np.array([[0,0],[0,0]])
self.compassTablePeaks=np.array([[0,0],[0,0]])
def cvxhull_area(self):
'''
Calculate the convexhull area such that:
A=0.5*Sumfrom 0 to N-1 of {xn+1*yn-xn*yn+1}
Pn(xn,yn) and PN=P0
see http://en.wikipedia.org/wiki/Polygon#Area_and_centroid
'''
#print "cvxhull called"
binIm=self.particuleImage>0
area=np.sum(binIm[:,:]==True)
print "area",area
#print binIm.dtype
contour=mahotas.bwperim(binIm)
#print "contour",contour.dtype
pointlist=mahotas.polygon.convexhull(contour)
N=len(pointlist)
fP=pointlist[0]
#duplicate the first point P0
#at the end such that PointN=Point0
pointlist.append(fP)
s=0
#compute the sum from 0 to N-1
for i in range(0,N-1):
cx=pointlist[i][0]#x of the current point
cy=pointlist[i][1]#y of the current point
#print "Point",i," x=",cx," y=",cy
nx=pointlist[i+1][0]#x of the next point
ny=pointlist[i+1][1]#y of the next point
#print "Point suiv",i+1," x=",nx," y=",ny
det=nx*cy-cx*ny
s=s+det
#print "det:",det," S=",s
CvxhParticleArea_ratio=area/(0.5*abs(s))
return 0.5*abs(s),CvxhParticleArea_ratio
def getVerticalImage(self):
return self.rotatedIm
def getcompassTableDerivative(self):
return self.compassTablePeaks
def orientationByErosion(self,step):
''' Performs successive rotations of a binary particle from 0 to 180 by "step"
Then at each angle performs:
*successives erosions
*pixel counts
Find relationship between angle and the minimal erosion to destroy the particle
An anisotropic particle such a long chromosome should have a
principal orientation, overlapping chromosomes should have at least two
orientation.
'''
def compassTableDerivative(table):
deriv=np.zeros(table.shape)
x=table[0,:]
y=table[1,:]
for i in range(1,y.shape[0]-1):
deriv[1,i]=(y[i+1]-y[i])/((x[i+1]-x[i]))
deriv[0,i]=x[i]
print "delta y",y[i+1]-y[i]
print "delta x",x[i+1]-x[i]
print "der",deriv[1,i]
self.compassTablePeaks=np.copy(deriv)
return self.compassTablePeaks
def erodeVParticle(im):
'''Counts the number of vertical erosion neccessary
to destroy a particule
vline is a 3X1 structuring element
DestroyParticle returns the scalar number:nb_erosion
necessary to erode totally a binary particle'''
def countTruePix(bim):
'''bim must be a flat binary (True,False) array
counts the occurence of True (pixel) in a binary array/image '''
return np.sum(bim[:,:]==True)
vline=np.array([[0,1,0],
[0,1,0],
[0,1,0]])
#threshold the image
binaryIm=(im>0)
pixcount=countTruePix(binaryIm)
nb_erosion=1
while pixcount>0:
erodedIm=nd.binary_erosion(binaryIm,structure=vline,iterations=nb_erosion)
nb_erosion=nb_erosion+1
pixcount=countTruePix(erodedIm)
#erodedParticle work should be done here
return nb_erosion
#
#Analysis of the particle orientation starts here
#
maxRotationNubr=180/step
compassTable=np.zeros((2,maxRotationNubr+1),dtype=np.uint8)#rotation angle from 0 to 180 by step of 5°
#Store
#print compassTable.shape
i=0
angle=i*step
rotatedIm=self.particuleImage
#print self.particuleImage.shape
while i<=maxRotationNubr:
maxErosion=erodeVParticle(rotatedIm)
compassTable[0,i]=angle
compassTable[1,i]=maxErosion
#print "i=",i," angle=",angle," -erosion max:",maxErosion
rotatedIm=nd.rotate(self.particuleImage,angle)
i=i+1
angle=i*step
majorAngle=PeakByModalValue(compassTable)
self.rotatedIm=nd.rotate(self.particuleImage,majorAngle)
self.rotatedFlag=True
return majorAngle,compassTable,self.rotatedIm
user=os.path.expanduser("~")
#modify the path to your image
workdir=os.path.join(user,"Applications","ImagesTest","jp","Jpp48","13","DAPI","particules")
file="part14.png"
complete_path=os.path.join(workdir,file)
if __name__ == "__main__":
im=readmagick.readimg(complete_path)
hip=pymorph.subm(im,nd.gaussian_filter(im,3))
#print im.dtype#print uint16
p1=particle(im)
print "particle 1",p1.cvxhull_area()
p2=particle(hip)
print "particle 1 hi pass",p2.cvxhull_area()
theta1,rosedesvents,VImage=p1.orientationByErosion(10)
theta2,rdv,VHip=p2.orientationByErosion(10)
x=rosedesvents[0,:]
y=rosedesvents[1,:]
xh=rdv[0,:]
yh=rdv[1,:]
pylab.subplot(321,frameon=False, xticks=[], yticks=[])
pylab.gray()
pylab.imshow(im)
pylab.subplot(322,frameon=False, xticks=[], yticks=[])
pylab.imshow(hip)
pylab.subplot(323)
pylab.title("Resistance to vertical erosion")
pylab.ylabel("nb erosion")
pylab.xlabel("rotation angle")
pylab.plot(x,y,xh,yh)
pylab.subplot(324,frameon=False, xticks=[], yticks=[])
pylab.imshow(VImage)
pylab.subplot(326,frameon=False, xticks=[], yticks=[])
pylab.imshow(VHip)
pylab.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment