Hausdorff Distance

The Hausdorff distance measures how far two sets are from each other by their worst-placed point: it is the largest distance from any point of either set to the nearest point of the other. Two shapes are within Hausdorff distance $\varepsilon$ exactly when each lies inside the other grown by $\varepsilon$. In computer vision it compares segmentation boundaries with a reference, as in medical image segmentation and challenges such as BraTS; it compares point sets and surfaces, as the Chamfer distance article describes; and it is an edge-based matching criterion: the mutual information registration article lists Huttenlocher’s Hausdorff distance among the methods that compare edge maps, and the Chamfer Distance article describes Huttenlocher, Klanderman and Rucklidge’s (1993) paper. This article covers the distance as a general set distance, its directed, percentile, modified and partial forms, the conventions with which libraries compute it on pixel grids, and matching under translation with a distance transform.

Formula

Let $A$ and $B$ be non-empty sets of points and $d$ a distance between points, the Euclidean distance unless stated. The directed Hausdorff distance from $A$ to $B$ is

$$h(A, B) = \sup_{a \in A} \, \inf_{b \in B} \, d(a, b),$$

and the Hausdorff distance is the larger of the two directions,

$$H(A, B) = \max\big\{ h(A, B),\; h(B, A) \big\}.$$

Ribera et al. (2019) state it in this sup–inf form (their Eq. 1) and note that on a bounded pixel grid the suprema and infima are attained, so they become maxima and minima. For finite sets, such as the pixels of two masks,

$$h(A, B) = \max_{a \in A} d(a, B), \qquad d(a, B) = \min_{b \in B} \|a - b\|.$$

The value is a length, in pixels or, with the voxel spacing applied, in millimetres. It ranges from 0, for identical closed sets, to the diameter of the image or volume; smaller is better. On closed bounded sets $H$ is a metric: it is non-negative, zero only for equal sets, symmetric, and satisfies the triangle inequality (Ribera et al. list the four properties, citing Attouch et al.; the Chamfer Distance article attributes the statement to Huttenlocher et al., 1993). The directed distance is not symmetric: $h(A, B) = 0$ only says that $A$ lies inside $B$.

Writing $B_\varepsilon = \{x : d(x, B) \le \varepsilon\}$ for $B$ dilated by a disc of radius $\varepsilon$, it follows from the definition that $h(A, B) \le \varepsilon$ exactly when $A \subseteq B_\varepsilon$, so that $H(A, B)$ is the smallest $\varepsilon$ for which each set lies inside the other’s dilation (our derivation). This is the link to mathematical morphology, and it is why a distance transform computes the directed distance directly: $h(A, B)$ is the largest value of the distance map of $B$ read at the points of $A$.

Intuition

The directed distance asks how far the farthest point of $A$ has to travel to reach $B$. Every point of $A$ finds its nearest neighbour in $B$, and only the worst of those distances counts; the Chamfer distance averages the same nearest-neighbour distances instead. The symmetric version asks the question both ways, so it penalizes a prediction that reaches where the reference does not (a false-positive spur) and a reference part the prediction misses (a false-negative gap) alike.

Because only the maximum counts, the Hausdorff distance is a worst-case guarantee: every point of each set is within $H$ of the other set. The same property makes it sensitive to a single stray point. Ribera et al. call this a high sensitivity to outliers and illustrate two point configurations with the same Hausdorff distance despite very different typical errors (their Figure 2); the Medical Image Segmentation and Chamfer Distance articles show one far voxel or point setting the value.

Interpretation

What the number means depends on which points enter it.

• Boundaries or filled regions. Segmentation evaluation compares boundary points: the outlines of two masks in 2D, their surfaces in 3D, as in the definition Medical Image Segmentation gives. scikit-image’s hausdorff_distance(image0, image1) instead uses every nonzero pixel of each image (its source takes np.nonzero of both inputs). For two discs the two choices give the same Hausdorff distance, $\delta + |r_1 - r_2|$ for radii $r_1, r_2$ and centres $\delta$ apart (our geometry), but not when one region has a hole: the filled-set distance is the radius of the hole, the largest distance from a point of the hole to the surrounding region, while the boundary distance runs from the hole’s edge to the outer edge. The example below shows both.
• Directed values. $h(\text{prediction}, \text{reference})$ is large when the prediction has a false-positive part far from the reference; $h(\text{reference}, \text{prediction})$ is large when part of the reference is missed. MONAI’s directed=True returns the first only. In matching, the forward distance $h(\text{model}, \text{image})$ asks whether every model point is near some image feature, and extra image features cannot increase it.
• Units and resolution. On a pixel grid, rasterization alone moves the boundary by up to about a pixel, so differences below a pixel say little; with anisotropic voxels the spacing is part of the metric, as Medical Image Segmentation shows.
• Scale. The Hausdorff distance is not normalized by object size: 5 mm is a small error for a liver and a large one for a 10 mm lesion.

When to Use It

Use the Hausdorff distance when the worst boundary error matters, for example a margin that a treatment plan must respect, and when the sets are reasonably clean, since one outlier decides the value. Report it alongside an overlap score such as the Dice coefficient, which counts pixels without regard to where errors lie, and alongside an average distance, which describes the typical error. Metrics Reloaded recommends HD or a percentile of HD as boundary metrics when distance errors should be penalized with a focus on outliers, as Medical Image Segmentation summarizes.

Use a percentile form (HD95) when isolated outliers should not decide the score, a partial (ranked) form when part of a model may be legitimately missing, as with occlusion in matching, and an average form when the whole contour matters. For point clouds sampled from surfaces, where sampling density varies, the Chamfer distance and the F-score at a threshold answer related questions, as that article explains.

Failure Cases

