Retinal vessels: tiling, and a metric Dice cannot replace

Why a vessel cannot be resized

FIVES images are 2048 x 2048, and the obvious move is to resize them to something a U-Net is comfortable with. It is the wrong move here, and the reason is a measurement rather than a preference.

A FIVES vessel is about nine pixels wide - measured as mask area over skeleton length, over six random masks, on the engine side where the skeleton lives. Resized to 256 the same vessels are about 1.2 pixels. That matters because of what happens at one pixel: shift a one-pixel structure by a single pixel and it overlaps the truth nowhere at all, so Dice, IoU and clDice alike read zero and none of them can tell a near miss from a total failure. Tiling keeps the data out of that regime, which is what splits is for.

(The nine-pixel figure corrects the one-to-five this dataset is usually described with. One to five is what remains after a resize - a property of the preprocessing, not of the retina.)

The data

One 1.7 GB download from figshare 19688169, CC BY 4.0. The authors’ own split is already balanced across the four groups - 150 train and 50 test of each - so it is used rather than replaced. Point PLATYPUS_FIVES at the unpacked directory.

The dataset is laid out as Original/ and Ground truth/ with matching file names, which is not a layout platypus discovers on its own, so the pairs are written into a CSV:

root <- file.path(
  Sys.getenv("PLATYPUS_FIVES"),
  "FIVES A Fundus Image Dataset for AI-based Vessel Segmentation"
)
stopifnot(dir.exists(root))

write_split <- function(split, path) {
  # `pattern` rather than everything in the directory: the archive ships a `Thumbs.db`
  # beside the 600 training images, and a Windows thumbnail cache is not a fundus
  # photograph. Filtering by extension is also what keeps the pairing below honest.
  images <- sort(list.files(
    file.path(root, split, "Original"), pattern = "\\.png$", full.names = TRUE
  ))
  frame <- data.frame(
    # The key is what `evaluate_cases()` reports and what `group_by` reads. Without it a
    # case is called `row 55`, which names nothing on disk and cannot be grouped.
    key = tools::file_path_sans_ext(basename(images)),
    images = images,
    masks = file.path(root, split, "Ground truth", basename(images))
  )
  stopifnot(all(file.exists(frame$masks)))
  write.csv(frame, path, row.names = FALSE)
  nrow(frame)
}

splits <- file.path(tempdir(), c(train = "train.csv", test = "test.csv"))
names(splits) <- c("train", "test")
c(train = write_split("train", splits[["train"]]), test = write_split("test", splits[["test"]]))
#> train  test 
#>   600   200

The disease is in the file name - _N, _A, _D, _G for normal, AMD, diabetic retinopathy and glaucoma - which is what the grouping below reads.

Two things that cost time if taken on trust rather than checked: the archive ships a Thumbs.db beside the training images, so the listing is filtered by extension; and the masks are paired by name, which stopifnot above asserts rather than assumes.

The configuration

spec <- platypus_spec(
  task = "semantic_segmentation",
  data = segmentation_data(
    train = splits[["train"]],
    validation = splits[["test"]],
    colormap = list(c(0, 0, 0), c(255, 255, 255)),
    mode = "config_file"
  ),
  models = list(
    u_net(
      name = "vessels",
      input_shape = c(512, 512),
      # 2048 read whole and cut 4 x 4. This is the field the vignette exists for.
      splits = c(4, 4),
      blocks = 4,
      filters = 16,
      batch_size = 4,
      epochs = 15,
      loss = loss_cce_dice(),
      metrics = list(
        metric_dice(include_background = FALSE),
        metric_iou(include_background = FALSE),
        # Measured on the reassembled image by `evaluate_cases()`, not per tile.
        metric_cldice(include_background = FALSE)
      ),
      optimizer = optimizer_adam(learning_rate = 1e-3),
      augmentation = list(
        augmentation_step("HorizontalFlip", p = 0.5),
        augmentation_step("VerticalFlip", p = 0.5)
      )
    )
  ),
  seed = 1
)
spec
#> platypus specification
#>   data        : /tmp/RtmpLgvYTS/train.csv 
#>   classes     : 2 
#>   models      : 1 
#>     vessels        u_net            512x512, tiled 4x4, cce_dice, 15 epochs

Two things about that worth saying out loud.

metric_cldice() sits beside metric_dice() even though the run is tiled, and the engine reassembles each case before measuring it. It will not appear in the epoch log, and that is deliberate rather than an omission: an epoch averages its metrics over batches, which is exact for Dice - a ratio of sums - and meaningless for a skeleton, because a tile severs every vessel crossing its edge. So there is no val_cldice to watch, and a callback that tried would be refused when the specification is read.

Fifteen epochs rather than the forty this model wants, because a vignette that takes seven hours to rebuild is a vignette nobody rebuilds. The shape of the result is the same; the numbers are a little worse than the model can do.

Training

