def getCellAngle(pix):
[docs] """ 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 getAngleField(sigma):
[docs] """ Get field with the cell angles
:param sigma: numpy array with cell id's
:return: numpy array with cell angles in radians
"""
angles = np.zeros(sigma.shape)
for id in np.unique(sigma):
if id == 0:
continue
angles[sigma==id] = getCellAngle(np.where(sigma==id))
return angles
def getOrderParameterFromGrid(sigma,angles,r):
[docs] """ Calculate order parameter for a morphology using the cpm grid data. Based on the given radius the order parameters is either calculated locally :func:`~OpenSimUtils.AnalysisUtils.getLocalOrderParameter` or globally :func:`~OpenSimUtils.AnalysisUtils.getGlobalOrderParameter`.
:param sigma: numpy array with cell id's
:param angles: numpy array with cell angles (radians)
:param r: radius of neighborhood
:type r: int
.. seealso:: :func:`~OpenSimUtils.AnalysisUtils.getLocalOrderParameterFromGrid`, :func:`~OpenSimUtils.AnalysisUtils.getGlobalOrderParameterFromGrid`
"""
(nx,ny) = sigma.shape
if r < np.sqrt(nx**2+ny**2):
return getLocalOrderParameterFromGrid(sigma,angles,r)
else:
return getGlobalOrderParameterFromGrid(sigma,angles)
def getLocalOrderParameterFromGrid(sigma,angles,r):
[docs] """ Calculate local order parameter.
:param sigma: numpy array with cell id's
:param angles: numpy array with cell angles (radians)
:param r: radius of neighborhood
:type r: int
:return: local order parameter
.. seealso:: :func:`~OpenSimUtils.AnalysisUtils.getDirector`
"""
s = 0
n = 0
for id in np.unique(sigma):
if id == 0:
continue
a = angles[np.where(sigma==id)][0]
A = getDirector(getCoM(np.where(sigma==id)),r,sigma,angles)
dA = abs(a-A) if abs(a-A) < 90 else 180-abs(a-A)
s += (3*np.power(np.cos((2*np.pi*dA)/360.0),2)-1)/2.0
n += 1
return s/n
def getGlobalOrderParameterFromGrid(sigma,angles):
[docs] """ Calculate global order parameter.
:param sigma: numpy array with cell id's
:param angles: numpy array with cell angles (radians)
:return: global order parameter
"""
a = angles[np.where(sigma>0)]
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)
def getDirector(com,r,sigma,angles):
[docs] """ Find the director of the center of mass of a cell.
:param com: center of mass of the cell (x,y)
:param r: radius of neighborhood
:type r: number
:param sigma: numpy array with cell id's
:param angles: numpy array with cell angles (radians)
:return: director (radians)
"""
# Because we only need the pixels within a radius r from com, we create new
# array sigma and angles that only contain pixels within 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 between all pixels and com
d = np.sqrt(np.power(x-com[0],2)+np.power(y-com[1],2))
# find all angles at pixels within radius r and sigma > 0
sa = sigma[d<=r]
a = angles[d<=r]
an = a[sa>0]
if len(an) > 0:
return np.sum(an)/len(an)
else:
return np.sum(an)
def getRelativeDirField(sigma,r):
[docs] """ 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: numpy array with cell id's
:param r: radius of neighborhood
:type r: int
:return: field with relative director values
.. seealso:: :func:`~OpenSimUtils.AnalysisUtils.getDirector`, :func:`~OpenSimUtils.AnalysisUtils.getAngleField`
"""
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] = getDirector((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
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
def getCompactness(sigma,minval=0):
[docs] """
Calculate compactness of a morphology: :math:`\\frac{A_{\\text{area}}}{A_{\\text{convex hull}}}` .
:param sigma: numpy array with cell id's
:param minval: minimum cell id for non medium pixels
:type minval: int
:return: compactness
"""
# 1: create negative image
(nx,ny) = sigma.shape
imneg = np.uint8(im>sigma)
# 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)
def getLCC(sigma):
[docs] """ Find largest connected component of an image
:param sigma: numpy array with cell id's
:return: image with largest connected component
"""
lab,nr = labeled.label(sigma)
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 getCellClusters(field,sigma,th=15,minlabsize=50,opendisk=1,mincellsize=.25):
[docs] """ 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:`~OpenSimUtils.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:`~OpenSimUtils.AnalysisUtils.Cluster` instances
.. seealso:: :func:`~OpenSimUtils.AnalysisUtils._getBlobs`, :class:`~OpenSimUtils.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
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([[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 calcMSDTransForCellTC(com):
[docs] """ 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)
def calcMSDRotForCellTC(vecset):
[docs] """ 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)
class Cluster:
[docs] """ 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
def addCell(self,cellid):
[docs] """ Add cell to cluster
:param cellid: id of cell
"""
if cellid not in self.cells:
self.cells.append(cellid)
def removeCell(self,cellid):
[docs] """ 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
def getClusterSize(self):
[docs] """ Calculate number of cell in cluster
:return: number of cells in cluster
"""
return len(self.cells)
class ClusterCellTC():
[docs] def __init__(self,id):
self.id = id
self.clusterId = np.array([],dtype=np.int)
self.clusterSize = np.array([],dtype=np.int)
self.time = np.array([],dtype=np.int)
self.laxis = np.empty((0,2),dtype=np.float)
self.com = np.empty((0,2),dtype=np.float)
def addTimeStep(self,t,pix,cid,csz):
[docs] self.time = np.append(self.time,(int(t)))
self.com = np.vstack((self.com,getCoM(pix)))
self.laxis = np.vstack((self.laxis,getCellOrientation(pix)))
self.clusterId = np.append(self.clusterId,(cid))
self.clusterSize = np.append(self.clusterSize,(csz))