Newer
Older
Base functions to extract information from a single cell.
These functions are automatically read by extractor.py, and
so can only have the cell_mask and trap_image as inputs. They
must return only one value.
They assume that there are no NaNs in the image.
We use the module bottleneck when it performs faster than numpy:
- median
- values containing NaNs (but we make sure this never happens).
import math
import typing as t
import bottleneck as bn
import numpy as np
from scipy import ndimage
def area(cell_mask) -> int:
def eccentricity(cell_mask) -> float:
Find the eccentricity using the approximate major and minor axes.
min_ax, maj_ax = min_maj_approximation(cell_mask)
return np.sqrt(maj_ax**2 - min_ax**2) / maj_ax
def mean(cell_mask, trap_image) -> float:
Find the mean of the pixels in the cell.
return np.mean(trap_image[cell_mask])
def total(cell_mask, trap_image) -> float:
"""
Find the sum of the pixels in the cell.
Parameters
----------
cell_mask: 2d array
Segmentation mask for the cell.
trap_image: 2d array
"""
return np.sum(trap_image[cell_mask])
def total_squared(cell_mask, trap_image) -> float:
"""
WARNING: produces overflow error when converted to float16
Find the sum of the square of the pixels in the cell.
For finding variances.
Parameters
----------
cell_mask: 2d array
Segmentation mask for the cell.
trap_image: 2d array
"""
return np.sum(trap_image[cell_mask] ** 2)
def median(cell_mask, trap_image) -> int:
Find the median of the pixels in the cell.
return bn.median(trap_image[cell_mask])
def max2p5pc(cell_mask, trap_image) -> float:
Find the mean of the brightest 2.5% of pixels in the cell.
npixels = np.sum(cell_mask)
n_top = int(np.ceil(npixels * 0.025))
# sort pixels in cell and find highest 2.5%
top_values = bn.partition(pixels, len(pixels) - n_top)[-n_top:]
return np.mean(top_values)
def max5px_median(cell_mask, trap_image) -> float:
Estimate the degree of localisation.
Find the mean of the five brightest pixels in the cell divided by the
median of all pixels.
pixels = trap_image[cell_mask]
if len(pixels) > 5:
top_values = bn.partition(pixels, len(pixels) - 5)[-5:]
# find mean of five brightest pixels
max5px = np.mean(top_values)
med = np.median(pixels)
if med == 0:
return np.nan
else:
return max5px / np.median(pixels)
Find the standard deviation of the values of the pixels in the cell.
return np.std(trap_image[cell_mask])
def volume(cell_mask) -> float:
Estimate the volume of the cell.
Assumes the cell is an ellipsoid with the mask providing
a cross-section through its median plane.
min_ax, maj_ax = min_maj_approximation(cell_mask)
Parameters
----------
cell_mask: 2D array
Segmentation mask for the cell
"""
padded = np.pad(cell_mask, 1, mode="constant", constant_values=0)
nearest_neighbor = (
ndimage.morphology.distance_transform_edt(padded == 1) * padded
)
return 4 * np.sum(nearest_neighbor)
Estimate the volume of the cell.
Assumes the cell is a sphere with the mask providing
a cross-section through its median plane.
Parameters
----------
cell_mask: 2d array
Segmentation mask for the cell
"""
total_area = area(cell_mask)
r = math.sqrt(total_area / np.pi)
return (4 * np.pi * r**3) / 3
def min_maj_approximation(cell_mask) -> t.Tuple[int]:
Find the lengths of the minor and major axes of an ellipse from a cell mask.
Parameters
----------
cell_mask: 3d array
Segmentation masks for cells
"""
# pad outside with zeros so that the distance transforms have no edge artifacts
padded = np.pad(cell_mask, 1, mode="constant", constant_values=0)
# get the distance from the edge, masked
nn = ndimage.morphology.distance_transform_edt(padded == 1) * padded
# get the distance from the top of the cone, masked
dn = ndimage.morphology.distance_transform_edt(nn - nn.max()) * padded
# get the size of the top of the cone (points that are equally maximal)
cone_top = ndimage.morphology.distance_transform_edt(dn == 0) * padded
# minor axis = largest distance from the edge of the ellipse
min_ax = np.round(np.max(nn))
# major axis = largest distance from the cone top
# + distance from the center of cone top to edge of cone top
maj_ax = np.round(np.max(dn) + np.sum(cone_top) / 2)
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
def moment_of_inertia(cell_mask, trap_image):
"""
Find moment of inertia - a measure of homogeneity.
From iopscience.iop.org/article/10.1088/1742-6596/1962/1/012028
which cites ieeexplore.ieee.org/document/1057692.
"""
# set pixels not in cell to zero
trap_image[~cell_mask] = 0
x = trap_image
if np.any(x):
# x-axis : column=x-axis
columnvec = np.arange(1, x.shape[1] + 1, 1)[:, None].T
# y-axis : row=y-axis
rowvec = np.arange(1, x.shape[0] + 1, 1)[:, None]
# find raw moments
M00 = np.sum(x)
M10 = np.sum(np.multiply(x, columnvec))
M01 = np.sum(np.multiply(x, rowvec))
# find centroid
Xm = M10 / M00
Ym = M01 / M00
# find central moments
Mu00 = M00
Mu20 = np.sum(np.multiply(x, (columnvec - Xm) ** 2))
Mu02 = np.sum(np.multiply(x, (rowvec - Ym) ** 2))
# find invariants
Eta20 = Mu20 / Mu00 ** (1 + (2 + 0) / 2)
Eta02 = Mu02 / Mu00 ** (1 + (0 + 2) / 2)
# find moments of inertia
moi = Eta20 + Eta02
return moi
else:
return np.nan
def ratio(cell_mask, trap_image):
"""Find the median ratio between two fluorescence channels."""
if trap_image.ndim == 3 and trap_image.shape[-1] == 2:
fl_0 = trap_image[..., 0][cell_mask]
fl_1 = trap_image[..., 1][cell_mask]
if np.any(fl_1 == 0):
div = np.nan
else:
div = np.median(fl_0 / fl_1)
else:
div = np.nan
return div
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
def centroid(cell_mask):
"""Find the cell's centroid."""
weights_c = np.arange(1, cell_mask.shape[1] + 1, 1).reshape(
1, cell_mask.shape[1]
)
weights_v = np.arange(1, cell_mask.shape[0] + 1, 1).reshape(
cell_mask.shape[0], 1
)
# moments
M00 = np.sum(cell_mask)
M10 = np.sum(np.multiply(cell_mask, weights_c))
M01 = np.sum(np.multiply(cell_mask, weights_v))
# centroid
Xm = M10 / M00
Ym = M01 / M00
return (Xm, Ym)
def centroid_x(cell_mask):
"""Return x coordinate of a cell's centroid."""
return centroid(cell_mask)[0]
def centroid_y(cell_mask):
"""Return y coordinate of a cell's centroid."""
return centroid(cell_mask)[1]