fit <- platypus_fit(spec)
fit
#> platypus fit
#>   models   : vessels 
#>   device   : NVIDIA GeForce GTX 1070 (torch 2.7.1+cu126) 
#>   epochs   : vessels=15 
#>   next     : evaluate(fit), predict(fit, "vessels")

plot(fit)

plot of chunk history

One number hides the finding

evaluate(fit)
#>     model architecture loss_function parameters  load_shape epochs_run
#> 1 vessels        u_net      cce_dice    1942594 2048 x 2048         15
#>   best_epoch       loss foreground      dice       iou dice_worst  iou_worst
#> 1         NA 0.07707521  0.0701509 0.8725505 0.7859755   0.100926 0.05315983

best_epoch is empty, and it is worth a word because the column is new in the engine’s 0.8.0a73 and reads as a gap otherwise. It names the epoch whose weights the run ended on, which only a callback that replaces weights can answer - callback_early_stopping() with its restore_best default, or one of the two averagers. The specification above asks for no callbacks at all, so nothing selected an epoch and this run ended on its last: epochs_run already says which.

That is the whole test set as one row, and it is the least informative view of it. The distribution is the finding, and summary() gives it - including which case is the worst, because the first thing anyone does with a low score is go and look at that image:

per_case <- evaluate_cases(fit)
summary(per_case)
#>   metric   n      mean        sd    median        min       max worst_case
#> 1   dice 200 0.8725505 0.1075551 0.9028940 0.10092598 0.9557962      122_G
#> 2    iou 200 0.7859755 0.1299171 0.8229785 0.05315983 0.9153351      122_G
#> 3 cldice 200 0.7892405 0.1445759 0.8286990 0.15440816 0.9353433      122_G

Now by the disease in the file name. group_by pools a group’s cases into one score, the way a patient’s slices would be pooled:

evaluate_cases(fit, group_by = "_([NADG])$")
#> platypus scores by group (4 rows)
#>  group      dice       iou    cldice
#>      A 0.9129180 0.8397875 0.8297471
#>      D 0.8900999 0.8019639 0.7770124
#>      G 0.8916201 0.8044354 0.7302487
#>      N 0.8843446 0.7926683 0.8199539

Pooling answers “how good is it on this group”. It does not show the spread, and the spread is where this dataset is interesting, so that comes off the per-case rows:

diseases <- c(N = "normal", A = "AMD", D = "diabetic retinopathy", G = "glaucoma")
per_case$disease <- diseases[sub(".*_([NADG])$", "\\1", per_case$case)]

spread <- function(x) c(n = length(x), mean = mean(x), sd = sd(x), worst = min(x))
round(t(vapply(split(per_case$dice, per_case$disease), spread, numeric(4))), 4)
#>                       n   mean     sd  worst
#> AMD                  50 0.9065 0.0515 0.7030
#> diabetic retinopathy 50 0.8753 0.0756 0.6137
#> glaucoma             50 0.8267 0.1803 0.1009
#> normal               50 0.8817 0.0520 0.6816

Compare that table with the pooled one above and the reason for having both becomes concrete. Pooled, glaucoma and diabetic retinopathy are within a few thousandths of each other. Per case, glaucoma’s mean is several points lower and its spread more than twice as wide, with a worst case far below its own group’s mean. (No figure is quoted here on purpose: the table is directly above, and a number repeated in prose is a number that goes stale the next time this is rebuilt.) Pooling hid it because Dice is a ratio of sums: fifty images’ counts added together are dominated by the ones with the most vessel to find, and a photograph the model failed on contributes little to the numerator and little to the denominator.

Before concluding that the model is worse at the disease whose row is worst, check whether it is the group whose photographs are worse. FIVES ships a Quality Assessment.xlsx grading every image for illumination, blur and contrast; nothing in this pipeline reads it, and when it is joined to these scores the ordering above changes substantially. The engine’s measurement of that is in measurements/fives/quality.py in the project’s workbench.

That is the general lesson rather than a fact about retinas: a per-group table is worth having precisely because it invites the question why is this group different, and the first answer is rarely the right one.

clDice beside Dice

round(c(
  correlation = cor(per_case$dice, per_case$cldice),
  mean_difference = mean(per_case$cldice - per_case$dice),
  sd_of_difference = sd(per_case$cldice - per_case$dice)
), 4)
#>      correlation  mean_difference sd_of_difference 
#>           0.8581          -0.0833           0.0760

Read those two numbers together, because they say different things. The correlation is high: the two metrics rank the cases much the same way, so clDice is not finding a different set of bad images. The mean difference is negative and not small - clDice sits several points below Dice - which says the vessels this model draws are less connected than their pixel overlap suggests. It is getting the right pixels and breaking the lines.

That gap is the metric’s peel count rather than the model’s training, and the way to tell is to measure one set of weights two ways. metric_cldice() above takes its default iterations = 5; the engine’s article measures the same metric at iterations: 12, a number it argues for against vessels nine pixels across. Those are two measurements rather than one, so reading this page’s clDice against the article’s says something about the metric before it says anything about either model.

