Segmentation

Region Growing

Segmenting an image by growing regions from seed pixels for as long as a homogeneity criterion holds, as in seeded region growing, or by splitting non-uniform quadtree blocks and merging similar neighbours.

beginner

Region growing segments an image from the inside out. It starts from one or more seed pixels and repeatedly adds neighbouring pixels that are similar enough to the region, judged by a homogeneity criterion such as the difference between a pixel’s gray value and the region’s mean, until no neighbour qualifies. Its counterpart, split-and-merge, works from the outside in: it splits the image into quadrants while a block is not uniform, then merges adjacent blocks that are similar. Both are among the oldest region-based methods of image segmentation. Seeded growing with a fixed tolerance is the flood fill that scikit-image’s documentation compares to the paint-bucket tool of raster graphics programs, and medical imaging toolkits such as ITK ship several region growing filters.

The methods date from the beginnings of computer vision. Brice and Fennema (1970) grew regions by merging; quadtree splitting followed in the early 1970s, and Horowitz and Pavlidis (1976) combined splitting and merging (Haralick & Shapiro, 1985; Szeliski, 2022). Haralick and Shapiro (1985) sorted the region growing methods of their time into single linkage, hybrid linkage and centroid linkage schemes, a classification this article follows. Adams and Bischof (1994) gave a widely cited form, seeded region growing (SRG), in which several seeded regions compete for pixels in order of similarity and no homogeneity threshold is needed. The Adams and Bischof paper could not be read for this article (the publisher’s site does not allow automated access, and no authors’ copy was found), so its description here follows a 1996 dissertation by Bub that restates the algorithm in its own notation. The Horowitz and Pavlidis and the Brice and Fennema papers are likewise described through Haralick and Shapiro’s survey and through Szeliski’s textbook.

Problem

Partition an image, or extract one object from it, into connected regions whose pixels share some property, usually gray level or colour, so that each region is uniform and adjacent regions differ. Haralick and Shapiro (1985) state these two demands and note that they conflict: strictly uniform regions are full of small holes and have ragged boundaries, while insisting on large differences between neighbours merges regions and loses boundaries. Region growing adds a spatial constraint that thresholding lacks: a pixel joins a region only through a chain of neighbours that are also in it.

Inputs and Outputs

  • Input: a grayscale or multichannel image, a neighbourhood (4- or 8-connected in 2D, 6- or 26-connected in 3D), a homogeneity criterion, and, for seeded methods, seed pixels or seed areas with their labels. Single-region growing and split-and-merge also take a threshold; SRG takes none.
  • Output: a binary mask of the grown region, or a label image in which each pixel carries the label of its region; some variants also mark pixels that touch two regions as boundary pixels.

Intuition

Pour paint on a seed and let it spread to every neighbouring pixel that looks like the paint’s colour. If “looks like” means “within TT of the seed’s value”, this is the paint-bucket tool, and the result depends only on the seed and TT. If it means “within TT of the average of what is already painted”, the region’s notion of its own colour adapts as it grows: it can follow a slow change of brightness across an object, but it can also drift across a weak boundary into the background, one similar pixel at a time.

SRG removes the threshold by making regions compete. Every seeded region proposes its most similar unlabelled neighbour, the most similar of all proposals is accepted, and the process repeats until every pixel has been assigned. A pixel on a weak boundary is not asked “is it similar enough?” but “which region is it more similar to?”, and it is decided only after the clearer pixels on both sides.

Split-and-merge needs no seeds. It asks whether a block is uniform; if not, it cuts it into four and asks again. Because cuts follow the quadtree and not the image, adjacent blocks of the same object end up in different branches, and a merge step joins neighbours whose statistics agree.

Algorithm

Single-region growing from a seed. With a threshold TT and a seed p0p_0:

  1. Set R={p0}R = \{p_0\} and put the neighbours of p0p_0 in a queue.
  2. Take a pixel qq from the queue. If it satisfies the criterion (for example ∣g(q)−gˉR∣≤T|g(q) - \bar g_R| \le T, with gˉR\bar g_R the current mean of RR), add it to RR, update the region statistics, and queue its neighbours that are not yet in RR or in the queue.
  3. Repeat until the queue is empty.

With a criterion that compares every pixel to the seed’s value, this is a flood fill. scikit-image’s flood and flood_fill implement it with a breadth-first queue and fill connected pixels within plus or minus tolerance of the seed’s value, inclusive; with no tolerance, only pixels equal to the seed’s value (scikit-image 0.25 documentation and source).

Seeded region growing (Adams and Bischof, as described by Bub, 1996). The seeds are grouped into nn sets A1,…,AnA_1, \dots, A_n (Bub writes S1,…,SnS_1, \dots, S_n, and AiA_i in his definition of δ\delta); a set may be one pixel or a group of pixels. Bub used 8-connected neighbours; the example below uses 4-connected ones.

  1. Label each seed set and compute its mean.
  2. Put every unlabelled pixel that touches a labelled region into a priority queue, keyed by δ\delta, its difference from the region it adjoins (defined below).
  3. Remove the pixel with the smallest δ\delta. If its labelled neighbours all carry one label, give it that label and update the region’s mean. If they carry several labels, either give it the label that minimises δ\delta or flag it as a boundary pixel.
  4. Queue its unlabelled neighbours that are not already queued, keyed by their δ\delta at that moment.
  5. Repeat until the queue is empty: every pixel has then been allocated.

Bub calls the queue a sequentially sorted data structure and notes that what makes the algorithm distinct is that pixels are added in its order.

Split-and-merge.

  1. Start with the whole image (or with square blocks of a chosen size) as the initial segments.
  2. Split every segment that fails the homogeneity test into four quadrants; repeat on the quadrants, down to a minimum block size. The leaves form a quadtree.
  3. Merge adjacent segments, including leaves from different branches, whose statistics are similar enough; repeat until no pair qualifies.

Haralick and Shapiro credit the splitting algorithms to Robertson and to Klinger, both in 1973, and the split-and-merge strategy, which joins adjacent quarters from different split segments, to Horowitz and Pavlidis (1976). Szeliski (2022) describes the Horowitz–Pavlidis method as starting from a medium-grain quadtree segmentation and allowing both merging and splitting operations.

Mathematical Formulation

Homogeneity criteria

Write g(p)g(p) for the gray value of pixel pp, RR for a region with NN pixels, mean Xˉ\bar X and scatter S2=∑p∈R(g(p)−Xˉ)2S^2 = \sum_{p \in R} (g(p) - \bar X)^2. Haralick and Shapiro (1985) distinguish criteria by what a candidate pixel is compared with.

  • Single linkage: two neighbouring pixels are joined when they are similar enough, in the simplest case when ∣g(p)−g(q)∣≤T|g(p) - g(q)| \le T, and segments are the connected components of the resulting graph. Bryant normalized the difference by 2\sqrt{2} times the root mean square of neighbouring pixel differences over the image; for vector-valued pixels a norm of the difference vector is used. Because one leaking link merges two regions, Haralick and Shapiro call these schemes the simplest and the most prone to unwanted merges, a problem they call chaining.
  • Hybrid linkage: each pixel gets a property vector computed from its K×KK \times K neighbourhood, and pixels are linked when their neighbourhoods are similar. One scheme links neighbouring pixels that are both non-edge pixels of an edge operator, so regions are the connected components of the non-edge pixels.
  • Centroid linkage: a pixel is compared with the mean of an existing, not necessarily complete, neighbouring segment, and the mean is updated when it joins. A statistical form tests whether y=g(q)y = g(q) could come from the population of RR. Under the assumption that the pixels of RR and yy are independent and identically distributed normal variables,
T=[(N−1)NN+1 (y−Xˉ)2S2]1/2T = \left[ \frac{(N-1)N}{N+1} \, \frac{(y - \bar X)^2}{S^2} \right]^{1/2}

has a tt distribution with N−1N - 1 degrees of freedom, and yy joins RR when TT is below the critical value tN−1(α)t_{N-1}(\alpha) for a user-chosen significance level α\alpha. Since the critical value falls as the degrees of freedom grow, a pixel must be closer to the mean of a large region than of a small one to join it, for equal scatters. A new one-pixel region has zero scatter, so Haralick and Shapiro add a prior variance VV and an initial size N0N_0, starting each region with scatter N0VN_0 V and size N0N_0.

The flood-fill criterion ∣g(q)−g(p0)∣≤T|g(q) - g(p_0)| \le T compares every candidate with the seed alone. It follows from this definition that the result is the connected component containing p0p_0 of the band {p:g(p0)−T≤g(p)≤g(p0)+T}\{p : g(p_0) - T \le g(p) \le g(p_0) + T\}: a double threshold followed by connected components, independent of the order in which pixels are visited.

The SRG ordering

Let TT now denote the set of unlabelled pixels that border at least one region,

T={x∉⋃iAi  :  N(x)∩⋃iAi≠∅},T = \Big\{ x \notin \textstyle\bigcup_i A_i \; : \; N(x) \cap \bigcup_i A_i \neq \emptyset \Big\},

with N(x)N(x) the neighbours of xx. For x∈Tx \in T adjoining region Ai(x)A_{i(x)}, Bub gives as the simplest measure

δ(x)=∣g(x)−mean⁡y∈Ai(x)g(y)∣,\delta(x) = \Big| g(x) - \operatorname*{mean}_{y \in A_{i(x)}} g(y) \Big| ,

and each step allocates the pixel z∈Tz \in T with δ(z)=min⁡x∈Tδ(x)\delta(z) = \min_{x \in T} \delta(x). When N(x)N(x) meets several regions, i(x)i(x) is the one that minimises δ(x)\delta(x), or xx is put in a set of boundary pixels. Bub stresses that δ\delta can be replaced by any distance suited to the image’s features. The criterion contains no threshold: the seeds alone decide how many regions there are and which features they represent.

Split-and-merge tests

The simplest test calls a block uniform when the difference between its largest and smallest gray values is small enough; Fukada split until the sample variance was small enough (Haralick & Shapiro, 1985). Chen and Pavlidis, as Haralick and Shapiro report, required instead that the mean of each quarter not differ significantly from the block’s mean, ∣Xˉi⋅−Xˉ⋅⋅∣≤ϵ|\bar X_{i\cdot} - \bar X_{\cdot\cdot}| \le \epsilon for i=1,…,4i = 1, \dots, 4, assuming equal and known variances. Without assuming the variance known, Haralick and Shapiro give an FF test of the hypothesis that the quarters have identical means and variances: with KK pixels per quarter, quarter means Xˉi⋅\bar X_{i\cdot} and grand mean Xˉ⋅⋅\bar X_{\cdot\cdot},

F=K∑i=14(Xˉi⋅−Xˉ⋅⋅)2/3∑i=14∑k=1K(Xik−Xˉi⋅)2/(4(K−1))∼F3, 4(K−1),F = \frac{K \sum_{i=1}^{4} (\bar X_{i\cdot} - \bar X_{\cdot\cdot})^2 / 3}{\sum_{i=1}^{4} \sum_{k=1}^{K} (X_{ik} - \bar X_{i\cdot})^2 / \big(4(K-1)\big)} \sim F_{3,\,4(K-1)} ,

and the block is declared not uniform when FF is too high. For merging, Muerle and Allen compared the gray-value distributions of adjacent regions with the Kolmogorov–Smirnov test.

Parameters

  • Seeds. Their positions and number decide what is segmented: in SRG every region comes from a seed set, and an object without a seed is absorbed by its neighbours. Single pixels are fragile in noise. Bub reports that Adams suggested small seed areas instead of single pixels, and found groups of seed points necessary for noisy CT images.
  • Homogeneity threshold TT (single-region growing). Too small and the region stops short, leaving holes and a ragged boundary; too large and it leaks through the weakest point of its boundary. For a statistical criterion, the significance level α\alpha replaces TT.
  • Criterion. Seed value, region mean, a statistical test, or a texture or colour distance; comparisons with the region mean adapt to slow changes, comparisons with the seed do not.
  • Connectivity. 8-connected growth crosses diagonal gaps that 4-connected growth does not, and so leaks more easily.
  • Split-and-merge thresholds and minimum block size. The split test sets how finely blocks are cut, the merge threshold how far merging goes, and the minimum block size how blocky boundaries can be.

Complexity

A flood fill visits each pixel of the region once and examines its neighbours, so, by construction of its queue, it is linear in the region’s size. SRG keeps a priority queue of boundary pixels; with a binary heap, NN insertions and removals cost O(Nlog⁡N)O(N \log N) for an NN-pixel image (our reasoning from the heap’s costs). Bub found a sorted linked list too slow, its insertions costing O(n2)O(n^2) in total, and used a Fibonacci heap. For split-and-merge, Haralick and Shapiro warned that the data structures for images larger than 512×512512 \times 512 are very large, so that on virtual-memory computers of the time paging could dominate; Browning and Tanimoto split and merged within separate tiles and then merged across tile borders. Single-pass centroid linkage, by contrast, needs only the statistics of regions active in the current row, which Haralick and Shapiro bound by the number of pixels in a row.

Implementation

The example compares four methods on a 128×128128 \times 128 synthetic image with known ground truth: a background whose brightness ramps from 0.25 on the left to 0.45 on the right, a bright disk A (mean 0.75) with smooth internal texture, and a dim disk B (0.58) on the bright side of the ramp, only 0.15 to 0.23 above its local background, all with Gaussian noise of standard deviation 0.05. It runs (1) scikit-image’s flood and (2) single-seed centroid linkage at thresholds from 0.04 to 0.16, from five seeds in disk B: its centre and points 12 pixels up, down, left and right. It compares centroid linkage in first-in-first-out (FIFO) and best-first order, runs (3) SRG from 3×33 \times 3 seed areas in each disk and one background seed set in the four corners, with both rules for pixels that touch two regions, and (4) split-and-merge with two split tests: the range test and a standard-deviation test. Masks are scored by intersection over union (IoU) with the true disk and by boundary error, the average symmetric distance in pixels between predicted and true boundary pixels, which also counts the boundaries of holes. Split-and-merge has no seeds, so each disk is scored against the region that matches it best. The coins image bundled with scikit-image has no ground truth, so on it the script reports only quantities that need none: flood areas against tolerance, and region counts. It uses NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2, and runs in under ten seconds.

import heapq
from collections import deque
import numpy as np
from scipy import ndimage as ndi
from skimage import data
from skimage.filters import threshold_otsu
from skimage.morphology import flood

N4 = ((-1, 0), (1, 0), (0, -1), (0, 1))


def make_scene(seed):
    """128 x 128 scene: illumination ramp, bright textured disk A, dim disk B."""
    rng = np.random.default_rng(seed)
    H = W = 128
    yy, xx = np.mgrid[0:H, 0:W]
    gt = np.zeros((H, W), int)
    gt[(yy - 40) ** 2 + (xx - 38) ** 2 <= 22 ** 2] = 1          # disk A, left
    gt[(yy - 84) ** 2 + (xx - 90) ** 2 <= 24 ** 2] = 2          # disk B, on the bright side
    img = 0.25 + 0.20 * xx / (W - 1)                             # background ramp 0.25 -> 0.45
    tex = ndi.gaussian_filter(rng.normal(0, 1, (H, W)), 2)
    img = np.where(gt == 1, 0.75 + 0.05 * tex / tex.std(), img)  # textured bright disk
    img = np.where(gt == 2, 0.58, img)                           # weak-contrast disk
    return img + rng.normal(0, 0.05, (H, W)), gt


def iou(a, b):
    return (a & b).sum() / (a | b).sum()


def bdist(pred, gt):
    """Average symmetric distance (pixels) between the boundaries of two masks."""
    if not pred.any():
        return np.inf
    bp, bg = pred & ~ndi.binary_erosion(pred), gt & ~ndi.binary_erosion(gt)
    return 0.5 * (ndi.distance_transform_edt(~bg)[bp].mean()
                  + ndi.distance_transform_edt(~bp)[bg].mean())


def grow_mean(img, seed, T, order="fifo", r=0):
    """Centroid linkage from a seed (a (2r+1)^2 seed area): add a neighbor if |g - mean(R)| <= T."""
    H, W = img.shape
    inR = np.zeros((H, W), bool)
    queued = np.zeros((H, W), bool)
    y0, x0 = seed
    inR[y0 - r:y0 + r + 1, x0 - r:x0 + r + 1] = True
    queued[inR] = True
    s, n = float(img[inR].sum()), int(inR.sum())
    q, count = deque() if order == "fifo" else [], 0

    def push(y, x):
        nonlocal count
        for dy, dx in N4:
            v, u = y + dy, x + dx
            if 0 <= v < H and 0 <= u < W and not queued[v, u]:
                queued[v, u] = True
                if order == "fifo":
                    q.append((v, u))
                else:
                    heapq.heappush(q, (abs(img[v, u] - s / n), count, v, u)); count += 1

    for y, x in zip(*np.nonzero(inR)):
        push(y, x)
    while q:
        v, u = q.popleft() if order == "fifo" else heapq.heappop(q)[2:]
        if abs(img[v, u] - s / n) <= T:                      # test against the current mean
            inR[v, u] = True; s += img[v, u]; n += 1
            push(v, u)
        else:
            queued[v, u] = False                             # may be re-tested later
    return inR


def srg(img, seeds, mode="min", nbrs=N4):
    """Seeded region growing (Adams & Bischof, as described by Bub, 1996).
    seeds: int array, 0 = unlabeled, 1..K = seed sets. Pixels are taken from a
    priority queue ordered by delta = |g - mean of adjoining region| at insertion.
    mode 'min': a pixel touching several regions joins the one with the smallest
    delta; mode 'boundary': it gets the boundary label -1."""
    H, W = img.shape
    lab = seeds.copy()
    K = lab.max()
    S = np.array([0.0] + [img[lab == k].sum() for k in range(1, K + 1)])
    N = np.array([1] + [(lab == k).sum() for k in range(1, K + 1)])
    inq = np.zeros((H, W), bool)
    heap, count = [], 0

    def adj(v, u):
        return {lab[v + dy, u + dx] for dy, dx in nbrs
                if 0 <= v + dy < H and 0 <= u + dx < W and lab[v + dy, u + dx] > 0}

    def push(y, x):
        nonlocal count
        for dy, dx in nbrs:
            v, u = y + dy, x + dx
            if 0 <= v < H and 0 <= u < W and lab[v, u] == 0 and not inq[v, u]:
                d = min(abs(img[v, u] - S[k] / N[k]) for k in adj(v, u))
                heapq.heappush(heap, (d, count, v, u)); count += 1
                inq[v, u] = True

    for y, x in zip(*np.nonzero(lab)):
        push(y, x)
    while heap:
        _, _, y, x = heapq.heappop(heap)
        ks = adj(y, x)
        if len(ks) > 1 and mode == "boundary":
            lab[y, x] = -1
            continue
        k = min(ks, key=lambda k: abs(img[y, x] - S[k] / N[k]))
        lab[y, x] = k; S[k] += img[y, x]; N[k] += 1
        push(y, x)
    return lab


def split_merge(img, Ts, Tm, pred="range", min_size=2):
    """Quadtree split while a block is not homogeneous, then greedy merging of
    adjacent regions in order of mean difference, while the difference < Tm."""
    lab = np.zeros(img.shape, int)
    leaves = []

    def homog(b):
        return (b.max() - b.min() <= Ts) if pred == "range" else (b.std() <= Ts)

    def rec(y0, x0, h, w):
        b = img[y0:y0 + h, x0:x0 + w]
        if homog(b) or h <= min_size or w <= min_size:
            leaves.append((y0, x0, h, w)); return
        h2, w2 = h // 2, w // 2
        rec(y0, x0, h2, w2); rec(y0, x0 + w2, h2, w - w2)
        rec(y0 + h2, x0, h - h2, w2); rec(y0 + h2, x0 + w2, h - h2, w - w2)

    rec(0, 0, *img.shape)
    for i, (y0, x0, h, w) in enumerate(leaves):
        lab[y0:y0 + h, x0:x0 + w] = i
    n = len(leaves)
    S = ndi.sum(img, lab, np.arange(n)); C = np.bincount(lab.ravel(), minlength=n).astype(float)
    pairs = set()
    for a, b in ((lab[:-1], lab[1:]), (lab[:, :-1], lab[:, 1:])):
        m = a != b
        pairs |= set(zip(np.minimum(a[m], b[m]).tolist(), np.maximum(a[m], b[m]).tolist()))
    nb = [set() for _ in range(n)]
    for a, b in pairs:
        nb[a].add(b); nb[b].add(a)
    parent = list(range(n)); ver = [0] * n
    heap = [(abs(S[a] / C[a] - S[b] / C[b]), a, b, 0, 0) for a, b in pairs]
    heapq.heapify(heap)
    while heap:
        d, a, b, va, vb = heapq.heappop(heap)
        if d >= Tm:
            break
        if parent[a] != a or parent[b] != b or ver[a] != va or ver[b] != vb:
            continue                                          # stale entry
        parent[b] = a; S[a] += S[b]; C[a] += C[b]; ver[a] += 1
        nb[a] |= nb[b]; nb[a] -= {a, b}
        for c in nb[b]:
            nb[c].discard(b)
            if c != a:
                nb[c].add(a)
        for c in nb[a]:
            heapq.heappush(heap, (abs(S[a] / C[a] - S[c] / C[c]), a, c, ver[a], ver[c]))
    root = np.array(parent)
    for i in range(n):                                        # path compression
        r = i
        while root[r] != r:
            r = root[r]
        root[i] = r
    return root[lab], len(leaves)


def best_match(pred, gt_mask):
    """Region of a seedless segmentation that best matches a true object."""
    ids = np.unique(pred[gt_mask])
    j = max(ids, key=lambda j: iou(pred == j, gt_mask))
    return pred == j


def seeds_for(shape, b_seed):
    s = np.zeros(shape, int)
    s[39:42, 37:40] = 1                                       # 3 x 3 seed area in disk A
    s[b_seed[0] - 1:b_seed[0] + 2, b_seed[1] - 1:b_seed[1] + 2] = 2   # disk B
    for y, x in ((1, 1), (1, 124), (124, 1), (124, 124)):     # one background set: 4 corners
        s[y:y + 3, x:x + 3] = 3
    return s


img, gt = make_scene(0)
B = gt == 2
b_seeds = [(84, 90), (72, 90), (96, 90), (84, 78), (84, 102)]   # center, 12 px up/down/left/right
print("disk B seed values:", np.round([img[s] for s in b_seeds], 3))

# 1. Homogeneity threshold and seed: flood (fixed seed value) vs centroid linkage (running mean)
for T in (0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16):
    fl = [flood(img, s, connectivity=1, tolerance=T) for s in b_seeds]
    cl = [grow_mean(img, s, T) for s in b_seeds]
    print(f"T={T:.2f} flood IoU {np.round([iou(m, B) for m in fl], 2)} "
          f"bd {np.round([bdist(m, B) for m in fl], 1)} | centroid IoU {np.round([iou(m, B) for m in cl], 2)} "
          f"bd {np.round([bdist(m, B) for m in cl], 1)}")

# 2. Order dependence: centroid linkage FIFO vs best-first; flood has no order effect
for T in (0.08, 0.12):
    a, b = grow_mean(img, b_seeds[0], T, "fifo"), grow_mean(img, b_seeds[0], T, "best")
    print(f"T={T:.2f} centroid linkage, FIFO vs best-first: {(a != b).sum()} pixels differ "
          f"(sizes {a.sum()} / {b.sum()})")

# 3. Seeded region growing with competing seeds, for each disk-B seed, both boundary rules
for mode in ("min", "boundary"):
    res = []
    for s in b_seeds:
        lab = srg(img, seeds_for(img.shape, s), mode)
        res.append((iou(lab == 2, B), bdist(lab == 2, B), iou(lab == 1, gt == 1), (lab == -1).sum()))
    r = np.array(res)
    print(f"SRG ({mode}): disk B IoU {np.round(r[:, 0], 3)} bd {np.round(r[:, 1], 2)}; "
          f"disk A IoU {np.round(r[:, 2], 3)}; boundary pixels {r[:, 3].astype(int)}")
lab_a = srg(img, seeds_for(img.shape, b_seeds[0]))
lab_b = srg(img, seeds_for(img.shape, b_seeds[0]), nbrs=N4[::-1])
print(f"SRG, neighbor order reversed: {(lab_a != lab_b).sum()} pixels change label")

# 4. Split-and-merge, two homogeneity predicates (range: H&S; std: variance test)
for pred, Ts in (("range", 0.20), ("std", 0.06)):
    sm, nleaf = split_merge(img, Ts, 0.07, pred)
    mA, mB = best_match(sm, gt == 1), best_match(sm, B)
    print(f"split-and-merge ({pred} <= {Ts}, merge < 0.07): {nleaf} leaves -> {len(np.unique(sm))} regions; "
          f"disk A IoU {iou(mA, gt == 1):.3f} bd {bdist(mA, gt == 1):.2f}; "
          f"disk B IoU {iou(mB, B):.3f} bd {bdist(mB, B):.2f}")

# 5. Spread over noise seeds 0-9
rows = []
for k in range(10):
    im, g = make_scene(k)
    Bk = g == 2
    lab = srg(im, seeds_for(im.shape, b_seeds[0]))
    f = flood(im, b_seeds[0], connectivity=1, tolerance=0.10)
    c = grow_mean(im, b_seeds[0], 0.10)
    c3 = grow_mean(im, b_seeds[0], 0.10, r=1)                  # 3 x 3 seed area
    rows.append((iou(lab == 2, Bk), bdist(lab == 2, Bk), iou(f, Bk), iou(c, Bk), iou(c3, Bk), im[b_seeds[0]]))
r = np.array(rows)
print(f"noise seeds 0-9: SRG disk B IoU {r[:, 0].min():.3f}-{r[:, 0].max():.3f}, "
      f"bd {r[:, 1].min():.2f}-{r[:, 1].max():.2f}; flood T=0.10 IoU {r[:, 2].min():.2f}-{r[:, 2].max():.2f}; "
      f"centroid T=0.10 IoU {r[:, 3].min():.2f}-{r[:, 3].max():.2f}")
print("  seed pixel value     ", np.round(r[:, 5], 3))
print("  flood IoU            ", np.round(r[:, 2], 2))
print("  centroid IoU         ", np.round(r[:, 3], 2))
print("  centroid, 3x3 seed   ", np.round(r[:, 4], 2))

# 6. coins: no ground truth, so only label-free measurements
c = data.coins().astype(float)
for s in ((125, 45), (118, 40), (130, 52)):
    areas = [int(flood(c, s, connectivity=1, tolerance=t).sum()) for t in range(10, 101, 10)]
    print(f"coins flood from {s} (value {c[s]:.0f}), tolerance 10..100: areas {areas}")
sm_ = ndi.gaussian_filter(c, 2)
coin_seeds, n_coin = ndi.label(sm_ >= 170)
seeds = coin_seeds.copy()
seeds[sm_ <= 50] = n_coin + 1                                  # one background seed set
lab = srg(c, seeds)
lab_rev = srg(c, seeds, nbrs=N4[::-1])
coin = (lab >= 1) & (lab <= n_coin)
cc, n_cc = ndi.label(coin)
sizes = np.bincount(cc.ravel())[1:]
otsu = c > threshold_otsu(c)
print(f"coins SRG: {n_coin} coin seed components; coin label covers {coin.mean():.3f} of the image "
      f"in {n_cc} components ({(sizes >= 200).sum()} with >= 200 px, sizes {sizes.min()}-{sizes.max()}); "
      f"agreement with Otsu mask {(coin == otsu).mean():.3f}")
coin_rev = (lab_rev >= 1) & (lab_rev <= n_coin)
print(f"coins SRG, neighbor order reversed: {(lab != lab_rev).sum()} pixels change label, "
      f"{(coin != coin_rev).sum()} change between coin and background")
for pred, Ts in (("range", 60), ("std", 15)):
    sm, nleaf = split_merge(c, Ts, 20, pred)
    sizes = np.bincount(np.unique(sm, return_inverse=True)[1].ravel())
    print(f"coins split-and-merge ({pred} <= {Ts}, merge < 20): {nleaf} leaves -> {len(sizes)} regions, "
          f"{(sizes >= 200).sum()} with >= 200 px")

Output:

disk B seed values: [0.607 0.585 0.568 0.573 0.573]
T=0.04 flood IoU [0.08 0.08 0.05 0.14 0.08] bd [13.5 12.7 13.2 10.3 13.2] | centroid IoU [0.   0.06 0.17 0.12 0.21] bd [22.6 15.  10.3 11.2  9.4]
T=0.06 flood IoU [0.7  0.76 0.75 0.76 0.76] bd [3.9 3.9 3.7 3.7 3.7] | centroid IoU [0.76 0.77 0.76 0.77 0.76] bd [3.9 3.8 3.9 3.8 3.9]
T=0.08 flood IoU [0.84 0.9  0.88 0.89 0.89] bd [3.6 3.3 3.4 3.4 3.4] | centroid IoU [0.9 0.9 0.9 0.9 0.9] bd [3.3 3.3 3.4 3.3 3.3]
T=0.10 flood IoU [0.92 0.95 0.94 0.95 0.95] bd [3.1 2.8 2.7 2.6 2.6] | centroid IoU [0.95 0.95 0.95 0.95 0.95] bd [2.6 2.7 2.7 2.7 2.7]
T=0.12 flood IoU [0.96 0.97 0.96 0.97 0.97] bd [2.3 1.8 2.  1.7 1.7] | centroid IoU [0.97 0.97 0.97 0.97 0.97] bd [1.6 1.6 1.6 1.6 1.6]
T=0.14 flood IoU [0.98 0.97 0.81 0.95 0.95] bd [1.4 0.6 3.9 1.3 1.2] | centroid IoU [0.97 0.97 0.97 0.97 0.97] bd [0.8 0.9 0.8 0.9 1. ]
T=0.16 flood IoU [0.98 0.8  0.48 0.51 0.51] bd [ 0.5  4.2 11.5 11.3 11.3] | centroid IoU [0.15 0.15 0.15 0.15 0.15] bd [17.6 17.8 18.  17.6 18.1]
T=0.08 centroid linkage, FIFO vs best-first: 13 pixels differ (sizes 1622 / 1621)
T=0.12 centroid linkage, FIFO vs best-first: 2 pixels differ (sizes 1785 / 1787)
SRG (min): disk B IoU [0.992 0.991 0.991 0.992 0.991] bd [0.1 0.1 0.1 0.1 0.1]; disk A IoU [1. 1. 1. 1. 1.]; boundary pixels [0 0 0 0 0]
SRG (boundary): disk B IoU [0.989 0.99  0.988 0.99  0.989] bd [0.14 0.13 0.15 0.13 0.14]; disk A IoU [0.978 0.978 0.978 0.978 0.978]; boundary pixels [264 264 264 264 264]
SRG, neighbor order reversed: 0 pixels change label
split-and-merge (range <= 0.2, merge < 0.07): 1909 leaves -> 85 regions; disk A IoU 0.670 bd 3.65; disk B IoU 0.956 bd 0.52
split-and-merge (std <= 0.06, merge < 0.07): 478 leaves -> 53 regions; disk A IoU 0.856 bd 2.38; disk B IoU 0.952 bd 0.57
noise seeds 0-9: SRG disk B IoU 0.958-0.997, bd 0.04-0.59; flood T=0.10 IoU 0.00-0.94; centroid T=0.10 IoU 0.00-0.96
  seed pixel value      [0.607 0.647 0.596 0.665 0.494 0.566 0.618 0.557 0.667 0.621]
  flood IoU             [0.92 0.75 0.93 0.55 0.   0.94 0.89 0.9  0.07 0.89]
  centroid IoU          [0.95 0.96 0.94 0.96 0.   0.96 0.95 0.95 0.95 0.95]
  centroid, 3x3 seed    [0.95 0.96 0.94 0.96 0.94 0.96 0.95 0.95 0.95 0.95]
coins flood from (125, 45) (value 179), tolerance 10..100: areas [4, 729, 977, 1084, 1211, 1293, 1322, 1335, 1357, 35094]
coins flood from (118, 40) (value 157), tolerance 10..100: areas [6, 90, 412, 768, 989, 1098, 1192, 36017, 47099, 66178]
coins flood from (130, 52) (value 148), tolerance 10..100: areas [1, 3, 7, 32, 797, 1015, 35112, 46129, 63996, 79918]
coins SRG: 54 coin seed components; coin label covers 0.292 of the image in 22 components (22 with >= 200 px, sizes 1081-3110); agreement with Otsu mask 0.887
coins SRG, neighbor order reversed: 0 pixels change label, 0 change between coin and background
coins split-and-merge (range <= 60, merge < 20): 9244 leaves -> 1312 regions, 38 with >= 200 px
coins split-and-merge (std <= 15, merge < 20): 9406 leaves -> 1381 regions, 39 with >= 200 px

The threshold. Both single-seed methods have a narrow working range on disk B. At T≤0.06T \le 0.06 they reach IoU 0.77 at best: noise pushes many disk pixels outside the band, the region stops at them, and holes and ragged edges give boundary errors of 3.7 to 22.6 pixels. Raising TT closes the holes, and centroid linkage reaches IoU 0.97 with a boundary error of about one pixel or less at T=0.14T = 0.14. At 0.16 it leaks into the background from every seed (IoU 0.15). The local contrast of 0.15 is about three noise standard deviations, so a criterion loose enough to admit the disk’s noisy pixels also admits the background’s brightest ones, and one chain of them is enough.

The seed. The flood fill compares everything with the seed’s own noisy value. Its five seeds lie between 0.568 and 0.607, and at T=0.14T = 0.14 and 0.160.16 their results range from IoU 0.98 to 0.48; only the brightest seed, whose band stops above most of the background, avoids the leak at 0.16. Centroid linkage replaces the seed value by a running mean that soon averages over many pixels, and from T=0.06T = 0.06 to 0.140.14 its five results agree to within 0.01 IoU. Over noise seeds 0–9, a single-pixel seed failed outright when its value was an outlier: the seed at 0.494 reached IoU 0.00 under both methods, and the flood from 0.667 reached only 0.07. Starting centroid linkage from the 3×33 \times 3 seed area removed the failure (IoU 0.94–0.96 for all ten), as the advice to use seed areas leads one to expect.

Competition instead of a threshold. SRG segmented disk B with IoU 0.991–0.992 and a boundary error of 0.1 pixels from all five seed positions, and disk A perfectly, with no threshold to choose; over noise seeds 0–9 disk B scored 0.958–0.997. A pixel near the weak boundary is allocated only when it is the most similar candidate left, and then to whichever region it resembles more, so leakage would require a background pixel to look more like the disk than the disk’s own unallocated pixels do. With the boundary-label rule, 264 pixels touching two regions are left unlabelled, which costs disk A about 2% of its IoU in this scoring.

Order. FIFO and best-first centroid linkage differed by 13 pixels at T=0.08T = 0.08 and 2 at T=0.12T = 0.12, out of about 1,600 to 1,800: because the mean changes as pixels join, a pixel tested early can be rejected and a pixel tested late accepted. The flood fill has no such dependence. Reversing the order in which SRG visits neighbours changed no labels, on the synthetic image or on the 8-bit coins image, where equal δ\delta values are common.

Split-and-merge. Without seeds, split-and-merge found disk B with IoU 0.95–0.96, but its boundary error of about 0.5 pixels is five times SRG’s, because boundaries follow block edges. No single region matched the textured disk A well (IoU 0.67 with the range test, 0.86 with the standard-deviation test): the range of a block grows with texture and noise, so the range test split the image into 1,909 leaves, against 478 for the standard deviation. The block-edge boundaries are the squarish boundaries Haralick and Shapiro describe.

Coins. The flood areas show the leak directly. From the brightest seed (179), the region grows slowly from 729 pixels at tolerance 20 to 1,357 at 90, about the size of a coin, then jumps to 35,094 at 100 as it escapes into the background. The two darker seeds in the same coin leak at tolerances 80 and 70: the tolerance that works depends on where the user clicks. SRG, seeded by 54 bright components and a set of dark background pixels, produced 22 coin-sized regions of 1,081 to 3,110 pixels for the 24 coins in the image (our count). Since a region can only grow from a seed, a coin with no pixel above the seed level cannot get one. The SRG mask and Otsu’s global threshold disagree on 11% of pixels; without ground truth this example cannot say which is right, although uneven illumination, as in coins, is a weakness of a single global threshold that the thresholding article discusses. Split-and-merge produced more than 1,300 regions, 38 or 39 of them with at least 200 pixels, far more large regions than objects.

The code shown is the third version of the script. After the first run, the coins flood areas were reformatted from NumPy integer representations to plain integers; the noise-seed loop was extended to print each seed’s value and the per-seed IoUs, because the first version printed only ranges, which hid why single-pixel methods reached IoU 0.00; and the centroid-linkage run from a 3×33 \times 3 seed area was added, with a seed-radius argument whose default reproduces the earlier runs, to separate the effect of seed areas from competition between seeds. The size range of the SRG coin regions was also added to the printout. After the second run, the reversed-order SRG run on coins was added, because the synthetic image, being floating point, has no ties in δ\delta. No other line of the output changed between versions.

Properties and Behavior

  • Connected, closed regions. Every grown region is connected by construction, and SRG assigns every pixel. Pixels enclosed by a region, such as the dark relief inside a coin, can only be reached through it and so join it; a global threshold leaves them as holes.
  • Leakage and chaining. A region grows wherever its criterion is met along some path, so a single gap in a boundary is enough to let it flood a neighbouring area. Haralick and Shapiro (1985) found single linkage the most prone to such unwanted merges, hybrid and centroid linkage better, and split-and-merge less subject to them, at the cost of memory and blocky boundaries.
  • Order dependence. Whenever the criterion depends on the region’s current state, the result can depend on the order of processing. Haralick and Shapiro call this dependence on the order in which pixels and regions are examined inherent to region growing schemes, adding that the differences caused by scan order are usually minor. They point out that a single left-to-right, top-to-bottom pass cannot make the two arms of a V-shaped region one segment, so the pass must be followed by a merging pass, for which they used the two-pass label propagation of the connected-components algorithm of Lumia et al. A criterion fixed in advance, as in the flood fill, removes the dependence. In SRG, by our reading of Bub’s description, each pixel’s label is decided when it leaves the queue, from the regions it touches at that moment, but its place in the queue was set by δ\delta when it entered, which can be out of date by then; ties between equal δ\delta values are broken by the queue. In the example neither effect changed a label.
  • Seed dependence. Seeded methods inherit the seeds’ errors: a seed on noise or on a boundary starts from a misleading value, and in SRG an object without a seed is absorbed. Bub writes that the choice of seeds embeds the decision of what is a feature of interest.

Limitations

  • Seeds must be supplied. By a user, by a threshold, or by a detector; Bub’s dissertation devotes much of a chapter to choosing seeds automatically for brain CT from histogram peaks and thresholded blocks.
  • Gray-level homogeneity. The criteria assume that a region is roughly uniform in the measured property. Textured objects and smooth shading break that assumption unless the criterion uses texture features or region models; Pavlidis, as Haralick and Shapiro report, merged regions approximated by fitted functions instead of constants.
  • Weak boundaries. Where an object’s contrast with its surroundings is within the noise, single-region growing has no threshold that both fills the object and stops at its edge, as disk B shows.
  • Blocky boundaries in split-and-merge. Quadtree splits produce squarish, slightly artificial boundaries, as Haralick and Shapiro put it, which merging does not fully remove.
  • No semantics. Haralick and Shapiro end their survey by noting that whether two segments belong together can depend on the object parts they represent, not on gray-level distributions, and recommend semantic information as the area for further study.

Variants

  • Merging from small regions. Brice and Fennema (1970) started from segments of pixels with identical intensity and merged adjacent regions if a significant fraction of their common border had a small intensity difference across it (Haralick & Shapiro, 1985); Szeliski (2022) describes their dual grid of boundaries between pixels and merging based on relative boundary lengths and edge strengths. The Felzenszwalb–Huttenlocher algorithm is a later graph-based merging method, and merged small regions are one way to produce superpixels.
  • Statistical centroid linkage. Variants differ in their test: the absolute difference from the mean of a neighbouring region (Nagy and Tolaba), a tt test on the difference between the pixel and the nearest region (Gupta, Kettig, Landgrebe and Wintz), an FF test for multiband data (Kettig and Landgrebe), or a statistic that demands larger differences when a region’s scatter is high (Levine and Shaheen), all as Haralick and Shapiro summarize them.
  • Confidence intervals and isolated seeds. The Insight Toolkit (ITK) 5.4.0 offers several region growing filters. ConnectedThreshold keeps pixels connected to a seed within a lower and an upper threshold; NeighborhoodConnected requires every pixel of a neighbourhood to lie in the range. ConfidenceConnected takes the mean and variance around the seed, grows the region within the mean plus or minus a multiplier times the standard deviation, then re-estimates both from the whole region and repeats for a set number of iterations; its documentation notes that a multiplier of 2.5 would capture 99% of samples if intensities were Gaussian. IsolatedConnected uses a binary search to find the threshold that includes one set of seeds and excludes another.
  • Seeds as markers in other algorithms. Marker-controlled watershed also floods from seeds with a priority queue, but the priority is the fixed value of a relief image, usually the gradient magnitude, while in SRG it is the distance to a region’s current mean. Graph cuts use seeds as hard constraints and find a globally optimal binary segmentation under region and boundary terms; interactive segmentation tools likewise start from user clicks or strokes.
  • Pyramid linking. Burt, Hong and Rosenfeld linked each pixel to a parent in a reduced-resolution pyramid and iterated parent means and child links, a spatial clustering that Haralick and Shapiro classify separately from region growing.

Medical Imaging

Region growing has a long record in medical image analysis, where an operator can place a seed in an organ or lesion. Litjens et al. (2017), surveying deep learning in medical image analysis, describe the field’s systems from the 1970s to the 1990s as chains of low-level pixel processing, such as edge and line detectors and region growing, with mathematical modelling, built into rule-based programs for particular tasks.

Bub (1996) segmented brain CT with SRG after extracting the intracranial volume and suppressing noise, choosing seed groups automatically. Against a qualified operator’s manual outlines, the SRG volume of a cyst over 41 slices was 56.056 cm³ against 56.184 cm³, a difference of 0.228%; for part of the ventricular system over 20 slices it was 19.47 cm³ against 19.71 cm³, which the dissertation reports as a 0.24% difference, although these volumes differ by about 1.2%. On simulated phantoms, the abstract reports a maximum volume error of 15% for objects larger than 741.51 mm³ with a mean gray-level distance of 10 from any other object. These are single image sets compared with one operator.

Region growing also appears as a component of larger systems. In the LUNA16 challenge on pulmonary nodule detection in CT, the M5L system’s Multi Opening and Threshold CAD obtained the lung volume by 3D region growing, with procedures to exclude the trachea and separate the lungs, before detecting nodule candidates (Setio et al., 2017). In the BRATS brain tumour benchmark, annotators outlined structures in every third axial slice and interpolated the segmentation with what Menze et al. (2015) describe as morphological operators (region growing) in the 3D Slicer software, then corrected the result by hand. Learned models are now the main approach to medical image segmentation, as that article describes; in these sources region growing serves as a step for delineating a high-contrast structure or for speeding up manual annotation.

Related

  • Watershed

    A morphological segmentation algorithm that treats an image, usually its gradient magnitude, as a relief and floods it from its minima or from markers, so that regions are catchment basins and boundaries are the lines where floods meet.

  • Graph-Based Image Segmentation

    The Felzenszwalb–Huttenlocher algorithm, which scans the edges of a pixel graph in order of increasing weight and merges two regions unless the edge between them exceeds one region's internal variation plus a size-dependent tolerance set by a scale parameter k.

  • Superpixels

    An oversegmentation of an image into a few hundred compact, roughly uniform regions that follow image boundaries, used in place of pixels to make later processing cheaper.

  • Graph Cuts

    Energy minimization by minimum s–t cut, which labels the pixels of an image optimally for binary energies with submodular pairwise terms and approximately, through expansion and swap moves, for many labels.

  • Medical Image Segmentation

    Delineating organs, tumors, lesions and other structures voxel by voxel in CT, MRI and other medical images; the 2D, 2.5D and 3D approaches, losses for very small targets, voxel spacing and anisotropy, metrics such as Dice and surface distances with their pitfalls, and benchmarks such as the Medical Segmentation Decathlon.

  • Computed Tomography (CT)

    An X-ray imaging method that measures attenuation along many lines through the body and reconstructs cross-sectional images from them, with each voxel calibrated in Hounsfield units.

References

  1. Brice, C. R. & Fennema, C. L. (1970). Scene Analysis Using Regions. Artificial Intelligence, 1(3–4), 205–226.
  2. Horowitz, S. L. & Pavlidis, T. (1976). Picture Segmentation by a Tree Traversal Algorithm. Journal of the ACM, 23(2), 368–388.
  3. Haralick, R. M. & Shapiro, L. G. (1985). Image Segmentation Techniques. Computer Vision, Graphics, and Image Processing, 29(1), 100–132 (read in the scan on Haralick's website).
  4. Adams, R. & Bischof, L. (1994). Seeded Region Growing. IEEE Transactions on Pattern Analysis and Machine Intelligence, 16(6), 641–647.
  5. Bub, A. M. (1996). Segmentation of Brain X-Ray CT Images Using Seeded Region Growing. MSc (Eng) dissertation, University of Cape Town.
  6. Menze, B. H., Jakab, A., Bauer, S., Kalpathy-Cramer, J., Farahani, K., Kirby, J., Burren, Y., Porz, N., Slotboom, J., Wiest, R., et al. (2015). The Multimodal Brain Tumor Image Segmentation Benchmark (BRATS). IEEE Transactions on Medical Imaging, 34(10), 1993–2024 (read in the accepted version, HAL hal-00935640v2).
  7. Litjens, G., Kooi, T., Ehteshami Bejnordi, B., Setio, A. A. A., Ciompi, F., Ghafoorian, M., van der Laak, J. A. W. M., van Ginneken, B. & Sánchez, C. I. (2017). A survey on deep learning in medical image analysis. Medical Image Analysis, 42, 60–88 (read in arXiv:1702.05747v2).
  8. Setio, A. A. A., Traverso, A., de Bel, T., Berens, M. S. N., van den Bogaard, C., Cerello, P., Chen, H., Dou, Q., Fantacci, M. E., Geurts, B., et al. (2017). Validation, comparison, and combination of algorithms for automatic detection of pulmonary nodules in computed tomography images: The LUNA16 challenge. Medical Image Analysis, 42, 1–13 (read in the authors' post-print, University of Pisa repository).
  9. Szeliski, R. (2022). Computer Vision: Algorithms and Applications (2nd ed.). Springer. Section 7.5 (read in the author's final draft of 30 September 2021).
  10. scikit-image documentation (version 0.25): skimage.morphology.flood and flood_fill.
  11. Insight Toolkit (ITK) 5.4.0 source: RegionGrowing module (ConnectedThreshold, NeighborhoodConnected, ConfidenceConnected and IsolatedConnected image filters).