Nuclei, five years later: what a mean Dice hides

In 2020 I trained a U-Net on the 2018 Data Science Bowl nuclei with the first version of this package, watched a Dice score climb during training, and finished the post with four predicted masks on the competition’s test images. Those images have no masks. The pictures looked good, and I had no way of knowing where the model was wrong – so I never asked.

This post asks. Nothing in it trains: the model is the published dsbowl-unet from the first post of the series, scored on the images it was held out from. Scoring all of them took about twelve seconds on a processor.

The images it was held out from

The data is BBBC038v1 from the Broad Institute – the Data Science Bowl’s images under a CC0 waiver. Its stage1_train holds 670 annotated images, one directory each, with one mask file per nucleus.

data_dir <- "stage1_train"   # BBBC038v1, unzipped

data_dir = "stage1_train"   # BBBC038v1, unzipped

The weights were measured on an 80/20 split of all 670 at seed 1, so that split can be made again. split_dataset() copies nothing; it writes CSVs listing the files:

library(rplatypus)

split <- split_dataset(data_dir, tempfile("splits"),
                       fractions = c(0.8, 0.2), seed = 1)
split
#> platypus split
#>   train        536 samples
#>   validation   134 samples
#>   every sample its own group
#>   seed 1 - the same inputs give this split again

import tempfile
import pyplatypus as pp

split = pp.split_dataset(data_dir, tempfile.mkdtemp(),
                         fractions=(0.8, 0.2), seed=1)
split["samples"]
#> {'train': 536, 'validation': 134, 'test': 0}

Making the split again matters more than it looks. My own copy of this dataset has a directory called stage1_validation with exactly 134 images in it – the same count – and 101 of them are in the training part of the split the weights were measured on. Scored there, the model would be graded partly on images it learnt from, and the count would never tell you.

One number

spec <- platypus_spec(
  data = segmentation_data(split, colormap = binary_colormap),
  models = list(
    u_net("nuclei", input_shape = c(256, 256), blocks = 4, filters = 16,
          weights = "dsbowl-unet", fit = FALSE,
          metrics = list(metric_dice(include_background = FALSE)))
  )
)
fit <- platypus_fit(spec, verbose = FALSE)

cases <- evaluate_cases(fit)
summary(cases)
#>   metric   n      mean         sd    median       min       max
#> 1   dice 134 0.9205236 0.05322748 0.9313503 0.7239583 0.9860661
#>                                                         worst_case
#> 1 7f38885521586fc6011bef1314a9fb2aa1e4935bd581b2991e1d963395eab770

engine = pp.build_engine(pp.from_dict({
    "task": "semantic_segmentation",
    "data": {
        "train_path": str(split["paths"]["train"]),
        "validation_path": str(split["paths"]["validation"]),
        "colormap": [[0, 0, 0], [255, 255, 255]],
        "mode": "config_file",
    },
    "models": [{
        "name": "nuclei",
        "architecture": "u_net",
        "input_shape": [256, 256],
        "blocks": 4,
        "filters": 16,
        "weights": "dsbowl-unet",
        "fit": False,
        "metrics": [{"name": "dice", "include_background": False}],
    }],
}))
engine.fit()

cases = engine.evaluate_cases("nuclei")
pp.summarise_cases(cases)
#> [{'metric': 'dice', 'n': 134, 'mean': 0.9205236150257623,
#>   'sd': 0.05322747958922138, 'median': 0.9313502907752991,
#>   'min': 0.7239583134651184, 'max': 0.9860661029815674,
#>   'worst_case':
#>     '7f38885521586fc6011bef1314a9fb2aa1e4935bd581b2991e1d963395eab770'}]

A mean Dice of 0.9205 and a worst image at 0.7240 – the model card’s numbers to four digits, which is the check that this is the card’s split. include_background = FALSE keeps the background out of the average; with it in, a model that found no nuclei at all would still score around half.

0.92 is a good number, and it is true. The question is what it is an average of.

One row per image

evaluate() gives one row per model. evaluate_cases() gives one row per image, and that is the table worth reading:

head(cases[order(cases$dice), ], 5)
#> platypus scores by case (5 rows)
#>                                                              case      dice
#>  7f38885521586fc6011bef1314a9fb2aa1e4935bd581b2991e1d963395eab770 0.7239583
#>  29ea4f6eb4545f43868a9b40a60000426bf8dfd9d062546656a37bd2a2aaf9ec 0.7531800
#>  8f94a80b95a881d0efdec36affc915dca9609f4cba8134c4a91b219d418778aa 0.7541562
#>  a6515d73077866808ad4cb837ecdac33612527b8a1041e82135e40fce2bb9380 0.7588775
#>  2a1a294e21d76efd0399e4eb321b45f44f7510911acd92c988480195c5b4c812 0.7702851

