Segmenting volumes: CT, spacing and per-patient scores

Synthetic spheres with known volumes, a split that keeps patients whole, and the question Dice cannot answer.

The other articles work on two-dimensional pictures. This one is about volumes, which is where medical imaging mostly lives: a CT or an MRI is a stack of slices, and three things that do not exist in 2D decide whether a result means anything.

Voxels have a physical size, and it differs between scans. Slices belong to a patient, so a split that ignores who they came from measures memory rather than generalisation. And a mask has a position in space, without which it cannot be laid over the scan it describes.

The data here is synthetic, generated by the code below, so every line of this page runs without downloading anything or registering anywhere. That is the reason it is generated inline rather than by a script in examples/: everything else in that directory is a program you point at real data, and a page whose data arrives with it is the one page here that needs nothing.

Where the numbers come from. One run of exactly the code below, on this machine: 24 scans, 12 epochs, and every table transcribed from its output. Nothing is determined - this specification does not ask for deterministic - so a second run moves the per-case scores by about 0.03 of Dice. The shape of every finding survives that; the third decimal does not.

This page mirrors vignette("volumes") in the R package, which is the source for the prose. The walk is the same and the route is not.

Making something to segment

Twenty-four scans from twelve patients, two each. Each holds one spherical “lesion” in air, and the scans differ in three ways that are ordinary in clinical data: slice thickness, half at 1 mm and half at 2.5 mm, so the same anatomy occupies a different number of voxels; lesion size, 4 to 12 mm in radius; and noise.

The variety is the point. A synthetic dataset where every case is equally easy would make the per-patient section below look like a formality, and it is not one.

import math
from pathlib import Path

import nibabel as nib
import numpy as np

rng = np.random.default_rng(1)

def make_scan(root: Path, patient: int, study: int, slice_mm: float, radius_mm: float, noise: float) -> str:
    shape, spacing = (48, 48, 32), (1.0, 1.0, slice_mm)
    grid = np.stack(np.meshgrid(*[np.arange(s) for s in shape], indexing="ij"), axis=-1)
    # Millimetres from the centre, per voxel - which is what makes the sphere a sphere in
    # the patient rather than in the array.
    offset = (grid - (np.array(shape) - 1) / 2) * np.array(spacing)
    # Move the lesion per patient, so the model cannot learn one position and be done.
    offset[..., 0] -= (patient - 6) * 1.5
    distance = np.sqrt((offset ** 2).sum(axis=-1))

    lesion = (distance <= radius_mm).astype(np.int16)
    # Hounsfield units: air is -1000, soft tissue about 40, plus scanner noise.
    scan = np.where(lesion == 1, 40.0, -1000.0) + rng.normal(0, noise, shape)

    name = f"patient{patient:02d}_study{study}"
    for folder in ("images", "masks"):
        (root / name / folder).mkdir(parents=True, exist_ok=True)
    affine = np.diag([*spacing, 1.0])
    nib.save(nib.Nifti1Image(scan.astype(np.float32), affine), str(root / name / "images" / "ct.nii.gz"))
    nib.save(nib.Nifti1Image(lesion, affine), str(root / name / "masks" / "lesion.nii.gz"))
    return name

Keeping the radius is the part worth copying into a real study where it is possible: the volume of a sphere is known exactly, which makes it possible to ask later whether the model got the size right and not only the overlap. That question is the last third of this page and no overlap metric answers it.

Two scans, the same number of slices, and one covers 32 mm of patient where the other covers 80 mm. volume_spacing() and volume_shape() read that out of the files themselves, through the same reader training uses:

from pyplatypus.data.volumes import volume_shape, volume_spacing

volume_spacing("scans/patient01_study1/images/ct.nii.gz")   # (1.0, 1.0, 2.5)
volume_shape("scans/patient01_study1/images/ct.nii.gz")     # (48, 48, 32)

A voxel in the first is two and a half times the volume of a voxel in the second. Nothing in the arrays says so, which is the whole reason this pair of functions exists.

Splitting by patient, not by scan

