What Dice Misses

Why look beyond Dice

My master’s research focuses on improving boundary detection in U-Net architectures for medical image segmentation. A comparison reduced to a single Dice column leaves the part I want to measure unanswered: where the boundary failed and by how many millimetres.

In a comparison of U-Net architectures I published here, the standard U-Net achieved the highest mean Dice. Inspecting the masks, however, separated very different errors: leakage into the background, a peripheral halo, a missing lobe, a discontinuity, and failure on a small target. All of them change overlap, but the number does not say which one occurred, where it occurred, or how far the predicted boundary was from the reference.

Dice is simple, avoids counting the many true negatives in the background, and answers the question it was designed to answer: how much do two regions overlap? Its answer becomes incomplete when it represents segmentation quality on its own. The medical image validation literature documents this mismatch between the selected metric and the actual interest of an application.^{1 2}

I want to mark the edge of its field of view and identify the questions that require other measures.

What Dice measures

An overlap statistic

Let \(G\) be the binary reference mask and \(P\) the predicted mask. Proposed by Lee Dice in 1945 to measure ecological association,³ the coefficient takes the following form in segmentation:

\[ \operatorname{Dice}(P,G) = \frac{2|P\cap G|}{|P|+|G|} = \frac{2TP}{2TP+FP+FN}. \]

When at least one mask is nonempty, the value lies between zero and one: it is one for identical nonempty masks and zero when there is no intersection. When \(P\) and \(G\) are nonempty, it is also the pixel- or voxel-level \(F_1\): the harmonic mean of precision and recall. Three properties explain much of its popularity.

First, true negatives are absent from the formula. In an image where the anatomy of interest occupies a small region, correctly classified background does not dominate the result as it would dominate accuracy. Second, the measure is symmetric: swapping reference and prediction leaves it unchanged. Third, its geometric interpretation is direct and does not depend on the physical unit of the image.

The same formula omits error direction and location. False positives and false negatives enter the same denominator, so Dice does not show whether a model tends to expand or contract a structure. Coordinates are absent from the formula: a false positive touching the boundary and one thirty millimetres away carry the same weight. Dividing by the sum of the volumes also makes the penalty relative to size, so the same absolute error consumes a larger fraction of a small structure.¹

Dice and IoU tell the same story on a different scale

Intersection over union, or IoU, is

\[ \operatorname{IoU}(P,G)=\frac{|P\cap G|}{|P\cup G|} =\frac{TP}{TP+FP+FN}. \]

For a pair of masks with a nonempty union, Dice and IoU are linked exactly by

\[ \operatorname{Dice}=\frac{2\operatorname{IoU}}{1+\operatorname{IoU}}, \qquad \operatorname{IoU}=\frac{\operatorname{Dice}}{2-\operatorname{Dice}}. \]

A mask with higher Dice therefore also has higher IoU. Reporting both does not provide two independent geometric views; it provides a monotonic transformation of the same counts.⁴ There is one qualification: because the transformation is nonlinear, means of per-case Dice and IoU need not preserve differences or rankings after aggregation. That is an effect of the statistical summary, not new spatial information.

Metric and loss function should also be kept separate. The Dice above compares binary masks after a decision rule. A soft Dice loss operates on probabilities, depends on reduction across classes and batches, and usually includes a smoothing term. Training with one of those versions does not settle how final evaluation should be performed; it only defines part of the optimization problem.

The same Dice, different errors

Consider a \(64\times64\)-pixel square reference with area \(4096\). I constructed three predictions, each also containing \(4096\) pixels and sharing the same \(3968\)-pixel intersection with the reference:

1. Shift: the whole square moves two pixels to the right.
2. Island: a strip of \(128\) pixels is removed from the boundary and placed far from the object, creating a second component.
3. Hole: \(128\) pixels are removed from the interior and restored as a strip attached to the outer boundary.

By construction, all three predictions have

\[ \operatorname{Dice}=\frac{2\cdot3968}{4096+4096}=0.96875, \qquad \operatorname{IoU}=\frac{3968}{4224}\approx0.93939. \]