worst = sorted(cases, key=lambda row: row["dice"])[:5]
[(row["case"][:16], round(row["dice"], 4)) for row in worst]
#> [('7f38885521586fc6', 0.724), ('29ea4f6eb4545f43', 0.7532),
#>  ('8f94a80b95a881d0', 0.7542), ('a6515d7307786680', 0.7589),
#>  ('2a1a294e21d76efd', 0.7703)]

Those names say nothing, so give each image a label. The Data Science Bowl mixed three kinds of microscopy, and two numbers tell them apart well enough: whether the picture has colour, and whether its background is light or dark.

image_kind <- function(case) {
  path <- list.files(file.path(data_dir, case, "images"), full.names = TRUE)
  pixels <- png::readPNG(path[1])[, , 1:3]
  colour <- mean(abs(pixels[, , 1] - pixels[, , 2])) +
    mean(abs(pixels[, , 2] - pixels[, , 3]))
  if (colour > 0.02) return("stained")
  if (mean(pixels) > 0.5) "brightfield" else "fluorescence"
}
cases$kind <- vapply(cases$case, image_kind, character(1))

aggregate(dice ~ kind, cases, function(x) {
  round(c(n = length(x), mean = mean(x), min = min(x)), 3)
})
#>           kind  dice.n dice.mean dice.min
#> 1  brightfield   3.000     0.749    0.724
#> 2 fluorescence 105.000     0.935    0.753
#> 3      stained  26.000     0.883    0.805

from pathlib import Path
import numpy as np
from PIL import Image

def image_kind(case):
    path = next((Path(data_dir) / case / "images").glob("*.png"))
    pixels = np.asarray(Image.open(path).convert("RGB")) / 255
    colour = (np.abs(pixels[..., 0] - pixels[..., 1]).mean()
              + np.abs(pixels[..., 1] - pixels[..., 2]).mean())
    if colour > 0.02:
        return "stained"
    return "brightfield" if pixels.mean() > 0.5 else "fluorescence"

for row in cases:
    row["kind"] = image_kind(row["case"])

kinds = sorted({row["kind"] for row in cases})
dice = {k: [r["dice"] for r in cases if r["kind"] == k] for k in kinds}
{k: (len(d), round(float(np.mean(d)), 3), round(min(d), 3))
 for k, d in dice.items()}
#> {'brightfield': (3, 0.749, 0.724), 'fluorescence': (105, 0.935, 0.753),
#>  'stained': (26, 0.883, 0.805)}

library(ggplot2)
ggplot(cases, aes(dice, kind)) +
  geom_vline(xintercept = mean(cases$dice), linetype = "dashed",
             colour = "grey50") +
  geom_jitter(height = 0.15, width = 0, alpha = 0.7, size = 2) +
  labs(x = "Dice, per image", y = NULL) +
  theme_minimal(base_size = 13)

import matplotlib.pyplot as plt

rng = np.random.default_rng(1)
fig, ax = plt.subplots(figsize=(7, 3.2))
ax.axvline(np.mean([r["dice"] for r in cases]), linestyle="--", color="grey")
for i, k in enumerate(kinds):
    jitter = rng.uniform(-0.15, 0.15, len(dice[k]))
    ax.scatter(dice[k], i + jitter, alpha=0.7)
ax.set_yticks(range(len(kinds)), kinds)
ax.set_xlabel("Dice, per image")

There are three brightfield images among the 134, and they score 0.72, 0.75 and 0.77 – three of the five worst. Fluorescence averages 0.93 over 105 images, and it is most of what the mean measures. Stained tissue sits between the two, at 0.88.

predictions <- predict(fit, split = "validation")

files_of <- function(case, what) {
  list.files(file.path(data_dir, case, what), full.names = TRUE)
}
truth_of <- function(case) {
  masks <- lapply(files_of(case, "masks"), function(f) {
    read_masks(f, binary_colormap, size = c(256, 256))[1, , ]
  })
  unite_masks(masks)
}

middle <- which.min(abs(cases$dice - median(cases$dice)))
shown <- c(which(cases$kind == "brightfield"), middle)
images <- read_images(
  vapply(cases$case[shown], function(k) files_of(k, "images")[1], ""),
  size = c(256, 256)
)
truth <- simplify2array(lapply(cases$case[shown], truth_of))