• One stray point. A single false-positive voxel far from the target sets the Hausdorff distance, and so does one erroneous voxel in the reference; BraTS adopted HD95 because a few small false-positive regions can dominate the maximum, as Tumor Segmentation reports.
• A large error the percentile keeps. HD95 discards only the worst 5% of the boundary. When a discrepancy at one distance covers more of the boundary than that, HD95 equals HD; in the example, a hole whose edge is 15/55 of the annulus boundary keeps the analytic HD95 at the full distance, and every boundary convention’s HD95 within 0.3 pixel of its HD.
• Filled sets hide what boundaries show. On filled regions, interior points that lie inside the other region contribute zero, so the filled-set distance can be much smaller than the boundary distance, as with the hole below (15 against 25 pixels analytically).
• Empty sets. The distance to an empty set is undefined. scikit-image returns 0 when both images are empty and infinity when only one is; DeepMind’s compute_robust_hausdorff returns infinity when either surface is empty; MONAI returns NaN when both are empty and infinity when one is, with a warning (MONAI 1.6.1 source). A method that misses a lesion therefore gets an infinite, NaN or, after a challenge’s penalty, a fixed value, and how such cases are aggregated decides rankings, as BraTS and Medical Image Segmentation describe.
• Convention differences. Border extraction, surface weighting and the percentile rule differ between implementations and change values by a fraction of a pixel to about a pixel on the same masks (the example below).
• Small structures. With few boundary pixels, a percentile falls between widely spaced distances, and a one-pixel change in the mask is a large fraction of the object.

Variants

Directed and symmetric

The directed distance $h(A, B)$ is used when only one direction is meaningful, for example to ask only whether a prediction strays outside the reference; MONAI offers it as directed=True. The symmetric $H$ is what scikit-image and DeepMind’s code return, and MONAI’s default. BOP’s symmetry test, in the object pose estimation article, compares an object’s vertices with their transformed copy by a Hausdorff distance.

Percentile Hausdorff distance (HD95)

Replacing each directed maximum by its 95th percentile, and taking the larger of the two, gives HD95; Medical Image Segmentation describes it as a form of the quantile method of Huttenlocher et al., following Taha and Hanbury (2015). Implementations differ in how the percentile is taken:

• Interpolation. MONAI’s HausdorffDistanceMetric(percentile=95) applies torch.quantile to each direction’s distances and takes the maximum; torch.quantile, like NumPy’s default np.percentile, interpolates linearly between order statistics.
• Inverse CDF with area weights. DeepMind’s compute_robust_hausdorff sorts the distances, accumulates each surface element’s length (2D) or area (3D), normalizes the sum, and takes the first distance whose cumulative weight reaches the percentile (np.searchsorted), without interpolation.
• Pooling. The panoptica toolkit in BraTS’s evaluation code takes one percentile of both directions’ distances pooled, which the BraTS article compares with the directed-maximum form.

Partial (ranked) Hausdorff distance

For matching, where part of the model may be occluded, the maximum is replaced by a ranked value. With $|A| = n$,

$$h_K(A, B) = \operatorname*{K^{th}}_{a \in A} \, d(a, B),$$

the $K$-th smallest of the nearest-neighbour distances, often written as a fraction $f = K/n$; $K = n$ gives the directed distance. The Chamfer Distance article reports that Huttenlocher et al. extended the distance to compare portions of the sets (their Section III); the ranked form above is the formulation used here. The ranked value tolerates up to $n - K$ points that are missing from $B$: if at least $K$ points of $A$ lie exactly on $B$, $h_K(A, B) = 0$ whatever the other points do. In the forward direction from model to image it also ignores clutter, since added image points can only lower distances. A percentile Hausdorff distance is the same idea with $K$ set by a percentage; MONAI’s directed=True, percentile=80 computes $h_{0.8}$ up to its interpolation. Two of the same authors, Huttenlocher and Rucklidge, wrote an earlier Cornell technical report (1992), A Multi-Resolution Technique for Comparing Images Using the Hausdorff Distance.

Modified Hausdorff distance

scikit-image implements the modified Hausdorff distance (MHD) of Dubuisson and Jain (1994) as hausdorff_distance(..., method='modified'), which, per its docstring and source, computes the mean nearest-neighbour distance in each direction and returns the larger:

$$\mathrm{MHD}(A, B) = \max\Big\{ \tfrac{1}{|A|} \textstyle\sum_{a \in A} d(a, B),\; \tfrac{1}{|B|} \sum_{b \in B} d(b, A) \Big\}.$$

The docstring states that Dubuisson and Jain found it to perform better than the directed Hausdorff distance. Medical Image Segmentation notes that Taha and Hanbury’s “average distance” takes the same larger of the two directed means.

Average Hausdorff distances

Averaging instead of maximizing gives a family whose members carry similar names:

• the pooled average symmetric surface distance (ASSD), which divides the sum of both directions’ distances by the total number of boundary points, defined in Medical Image Segmentation;
• the larger of the two directed means, the MHD above;
• the sum of the two directed means, which Ribera et al. (2019) call the average Hausdorff distance (their Eq. 4). They note that it keeps non-negativity, the identity property and symmetry but loses the triangle inequality, and that it is differentiable with respect to the points.

For the same pair of sets these three values differ, so the formula must be stated with the number.

Weighted Hausdorff distance as a loss

Ribera et al. (2019) made the average Hausdorff distance trainable for locating objects without bounding boxes. Their weighted Hausdorff distance replaces the predicted point set by a per-pixel probability map from a fully convolutional network: the first term weights each pixel’s distance to the nearest ground-truth point by its probability, and the second replaces the minimum over predicted points by a generalized mean that penalizes low probabilities near each ground-truth point (their Eqs. 5–7). With binary probabilities, the minimum in place of the generalized mean and $\epsilon = 0$, it reduces to the average Hausdorff distance.