The published forty-epoch weights settle how much is which, because they can be re-scored without being retrained. On their own test split, at this page’s own five peels, they give a correlation of 0.856 and a mean difference of -0.082 - this page’s own 0.8581 and -0.0833, to within 0.002 on both. At twelve peels the same weights give 0.973 and -0.001.

So of the distance between this page’s figures and the article’s, about 0.085 is the peel count and 0.0015 is the forty epochs. A large negative gap here is a skeleton that was not thinned, rather than a model breaking lines. Twelve peels is chosen so that the widest vessels in the set reach a centreline and not only the average one; at five they do not, and what is left of them is counted as skeleton - which makes part of this number a measure of thickness, the one thing clDice exists to ignore.

The article’s own table was measured before clDice’s definition changed, and the two re-scores quoted above were not. 0.8.0a72 made the soft skeleton erode with a four-neighbour cross rather than a 3 x 3 square, because that is what the reference implementation does; a cross erodes less, leaves a thinner skeleton, and scores lower. How much lower depends on the prediction - measured on ten of these masks with the prediction thinned, between 0.004 and 0.074, and on the published weights at twelve peels, 0.004 - so a clDice from one definition and a clDice from the other are not a difference worth reading. Everything this document prints, and both re-scores above, come from the pinned engine. The article’s forty-epoch table and every clDice figure in the tables below do not.

The worst case, drawn

A table says one case scored far below the rest. A picture says why — and this is the panel that found the quality finding above, because the worst glaucoma photograph turns out to be nearly black.

Drawing it needs two images out of two hundred, and the whole split at 2048 x 2048 is about 6.7 GB, so predict() is asked for the two by name. It keeps one answer at a time and returns the pair in the order asked for:

chosen <- per_case$case[c(which.max(per_case$dice), which.min(per_case$dice))]
masks <- predict(fit, split = "validation", cases = chosen)

at_source <- function(what) file.path(root, "test", what, paste0(chosen, ".png"))

plot_masks(
  images = read_images(at_source("Original"), size = NULL, channels = 3),
  # One `height x width` matrix per case, stacked the way plot_masks wants them:
  # image x height x width.
  prediction = aperm(simplify2array(masks), c(3, 1, 2)),
  truth = read_masks(at_source("Ground truth"), colormap = binary_colormap, size = NULL),
  labels = sprintf("%s  dice %.3f", chosen, per_case$dice[match(chosen, per_case$case)])
)

plot of chunk worst

Two of those choices, measured afterwards

loss_cce_dice() and splits = c(4, 4) are the two decisions in the specification above that read as taste. Both have since been measured on this dataset - after this vignette was written, with the engine, and with every arm kept rather than described - and one of the answers contradicts the section this document opens with.

The tables below are transcriptions. The runs live in the project’s workbench under measurements/fives/, and each table is printed by a script that reads the stored records back: losses-read.py, boundary-read.py and tiling-read.py. Nothing in this document runs them, and no number below was produced by the knit you are reading.

Every arm below was measured between 2026-10-16 and 2026-10-22, which the clDice - dice column makes matter. clDice’s definition moved in 0.8.0a72, after all of them had run, and the arms’ checkpoints were not kept - so that column cannot be re-read under the current definition without retraining fifteen models, and it is left as it was measured rather than adjusted. Within each table the comparison is sound, because every arm in it was scored the same way; it is reading one of these columns against this document’s own clDice that is not. The Dice columns are unaffected.

The setting is not this vignette’s, which is the first thing to know before comparing anything to the numbers above. Every arm trained on 148 of the 600 training retinas - every second one within each disease group, so the four stay balanced - for up to 40 epochs with early stopping on val_dice at patience 8, and every arm was scored case by case on the same 200 test images. So no absolute Dice below is comparable to the one this document printed; what the tables are for is the differences, each paired case against case.

A difference is read against a measured floor rather than against zero. The baseline ran twice with different seeds, and the floor is the larger of that pair’s mean difference and that mean’s own standard error: 0.0009 over all 200 cases and 0.0022 over the 67 that FIVES’s own quality assessment flags - the grading the section above says nothing in this pipeline reads. The x columns are multiples of that floor, and a row at 2x says nothing.

Notes

The card. This was precomputed on a GTX 1070, which is a Pascal-generation GPU and needs the last PyTorch built against CUDA 12. platypus_use_torch("pascal") asks for it, and it has to be called before anything starts the engine. Without it this machine trains on the processor and says so. platypus_device() reports what yours can see.

The cost. Fifteen epochs over 600 retinas at 2048, cut sixteen ways, takes a few hours on that card; scoring the 200 test images - including reassembling each one for clDice - takes about two minutes. Tiling is what makes the first number bearable: a tiled run reuses each decoded image for all of its tiles, which is worth roughly a factor of two.

Asking for a few cases is cheap in memory, not in time. cases = never holds more than one answer, which is what makes the figure above possible at all - but the engine’s stream is sequential, so two cases out of two hundred still run the model over two hundred. Pass a function to each = instead when you want to walk the split and keep a number rather than a mask from each image.