Watershed

The watershed transform segments an image by reading it as a landscape. Gray values become altitudes; water poured on the surface collects in valleys, the catchment basins, and the ridges that separate neighboring basins are the watershed lines. On the gradient magnitude, basins are regions of nearly uniform intensity and lines follow the strongest edges; on the negated distance transform of a binary mask, basins split touching blobs into separate objects. The method comes from mathematical morphology and is one of the classical image segmentation methods, still used to separate touching cells, nuclei and grains and to build region hierarchies and superpixels.

As Vincent and Soille (1991) recount, the watershed entered image analysis through Digabel and Lantuéjoul, who applied it to stacks of thresholded binary images. Beucher and Lantuéjoul (1979) extended it to grayscale images and defined contours as the watersheds of the gradient modulus. Vincent and Soille (1991) made it fast by simulating immersion with sorted pixels and first-in-first-out (FIFO) queues. Meyer’s flooding algorithm with a hierarchical (priority) queue, published in French in 1992 and presented in English by Beucher and Meyer (1993), is the algorithm OpenCV implements.

Problem

Given an image $f$, partition its pixels into regions, one per significant valley of $f$, with boundaries on the crests between valleys. The relief is chosen so that objects are valleys: the gradient magnitude, the inverted intensity of bright blobs, or the negated distance to the background of a binary mask.

Inputs and Outputs

The input is a grayscale relief $f$, a connectivity (4 or 8 neighbors in 2D), and optionally a marker image of labeled seeds and a mask. The output is a label image: each pixel carries the label of its basin or marker, and some algorithms give pixels on watershed lines a special label.

Intuition

Pierce a hole at the bottom of every valley and lower the landscape slowly into a lake. Water rises through the holes and fills each basin from its lowest point. Where water from two basins would merge, build a dam. When the surface is submerged, the dams are the watershed lines and enclose the catchment basins (Vincent & Soille, 1991). Equivalently, a drop of water runs down the steepest slope to one minimum; a basin is the set of pixels whose drops reach the same minimum.

Flooding from markers pierces holes only in chosen regions. A basin without a marker fills by overflowing from a neighbor at the lowest pass, and dams separate only floods from different markers (Beucher & Meyer, 1993).

Algorithm

Immersion (Vincent & Soille, 1991).

1. Sort the pixels by gray value, in linear time for integer images, with a histogram-based distributive sort.
2. For each level $h$ in increasing order, put pixels at level $h$ that touch a labeled basin in a FIFO queue and extend the basins into level $h$ by breadth-first propagation; pixels equidistant from two basins become watershed pixels.
3. Give each connected group of still-unreached pixels at level $h$, a new minimum, a new label.

Marker-based flooding with a priority queue (Meyer; Beucher & Meyer, 1993).

1. Label the markers; insert their boundary pixels into a priority queue with priority $f(p)$.
2. Repeatedly extract the pixel with the lowest priority, taking the earliest-inserted pixel among equals.
3. Give each unlabeled neighbor the label of the extracted pixel and insert it with priority $f$ of that neighbor.
4. Stop when the queue is empty.

This version produces touching basins without lines; a variant marks pixels where different floods meet. FIFO order among equal priorities splits plateaus by distance to their lower border. With the regional minima of $f$ as markers, the result is the plain watershed of $f$; with other markers, it is the watershed of $f$ modified so that the markers are its only minima.

Mathematical Formulation

By immersion. Let $T_h = \{p : f(p) \le h\}$ be the threshold set at level $h$, $\min_h$ the regional minima at altitude $h$, and $IZ_A(B)$ the geodesic influence zones of the connected components of $B$ inside $A$, the pixels of $A$ closer, by paths inside $A$, to one component than to any other. Vincent and Soille define the basins by the recursion

$$X_{h_{\min}} = T_{h_{\min}}, \qquad X_{h+1} = \min\nolimits_{h+1} \cup IZ_{T_{h+1}}(X_h),$$

and the watershed as the pixels outside $X_{h_{\max}}$.

By topographic distance. Meyer (1994), as Roerdink and Meijster (2000) present it, defined basins through a distance that integrates slope along a path. For a smooth $f$,

$$T_f(p, q) = \inf_{\gamma} \int_\gamma \|\nabla f(\gamma(s))\|\, ds ,$$

over paths $\gamma$ from $p$ to $q$. The basin of a minimum $m_i$ is the set of points topographically closer to it than to any other minimum,

$$CB(m_i) = \{ x : f(m_i) + T_f(x, m_i) < f(m_j) + T_f(x, m_j) \ \ \forall j \ne i \},$$

and the watershed is the complement of all basins. Paths of steepest descent are the shortest paths of $T_f$, which formalizes the drop-of-water picture; the digital version uses a lower slope per pixel and a lower completion of plateaus (Roerdink & Meijster, 2000).

The definitions do not agree exactly. Roerdink and Meijster gave small images where immersion and topographic distance disagree, where results depend on connectivity, where watersheds are thick, and where ordered queues place plateau lines by processing order; they also found that the Vincent–Soille algorithm does not adhere exactly to its own definition.

Parameters

• The relief $f$. The most consequential choice: gradient magnitude for contrast boundaries, the negated distance transform for shape-based splitting.
• Pre-smoothing. A Gaussian filter or a morphological filter before the gradient removes some noise minima.
• Markers. Their number sets the number of regions; for distance-transform maxima, the window size and minimum distance decide when an elongated blob splits.
• Minima suppression. Filling minima shallower than a height $h$ (the h-minima transform) removes shallow noise basins.
• Connectivity. 4- or 8-neighbors change both the minima and the lines.

Complexity

Vincent and Soille’s algorithm is linear in the number of pixels $N$: two scans to sort and about three visits per pixel during flooding; they reported about 2.5 seconds for a $256 \times 256$ image on a Sun SPARCstation 1. A hierarchical queue with one FIFO per gray level processes each pixel once (Beucher & Meyer, 1993), also linear for integer images; a binary heap, as below, costs $O(N \log N)$. Memory is $O(N)$.

Implementation

The NumPy and SciPy code builds a $96 \times 96$ image of four bright disks, three overlapping, with Gaussian noise ($\sigma = 0.15$) and known instance masks. The watershed function is a marker-based priority flood: markers enter the queue first, so they act as the only minima, and a counter makes ties first in, first out. The script then (1) thresholds and labels connected components, (2) splits the mask by flooding the negated Euclidean distance transform from its local maxima, (3) floods the gradient magnitude from all its regional minima, at three smoothing scales, and (4) floods the gradient from one marker per object plus a background marker. Each object is scored by its best-matching region’s IoU.

import heapq
import numpy as np
from scipy import ndimage as ndi

SEED = 0
rng = np.random.default_rng(SEED)

# A 96 x 96 image of four bright disks, three of them overlapping, on a dark
# background with Gaussian noise. Ground truth: each pixel belongs to the disk
# whose center is nearest among the disks that contain it.
H = W = 96
yy, xx = np.mgrid[0:H, 0:W]
disks = [(30, 28, 14), (30, 52, 13), (52, 40, 12), (70, 72, 13)]   # (row, col, radius)
d2 = np.stack([(yy - r)**2 + (xx - c)**2 for r, c, _ in disks])
inside = d2 <= np.array([rad for *_, rad in disks])[:, None, None]**2
gt = np.where(inside.any(0), np.argmin(np.where(inside, d2, np.inf), axis=0) + 1, 0)
img = np.where(gt > 0, 0.8, 0.2) + rng.normal(0, 0.15, (H, W))

def watershed(f, markers, mask=None):
    """Flood f from labeled markers, lowest values first, with a priority queue."""
    lab = markers.copy()
    mask = np.ones(f.shape, bool) if mask is None else mask
    heap, count = [], 0                                       # count breaks ties: first in, first out
    for y, x in zip(*np.nonzero(lab)):                        # markers enter first: they act as
        heapq.heappush(heap, (-np.inf, count, y, x)); count += 1  # the only minima of f
    while heap:
        _, _, y, x = heapq.heappop(heap)
        for dy, dx in ((-1, 0), (1, 0), (0, -1), (0, 1)):     # 4-neighbors
            v, u = y + dy, x + dx
            if 0 <= v < f.shape[0] and 0 <= u < f.shape[1] and mask[v, u] and lab[v, u] == 0:
                lab[v, u] = lab[y, x]                         # the first flood to arrive wins
                heapq.heappush(heap, (f[v, u], count, v, u)); count += 1
    return lab

def match(pred, gt):                                          # best IoU of each true object
    return [max((((pred == j) & (gt == i)).sum() / ((pred == j) | (gt == i)).sum()
                 for j in np.unique(pred[pred > 0])), default=0.0) for i in range(1, gt.max() + 1)]

# 1. Smooth, threshold, and label connected components: touching disks merge.
mask = ndi.gaussian_filter(img, 1.0) > 0.5
cc, n_cc = ndi.label(mask)
print(f"connected components: {n_cc} (true objects: {gt.max()}), "
      f"best IoU per object {np.round(match(cc, gt), 3)}")

# 2. Distance transform; markers at its local maxima (7 x 7 window, distance > 4).
dist = ndi.distance_transform_edt(mask)
peaks = (dist == ndi.maximum_filter(dist, size=7)) & (dist > 4)
markers, n_mark = ndi.label(peaks)
ws = watershed(-dist, markers, mask)
print(f"distance-transform watershed: {n_mark} markers -> {len(np.unique(ws)) - 1} objects, "
      f"best IoU per object {np.round(match(ws, gt), 3)}")

# 3. Watershed of the gradient magnitude: every regional minimum floods.
def grad(f):
    return np.hypot(ndi.sobel(f, 0), ndi.sobel(f, 1))

for sigma in (0.0, 1.0, 2.0):
    g = grad(ndi.gaussian_filter(img, sigma) if sigma else img)
    cross = ndi.generate_binary_structure(2, 1)                    # 4-neighborhood, as in flooding
    minima, n_min = ndi.label(g == ndi.minimum_filter(g, footprint=cross))
    n_reg = len(np.unique(watershed(g, minima)))
    print(f"gradient watershed, smoothing sigma={sigma}: {n_min} minima -> {n_reg} regions")

# 4. Marker-controlled watershed of the gradient: one marker per object from
#    step 2, plus one background marker (far from any foreground pixel).
g = grad(ndi.gaussian_filter(img, 1.0))
bg = ndi.distance_transform_edt(~mask) > 3
mk = markers.copy()
mk[bg] = n_mark + 1
seg = watershed(g, mk)
seg[seg == n_mark + 1] = 0                                    # background label -> 0
fg_iou = ((seg > 0) & (gt > 0)).sum() / ((seg > 0) | (gt > 0)).sum()
print(f"marker-controlled gradient watershed: {n_mark + 1} markers -> {len(np.unique(seg))} regions, "
      f"foreground IoU {fg_iou:.3f}, best IoU per object {np.round(match(seg, gt), 3)}")

Output:

connected components: 2 (true objects: 4), best IoU per object [0.381 0.328 0.281 0.985]
distance-transform watershed: 4 markers -> 4 objects, best IoU per object [0.985 0.985 0.969 0.985]
gradient watershed, smoothing sigma=0.0: 1661 minima -> 1661 regions
gradient watershed, smoothing sigma=1.0: 1164 minima -> 1164 regions
gradient watershed, smoothing sigma=2.0: 540 minima -> 540 regions
marker-controlled gradient watershed: 5 markers -> 5 regions, foreground IoU 0.972, best IoU per object [0.325 0.631 0.417 0.97 ]

Thresholding finds the foreground, but the three overlapping disks form one component (IoU 0.28–0.38 each). The distance transform turns each disk into a cone peaked at its center, with saddles at the necks between disks; flooding the negated distance from its four peaks recovers all four objects with IoU 0.969–0.985. Beucher and Meyer (1993) separated overlapping coffee beans this way.

Flooding the gradient from all its minima gives 1,661 regions for four objects: every noise dip is a regional minimum. Smoothing reduces the count (1,164 at $\sigma = 1$, 540 at $\sigma = 2$) but nowhere near the right number. With five markers, the gradient watershed finds the object–background boundary well (foreground IoU 0.972) and the isolated disk with IoU 0.97, but the boundaries among the three overlapping disks land where noise puts them, because two disks of equal brightness have no gradient crest between them. The relief must contain the boundary to be found.

Over noise seeds 0–9, the distance-transform watershed always found four objects (IoU 0.955–0.991); the minima counts ranged over 1,661–1,715, 1,143–1,188 and 490–565; and in step 4 the overlapping disks scored from 0.325 to 0.941, while foreground IoU stayed at 0.972–0.979. Given the step 4 markers imposed as minima (marker pixels set below every gradient value), scikit-image 0.25.2’s watershed returned identical labels for every seed.

In libraries, skimage.segmentation.watershed(image, markers, mask=...) floods with a priority queue ordered by value and then entry time, and labels optional watershed lines 0. OpenCV’s cv2.watershed(image, markers) implements Meyer’s algorithm on an 8-bit 3-channel image and 32-bit markers, setting boundary pixels to $-1$.

Properties and Behavior

• Closed contours without a threshold. Watershed lines are closed, unlike thresholded gradient maps, which give contours either open or too thick (Beucher & Lantuéjoul, 1979).
• Over-segmentation. Each regional minimum produces a basin, and noise and texture create many minima, so the plain transform oversegments severely, even after filtering (Vincent & Soille, 1991; Beucher & Meyer, 1993).
• Nonlocal. Basin membership cannot be decided locally, which complicates parallel implementation (Roerdink & Meijster, 2000).

Limitations

• Good markers need knowledge of object size, shape or brightness. For very varied samples, Vincent and Soille (1991) found robust markers almost impossible to extract and merged basins instead.
• Distance-transform splitting assumes roughly convex objects: elongated objects produce several maxima and split, and deep overlaps produce one maximum and merge.
• Gradient watersheds are sensitive to noise and texture, and know nothing about object categories.

Variants

• Marker-controlled watershed. Markers are imposed as the only minima by homotopy modification (Beucher & Meyer, 1993), or the image is flooded directly from them; markers may come from the distance transform, h-minima, top-hats or user clicks.
• Hierarchical watersheds and waterfalls. The basins of an oversegmented watershed are merged by a watershed on a graph of adjacent regions valued by their boundaries; repeating this gives a hierarchy (Beucher & Meyer, 1993), and Beucher (1994) describes the related waterfall algorithm.
• Watersheds on edge-weighted graphs. Cousty et al. (2009) defined watershed cuts on edge-weighted graphs, proved their optimality in terms of minimum spanning forests, and gave a linear-time algorithm. Couprie et al. (2011) placed the spanning-forest watershed in a common energy framework with graph cuts, random walker and shortest paths, and introduced the power watershed.
• Compact watershed. Neubert and Protzel (2014) added a weighted Euclidean distance to the seed to the flooding priority, like SLIC’s compactness term, so grid seeds grow into regular superpixels in about 10 ms per $481 \times 321$ image; scikit-image exposes it as compactness.
• Other domains. Immersion extends to other grids, $n$-dimensional images and graphs (Vincent & Soille, 1991).