Graph-Based Image Segmentation

Graph-based image segmentation, in the sense used here, is the algorithm of Pedro Felzenszwalb and Daniel Huttenlocher, published in the International Journal of Computer Vision in 2004 and referred to by its authors’ names (Felzenszwalb–Huttenlocher, FH below) or by the paper’s title, “Efficient Graph-Based Image Segmentation”. It treats the image as a graph whose nodes are pixels and whose edges join neighbouring pixels, weighted by how different they are. It visits the edges from the smallest weight to the largest, as Kruskal’s minimum spanning tree algorithm does, and merges the two regions an edge connects unless the edge is clearly larger than the variation already inside one of them. One parameter, $k$, sets how much larger counts as clearly larger for small regions. The algorithm runs in $O(m \log m)$ time for $m$ edges, and its output satisfies two global properties that the paper defines and proves, even though every decision is greedy. It is one of the classical methods of image segmentation, and it became a standard source of superpixels and of the starting regions for selective search.

The definitions, proofs and experiments below come from the authors’ manuscript of the journal paper on the first author’s Brown University page; it carries their affiliations at the MIT Artificial Intelligence Lab and Cornell University, and the Springer version could not be read, so page-level details may differ from the published text. The authors had presented an earlier conference paper, “Image Segmentation Using Local Variation”, at CVPR 1998; it was not read for this article.

Problem

Felzenszwalb and Huttenlocher wanted a segmentation method that could be used as widely as edge detection: fast enough to run in nearly linear time, and precise enough that one could say what its output means. They illustrated the second requirement with a synthetic image (their Figure 1) that most people see as three regions: an intensity ramp on the left, a constant region on the right, and a rectangle of high-variability noise inside the constant region. Two observations follow from it.

• Variation alone is not evidence of a boundary. Intensity changes a lot inside the ramp and inside the noisy rectangle, so a method that assumes nearly constant regions splits both.
• Purely local decisions cannot work. The intensity step between the ramp and the constant region is smaller than many differences between neighbouring pixels inside the noisy rectangle. Any fixed threshold on neighbour differences either splits the rectangle or merges the ramp with the constant region.

The paper reviews why earlier methods fail on this image. Zahn’s (1971) method, as the paper describes it, breaks the large-weight edges of the minimum spanning tree, which is a fixed threshold in disguise. Splitting and merging regions by a uniformity criterion either splits the noisy region or merges it with its surroundings. Feature-space clustering assumes piecewise constant images. Normalized cuts captures non-local structure, but the authors note that it characterizes each cut rather than the final segmentation, leads to an NP-hard problem whose approximation error is not well understood, and was too slow for many applications. FH answers with an adaptive criterion: the evidence for a boundary is the difference across it measured relative to the variation inside the regions on either side.

Inputs and Outputs

The input is an undirected graph $G = (V, E)$ with non-negative edge weights $w(v_i, v_j)$ that measure dissimilarity. For images, $V$ is the set of pixels. The paper builds two kinds of graph:

• a grid graph, joining each pixel to its 8-connected neighbours with weight $|I(p_i) - I(p_j)|$, the absolute intensity difference after slight Gaussian smoothing;
• a nearest-neighbour graph, mapping each pixel to a point $(x, y, r, g, b)$ and joining it to its ten nearest neighbours in that 5-dimensional feature space, with Euclidean distance as the weight.

The output is a partition $S$ of $V$ into components, each connected in the graph. With the grid graph every region is connected in the image; with the nearest-neighbour graph a region can consist of pieces that are far apart in the image but similar in colour. There is no class label, no score and no hierarchy, only one flat partition per value of $k$.

Intuition

Think of each region as having a typical variation: the largest step needed to walk across it without leaving it. A smooth sky has a small one; grass or a patterned shirt has a large one. An edge between two regions is evidence of a boundary when it is larger than the typical variation of at least one of them, so a flat region keeps its border with a textured neighbour. A step of 20 grey levels is a boundary between two flat regions but not between two patches of the same grass.

Processing edges from the smallest to the largest makes this cheap. When an edge is reached, every smaller edge has already been examined, so the regions it connects are already formed at the scale of that edge, and their typical variation is simply the largest edge that has merged anything inside them. A single pixel has no variation yet, so the method adds a tolerance that is large for small regions and shrinks as they grow; this lets pixels group before their variation can be estimated, and it is where the scale parameter $k$ enters.

Algorithm

The paper’s Algorithm 1, for a graph with $n$ vertices and $m$ edges:

1. Sort the edges by non-decreasing weight, $o_1, \dots, o_m$.
2. Start with every vertex in its own component.
3. For $q = 1, \dots, m$: let $o_q = (v_i, v_j)$ join components $C_i$ and $C_j$. If $C_i \neq C_j$ and $w(o_q) \le \mathrm{MInt}(C_i, C_j)$, merge the two components; otherwise do nothing.
4. Return the components.