Related distances under other names

• PyTorch3D’s point_reduction="max" returns the square of $H$ with its default squared norm, as the Chamfer Distance article explains.
• BOP’s MSSD is a pose error: the largest deviation of a model vertex between the estimated and the true pose, minimized over the object’s symmetries, as robotic manipulation and object pose estimation define it. It takes a maximum over corresponding vertices, not over nearest-neighbour distances.
• PointNet’s universal-approximation theorem is stated for functions continuous with respect to the Hausdorff distance between point sets.

Computation

For finite sets of sizes $p$ and $q$, brute force computes all $pq$ distances. Faster methods exploit structure:

• Early break. SciPy’s directed_hausdorff(u, v, rng=0) implements the algorithm of Taha and Hanbury’s TPAMI paper (2015): it shuffles both point sets (with a fixed seed by default) and, for each point of u, stops scanning v as soon as a distance below the current maximum is found, because that point can no longer raise the maximum. The worst case stays $O(pq)$; the docstring reports that the authors showed the average runtime to be close to linear in the size of u. The v1.16.3 source carries a comment that the paper’s update condition is incorrect when the inner loop breaks early, and uses a corrected one. The function returns the distance and the indices of one pair of points that attains it.
• Nearest-neighbour trees. scikit-image builds a cKDTree on each set and queries the other; MONAI does the same on CPU. The Chamfer Distance article gives Huttenlocher et al.’s $O((p+q)\log(p+q))$ bound for planar point sets.
• Distance-transform lookup. On a grid, one exact Euclidean distance transform of the complement of $B$ gives $d(x, B)$ for every pixel $x$; reading it at the points of $A$ gives all directed distances. DeepMind’s code and the Medical Image Segmentation example work this way, with the voxel spacing passed to the transform.
• Matching under translation. The forward distance at a translation $t$ is $h_K(A + t, B)$, read from the same distance map of $B$ at the shifted model points. One transform serves every translation, and each candidate then costs $|A|$ lookups and a partial sort. Template matching scores translations by intensity differences or correlation instead; chamfer matching sums the same distance-map values instead of ranking them, as the Chamfer Distance article describes for Barrow et al.

Implementations

Conventions read from the source at the versions given (as of November 2026):

Implementation	Points used	Distance	Percentile	Empty sets
scikit-image 0.25.2 hausdorff_distance	all nonzero pixels of each image	k-d tree, Euclidean, pixel units	none; method='modified' gives the larger directed mean	0 if both empty, inf if one
SciPy 1.16.3 directed_hausdorff	the given point arrays, one direction	early-break exact search, Euclidean	none	not handled specially
DeepMind surface-distance 1f805ce	corner-grid contour (2D) or surface (3D) elements from 2×2(×2) neighbour codes, weighted by length or area	EDT with spacing in mm	inverse CDF on the cumulative weights (compute_robust_hausdorff)	inf
MONAI 1.6.1 HausdorffDistanceMetric	binary-erosion XOR edges, unweighted	k-d tree on CPU, EDT on GPU, with spacing; chessboard and taxicab options	torch.quantile (linear) per direction, then the maximum	NaN if both empty, inf if one, with a warning

In MONAI 1.6.1, compute_hausdorff_distance calls get_edge_surface_distance without its use_subvoxels option, so the Hausdorff metric always uses erosion edges, even when a spacing is given (our reading of hausdorff_distance.py and utils.py). The pooled HD95 of the BraTS evaluation code, panoptica, is described in the BraTS article.

Example

The script (NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2; about 2 seconds) has two parts. Part A computes the Hausdorff distance (HD), HD95 and the modified distance (MHD) of three pairs of rasterized discs with every convention above and compares them with the values of the continuous shapes. A disc is the closed set of pixel centres within its radius:

• P1: a disc of radius 40 against the same disc with a centred hole of radius 15, an annulus;
• P2: discs of radius 40 and 35 with centres 6 pixels apart;
• P3: discs of radius 6 and 4 with centres 2 pixels apart.

The analytic values are closed forms for P1. For two discs, the filled-set and the boundary Hausdorff distances are both the centre distance plus the difference of radii ($6 + 5 = 11$ and $2 + 2 = 4$, our geometry); the code takes the filled value from this formula and the boundary value as the largest quadrature sample, which matches it at the printed precision. The other P2 and P3 values (HD95 and MHD) come from quadrature on the exact circles and discs; HD95 for curves weights the boundary by length. The conventions are:

• scikit-image on all nonzero pixels;
• pixel borders by erosion with SciPy’s default cross, which keeps pixels with a background 4-neighbour and gives an 8-connected border (MONAI’s convention, with linear percentiles), and by erosion with a $3 \times 3$ square, which keeps pixels with a background 8-neighbour and gives a 4-connected border;
• a port of DeepMind’s 2D surface extraction, with its contour-length weights, robust Hausdorff and average surface distance.

The columns give HD; HD95 with linear interpolation and with the inverse CDF, both with unit weights; DeepMind’s length-weighted HD95; and MHD (for the DeepMind row, the larger of its two length-weighted average surface distances). The script checks SciPy, scikit-image and the distance-transform lookups against brute force, DeepMind’s rule with unit weights against NumPy’s inverted_cdf, and the port’s contour length on a $10 \times 10$ square against $4s - 4 + 2\sqrt2$.

import os
os.environ["OMP_NUM_THREADS"] = os.environ["OPENBLAS_NUM_THREADS"] = os.environ["MKL_NUM_THREADS"] = "1"
import numpy as np
from scipy import ndimage as ndi
from scipy.spatial import cKDTree
from scipy.spatial.distance import directed_hausdorff
from skimage.metrics import hausdorff_distance
from skimage.draw import polygon, line

