"""
Copyright © 2023 Howard Hughes Medical Institute, Authored by Carsen Stringer and Marius Pachitariu.
"""
import numpy as np
from . import utils, dynamics
from numba import jit
from scipy.optimize import linear_sum_assignment
from scipy.ndimage import convolve, mean
[docs]def mask_ious(masks_true, masks_pred):
"""Return best-matched masks."""
iou = _intersection_over_union(masks_true, masks_pred)[1:, 1:]
n_min = min(iou.shape[0], iou.shape[1])
costs = -(iou >= 0.5).astype(float) - iou / (2 * n_min)
true_ind, pred_ind = linear_sum_assignment(costs)
iout = np.zeros(masks_true.max())
iout[true_ind] = iou[true_ind, pred_ind]
preds = np.zeros(masks_true.max(), "int")
preds[true_ind] = pred_ind + 1
return iout, preds
[docs]def boundary_scores(masks_true, masks_pred, scales):
"""
Calculate boundary precision, recall, and F-score.
Args:
masks_true (list): List of true masks.
masks_pred (list): List of predicted masks.
scales (list): List of scales.
Returns:
tuple: A tuple containing precision, recall, and F-score arrays.
"""
diams = [utils.diameters(lbl)[0] for lbl in masks_true]
precision = np.zeros((len(scales), len(masks_true)))
recall = np.zeros((len(scales), len(masks_true)))
fscore = np.zeros((len(scales), len(masks_true)))
for j, scale in enumerate(scales):
for n in range(len(masks_true)):
diam = max(1, scale * diams[n])
rs, ys, xs = utils.circleMask([int(np.ceil(diam)), int(np.ceil(diam))])
filt = (rs <= diam).astype(np.float32)
otrue = utils.masks_to_outlines(masks_true[n])
otrue = convolve(otrue, filt)
opred = utils.masks_to_outlines(masks_pred[n])
opred = convolve(opred, filt)
tp = np.logical_and(otrue == 1, opred == 1).sum()
fp = np.logical_and(otrue == 0, opred == 1).sum()
fn = np.logical_and(otrue == 1, opred == 0).sum()
precision[j, n] = tp / (tp + fp)
recall[j, n] = tp / (tp + fn)
fscore[j] = 2 * precision[j] * recall[j] / (precision[j] + recall[j])
return precision, recall, fscore
[docs]def aggregated_jaccard_index(masks_true, masks_pred):
"""
AJI = intersection of all matched masks / union of all masks
Args:
masks_true (list of np.ndarrays (int) or np.ndarray (int)):
where 0=NO masks; 1,2... are mask labels
masks_pred (list of np.ndarrays (int) or np.ndarray (int)):
np.ndarray (int) where 0=NO masks; 1,2... are mask labels
Returns:
aji (float): aggregated jaccard index for each set of masks
"""
aji = np.zeros(len(masks_true))
for n in range(len(masks_true)):
iout, preds = mask_ious(masks_true[n], masks_pred[n])
inds = np.arange(0, masks_true[n].max(), 1, int)
overlap = _label_overlap(masks_true[n], masks_pred[n])
union = np.logical_or(masks_true[n] > 0, masks_pred[n] > 0).sum()
overlap = overlap[inds[preds > 0] + 1, preds[preds > 0].astype(int)]
aji[n] = overlap.sum() / union
return aji
[docs]def average_precision(masks_true, masks_pred, threshold=[0.5, 0.75, 0.9]):
"""
Average precision estimation: AP = TP / (TP + FP + FN)
This function is based heavily on the *fast* stardist matching functions
(https://github.com/mpicbg-csbd/stardist/blob/master/stardist/matching.py)
Args:
masks_true (list of np.ndarrays (int) or np.ndarray (int)):
where 0=NO masks; 1,2... are mask labels
masks_pred (list of np.ndarrays (int) or np.ndarray (int)):
np.ndarray (int) where 0=NO masks; 1,2... are mask labels
Returns:
ap (array [len(masks_true) x len(threshold)]):
average precision at thresholds
tp (array [len(masks_true) x len(threshold)]):
number of true positives at thresholds
fp (array [len(masks_true) x len(threshold)]):
number of false positives at thresholds
fn (array [len(masks_true) x len(threshold)]):
number of false negatives at thresholds
"""
not_list = False
if not isinstance(masks_true, list):
masks_true = [masks_true]
masks_pred = [masks_pred]
not_list = True
if not isinstance(threshold, list) and not isinstance(threshold, np.ndarray):
threshold = [threshold]
if len(masks_true) != len(masks_pred):
raise ValueError(
"metrics.average_precision requires len(masks_true)==len(masks_pred)")
ap = np.zeros((len(masks_true), len(threshold)), np.float32)
tp = np.zeros((len(masks_true), len(threshold)), np.float32)
fp = np.zeros((len(masks_true), len(threshold)), np.float32)
fn = np.zeros((len(masks_true), len(threshold)), np.float32)
n_true = np.array(list(map(np.max, masks_true)))
n_pred = np.array(list(map(np.max, masks_pred)))
for n in range(len(masks_true)):
#_,mt = np.reshape(np.unique(masks_true[n], return_index=True), masks_pred[n].shape)
if n_pred[n] > 0:
iou = _intersection_over_union(masks_true[n], masks_pred[n])[1:, 1:]
for k, th in enumerate(threshold):
tp[n, k] = _true_positive(iou, th)
fp[n] = n_pred[n] - tp[n]
fn[n] = n_true[n] - tp[n]
ap[n] = tp[n] / (tp[n] + fp[n] + fn[n])
if not_list:
ap, tp, fp, fn = ap[0], tp[0], fp[0], fn[0]
return ap, tp, fp, fn
@jit(nopython=True)
def _label_overlap(x, y):
"""Fast function to get pixel overlaps between masks in x and y.
Args:
x (np.ndarray, int): Where 0=NO masks; 1,2... are mask labels.
y (np.ndarray, int): Where 0=NO masks; 1,2... are mask labels.
Returns:
overlap (np.ndarray, int): Matrix of pixel overlaps of size [x.max()+1, y.max()+1].
"""
# put label arrays into standard form then flatten them
# x = (utils.format_labels(x)).ravel()
# y = (utils.format_labels(y)).ravel()
x = x.ravel()
y = y.ravel()
# preallocate a "contact map" matrix
overlap = np.zeros((1 + x.max(), 1 + y.max()), dtype=np.uint)
# loop over the labels in x and add to the corresponding
# overlap entry. If label A in x and label B in y share P
# pixels, then the resulting overlap is P
# len(x)=len(y), the number of pixels in the whole image
for i in range(len(x)):
overlap[x[i], y[i]] += 1
return overlap
def _intersection_over_union(masks_true, masks_pred):
"""Calculate the intersection over union of all mask pairs.
Parameters:
masks_true (np.ndarray, int): Ground truth masks, where 0=NO masks; 1,2... are mask labels.
masks_pred (np.ndarray, int): Predicted masks, where 0=NO masks; 1,2... are mask labels.
Returns:
iou (np.ndarray, float): Matrix of IOU pairs of size [x.max()+1, y.max()+1].
How it works:
The overlap matrix is a lookup table of the area of intersection
between each set of labels (true and predicted). The true labels
are taken to be along axis 0, and the predicted labels are taken
to be along axis 1. The sum of the overlaps along axis 0 is thus
an array giving the total overlap of the true labels with each of
the predicted labels, and likewise the sum over axis 1 is the
total overlap of the predicted labels with each of the true labels.
Because the label 0 (background) is included, this sum is guaranteed
to reconstruct the total area of each label. Adding this row and
column vectors gives a 2D array with the areas of every label pair
added together. This is equivalent to the union of the label areas
except for the duplicated overlap area, so the overlap matrix is
subtracted to find the union matrix.
"""
overlap = _label_overlap(masks_true, masks_pred)
n_pixels_pred = np.sum(overlap, axis=0, keepdims=True)
n_pixels_true = np.sum(overlap, axis=1, keepdims=True)
iou = overlap / (n_pixels_pred + n_pixels_true - overlap)
iou[np.isnan(iou)] = 0.0
return iou
def _true_positive(iou, th):
"""Calculate the true positive at threshold th.
Args:
iou (float, np.ndarray): Array of IOU pairs.
th (float): Threshold on IOU for positive label.
Returns:
tp (float): Number of true positives at threshold.
How it works:
(1) Find minimum number of masks.
(2) Define cost matrix; for a given threshold, each element is negative
the higher the IoU is (perfect IoU is 1, worst is 0). The second term
gets more negative with higher IoU, but less negative with greater
n_min (but that's a constant...).
(3) Solve the linear sum assignment problem. The costs array defines the cost
of matching a true label with a predicted label, so the problem is to
find the set of pairings that minimizes this cost. The scipy.optimize
function gives the ordered lists of corresponding true and predicted labels.
(4) Extract the IoUs from these pairings and then threshold to get a boolean array
whose sum is the number of true positives that is returned.
"""
n_min = min(iou.shape[0], iou.shape[1])
costs = -(iou >= th).astype(float) - iou / (2 * n_min)
true_ind, pred_ind = linear_sum_assignment(costs)
match_ok = iou[true_ind, pred_ind] >= th
tp = match_ok.sum()
return tp
[docs]def flow_error(maski, dP_net, device=None):
"""Error in flows from predicted masks vs flows predicted by network run on image.
This function serves to benchmark the quality of masks. It works as follows:
1. The predicted masks are used to create a flow diagram.
2. The mask-flows are compared to the flows that the network predicted.
If there is a discrepancy between the flows, it suggests that the mask is incorrect.
Masks with flow_errors greater than 0.4 are discarded by default. This setting can be
changed in Cellpose.eval or CellposeModel.eval.
Args:
maski (np.ndarray, int): Masks produced from running dynamics on dP_net, where 0=NO masks; 1,2... are mask labels.
dP_net (np.ndarray, float): ND flows where dP_net.shape[1:] = maski.shape.
Returns:
flow_errors (np.ndarray, float): Mean squared error between predicted flows and flows from masks.
dP_masks (np.ndarray, float): ND flows produced from the predicted masks.
"""
if dP_net.shape[1:] != maski.shape:
print("ERROR: net flow is not same size as predicted masks")
return
# flows predicted from estimated masks
dP_masks = dynamics.masks_to_flows(maski, device=device)
# difference between predicted flows vs mask flows
flow_errors = np.zeros(maski.max())
for i in range(dP_masks.shape[0]):
flow_errors += mean((dP_masks[i] - dP_net[i] / 5.)**2, maski,
index=np.arange(1,
maski.max() + 1))
return flow_errors, dP_masks