$\mathrm{MInt}$, the minimum internal difference, is defined in the next section. The components are kept in a disjoint-set forest (union-find) with union by rank and path compression. For each component the implementation stores only its size and its internal difference; when two components merge, the size is the sum and the internal difference is the weight of the edge that caused the merge.

Released implementations add a step the paper does not have. The authors’ C++ code, scikit-image and OpenCV all make a second pass over the sorted edges and merge any two adjacent components of which at least one is smaller than a minimum size; the Parameters section covers it.

Mathematical Formulation

The internal difference of a component $C$ is the largest weight in its minimum spanning tree:

$$\mathrm{Int}(C) = \max_{e \in \mathrm{MST}(C, E)} w(e).$$

The paper’s intuition is that $C$ stays connected only once edges of weight $\mathrm{Int}(C)$ are allowed. The difference between two components is the smallest edge joining them, or infinity if none does:

$$\mathrm{Dif}(C_1, C_2) = \min_{v_i \in C_1,\ v_j \in C_2,\ (v_i, v_j) \in E} w(v_i, v_j).$$

The pairwise comparison predicate says there is evidence for a boundary when the difference exceeds the minimum internal difference:

$$D(C_1, C_2) = \begin{cases} \text{true} & \text{if } \mathrm{Dif}(C_1, C_2) > \mathrm{MInt}(C_1, C_2), \\ \text{false} & \text{otherwise,} \end{cases}$$ $$\mathrm{MInt}(C_1, C_2) = \min\bigl(\mathrm{Int}(C_1) + \tau(C_1),\ \mathrm{Int}(C_2) + \tau(C_2)\bigr), \qquad \tau(C) = \frac{k}{|C|}.$$

Because of the minimum, a boundary needs to be large relative to only one of the two regions: a flat region next to a textured one keeps its boundary even though the textured region’s internal difference would absorb the step. $\tau$ is the threshold function: a single pixel has $\mathrm{Int} = 0$, which says nothing about local variation, so $\tau$ demands stronger evidence for a boundary around small components. The paper stresses that $k$ is not a minimum component size; a small component survives whenever the step to its neighbours is large enough. Any non-negative function of a component can replace $k/|C|$ without changing the algorithm’s guarantees; the authors suggest, for example, a $\tau$ that is large for components of an unwanted shape (long and thin, measured by perimeter over area) so that they are merged aggressively.

Too fine and too coarse. A segmentation is too fine if some pair of neighbouring regions has no evidence for a boundary between them (Definition 1). A segmentation $T$ is a refinement of $S$ if each component of $T$ lies inside a component of $S$, and a proper refinement if $T \neq S$; $S$ is too coarse if it has a proper refinement that is not too fine (Definition 2). A segmentation that is neither always exists (Property 1: start from the single component and keep refining), but it need not be unique.

The paper then proves, for the segmentation $S$ returned by Algorithm 1:

• Lemma 1. If an edge joins two distinct components and does not merge them, at least one of the two is already final. Consequently the edge that merges two components is the minimum edge between them, and the merging edges are exactly those Kruskal’s algorithm would choose for each component’s minimum spanning tree. This is why storing the merging edge’s weight keeps $\mathrm{Int}$ exact.
• Theorem 1. $S$ is not too fine.
• Theorem 2. $S$ is not too coarse.
• Theorem 3. $S$ does not depend on which non-decreasing order is used when edges have equal weights.

The appendix shows how delicate this is. Replacing the minimum in $\mathrm{Dif}$ by a quantile, such as the median edge between two regions, which would be more robust to a single low-weight edge, makes finding a segmentation that is neither too coarse nor too fine NP-hard; the proof reduces the minimum ratio cut problem to it.

In the nearest-neighbour graph the quantities have a geometric reading: if every feature point is replaced by a ball of radius $r$, the balls of $C$ form one connected volume exactly when $r \ge \mathrm{Int}(C)/2$, and $\mathrm{Dif}(C_1, C_2)$ is the dilation needed to connect the two sets. The authors relate this to the mean-shift clustering of Comaniciu and Meer, which dilates points by a fixed radius after its own transformation of the data, whereas FH chooses the radius adaptively per component (see mean shift).

Parameters

Scale $k$. $k$ sets a scale of observation: larger values prefer larger components. The tolerance $k/|C|$ is measured in units of edge weight, so a component of $k$ pixels tolerates edges one intensity level above its internal difference (with weights on the 0–255 scale), and a component of $k/10$ pixels tolerates ten levels; this is arithmetic from the definition, not a statement in the paper. The paper used two settings throughout: $k = 150$ for the $128 \times 128$ images of the COIL object database and $k = 300$ for images of $320 \times 240$ pixels or larger. Because $|C|$ counts pixels, the same object at a lower resolution has fewer pixels and therefore a larger tolerance, so a given $k$ yields coarser segmentations of smaller images.