# ---------------------------------------------------------------------------
# DeepMind surface-distance, 2-D path, ported from surface_distance/metrics.py and
# lookup_tables.py at commit 1f805ce (functions compute_surface_distances,
# compute_average_surface_distance, compute_robust_hausdorff,
# create_table_neighbour_code_to_contour_length).
# Copyright 2018 Google Inc. All Rights Reserved.
# Licensed under the Apache License, Version 2.0; modified (ported): 2-D only,
# shortened, scipy.ndimage.correlate in place of the deprecated ndimage.filters.
# Licence: https://github.com/google-deepmind/surface-distance/blob/1f805cea446805e559badc80dd27f34175e4173d/LICENSE
def dm_contour_length_table(spacing):
    t = np.zeros(16)
    vert, horiz, diag = spacing[0], spacing[1], 0.5 * np.hypot(spacing[0], spacing[1])
    for code, length in [("0001", diag), ("0010", diag), ("0011", horiz), ("0100", diag),
                         ("0101", vert), ("0110", 2 * diag), ("0111", diag), ("1000", diag),
                         ("1001", 2 * diag), ("1010", vert), ("1011", diag), ("1100", horiz),
                         ("1101", diag), ("1110", diag)]:
        t[int(code, 2)] = length
    return t

def dm_surface_distances(mask_gt, mask_pred, spacing=(1.0, 1.0)):
    table = dm_contour_length_table(spacing)
    kernel = np.array([[8, 4], [2, 1]])
    idx = np.argwhere(mask_gt | mask_pred)
    lo, hi = idx.min(0), idx.max(0)
    crops = []
    for m in (mask_gt, mask_pred):
        c = np.zeros(hi - lo + 2, np.uint8)
        c[0:-1, 0:-1] = m[lo[0]:hi[0] + 1, lo[1]:hi[1] + 1]
        crops.append(c)
    codes = [ndi.correlate(c, kernel, mode="constant", cval=0) for c in crops]
    borders = [(k != 0) & (k != 15) for k in codes]
    dist = [ndi.distance_transform_edt(~b, sampling=spacing) for b in borders]
    areas = [table[k] for k in codes]
    d_gt_pred, d_pred_gt = dist[1][borders[0]], dist[0][borders[1]]
    a_gt, a_pred = areas[0][borders[0]], areas[1][borders[1]]
    sort = lambda d, a: tuple(np.array(sorted(zip(d, a))).T)
    (d_gt_pred, a_gt), (d_pred_gt, a_pred) = sort(d_gt_pred, a_gt), sort(d_pred_gt, a_pred)
    return dict(distances_gt_to_pred=d_gt_pred, distances_pred_to_gt=d_pred_gt,
                surfel_areas_gt=a_gt, surfel_areas_pred=a_pred, borders=borders)

def dm_average_surface_distance(sd):
    return (np.sum(sd["distances_gt_to_pred"] * sd["surfel_areas_gt"]) / np.sum(sd["surfel_areas_gt"]),
            np.sum(sd["distances_pred_to_gt"] * sd["surfel_areas_pred"]) / np.sum(sd["surfel_areas_pred"]))

def dm_robust_hausdorff(sd, percent):
    out = []
    for d, a in ((sd["distances_gt_to_pred"], sd["surfel_areas_gt"]),
                 (sd["distances_pred_to_gt"], sd["surfel_areas_pred"])):
        idx = np.searchsorted(np.cumsum(a) / np.sum(a), percent / 100.0)
        out.append(d[min(idx, len(d) - 1)])
    return max(out)
# --------------------------------------------------------------------------- end of port

def disc(shape, c, r):
    yy, xx = np.indices(shape)
    return (yy - c[0]) ** 2 + (xx - c[1]) ** 2 <= r ** 2       # closed disc, pixel centres

border_cross = lambda m: m & ~ndi.binary_erosion(m)                          # 4-neighbour test: 8-connected border
border_square = lambda m: m & ~ndi.binary_erosion(m, np.ones((3, 3), bool))   # 8-neighbour test: 4-connected border

def brute_nearest(P, Q, chunk=2000):
    """Distance from each row of P to the nearest row of Q, by brute force."""
    return np.concatenate([np.sqrt(((P[i:i + chunk, None, :] - Q[None]) ** 2).sum(-1)).min(1)
                           for i in range(0, len(P), chunk)])

def gap_bound_ok(d, q=0.95):
    """|linear - inverted_cdf| percentile is at most the gap between the order statistics around rank q."""
    s = np.sort(d); n = len(s)
    i = int(np.floor(q * (n - 1)))
    ranks = [i, min(i + 1, n - 1), int(np.ceil(q * n)) - 1]
    return abs(np.percentile(d, 95) - np.percentile(d, 95, method="inverted_cdf")) <= s[max(ranks)] - s[min(ranks)] + 1e-12

