Graph Cuts

Graph cuts minimize an energy over pixel labelings by computing a minimum cut in a graph built from the image: pixels are nodes, two terminal nodes stand for the labels, and edge weights encode how well each pixel fits each label and what it costs neighbors to disagree. For two labels and a broad class of energies the cut is the exact global minimum; for many labels, sequences of binary cuts give approximations, with a guaranteed bound for expansion moves. The method has been used for interactive segmentation, stereo, image restoration and texture synthesis, and is the optimizer inside GrabCut. The image segmentation hub places it among the classical methods.

Problem

Many vision problems assign a label $L_p$ to every pixel $p$ (object or background, a disparity, an intensity) that fits the data and agrees with its neighbors, so that noise does not create isolated wrong pixels. Both demands go into one energy:

$$E(L) = \sum_{p \in \mathcal{P}} D_p(L_p) + \lambda \sum_{\{p,q\} \in \mathcal{N}} V_{pq}(L_p, L_q),$$

where $\mathcal{P}$ is the set of pixels, $\mathcal{N}$ the set of neighboring pairs, $D_p$ the data term, $V_{pq}$ the smoothness term, and $\lambda \ge 0$ balances them (Boykov and Jolly put the weight on the data term instead).

If $D_p(l) = -\log P(I_p \mid l)$ is the negative log-likelihood of the observed intensity and the smoothness term is the negative log of a Markov random field (MRF) prior favoring equal neighbors, minimizing $E$ gives the maximum a posteriori (MAP) labeling. With two labels and $V_{pq}(a,b) = [a \ne b]$ the prior is the Ising model; with more labels, the Potts model. Greig, Porteous and Seheult (1989) showed that the exact MAP estimate of a degraded binary image can be computed with a maximum-flow algorithm, and their experiments suggested that simulated annealing with practicable temperature schedules can give poor approximations to it.

Inputs and Outputs

• Input: data costs $D_p(l)$, smoothness costs $V_{pq}$, a neighborhood system (4- or 8-connected in 2D, 6- or 26-connected in 3D), $\lambda$, and, for interactive use, seed pixels with fixed labels.
• Output: a label per pixel: a global minimum of $E$ for binary problems with submodular pairwise terms, a local minimum with respect to large moves for multi-label problems.

Intuition

Think of the graph as pipes. Each pixel has a pipe from the source, wide if the pixel looks like foreground, a pipe to the sink, wide if it looks like background, and pipes to its neighbors, wide if they are similar. The narrowest set of pipes that disconnects source from sink separates foreground from background. Cutting between similar neighbors is expensive, so the cut follows image edges and ignores isolated noisy pixels unless their evidence is strong.

Algorithm

For labels 1 (foreground, source side) and 0 (background, sink side):

1. Create a node per pixel, plus a source $s$ and a sink $t$.
2. Add t-links (terminal links): $s \to p$ with capacity $D_p(0)$ and $p \to t$ with capacity $D_p(1)$, shifted by a per-pixel constant if needed to make them nonnegative.
3. Add n-links (neighbor links): $p \to q$ and $q \to p$, each with capacity $\lambda w_{pq}$, for a weighted Potts term $V_{pq}(a,b) = w_{pq}[a \ne b]$.
4. Compute a maximum flow from $s$ to $t$.
5. Label 1 the nodes reachable from $s$ in the residual graph (along edges with spare capacity), and 0 the rest.

An $s$–$t$ cut splits the nodes into a set $S \ni s$ and a set $T \ni t$ and costs the total capacity of edges from $S$ to $T$. A pixel in $T$ severs $s \to p$ and pays $D_p(0)$, a pixel in $S$ severs $p \to t$ and pays $D_p(1)$, and each neighboring pair split by the cut pays $\lambda w_{pq}$. So every labeling is a cut whose cost equals its energy. Ford and Fulkerson (1956) proved that the maximum flow value equals the minimum cut capacity; step 5 recovers that cut from the flow.

Mathematical Formulation

Which binary energies a cut can minimize

Kolmogorov and Zabih (2004) showed that an energy of binary variables made of terms on one and two variables can be represented by a graph cut exactly when every pairwise term satisfies

$$V_{pq}(0,0) + V_{pq}(1,1) \le V_{pq}(0,1) + V_{pq}(1,0).$$