Smoothing $\sigma$. The paper always smooths with a Gaussian of $\sigma = 0.8$ before computing weights, which it describes as producing no visible change while removing digitization artifacts; its nearest-neighbour example on a $40 \times 32$ synthetic image used $\sigma = 0$.

Colour. The paper segments colour images three times, once per red, green and blue plane, and puts two neighbouring pixels in the same region only if they share a component in all three. It notes that one run on a graph weighted by colour distance is the alternative, but that intersecting the three planes gave better results in its experiments. The released implementations do the alternative: the authors’ C++ code, scikit-image and OpenCV all build one graph with the Euclidean distance between smoothed colour vectors.

Minimum size. The paper has no minimum size, but the authors’ C++ code takes a third argument, a minimum component size “enforced by post-processing”, and so do the libraries:

Source	$\sigma$	$k$	Minimum size	Neighbourhood	Merge test
Paper (2004)	0.8	150 or 300	none	8-connected	$w \le \mathrm{MInt}$
Authors’ C++ code, README “typical parameters”	0.5	500	20	8-connected	$w \le \mathrm{MInt}$
scikit-image 0.25.2 felzenszwalb defaults	0.8	1 (scale)	20	8-connected	$w < \mathrm{MInt}$
OpenCV 4.10.0 createGraphSegmentation defaults	0.5	300	100	4-connected	$w \le \mathrm{MInt}$
Selective search Matlab code (bundles the authors’ C++ code)	0.8	50–300 in the paper’s modes	equal to $k$	8-connected	$w \le \mathrm{MInt}$

The authors’ web page shows example results at $\sigma = 0.5$ with $k = 500$ and minimum 50, and with $k = 1000$ and minimum 100. scikit-image calls $k$ scale; it converts the image to floating point in $[0, 1]$ and divides scale by 255, which its source comments say makes it behave like the reference implementation, so scale has the same meaning as the paper’s $k$ for 8-bit images. By that reading its default of 1 is a very small scale, which should give a very fine oversegmentation.

Complexity

Sorting takes $O(m \log m)$ time, or linear time with counting sort when the weights are integers. The merging pass takes $O(m\,\alpha(m))$, where $\alpha$ is the inverse Ackermann function: each edge needs at most three disjoint-set operations, and $\mathrm{MInt}$ is computed in constant time from the stored sizes and internal differences. A grid graph, or a graph with a fixed number of nearest neighbours per pixel, has $m = O(n)$ edges for $n$ pixels, so the whole algorithm runs in $O(n \log n)$; the abstract describes this as nearly linear in the number of graph edges. Memory is linear: the implementation holds the edge list plus a parent, rank, size and internal difference per pixel.

The paper reports that the method takes a small fraction of a second on its $320 \times 240$ Figure 1 image and generally a fraction of a second, without naming the hardware, and that it had been used in large-scale image database applications. Guimarães et al. (2017) timed 0.157 s per segmentation on a $321 \times 481$ image on a 2.50 GHz Xeon. For the nearest-neighbour graph, the authors found the neighbours with the approximate nearest-neighbour library ANN, which they found fast in practice for several hundred thousand 5-dimensional points.

Implementation

The example implements Algorithm 1 in NumPy and plain Python, checks it against skimage.segmentation.felzenszwalb, and then uses it to look at the behaviour the paper describes. The implementation builds the 8-connected graph in the same edge order as scikit-image and sorts with the same np.argsort, so that equal weights are visited in the same order. Two details are left open by the sources and are printed both ways: whether the merge test is strict (scikit-image uses $<$, the paper and the C++ code $\le$), and whether thresholds are held in single precision (scikit-image’s Cython code stores $\mathrm{Int} + \tau$ in float32 variables). Data: the bundled coins (grey, $303 \times 384$) and chelsea (colour, $300 \times 451$) images and a synthetic image with known regions. Versions: NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2; the whole script runs in about 15 seconds.

import numpy as np
from scipy import ndimage as ndi
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import connected_components
from skimage import data, transform
from skimage.segmentation import felzenszwalb
from skimage.util import img_as_float64