# ---- analytic values of the continuous shapes (quadrature on exact circles and discs)
def analytic_offset_discs(c1, r1, c2, r2, m=720_000, nr=1500, nt=1500):
    th = (np.arange(m) + 0.5) * 2 * np.pi / m
    u = np.stack([np.cos(th), np.sin(th)], 1)
    f12 = np.abs(np.linalg.norm(c1 + r1 * u - c2, axis=1) - r2)     # circle 1 -> circle 2
    f21 = np.abs(np.linalg.norm(c2 + r2 * u - c1, axis=1) - r1)
    rho = (np.arange(nr) + 0.5) / nr
    tt = (np.arange(nt) + 0.5) * 2 * np.pi / nt
    R, T = np.meshgrid(rho, tt, indexing="ij")
    w = R / R.sum()                                                  # area weights of the polar grid
    def filled_mean(ca, ra, cb, rb):
        p = np.stack([ca[0] + ra * R * np.cos(T), ca[1] + ra * R * np.sin(T)], -1)
        return (np.maximum(0, np.linalg.norm(p - cb, axis=-1) - rb) * w).sum()
    delta = np.linalg.norm(c1 - c2)
    return dict(hd_filled=max(r1 + delta - r2, r2 + delta - r1, 0.0),
                mhd_filled=max(filled_mean(c1, r1, c2, r2), filled_mean(c2, r2, c1, r1)),
                hd=max(f12.max(), f21.max()), hd95=max(np.quantile(f12, 0.95), np.quantile(f21, 0.95)),
                mhd=max(f12.mean(), f21.mean()))

pairs = {}
A = disc((128, 128), (64, 64), 40); B = A & ~disc((128, 128), (64, 64), 15)
pairs["P1  disc r=40 vs the same disc with a centred hole r=15"] = (A, B, dict(
    hd_filled=15.0, mhd_filled=2 * (15 ** 3 / 6) / 40 ** 2,          # mean of (15 - rho)+ over the disc
    hd=25.0, hd95=25.0, mhd=15 * 25 / 55))                          # hole circle is 15/55 of the boundary
c1, c2 = np.array([64.0, 60.0]), np.array([64.0, 66.0])
pairs["P2  disc r=40 vs disc r=35, centres 6 px apart"] = (
    disc((128, 128), c1, 40), disc((128, 128), c2, 35), analytic_offset_discs(c1, 40, c2, 35))
c1, c2 = np.array([20.0, 19.0]), np.array([20.0, 21.0])
pairs["P3  disc r=6 vs disc r=4, centres 2 px apart"] = (
    disc((40, 40), c1, 6), disc((40, 40), c2, 4), analytic_offset_discs(c1, 6, c2, 4))

print("Part A: Hausdorff conventions on rasterized discs (pixel units)")
checks = dict(scipy=0.0, skimage=0.0, edt=0.0, dm_unit=0.0, gap_ok=True)
for name, (A, B, an) in pairs.items():
    print(f"\n{name}   [pixels: {A.sum()} and {B.sum()}; B inside A: {bool((B & ~A).sum() == 0)}]")
    print(f"{'convention':34s}{'HD':>8s}{'HD95 lin':>10s}{'HD95 icdf':>10s}{'HD95 len':>10s}{'MHD':>8s}")
    print(f"{'analytic, filled sets':34s}{an['hd_filled']:8.3f}{'-':>10s}{'-':>10s}{'-':>10s}{an['mhd_filled']:8.3f}")
    print(f"{'analytic, boundary curves':34s}{an['hd']:8.3f}{'-':>10s}{'-':>10s}{an['hd95']:10.3f}{an['mhd']:8.3f}")
    # scikit-image: all nonzero pixels
    pa, pb = np.argwhere(A).astype(float), np.argwhere(B).astype(float)
    dab, dba = brute_nearest(pa, pb), brute_nearest(pb, pa)
    hs, hm = hausdorff_distance(A, B), hausdorff_distance(A, B, method="modified")
    checks["skimage"] = max(checks["skimage"], abs(hs - max(dab.max(), dba.max())),
                            abs(hm - max(dab.mean(), dba.mean())))
    print(f"{'scikit-image, all nonzero pixels':34s}{hs:8.3f}{'-':>10s}{'-':>10s}{'-':>10s}{hm:8.3f}")
    # pixel borders, distances by EDT lookup, unit weights
    for label, border in (("border, cross erosion (8-conn.)", border_cross),
                          ("border, 3x3 erosion (4-conn.)", border_square)):
        ba, bb = border(A), border(B)
        d1, d2 = ndi.distance_transform_edt(~bb)[ba], ndi.distance_transform_edt(~ba)[bb]
        qa, qb = np.argwhere(ba).astype(float), np.argwhere(bb).astype(float)
        checks["edt"] = max(checks["edt"], np.abs(d1 - brute_nearest(qa, qb)).max(),
                            np.abs(d2 - brute_nearest(qb, qa)).max())
        if label.startswith("border, cross"):
            for rng in (0, 1):
                s12, s21 = directed_hausdorff(qa, qb, rng=rng)[0], directed_hausdorff(qb, qa, rng=rng)[0]
                checks["scipy"] = max(checks["scipy"], abs(s12 - d1.max()), abs(s21 - d2.max()))
        lin = max(np.percentile(d1, 95), np.percentile(d2, 95))
        icdf = max(np.percentile(d1, 95, method="inverted_cdf"), np.percentile(d2, 95, method="inverted_cdf"))
        checks["gap_ok"] &= gap_bound_ok(d1) and gap_bound_ok(d2)
        print(f"{label:34s}{max(d1.max(), d2.max()):8.3f}{lin:10.3f}{icdf:10.3f}{'-':>10s}"
              f"{max(d1.mean(), d2.mean()):8.3f}")
    # DeepMind corner-grid contour, contour-length weights
    sd = dm_surface_distances(A, B)
    g, p = sd["distances_gt_to_pred"], sd["distances_pred_to_gt"]
    unit = dict(sd, surfel_areas_gt=np.ones_like(g), surfel_areas_pred=np.ones_like(p))
    icdf = max(np.percentile(g, 95, method="inverted_cdf"), np.percentile(p, 95, method="inverted_cdf"))
    checks["dm_unit"] = max(checks["dm_unit"], abs(dm_robust_hausdorff(unit, 95) - icdf))
    checks["gap_ok"] &= gap_bound_ok(g) and gap_bound_ok(p)
    lin = max(np.percentile(g, 95), np.percentile(p, 95))
    print(f"{'DeepMind corner-grid contour':34s}{dm_robust_hausdorff(sd, 100):8.3f}{lin:10.3f}{icdf:10.3f}"
          f"{dm_robust_hausdorff(sd, 95):10.3f}{max(dm_average_surface_distance(sd)):8.3f}")