They called such terms regular; it is submodularity for binary functions: agreeing costs no more in total than disagreeing. Their construction splits each term into a constant, two single-variable parts that become t-links, and one n-link of capacity $V_{pq}(0,1) + V_{pq}(1,0) - V_{pq}(0,0) - V_{pq}(1,1)$, nonnegative exactly when the inequality holds. They proved that energies with a non-regular pairwise term are NP-hard to minimize in general, and extended the result to terms on three variables. Potts terms satisfy the inequality, so binary segmentation is solved exactly.

Interactive segmentation

In Boykov and Jolly’s (2001) interactive segmentation of images and volumes, pixels the user marks as “object” or “background” are seeds, imposed as hard constraints: an object seed gets a t-link to the source of capacity $K = 1 + \max_p \sum_q B_{pq}$, more than all its n-links together, and none to the sink, so no minimum cut puts it on the wrong side. The cut is then the global minimum among all segmentations consistent with the seeds. Regional penalties were negative log-likelihoods under intensity histograms of the seed pixels, and the boundary penalty an ad hoc function

$$B_{pq} \propto \exp\!\left( -\frac{(I_p - I_q)^2}{2\sigma^2} \right) \cdot \frac{1}{\operatorname{dist}(p,q)},$$

large between similar pixels and small across edges, with $\sigma$ estimated as the camera noise. They used 8-connected neighborhoods in 2D and 26-connected in 3D, and segments could have any topology. A new seed changes only two t-links, and the new maximum flow can start from the previous one, which kept corrections interactive.

Multiple labels: expansion and swap moves

With more than two labels, minimizing even the Potts energy is NP-hard (Boykov, Veksler & Zabih, 2001). Their algorithms make large moves, each a binary problem solved by a minimum cut:

• An α-β swap lets any subset of pixels labeled α or β exchange those labels; it requires $V$ to be a semimetric ($V(a,b) = V(b,a) \ge 0$, zero exactly when $a = b$).
• An α-expansion lets any subset of pixels switch to α; it requires a metric, adding the triangle inequality, which Kolmogorov and Zabih later showed makes each expansion’s binary energy regular.

The algorithm cycles over labels (or label pairs), accepting moves that lower the energy, and stops after a cycle without improvement; in the authors’ experiments that took a few cycles. With nonnegative data terms, as they assume, any labeling $\hat{f}$ that no expansion move improves satisfies

$$E(\hat{f}) \le 2c\, E(f^*), \qquad c = \frac{\max_{a \ne b} V(a,b)}{\min_{a \ne b} V(a,b)},$$

where $f^*$ is the global minimum; for Potts, $c = 1$ and the bound is a factor of two. A local minimum for swap moves can be arbitrarily far from the global one.

Parameters

• Smoothness weight $\lambda$. Zero gives independent per-pixel decisions; too large a value deletes thin structures and eventually the object (see the example). In Boykov and Jolly’s synthetic test, where their $\lambda$ weights the regional term, one result held for $\lambda$ from 7 to 43, so good values can form a broad range.
• Data model. Histograms, Gaussians or Gaussian mixtures, from seeds or earlier data.
• Contrast scale $\sigma$ and neighborhood size, which trades memory for less grid bias.

Complexity

The graph has $n = |\mathcal{P}| + 2$ nodes and $m = 2|\mathcal{P}| + 2|\mathcal{N}|$ edges, so memory is linear in pixels times neighborhood size. Max-flow algorithms either augment flow along source-to-sink paths or, like Goldberg and Tarjan’s (1988) push–relabel method, push excess flow locally. Boykov and Kolmogorov (2004) list worst-case bounds of $O(mn^2)$ for Dinic’s augmenting-path algorithm and $O(n^3)$ and $O(n^2\sqrt{m})$ for two push–relabel variants. Their own augmenting-path algorithm grows search trees from both terminals and reuses them rather than rebuilding them. Its only bound, $O(mn^2|C|)$ with $|C|$ the minimum cut cost, is not polynomial, yet in their tests it was 2–5 times faster than the others on most vision graphs and consistently several times faster on 2D grids; on 3D segmentation graphs push–relabel was comparable.

Implementation

The NumPy and SciPy code builds the graph for a binary Potts energy on a 4-connected grid and solves it with scipy.sparse.csgraph.maximum_flow (SciPy 1.16, Dinic’s algorithm by default), which requires integer capacities, so costs are scaled by 1,000 and rounded. It first checks optimality on a $3 \times 4$ image with random costs against all $2^{12} = 4{,}096$ labelings, then segments a noisy $32 \times 32$ disk with a thin bar (2 by 10 pixels), with Gaussian negative log-likelihoods as data terms.