def grid_graph(img, sigma):
    """8-connected grid graph; weights = Euclidean distance between smoothed pixels."""
    img = np.atleast_3d(img).astype(np.float64)
    if sigma > 0:
        img = ndi.gaussian_filter(img, [sigma, sigma, 0], mode="reflect", truncate=4.0)
    h, w = img.shape[:2]
    idx = np.arange(h * w).reshape(h, w)
    pairs = [  # (a, b) slices: right, down, down-right, up-right (scikit-image's order)
        ((slice(None), slice(1, None)), (slice(None), slice(None, -1))),
        ((slice(1, None), slice(None)), (slice(None, -1), slice(None))),
        ((slice(1, None), slice(1, None)), (slice(None, -1), slice(None, -1))),
        ((slice(None, -1), slice(1, None)), (slice(1, None), slice(None, -1))),
    ]
    ea, eb, wt = [], [], []
    for sa, sb in pairs:
        d = img[sa] - img[sb]
        wt.append(np.sqrt(np.sum(d * d, axis=-1)).ravel())
        ea.append(idx[sa].ravel())
        eb.append(idx[sb].ravel())
    return h * w, np.concatenate(ea), np.concatenate(eb), np.concatenate(wt)


def segment_graph(n, ea, eb, wt, k, order=None, min_size=1, strict=True, f32=True):
    """Felzenszwalb-Huttenlocher Algorithm 1 (Kruskal-style merging), then min-size pass."""
    if order is None:
        order = np.argsort(wt)                 # default sort kind, as scikit-image
    a, b, w = ea[order].tolist(), eb[order].tolist(), wt[order].tolist()
    parent = list(range(n))
    size = [1] * n
    internal = [0.0] * n                       # Int(C): largest MST edge so far

    def find(x):
        while parent[x] != x:
            parent[x] = parent[parent[x]]
            x = parent[x]
        return x

    for i in range(len(w)):
        ra, rb = find(a[i]), find(b[i])
        if ra == rb:
            continue
        ta = internal[ra] + k / size[ra]       # Int(C) + tau(C), tau = k / |C|
        tb = internal[rb] + k / size[rb]
        if f32:
            ta, tb = float(np.float32(ta)), float(np.float32(tb))
        mint = min(ta, tb)                     # MInt(C1, C2)
        if (w[i] < mint) if strict else (w[i] <= mint):
            if size[ra] < size[rb]:
                ra, rb = rb, ra
            parent[rb] = ra
            size[ra] += size[rb]
            internal[ra] = w[i]                # the merging edge is the new MST maximum
    for i in range(len(w)):                    # post-processing: absorb small components
        ra, rb = find(a[i]), find(b[i])
        if ra != rb and (size[ra] < min_size or size[rb] < min_size):
            if size[ra] < size[rb]:
                ra, rb = rb, ra
            parent[rb] = ra
            size[ra] += size[rb]
    roots = np.array([find(x) for x in range(n)])
    return np.unique(roots, return_inverse=True)[1]


def fh_segment(img, k, sigma=0.8, min_size=1, strict=True, f32=True):
    """Image in [0, 1]; k on the 0-255 scale, as in the paper and scikit-image."""
    img = img_as_float64(img)
    n, ea, eb, wt = grid_graph(img, sigma)
    lab = segment_graph(n, ea, eb, wt, k / 255.0, None, min_size, strict, f32)
    return lab.reshape(img.shape[:2])


def same_partition(p, q):
    p, q = p.ravel(), q.ravel()
    joint = np.unique(p.astype(np.int64) * (q.max() + 1) + q).size
    return joint == np.unique(p).size == np.unique(q).size


def n_regions(lab):
    return np.unique(lab).size


coins, chelsea = data.coins(), data.chelsea()

print("A. NumPy implementation against skimage.segmentation.felzenszwalb (sigma 0.8)")
for name, img in [("coins", coins), ("chelsea", chelsea)]:
    for min_size in (1, 20):
        for k in (100, 300, 1000):
            ref = felzenszwalb(img, scale=k, sigma=0.8, min_size=min_size,
                               channel_axis=None if img.ndim == 2 else -1)
            mine32 = fh_segment(img, k, 0.8, min_size, strict=True, f32=True)
            mine64 = fh_segment(img, k, 0.8, min_size, strict=True, f32=False)
            le = fh_segment(img, k, 0.8, min_size, strict=False, f32=True)
            print(f"  {name:8s} min_size={min_size:2d} k={k:5d}: regions {n_regions(ref):5d} | "
                  f"same as skimage: float32 {same_partition(ref, mine32)}, "
                  f"float64 {same_partition(ref, mine64)} | "
                  f"'<=' instead of '<' same: {same_partition(mine32, le)}")