plot_masks(
  images,
  prediction = predictions[shown, , ],
  truth      = aperm(truth, c(3, 1, 2)),
  colormap   = binary_colormap,
  labels     = sprintf("%s, Dice %.2f", cases$kind[shown], cases$dice[shown])
)

predictions = engine.predict("nuclei", "validation")
data = engine.dataset(engine.runs["nuclei"].spec, "validation")

middle = np.median([r["dice"] for r in cases])
shown = [i for i, r in enumerate(cases) if r["kind"] == "brightfield"]
shown.append(min(range(len(cases)),
                 key=lambda i: abs(cases[i]["dice"] - middle)))

pp.plot_masks(
    np.stack([data[i][0] for i in shown]),
    prediction=predictions[shown],
    truth=np.stack([data[i][1] for i in shown]),
    labels=[f"{cases[i]['kind']}, Dice {cases[i]['dice']:.2f}"
            for i in shown],
)

The kind of microscope, or the size of a nucleus?

The pictures suggest a second reading. The brightfield images are 1024 by 1024 pixels, and the model sees every image at 256 by 256 – so a nucleus sixteen pixels across arrives as four. Measure the average nucleus on the model’s grid:

nucleus_area <- function(case) {
  files <- list.files(file.path(data_dir, case, "masks"), full.names = TRUE)
  areas <- vapply(files, function(f) {
    sum(read_masks(f, binary_colormap) == 2)
  }, numeric(1))
  size <- dim(read_masks(files[1], binary_colormap))[2:3]
  mean(areas) * 256^2 / prod(size)
}
cases$area <- vapply(cases$case, nucleus_area, numeric(1))

aggregate(area ~ kind, cases, function(x) round(range(x)))
#>           kind area.1 area.2
#> 1  brightfield     11     17
#> 2 fluorescence     69    957
#> 3      stained    137   1893

def nucleus_area(case):
    files = sorted((Path(data_dir) / case / "masks").glob("*.png"))
    masks = [np.asarray(Image.open(f)) > 0 for f in files]
    height, width = masks[0].shape[:2]
    return np.mean([m.sum() for m in masks]) * 256**2 / (height * width)

for row in cases:
    row["area"] = nucleus_area(row["case"])

areas = {k: [r["area"] for r in cases if r["kind"] == k] for k in kinds}
{k: (round(min(a)), round(max(a))) for k, a in areas.items()}
#> {'brightfield': (11, 17), 'fluorescence': (69, 957),
#>  'stained': (137, 1893)}

On every other image an average nucleus covers at least 69 pixels of the grid. On the three brightfield images it covers eleven to seventeen. So the three images are both the only brightfield in the split and the only ones where the nuclei are this small, and three images cannot tell those two explanations apart. Training does not help either: of the 536 training images, 13 are brightfield, all of them 1024 by 1024 as well.

Scoring at the files’ own resolution does not separate them either, because the detail was lost before the score was taken – the prediction was made at 256:

summary(evaluate_cases(fit, at = "source"))
#>   metric   n     mean         sd    median       min       max
#> 1   dice 134 0.920808 0.05371718 0.9313503 0.7164674 0.9860661
#>                                                         worst_case
#> 1 7f38885521586fc6011bef1314a9fb2aa1e4935bd581b2991e1d963395eab770

pp.summarise_cases(engine.evaluate_cases("nuclei", at="source"))
#> [{'metric': 'dice', 'n': 134, 'mean': 0.920808045276955,
#>   'sd': 0.053717183966230486, 'median': 0.9313502907752991,
#>   'min': 0.7164673805236816, 'max': 0.9860661029815674,
#>   'worst_case':
#>     '7f38885521586fc6011bef1314a9fb2aa1e4935bd581b2991e1d963395eab770'}]

What would separate them is a model that sees these images at a resolution where a nucleus is as large as it is everywhere else – tiling instead of resizing, which is a post of its own later in the series.

What to take from this

The mean was not wrong. 0.92 is what this model scores on these 134 images. But it is mostly an average over 105 fluorescence images, and the three images it does worst on share two properties at once, either of which could be the cause. None of that is visible in one number; all of it is visible in a table with one row per image and a column you add yourself.

If you report a segmentation score, report the worst case beside the mean, and look at what the worst cases have in common. Here it took one function call and a few lines of R or Python.

The series so far

1. rplatypus rebuilt: deep learning for medical images in R – why the 2020 package had to be rewritten, how it now installs its own Python, and a first mask from the same published weights used here.
2. Nuclei, five years later – this post.

The next one stays with segmentation but steps back from a single dataset: what rplatypus offers for it – four architectures and how they differ, the losses and metrics to choose from, and the YAML file that describes a whole experiment, field by field.