Each of these patients has two studies. Split the twenty-four scans at random and both studies of some patient land on opposite sides: the model meets that patient in training and is then tested on the same anatomy.

How much that costs is a separate question, and it has now been measured rather than asserted. This paragraph used to end “Dice comes out several points too high”, which nothing had ever checked. On MSD Spleen - 41 patients, 1043 slices - a slice-level split put 29 of 29 development patients on both sides, and the model it produced overstated its own score by +0.035 against patients it had never seen. A patient-level split on the same data overstated itself by +0.039. So the inflation is real and it is not the leak: any validation set differs from any test set by a few points because patients differ, and on this data the leak added nothing on top of that.

That is one dataset, one architecture and twelve epochs, which is not evidence that the mistake is harmless where a model has room to memorise. The reason to keep patients whole is that the mechanism is certain and invisible in the result, not that the size is known.

split_dataset() takes a pattern saying which part of a name is the patient, and keeps every patient whole:

import pyplatypus as pp

split = pp.split_dataset(
    "scans", "splits",
    group_by=r"^(patient\d+)_",
    fractions=(0.5, 0.25, 0.25),
    seed=1,
)

Twelve scans to train on, six to validate, six to test - six patients, three and three, with no patient on two sides. A pattern that matches nothing is an error rather than a fallback, which is the point: a silent fallback to one group per scan would reintroduce exactly the bug the argument exists to prevent.

The specification

Three volume-specific fields carry the weight. labels says the mask is a label map rather than a picture - which is how every volume format stores masks, and the reason colormap is not the only way to name classes. window maps Hounsfield units to 0-1 through a fixed window, so the same tissue is the same brightness in every scan. And target_spacing resamples everything to a common voxel size before cropping to the model’s shape.

task: semantic_segmentation
data:
  train_path: splits/train.csv
  validation_path: splits/validation.csv
  mode: config_file
  labels: [0, 1]
  window: soft_tissue
  target_spacing: [1, 1, 1]
models:
  - name: lesion
    architecture: u_net
    # Three entries, and that is the only thing making this a 3D run.
    input_shape: [32, 32, 32]
    channels: 1
    blocks: 3
    filters: 16
    epochs: 12
    batch_size: 2
    loss:
      name: dice
    metrics:
      - name: dice
        include_background: false
    augmentation:
      - name: HorizontalFlip
        params:
          p: 0.5
      - name: CubicSymmetry
        params:
          p: 0.5
seed: 1

Nothing else in the specification is different from a 2D run - the same losses, the same metrics, the same tiling if you want it.

window: soft_tissue is worth a sentence of its own. Scaling each scan by its own extremes is the obvious alternative and it throws away the thing Hounsfield units are for: air is -1000 on every scanner ever built, so a fixed window maps the same tissue to the same number in every scan, and one bright voxel - an implant, a marker - cannot rescale everything else.

Not every augmentation works on a volume. albumentations supports them unevenly, so ask rather than guess:

len(pp.available_transforms())           # 118 here
len(pp.available_transforms(rank=3))     # 91 here

Those two counts are a property of the installed albumentations rather than of this package, which is why the difference is not quoted as a fact about pyplatypus. A transform from the difference is refused by name when the run starts, rather than raising something from inside the library an hour in.

Training

spec = pp.from_yaml("volumes.yaml")
engine = pp.build_engine(spec)
history = engine.fit()["lesion"]

1,401,282 parameters at blocks: 3 and filters: 16, twelve epochs. A 3D U-Net grows quickly: the same blocks and filters in 2D is a fraction of that, because every convolution kernel gains a dimension.

Scores per patient, not per dataset

evaluate() gives one row per model - one number for the whole validation set:

model	architecture	loss_function	parameters	load_shape	epochs_run	loss	dice
lesion	u_net	dice	1401282	32 x 32 x 32	12	0.3171	0.9250

load_shape is in that table because a metric is only comparable to another measured on the same grid. With one model it is a note; with two it is the difference between a comparison and two unrelated numbers, and evaluate() says so out loud when the rows disagree.