import itertools
import numpy as np
from scipy.sparse import csr_array
from scipy.sparse.csgraph import maximum_flow, breadth_first_order

SCALE = 1000                                       # SciPy's max-flow needs integer capacities

def graph_cut(D0, D1, lam):
    """Minimize sum_p D_{L_p}(p) + lam * sum_{p~q} [L_p != L_q] on a 4-connected grid.
    Returns the label map (1 = source side = foreground) and the max-flow value."""
    H, W = D0.shape
    n = H * W
    s, t = n, n + 1
    idx = np.arange(n).reshape(H, W)
    rows, cols, caps = [], [], []
    def add(u, v, c):
        rows.append(u.ravel()); cols.append(v.ravel()); caps.append(np.ravel(c))
    add(np.full(n, s), idx, D0)                    # t-link s->p, cut if p is background
    add(idx, np.full(n, t), D1)                    # t-link p->t, cut if p is foreground
    for a, b in [(idx[:, :-1], idx[:, 1:]), (idx[:-1, :], idx[1:, :])]:
        w = np.full(a.size, lam)                   # n-links in both directions (Potts)
        add(a, b, w); add(b, a, w)
    r, c = np.concatenate(rows), np.concatenate(cols)
    cap = np.rint(np.concatenate(caps) * SCALE).astype(np.int32)
    G = csr_array((cap, (r, c)), shape=(n + 2, n + 2))
    res = maximum_flow(G, s, t)
    residual = G - res.flow                        # residual capacity of every edge
    residual.data[residual.data < 0] = 0
    residual.eliminate_zeros()
    reach = breadth_first_order(residual, s, return_predecessors=False)
    L = np.zeros(n + 2, int); L[reach] = 1         # reachable from s = source side
    return L[:n].reshape(H, W), res.flow_value / SCALE

def energy(L, D0, D1, lam):                        # same rounding as the graph
    q = lambda x: np.rint(np.asarray(x) * SCALE)
    data = np.where(L == 1, q(D1), q(D0)).sum()
    cuts = (L[:, 1:] != L[:, :-1]).sum() + (L[1:, :] != L[:-1, :]).sum()
    return (data + cuts * q(lam)) / SCALE

# 1. Optimality check on a 3 x 4 image: brute force over all 2^12 = 4,096 labelings.
rng = np.random.default_rng(1)
D0, D1 = rng.uniform(0, 2, (3, 4)), rng.uniform(0, 2, (3, 4))
for lam in [0.3, 1.0]:
    L, flow = graph_cut(D0, D1, lam)
    best = min(energy(np.array(b).reshape(3, 4), D0, D1, lam)
               for b in itertools.product([0, 1], repeat=12))
    print(f"lambda={lam}: max flow {flow:.3f}, energy of cut labeling {energy(L, D0, D1, lam):.3f}, "
          f"brute-force minimum {best:.3f}")

# 2. A noisy 32 x 32 binary image: a disk with a thin bar, foreground 0.7, background 0.3.
rng = np.random.default_rng(0)
H = W = 32
yy, xx = np.mgrid[0:H, 0:W]
gt = ((yy - 16)**2 + (xx - 12)**2 <= 7**2).astype(int)
bar = np.zeros_like(gt); bar[15:17, 19:29] = 1     # bar, 2 pixels wide and 10 long
gt |= bar
mu0, mu1, sigma = 0.3, 0.7, 0.35
img = np.where(gt == 1, mu1, mu0) + rng.normal(0, sigma, (H, W))
D0 = (img - mu0)**2 / (2 * sigma**2)               # negative log-likelihoods,
D1 = (img - mu1)**2 / (2 * sigma**2)               # shared constant dropped

def report(name, L):
    iou = (L & gt).sum() / (L | gt).sum()
    print(f"{name:<18} errors {int((L != gt).sum()):4d}  IoU {iou:.3f}  "
          f"foreground {int(L.sum()):4d} px  bar pixels kept {int((L & bar).sum()):2d}/20")

print(f"ground truth: {int(gt.sum())} foreground pixels")
report("threshold 0.5", (img > 0.5).astype(int))
for lam in [0.0, 0.25, 0.5, 1.0, 1.5, 2.0]:
    report(f"graph cut lam={lam}", graph_cut(D0, D1, lam)[0])