sq = np.zeros((20, 20), bool); sq[5:15, 5:15] = True
sd_sq = dm_surface_distances(sq, sq)
sd_c = dm_surface_distances(pairs["P1  disc r=40 vs the same disc with a centred hole r=15"][0],
                            pairs["P1  disc r=40 vs the same disc with a centred hole r=15"][0])
print(f"\nchecks: SciPy directed_hausdorff vs brute force {checks['scipy']:.1e}; scikit-image vs brute force "
      f"{checks['skimage']:.1e}; EDT lookup vs brute force {checks['edt']:.1e}")
print(f"checks: DeepMind rule with unit weights vs inverted_cdf {checks['dm_unit']:.1e}; "
      f"|linear - inverted_cdf| within the order-statistic gap: {checks['gap_ok']}")
print(f"checks: DeepMind contour length, 10 x 10 square {sd_sq['surfel_areas_gt'].sum():.6f} "
      f"(4s - 4 + 2 sqrt 2 = {36 + 2 * np.sqrt(2):.6f}); disc r=40 {sd_c['surfel_areas_gt'].sum():.3f} "
      f"(2 pi r = {2 * np.pi * 40:.3f})")

Output:

Part A: Hausdorff conventions on rasterized discs (pixel units)

P1  disc r=40 vs the same disc with a centred hole r=15   [pixels: 5025 and 4316; B inside A: True]
convention                              HD  HD95 lin HD95 icdf  HD95 len     MHD
analytic, filled sets               15.000         -         -         -   0.703
analytic, boundary curves           25.000         -         -    25.000   6.818
scikit-image, all nonzero pixels    15.033         -         -         -   0.749
border, cross erosion (8-conn.)     24.166    24.042    24.042         -   6.706
border, 3x3 erosion (4-conn.)       24.000    23.707    23.707         -   6.694
DeepMind corner-grid contour        25.000    25.000    25.000    25.000   6.772

P2  disc r=40 vs disc r=35, centres 6 px apart   [pixels: 5025 and 3853; B inside A: False]
convention                              HD  HD95 lin HD95 icdf  HD95 len     MHD
analytic, filled sets               11.000         -         -         -   1.046
analytic, boundary curves           11.000         -         -    10.936   5.454
scikit-image, all nonzero pixels    11.000         -         -         -   1.118
border, cross erosion (8-conn.)     11.000    10.770    10.770         -   5.340
border, 3x3 erosion (4-conn.)       11.000    10.440    10.440         -   5.129
DeepMind corner-grid contour        11.000    10.440    10.440    10.770   5.147

P3  disc r=6 vs disc r=4, centres 2 px apart   [pixels: 113 and 49; B inside A: True]
convention                              HD  HD95 lin HD95 icdf  HD95 len     MHD
analytic, filled sets                4.000         -         -         -   0.924
analytic, boundary curves            4.000         -         -     3.982   2.168
scikit-image, all nonzero pixels     4.123         -         -         -   1.179
border, cross erosion (8-conn.)      4.123     4.055     4.123         -   2.114
border, 3x3 erosion (4-conn.)        4.123     3.941     4.000         -   2.014
DeepMind corner-grid contour         4.123     4.000     4.000     4.000   1.937

checks: SciPy directed_hausdorff vs brute force 0.0e+00; scikit-image vs brute force 0.0e+00; EDT lookup vs brute force 0.0e+00
checks: DeepMind rule with unit weights vs inverted_cdf 0.0e+00; |linear - inverted_cdf| within the order-statistic gap: True
checks: DeepMind contour length, 10 x 10 square 38.828427 (4s - 4 + 2 sqrt 2 = 38.828427); disc r=40 266.593 (2 pi r = 251.327)

Every check held, printing zero or True: SciPy’s early-break search, scikit-image’s k-d trees and the distance-transform lookups agree exactly with brute force, the DeepMind rule with unit weights is NumPy’s inverse CDF, and the port gives the exact square perimeter. The containment flags (True, False, True) are as constructed, the pixel counts are within a few pixels of the discs’ areas, and two runs gave identical output. The port’s contour around the radius-40 disc is about 6% longer than the circle (266.593 against 251.327); since the weights are normalized, a uniform excess would not move the percentile.

Filled sets against boundaries (P1). scikit-image reports 15.033, the radius of the hole: the centre pixel’s nearest annulus pixel is at $\sqrt{226}$. Every boundary convention reports 24.0 to 25.0, the distance from the hole’s edge to the outer edge, and the MHD rises from 0.749 on filled sets to about 6.7 on boundaries, because on filled sets every pixel of the annulus and most of the disc lie inside the other set and contribute zero. The prediction (H1) that scikit-image’s HD would lie in $[15.0, 15.1]$ with MHD at most 1.0, and every boundary HD at least 22 with MHD at least 6.0, held. It was written with a planning computation of the same effect on a larger disc in view, so this part illustrates the construction rather than testing it. The filled-set MHD was smaller than every boundary MHD on all three pairs (H3 held: 0.749, 1.118 and 1.179 against at least 6.694, 5.129 and 1.937).

