from mahotas import labeled
from mahotas import _convex
from mahotas.thin import thin
import numpy as np
import pymorph as m
from pylab import imsave,imshow,show
import copy
import Image,ImageDraw,ImageFont
from scipy.spatial.distance import squareform,pdist
norm = lambda v: np.array(v)/np.linalg.norm(np.array(v))
angle = lambda x,y: np.arccos(np.dot(x,y)/(np.linalg.norm(x)*np.linalg.norm(y)))
[docs]class Cluster:
""" Container for a cell cluster
:param id: cluster id
:ivar cells: list of ids of the cells in the clusters
"""
def __init__(self,id):
self.id = id
self.cells = []
self.cellinterface = -1
self.ecminterface = -1
self.perimeter = -1
[docs] def addCell(self,cellid):
""" Add cell to cluster
:param cellid: id of cell
"""
if cellid not in self.cells:
self.cells.append(cellid)
[docs] def removeCell(self,cellid):
""" Remove cell from cluster
:param cellid: id of cell
"""
if cellid in self.cells:
self.cells.remove(cellid)
if len(self.cells) == 0:
return True
else:
return False
[docs] def getClusterSize(self):
""" Calculate number of cell in cluster
:return: number of cells in cluster
"""
#~ print len(self.cells)
return len(self.cells)
def realconvexhull(bwimg):
"""
Calcuate the real convex hull of an image. The default convex hull algorithms, part of mahotas, returns the coordinates of the pixels on the hull. Therefore, that convex hull intersects with the outer pixels. Here we calculate the hull over a set containing the 4 border points of each pixel.
:param bwimg: black-and-white image with 1 for occupied pixels and 0 at empty pixels
:type bwimg: numpy array
:return: numpy array of points in the hull
"""
# based on convexhull from mahotas.polygon
# extra points are added to prevent intersection between the hull and pixel centers
Y,X = np.where(bwimg)
ye = np.array([[2*y-1,2*y+1,2*y+1,2*y-1] for y in Y]).flatten()
xe = np.array([[2*x-1,2*x-1,2*x+1,2*x+1] for x in X]).flatten()
ind = np.lexsort((xe,ye))
P = list(zip(xe[ind],ye[ind]))
if len(P) <= 3:
return P
return np.array(_convex.convexhull(P))/2.0
[docs]def getCompactness(im,minval=0):
"""
Calculate compactness of a morphology: :math:`\\frac{A_{\\text{area}}}{A_{\\text{convex hull}}}`
:param im: image.
:type im: numpy array
:param minval: minimum value for valid pixels.
:type minval: int
:return: compactness
"""
# 1: create negative image
(nx,ny) = im.shape
imneg = np.uint8(im>minval)
# 2: find borders of negative image, these are all pixels where a neighbor is not cell
# thus, all borders are 2 pixels thick!
bim = labeled.borders(np.asarray(imneg))
# 3: get convex hull with intersection between bim and imneg -> the boundary pixels of the cells
# ph is an array of points on the hull
# Note that the input are only the pixels on the border of the image, thus speeds up the
# calculations A LOT.
ph = realconvexhull(imneg&bim)
# 4: calculate area of hull according to http://en.wikipedia.org/wiki/Polygon#Area_and_centroid
# note that first point on the hull is added to the end for this calculation
n = len(ph)
ph = np.vstack((ph,ph[0,:]))
A = -0.5*np.sum(ph[0:n,0]*ph[1:n+1,1]-ph[1:n+1,0]*ph[0:n,1])
# return compactness
return np.sum(imneg)/float(A)
[docs]def getLCC(im):
""" Find largest connected component of an image
:param im: input image
:type im: numpy array
:return: image with only largest component
"""
lab,nr = labeled.label(im)
sizes = np.bincount(np.reshape(lab,(lab.shape[0]*lab.shape[1])))
return lab==np.where(sizes==max(sizes[1:len(sizes)]))[0][0]
def getNodes(im,closedisk=5,thinloop=10):
"""
Find the branchpoints in the image. This is done in a number of steps:
1. Prepare image for skeleotanization with a closing operation that closes small holes and removes rough edges.
2. Create a skeleton with the thinning method of mahotas.
3. Find all the points in the skeleton that have more than 2 first order neighbors in the skeleton.
:param im: input image
:type im: numpy array
:param closedisk: size of the disk used for closing (step 1).
:type closedisk: int
:param thinloop: number of thinning iterations (step 2).
:type thinloop: int
:return: dictionary with nodes: {nodeid : (x,y)}, and skeleton image (nodelist,imout)
"""
#--- Make skeleton ---#
# 1) closing operation to get rid of rough edges
imout = m.close(im,m.sedisk(closedisk))
# 2) acutal thining
imout = m.thin(thin(imout), m.endpoints('loop'),thinloop)
# 3) find branchpoints
nodedict = {}
nz = np.where(imout==1)
for i in range(len(nz[0])):
x = nz[0][i]
y = nz[1][i]
nb = _getNNeighbors(x,y,imout)
if nb > 2:
nodedict[len(nodedict)] = (x,y)
return nodedict,imout
def mergeNodesAndEdges(nodes,edges,d=6):
""" Merge nodes that are too close and redefine the edges for these nodes.
:param nodes: dictionary with {nodeid : n(x,y)}.
:param edges: list of edges: [(nodeid_source,nodeid_target)].
:param d: minimum distance between nodes.
:type d: number
:return: dictionary with nodes: {nodeid : (x,y)} and list of edges
"""
nnodes = copy.deepcopy(nodes)
nedges = copy.deepcopy(edges)
D = squareform(pdist(np.array(nodes.values())))
visited = []
for i,node in nodes.iteritems():
if i in visited:
continue
nn = np.where(D[:,i]<d)[0]
if len(nn) > 1:
sx = 0.
sy = 0.
for n in nn:
visited.append(n)
sx += nodes[n][0]
sy += nodes[n][1]
newid = max(nnodes.keys())+1
nnodes[newid] = (int(sx/len(nn)),int(sy/len(nn)))
for n in nn:
if n in nnodes:
del nnodes[n]
remedges = []
newedges = []
for edge in nedges:
if (edge[0] in nn) or (edge[1] in nn):
if (edge[0] in nn) and (edge[1] not in nn):
newedges.append([edge[1],newid])
elif (edge[0] not in nn) and (edge[1] in nn):
newedges.append([edge[0],newid])
remedges.append(copy.deepcopy(edge))
for edge in remedges:
nedges.remove(edge)
for edge in newedges:
if edge not in nedges:
nedges.append(edge)
return nnodes,nedges
def purgeNodes(nodes,edges):
""" Remove nodes with only 2 edges.
:param nodes: dictionary with {nodeid : n(x,y)}.
:param edges: list of edges: [(nodeid_source,nodeid_target)].
:return: dictionary with nodes: {nodeid : (x,y)} and list of edges
"""
nnodes = copy.deepcopy(nodes)
nedges = copy.deepcopy(edges)
for id,node in nodes.iteritems():
degree = np.bincount((np.array(nedges)).flatten())
if id >= len(degree):
continue
if degree[id] == 2:
del nnodes[id]
pos = np.where(np.array(nedges)==id)
e0 = copy.deepcopy(nedges[pos[0][0]][:])
e1 = copy.deepcopy(nedges[pos[0][1]][:])
nedges.append([nedges[pos[0][0]][1-pos[1][0]],nedges[pos[0][1]][1-pos[1][1]]])
nedges.remove(e0)
nedges.remove(e1)
return nnodes,nedges
def getEndPoints(skel,nodes,d=6):
"""
Find the endpoints of a skeleton images. Endpoints are found with an iteration of the following steps:
1. Cut a given node N from the graph.
2. Label the resulting image to find all the connected components.
3. Iterate over the connected components and find the connected components without a node:
a. if a connected component has at least one node, an edge is connecting N to the rest of the image; no endpoint;
b. if a connected component has no nodes, it may be an edge to an endpoint.
4. For a connected component without nodes: check if there is a pixel in the neighbor that has only 1 first order neighbor and is further away than d.
5. One, or more endpoints can be found for each connected component, only save the endpoint that is farthest away from node N`.
:param skel: skeleton image
:type skel: numpy array
:param nodes: dictionary with {nodeid : (x,y)}
:param d: minimal distance between two branchpoints (4).
:type d: number
:return: list of endpoint positions
"""
im = copy.deepcopy(skel)
nof = 0
endpoints = []
edges = []
nof = 0
for nid,node in nodes.iteritems():
tmp = copy.deepcopy(im)
_delNode(tmp,node,w=3)
n = np.array([[1,1,1],[1,1,1],[1,1,1]])
lab,nr = labeled.label(tmp,n)
for i in range(1,nr+1):
noNodes = True
labim = lab==i
for n in nodes.values():
if (labim[n[0],n[1]] == 1):
noNodes = False
break
if noNodes:
nof += 1
ep = _findEndPoint(labim,nodes,d=6)
if len(ep) > 0:
dx = np.sqrt((node[0]-ep[0][0])**2+(node[1]-ep[0][1])**2)
idx = 0
for j,e in enumerate(ep):
if j == 0:
continue
if np.sqrt((node[0]-e[0])**2+(node[1]-e[1])**2) > dx:
dx = np.sqrt((node[0]-e[0])**2+(node[1]-e[1])**2)
idx = j
endpoints.append(ep[idx])
return endpoints
def _findEndPoint(im,nodes,d=6):
""" Find end points in an image. End points are defined as pixels with only 1 first order neighbor. This function is used by getEndPoints and is not intended for direct use. """
nz = np.where(im==1)
points = []
for i in range(len(nz[0])):
x = nz[0][i]
y = nz[1][i]
nb = _getNNeighbors(x,y,im)
if nb == 1:
isNode = False
for node in nodes.values():
if (np.sqrt((node[0]-x)**2+(node[1]-y)**2) < d):
isNode = True
break
if not isNode:
points.append((x,y))
return points
def getBlobs(im,minsize=10,opendisk=0):
""" Find blobs in an image.
:param im: input image
:type im: numpy array
:param minsize: minimal size of detected blobs
:type minsize: int
:param opendisk: disk size for opening operation (ommitted when opendisk <= 0)
:type opendisk: int
:return: identifiers for areas with blobs > minsize and labeled image (lac,lab)
"""
if opendisk > 0:
im = m.open(im.astype(np.int),m.sedisk(opendisk))
#~ n = np.array([[1,1,1],[1,1,1],[1,1,1]])
n = np.array([[0,1,0],[1,1,1],[0,1,0]])
lab,nr = labeled.label(im,n)
sizes = np.bincount(np.reshape(lab,(lab.shape[0]*lab.shape[1])))
lac = np.nonzero(sizes>minsize)[0]
return lac,lab
def getMeshes(im,minsize=20):
""" Find meshes (white areas) in an image.
:param im: input image
:type im: numpy array
:param minsize: minimal size of detected meshes
:type minsize: int
:return: list of mesh sizes
"""
lab,nr = labeled.label(im==0)
sizes = np.bincount(np.reshape(lab,(lab.shape[0]*lab.shape[1])))
# we are looking at white areas, this also includes the area around the image (this will be label=0)
lac = np.nonzero(sizes[1:nr+1]>minsize)[0]
return [sizes[i+1] for i in lac[1:len(lac)]]
def getEdges(nodedict,skel,w=3):
""" Find edges that connect the nodes in a skeleton. Edges are detected per node N with the following steps:
1. Add N to the skeleton and remove all other nodes;
2. Label image and continue with subimage that contains N;
3. All nodes that intersect with the image from (2) are connected to N;
:param nodedict: dictionary with {nodeid : (x,y)}
:param skel: skeleton image
:type skel: numpy array
:param w: node width that will be cut out of the first image
:type w: int
:returns: list of edges, values correspond with indices of nodelist
"""
# node id is position in list
edges = []
#~ for i,node in enumerate(nodelist):
for i,node in nodedict.iteritems():
connected = _getEdgesForNode(i,nodedict,skel,w)
for nd in connected:
e = [i,nd]
e.sort()
if e not in edges:
edges.append(e)
elif connected.count(nd) > 1:
edges.append(e)
return edges
def getNofEdges(nodes,skel,w=3):
""" Find number of edges for each node in nodelist.
:param nodes: dictionary with {nodeid : (x,y)}
:param skel: skeleton image
:type skel: numpy array
:param w: node width that will be cut out of the first image
:type w: int
:returns: dict with index of nodelist and number of edges
"""
return dict((n,_getNofEdgesForNode(n,nodes,skel,w)) for n,node in nodes.iteritems())
def _getNofEdgesForNode(idx,nodes,skel,w):
""" Collect all edges for a single node. This function is not intended for usage from outside but is called by getEdges """
# 1) create image of skeleton with node
im = copy.deepcopy(skel)
for i,node in nodes.iteritems():
if i is not idx:
_delNode(im,node,w=w)
_addNode(im,nodes[idx],w=w)
# 2) label image
n = np.array([[1,1,1],[1,1,1],[1,1,1]])
lab,nr = labeled.label(im,n)
nodeim = lab==lab[nodes[idx][0],nodes[idx][1]]
# 3) remove node from image
_delNode(nodeim,nodes[idx],w)
# 4) label image
lab,nr = labeled.label(nodeim,n)
return nr
def _getEdgesForNode(idx,nodedict,skel,w):
""" Collect all edges for a single node. This function is not intended for usage from outside but is called by getEdges """
#~ print idx
# 1) create image of skeleton with node
im = copy.deepcopy(skel)
for i,node in nodedict.iteritems():
if i is not idx:
_delNode(im,node,w=w)
_addNode(im,nodedict[idx],w=w)
#~ imshow(im)
#~ show()
# 2) label image
n = np.array([[1,1,1],[1,1,1],[1,1,1]])
lab,nr = labeled.label(im,n)
nodeim = lab==lab[nodedict[idx][0],nodedict[idx][1]]
#~ imshow(nodeim)
#~ show()
# 3) put nodes on labeled image with node idx and label again
for i,node in nodedict.iteritems():
if i is not idx:
_addNode(nodeim,node,w=w)
# 4) find nodes in nodelist that are in labeled image with idx
lab,nr = labeled.label(nodeim,n)
connim = lab==lab[nodedict[idx][0],nodedict[idx][1]]
#~ imshow(connim)
#~ show()
conn = []
for i,node in nodedict.iteritems():
if i == idx:
continue
if connim[node[0],node[1]] == 1:
conn.append(i)
return conn
def _addNode(im,pos,w=1):
# represent node as 3x3 block
nx = pos[0]+np.array([i for i in range(-w,w+1) for j in range(-w,w+1)])
ny = pos[1]+np.array([j for i in range(-w,w+1) for j in range(-w,w+1)])
#~ ny = pos[1]+np.array([-1,0,1,-1,1,0,-1,0,1])
pos = np.where(nx>=0)[0]
pos = pos[np.where(nx[pos]<im.shape[0])[0]]
pos = pos[np.where(ny[pos]>=0)[0]]
pos = pos[np.where(ny[pos]<im.shape[0])[0]]
im[nx[pos],ny[pos]] = 1
def _delNode(im,pos,w=1):
# represent node as 3x3 block
nx = pos[0]+np.array([i for i in range(-w,w+1) for j in range(-w,w+1)])
ny = pos[1]+np.array([j for i in range(-w,w+1) for j in range(-w,w+1)])
pos = np.where(nx>=0)[0]
pos = pos[np.where(nx[pos]<im.shape[0])[0]]
pos = pos[np.where(ny[pos]>=0)[0]]
pos = pos[np.where(ny[pos]<im.shape[0])[0]]
im[nx[pos],ny[pos]] = 0
def _getNNeighbors(x,y,im):
nx = x+np.array([-1,0,1,0])
ny = y+np.array([0,-1,0,1])
pos = np.where(nx>=0)[0]
pos = pos[np.where(nx[pos]<im.shape[0])[0]]
pos = pos[np.where(ny[pos]>=0)[0]]
pos = pos[np.where(ny[pos]<im.shape[0])[0]]
return np.sum(im[nx[pos],ny[pos]])
def _getNNeighbors2(x,y,im):
nx = x+np.array([-1,1,1,-1])
ny = y+np.array([-1,-1,1,1])
pos = np.where(nx>=0)[0]
pos = pos[np.where(nx[pos]<im.shape[0])[0]]
pos = pos[np.where(ny[pos]>=0)[0]]
pos = pos[np.where(ny[pos]<im.shape[0])[0]]
return np.sum(im[nx[pos],ny[pos]])
[docs]def getDirector(com,r,sigma,cells,weighted=False):
""" Find the director of the center of mass of a cell.
This implementation is best when the director is only calculated per cell, as done for the calculation of the order parameter. When the director is calculated per pixel use getDirectorPF.
:param com: center of mass of the cell [x,y]
:param r: radius of neighborhood
:type r: int
:param sigma: CPM grid
:type sigma: numpy array
:param cells: dictionary with cell id's as keys and :class:`~SimUtils.Cell.Cell` instances as values
:type cells: dict
:param weighted: average over cells (weighted=False) or over pixels (weighted=True)
:return: director (degrees)
"""
# Cut part we need out of sigma and create x- and y-grid in this range
xmin = 0 if (com[0]-r) < 0 else int(np.floor(com[0]-r))
xmax = sigma.shape[0] if (com[0]+r) > sigma.shape[0] else int(np.ceil(com[0]+r))
ymin = 0 if (com[1]-r) < 0 else int(np.floor(com[1]-r))
ymax = sigma.shape[1] if (com[1]+r) > sigma.shape[1] else int(np.ceil(com[1]+r))
sigma = sigma[xmin:xmax,ymin:ymax]
(x,y) = np.mgrid[xmin:xmax,ymin:ymax]
# calculate distances on the grid to find pixels within the circle
d = np.sqrt(np.power(x-com[0],2)+np.power(y-com[1],2))
neighborhood = np.unique(sigma[d<=r])
neighborhood = neighborhood[neighborhood>0]
# calculate average angle
a = []
if weighted:
w = np.bincount(sigma[d<=r])/(1.0*np.sum(b[neighborhood]))
else:
w = np.ones((len(neighborhood)))
for id in neighborhood:
a.append(cells[id].getAngle())
a = np.array(a)
a[a>90] = 180-a
return np.average(a,weights=w)
## Find the director of the com of a cell.
def getDirectorPF(com,r,sigma,angles,weighted=False):
""" Find the director of the center of mass of a cell.
This implementation is faster (for some reason) if the director for every pixel is calculated, if the director per cell is calculated use getDirector.
:param com: center of mass of the cell
:param r: radius of neighborhood
:param sigma: CPM grid
:type sigma: numpy array
:param angles: cell angle (degrees) at every coordinate (-1 for medium)
:type angles: numpy array
:param weighted: average over cells (weighted=False) or over pixels (weighted=True)
:return: director (degrees)
"""
# Cut part we need out of sigma and create x- and y-grid in this range
xmin = 0 if (com[0]-r) < 0 else int(np.floor(com[0]-r))
xmax = sigma.shape[0] if (com[0]+r) > sigma.shape[0] else int(np.ceil(com[0]+r))
ymin = 0 if (com[1]-r) < 0 else int(np.floor(com[1]-r))
ymax = sigma.shape[1] if (com[1]+r) > sigma.shape[1] else int(np.ceil(com[1]+r))
sigma = sigma[xmin:xmax,ymin:ymax]
angles = angles[xmin:xmax,ymin:ymax]
(x,y) = np.mgrid[xmin:xmax,ymin:ymax]
# calculate distances on the grid to find pixels within the circle
d = np.sqrt(np.power(x-com[0],2)+np.power(y-com[1],2))
sa = sigma[d<=r]
a = angles[d<=r]
an = a[sa>0]
if weighted:
if len(an) > 0:
return np.sum(an)/len(an)
else:
return np.sum(an)
else:
return np.mean([a[np.where(sa==cid)][0] for cid in np.unique(sa[np.where(sa>0)])])
# http://en.wikipedia.org/wiki/Liquid_crystal#Order_parameter
def getOrderParameter(cells,sigma,r,weighted=False):
""" Calculate order parameter for a morphology. Based on the given radius the order parameters is either calculated locally :func:`~SimUtils.AnalysisUtils.getLocalOrderParameter` or globally :func:`~SimUtils.AnalysisUtils.getGlobalOrderParameter`.
:param cells: cell objects
:type cells: dict
:param sigma: CPM grid
:type sigma: numpy array
:param r: radius of neighborhood
:type r: int
:param weighted: average over cells (weighted=False) or over pixels (weighted=True)
:type weighted: bool
.. seealso:: :func:`~SimUtils.AnalysisUtils.getLocalOrderParameter`, :func:`~SimUtils.AnalysisUtils.getGlobalOrderParameter`
"""
(nx,ny) = sigma.shape
rmax = np.sqrt(nx**2+ny**2)
if r < rmax:
return getLocalOrderParameter(cells,sigma,r,weighted=False)
else:
#~ print 'calculate global'
return getGlobalOrderParameter(cells,sigma,weighted=False)
[docs]def getOrderParameterFromGrid(sigma,angles,r,weighted=False):
""" Calculate order parameter for a morphology using the cpm grid data. Based on the given radius the order parameters is either calculated locally :func:`~SimUtils.AnalysisUtils.getLocalOrderParameter` or globally :func:`~SimUtils.AnalysisUtils.getGlobalOrderParameter`.
:param sigma: CPM grid
:type sigma: numpy array
:param angles: array of cell angles
:type angles: numpy array
:param r: radius of neighborhood
:type r: int
:param weighted: average over cells (weighted=False) or over pixels (weighted=True)
:type weighted: bool
.. seealso:: :func:`~SimUtils.AnalysisUtils.getLocalOrderParameterFromGrid`, :func:`~SimUtils.AnalysisUtils.getGlobalOrderParameterFromGrid`
"""
(nx,ny) = sigma.shape
if r < np.sqrt(nx**2+ny**2):
return getLocalOrderParameterFromGrid(sigma,angles,r,weighted)
else:
return getGlobalOrderParameterFromGrid(sigma,angles,weighted)
[docs]def getLocalOrderParameterFromGrid(sigma,angles,r,weighted=False):
""" Calculate local order parameter for a morphology using data from the cpm grid.
:param sigma: CPM grid
:type sigma: numpy array
:param angles: array of cell angles
:type angles: numpy array
:param r: radius of neighborhood
:type r: radius of neighborhood
:param weighted: average over cells (weighted=False) or over pixels (weighted=True)
:type weighted: bool
:return: local order parameter
.. seealso:: :func:`~SimUtils.AnalysisUtils.getDirector`
"""
s = 0
n = 0
for id in np.unique(sigma):
if id == 0:
continue
a = angles[np.where(sigma==id)][0]
#~ print np.min(a),np.max(a)
A = getDirectorPF(getCoM(np.where(sigma==id)),r,sigma,angles,weighted=weighted)
#~ print np.min(A),np.max(A)
#~ dA = abs(a-A) if abs(a-A) < 90 else 180-abs(a-A)
dA = abs(a-A) if abs(a-A) < 0.5*np.pi else np.pi-abs(a-A)
#~ s += (3*np.power(np.cos((2*np.pi*dA)/360.0),2)-1)/2.0
s += (3*np.power(np.cos(dA),2)-1)/2.0
n += 1
return s/n
[docs]def getGlobalOrderParameterFromGrid(sigma,angles,weighted=False):
""" Calculate local order parameter for a morphology using data from the cpm grid.
:param sigma: CPM grid
:type sigma: numpy array
:param angles: array of cell angles
:type angles: numpy array
:param weighted: average over cells (weighted=False) or over pixels (weighted=True)
:type weighted: bool
:return: global order parameter
"""
if weighted:
a = angles[np.where(sigma>0)]
else:
a = np.unique(angles[np.where(sigma>0)])
A = np.mean(a)
dA = np.abs(a-A)
dA[dA > .5*np.pi] = np.abs(dA[dA > .5*np.pi]-np.pi)
return np.mean((3*np.power(np.cos(dA),2)-1)/2.0)
#~ return np.mean((3*np.power(np.cos((2*np.pi*dA)/360.0),2)-1)/2.0)
# http://en.wikipedia.org/wiki/Liquid_crystal#Order_parameter
def getLocalOrderParameter(cells,sigma,r,weighted=False):
""" Calculate local order parameter for a morphology
:param cells: cell objects
:type cells: dict
:param sigma: CPM grid
:type sigma: numpy array
:param r: radius of neighborhood
:type r: int
:param weighted: average over cells (weighted=False) or over pixels (weighted=True)
:type weighted: bool
:return: local order parameter
.. seealso:: :func:`~SimUtils.AnalysisUtils.getDirector`
"""
n = 0
s = 0
angles = np.zeros(sigma.shape)
adict = dict((id,cell.getAngle()) for id,cell in cell.iteritems())
for id,cell in cells.iteritems():
if cell.type == 1:
continue
if cell.id in sigma:
angles[sigma==s] = adict[id]
for id,cell in cells.iteritems():
if cell.type == 1:
continue
a = adict[id]
if a > 90:
a = a-180
#~ A = getDirector(cell.getCoM(),r,sigma,cells,weighted)
A = getDirectorPF(cell.getCoM(),r,sigma,cells,angles,weighted=True)
dA = abs(a-A) if abs(a-A) < 90 else 180-abs(a-A)
#~ print a,A,dA
#~ print (2*np.pi*dA)/360.0
s += (3*np.power(np.cos((2*np.pi*dA)/360.0),2)-1)/2.0
n += 1
return s/n
# http://en.wikipedia.org/wiki/Liquid_crystal#Order_parameter
def getGlobalOrderParameter(cells,sigma,weighted=False):
""" Calculate local order parameter for a morphology
:param cells: cell objects
:type cells: dict
:param sigma: CPM grid
:type sigma: numpy array
:param weighted: average over cells (weighted=False) or over pixels (weighted=True)
:type weighted: bool
:return: global order parameter
"""
if weighted:
a = np.zeros(sigma.shape)
for id,cell in cells.iteritems():
if cell.type == 1:
continue
if cell.id in sigma:
a[sigma==s] = cell.getAngle()
else:
a = []
for id,cell in cells.iteritems():
if cell.type == 1:
continue
if cell.id in sigma:
a.append(cell.getAngle())
a = np.array(a)
a[a>90] = a[a>90]-180
A = np.mean(a)
dA = np.abs(a-A)
dA[dA > 90] = np.abs(dA[dA > 90]-180)
return np.mean((3*np.power(np.cos((2*np.pi*dA)/360.0),2)-1)/2.0)
[docs]def getCellAngle(pix):
""" Calculate angle of a cell on interval :math:`[0,\pi]`
:param pix: cell coordinates ([x1,...,xn],[y1,...,yn])
:return: angle in radians on interval :math:`[0,\pi]`
"""
l = getCellOrientation(pix)
a = np.arctan2(l[1],l[0])
if a < 0:
a += np.pi
return a
def getCellAngleFromOrientation(l):
""" Calculate angle of a cell using the cell orientation on interval :math:`[0,\pi]`
:param l: cell orientation [x,y]
:return: angle in radians on interval :math:`[0,\pi]`
"""
a = np.arctan2(l[1],l[0])
if a < 0:
#~ if a < 1:
a += np.pi
return a
def getCoM(pix):
""" Calculate center of mass of a cell
:param pix: cell coordinates ([x1,...,xn],[y1,...,yn])
:return: center of mass (x,y)
"""
try:
a = len(pix[0])
cx = 1/float(a)*(np.sum(pix[0]))
cy = 1/float(a)*(np.sum(pix[1]))
except:
print pix
return (cx,cy)
def getLargestCellComponent(sigma,cellid):
""" Identify connected components of a cell and return the pixels for the largest connected component.
:param sigma: CPM grid
:type sigma: numpy array
:param cellid: cell identifier
:type cellid: int
:return: tuple with lists of x and y pixels of the largest component
"""
imBW = np.uint8(sigma==cellid)
se = np.array([
[1, 1, 1],
[1, 1, 1],
[1, 1, 1]])
lab,nr = labeled.label(imBW,Bc=se)
sizes = np.bincount(np.reshape(lab,(lab.shape[0]*lab.shape[1])))
lidx = np.where(sizes==max(sizes[1:len(sizes)]))[0][0]
return np.where(lab==lidx)
def getCellInertiaTensor(pix):
""" Get inertia tensor for a cell
:param pix: cell coordinates ([x1,...,xn],[y1,...,yn])
:return: inertia tensor [[Ixx,Ixy],[Ixy,Iyy]]
"""
a = len(pix[0])
cx = 1/float(a)*(np.sum(pix[0]))
cy = 1/float(a)*(np.sum(pix[1]))
Ixx = np.sum(np.power(pix[0]-cx,2))
Iyy = np.sum(np.power(pix[1]-cy,2))
Ixy = -np.sum((pix[0]-cx)*(pix[1]-cy))
return np.array([[Ixx,Ixy],[Ixy,Iyy]])
def getCellAreas(sigma):
""" Calculate areas of cells.
:param sigma: CPM grid
:type sigma: numpy array
:return: list of cell areas
"""
a = []
for cid in np.unique(sigma):
if cid == 0:
continue
a.append(np.sum(sigma==cid))
return np.array(a)
def getCellLengths(sigma):
""" Calculate length of each cell.
:param sigma: CPM grid
:type sigma: numpy array
:return: list of cell length
"""
l = []
for cid in np.unique(sigma):
if cid == 0:
continue
l.append(getCellLength(np.where(sigma==cid)))
return np.array(l)
def getCellComponents(sigma):
""" Get all connected components of cells.
:param sigma: CPM grid
:type sigma: numpy array
:return: list of connected components
"""
c = []
for cid in np.unique(sigma):
if cid == 0:
continue
[lac,lab] = getBlobs(sigma==cid,minsize=1)
c.append(len(lac)-1)
return np.array(c)
def getCellLength(pix):
""" Calculate length of a cell using the inertia tensor.
:param pix: cell coordinates ([x1,...,xn],[y1,...,yn])
:return: cell length
.. seealso:: :func:`~SimUtils.AnalysisUtils.getCellInertiaTensor`
"""
a = len(pix[0])
[V,D] = np.linalg.eig(getCellInertiaTensor(pix))
return 4.0*np.sqrt(max(abs(V))/float(a))
def getCellOrientation(pix):
""" Calculate orientation of a cell. The orientation is the eigenvector corresponding to the largest eigenvalue of the cells' inertia tensor.
:param pix: cell coordinates ([x1,...,xn],[y1,...,yn])
:return: unit vector of the cell orientation
.. seealso:: :func:`~SimUtils.AnalysisUtils.getCellInertiaTensor`
"""
[V,D] = np.linalg.eig(getCellInertiaTensor(pix))
return [D[V==max(V)][0][0],D[V==max(V)][0][1]]
def getCellOrientationDiff(v1,v2):
""" Calculate angle between to cell orientations
:param v1: orientation vector of the first cell
:param v2: orientation vector of the second cell
:return: angle in radians on interval :math:`[0,\pi]`
"""
da = np.arccos(v1[0]*v2[0]+v1[1]*v2[1])
if da > .5*np.pi:
v1[0] = -1*v1[0]
v1[1] = -1*v1[1]
da = np.arccos(v1[0]*v2[0]+v1[1]*v2[1])
return da
def calcPerimeter(im):
""" Calculate perimeter length
:param im: boolean image; 1 for occupied pixels and 0 for empty pixels.
:type im: numpy array
:return: perimeter length
.. seealso:: :func:`~SimUtils.AnalysisUtils.getPerimeter`
"""
return np.sum(getPerimeter(im))
def getPerimeter(im):
""" Calculate perimeter of the object described by the non-zero pixels. This is done by finding the border pixels using using :func:`mahotas.labeled.borders` which returns an image with 1 at border pixels and zero at all other pixels. Note that :func:`mahotas.labeled.borders` finds border pixels at both sides of the interface, here we only use the pixels on the outside (im==0).
:param im: boolean image; 1 for occupied pixels and 0 for empty pixels.
:type im: numpy array
:return: image with 1 at border pixels and 0 at all other pixels.
"""
imtemp = copy.deepcopy(im)
imtemp[imtemp>0] = 1
bimall = labeled.borders(imtemp)
# labeled.borders returns pixels at both sides of the interface, we only use the inner perimeter.
bimall[im==0] = 0
return bimall
def _hasMediumNeighbor(pos,tau,medium):
return np.any(tau[pos[0]-1:pos[0]+2,pos[1]-1:pos[1]+2]==medium)
def _countMediumNeighbor(pos,tau,medium):
return np.sum(tau[pos[0]-1:pos[0]+2,pos[1]-1:pos[1]+2]==medium)
def countMediumNeighbor(cellim,tau,medium=0):
""" Count number of pixels with medium neighbors and the total number of medium pixels that are seen by the border pixels.
:param cellim: boolean image; 1 for occupied pixels and 0 for empty pixels.
:type cellim: numpy array
:param tau: array with cell types
:type tau: numpy array
:param medium: type of medium
:type medium: int
:return: number of cells that has at least one medium neighbor, total number of medium neighbors that pixels can be copied to.
"""
per = getPerimeter(cellim)
pix = np.where(per)
ncellsWithMedium = 0
totalMediumNeighbor = 0
for i in range(len(pix[0])):
n = _countMediumNeighbor([pix[0][i],pix[1][i]],tau,medium)
if n > 0:
ncellsWithMedium += 1
totalMediumNeighbor += n
return ncellsWithMedium,totalMediumNeighbor
def calcInterfaceTypes(im,tau):
""" Calculate perimeter of the object described by the non-zero pixels. This is done by finding the border pixels using using :func:`mahotas.labeled.borders` which returns an image with 1 at border pixels and zero at all other pixels. Note that :func:`mahotas.labeled.borders` finds border pixels at both sides of the interface, here we only use the pixels on the outside (im==0).
:param im: boolean image; 1 for occupied pixels and 0 for empty pixels.
:type im: numpy array
:param tau: array with cell types
:type tau: numpy array
:return: list with the type of each border pixel
"""
imtemp = copy.deepcopy(im)
imtemp[imtemp>0] = 1
bimall = labeled.borders(imtemp)
# labeled.borders returns pixels at both sides of the interface, we only use the outer perimeter.
bimall[im>0] = 0
return tau[bimall]
[docs]def calcClusterInterfaces(sigma,clusters,cells):
"""
Calculate length of cell-cluster interface and ECM-cluster interface.
For this we use the folowing perimeters:
- :math:`p_{\\text{M}}` = perimeter morphology
- :math:`p_{\\text{C}}` = perimeter cluster
- :math:`p_{\\text{NC}}` = perimeter of morppology WITHOUT the cluster
Not the total interface between the cluster and other cells is: :math:`p_\\text{cell,cluster} = \\frac{p_\\text{C}+p_\\text{NC}-p_\\text{M}}{2}\\;` and the interface between the cluster and ECM is: :math:`p_{\\text{ecm,cluster}} = p_\\text{C}-p_\\text{cell,cluster}`.
:param sigma: CPM grid
:type sigma: numpy array
:param clusters: dictionary with cluster id as key and :class:Cluster instaces as values
:param cells: list of :class:`~SimUtils.Cell.ClusterCell` instances
"""
# calcuate pM
pM = calcPerimeter(sigma)
# calculate perimeters per cluster
for id,c in clusters.iteritems():
imC = np.zeros_like(sigma)
imNC = copy.deepcopy(sigma)
imNC[imNC>0] = 1
for cid in c.cells:
imC[cells[cid].pix] = 1
imNC[cells[cid].pix] = 0
c.perimeter = calcPerimeter(imC)
c.cellinterface = .5*(c.perimeter+calcPerimeter(imNC)-pM)
c.ecminterface = c.perimeter-c.cellinterface
[docs]def getRelativeDirField(sigma,r,cells=None):
""" Calculate field with relative director for each pixel. The relative director is the difference to the angle of the cell at that pixel and the relative director on the pixel. Pixels with high values represent unordered areas, such as branchpoints.
:param sigma: CPM grid
:type sigma: numpy array
:param cells: dict with cell identifiers as keys and :class:`~SimUtils.Cell.ClusterCell` instances as values
:param r: radius used for calculating the director
:type r: int
:return: field with relative director values
.. seealso:: :func:`~SimUtils.AnalysisUtils.getDirectorPF`
.. seealso:: :func:`~SimUtils.AnalysisUtils.getDirectorPF`
"""
angles = getAngleField(sigma)
dir = np.zeros(sigma.shape)
# calculate director at every pixel
(nx,ny) = sigma.shape
for i in range(nx):
for j in range(ny):
if sigma[i,j] > 0:
dir[i,j] = getDirectorPF((i,j),r,sigma,angles,weighted=True)
# calcuate local difference between director and cell angles
a = copy.deepcopy(angles)
field = np.abs(dir-a)
field[np.where(field>0.5*np.pi)] = np.pi-field[np.where(field>0.5*np.pi)]
return field
[docs]def getPixelClusters(field,sigma,th=15,minlabsize=50,opendisk=1):
""" Get clusters of pixels for a single morphology. Clusters are detected with the following steps:
1) Remove all pixels that have a value in field higher than a given threshold.
2) Detect blobs in remaining image with a labeling algorith:
a. an opening operation may be performed before labeling.
b. areas smaller than a given size are ignored;
:param field: numpy array with values on which data is seperated
:param sigma: CPM grid
:param th: threshold value for step 1
:param minlabsize: labelled areas smaller than this value are ignored (2b)
:param opendisk: disk size for opening operation (2a)
:param mincellsize: minimal fraction of the cell that must be on the labelled area to be added to the cluster
:return: dictionary: {clusterid : list of cluster pixels}
.. seealso:: :func:`~SimUtils.AnalysisUtils.getBlobs`
"""
cfield = np.zeros_like(field)
#~ print np.unique(field)
cfield[field<th] = 120
cfield[sigma==0] = 0
labnum,lab = getBlobs(cfield,minsize=minlabsize,opendisk=opendisk)
clusters = {}
for n in labnum:
if n == 0:
continue
clusters[n] = np.where(lab==n)
return clusters
[docs]def getCellClusters(field,sigma,th=15,minlabsize=50,opendisk=1,mincellsize=.25):
""" Get clusters for a single morphology. Clusters are detected with the following steps:
1) Remove all pixels that have a value in field higher than a given threshold.
2) Detect blobs in remaining image with a labeling algorith:
a. an opening operation may be performed before labeling.
b. areas smaller than a given size are ignored;
3) Map each blob on the CPM grid:
a. at least a given fraction of the cell must be on the labeled area.
4) Check for cells in multiple clusters:
a. remove cell from all but biggest cluster;
b. remove cluster if it is empty after (a).
:param field: numpy array with values on which data is seperated
:param cells: dict with cell identifiers as keys and :class:`~SimUtils.Cell.ClusterCell` instances as values
:param sigma: CPM grid
:param th: threshold value for step 1
:param minlabsize: labelled areas smaller than this value are ignored (2b)
:param opendisk: disk size for opening operation (2a)
:param mincellsize: minimal fraction of the cell that must be on the labelled area to be added to the cluster
:return: dictionary with cluster id as key and :class:`~SimUtils.AnalysisUtils.Cluster` instances
.. seealso:: :func:`~SimUtils.AnalysisUtils.getBlobs`, :class:`~SimUtils.AnalysisUtils.Cluster`
"""
cfield = np.zeros_like(field)
cfield[field<th] = 120
cfield[sigma==0] = 0
celldict = dict((cellid,[]) for cellid in np.unique(sigma[sigma>0]))
labnum,lab = getBlobs(cfield,minsize=minlabsize,opendisk=opendisk)
clusters = {}
for n in labnum:
c = Cluster(n)
if n == 0:
continue
cellpos = list(sigma[np.where(lab==n)])
for cellid in np.unique(sigma[sigma>0]):
if cellpos.count(cellid) >= mincellsize*np.sum(sigma==cellid):
celldict[cellid].append(n)
c.addCell(cellid)
if c.getClusterSize() > 0:
clusters[n] = c
nolab = 0
for cellid,labels in celldict.iteritems():
if len(labels) == 0:
nolab += 1
elif len(labels) > 1:
lc = [clusters[lab].getClusterSize() for lab in labels]
cmax = labels[lc.index(max(lc))]
for lab in labels:
if lab is not cmax:
empty = clusters[lab].removeCell(cellid)
if empty:
del clusters[lab]
return clusters
[docs]def getAngleField(sigma):
""" Get field with the cell angles
:param sigma: CPM grid
:return: field with angles
"""
angles = np.zeros(sigma.shape)
for id in np.unique(sigma):
if id == 0:
continue
angles[sigma==id] = getCellAngle(np.where(sigma==id))
#~ angles[angles>90] = angles[angles>90]-180
return angles
[docs]def calcMSDTransForCellTC(com):
""" Calculate the translational MSD for a single object.
:param com: list of centers of mass at each time step
:return: list with MSD for each time step
"""
return np.power(com-com[0,:],2).sum(axis=1)
[docs]def calcMSDRotForCellTC(vecset):
""" Calculate the rotational MSD for a single object.
:param vecset: list of orientation vectors for each time step
:return: list with MSD for each time step
"""
# 1) Calculate angles between two vectors of consecutive time steps: $a$, this is used to find the best fitting orientation of each cell:
a = np.array([angle(vecset[i],vecset[i-1]) for i in range(1,len(vecset))])
# Find all indices in $a$ where the $a > \frac{1}{2}\pi$; here the opposite orientation should be used
# Rotate the vector corresponding to that angle 180 degrees.
# Update angles that are calculated with that vector.
# Repeat these steps untill all angles are correct.
# Now we have a set of vectors which in which all cells have the most likely orientation: vecset
while len(np.where(a>.5*np.pi)[0]):
i = np.where(a>.5*np.pi)[0][0]
vecset[i+1] = -1*vecset[i+1]
a[i] = angle(vecset[i+1],vecset[i])
if i+2 < len(vecset):
a[i+1] = angle(vecset[i+2],vecset[i+1])
# 2) Calculate the angle between the possitive x-axis and each vector
xa = np.arctan2(vecset[:,1],vecset[:,0])
# Remap these angles from the [0,-pi] domain to the [pi,2pi] domain.
xa[np.where(xa<0)] = 2*np.pi+xa[np.where(xa<0)]
# 3) Calculate differences between consecutive angles
# calculate differences between angles
da = xa[1:len(xa)] - xa[0:-1]
# detect transition over 2*np.pi and correct the angular differences
da[np.where(da<-1*np.pi)] = 2*np.pi+da[np.where(da<-1*np.pi)]
da[np.where(da>1*np.pi)] = -2*np.pi+da[np.where(da>1*np.pi)]
# summarize all angular differences to get a(t)-a(0) = sum(a(t)-a(t-1)) for eacht t
# 4) Calculate cumulative sum of the angles and square
sda = np.cumsum(da)
return np.power(np.concatenate(([0],sda)),2)
def getNeighborsForCell(sigma,cellid):
""" Get list of cell id's of the cells that are neighbors of cellid
:param sigma: CPM grid
:param cellid: id of the cell for which the neighbors are returned
:return: list with cell id's of the neighbors of cellid
"""
neighbormap = labeled.borders(sigma==cellid)
nb = np.unique(sigma[neighbormap])
return nb[(nb>0)&(nb!=cellid)]
def getNeighborsFromSigma(sigma):
""" Get neighbors based on the CPM grid
:param sigma: CPM grid
:return: dictionary with cellid:[neighbors]
"""
neighbors = {}
for cellid in np.unique(sigma):
if cellid == 0:
continue
neighbors[cellid] = getNeighborsForCell(sigma,cellid)
return neighbors
class TipFinder:
""" Class to find cells on the end points of sprouts
:param sigma: CPM grid
:param tau: type grid
:param nodes: dictionary with {nodeid : (x,y)}
:param edges: list with edges [(nodeid_soure,nodeid_target)]
"""
def __init__(self,sigma,tau,nodes,edges):
self._buildCellGraph(sigma,tau)
self._buildSkeletonGraph(sigma,nodes,edges)
def _buildCellGraph(self,sigma,tau):
from pygraph.classes.graph import graph
self.cg = graph()
for cellid in np.unique(sigma):
if cellid == 0:
continue
type = np.unique(tau[sigma==cellid])[0]
pos = getCoM(np.where(sigma==cellid))
attr = {'type':type,'pos':pos,'size':1,'color':(0,0,0)}
self.cg.add_node(cellid,attrs=[attr])
edges = getNeighborsFromSigma(sigma)
for id,neighbors in edges.iteritems():
for n in neighbors:
if not self.cg.has_edge((id,n)):
self.cg.add_edge((id,n))
def _buildSkeletonGraph(self,sigma,nodes,edges):
from pygraph.classes.graph import graph
self.sg = graph()
idmap = {}
for i,n in nodes.iteritems():
id = int(sigma[n[0],n[1]])
if id==0:
ids = np.bincount(sigma[(n[0]-5):(n[0]+5),(n[1]-5):(n[1]+5)].flatten())
ids = ids[1:len(ids)]
id = np.where(ids==np.max(ids))[0][0]+1
attr = {'pos':self.cg.node_attributes(id)[0]['pos'],'size':1,'color':(0,0,0)}
idmap[i] = id
if not self.sg.has_node(id):
self.sg.add_node(id,attrs=[attr])
for e in edges:
id1 = idmap[e[0]]
id2 = idmap[e[1]]
if not (id1 == id2):
if not self.sg.has_edge((id1,id2)):
#~ print (id1,id2)
self.sg.add_edge((id1,id2))
#~ print '-----'
def findTips(self,source,guess):
""" Find cells at the tip of the branch from source to guess
- Calculate the shortest path from source to all other nodes.
- Set the reference distance to the distance from source to guess (d) and add guess to tip nodes.
- Visit all neighbors of guess and calculate the shortest path from neighbor to source (spn).
* Check if we stumbled upon another node in the skeleton graph. If this happened we are walking over the wrong edge.
* If spn > d:
+ replace tip nodes with neighbor;
+ set d to spn;
+ clear list of nodes to visit next iteration (nnextNodes);
+ add all unvisited neighbors of neighbor nnextNodes.
* If spn = d:
+ add neighbor to tip nodes;
+ add all unvisited neighbors of neighbor nnextNodes.
:return: list of cell identifiers of the cells at the end point
"""
from pygraph.algorithms.minmax import shortest_path
sp = shortest_path(self.cg,source)[1]
if guess not in sp:
return []
d = sp[guess]
nextNodes = self.cg.neighbors(guess)
if len(nextNodes) < 2:
print source,guess,nextNodes
prevNodes = [guess]
tipNodes = [guess]
endpoints = []
wrongPath = False
while len(nextNodes) > 0:
nnextNodes = []
for n in nextNodes:
prevNodes.append(n)
for n in nextNodes:
if (n in self.sg.nodes()) and ((n is not guess) and (n is not source)):
#~ print 'I should not be here, break and see wat whe got'
wrongPath = True
break
if sp[n] > d:
d = sp[n]
nnextNodes = []
tipNodes = [n]
nadd = 0
for nn in self.cg.neighbors(n):
if nn not in prevNodes:
nnextNodes.append(nn)
nadd+=1
if nadd == 0:
endpoints.append(n)
elif sp[n] == d:
tipNodes.append(n)
nadd = 0
for nn in self.cg.neighbors(n):
if nn not in prevNodes:
nnextNodes.append(nn)
nadd+=1
if nadd == 0:
endpoints.append(n)
nextNodes = nnextNodes
if wrongPath:
return endpoints
else:
return tipNodes
def findAllCellsOnTips(self):
""" Find all cells at end points in the graph.
Iterate over all edges:
- If edge connects node of order 1 with another node
* the node with order one is used as the initial guess for the tip node;
* the other node is the source node, from where we start looking for the tip node;
* if both nodes have order 1, we look for the tip twice, once with the first node as source and the second as guess, and vice versa;
-
- for all edges that connect an end point (node with degree 1) to another node
:return: dictionary with (sigma,tau)
.. seealso:: :func:`~SimUtils.AnalysisUtils.TipFinder.findTips`
"""
tips = []
nt = 0
vsources = []
vguesses = []
for n1,n2 in self.sg.edges():
#~ print n1,n2
if ((n1 in vsources) and (n2 in vguesses)) or ((n2 in vsources) and (n1 in vguesses)):
continue
guess = None
source = None
twoWay = False
tipNodes = []
if (self.sg.node_order(n1) == 1) and (self.sg.node_order(n2) > 1):
guess = n1
source = n2
#~ tipNodes = self.findTips(n2,n1)
elif (self.sg.node_order(n2) == 1) and (self.sg.node_order(n1) > 1):
guess = n2
source = n1
#~ tipNodes = self.findTips(n1,n2)
elif (self.sg.node_order(n1) == 1) and (self.sg.node_order(n2) == 1):
guess = n1
source = n2
twoWay = True
if (guess is not None) and (source is not None):
if not twoWay:
tipNodes = self.findTips(source,guess)
vsources.append(source)
vguesses.append(guess)
else:
tipNodes = []
t1 = self.findTips(source,guess)
t2 = self.findTips(guess,source)
for t in t1:
tipNodes.append(t)
for t in t2:
tipNodes.append(t)
vsources.append(source)
vguesses.append(source)
vsources.append(guess)
vguesses.append(guess)
for n in tipNodes:
if n not in tips:
tips.append(n)
tipdict = {}
for t in tips:
if t not in tipdict:
tipdict[t] = self.cg.node_attributes(t)[0]['type']
return tipdict