One number for six scans is the least informative view of them. evaluate_cases() with group_by pools each patient’s studies into one score, the way a radiologist would read a patient rather than a file:

group	dice
patient06	0.9151
patient11	0.9148
patient10	0.9448

and without group_by it is a row per scan, with summarise_cases() turning those into the distribution:

metric	n	mean	sd	median	min	max	worst_case
dice	6	0.9250	0.0155	0.9154	0.9142	0.9476	patient11_study1

The minimum and the name beside it are the two numbers worth reading. A mean of 0.925 says the model works; worst_case says which scan to open first, and on a real dataset that is where the hour goes.

Looking at a slice

A volume has no single picture, so plot_masks() asks which slice to draw and refuses a volume without one rather than picking the middle:

masks = engine.predict("lesion", split="validation", space="source")
figure = pp.plot_masks(images, prediction=np.stack(masks), truth=truth, slice=16)

The refusal is the useful part. A middle slice is a reasonable default for a sphere in the centre of a scan and a terrible one for anything real - a lesion in the apex of a lung is not at depth // 2, and a page of pictures that quietly all missed it looks exactly like a page of pictures where the model found nothing.

Millilitres, not voxels

The size of a finding is reported in millilitres. A voxel count means nothing outside the scanner that produced it - the same lesion is 4,000 voxels on one machine and 12,000 on another - and resampling is what makes the conversion trivial: at target_spacing: [1, 1, 1] a voxel is exactly one cubic millimetre.

# On the model's grid, where target_spacing has made a voxel one cubic millimetre.
prediction = engine.predict("lesion", split="validation", space="model")
millilitres = [float((case.argmax(axis=-1) == 1).sum()) / 1000 for case in prediction]

case	predicted ml	true ml	error
patient06_study1	0.332	0.268	+24%
patient06_study2	0.332	0.268	+24%
patient11_study1	0.589	0.524	+12%
patient11_study2	0.572	0.524	+9%
patient10_study1	1.013	0.905	+12%
patient10_study2	1.024	0.905	+13%

Dice 0.914 to 0.948, and volumes 9% to 24% too large. That is not a contradiction. Dice is dominated by the interior of an object and forgiving about its boundary, and for a small lesion the boundary is most of it: a model that draws every lesion one voxel too wide scores well and reports a volume that is wrong by a quarter. If the question is is there a lesion, the Dice is the relevant number. If it is has it grown since March, this table is, and neither follows from the other.

Which grid the counting happens on is not a detail, and getting it wrong is easy. The line above counts on the model’s grid, where target_spacing has made every voxel one cubic millimetre. Count the same lesions on the source grid and half these scans have voxels of 1 x 1 x 2.5 mm, so the identical count means two and a half times the volume. Writing this page, counting source-grid voxels as though they were cubic millimetres reported patient11 as 53% too small where the honest figure is 12% too large - a number that contradicted its own Dice, and the only reason it was caught is that Dice 0.91 and a volume off by half cannot both be true. That is the section’s own lesson arriving uninvited.

Writing a mask back onto its scan

save_volumes() writes a prediction as NIfTI carrying the affine of the scan it came from. Without that affine a mask cannot be laid over its scan by any viewer or registration tool - it either lands in the wrong place or is refused, and nothing about the file looks wrong.

Carrying the affine means the mask has to be on the same grid as the scan, and after resampling it is not: these scans are 48 x 48 x 32 voxels while the model works at 32 x 32 x 32. space="source" maps each prediction back out of the model’s grid onto the scan’s, undoing the resampling and the crop:

masks = engine.predict("lesion", split="validation", space="source")
len(masks)        # 6
masks[0].shape    # (48, 48, 32, 2)

A list rather than an array, because scans differ in size and a stack needs one shape. That is the API being explicit rather than resizing them to match, which is how a mask ends up describing anatomy it was not computed from.

Every scan here is called ct.nii.gz inside its own directory, so names taken from the references would collide and every mask would be written to one path - an error rather than five silently lost files. The case keys are what to pass instead:

written = pp.save_volumes(masks, "predictions", reference=images, cases=keys)
[path.name for path in written]
# ['patient06_study1.nii.gz', 'patient06_study2.nii.gz', ...]

and the written masks carry the scans’ geometry, which is what makes them open on top of them - volume_spacing() says (1.0, 1.0, 1.0) and volume_shape() says (48, 48, 32) for patient06_study1, which is that scan’s own grid rather than the model’s.

Two things are lost on the way back and it is worth knowing which. Interpolation softens a boundary, so the mask is slightly smoother than the model drew it. And where reading cropped anatomy away, the mask returns padded with background - which means not examined, not nothing there. If that distinction matters, give the model an input_shape that covers the anatomy.

DICOM, which is how the data usually arrives

The scans above are NIfTI because that is what research datasets ship. Data out of a hospital arrives as a directory of DICOM files per series, and this package reads those as volumes too - sorting them by their position in the patient rather than by filename, which is the difference between a scan and a shuffled deck.

Point the specification at directories of slices and nothing else changes. Before training, run the export past describe_series():

from pyplatypus.data.dicom_series import describe_series

for directory in sorted(Path("dicom_export").iterdir()):
    report = describe_series(directory)

It reads no pixels - every check runs on the headers - so it stays usable on a hundred gigabytes of data, and it is the first thing worth running on a new export. An archive of a hundred cases will contain a few with a missing slice, duplicated files, or two series sitting in one folder. A missing slice does not shorten a volume: it moves everything past the gap, and no metric would show it.

read_dicom_series() is the reader itself, looks_like_dicom_series() the test that a directory is one, and series_spacing() the slice spacing the headers imply. Unlike NIfTI, which is reliably named, DICOM is detected by content, because exports so often carry no useful extension.

What this leaves out

Resampling is offered, not imposed: drop target_spacing and the scans are cropped or padded to input_shape at their own voxel sizes, which is the right choice when they already agree.

What it costs to impose it, and the trap in doing so, were measured on MSD Spleen - 41 abdominal CTs whose slice thickness runs 1.5 mm to 8 mm and whose voxels differ in volume by 5.9 times, which is the variation target_spacing exists for. Three arms, one input_shape of 64 cubed, every one scored against the truth as stored:

what the model reads	dice	worst case	volume error
the whole scan, squashed into the shape	0.709	0.014	87 ml
`target_spacing: [1.5, 1.5, 1.5]`	0.000	0.000	231 ml
`target_spacing: [7.5, 7.5, 7.5]`	0.822	0.770	29 ml

of a spleen averaging 234 ml.

The middle row is the trap and it is total. Resampling does not only change the scale: it changes what input_shape means. Without it, the shape is the whole scan squashed to fit. With it, the shape is a box of so many millimetres taken from the middle - 64 voxels at 1.5 mm is a 96 mm cube - and a spleen sits up to 165 mm off the centre of the scan and is up to 143 mm across. The organ is outside the box in every case, so the model finds nothing at all, three runs out of three.

Choose the field of view in millimetres first, then the spacing and the shape that cover it. Here that is a 480 mm box, and it is worth +0.11 of Dice over resizing and a third of the volume error - 14% of the organ against 30%.

The difference in reliability is larger than the difference in the mean: across three runs the resampled arm’s worst case stayed between 0.72 and 0.77, while the resized arm’s was 0.00, 0.58 and 0.01. It does not score a little lower; it loses a whole scan, a different one each time. Which scan is not explained by its geometry - the two volumes covering the most patient came through fine - and that is an open question rather than a tidy ending.

A coarse 7.5 mm is a poor way to see a boundary, so the arm this did not run is a finer spacing with splits, which is how a large field of view is covered without a large input. #201.

Several volumes per sample - one per modality, as BraTS ships - works through channels_from, where the order is stated rather than inferred.

And one thing above is smaller than it looks: a validation set of six scans from three patients is enough to show a distribution and nowhere near enough to measure one. The spread in the table above is a demonstration of a mechanism, not an estimate of anything.