Border conventions (H2). The erosion borders lie inside each object, within one pixel (cross) or $\sqrt2$ (square) of its circle, and the hole’s border lies outside the hole; both offsets narrow the gap between the hole’s edge and the outer edge, so the erosion rows fall below 25 (24.166 and 24.000). DeepMind’s contour runs through points between inside and outside pixel centres, within $\sqrt2/2$ of the circle, and gives 25.000. On P2 every convention gives exactly 11.000, the analytic value and, for the filled sets and the erosion borders, the lower bound that the construction fixes on the line of the centres: the pixel 46 from the smaller disc’s centre belongs to the large disc and to both of its borders, and every pixel of the small disc lies within 35 of its centre, so that pixel is at least 11 from the other set. Lattice distances are square roots of integers, so the printed 11.000 is 11 itself. On P3 every convention gives $\sqrt{17} = 4.123$ against the analytic 4, an error of 3% at this size. The prediction that each value would lie in an interval derived from these bounds (for example $[23.0, 25.1]$ for the cross border on P1 and $[11.0, 12.1]$ on P2) held for all eleven values (the three boundary conventions on P1, and those three and scikit-image on P2 and P3), three of them at the lower limit 11.0.

Percentile rules (H4, H5). On the large shapes, linear interpolation and the inverse CDF gave the same HD95 to three decimals for every surface; in our reading, ties among the lattice distances make the order statistics on both sides of the 95% rank equal. The weighting mattered more than the rule: on P2, DeepMind’s length weights move its HD95 from 10.440 ($\sqrt{109}$) with unit weights to 10.770 ($\sqrt{116}$), a jump between two lattice distances. The prediction that on P1 and P2 the rules for each surface would differ by at most 0.5 (H4a) held: the largest spread was 0.330. The prediction that the largest spread on the small P3 would exceed the largest on both P1 and P2 (H4b) failed: the spreads on P3 were 0.068 (cross border, 4.055 against 4.123), 0.059 (square border) and 0 (DeepMind), more than P1’s 0 but less than P2’s 0.330. On P3 linear interpolation produced values between lattice distances (4.055 and 3.941) that no pair of pixels attains. On P1 every HD95 stayed at least 22 (H5 held: 23.707 to 25.000): the hole’s edge is more than a quarter of the annulus boundary, so the 95th percentile keeps it. None of the HD95 values reaches the continuous HD95 of P2 (10.936) more closely than about 0.17; all lie below it.

Part B matches a template’s edge points under translation, in a scene of our own. The template is the border of a filled arrow, 131 edge points. The image is $192 \times 192$ with the template’s edges at translation (61, 83). The partial distance $h_K$ of the shifted template to the image edges is computed for all 21,744 translations from one distance transform, for $K = 100\%$ (the directed Hausdorff distance) and $K = 80\%$ (the 105th smallest distance), and the minimizer is taken as the match. The conditions are:

• clean;
• 30 random line segments of clutter;
• the right-most 15% of the template’s edge points removed from the image, as if occluded (20 points);
• that occlusion with clutter;
• 30% occlusion (39 points) with clutter, more than the 20% that $K = 80\%$ tolerates.

Clutter is drawn for seeds 0 to 4. The table gives the number of runs whose minimizer is within 1 pixel of the truth, the largest localization error, the ranges of $h_K$ at the truth and at the minimum, and the largest number of tied minima. The ranked values are checked against a k-d tree and the $K = 100\%$ value against SciPy.

# ---------------------------------------------------------------------------
# Part B: partial (ranked) directed Hausdorff matching of a template under translation
H = W = 192
T0 = np.array([61, 83])                                   # true translation of the template (row, col)
tm = np.zeros((50, 43), bool)
rr, cc = polygon(np.array([1, 1, 21, 21, 48, 21, 21]), np.array([11, 31, 31, 41, 21, 1, 11]), tm.shape)
tm[rr, cc] = True                                          # an arrow pointing down
T = np.argwhere(border_cross(tm))                          # template points (model edge pixels)
n = len(T)
order = np.lexsort((T[:, 0], -T[:, 1]))                    # right-most columns first: an occluder from the right

def make_scene(occluded_frac, clutter_seed):
    m = int(round(occluded_frac * n))
    visible = np.ones(n, bool); visible[order[:m]] = False
    E = np.zeros((H, W), bool)
    E[T[visible, 0] + T0[0], T[visible, 1] + T0[1]] = True
    if clutter_seed is not None:
        rng = np.random.default_rng(clutter_seed)
        for _ in range(30):                                 # 30 random straight edge segments
            y0, x0 = rng.uniform(0, [H, W]); a = rng.uniform(0, np.pi); L = rng.uniform(10, 40)
            y1, x1 = np.clip([y0 + L * np.sin(a), x0 + L * np.cos(a)], 0, [H - 1, W - 1])
            r, c = line(int(y0), int(x0), int(y1), int(x1)); E[r, c] = True
    return E, visible, m

def ranked_map(D, frac):
    """h_f(T + t, E) for every translation t: the ceil(frac*n)-th smallest distance from template points to edges."""
    k = int(np.ceil(frac * n)) - 1
    th, tw = T.max(0) + 1
    xs = np.arange(W - tw + 1)
    out = np.empty((H - th + 1, len(xs)))
    for dy in range(H - th + 1):
        v = D[T[:, 0][:, None] + dy, T[:, 1][:, None] + xs[None, :]]
        out[dy] = np.partition(v, k, axis=0)[k]
    return out

conditions = [("clean", 0.0, False), ("clutter", 0.0, True), ("occlusion 15%", 0.15, False),
              ("occlusion 15% + clutter", 0.15, True), ("occlusion 30% + clutter", 0.30, True)]
seeds = [0, 1, 2, 3, 4]
print(f"\nPart B: template of {n} edge points, true translation ({T0[0]}, {T0[1]}), {len(seeds)} clutter seeds "
      f"{seeds}, 30 clutter segments")