Output:

lambda=0.3: max flow 10.455, energy of cut labeling 10.455, brute-force minimum 10.455
lambda=1.0: max flow 10.748, energy of cut labeling 10.748, brute-force minimum 10.748
ground truth: 168 foreground pixels
threshold 0.5      errors  270  IoU 0.304  foreground  338 px  bar pixels kept 15/20
graph cut lam=0.0  errors  270  IoU 0.304  foreground  338 px  bar pixels kept 15/20
graph cut lam=0.25 errors  148  IoU 0.477  foreground  250 px  bar pixels kept 14/20
graph cut lam=0.5  errors   53  IoU 0.735  foreground  179 px  bar pixels kept 13/20
graph cut lam=1.0  errors   36  IoU 0.794  foreground  146 px  bar pixels kept  6/20
graph cut lam=1.5  errors   58  IoU 0.657  foreground  112 px  bar pixels kept  0/20
graph cut lam=2.0  errors  168  IoU 0.000  foreground    0 px  bar pixels kept  0/20

The maximum flow, the energy of the cut’s labeling and the brute-force minimum agree. With $\lambda = 0$ the graph cut is thresholding at 0.5, the per-pixel MAP decision for two equal-variance Gaussians: 270 of 1,024 pixels are wrong and the foreground IoU is 0.304. At $\lambda = 1$ only 36 are wrong (IoU 0.794).

Larger $\lambda$ overrides correct evidence. It follows from the energy that a foreground pixel gains on average $(\mu_1 - \mu_0)^2 / 2\sigma^2 \approx 0.65$ in data cost, while each unit of boundary costs $\lambda$. A bar two pixels wide gains about 1.3 per unit length but pays $2\lambda$ for its two sides, so it should survive only below $\lambda \approx 0.65$; it is mostly gone at 1 and entirely at 1.5. The disk alone (149 pixels, 60 boundary edges) gains about 97 against $60\lambda$, so above $\lambda \approx 1.6$ the empty labeling is cheaper in expectation, and at $\lambda = 2$ the object vanishes. Over noise seeds 0–9, IoU at $\lambda = 1$ ranged from 0.66 to 0.85 against 0.27–0.31 for thresholding, and at $\lambda = 2$ the empty labeling won in 6 of the 10 runs.

Limitations

• Shrinking bias. Charging for boundary length favors short boundaries, cutting off thin structures and shrinking objects when the data term is weak; Boykov and Jolly showed an object shrinking with the boundary term alone.
• Restricted energies. Exact solutions need submodular binary energies or, as Boykov, Veksler and Zabih note of earlier multi-label constructions, convex $V$, which over-penalizes discontinuities; the rest is approximate.
• Grid bias. It follows from counting cut edges that a 4-connected grid measures boundary length in the city-block metric, so diagonal boundaries cost more; larger neighborhoods reduce this.
• Fixed appearance models. The data term must be known before the cut; seeds supply it, and GrabCut re-estimates it.
• Large graphs. 3D graphs are large; Boykov and Jolly reported minutes for $512 \times 512 \times 50$ volumes, against under a second for most 2D images up to $512 \times 512$ (on a 333 MHz Pentium III).

Variants

• GrabCut (Rother, Kolmogorov & Blake, 2004) starts from a rectangle around the object, models foreground and background colors with full-covariance Gaussian mixtures (typically 5 components each), and alternates between assigning pixels to components, refitting the mixtures, and a minimum cut. Each step lowers the energy, so it converges to a local minimum; border matting then estimates transparency along the boundary.
• Stereo and multi-view. Boykov, Veksler and Zabih applied their moves to stereo matching; in the Scharstein–Szeliski evaluation graph cuts were the most accurate of the tested optimizers and among the slowest, and semi-global matching approximates such global energies much faster. Volumetric graph cuts appear in multi-view stereo.
• Restoration. Greig et al.’s noisy binary images are a special case of image denoising; Boykov, Veksler and Zabih restored images with many intensity labels.
• Seam finding. Graphcut textures (Kwatra et al., 2003) stitch copied patches along minimum-cut seams, which, unlike dynamic programming, works in any dimension, including video.
• Normalized cuts partition a similarity graph without terminals, normalizing each cut by the segments’ total connections to counter the minimum cut’s bias toward small segments; the exact problem is NP-complete and solved approximately.