Distance Transforms & Skeletonization

"Every pixel in this shape keeps asking me how far it is from the edge. I am the ridge where the answer stops growing. Some call it the medial axis. I call it the only place with a view of both coasts."

A Deeply Centered Medial Axis

The previous section turned blobs into rows of measurements, but its properties were global per object: one area, one centroid, one eccentricity. Many industrial questions are local: is this weld seam at least 3 mm wide everywhere? Where exactly should the gripper pinch this part? Where do two touching tablets actually meet? All of these reduce to distances from boundaries, and computing them pixel by pixel against every boundary point would be ruinously quadratic. The distance transform computes them all in two sweeps. From it flows this section's second hero: the skeleton, which compresses an elongated object onto its centerline so that length and width become countable.

1. The Distance Transform, From Scratch Intermediate

The distance transform of a foreground set $A$ assigns each foreground pixel its distance to the nearest background pixel:

$$ \mathrm{DT}(p) \;=\; \min_{q \,\notin\, A} \; d(p, q), $$

under one of the grid metrics from Section 6.1: city-block, chessboard, or Euclidean. Background pixels get 0; values climb toward the interior like elevation rising away from a coastline. The classic chamfer algorithm computes an excellent Euclidean approximation in exactly two raster sweeps, by exploiting a simple fact: the distance at a pixel is at most the distance at a neighbor plus the cost of the step from it. Sweeping top-left to bottom-right propagates distances from above and the left; the mirrored return sweep finishes the job. Integer step costs of 3 (edge) and 4 (diagonal), divided by 3 at the end, approximate the true $1$ and $\sqrt{2} \approx 1.333\ldots \times 3/3$ remarkably well, resolving the metric sandwich of Section 6.1's Exercise 6.1.3:

import numpy as np

def chamfer_dt(mask):
    """Two-pass (3,4)-chamfer distance transform of a boolean mask."""
    INF = 10**9
    h, w = mask.shape
    d = np.where(mask, INF, 0).astype(np.int64)
    for y in range(h):                       # forward sweep: causal neighbors
        for x in range(w):
            if d[y, x]:
                if x:          d[y, x] = min(d[y, x], d[y, x-1] + 3)
                if y:
                    d[y, x] = min(d[y, x], d[y-1, x] + 3)
                    if x:      d[y, x] = min(d[y, x], d[y-1, x-1] + 4)
                    if x < w-1: d[y, x] = min(d[y, x], d[y-1, x+1] + 4)
    for y in range(h - 1, -1, -1):           # backward sweep: the mirror image
        for x in range(w - 1, -1, -1):
            if d[y, x]:
                if x < w-1:    d[y, x] = min(d[y, x], d[y, x+1] + 3)
                if y < h-1:
                    d[y, x] = min(d[y, x], d[y+1, x] + 3)
                    if x < w-1: d[y, x] = min(d[y, x], d[y+1, x+1] + 4)
                    if x:      d[y, x] = min(d[y, x], d[y+1, x-1] + 4)
    return d / 3.0

disk = (np.hypot(*np.mgrid[-50:51, -50:51]) <= 40)
print(f"chamfer max = {chamfer_dt(disk).max():.2f}  (true inradius 40)")
# Output:
# chamfer max = 40.33  (true inradius 40)

2. Reading the Map: Thickness, Inradius, Clearance Beginner

A surprising number of measurements are single NumPy expressions on the distance map. The global maximum is the inradius, the radius of the largest inscribed disk, and its location is the deepest point of the object: the natural pick-point for a suction gripper and the most defensible "center" of a non-convex shape (the centroid of a crescent lies outside the crescent; the DT maximum never embarrasses you that way). Along a structure of roughly constant width, twice the DT value at the centerline is the local width. And running the transform on the background measures clearances: how close each free pixel is to the nearest obstacle, the standard costmap ingredient in robot path planning.

import cv2
import numpy as np

mask = cv2.imread("gasket_mask.png", cv2.IMREAD_GRAYSCALE)
dist = cv2.distanceTransform(mask, cv2.DIST_L2, cv2.DIST_MASK_PRECISE)

inradius = float(dist.max())
cy, cx = np.unravel_index(np.argmax(dist), dist.shape)
print(f"largest inscribed disk: radius {inradius:.1f} px at ({cx},{cy})")
print(f"minimum wall thickness: {2 * dist[dist > 0].min():.1f} px")
# Representative output:
# largest inscribed disk: radius 47.3 px at (312,198)
# minimum wall thickness: 2.0 px

3. Separating Touching Objects: Peaks as Seeds Intermediate

Here is the distance transform's most celebrated trick. When two convex objects touch, their merged silhouette has one component but two interior summits in the distance map, one per object, separated by a saddle at the contact neck. Find the summits, and you have one seed per true object; grow the seeds back out (flooding downhill on the negated distance map) and the contact line falls exactly along the saddle. The growing step is the watershed algorithm, treated properly with the rest of classical segmentation in Chapter 11; here we use it as a black box, seeded by DT peaks. Figure 6.5.1 visualizes the elevation picture, and the code that follows is the standard recipe, worth keeping within arm's reach forever:

import numpy as np
from scipy import ndimage as ndi
from skimage.feature import peak_local_max
from skimage.segmentation import watershed

# mask: bool array, touching tablets merged into single components
dist = ndi.distance_transform_edt(mask)

# One peak per tablet: enforce a minimum separation ~ tablet radius.
coords = peak_local_max(dist, min_distance=18, labels=mask)
seeds = np.zeros(dist.shape, dtype=bool)
seeds[tuple(coords.T)] = True
markers, n_seeds = ndi.label(seeds)

labels = watershed(-dist, markers, mask=mask)     # flood downhill from seeds
print(f"components before: {ndi.label(mask)[1]}   after: {labels.max()}")
# Representative output on a tray of 24 tablets with 7 contacts:
# components before: 17   after: 24

The recipe has one tuning parameter that matters, min_distance, and it should come from physics: just under the smallest expected object radius. Set it too low and a single elongated object yields two summits and gets cut in half; too high and genuinely separate small objects share a seed and stay merged. For objects far from convex, DT peaks stop being reliable seeds at all, which is the cue to graduate to the full marker strategies of Chapter 11.

4. Skeletons: The Shape, Minus Its Width Intermediate

Now follow the distance map's ridge. Harry Blum's medial axis (1967) is the set of interior points with two or more nearest boundary points, equivalently the centers of maximal inscribed disks, equivalently the ridge crest of the distance transform. Blum's metaphor is unimprovable: set fire to the entire boundary at once and let the flame front advance inward at uniform speed; the medial axis is where flame fronts meet. For an elongated object (wire, crack, vessel, road) the axis is a one-pixel-wide centerline that preserves the object's topology, and paired with the distance values along it, the axis carries the object's complete width profile: shape equals skeleton plus radius function, almost exactly.

In practice two related operations coexist, and scikit-image ships both. skimage.morphology.medial_axis(mask, return_distance=True) computes the ridge definition and hands back the distance map with it. skimage.morphology.skeletonize(mask) instead runs morphological thinning (the Zhang-Suen style iterative peeling built from Section 6.3's hit-or-miss templates, removing boundary pixels whose deletion cannot break connectivity) and generally produces a smoother, less whiskery centerline. The whiskers deserve their warning label: the medial axis is exquisitely sensitive to boundary noise, since every small bump on the contour spawns a spur branch. A light opening before skeletonizing, or pruning of short branches after, is standard hygiene.

import numpy as np
from skimage.morphology import medial_axis, skeletonize

# crack_mask: bool mask of a crack traced on a concrete surface
skel, dist = medial_axis(crack_mask, return_distance=True)
skel2 = skeletonize(crack_mask)            # thinning variant: fewer spurs

length_px = int(skel2.sum())               # centerline pixel count
widths = 2.0 * dist[skel]                  # local width along the axis
print(f"crack length ~ {length_px} px")
print(f"width: mean {widths.mean():.1f} px, "
      f"max {widths.max():.1f} px at the widest point")
# Representative output:
# crack length ~ 1841 px
# width: mean 4.6 px, max 11.2 px at the widest point