print(f"{'condition':26s}{'K':>5s}{'removed':>8s}{'found':>7s}{'max err px':>11s}{'h at truth':>17s}"
      f"{'h at minimum':>17s}{'ties':>6s}")
cb = dict(edt=0.0, scipy=0.0, data=True)
for cname, f, clutter in conditions:
    for frac in (1.0, 0.8):
        errs, h_true, h_min, ties = [], [], [], []
        for s in (seeds if clutter else [None]):
            E, visible, m = make_scene(f, s)
            D = ndi.distance_transform_edt(~E)
            hm = ranked_map(D, frac)
            t = np.array(np.unravel_index(hm.argmin(), hm.shape))
            errs.append(np.linalg.norm(t - T0)); h_true.append(hm[tuple(T0)]); h_min.append(hm.min())
            ties.append(int((hm == hm.min()).sum()))
            E_pts = np.argwhere(E).astype(float)
            k = int(np.ceil(frac * n)) - 1
            for tt in (T0, t):
                bf = np.sort(cKDTree(E_pts).query(T + tt)[0])[k]
                cb["edt"] = max(cb["edt"], abs(bf - hm[tuple(tt)]))
            if frac == 1.0:
                cb["scipy"] = max(cb["scipy"], abs(directed_hausdorff((T + T0).astype(float), E_pts)[0] - hm[tuple(T0)]))
            cb["data"] &= bool(E[T[visible, 0] + T0[0], T[visible, 1] + T0[1]].all()) and m == visible.size - visible.sum()
            if not clutter:
                cb["data"] &= not E[T[~visible, 0] + T0[0], T[~visible, 1] + T0[1]].any()
        found = sum(e <= 1.0 for e in errs)
        rng_s = lambda v: f"{min(v):.3f}-{max(v):.3f}" if len(v) > 1 else f"{v[0]:.3f}"
        print(f"{cname:26s}{int(frac * 100):4d}%{m:8d}{found:4d}/{len(errs):<2d}{max(errs):11.3f}"
              f"{rng_s(h_true):>17s}{rng_s(h_min):>17s}{max(ties):6d}")
print(f"checks: ranked value by EDT vs brute force {cb['edt']:.1e}; K=100% vs SciPy directed_hausdorff "
      f"{cb['scipy']:.1e}; scene data as planned: {cb['data']}")

Output:

Part B: template of 131 edge points, true translation (61, 83), 5 clutter seeds [0, 1, 2, 3, 4], 30 clutter segments
condition                     K removed  found max err px       h at truth     h at minimum  ties
clean                      100%       0   1/1       0.000            0.000            0.000     1
clean                       80%       0   1/1       0.000            0.000            0.000     1
clutter                    100%       0   5/5       0.000      0.000-0.000      0.000-0.000     1
clutter                     80%       0   5/5       0.000      0.000-0.000      0.000-0.000     1
occlusion 15%              100%      20   0/1       5.000           10.050            6.000     6
occlusion 15%               80%      20   1/1       0.000            0.000            0.000     1
occlusion 15% + clutter    100%      20   0/5       5.099     6.083-10.050      3.162-6.000     6
occlusion 15% + clutter     80%      20   5/5       0.000      0.000-0.000      0.000-0.000     1
occlusion 30% + clutter    100%      39   0/5       5.385     6.083-10.050      5.099-8.246    19
occlusion 30% + clutter     80%      39   1/5       3.000      2.000-3.606      1.414-3.000     4
checks: ranked value by EDT vs brute force 0.0e+00; K=100% vs SciPy directed_hausdorff 0.0e+00; scene data as planned: True

The three checks in the last line held. In the clean and cluttered scenes both forms find the template exactly with value 0, as predicted from the construction (HB1 held): at the truth every template point lies on an edge pixel, and clutter adds edges, which can only lower forward distances. This prediction, too, was written with a planning run of a clean match in view and is an illustration.

Occlusion separates the two forms. With 15% of the outline missing and no clutter, the directed Hausdorff distance at the truth is 10.050, set, in our reading, by the occluded right-hand corner of the arrowhead, while six other translations reach 6.000. Moving the template toward the visible edges trades the occluded points’ distance against the visible points’ displacement, so the minimizer of the maximum lies 5 pixels from the truth (the printed error is that of the first of the six tied minima, none of which is the truth, since $10.050 > 6.000$). The ranked form with $K = 80\%$ still has value 0 at the truth, because 111 of the 131 points, more than 105, are on edges, and finds it exactly. The predictions for this condition (HB2) and for occlusion with clutter (HB3: $K = 80\%$ exact in all five seeds, $K = 100\%$ positive at the truth and more than 1 pixel off in all five) held. Clutter lowers the directed distance at the minimum (3.162 to 6.000) but does not bring it back to the truth.

With 30% occluded, only 92 points remain, fewer than the 105 that $K = 80\%$ ranks, so its value at the truth is positive in every seed (2.000 to 3.606; HB4 held). It then found the template in one of the five seeds, and the largest error was 3 pixels; the directed form found it in none of the five, with errors up to 5.385 pixels. The ranked distance works as long as $K$ does not exceed the number of model points that are actually visible. In our reading, this is the trade-off in choosing $K$: a smaller $K$ tolerates more occlusion but also lets a match succeed on a smaller part of the model.

Related Metrics

• Chamfer distance: the mean of the same nearest-neighbour distances, for point clouds and surfaces.
• ASSD, MASD and the normalized surface distance (surface Dice): averaged and thresholded boundary distances, defined in medical image segmentation.
• Dice coefficient and intersection over union: overlap scores that ignore where errors lie, reported together with a boundary distance.
• The distance transform: the map from which directed distances, chamfer matching and partial Hausdorff matching are read.
• Similarity and distance measures, including the Earth Mover’s Distance, for comparing sets and distributions.