The code below reproduces the example using only NumPy. It uses a four-neighbour inner contour, distances between pixel centres, the maximum of the two directional percentiles for HD95, and a pooled mean for ASSD. Surface Dice is an illustrative approximation that gives every boundary pixel the same weight. The code also materializes all pairwise distances, so it is neither a reference implementation nor a solution for real medical volumes.

import numpy as np


def boundary(mask: np.ndarray) -> np.ndarray:
    """Four-neighbour inner boundary of a 2-D binary mask."""
    m = np.pad(mask, 1, constant_values=False)
    interior = (
        m[1:-1, 1:-1]
        & m[:-2, 1:-1]
        & m[2:, 1:-1]
        & m[1:-1, :-2]
        & m[1:-1, 2:]
    )
    return mask & ~interior


def directed_distances(source: np.ndarray, target: np.ndarray) -> np.ndarray:
    """Euclidean distance from each source-boundary pixel to the target."""
    a = np.argwhere(boundary(source))
    b = np.argwhere(boundary(target))
    squared = ((a[:, None, :] - b[None, :, :]) ** 2).sum(axis=2)
    return np.sqrt(squared.min(axis=1))


def count_components(mask: np.ndarray) -> int:
    """Number of four-connected foreground components."""
    remaining = set(map(tuple, np.argwhere(mask)))
    total = 0
    while remaining:
        total += 1
        stack = [remaining.pop()]
        while stack:
            row, column = stack.pop()
            for neighbour in (
                (row - 1, column),
                (row + 1, column),
                (row, column - 1),
                (row, column + 1),
            ):
                if neighbour in remaining:
                    remaining.remove(neighbour)
                    stack.append(neighbour)
    return total


def count_holes(mask: np.ndarray) -> int:
    """Background components that do not reach the image border."""
    padded_background = np.pad(~mask, 1, constant_values=True)
    return count_components(padded_background) - 1


def metrics(reference: np.ndarray, prediction: np.ndarray, tolerance=2.0):
    intersection = np.logical_and(reference, prediction).sum()
    union = np.logical_or(reference, prediction).sum()
    ref_to_pred = directed_distances(reference, prediction)
    pred_to_ref = directed_distances(prediction, reference)
    pooled = np.concatenate((ref_to_pred, pred_to_ref))
    return {
        "Dice": 2 * intersection / (reference.sum() + prediction.sum()),
        "IoU": intersection / union,
        "HD95": max(np.percentile(ref_to_pred, 95), np.percentile(pred_to_ref, 95)),
        "ASSD": pooled.mean(),
        "SD@2": np.mean(pooled <= tolerance),
        "RVE_pct": 100 * (prediction.sum() - reference.sum()) / reference.sum(),
        "components": count_components(prediction),
        "holes": count_holes(prediction),
    }


reference = np.zeros((192, 192), dtype=bool)
reference[48:112, 48:112] = True

shift = np.zeros_like(reference)
shift[48:112, 50:114] = True

island = reference.copy()
island[48:112, 110:112] = False
island[144:152, 144:160] = True

hole = reference.copy()
hole[76:84, 72:88] = False
hole[48:112, 112:114] = True

for name, prediction in {
    "shift": shift,
    "island": island,
    "hole": hole,
}.items():
    values = metrics(reference, prediction)
    text = " ".join(f"{key}={value:.3f}" for key, value in values.items())
    print(name, text)

# shift  Dice=0.969 IoU=0.939 HD95=2.000  ASSD=1.000 SD@2=1.000 RVE_pct=0.000 components=1.000 holes=0.000
# island Dice=0.969 IoU=0.939 HD95=57.347 ASSD=4.887 SD@2=0.919 RVE_pct=0.000 components=2.000 holes=0.000
# hole   Dice=0.969 IoU=0.939 HD95=27.000 ASSD=2.626 SD@2=0.914 RVE_pct=0.000 components=1.000 holes=1.000

The table adds boundary distances and topology to the overlap results:

The example was designed to isolate geometry, not to represent clinical risk. Dice and IoU are identical and high, while volume error is zero in all three cases. HD95 reacts to the remote island. Component counting records the prediction’s two objects, while hole counting captures a change that component count would miss.

Surface Dice at two pixels considers the entire shift acceptable because no part of its boundary lies beyond the tolerance. Between the two structural defects, it gives the island a slightly higher value than the hole, whereas HD95 judges the island far worse because of the remote component. One metric counts the fraction of surface within tolerance; the other summarizes the tail of the distances. The ordering depends on the property selected for the task: acceptable boundary extent, severe local error, connectivity, or something else.

The relative price of a pixel

Object size creates another effect. If a reference has area \(n\) and an erosion removes \(e\) pixels, the prediction is contained in the reference and \(|P|=|P\cap G|=n-e\). In that case,

\[ \operatorname{Dice}=\frac{2(n-e)}{2n-e}. \]

The ratio \(e/n\) tends to be larger for small structures, where the boundary accounts for a greater fraction of the area. Applying one four-neighbour erosion to two disks on the same grid makes the effect visible:

def disk(radius: int, size=128) -> np.ndarray:
    y, x = np.ogrid[:size, :size]
    center = (size - 1) / 2
    return (y - center) ** 2 + (x - center) ** 2 <= radius**2


def four_neighbour_erosion(mask: np.ndarray) -> np.ndarray:
    m = np.pad(mask, 1, constant_values=False)
    return (
        m[1:-1, 1:-1]
        & m[:-2, 1:-1]
        & m[2:, 1:-1]
        & m[1:-1, :-2]
        & m[1:-1, 2:]
    )


for radius in (10, 40):
    reference = disk(radius)
    prediction = four_neighbour_erosion(reference)
    intersection = np.logical_and(reference, prediction).sum()
    dice = 2 * intersection / (reference.sum() + prediction.sum())
    removed = reference.sum() - prediction.sum()
    print(f"radius={radius:2d} area={reference.sum():4d} removed={removed:3d} Dice={dice:.3f}")

# radius=10 area= 316 removed= 56 Dice=0.903
# radius=40 area=5024 removed=224 Dice=0.977

The boundary operation is the same in both cases, but the Dice drop is much larger for the small disk. The result is consistent with the coefficient’s definition of overlap, but it prevents the same Dice difference from being interpreted as the same geometric difference across structures of different sizes.

Distance to the boundary

To evaluate the contour, we first need to represent it. Let \(\partial P\) and \(\partial G\) be the surfaces of the prediction and reference, and let

\[ d(x,S)=\inf_{y\in S}\lVert x-y\rVert_2 \]

be the distance from point \(x\) to the closest point on surface \(S\). The two directed collections,

\[ D_{P\to G}={d(p,\partial G):p\in\partial P}, \qquad D_{G\to P}={d(g,\partial P):g\in\partial G}, \]

retain information that pixel counts discard. Different summaries of those collections answer different questions.⁴

Hausdorff and HD95

The symmetric Hausdorff distance takes the largest error in either direction:

\[ HD(P,G)=\max\bigl(\sup D_{P\to G},\sup D_{G\to P}\bigr). \]

It finds the point of greatest disagreement, so it reacts strongly to a remote island or a single stray pixel. That sensitivity may expose a serious failure or annotation noise. Under a directional, area-weighted convention, replacing the maximum with the 95th percentile discards the upper 5% of each distance distribution; other conventions truncate a different tail. A distant false component may not appear in HD95 if it occupies a small enough fraction of the surface.

Even the name “HD95” does not specify an implementation. One can take the maximum of the directional percentiles, as in the example above, or pool both collections before taking the percentile; one can measure between voxel centres or mesh elements; one can weight by physical surface area or not. A comparison of five libraries found systematic differences caused by choices of this kind.⁵ The percentile, boundary extraction, weighting, connectivity, library, and version therefore belong in the experimental method.

Average symmetric surface distance

One convention for average symmetric surface distance, or ASSD, pools both directions and weights each surface element:

\[ \operatorname{ASSD}_{\mathrm{pool}}= \frac{ \sum_{p\in\partial P} a_p d(p,\partial G) {}+\sum_{g\in\partial G} a_g d(g,\partial P) }{ \sum_{p\in\partial P}a_p+\sum_{g\in\partial G}a_g }, \]

where \(a_p\) and \(a_g\) are lengths in 2-D or areas in 3-D. It describes typical contour separation and is less dominated by a single extreme than Hausdorff. A severe local failure, however, can be diluted among thousands of elements close to the reference.

Another convention also called ASSD first computes a mean in each direction and then averages the two. The formulas differ when the surfaces have unequal areas. The illustrative code used the pooled form with \(a=1\); a paper should state the convention and preserve physical surface areas.

Surface Dice

Surface Dice at tolerance \(\tau\) asks what fraction of both surfaces lies within a distance deemed acceptable from the other. With \(A(\cdot)\) denoting surface area,

\[ SD_{\tau}(P,G)= \frac{ A\bigl({p\in\partial P:d(p,\partial G)\le\tau}\bigr) {}+A\bigl({g\in\partial G:d(g,\partial P)\le\tau}\bigr) }{A(\partial P)+A(\partial G)}. \]

The metric was proposed for delineating organs at risk in radiotherapy, with structure-specific tolerances estimated from variation between specialists.⁶ The Medical Segmentation Decathlon also combined volumetric Dice with normalized surface Dice, using different physical tolerances by task and anatomy, selected through clinical feedback.⁷

The parameter \(\tau\) formalizes contour tolerance and must come from the task: annotation variability, editing tolerance, acquisition resolution, or a criterion defined with specialists. It also creates a cutoff: points just below and above \(\tau\) receive different outcomes, and the metric does not say how far beyond the limit an error lies.

In three-dimensional medical images, all these distances must respect physical spacing. In a volume with \(0.8\times0.8\times5.0\) mm voxels, one step in the axial plane does not equal one step between slices. Reporting “3 voxels” without a direction, or computing on a resampled array without recording the transformation, answers a different question. For surface Dice, counting boundary voxels is also insufficient: surface elements should be weighted by their physical length or area.⁶

What boundary metrics also leave out

Replacing Dice with one column of HD95 does not solve the underlying problem. Contour distance remains only one family of answers; error direction, object count, topology, and volume may remain invisible.

Error direction and volume

Precision and recall separate the two terms that Dice combines:

\[ \operatorname{Precision}=\frac{TP}{TP+FP}, \qquad \operatorname{Recall}=\frac{TP}{TP+FN}. \]

If missing tissue and including background have unequal costs, that asymmetry should appear in the report or in a prespecified weighted measure. Signed relative volume error,

\[ RVE=\frac{V_P-V_G}{V_G}, \]

shows a tendency toward oversegmentation (positive) or undersegmentation (negative); multiply the result by 100 to express it as a percentage. Precision is undefined for an empty prediction; recall and RVE are undefined for an empty reference unless a declared convention is imposed. RVE should only answer a volume question: two disjoint masks can have \(RVE=0\). If volumetry is the outcome, absolute error in millilitres or cubic centimetres and patient-level bias are more interpretable measures; Dice is a poor proxy.²

Objects and topology

In multifocal segmentation, a Dice computed over the full volume implicitly weights each lesion by its voxel count. One large, correctly segmented mass can hide several missing small lesions. Object-level evaluation must first match predicted and reference components using a declared rule; it can then report lesion-wise recall, precision or false positives per scan, and delineation quality only for detected lesions. Recent results on metastases in PET/CT document how global metrics favour larger components.⁸

For vessels, airways, neurons, and other network structures, connectivity may matter more than area. Centerline Dice, or clDice, intersects morphological skeletons with the masks. Here \(S_P\) and \(S_G\) denote the prediction and reference skeletons, not the surfaces used in the previous section:

\[ T_{\mathrm{prec}}=\frac{|S_P\cap G|}{|S_P|}, \qquad T_{\mathrm{sens}}=\frac{|S_G\cap P|}{|S_G|}, \qquad clDice=2\frac{T_{\mathrm{prec}}T_{\mathrm{sens}}}{T_{\mathrm{prec}}+T_{\mathrm{sens}}}. \]

Those properties concern tubular structures under specific assumptions.⁹ For compact tumours or organs, component count, hole count, Euler characteristic, Betti numbers, or an application-specific rule may be more appropriate. The relevant topology must be defined before its measure is selected.

When implementation changes the question

A metric value also depends on the protocol that produced it. Each case below changes the quantity being estimated.

1. Empty masks. If both reference and prediction are empty, the Dice formula gives \(0/0\); depending on the convention, the result may be one, zero, or NaN. If only one mask is empty, surface metrics do not have two surfaces to compare. In datasets where the structure may be absent, presence or absence should be assessed as a classification task, with delineation assessed in the applicable cases. Each metric’s result for empty masks must be specified in advance. Silently removing these cases also removes false positives and can inflate the result.¹⁰
2. Decision threshold. Turning probabilities into a mask requires a threshold. Tuning it on the test set transfers information from the final set to the model; the threshold and any post-processing should be selected on validation data and frozen before evaluation.
3. Resampling. Categorical masks require suitable interpolation, usually nearest neighbour. Orientation, origin, and spacing must remain aligned; resampling can smooth or move boundaries and alter precisely the metrics under study.
4. Two or three dimensions. Averaging Dice per slice, computing one whole-volume mask, and pooling every voxel in the dataset estimate three different quantities. The unit should follow the use case: slice, scan, structure, or patient.
5. Aggregation. A global Dice obtained by concatenating all voxels gives more weight to cases and classes with larger foreground masks, that is, larger \(|P|+|G|\). Computing by patient before summarizing respects the hierarchy of the data; classes and structures should also remain separate when an average would conceal a specific failure.²
6. Uncertainty and hard cases. A mean without a distribution does not show tails or heterogeneity. The median and interquartile range, a patient-level bootstrap confidence interval, and strata defined before analysis—such as structure size and acquisition centre—show where a gain occurred. Qualitative cases should follow a selection rule, not convenience.
7. Executable definition. The name, formula, parameters, empty-case handling, unit, library, and version should be recorded. HD95 from two libraries should not be presumed identical.⁵

A minimum evaluation protocol

The metric set should follow the decision that segmentation will support. An initial configuration, to be narrowed or expanded for the anatomy and use case, can be organized by the questions below.

I would begin with the analysis unit and presence rule. I would then choose one overlap measure and add only the families tied to the intended use. Contour tolerance would be prespecified from clinical acceptability, inter-annotator variation, or an external task protocol without examining test predictions. Decision threshold and post-processing would be selected using training and validation data; empty-case rules would also be frozen before testing. Results would be computed by patient and class, with uncertainty and size stratification. Finally, I would publish the implementation or, at minimum, its executable definition and tests on known synthetic masks.

The result is an error profile distributed across several measures. Metrics Reloaded similarly recommends selecting metrics from a problem fingerprint and combining complementary families.² The table becomes wider, but it shows why two models differ.

Synthesis

The three-mask experiment shows how I would frame a comparison. The predictions tied exactly on Dice, IoU, and volume, but the additional measures separated different properties: HD95 exposed the distance to the remote component, surface Dice quantified the fraction within tolerance, and the topological descriptors recorded components and holes. An ordering becomes meaningful only after the task’s property of interest has been stated.

For my problem, a claim that an architecture “improves boundaries” should answer at least three questions on the same test set: how much overlap changed, whether surface distance fell and by how many millimetres, and which structure sizes benefited. If there are multiple lesions, I would add how many were actually found. A higher Dice answers only the first.

Two questions remain open. Do methods designed to refine contours reduce physical distances without losing lesion-wise recall? And do apparent gains persist after stratifying by size, or are they concentrated in the large structures for which Dice is more forgiving?

References