A. NumPy implementation against skimage.segmentation.felzenszwalb (sigma 0.8)
  coins    min_size= 1 k=  100: regions  1252 | same as skimage: float32 True, float64 True | '<=' instead of '<' same: True
  coins    min_size= 1 k=  300: regions   510 | same as skimage: float32 True, float64 True | '<=' instead of '<' same: True
  coins    min_size= 1 k= 1000: regions   129 | same as skimage: float32 True, float64 True | '<=' instead of '<' same: True
  coins    min_size=20 k=  100: regions   295 | same as skimage: float32 True, float64 True | '<=' instead of '<' same: True
  coins    min_size=20 k=  300: regions    74 | same as skimage: float32 True, float64 True | '<=' instead of '<' same: True
  coins    min_size=20 k= 1000: regions    52 | same as skimage: float32 True, float64 True | '<=' instead of '<' same: True
  chelsea  min_size= 1 k=  100: regions  1145 | same as skimage: float32 True, float64 True | '<=' instead of '<' same: True
  chelsea  min_size= 1 k=  300: regions   643 | same as skimage: float32 True, float64 True | '<=' instead of '<' same: True
  chelsea  min_size= 1 k= 1000: regions   343 | same as skimage: float32 True, float64 True | '<=' instead of '<' same: True
  chelsea  min_size=20 k=  100: regions   287 | same as skimage: float32 True, float64 True | '<=' instead of '<' same: True
  chelsea  min_size=20 k=  300: regions   128 | same as skimage: float32 True, float64 True | '<=' instead of '<' same: True
  chelsea  min_size=20 k= 1000: regions    51 | same as skimage: float32 True, float64 True | '<=' instead of '<' same: True

The plain implementation reproduces scikit-image’s partition exactly in all twelve cases, with and without the minimum-size pass, for a grey and a colour image, so the library is the paper’s Algorithm 1, on one colour-distance graph, plus the post-processing step. The two open details make no difference here: an edge weight has to equal $\mathrm{Int}(C) + k/|C|$ exactly for $<$ and $\le$ to disagree, or lie within single-precision rounding of it for float32 to matter, and on these images no such case changes the partition. The minimum size changes the result far more than either: at $k = 100$ it reduces the coins segmentation from 1,252 to 295 regions.

The rest of the script continues in the same file. Part B builds an image after the paper’s Figure 1, $240 \times 160$ pixels: a ramp from 0.1 to 0.5 on the left half, a constant 0.55 on the right, and a $60 \times 60$ square of uniform noise in $[0.7, 1.0]$ inside the constant half. It scores FH against a fixed threshold on edge weights (the connected components of all edges with $w \le t$, the kind of local rule the paper argues against), by the best intersection over union any predicted region reaches with each true region. Parts C–E use scikit-image, now shown to be equivalent, for the effect of $k$, the minimum size and the resolution, for nesting across scales, and for Theorem 3.

