In [1]:
"""
Created on Tue Oct 23 08:06:37 2012
@author: Jean-Patrick
"""
import os
import numpy as np
import matplotlib.pyplot as plt
import skimage.io as sio
import skimage as sk
import skimage.viewer as viewer
from skimage.draw import ellipse
from skimage.measure import find_contours, approximate_polygon, \
subdivide_polygon
import cmath
from math import pi
import itertools
def addpixelborder(im):
x,y = im.shape
dy = np.zeros((y,))
im = np.vstack((im,dy))
im = np.vstack((dy,im))
im = np.transpose(im)
x,y = im.shape
dy = np.zeros((y,))
im = np.vstack((im,dy))
im = np.vstack((dy,im))
im = np.transpose(im)
return im
def angle(p1,p2,p3):
'''get the angle between the vectors p2p1,p2p3 with p a point
'''
z1=complex(p1[0],p1[1])
z2=complex(p2[0],p2[1])
z3=complex(p3[0],p3[1])
v12=z1-z2
v32=z3-z2
#print v12,v32
angle=(180/pi)*cmath.phase(v32/v12)
return angle
def angleBetweenSegments(contour):
angles = []
points = list(contour)
points.pop()
#print 'contour in def',len(contour)
l = len(points)
for i in range(len(points)+1):
#print 'i',i,' i-1',i-1,' (i+1)%l',(i+1)%l
prv_point=points[(i-1)%l]
mid_point=points[i % l]
nex_point=points[(i+1)%l]
angles.append(angle(prv_point,mid_point,nex_point))
return angles
def quadArea(A,B,C,D):
x1 = A[0]
y1 = A[1]
x2 = B[0]
y2 = B[1]
x3 = C[0]
y3 = C[1]
x4 = D[0]
y4 = D[1]
#sign +/- if direct/indirect quadrilateral
return 0.5*(x1*y2+x2*y3+x3*y4+x4*y1-x2*y1-x3*y2-x4*y3-x1*y4)
def maximalQuad(A,B,C,D):
maxarea = 0
for quad in itertools.permutations((A,B,C,D)):
area = quadArea(quad[0],quad[1],quad[2],quad[3])
if area > maxarea:
maxarea = area
maxquad = quad
return maxquad
In [18]:
user=os.path.expanduser("~")
#workdir=os.path.join(user,"QFISH","JPPAnimal","JPP52","11","DAPI","particles")
#image='part25.png'
#QFISH/CytoProject/Jpp48/1/DAPI/particles/part9.png
workdir=os.path.join(user,"QFISH","CytoProject","Jpp48","1","DAPI","particles")
image='part9.png'
complete_path=os.path.join(workdir,image)
#print complete_path
im = sio.imread(complete_path)
im = addpixelborder(im)
original_im = np.copy(im)
fig, ax1= plt.subplots(ncols=1, figsize=(10, 10))
#subplot(111,xticks=[],yticks=[])
contours = find_contours(im, 0.01)
p = subdivide_polygon(contours[0], degree=5, preserve_ends=True)
#coord = approximate_polygon(p, tolerance=1.5)
coord2 = approximate_polygon(p, tolerance=2)
ax1.imshow(im, interpolation='nearest', cmap=plt.cm.gray)
contour0=contours[0]
pf=contours[-1]
angles = angleBetweenSegments(coord2)
table = np.array([coord2[0][0],coord2[0][1],angles[0]])
for i in range(len(angles[1:])):
current = np.array([coord2[i][0],coord2[i][1],angles[i]])
table = np.vstack((table,current))
pos = table[np.where(np.logical_and(table[:,2]>80, table[:,2]<145))]
neg = table[np.where(np.logical_and(table[:,2]>-130, table[:,2]<0))]
flat = table[np.where(np.logical_and(table[:,2]>=-180, table[:,2]<=-140))]
for n, contour in enumerate(contours):
ax1.plot(contour[:, 1], contour[:, 0], linewidth=4)
ax1.plot(coord2[:, 1], coord2[:, 0], linewidth=1, color='r')
#ax1.plot(coord2[:, 1], coord2[:, 0], linewidth=1, color='y')
ax1.scatter(contour0[0][1],contour0[0][0],s=150, c='m',marker='o')
ax1.scatter(pf[0][1],pf[0][0],s=150, c='y',marker='*')
#print table[pos]
ax1.scatter(pos[:,1], pos[:,0], s=200, c='r')
ax1.scatter(neg[:,1], neg[:,0], s=200, c='y')
ax1.scatter(flat[:,1], flat[:,0], s=200, c='g')
plt.show()
neg is a table containing points coordinates with a negative angle between vectors starting from that point¶
In [4]:
neg[:,0:3]
Out[4]:
array([[ 117.03124939, 82.00000057, -114.97592621],
[ 98.00000052, 84.96875058, -74.09936894],
[ 89.00000065, 104.03124929, -97.90399103],
[ 116.99999939, 106.03124933, -102.89160021]])
In [5]:
neg[0,0:2]
Out[5]:
array([ 117.03124939, 82.00000057])
Let's make combinations of four points¶
In [6]:
coord = []
for r in range(0, neg[:,0:2].shape[0]):
coord.append(neg[r,0:2])
quad_it = itertools.combinations(coord,4)
q=next(quad_it)
q
Out[6]:
(array([ 117.03124939, 82.00000057]), array([ 98.00000052, 84.96875058]), array([ 89.00000065, 104.03124929]), array([ 116.99999939, 106.03124933]))
maxquad is the permutation of four points of maximal area, it corresponds to a direct quadrilateral:¶
In [7]:
maxQuad = maximalQuad(*q)
print maxQuad
print quadArea(*maxQuad)
(array([ 117.03124939, 82.00000057]), array([ 116.99999939, 106.03124933]), array([ 89.00000065, 104.03124929]), array([ 98.00000052, 84.96875058])) 504.500921106
In [8]:
from matplotlib.patches import Polygon
p = Polygon(maxQuad, alpha=0.5, color='b' )
fig, ax1= plt.subplots(ncols=1, figsize=(8, 8))
ax1.imshow(np.transpose(im), interpolation='nearest', cmap=plt.cm.gray)
ax1.scatter(neg[:,0], neg[:,1], s=100, c='y')
ax1.add_artist(p)
Out[8]:
<matplotlib.patches.Polygon at 0xb3603ac>
In [9]:
from skimage.draw import polygon
print im.shape
(164, 128)
With skimage, let's convert the quadrilateral coordinates into a quadrilateral image¶
In [10]:
overlapp = np.zeros(im.shape, dtype=int)
print overlapp.shape
x = np.asarray([int(r[0]) for r in maxQuad])
y = np.asarray([int(r[1]) for r in maxQuad])
print x
print y
rr,cc = polygon(x, y)
overlapp[rr,cc] = True
cut = np.logical_and(im>0,np.logical_not(overlapp))
subplot(121, xticks=[],yticks=[])
figsize(12,12)
imshow(cut, interpolation ='nearest')
subplot(122, xticks=[],yticks=[])
title('overlapp domain')
imshow(overlapp, interpolation ='nearest')
(164, 128) [117 116 89 98] [ 82 106 104 84]
Out[10]:
<matplotlib.image.AxesImage at 0xb1aad4c>
The overlapping domain is removed from the original image to yield chromosomal elements. Those elements are labelled, too small elements are merged to the overlapping domain:¶
In [11]:
import mahotas as mh
labcut = sk.morphology.label(cut, neighbors=4)
labsize = mh.labeled.labeled_size(labcut)
#print labsize
#print labsize < 2
#print np.where(labsize < 2)
smallimage = labcut==np.where(labsize < 2)
labcut = mh.labeled.remove_regions(labcut,4)
mh.labeled.relabel(labcut, inplace=True)
overlapp = (labcut.max()+1)* (smallimage+overlapp > 0)
figsize(4,4)
subplot(121,xticks=[],yticks=[])
title(str(mh.labeled.labeled_size(labcut)[1:]))
#print np.histogram(labcut, bins=1)
imshow(labcut, interpolation = 'nearest')
subplot(122,xticks=[],yticks=[])
title('overlapp lab:'+str(overlapp.max()))
imshow(smallimage+overlapp, interpolation = 'nearest')
Out[11]:
<matplotlib.image.AxesImage at 0xb62f22c>
In [12]:
from matplotlib import colors
testim = np.array([[0,1,1,0],[0,2,0,0],[0,0,3,3],[4,0,0,0]])
elements = [labcut==l for l in range(1,labcut.max()+1)]
#print len (elements)
figsize(6,6)
elements2=[]
mycmap = colors.ListedColormap(['black','red','cyan','blue','green','yellow'])
bounds = [0,1,2,3,4,5]
mynorm = colors.BoundaryNorm(bounds, mycmap.N+1)
for imlab in range(1,labcut.max()+1):
cur = labcut == imlab
cur = cur * imlab
#print imlab, cur.min(), cur.max()
elements2.append(cur)
subplot(2,2,imlab, xticks=[],yticks=[])
title(str(cur.max())+' label='+str(imlab))
#print cur.max()
imshow(cur, interpolation='nearest',cmap=mycmap, norm=mynorm)
Supposing that there are two chromosomes, let's make combinations of of two elements in the set of four elements. Two of those combinations are in a list called assemble:¶
In [13]:
candidates_it = itertools.combinations(elements2,2)
candidates = [c for c in candidates_it]
assemble = []
for c in candidates:
image = np.zeros(im.shape, dtype = int)
for i in c:
image = image + i
image = image+overlapp
#image = image > 0
assemble.append(image)
In [14]:
figsize(5,5)
convexhull00 = sk.morphology.convex_hull_image(assemble[0])
contour00 = mh.bwperim(convexhull00)
subplot(111, xticks=[],yticks=[])
imshow(2*contour00+assemble[0], interpolation = 'nearest')
Out[14]:
<matplotlib.image.AxesImage at 0xb5de04c>
In [15]:
def contourNegCorners(im, angle0=-130, angle1 = 0):
contours = find_contours(im, 0.01)
subpol = subdivide_polygon(contours[0], degree=5, preserve_ends=True)
coord = approximate_polygon(subpol, tolerance=1.5)
angles = angleBetweenSegments(coord)
table = np.array([coord[0][0],coord[0][1],angles[0]])
for i in range(len(angles[1:])):
current = np.array([coord2[i][0],coord2[i][1],angles[i]])
table = np.vstack((table,current))
neg = table[np.where(np.logical_and(table[:,2]>angle0, table[:,2]<angle1))]
return neg
Let's visualize all combinations of two elements+overlapp domain.¶
In [16]:
from mahotas.polygon import convexhull as mhcvxh
square = np.array([[1, 1, 1],[1, 1, 1], [1, 1, 1]])
print len(assemble)
n=1
contour = mh.bwperim(im > 0)
figsize(12,12)
chromosomes = {}
for im in assemble:
subplot(2,3,n,xticks=[],yticks=[])
cvxh = sk.morphology.convex_hull_image(im)
ratio =100 * ((1.0*np.sum(im>0)) / (1.0*np.sum(cvxh > 0)))
lab,_ = mh.label(im>0, Bc=square)
prop = sk.measure.regionprops(lab,properties=['Eccentricity','ConvexArea','Area'])
Ecc = prop[0]['Eccentricity']
Area = prop[0]['Area']
CVXArea = prop[0]['ConvexArea']
#print Ecc, str(100.0*Area/CVXArea)[0:5]
labim, nim = mh.label(im)
negcorners = contourNegCorners(im>0)
#
cont = sk.measure.find_contours(im,0.01)
pol = subdivide_polygon(cont[0], degree=5, preserve_ends=True)
coord_pol = approximate_polygon(pol, tolerance=3)
angles_pol = angleBetweenSegments(coord_pol)
#print angles_pol
table_im = np.array([coord_pol[0][0],coord_pol[0][1],angles_pol[0]])
for i in range(len(angles_pol[1:])):
current_im = np.array([coord_pol[i][0],coord_pol[i][1],angles_pol[i]])
table_im = np.vstack((table_im,current_im))
neg_dom = table_im[np.where(np.logical_and(table_im[:,2]>-110, table_im[:,2]<0))]
#Build a dict to filter assembled chromosomes
chromosomes[ratio] = im
title('cc='+str(nim)+' r='+str(ratio)[0:5]+' Ecc:'+str(100*Ecc)[0:4]+' corner='+str(len(neg_dom)))
imshow(im, interpolation = 'nearest',cmap=mycmap)
convexhull = mhcvxh(im>0)
plot(convexhull[:,1], convexhull[:,0],linewidth=1, color='w')
plot(coord_pol[:, 1], coord_pol[:, 0], linewidth=2, color='c')
if len(neg_dom)>0:
scatter(neg_dom[:,1],neg_dom[:,0],s=100,c='white')
#scatter(negcorners[:,1], negcorners[:,0], s=100, c='y')
n =n +1
6
Let's extract the two chromosomes from the chromosomes dictionary¶
In [25]:
chromkeys = chromosomes.keys()
chromkeys.sort()
twoHighestRatio = chromkeys[-2:]
r1=twoHighestRatio[0]
r2=twoHighestRatio[1]
chrom1=chromosomes[r1]
chrom2=chromosomes[r2]
figsize(10,10)
subplot(221,xticks=[],yticks=[])
imshow(chrom1, interpolation='nearest')
subplot(222,xticks=[],yticks=[])
segmented1 = (chrom1>0)*original_im
imshow(segmented1, interpolation='nearest', cmap=pylab.cm.gray)
subplot(223,xticks=[],yticks=[])
imshow(chrom2, interpolation='nearest')
subplot(224,xticks=[],yticks=[])
segmented2 = (chrom2>0)*original_im
imshow(segmented2, interpolation='nearest', cmap=pylab.cm.gray)
Out[25]:
<matplotlib.image.AxesImage at 0xca71cac>
In [ ]: