Distance Transform

A distance transform turns a binary image into a distance map: every pixel receives its distance to the nearest pixel of a chosen set, the features, which may be edge points, the background, or a single seed. The result is a grey-level image whose ridges run along the middle of shapes and whose level sets are the shapes grown or shrunk by a fixed amount. Distance transforms supply the relief that a watershed floods to split touching objects, the lookup table of chamfer matching, the thresholds that implement erosion and dilation by a disc in mathematical morphology, the medial axis of a shape, and the inside and outside distances from which a signed distance function is built on a grid.

This article covers the transforms of binary images: the two-pass chamfer transforms, which approximate the Euclidean distance with a small mask of integer steps, and exact Euclidean transforms. The generalized distance transform, in which every grid point carries a cost rather than being a feature or not, belongs to pictorial structures and deformable part models, whose articles give its lower-envelope algorithm and check it against brute force; it is summarized under Variants. The Chamfer distance between two point sets is a different object, a metric for comparing shapes, although it takes its name from the chamfering described here.

Problem

Let $F$ be the set of feature pixels of a digital image on a grid $G$, and $d$ a distance between grid points. The distance transform is

$$D_F(p) = \min_{q \in F} d(p, q), \qquad p \in G .$$

Computed directly, it compares every pixel with every feature, $O(|G|\,|F|)$ operations, which is impractical for a megapixel image with tens of thousands of edge points. Distance transform algorithms compute the same map in time proportional to the number of pixels, either exactly or by replacing $d$ with a metric that can be propagated through local neighbourhoods.

Which pixels are the features is a convention that differs between sources and libraries. Barrow et al. (1977) start from an array that is 0 on feature points and infinity elsewhere. SciPy’s distance_transform_edt and distance_transform_cdt give each nonzero (foreground) element its distance to the nearest zero (background) element, so the features are the zeros of the input (SciPy 1.16.3 docstrings). To measure distances from edge points with SciPy, the edge map is inverted first.

Inputs and Outputs

• Input: a binary image (any number of dimensions in SciPy), and the choice of metric: Euclidean, taxicab (city block, the $L_1$ distance), chessboard (the $L_\infty$ distance), or a chamfer mask. SciPy’s Euclidean transform also takes a sampling argument, the spacing of the grid along each axis, for anisotropic pixels or voxels; its chamfer transform has none.
• Output: the distance map, the same size as the input. SciPy returns float64 for the Euclidean transform and int32 for distance_transform_cdt.
• Optional output, the feature transform: for every pixel, the coordinates of a nearest feature (return_indices=True in SciPy). It partitions the grid into the discrete Voronoi cells of the features, and labelling each pixel with the label of its nearest feature is how scikit-image’s expand_labels grows labelled regions without letting them overlap (scikit-image 0.25.2 source). When two features are equally near, which one is returned depends on the implementation; the expand_labels docstring notes this.

Intuition

Distances obey a local rule. If a pixel’s nearest feature lies at distance $D$, then a neighbour one step away is at most $D$ plus the length of that step from the same feature. Starting from 0 on the features and “infinity” elsewhere, repeatedly replacing each pixel’s value by the minimum of its own value and its neighbours’ values plus the step lengths spreads the distances outward like a front, until nothing changes.

Sequential (raster) scans make this cheap. A forward pass from the top-left corner, row by row, lets each pixel take the minimum over its already-visited neighbours above and to the left, so distances propagate down and to the right. A backward pass from the bottom-right corner propagates them up and to the left. Barrow et al. (1977) compute their distances in two passes with exactly these two loops, and for every mask in the example below two passes reach the exact path metric. The distance produced is then the length of the shortest path through the grid made of the mask’s steps, not the straight-line distance, and the mask’s weights decide how far the two differ.

The Euclidean distance cannot be reproduced exactly by any fixed set of step lengths: a path metric built from finitely many steps has polygonal level sets, never circles (our reasoning; the error bounds below quantify the gap). Exact Euclidean algorithms instead use the fact that the squared distance is a sum over axes, and solve a one-dimensional problem along each axis in turn.

Algorithm

Two-pass chamfer transform (a mask of steps $(\Delta y, \Delta x)$ with integer weights $w$, split into a forward half, the steps to pixels visited earlier in raster order, and the mirrored backward half):

1. Set $D(p) = 0$ for features and $D(p) = \infty$ (a large integer) elsewhere. Pad the image so that mask steps never leave the array; padded pixels keep the large value.
2. Forward pass, rows top to bottom, each row left to right: $D(p) \leftarrow \min\big(D(p),\ \min_{(\Delta, w) \in \text{forward}} D(p + \Delta) + w\big)$.
3. Backward pass, rows bottom to top, each row right to left: the same update with the mirrored steps $-\Delta$.
4. Divide by the axial weight to express the result in pixels.

For a $3 \times 3$ mask the forward half holds the left, upper-left, upper and upper-right neighbours. Barrow et al. (1977) use increments 2 and 3, which they describe as approximating distances 1 and $\sqrt 2$; the 3-4 mask is another integer pair. A $5 \times 5$ mask adds the eight “knight’s move” steps $(\pm 1, \pm 2)$ and $(\pm 2, \pm 1)$; the 5-7-11 mask gives axial steps weight 5, diagonal steps 7 and knight steps 11. Weights 1 and 2 give the taxicab distance (a diagonal step costs two axial ones), and weights 1 and 1 the chessboard distance.

SciPy’s distance_transform_cdt runs this scheme with unit steps only. Its C routine makes one pass over the array, reverses the array and makes a second pass, each time taking the minimum over the first half of a binary $3 \times \dots \times 3$ structure plus 1; “taxicab” uses the structure with face neighbours only and “chessboard”, the default, the full $3 \times 3$ neighbourhood (SciPy 1.16.3 source). It therefore computes taxicab, chessboard or another unit-step metric, but not weighted masks such as 3-4 or 5-7-11, which the example implements in Python. Its neighbours outside the array are skipped, so the image border acts neither as feature nor as obstacle.

Exact Euclidean transform. SciPy’s distance_transform_edt computes the feature transform first, with the algorithm of Maurer, Qi and Raghavan (2003), to which a comment in the C source attributes it, and then the distances as the length of the vector from each pixel to its feature (SciPy 1.16.3 source). The Maurer et al. paper could not be read for this article, and what follows is read from the source. The routine works recursively over dimensions: it finds, along each line of the first axis, the nearest feature in that line, and then, along each line of each further axis, keeps only the candidate features whose Voronoi regions cross the line, discarding a candidate when the two on either side of it hide it, and reads off the nearest one at every position. Because the distances are computed from integer offsets, $\sqrt{\Delta y^2 + \Delta x^2}$, the result equals brute force exactly, as the example’s first check prints.

Felzenszwalb and Huttenlocher’s lower-envelope algorithm reaches the same goal along a different route: for the squared distance in one dimension, the transform is the lower envelope of parabolas rooted at the grid points, built in one pass, and a two-dimensional transform is a pass along every column followed by one along every row. The deformable part models article derives the intersection formula and implements it; the pictorial structures article attributes the lower-envelope method to “Distance Transforms of Sampled Functions” (Theory of Computing, 2012), which was not read for this article.

Mathematical Formulation

Chamfer distance to one feature. For a $3 \times 3$ mask with axial weight $a$ and diagonal weight $b$, $a \le b \le 2a$, take an offset $(x, y)$ with $x \ge y \ge 0$ (the other octants follow by symmetry). The cheapest path takes $y$ diagonal steps and $x - y$ axial ones, so the path length, divided by $a$, is

$$d_{a,b}(x, y) = x + \left(\tfrac{b}{a} - 1\right) y .$$

For the 5-7-11 mask, divided by 5, the cheapest path uses knight and axial steps below the direction $(2, 1)$, and knight and diagonal steps above it:

$$d_{5,7,11}(x, y) = \begin{cases} x + 0.2\,y, & 2y \le x, \\ 0.8\,x + 0.6\,y, & 2y > x . \end{cases}$$

These closed forms are our derivation; in the example, the two-pass transforms agree with them on every pixel of both test images (the “vs exact” column).

Error bounds (our derivation). Borgefors (1986) could not be read for this article, so the bounds below are derived here, with our normalization: each chamfer distance is divided by its axial weight, and the error is the relative error $e = (d_{\text{mask}} - d_E)/d_E$ against the Euclidean distance $d_E$. With $t = y/x \in [0, 1]$, the $3 \times 3$ mask gives

$$1 + e(t) = \frac{1 + c\,t}{\sqrt{1 + t^2}}, \qquad c = \frac{b}{a} - 1 .$$

Its derivative vanishes at $t = c$, where $1 + e = \sqrt{1 + c^2}$, the maximum; the minimum is at an end of the interval, $1 + e(0) = 1$ or $1 + e(1) = (1 + c)/\sqrt 2$. The same steps applied to each branch of the 5-7-11 formula give its band, with the extremes at $t = 1/5$ and $t = 1/2$; on the upper branch, the stationary point $t = 3/4$ has zero error.

Mask (divided by)	$c$	Lowest error	Highest error	Directions $y/x$ of the extremes
chessboard 1-1	0	$1/\sqrt 2 - 1 = -29.29\%$	0	1; 0
3-4 (3)	1/3	$\tfrac{4}{3\sqrt 2} - 1 = -5.72\%$	$\sqrt{10/9} - 1 = +5.41\%$	1; 1/3
2-3 (2), Barrow et al.	1/2	0	$\sqrt{5/4} - 1 = +11.80\%$	0; 1/2
taxicab 1-2	1	0	$\sqrt 2 - 1 = +41.42\%$	0; 1
5-7-11 (5)	—	$\tfrac{1.1}{\sqrt{1.25}} - 1 = -1.61\%$	$\sqrt{1.04} - 1 = +1.98\%$	1/2; 1/5

The 2-3 mask, whose $(1 + c)/\sqrt 2 = 1.061$ exceeds 1, overestimates in every direction except the axes. The 3-4 mask over- and underestimates by similar amounts, and the 5-7-11 mask narrows the band about threefold.

From one feature to many (our derivation). The chamfer transform of a feature set is $D_F(p) = \min_q d_{\text{mask}}(p, q)$. If $q^*$ is the Euclidean-nearest feature, $D_F(p) \le d_{\text{mask}}(p, q^*) \le (1 + e_{\max})\, d_E(p)$; if $q'$ is the feature nearest under the mask, $D_F(p) = d_{\text{mask}}(p, q') \ge (1 + e_{\min})\, d_E(p, q') \ge (1 + e_{\min})\, d_E(p)$. Every pixel of every image therefore lies within the one-feature band, although a particular image need not reach its ends. The absolute error is the relative error times the distance, so its bound grows in proportion to the distance from the features.

Parameters

• Metric or mask. Chessboard and taxicab are exact for their own metrics, but their unit discs are squares and diamonds, and their errors against the Euclidean distance reach 29% and 41%. Integer chamfer masks are a compromise between cost and accuracy: a larger mask has more steps per pixel and a narrower error band. An exact Euclidean transform removes the approximation altogether.
• Normalization. A chamfer result is in units of the axial weight; dividing by it expresses distances in pixels. Comparisons with a threshold can stay in integers, as in the example ($D \le 3R$ for the 3-4 mask).
• Grid spacing. For pixels or voxels that are not square, SciPy’s Euclidean transform scales each axis by its sampling entry; the medical image segmentation and organ segmentation articles use it for surface distances on voxel grids with unequal spacing.
• Feature set. What counts as a feature decides the meaning of the map: background pixels give each object pixel its depth inside the object; edge pixels give the distance to the nearest edge, as in chamfer matching; a single seed gives a distance from that point.

Complexity

The two-pass chamfer transform visits each pixel twice and looks at half the mask each time, $O(N k)$ for $N$ pixels and a mask of $k$ steps (our count from the loops), with memory for one integer array. The Maurer et al. algorithm runs in time linear in the number of pixels, as its title states, and SciPy stores an int32 feature transform with one array per dimension before computing float64 distances from it (SciPy 1.16.3 source). The lower-envelope algorithm is also linear, as the deformable part models article explains. SciPy’s distance_transform_bf is a brute-force routine: it compares every foreground element with every background element that touches the foreground (8-connected in two dimensions), which its docstring calls slow (SciPy 1.16.3 docstring and source).

Implementation

The NumPy and SciPy code compares five transforms with SciPy’s exact Euclidean transform on two binary images with known geometry. Scene A is a $121 \times 121$ image with a single feature pixel in the centre, so every pixel’s nearest feature is known and every direction up to the image edge is represented. Scene B is a $160 \times 160$ image with three shapes, a disc of radius 30, a rectangle with half-sides 42 and 14 rotated by 35° and a triangle, at least 12 pixels from the border; distances are measured inside the shapes to the nearest background pixel. The methods are a hand-written two-pass chamfer transform with Barrow et al.’s 2-3 mask, the 3-4 mask and the 5-7-11 mask, and SciPy’s distance_transform_cdt with the taxicab and chessboard metrics.

The script checks every transform against its definition, computed by brute force over all background pixels: the Euclidean transform against $\min_q \lVert p - q \rVert$, and each chamfer transform against the minimum of its closed-form path length from the section above (“vs exact” counts the pixels where they differ). It checks the feature transform, the agreement of the hand-written transform with SciPy’s on the unit masks, and that thresholding the exact distance map reproduces dilation (scene A, radius 40) and erosion (scene B, radius 6) by a disc footprint, the set of offsets with $\Delta y^2 + \Delta x^2 \le R^2$. For each method it prints the lowest and highest relative error in percent, the signed margins to the derived band (min − lo and hi − max; a negative margin would mean the error left the band), the numbers of under- and overestimating pixels, the largest absolute error in pixels, and the number of pixels on which thresholding the method’s map disagrees with the exact disc morphology (“disc diff”).

It ran with NumPy 2.3.5 and SciPy 1.16.3 in about 5 to 6 seconds, single-threaded; two runs printed identical output. Nothing in it is random.

import os
for v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[v] = "1"                                   # one BLAS thread, set before NumPy loads
import numpy as np
from scipy import ndimage as ndi

def chamfer(fg, mask):
    """Two-pass chamfer transform of the foreground pixels of a boolean image.
    mask: forward-half offsets (dy, dx, w) with integer weights; the backward pass uses
    the mirrored offsets. Returns integer distances (not divided by the axial weight)."""
    H, W = fg.shape
    p = max(max(abs(dy), abs(dx)) for dy, dx, _ in mask)
    BIG = 10 ** 9
    D = np.full((H + 2 * p, W + 2 * p), BIG, dtype=np.int64)    # padding never wins a minimum
    D[p:p + H, p:p + W] = np.where(fg, BIG, 0)
    D = D.tolist()
    for y in range(p, p + H):                                   # forward: top-left to bottom-right
        for x in range(p, p + W):
            if D[y][x]:
                D[y][x] = min([D[y][x]] + [D[y + dy][x + dx] + w for dy, dx, w in mask])
    for y in range(p + H - 1, p - 1, -1):                       # backward: bottom-right to top-left
        for x in range(p + W - 1, p - 1, -1):
            if D[y][x]:
                D[y][x] = min([D[y][x]] + [D[y - dy][x - dx] + w for dy, dx, w in mask])
    return np.array(D, dtype=np.int64)[p:p + H, p:p + W]

M3 = lambda a, b: [(0, -1, a), (-1, -1, b), (-1, 0, a), (-1, 1, b)]       # 3 x 3, forward half
M5711 = M3(5, 7) + [(-1, -2, 11), (-2, -1, 11), (-2, 1, 11), (-1, 2, 11)]  # 5 x 5, forward half

def path_metric(dy, dx, name):
    """Exact length of the cheapest mask path for an offset (the chamfer distance to one point)."""
    X = np.maximum(abs(dy), abs(dx)); Y = np.minimum(abs(dy), abs(dx))
    a, b, c = {"2-3": (2, 3, 0), "3-4": (3, 4, 0), "5-7-11": (5, 7, 11), "taxicab": (1, 2, 0),
               "chessboard": (1, 1, 0)}[name]
    if c:
        return np.where(2 * Y <= X, a * X + (c - 2 * a) * Y, (c - b) * X + (2 * b - c) * Y)
    return a * X + (b - a) * Y

def brute(fg, f):
    """min over background pixels q of f(p - q) at every foreground pixel p: the definition."""
    by, bx = np.nonzero(~fg); fy, fx = np.nonzero(fg)
    out = np.zeros(fg.shape)
    for i in range(0, len(fy), 512):
        out[fy[i:i + 512], fx[i:i + 512]] = f(fy[i:i + 512, None] - by, fx[i:i + 512, None] - bx).min(1)
    return out

# Scene A: a single background (feature) pixel. Scene B: three shapes, distances measured inside them.
fgA = np.ones((121, 121), bool); fgA[60, 60] = False
r, c = np.mgrid[0:160, 0:160].astype(float)
disc = (r - 48) ** 2 + (c - 46) ** 2 <= 30 ** 2
th = np.deg2rad(35)
rect = (abs((c - 104) * np.cos(th) + (r - 108) * np.sin(th)) <= 42) & \
       (abs(-(c - 104) * np.sin(th) + (r - 108) * np.cos(th)) <= 14)
side = lambda p, q, s: np.sign((q[1] - p[1]) * (r - p[0]) - (q[0] - p[0]) * (c - p[1])) == \
                       np.sign((q[1] - p[1]) * (s[0] - p[0]) - (q[0] - p[0]) * (s[1] - p[1]))
A_, B_, C_ = (14, 100), (14, 148), (60, 124)
fgB = disc | rect | (side(A_, B_, C_) & side(B_, C_, A_) & side(C_, A_, B_))
ys, xs = np.nonzero(fgB)
print(f"scene A: {int((~fgA).sum())} background pixel")
print(f"scene B: {ndi.label(fgB, np.ones((3, 3)))[1]} shapes, {fgB.sum()} foreground pixels, "
      f"margin to the border {min(ys.min(), xs.min(), 159 - ys.max(), 159 - xs.max())} px")

# Our derived bounds on the relative error (normalized by the axial weight), and the 3 x 3 / 5 x 5 masks.
METHODS = {  # name: (scale, low bound, high bound)
    "2-3": (2, 0.0, np.sqrt(5 / 4) - 1),
    "3-4": (3, (4 / 3) / np.sqrt(2) - 1, np.sqrt(10 / 9) - 1),
    "5-7-11": (5, 1.1 / np.sqrt(1.25) - 1, np.sqrt(1.04) - 1),
    "taxicab": (1, 0.0, np.sqrt(2) - 1),
    "chessboard": (1, 1 / np.sqrt(2) - 1, 0.0),
}
disc_fp = lambda R: (lambda g: g[0] ** 2 + g[1] ** 2 <= R * R)(np.mgrid[-R:R + 1, -R:R + 1])

for name_s, fg, op, R in [("A", fgA, "dilation", 40), ("B", fgB, "erosion", 6)]:
    dE, (iy, ix) = ndi.distance_transform_edt(fg, return_indices=True)
    yy, xx = np.indices(fg.shape)
    print(f"\nscene {name_s} {fg.shape}: EDT vs brute force max |diff| "
          f"{np.abs(dE - brute(fg, lambda dy, dx: np.sqrt(dy ** 2 + dx ** 2))).max():.1e}; "
          f"distance to returned feature max |diff| {np.abs(np.hypot(yy - iy, xx - ix) - dE).max():.1e}, "
          f"features on foreground {int(fg[iy, ix].sum())}")
    if op == "dilation":   # dilate the feature pixel by a disc / erode the shapes by a disc
        exact, thr_E = ndi.binary_dilation(~fg, structure=disc_fp(R)), dE <= R
    else:
        exact, thr_E = ndi.binary_erosion(fg, structure=disc_fp(R)), dE > R
    print(f"{op} by a disc of radius {R}: footprint vs EDT threshold, {int((exact != thr_E).sum())} pixels differ")
    D = {"2-3": chamfer(fg, M3(2, 3)), "3-4": chamfer(fg, M3(3, 4)), "5-7-11": chamfer(fg, M5711),
         "taxicab": ndi.distance_transform_cdt(fg, metric="taxicab"),
         "chessboard": ndi.distance_transform_cdt(fg, metric="chessboard")}
    print(f"hand two-pass vs SciPy cdt: taxicab {int((chamfer(fg, M3(1, 2)) != D['taxicab']).sum())}, "
          f"chessboard {int((chamfer(fg, M3(1, 1)) != D['chessboard']).sum())} pixels differ")
    print("method      vs exact  min err %  max err %  min-lo   hi-max  under  over  max|err| px  disc diff")
    nz = fg & (dE > 0)
    for m, (s, lo, hi) in METHODS.items():
        exact_m = brute(fg, lambda dy, dx: path_metric(dy, dx, m))
        rel = (D[m][nz] / s - dE[nz]) / dE[nz]
        thr_m = D[m] <= s * R if op == "dilation" else D[m] > s * R   # integer comparison
        print(f"{m:10s} {int((D[m] != exact_m).sum()):6d} {100 * rel.min():+10.4f} {100 * rel.max():+10.4f} "
              f"{rel.min() - lo:+8.1e} {hi - rel.max():+8.1e} {int((rel < 0).sum()):6d} {int((rel > 0).sum()):5d} "
              f"{np.abs(D[m][nz] / s - dE[nz]).max():10.3f} {int((thr_m != exact).sum()):8d}")

Output:

scene A: 1 background pixel
scene B: 3 shapes, 6247 foreground pixels, margin to the border 12 px

scene A (121, 121): EDT vs brute force max |diff| 0.0e+00; distance to returned feature max |diff| 0.0e+00, features on foreground 0
dilation by a disc of radius 40: footprint vs EDT threshold, 0 pixels differ
hand two-pass vs SciPy cdt: taxicab 0, chessboard 0 pixels differ
method      vs exact  min err %  max err %  min-lo   hi-max  under  over  max|err| px  disc diff
2-3             0    +0.0000   +11.8034 +0.0e+00 -1.4e-17      0 14400      8.038      704
3-4             0    -5.7191    +5.4093 -2.1e-17 -3.5e-17   3600 10680      4.853      296
5-7-11          0    -1.6130    +1.9804 -1.8e-16 -3.8e-17   8400  5840      1.212       64
taxicab         0    +0.0000   +41.4214 +0.0e+00 -5.6e-17      0 14400     35.147     1744
chessboard      0   -29.2893    +0.0000 +0.0e+00 +0.0e+00  14400     0     24.853     1536

scene B (160, 160): EDT vs brute force max |diff| 0.0e+00; distance to returned feature max |diff| 0.0e+00, features on foreground 0
erosion by a disc of radius 6: footprint vs EDT threshold, 0 pixels differ
hand two-pass vs SciPy cdt: taxicab 0, chessboard 0 pixels differ
method      vs exact  min err %  max err %  min-lo   hi-max  under  over  max|err| px  disc diff
2-3             0    +0.0000   +11.8034 +0.0e+00 +1.1e-16      0  4765      1.875       87
3-4             0    -5.7191    +5.4093 +6.9e-18 +7.6e-17   2623  1896      1.118        0
5-7-11          0    -1.6130    +1.9804 -1.8e-16 -3.8e-17   3615   912      0.397        0
taxicab         0    +0.0000   +41.4214 +0.0e+00 -5.6e-17      0  4765      3.787      239
chessboard      0   -29.2893    +0.0000 +0.0e+00 +0.0e+00   5191     0      8.017      835

The checks all pass. The data have the planned form, one feature pixel in scene A and three separate shapes 12 pixels from the border in scene B. The Euclidean transform matches brute force exactly on both scenes, and each returned feature is a background pixel at exactly the distance reported. Each of the five transforms equals its own path metric on every pixel, so the two raster passes compute the shortest mask path even with the knight steps of the 5-7-11 mask, and the hand-written unit masks agree with SciPy’s chamfer transform. Thresholding the exact map reproduces dilation and erosion by a disc with no pixel differing, which is how scikit-image’s isotropic_dilation and isotropic_erosion are implemented (scikit-image 0.25.2 source).

On scene A every method reaches both ends of its derived band, to rounding error ($10^{-16}$), because the extreme directions, offsets such as $(3, 1)$, $(2, 1)$ and $(5, 1)$, all occur in the image. On scene B every error stays within its band: the margins are zero or positive apart from values of order $10^{-16}$, the same rounding seen on scene A, so the derived bound held within its tolerance of $10^{-12}$. Scene B also reaches the ends of every band; in our reading, some pixels inside the shapes see their nearest background pixel along those same small offsets. The planned order of the largest absolute relative error, 5-7-11 (1.98%) below 3-4 (5.72%) below 2-3 (11.80%) below chessboard (29.29%) below taxicab (41.42%), held on both scenes. The planned sign behaviour also held on both scenes: the 2-3 mask never underestimates, while 3-4 and 5-7-11 both under- and overestimate; as derived above, taxicab never underestimates and chessboard never overestimates.

The absolute errors follow the distances. On scene A, where distances reach 84.9 pixels at the corners, the largest error is 4.85 pixels for 3-4 and 1.21 for 5-7-11; inside the shapes of scene B, where the largest distances, near the centre of the disc of radius 30, are about 30 pixels, they are 1.12 and 0.40. The morphology counts measure what the errors mean for a disc operation. On scene A the planned order held: dilation by a disc of radius 40 differs from the exact disc on 64 pixels for 5-7-11, 296 for 3-4 and 704 to 1,744 for the others. On scene B, eroding by a disc of radius 6, the planned hypothesis that 5-7-11 would disagree on fewer pixels than 3-4 failed: both disagree on no pixel, a tie at 0. Its other part, that 3-4 disagrees on fewer pixels than each of 2-3, taxicab and chessboard, held.

Properties and Behavior

• A chamfer transform is exact for its own metric. It computes the shortest-path distance on the graph whose edges are the mask’s steps, so it is a metric on the grid, and its errors are the differences between that metric and the Euclidean one, not numerical errors. The example’s “vs exact” column shows this: zero differences on every pixel.
• Anisotropy. The level sets of a chamfer transform around one point are polygons: squares for chessboard, diamonds for taxicab, octagons for $3 \times 3$ masks with $a < b < 2a$, and sixteen-sided figures for 5-7-11 (it follows from the closed forms, which are piecewise linear in the direction). Thresholding a chamfer map therefore grows or shrinks shapes by polygons rather than discs.
• Bounded relative error, growing absolute error. The relative error depends only on direction, so the absolute error grows with the distance, as the example’s two scenes show; a chamfer map is accurate near the features and less so far from them.
• Exact Euclidean maps on a grid. The Euclidean transform gives the distance between pixel centres. When the features are the pixels of a digitized boundary, the result is exact for those pixels, not for the continuous boundary they sample (our reading).
• Morphology by thresholding. Erosion of a shape by a disc of radius $r$ keeps the pixels whose exact distance to the background exceeds $r$, and dilation adds the pixels of the background within $r$ of the shape. Computing the distance map once gives these results for every radius, at a cost that does not grow with $r$, unlike erosion with a footprint of $(2r + 1)^2$ pixels (our reasoning). scikit-image’s isotropic morphology functions are implemented this way, and their docstrings say they are faster for large circular footprints (scikit-image 0.25.2 source).
• Uses of the map. The watershed article floods the negated Euclidean distance transform to separate overlapping disks, and cell segmentation uses the same recipe for touching nuclei. scikit-image’s medial_axis computes the Euclidean distance transform of the shape and can return it with the skeleton (scikit-image 0.25.2 source); the medial axis is the set of points with two or more nearest boundary points, as the signed distance function article puts it for surfaces. Barrow et al. (1977) chamfered the edge map of an image once and then scored each predicted placement of a map by summing the distance values at the predicted feature points, so that a candidate pose costs one lookup per point; the Chamfer distance article continues that history.

Limitations

• Chamfer masks are approximations of the Euclidean distance. Even 5-7-11 is off by up to about 2% in some directions, which matters when a distance is used as a measurement, for example a thickness, rather than as an ordering or a relief.
• The Euclidean distance has no exact local propagation. Exact algorithms need the separable, line-by-line structure of the squared distance, and they keep track of feature coordinates rather than distances alone; propagation with a fixed mask cannot be made exact by choosing weights, since its level sets stay polygons (our reasoning).
• Ties and discontinuities. The feature transform jumps where two features are equally near; along those lines, which become the medial axis of a shape, the nearest-feature label is not unique, and implementations break the tie in their own way.
• Binary input. The transforms here need a binary decision about what is a feature. A thresholded or noisy edge map passes its errors on, since a single spurious feature pixel sets the distance of everything around it to near zero.
• Distances in a non-convex domain. The transforms measure straight-line (or mask-path) distances across the whole grid. Distances measured within a shape, around its holes and bends, need a transform that constrains propagation to the shape (a geodesic distance, in our reading of the term); the methods here do not compute it.

Variants

• Unit-step transforms. Taxicab and chessboard distances, computed by two passes with unit steps, are the simplest. Barrow et al. (1977) cite Munson (an internal memo, Stanford Research Institute, 1973) and Rosenfeld for the process “known as chamfering”; the Rosenfeld entry of their reference list is Rosenfeld and Pfaltz, “Distance Functions on Digital Pictures” (Pattern Recognition, 1968), which was not read for this article.
• Weighted chamfer masks. Integer weights on $3 \times 3$ and larger masks, as in the example. Borgefors (1986) could not be read for this article; Borgefors (1988) used chamfer transforms for hierarchical edge matching, as the Chamfer distance article describes.
• Exact Euclidean transforms. Feature-transform methods such as that of Maurer, Qi and Raghavan (2003), used by SciPy, and the lower-envelope method of Felzenszwalb and Huttenlocher, which works on squared distances and extends to arbitrary costs.
• Generalized distance transforms. Replacing the binary feature indicator by an arbitrary cost function $f$ gives $\mathcal{D}_f(p) = \min_q \big(d(p, q) + f(q)\big)$, the min-convolution at the heart of efficient inference in pictorial structures and matching in deformable part models. Both articles treat it in full.
• Signed distance transforms. Subtracting the transform of a shape’s inside (each inside pixel’s distance to the outside) from that of its outside gives a signed distance on the grid, negative inside and positive outside. That is the convention of the signed distance function article, which notes that some systems flip the sign and covers signed fields, their eikonal property, and the fast marching and fast sweeping methods that fill a grid from a band of known distances.
• Higher dimensions and anisotropic grids. SciPy’s transforms work in any number of dimensions, and its Euclidean transform takes a grid spacing per axis, which volumetric medical images need.