print("\nB. Synthetic image after the paper's Figure 1 (ramp | constant with noisy square)")
rng = np.random.default_rng(0)
H, W = 160, 240
img = np.empty((H, W))
img[:, :W // 2] = np.linspace(0.1, 0.5, W // 2)[None, :]
img[:, W // 2:] = 0.55
gt = np.zeros((H, W), int)
gt[:, W // 2:] = 1
r0, c0 = H // 2 - 30, W // 2 + W // 4 - 30
img[r0:r0 + 60, c0:c0 + 60] = rng.uniform(0.7, 1.0, (60, 60))
gt[r0:r0 + 60, c0:c0 + 60] = 2


def best_iou(lab, gt):
    scores = []
    for g in range(3):
        G = gt == g
        cand = np.unique(lab[G])
        scores.append(max(np.sum(G & (lab == c)) / np.sum(G | (lab == c)) for c in cand))
    return np.array(scores)


for k in (100, 300, 1000):
    lab = fh_segment(img, k, 0.8, 1)
    s = best_iou(lab, gt)
    print(f"  FH k={k:5d}: regions {n_regions(lab):5d}, best IoU ramp/constant/noise "
          f"{s[0]:.3f}/{s[1]:.3f}/{s[2]:.3f}, mean {s.mean():.3f}")
n, ea, eb, wt = grid_graph(img, 0.8)
for t in (0.5, 1, 2, 4, 8, 16, 32):
    keep = wt <= t / 255.0
    g = coo_matrix((np.ones(keep.sum()), (ea[keep], eb[keep])), shape=(n, n))
    lab = connected_components(g, directed=False)[1].reshape(H, W)
    s = best_iou(lab, gt)
    print(f"  fixed threshold t={t:4}: regions {n_regions(lab):5d}, best IoU "
          f"{s[0]:.3f}/{s[1]:.3f}/{s[2]:.3f}, mean {s.mean():.3f}")

print("\nC. Scale parameter k and region count, coins (skimage, sigma 0.8)")
half = transform.rescale(coins, 0.5, anti_aliasing=True)
print("      k   min_size=1  min_size=20  half-res min_size=1")
for k in (10, 30, 100, 300, 1000, 3000):
    c1 = n_regions(felzenszwalb(coins, scale=k, sigma=0.8, min_size=1, channel_axis=None))
    c20 = n_regions(felzenszwalb(coins, scale=k, sigma=0.8, min_size=20, channel_axis=None))
    ch = n_regions(felzenszwalb(half, scale=k, sigma=0.8, min_size=1, channel_axis=None))
    print(f"  {k:5d}   {c1:10d}  {c20:11d}  {ch:19d}")

print("\nD. Are segmentations at increasing k nested? coins, k = 50, 100, ..., 1000, min_size=1")
ks = list(range(50, 1001, 50))
labs = [felzenszwalb(coins, scale=k, sigma=0.8, min_size=1, channel_axis=None) for k in ks]
counts = [n_regions(l) for l in labs]
rises = sum(c2 > c1 for c1, c2 in zip(counts, counts[1:]))
not_nested = 0
for fine, coarse in zip(labs, labs[1:]):
    pairs = np.unique(fine.ravel().astype(np.int64) * (coarse.max() + 1) + coarse.ravel()).size
    not_nested += pairs > n_regions(fine)
print(f"  region counts: {counts}")
print(f"  steps where the count rises: {rises} of {len(ks) - 1}")
print(f"  consecutive pairs where a finer region is split by the coarser one: "
      f"{not_nested} of {len(ks) - 1}")

print("\nE. Theorem 3: tie order, coins as integers (sigma 0, k = 300 on the 0-255 scale)")
n, ea, eb, wt = grid_graph(coins.astype(np.float64), 0)
print(f"  edges {wt.size}, distinct weights {np.unique(wt).size}")
for min_size in (1, 20):
    base = segment_graph(n, ea, eb, wt, 300.0, np.argsort(wt, kind="stable"), min_size)
    diff = 0
    for seed in (0, 1, 2):
        key = np.random.default_rng(seed).random(wt.size)
        lab = segment_graph(n, ea, eb, wt, 300.0, np.lexsort((key, wt)), min_size)
        diff += not same_partition(base, lab)
    print(f"  min_size={min_size:2d}: {n_regions(base)} regions; shuffled tie orders "
          f"giving a different partition: {diff} of 3")

B. Synthetic image after the paper's Figure 1 (ramp | constant with noisy square)
  FH k=  100: regions   183, best IoU ramp/constant/noise 0.017/0.948/0.900, mean 0.621
  FH k=  300: regions    31, best IoU ramp/constant/noise 0.983/0.948/0.951, mean 0.961
  FH k= 1000: regions    25, best IoU ramp/constant/noise 0.983/0.948/0.951, mean 0.961
  fixed threshold t= 0.5: regions  3242, best IoU 0.008/0.959/0.002, mean 0.323
  fixed threshold t=   1: regions  2415, best IoU 0.983/0.974/0.004, mean 0.654
  fixed threshold t=   2: regions  1350, best IoU 0.992/0.974/0.026, mean 0.664
  fixed threshold t=   4: regions   282, best IoU 1.000/0.985/0.814, mean 0.933
  fixed threshold t=   8: regions    23, best IoU 0.552/0.448/0.979, mean 0.660
  fixed threshold t=  16: regions     2, best IoU 0.552/0.448/1.000, mean 0.667
  fixed threshold t=  32: regions     1, best IoU 0.500/0.406/0.094, mean 0.333

C. Scale parameter k and region count, coins (skimage, sigma 0.8)
      k   min_size=1  min_size=20  half-res min_size=1
     10        10010          885                 3079
     30         3169          618                  871
    100         1252          295                  229
    300          510           74                  111
   1000          129           52                   54
   3000          120           26                   45

D. Are segmentations at increasing k nested? coins, k = 50, 100, ..., 1000, min_size=1
  region counts: [2081, 1252, 835, 638, 559, 510, 337, 314, 289, 252, 235, 220, 186, 176, 144, 135, 123, 126, 130, 129]
  steps where the count rises: 2 of 19
  consecutive pairs where a finer region is split by the coarser one: 19 of 19

E. Theorem 3: tie order, coins as integers (sigma 0, k = 300 on the 0-255 scale)
  edges 463349, distinct weights 204
  min_size= 1: 3040 regions; shuffled tie orders giving a different partition: 0 of 3
  min_size=20: 433 regions; shuffled tie orders giving a different partition: 3 of 3

The adaptive predicate. At $k = 300$ and $k = 1000$ FH finds all three regions, with IoUs between 0.948 and 0.983. The other regions do not affect this score; smoothing creates intermediate values along each boundary, which on grid graphs the paper finds produces narrow boundary regions. At $k = 100$ the ramp falls apart, and the arithmetic of $\tau$ explains why. The ramp rises by $0.4/119 \times 255 \approx 0.86$ grey levels per column and is constant down each column, so the zero-weight vertical edges come first in the order and make each column a 160-pixel component with $\mathrm{Int} = 0$. The horizontal edges of 0.86 then need $0.86 \le k/160$: false at $k = 100$ ($0.63$), true at $k = 300$ ($1.88$). The best IoU of 0.017 is that of a two-column strip with the 120-column ramp. Once two columns merge, their internal difference becomes 0.86 and every later merge across the ramp passes, which is why the ramp jumps from strips to one region between the two settings.

The fixed threshold does better than the plan for this example expected (“no fixed threshold exceeds about 0.75”): at $t = 4$ it scores 0.933 with 282 regions. The paper’s dichotomy still holds in the strict sense, since the square is split (its best region has IoU 0.814) while the ramp and the constant region stay apart, but the split is mild: in our reading, after smoothing with $\sigma = 0.8$, enough edges inside the square fall below $t = 4$ to connect most of it. Whether the paper’s own Figure 1 image would behave the same cannot be checked here. The window is narrow, though. At $t = 2$ the square shatters and at $t = 8$ the ramp and the constant region merge, so both score about 0.66, while FH gives the same IoUs at $k = 300$ and 1000 with about a tenth as many regions as the best threshold.

Scale, minimum size and resolution. The region count falls steeply with $k$, from about 10,000 to about 120 over $k = 10$ to 3,000, but not in proportion to it, and it barely changes from $k = 1000$ to 3,000. The minimum-size column shows why: most of the regions left at large $k$ are smaller than 20 pixels, since the minimum size removes them (120 regions become 26 at $k = 3000$). In our reading they survive because the step to their large neighbour exceeds that neighbour’s $\mathrm{Int} + k/|C|$, which is the paper’s point that $k$ is not a minimum size. At $k = 10$ the minimum size cuts the count by an order of magnitude. The half-resolution image gives fewer regions than the full-resolution one at every $k$, which in our reading follows from each object having a quarter of the pixels and so four times the tolerance; settings tuned at one image size do not transfer to another.

Not a hierarchy. On the dense grid of $k$ the count rises twice ($k = 850$ to 900 to 950), and in every one of the 19 consecutive pairs some region of the finer segmentation is split between regions of the coarser one. Raising $k$ does not only merge regions; it moves boundaries. In our reading, components grow earlier at a larger $k$, so their size and internal difference differ when later edges are examined, and a pixel strip that one component absorbed at a small $k$ can be taken by its neighbour at a larger one. Guimarães et al. (2017) describe the same behaviour as violations of the causality and location principles of multiscale analysis.

Theorem 3. The integer image has only 204 distinct weights among 463,349 edges, so ties are everywhere. Without post-processing, the three random tie orders give exactly the stable-sort partition, as Theorem 3 guarantees. With a minimum size of 20, all three differ: the post-processing pass merges small components in edge order, and the theorem does not cover it. In our reading, a library whose sort does not keep ties in a fixed order (NumPy’s default argsort kind is not stable) could therefore give slightly different minimum-size results across sort implementations.

Properties and Behavior

Guarantees. The output is neither too fine nor too coarse with respect to $D$ (Theorems 1 and 2) and independent of the tie order (Theorem 3). These are statements about consistency with the predicate, not about agreement with people: they say the greedy result could not be split or merged without contradicting $D$, not that $D$ captures perceptual boundaries. The guarantees hold for Algorithm 1 alone; the minimum-size post-processing, which most users run, is outside them.

Preserving detail where variation is low. The abstract names this as an important characteristic: detail is kept in low-variability regions and ignored in high-variability ones. The paper’s indoor scene keeps name tags and objects behind windows as separate regions while grouping a textured air-conditioning duct, clothing and furniture into large regions; in the street scene, specular reflections on a van are treated as internal variation. Slow intensity changes from lighting, as across the COIL objects, do not split a region.

Sensitivity to one weak edge. $\mathrm{Dif}$ is the minimum edge between two regions, so a single low-weight edge, for example a gap of similar colour through a faint boundary, is enough to merge them. The authors acknowledge this limitation and point out that the edge is compared only with minimum spanning tree edges inside each component; the robust alternative, a quantile, is NP-hard, as the Mathematical Formulation section shows. A related leak appears in the paper’s baseball scene, taken from the normalized cuts paper: a large grassy region absorbed part of the wall, which the authors attribute to the high variation inside the region and to a long, slow change in intensity between grass and wall, similar in magnitude to the folds in the players’ uniforms.

Boundary regions. On grid graphs, narrow regions one or two pixels wide appear along boundaries where pixels are halfway between the two sides in colour; the paper calls this common to any grid-graph method and suggests removing long thin regions whose colour is close to the average of their neighbours. The thin strips in the example above are the same effect.

Irregular shapes and no control of the count. The predicate has no term for shape or size beyond $\tau$, so regions are as irregular as the image and their number is set only indirectly by $k$, as scikit-image’s documentation also warns. The superpixels article compares FH with methods such as SLIC that give compact cells and control their number.

Accuracy on benchmarks. Arbeláez, Maire, Fowlkes and Malik (2011) counted FH among the three algorithms most widely used as sources of image segments, because of reasonable performance and public code, and evaluated it over a range of parameters on the Berkeley Segmentation Dataset. On the 200 test images of BSDS500 (their Tables 1 and 2, authors’ copy), FH reached a boundary F-measure of 0.61 at the optimal dataset scale (ODS) and 0.64 at the optimal image scale (OIS), with average precision 0.56, against 0.64, 0.68 and 0.56 for mean shift, 0.64, 0.68 and 0.45 for multiscale normalized cuts, and 0.73, 0.76 and 0.73 for the authors’ gPb-owt-ucm. On the region benchmarks its segmentation covering was 0.52 (ODS), 0.57 (OIS) and 0.69 (best over all parameters), against 0.54, 0.58 and 0.66 for mean shift.

Limitations

• Greedy and local in the end. The global properties are relative to $D$, and $D$ depends on a single minimum edge; leaks through weak boundaries are the typical failure.
• One flat partition per $k$. Partitions at different $k$ are not nested and the region count is not monotonic in $k$, which complicates choosing a scale and using the method for multiscale analysis.
• Resolution-dependent parameters. $k/|C|$ and the minimum size count pixels, so the same settings give different results at different image sizes.
• Weak cues. The grid-graph weights compare single smoothed pixels, so texture boundaries between two regions of similar mean colour are invisible to it, and textured regions are grouped only when their neighbour differences are small enough relative to $k$.
• Unspecified post-processing. The minimum-size pass that released code applies changes the result substantially, depends on the edge order, and is not part of the paper’s analysis.
• Colour handling differs from the paper. Results from common libraries come from one graph on colour distance, not the paper’s intersection of three per-channel segmentations, which the authors found to work better.

Variants

Nearest-neighbour graphs. The paper’s own second variant replaces the grid by the ten nearest neighbours in $(x, y, r, g, b)$. It gave very similar results to the grid graph on the street and baseball scenes, and it can group regions that are not spatially connected: the red flowers of a garden scene, disjoint in the foreground, mostly ended up in one region, and the lights of the Eiffel Tower at night formed one disconnected region.

Other threshold functions. Any non-negative $\tau$ keeps the algorithm and its guarantees, which the paper suggests using for shape preferences. The quantile version of $\mathrm{Dif}$, by contrast, loses tractability.

Hierarchical versions. Guimarães, Kenmochi, Cousty, Patrocinio and Najman (2017) built hGB, a hierarchy in which each level is meant to be the largest segmentation that is not too coarse at the corresponding scale, chosen among the partitions of a quasi-flat zones hierarchy (their algorithm is greedy and does not guarantee it), so that partitions are nested, the region count decreases with scale and any number of regions can be read off directly. Their FH baseline used 8-adjacency, Euclidean RGB weights and an area filter tuned on the training and validation folds of BSDS500. On the 200 test images, hGB’s average optimal region F-measure was 0.604 against 0.575 for FH under the OIS setting and 0.492 against 0.477 under ODS; for the boundary F-measure under OIS the paired difference was not significant (confidence interval $[-0.0179, 0.00717]$, $p = 0.21$). They computed the whole hierarchy in 1.37 s for a $321 \times 481$ image, against 7.84 s for 50 separate FH runs.

Starting regions for proposals. Selective search (Uijlings et al., 2013) uses FH to produce the initial regions that it merges hierarchically into region proposals; the authors chose it as, to their knowledge, the fastest publicly available method giving high-quality starting regions. It diversifies the starting regions by running FH in several colour spaces and at several values of $k$, and found that a lower $k$, giving smaller starting regions, gave better proposals at the cost of more boxes. Its released Matlab code uses $\sigma = 0.8$ and a minimum size equal to $k$. The proposals of selective search in turn fed the R-CNN detector.

Other graph methods. Graph cuts and normalized cuts also segment a pixel graph, but by optimizing a cut criterion, globally for a binary labelling with seeds or approximately by eigenvectors; FH optimizes nothing and instead guarantees properties of the greedy result. Its merging order is Kruskal’s, which is also the order of single-linkage clustering; in our reading FH is single linkage with a stopping rule that adapts to each region. The watershed transform is another fast way to oversegment an image, and Guimarães et al. note that it too is related to a minimum spanning tree.