Markov Random Field

A Markov random field (MRF), also called a Markov network or an undirected graphical model, is a joint probability distribution over many variables arranged on a graph, in which each variable, given its neighbours, is independent of all the others. In vision the variables are usually the pixels of an image or the labels assigned to them (object or background, a disparity, an intensity), and the graph is the pixel grid. The model says that a pixel’s value is predicted by the values around it, and by nothing further away once those are known.

Two results make the model usable. The Hammersley–Clifford theorem turns the local statement into a global formula, a Gibbs distribution $p(x) \propto \exp(-E(x))$ whose energy is a sum of terms over small groups of neighbours. And Bayes’ rule (see probability and statistics) turns an MRF prior and a per-pixel noise model into a posterior that is again an MRF, so that restoring or labelling an image becomes minimizing an energy of a data term plus a smoothness term. Geman and Geman (1984) built image restoration on these two facts. Szeliski et al. (2008) write that pixel-labelling problems are naturally posed as energy minimization, justified as MAP estimation of an MRF, and benchmark optimizers on stereo, stitching, segmentation and denoising energies. Graph cuts own the labelling energy and its exact minimization for two labels; this article covers the model behind it, its common forms, what can be computed from it and the other optimizers used on it, and compares them in an example.

Definition

Let $\mathcal{S}$ be a set of sites (pixels), and give every site $s$ a set of neighbours $\mathcal{N}(s) \subset \mathcal{S}$ , with $s \notin \mathcal{N}(s)$ and $t \in \mathcal{N}(s)$ exactly when $s \in \mathcal{N}(t)$ . A random field $X = \{X_s\}$ is a Markov random field with respect to this neighbourhood system if, for every site and every configuration $x$ ,

P(X_s = x_s \mid X_t = x_t,\ t \ne s) = P(X_s = x_s \mid X_t = x_t,\ t \in \mathcal{N}(s)).

This is Geman and Geman’s (1984, Eq. 1.1, p. 722) definition, written for the pixel lattice. The neighbourhood is part of the model: the usual choices on an image are the 4-connected grid (left, right, up, down) and the 8-connected grid. A clique is a set of sites that are all neighbours of one another; on the 4-connected grid the cliques are single pixels and horizontal or vertical pairs.

Intuition

Think of the image as a sheet of small magnets that each point up or down, where each magnet feels only its four nearest neighbours. Left to themselves, neighbouring magnets prefer to agree, so the sheet forms patches. Observing a noisy photograph adds a second pull on each magnet, towards the value the photograph suggests at that pixel. The most probable arrangement balances the two pulls: isolated disagreements, which noise creates, cost more than they explain, while a large region that the data support survives even though its boundary costs something. Geman and Geman (1984, p. 721) made this analogy with statistical mechanics explicit, treating grey levels and edges as states of atoms in a lattice, and borrowed from physics the idea of slowly lowering a temperature until the system settles into a low-energy state.

Formal Definition

Markov properties

Murphy (2023, §4.3.6) states three conditional-independence properties of an undirected graph $G$ . The global Markov property: sets of variables $A$ and $B$ are independent given a set $C$ whenever $C$ separates $A$ from $B$ in $G$ . The local property: each variable is independent of all others given its neighbours, its Markov blanket (the definition above). The pairwise property: two variables not joined by an edge are independent given all the rest. Global implies local implies pairwise, and for distributions that are strictly positive everywhere the three are equivalent (Murphy, 2023, §4.3.6.1). On the 4-connected grid the Markov blanket of a pixel is its four neighbours.

Gibbs distribution and the Hammersley–Clifford theorem

A Gibbs distribution on the graph is

p(x) = \frac{1}{Z} \prod_{c \in \mathcal{C}} \psi_c(x_c) = \frac{1}{Z} \exp\!\big(-E(x)\big), \qquad E(x) = \sum_{c \in \mathcal{C}} E_c(x_c), \qquad Z = \sum_x \prod_{c} \psi_c(x_c),

where $\mathcal{C}$ is the set of cliques, each potential $\psi_c = \exp(-E_c)$ is a non-negative function of the variables in clique $c$ , and the partition function $Z$ makes the distribution sum to 1 (Murphy, 2023, Eqs. 4.69–4.73). Low energy means high probability. Geman and Geman include a temperature, $\pi = e^{-U/T}/Z$ (1984, p. 731).

The Hammersley–Clifford theorem says that a strictly positive distribution with the conditional independences of $G$ can be written in this form over the (maximal) cliques of $G$ ; the converse, that every Gibbs distribution is Markov with respect to its graph, follows from the factorization (our remark, from the local conditional below). Murphy (2023, §4.3.1.1) notes that the theorem was never published and that it holds only for positive distributions, which rules out some models with hard constraints. Geman and Geman (1984, abstract and p. 722) call the result the Gibbs distribution–MRF equivalence and use it in both directions: the energy is the convenient way to specify a model, and the local conditionals are the convenient way to compute with it.

The local conditional of one site follows from the energy alone, because every clique not containing $s$ cancels between numerator and denominator:

p(x_s \mid x_{\mathcal{N}(s)}) = \frac{\exp\!\big(-\sum_{c \ni s} E_c(x_c)\big)}{\sum_{x_s'} \exp\!\big(-\sum_{c \ni s} E_c(x_s', x_{c \setminus s})\big)}.

$Z$ drops out, which is why sampling and greedy optimization can work site by site even though $Z$ itself is a sum over every configuration of the image.

Ising and Potts models

The Ising model has binary spins $x_s \in \{-1, +1\}$ and pairwise energy

E(x) = -J \sum_{s \sim t} x_s x_t - \alpha \sum_s x_s,

where $s \sim t$ runs over neighbouring pairs, $J$ is the coupling and $\alpha$ a unary bias, the “external field” (Murphy, 2023, Eqs. 4.78 and 4.82). With $J > 0$ neighbours prefer to agree. The Potts model generalizes it to $K$ labels with potential $e^{J}$ for equal neighbours and $e^0$ otherwise; for two labels it reduces to the Ising model with $J_{\text{Potts}} = 2 J_{\text{Ising}}$ (Murphy, 2023, §4.3.2.2). In the cost form used in vision, a penalty $\lambda [x_s \ne x_t]$ on unequal neighbours, the factor of 2 appears because $-J x_s x_t = 2J[x_s \ne x_t] - J$ for $\pm 1$ spins: a Potts penalty $\lambda$ is an Ising coupling $J = \lambda/2$ (our arithmetic). Szeliski et al. (2008, §2.1) likewise call the two-label Potts cost the Ising model.

The coupling, or equivalently the inverse temperature, controls what samples look like. On an infinite square lattice the Ising model has a phase transition at the critical coupling

J_c = \tfrac{1}{2} \log\big(1 + \sqrt{2}\big) \approx 0.44,

below which samples break into many small clusters and above which large clusters form; the $K$ -state Potts model’s transition is at $J_c = \log(1 + \sqrt{K})$ (Murphy, 2023, Eqs. 4.79 and 4.85, and Figs. 4.17–4.18, which show samples on either side). In our reading, samples from such priors are blobs, not images: an MRF prior on its own encodes local agreement, not shapes.

Gaussian Markov random fields

For continuous variables the pairwise model with quadratic energies is a multivariate Gaussian in information form, $p(x) \propto \exp(\eta^{\mathsf T} x - \tfrac{1}{2} x^{\mathsf T} \Lambda x)$ , with precision matrix $\Lambda = \Sigma^{-1}$ (Murphy, 2023, §4.3.5.1). A zero entry $\Lambda_{ij} = 0$ means there is no edge between $x_i$ and $x_j$ , and then $x_i$ and $x_j$ are independent given all the other variables. The sparsity pattern of the precision matrix is therefore the graph of the MRF, the reading that factor graphs use for the information matrix of a least-squares problem. A grid image with a quadratic smoothness prior and Gaussian noise is a Gaussian MRF whose MAP estimate solves a sparse linear system; Szeliski et al.’s (2008, §4.4 and §5) “House” denoising benchmark is of this kind; on it a hierarchically preconditioned conjugate-gradient method outperformed every other algorithm by a large margin in time, although after its real-valued solution was rounded to integer labels its energy was slightly higher than TRW-S’s.

The posterior of a hidden MRF

Images enter as observations of a hidden field. With a hidden labelling $x$ , an MRF prior $p(x)$ and an observation $y$ whose pixels are independent given $x$ ,

p(x \mid y) \propto p(x) \prod_s p(y_s \mid x_s),

and taking negative logarithms gives the energy

E(x) = \sum_s D_s(x_s) + \sum_{s \sim t} V_{st}(x_s, x_t), \qquad D_s(x_s) = -\log p(y_s \mid x_s),

up to a constant, the labelling energy that graph cuts minimize. Szeliski et al. (2008, Eqs. 1–3) write the energy of a pixel-labelling problem exactly this way and say that in the MRF framework the data energy comes from the negative log-likelihood of the measurement noise and the smoothness energy from the negative log of the prior. The posterior has the same cliques as the prior, plus single-site terms, so it is again an MRF: Geman and Geman (1984, abstract) state this for a range of degradations including blur, nonlinear deformations and additive or multiplicative noise, with a structure “akin to” the prior’s. Murphy (2023, §4.3.8, Fig. 4.26) gives the binary case, an Ising prior with Gaussian pixel noise, and notes that exact computation of its posterior marginals or MAP estimate is intractable for large grids. For the MAP estimate with two labels and $J > 0$ , however, a minimum cut finds the minimum exactly, as Graph Cuts reports from Greig, Porteous and Seheult (1989).

Pairwise terms used in vision

Szeliski et al. (2008, Eqs. 4–5) use smoothness costs of the form $w_{st} V(|x_s - x_t|)$ with the clipped monomial

V(\Delta x) = \min\big(|\Delta x|^k, V_{\max}\big), \qquad k \in \{1, 2\},

which gives the Potts model with $V_{\max} = 1$ , the truncated linear and truncated quadratic costs with larger $V_{\max}$ , and a convex quadratic cost without truncation. Truncation lets the field change abruptly at real discontinuities instead of smoothing across them. The weights $w_{st}$ can vary over the image; in their segmentation benchmarks the Potts weight falls with the colour difference of the two neighbours (their Eq. 10), so that label boundaries are cheaper along image edges. In their photomontage benchmarks the pairwise cost depends on the labels and the images themselves and measures how visible a seam would be.

Geman and Geman (1984, pp. 722, 732) handled discontinuities differently: they model the image as a pair $(F, L)$ of an intensity process $F$ on the pixels and an unobserved line process $L$ of edge elements on a dual lattice between them. When an edge element between two pixels is “on”, the bond between them is broken and their pairwise energy is set to zero, and the line process has its own energy on cliques of four edge sites (p. 732).

Properties

Locality makes sampling and local search cheap. Each site’s conditional depends only on its neighbours, so a site can be resampled or re-optimized in constant time. On the 4-connected grid the pixels of one colour of a checkerboard are independent given the other colour, so half the image can be updated in parallel (Murphy, 2023, §12.3.3); Geman and Geman (1984, p. 722) stressed this parallelism.
The partition function is intractable. $Z$ sums over all $L^n$ configurations; Geman and Geman (1984, p. 722) note that even a $64 \times 64$ binary image rules out direct search. Probabilities of whole configurations cannot be evaluated, only ratios, and learning the parameters is hard because the gradient of $\log Z$ is an expectation under the model (below).
Estimates depend on the loss. The MAP estimate, $\arg\max_x p(x \mid y)$ , is the minimum-energy labelling. The MPM estimate (maximizer of the posterior marginals) takes for each site the most probable value of its marginal $p(x_s \mid y)$ . Murphy (2023, §9.3.3.2) explains that they differ even when the marginals are exact: MPM minimizes the expected number of wrong sites (bit error), MAP the probability that the whole configuration is wrong, and MAP respects global constraints that MPM can break, while MPM can be more robust when such constraints are weak. A third estimate, the posterior mean, minimizes squared error; Geman and Geman (1984, p. 722) note that their sampler could compute it but do not pursue it.
Exact MAP is hard in general. With more than two labels, minimizing an energy whose smoothness costs include the Potts model is NP-hard, which Szeliski et al. (2008, §2.1) take to mean that most MRF energies used in vision are intractable. Special cases are exact: two labels with Ising costs by graph cuts, convex costs on ordered labels by other graph constructions, and tree-structured graphs by dynamic programming, as in pictorial structures.
The lowest energy is not the right answer. Szeliski et al. (2008, §6 and Table 1) found that on the Middlebury stereo benchmarks the ground-truth disparities had substantially higher energy than the labellings their optimizers found, and they report Meltzer et al.‘s finding that the global minimum was no more accurate than graph cuts’ or LBP’s results. For the models they considered, they conclude, better minimization is unlikely to give significantly more accurate labellings.

Inference algorithms

Iterated conditional modes (ICM) sets each site in turn to the value that minimizes the energy given its neighbours and repeats until nothing changes. Szeliski et al. (2008, §3.1), who attribute it to Besag (1986), note that it always converges, quickly, to a local minimum, that the result is very sensitive to the starting labelling, and that they start it “winner-take-all”, with each pixel’s lowest-cost label.
Gibbs sampling and simulated annealing. The Gibbs sampler, named by Geman and Geman (1984, p. 722), visits each site and draws its value from the local conditional at temperature $T$ ; their Theorem A (p. 731) shows that the distribution of the configuration then converges to the Gibbs distribution whatever the start. Lowering $T$ while sampling is annealing: Theorem B (p. 731) guarantees convergence to the set of global minima if $T(t) \to 0$ no faster than $N\Delta/\log t$ , where $N$ is the number of sites and $\Delta$ the range of the energy. The authors call this rate the main practical weakness: with $N = 20{,}000$ and $\Delta = 1$ it would take $e^{40{,}000}$ site visits to reach $T = 0.5$ . In practice they used $T(k) = C/\log(1 + k)$ over $k$ full sweeps, with $C = 3$ or 4, and obtained good results by lowering $T$ from about 4 to 0.5 over 300–1000 sweeps (p. 732). Murphy (2023, §12.9.1) reports the logarithmic guarantee and that exponential cooling, $T_{t+1} = \gamma T_t$ , is common in practice. The site update is the “heat bath” version of the Metropolis algorithm (Geman and Geman, 1984, p. 721). Graph cuts reports Greig, Porteous and Seheult’s (1989) finding that annealing with practicable schedules can give poor approximations to the exact binary MAP.
Graph cuts compute the exact minimum for binary submodular energies and strong local minima, with respect to expansion or swap moves, for many labels (graph cuts).
Loopy belief propagation (LBP) passes messages between neighbours as if the grid were a tree. Max-product (min-sum, in energies) messages approximate the MAP, and sum-product messages approximate the marginals (Szeliski et al., 2008, §3.3). On graphs with cycles it need not converge; Murphy (2023, §9.4.4.2) recommends damping, replacing each message by a weighted average of its new and old values, with a weight around 0.5 as standard practice. Semi-global matching has been interpreted as min-sum belief propagation, as that article reports.
Tree-reweighted message passing (TRW-S) is similar to LBP on the surface but also computes a lower bound on the energy that never decreases, which Szeliski et al. (2008, §3.4) used to measure how close every method came to the optimum.
Mean-field variational inference replaces the posterior by a fully factorized one; Murphy (2023, §4.3.8, Fig. 4.27) uses it to denoise a binary image under an Ising prior.

Szeliski et al. (2008) compared ICM, two LBP implementations, swap and expansion graph cuts and TRW-S on stereo, photomontage, binary segmentation and denoising energies. On every benchmark the best methods came within 1% of TRW-S’s lower bound, and often within 0.1%. On the binary segmentation benchmarks graph cuts and TRW-S reached the global minimum, and both LBP implementations came within 0.1% of it but never attained it (§5). Their strongest impression was how much better the modern methods were than ICM, which they took to reflect that the methods used until the late 1990s performed poorly, citing Boykov et al.’s similar comparison of graph cuts with simulated annealing. Expansion moves were the fastest of the competitive methods, TRW-S was consistently strong, and both LBP implementations did surprisingly poorly on many benchmarks (§6).

Learning

For an MRF with log-linear potentials, $p(x \mid \theta) \propto \exp\big(\sum_c \theta_c^{\mathsf T} \phi_c(x)\big)$ , the gradient of the average log-likelihood is the difference between the features’ average over the data and their expectation under the model, so the maximum-likelihood fit matches the two (moment matching); the model expectation, the derivative of $\log Z$ , needs samples from the model at every step (Murphy, 2023, §4.3.9.1–4.3.9.2), which is what makes maximum likelihood estimation expensive here. Pseudolikelihood, which Murphy credits to Besag (1975), maximizes instead the product of each site’s conditional given its neighbours, which needs only local normalization; it equals maximum likelihood for Gaussian MRFs, is consistent in the large-sample limit, and in experiments Murphy cites worked as well as exact maximum likelihood for fully observed Ising models (§4.3.9.3). Graph Cuts’ example shows how much a labelling depends on the smoothness weight $\lambda$ .

Examples

Same energy, different optimizers

The example builds one binary MRF energy on a synthetic image and minimizes it with five methods, so that the energies they reach can be compared with the exact minimum.

Scene: a $64 \times 64$ binary ground truth with a ring (centre at pixel (21, 21), $9.5 < r \le 14.5$ pixels) and ten $6 \times 6$ squares in two rows, intensity 0 for background and 1 for foreground, with independent Gaussian noise $\sigma = 0.8$ (noise seeds 0–4) and a noise-free copy.
Energy: $E(L) = \sum_p D_p(L_p) + \lambda \sum_{p \sim q} [L_p \ne L_q]$ on the 4-connected grid, with the Gaussian negative log-likelihood $D_p(l) = (I_p - \mu_l)^2/(2\sigma^2)$ (its constant dropped) and $\lambda = 0.6$ . The value of $\lambda$ was fixed from the construction before any code ran: a $6 \times 6$ square gains on average $36\,\Delta = 28.1$ in data cost, with $\Delta = (\mu_1 - \mu_0)^2/(2\sigma^2) = 0.78$ , against $24\lambda = 14.4$ for its boundary, and a lone wrong pixel costs $4\lambda = 2.4$ in disagreements. In spins the prior is an Ising model with coupling $J = \lambda/2 = 0.3$ , below $J_c \approx 0.44$ . All costs are multiplied by 1000 and rounded to integers once, because SciPy’s maximum flow needs integer capacities, and every method works on these same integers.
Exact minimum: a minimum $s$ – $t$ cut with scipy.sparse.csgraph.maximum_flow, built as in graph cuts.
ICM from the winner-take-all labelling (each pixel’s lower data cost), as Szeliski et al. start it. It updates the two colours of a checkerboard in turn, keeps a pixel’s label on ties, and stops after a sweep with no change. The visiting order is our choice; Besag’s paper was not read.
Annealing: a Gibbs sampler with checkerboard updates and Geman and Geman’s schedule $T(k) = 3/\log(1 + k)$ for $k = 1, \dots, 1000$ sweeps (1984, p. 732), from $T = 4.33$ down to 0.434, started from winner-take-all. Geman and Geman report the state after the last sweep; finishing with ICM, which is annealing at $T = 0$ , is a choice the sources leave open, so both are printed.
Loopy BP: min-sum messages, all updated at once and damped by 0.5, 200 iterations, each pixel taking its lower belief.
MPM: a Gibbs sampler at $T = 1$ , which samples the posterior of the model, with 200 sweeps discarded and 1000 recorded; each pixel takes the label it held in more than half of them.

Each method’s energy, its gap to the minimum (absolute, and relative to the minimum energy), and its pixel errors against the ground truth are printed. The relative gap depends on the constant dropped from the data term, so the absolute gap is printed beside it.

Six checks test the code and data: the scene’s pixel counts (C1); the maximum flow equals the energy of the labelling it returns (C2); the noise’s standard deviation is within 0.01 of 0.8 (C3); the winner-take-all labelling’s error fraction is within 0.015 of $\Phi(-(\mu_1 - \mu_0)/2\sigma) = 0.2660$ (C4); no method falls below the minimum (C5); and on noise-free data the exact labelling has no errors, ICM, annealing with the ICM finish and loopy BP reach the minimum exactly, and MPM has no errors (C6). The hypotheses were written before any script for the example ran: (H1) ICM’s relative gap is larger than loopy BP’s on every seed, after Szeliski et al.’s Discussion; (H2) loopy BP’s relative gap is below 0.001 on every seed, the 0.1% of Szeliski et al.’s binary segmentation benchmarks; (H3) annealing with the ICM finish has a smaller relative gap than ICM on every seed; (H4) the annealed final state, without the finish, lies above the minimum on every seed, since its last temperature is above zero; (H5) MPM has no more pixel errors than the exact MAP labelling on at least three of the five seeds, while its energy lies above the minimum on every seed. H5 rests on Murphy’s account of MPM as minimizing per-pixel error under the model; in our reading it was uncertain from the start, since the ring and squares are not drawn from the Ising prior the model assumes.

The script uses NumPy 2.3.5 and SciPy 1.16.3 and runs in about 2.5 seconds.

# Same energy, different optimizers: a binary MRF (Potts/Ising prior, Gaussian likelihood) on a 64 x 64 scene.
import os
os.environ["OPENBLAS_NUM_THREADS"] = "1"; os.environ["OMP_NUM_THREADS"] = "1"; os.environ["MKL_NUM_THREADS"] = "1"
import numpy as np
from scipy.sparse import csr_array
from scipy.sparse.csgraph import maximum_flow, breadth_first_order
from scipy.stats import norm

H = W = 64
SCALE = 1000                          # costs in integer thousandths (maximum_flow needs integers)
MU0, MU1, SIGMA, LAM = 0.0, 1.0, 0.8, 0.6
LQ = int(round(LAM * SCALE))
C_ANNEAL, K_ANNEAL = 3.0, 1000        # Geman & Geman's schedule T(k) = C / log(1 + k)
BURN, KEEP = 200, 1000                # MPM sampler at T = 1
LBP_ITERS, DAMP = 200, 0.5

yy, xx = np.mgrid[0:H, 0:W]
r2 = (yy - 21) ** 2 + (xx - 21) ** 2
ring = (r2 > 9.5 ** 2) & (r2 <= 14.5 ** 2)
squares = np.zeros((H, W), bool)
for r0 in (42, 52):
    for c0 in (6, 18, 30, 42, 54):
        squares[r0:r0 + 6, c0:c0 + 6] = True
gt = (ring | squares).astype(int)
black = (yy + xx) % 2 == 0            # checkerboard colours: each colour's pixels are independent given the other

def costs(img):
    D0 = np.rint((img - MU0) ** 2 / (2 * SIGMA ** 2) * SCALE).astype(np.int64)
    D1 = np.rint((img - MU1) ** 2 / (2 * SIGMA ** 2) * SCALE).astype(np.int64)
    return D0, D1

def energy(L, D0, D1):
    cuts = (L[:, 1:] != L[:, :-1]).sum() + (L[1:, :] != L[:-1, :]).sum()
    return int(np.where(L == 1, D1, D0).sum() + LQ * cuts)

def ones_around(L):                   # number of 4-neighbours labelled 1, and number of 4-neighbours
    P = np.pad(L, 1)
    n1 = P[:-2, 1:-1] + P[2:, 1:-1] + P[1:-1, :-2] + P[1:-1, 2:]
    Q = np.pad(np.ones_like(L), 1)
    return n1, Q[:-2, 1:-1] + Q[2:, 1:-1] + Q[1:-1, :-2] + Q[1:-1, 2:]

def local(L, D0, D1):                 # energy of each pixel's two labels given its neighbours (integers)
    n1, nn = ones_around(L)
    return D0 + LQ * n1, D1 + LQ * (nn - n1)

def maxflow(D0, D1):
    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(np.ravel(u)); cols.append(np.ravel(v)); caps.append(np.ravel(c))
    add(np.full(n, s), idx, D0); add(idx, np.full(n, t), D1)
    for a, b in [(idx[:, :-1], idx[:, 1:]), (idx[:-1, :], idx[1:, :])]:
        w = np.full(a.size, LQ); add(a, b, w); add(b, a, w)
    G = csr_array((np.concatenate(caps).astype(np.int32), (np.concatenate(rows), np.concatenate(cols))), shape=(n + 2, n + 2))
    res = maximum_flow(G, s, t)
    R = G - res.flow; R.data[R.data < 0] = 0; R.eliminate_zeros()
    L = np.zeros(n + 2, int); L[breadth_first_order(R, s, return_predecessors=False)] = 1
    return L[:n].reshape(H, W), int(res.flow_value)

def icm(L, D0, D1, max_sweeps=1000):
    L = L.copy()
    for sweep in range(1, max_sweeps + 1):
        changed = 0
        for colour in (black, ~black):
            e0, e1 = local(L, D0, D1)
            new = np.where(e1 < e0, 1, np.where(e0 < e1, 0, L))      # ties keep the current label
            m = colour & (new != L); changed += int(m.sum()); L[m] = new[m]
        if changed == 0:
            return L, sweep
    return L, max_sweeps

def gibbs_sweep(L, D0, D1, T, rng):
    for colour in (black, ~black):
        e0, e1 = local(L, D0, D1)
        p1 = 1.0 / (1.0 + np.exp((e1 - e0) / (SCALE * T)))
        L[colour] = (rng.random((H, W)) < p1)[colour]

def anneal(L, D0, D1, rng):
    L = L.copy()
    for k in range(1, K_ANNEAL + 1):
        gibbs_sweep(L, D0, D1, C_ANNEAL / np.log(1 + k), rng)
    return L

def mpm(L, D0, D1, rng):
    L = L.copy(); count = np.zeros((H, W))
    for _ in range(BURN):
        gibbs_sweep(L, D0, D1, 1.0, rng)
    for _ in range(KEEP):
        gibbs_sweep(L, D0, D1, 1.0, rng); count += L
    return (count / KEEP > 0.5).astype(int)

def lbp(D0, D1):
    D = np.stack([D0, D1], -1).astype(float)
    # m[d][i, j, l]: message arriving at pixel (i, j) from its neighbour in direction d (0 up, 1 down, 2 left, 3 right)
    m = np.zeros((4, H, W, 2))
    for _ in range(LBP_ITERS):
        h = D + m.sum(0)
        new = np.zeros_like(m)
        def send(a):                  # a[..., l_p] -> message over l_q: min(a[l_q], a[1 - l_q] + LQ), min-normalized
            out = np.minimum(a, a[..., ::-1] + LQ)
            return out - out.min(-1, keepdims=True)
        new[0, 1:] = send(h[:-1] - m[1, :-1])          # from the pixel above: it excludes what it received from below
        new[1, :-1] = send(h[1:] - m[0, 1:])
        new[2, :, 1:] = send(h[:, :-1] - m[3, :, :-1])
        new[3, :, :-1] = send(h[:, 1:] - m[2, :, 1:])
        new = DAMP * new + (1 - DAMP) * m
        change = np.abs(new - m).max(); m = new
    b = D + m.sum(0)
    return (b[..., 1] < b[..., 0]).astype(int), change

def run(img, seed):
    D0, D1 = costs(img)
    Lmin, flow = maxflow(D0, D1); Emin = energy(Lmin, D0, D1)
    wta = (D1 < D0).astype(int)
    L_icm, sweeps = icm(wta, D0, D1)
    L_ann = anneal(wta, D0, D1, np.random.default_rng(1000 + seed))
    L_annicm, _ = icm(L_ann, D0, D1)
    L_lbp, change = lbp(D0, D1)
    L_mpm = mpm(wta, D0, D1, np.random.default_rng(2000 + seed))
    rows = [("max-flow (exact)", Lmin), ("ICM from WTA", L_icm), ("annealing, final state", L_ann),
            ("annealing + ICM", L_annicm), ("loopy BP (min-sum)", L_lbp), ("MPM, Gibbs T=1", L_mpm)]
    out = {}
    for name, L in rows:
        E = energy(L, D0, D1)
        out[name] = (E, E - Emin, (E - Emin) / Emin, int((L != gt).sum()))
    return out, flow, Emin, energy(gt, D0, D1), sweeps, change

def show(out):
    print(f"  {'method':<24}{'energy':>12}{'gap':>11}{'rel. gap':>11}{'errors':>8}")
    for name, (E, g, r, e) in out.items():
        print(f"  {name:<24}{E / SCALE:12.3f}{g / SCALE:11.3f}{r:11.6f}{e:8d}")

print(f"scene: ring {int(ring.sum())} px, squares {int(squares.sum())} px, overlap {int((ring & squares).sum())}, "
      f"foreground {int(gt.sum())} of {H * W}")
clean = np.where(gt == 1, MU1, MU0).astype(float)
out, flow, Emin, Egt, sweeps, change = run(clean, 99)
print(f"noise-free data: max flow {flow / SCALE:.3f}, minimum energy {Emin / SCALE:.3f}, ICM sweeps {sweeps}, "
      f"LBP last change {change:.3e}")
show(out)
expected_err = norm.cdf(-(MU1 - MU0) / (2 * SIGMA))
res = {}
for seed in range(5):
    noise = np.random.default_rng(seed).normal(0, SIGMA, (H, W))
    img = clean + noise
    D0, D1 = costs(img)
    wta_err = int(((D1 < D0).astype(int) != gt).sum())
    out, flow, Emin, Egt, sweeps, change = run(img, seed)
    res[seed] = out
    print(f"seed {seed}: noise std {noise.std():.4f}; WTA error fraction {wta_err / (H * W):.4f} (expected {expected_err:.4f}); "
          f"max flow {flow / SCALE:.3f}; energy of ground truth {Egt / SCALE:.3f}; ICM sweeps {sweeps}; "
          f"LBP last change {change:.3e}, converged {change < 1e-3}")
    show(out)
    below = [n for n, v in out.items() if v[1] < 0]
    print(f"  methods below the max-flow energy: {below if below else 'none'}")
print("range over seeds 0-4:")
for name in res[0]:
    r = [res[s][name][2] for s in res]; e = [res[s][name][3] for s in res]
    print(f"  {name:<24} rel. gap {min(r):.6f} to {max(r):.6f}, errors {min(e)} to {max(e)}")

Output:

scene: ring 372 px, squares 360 px, overlap 0, foreground 732 of 4096
noise-free data: max flow 259.200, minimum energy 259.200, ICM sweeps 1, LBP last change 0.000e+00
  method                        energy        gap   rel. gap  errors
  max-flow (exact)             259.200      0.000   0.000000       0
  ICM from WTA                 259.200      0.000   0.000000       0
  annealing, final state       301.354     42.154   0.162631      34
  annealing + ICM              259.200      0.000   0.000000       0
  loopy BP (min-sum)           259.200      0.000   0.000000       0
  MPM, Gibbs T=1               259.200      0.000   0.000000       0
seed 0: noise std 0.7980; WTA error fraction 0.2712 (expected 0.2660); max flow 2203.680; energy of ground truth 2297.721; ICM sweeps 5; LBP last change 5.870e-51, converged True
  method                        energy        gap   rel. gap  errors
  max-flow (exact)            2203.680      0.000   0.000000     202
  ICM from WTA                2288.702     85.022   0.038582     425
  annealing, final state      2299.261     95.581   0.043373     314
  annealing + ICM             2222.907     19.227   0.008725     256
  loopy BP (min-sum)          2203.927      0.247   0.000112     191
  MPM, Gibbs T=1              2294.134     90.454   0.041047     352
  methods below the max-flow energy: none
seed 1: noise std 0.8023; WTA error fraction 0.2651 (expected 0.2660); max flow 2230.836; energy of ground truth 2319.260; ICM sweeps 5; LBP last change 2.018e-47, converged True
  method                        energy        gap   rel. gap  errors
  max-flow (exact)            2230.836      0.000   0.000000     226
  ICM from WTA                2306.610     75.774   0.033967     435
  annealing, final state      2343.688    112.852   0.050587     286
  annealing + ICM             2244.150     13.314   0.005968     228
  loopy BP (min-sum)          2231.214      0.378   0.000169     234
  MPM, Gibbs T=1              2323.307     92.471   0.041451     359
  methods below the max-flow energy: none
seed 2: noise std 0.8015; WTA error fraction 0.2727 (expected 0.2660); max flow 2210.714; energy of ground truth 2314.981; ICM sweeps 4; LBP last change 2.167e-50, converged True
  method                        energy        gap   rel. gap  errors
  max-flow (exact)            2210.714      0.000   0.000000     222
  ICM from WTA                2297.124     86.410   0.039087     455
  annealing, final state      2322.550    111.836   0.050588     332
  annealing + ICM             2223.478     12.764   0.005774     264
  loopy BP (min-sum)          2210.714      0.000   0.000000     222
  MPM, Gibbs T=1              2279.741     69.027   0.031224     355
  methods below the max-flow energy: none
seed 3: noise std 0.8017; WTA error fraction 0.2739 (expected 0.2660); max flow 2224.223; energy of ground truth 2316.155; ICM sweeps 6; LBP last change 6.980e-43, converged True
  method                        energy        gap   rel. gap  errors
  max-flow (exact)            2224.223      0.000   0.000000     230
  ICM from WTA                2312.870     88.647   0.039855     429
  annealing, final state      2337.943    113.720   0.051128     325
  annealing + ICM             2235.226     11.003   0.004947     237
  loopy BP (min-sum)          2224.223      0.000   0.000000     230
  MPM, Gibbs T=1              2300.589     76.366   0.034334     355
  methods below the max-flow energy: none
seed 4: noise std 0.7967; WTA error fraction 0.2754 (expected 0.2660); max flow 2200.864; energy of ground truth 2290.543; ICM sweeps 5; LBP last change 1.148e-44, converged True
  method                        energy        gap   rel. gap  errors
  max-flow (exact)            2200.864      0.000   0.000000     192
  ICM from WTA                2309.602    108.738   0.049407     450
  annealing, final state      2302.121    101.257   0.046008     285
  annealing + ICM             2216.818     15.954   0.007249     240
  loopy BP (min-sum)          2200.864      0.000   0.000000     192
  MPM, Gibbs T=1              2293.953     93.089   0.042297     354
  methods below the max-flow energy: none
range over seeds 0-4:
  max-flow (exact)         rel. gap 0.000000 to 0.000000, errors 192 to 230
  ICM from WTA             rel. gap 0.033967 to 0.049407, errors 425 to 455
  annealing, final state   rel. gap 0.043373 to 0.051128, errors 285 to 332
  annealing + ICM          rel. gap 0.004947 to 0.008725, errors 228 to 264
  loopy BP (min-sum)       rel. gap 0.000000 to 0.000169, errors 191 to 234
  MPM, Gibbs T=1           rel. gap 0.031224 to 0.042297, errors 352 to 359

All six checks held. The scene has the planned 372 ring and 360 square pixels without overlap (C1); the maximum flow equals the minimum energy in all six runs (C2); the noise’s standard deviation is 0.7980, 0.8023, 0.8015, 0.8017 and 0.7967 on seeds 0–4 (C3); the winner-take-all labelling misclassifies a fraction 0.2712, 0.2651, 0.2727, 0.2739 and 0.2754 of the pixels against the predicted 0.2660 (C4); no method’s energy is below the minimum on any seed, and on noise-free data every gap is zero or positive (C5); and on noise-free data the exact labelling has no errors and ICM, annealing with the ICM finish, loopy BP and MPM all reach it (C6). The annealed final state alone misses it on clean data, by 34 pixels: at its last temperature, 0.434, the sampler still draws higher-energy labels with positive probability.

ICM. H1 held: ICM’s relative gap is 0.033967–0.049407 on the five seeds, against at most 0.000169 for loopy BP. ICM stops after 5, 5, 4, 6 and 5 sweeps with 425–455 wrong pixels, roughly twice the exact labelling’s count on every seed (425 against 202 on seed 0, 429 against 230 on seed 3). It starts from the winner-take-all labelling, with about 27% of the pixels wrong (C4), and only makes changes that lower the energy one pixel at a time; in our reading it removes isolated errors but keeps clusters of wrong pixels in which each pixel, flipped alone, would raise the energy.

Annealing. H4 held: the state after the last sweep is 95.581, 112.852, 111.836, 113.720 and 101.257 energy units above the minimum on seeds 0–4. H3 held: finished with ICM, annealing lands 19.227, 13.314, 12.764, 11.003 and 15.954 above the minimum, a relative gap of 0.004947–0.008725 against ICM’s 0.033967–0.049407, smaller on every seed and still not the minimum. The schedule is the one Geman and Geman used in practice, not the rate their Theorem B requires for the guarantee.

Loopy BP. H2 held: the largest relative gap is 0.000169 (seed 1), below the 0.001 bar. Loopy BP converged on every seed, reached the exact minimum’s energy on seeds 2, 3 and 4, and was 0.247 and 0.378 above it on seeds 0 and 1. Szeliski et al. found that their two LBP implementations never attained the global minimum on their three binary segmentation benchmarks, which used Gaussian-mixture data costs and a contrast-sensitive Potts term on three colour images with user strokes; ours reached it on three of five synthetic instances.

MPM against MAP. H5 failed. Its error part failed on all five seeds: MPM has 352, 359, 355, 355 and 354 wrong pixels against the exact MAP labelling’s 202, 226, 222, 230 and 192, so MPM had no more errors than MAP on 0 of the 5 seeds, against the planned at least 3. Its energy part held: MPM lies 90.454, 92.471, 69.027, 76.366 and 93.089 above the minimum on seeds 0–4.

The ground truth is not the minimum. On every seed the true labelling has a higher energy than the minimum (2297.721 against 2203.680 on seed 0, 2290.543 against 2200.864 on seed 4), and the minimum has 192–230 wrong pixels. As with Szeliski et al.’s stereo ground truth, the right answer is not the most probable one under the model, so the exact optimizer’s errors belong to the model (the data term and $\lambda$ ), not to the optimizer.

Process note. The plan for this example listed nine departures from the batch plan’s description, before the pilot: (1) the “grid of squares” is two rows of five $6 \times 6$ squares, and the ring’s radii are not integers so that the noise-free ground truth has no one-pixel tips and is itself the minimum; (2) the absolute gap is printed beside the relative gap, because the relative gap depends on the data term’s dropped constant; (3) annealing is reported both as its final state and with an ICM finish, a choice the sources leave open; (4) the ground truth’s energy is printed for each seed, with no criterion; (5) checks C1–C6 were added; (6) ICM visits pixels in checkerboard order rather than a raster order; (7) the batch plan’s expectation that ICM ends with the largest gap is tested against loopy BP (H1) and against annealing with the ICM finish (H3), not against the raw annealed state, which is a sample at $T = 0.434$ , or against MPM, which does not minimize the energy; (8) the optional block of Gibbs samples from the Ising prior at two couplings was not run, and the prose describes Murphy’s figures instead; (9) the loopy BP and MPM settings were fixed in the plan before the pilot. The tolerances of two checks were set after a planning script had printed the values they test: C3 (0.01) and C4 (0.015), after the script had printed the same noise standard deviations and winner-take-all errors that the output above shows. C1’s ring count of 372 and C6’s zero errors for the exact labelling on noise-free data were likewise written after that script printed them; both are exact properties of the construction, not tolerances. That script ran only the exact max-flow minimum and printed none of the measures H1–H5 compare (no ICM, annealing, loopy BP or MPM result), so H1–H5 are tests, not illustrations. The pilot printed the same output as the first run; no change was made after it.

Common Misconceptions

“Markov means each pixel depends only on its neighbours.” Each pixel is conditionally independent of the rest given its neighbours. Its marginal correlations reach much further, through chains of neighbours; Geman and Geman (1984, p. 722) point out that with a hidden line process the marginal distribution of the intensities alone has a completely connected graph.
“The MRF energy can be turned into probabilities directly.” Only differences of energy can: $p(x) = e^{-E(x)}/Z$ needs $Z$ , a sum over all configurations. Local conditionals do not need it.
“Graph cuts solve MRFs.” They solve the binary submodular case exactly and approximate many-label metric cases (Graph Cuts); NP-hardness covers the rest (Szeliski et al., 2008).
“Simulated annealing finds the global minimum.” It is guaranteed to only with Geman and Geman’s logarithmic rate, which they call its main practical weakness; with the schedule they used in practice there is no such guarantee, and in the example both the annealed final state and its ICM finish stayed above the exact minimum on every seed.
“The MAP labelling is the best labelling.” MAP minimizes one loss; MPM and the posterior mean minimize others (Murphy, 2023, §9.3.3.2). And under a simplified model the ground truth can have a higher energy than the MAP labelling, as Szeliski et al. found for stereo and the example shows.
“Ising coupling and Potts penalty are the same number.” For two labels a Potts penalty $\lambda$ on unequal neighbours is an Ising coupling $\lambda/2$ (Murphy, 2023, §4.3.2.2).

Where It Is Used

Segmentation and labelling. Binary and multi-label image segmentation with Potts priors, solved by graph cuts, in interactive segmentation and GrabCut; normalized cuts instead partition a pixel-similarity graph without terminals, as Graph Cuts describes. A superpixel MRF that refines foreground masks ranked second overall at CDnet 2012, as background subtraction reports.
Stereo. Global stereo matching minimizes MRF energies over disparities; in the Scharstein–Szeliski evaluation global methods on 2D Markov random fields generally performed best, and semi-global matching approximates such energies along scanlines.
Restoration. Geman and Geman’s (1984) restoration of blurred and noisy images, binary denoising as in the example, and denoising and inpainting with truncated quadratic costs (Szeliski et al., 2008); image restoration places MAP estimation among its approaches.
Stitching. Photomontage seams chosen by minimizing a labelling energy whose pairwise term measures seam visibility (Szeliski et al., 2008), as in image stitching.
Part-based models. Pictorial structures are tree-structured MRFs over part locations, minimized exactly by dynamic programming.
Least-squares estimation. In SLAM and bundle adjustment the sparsity of the information matrix is the graph of a Gaussian MRF (factor graph).
Conditional models. Conditional random fields condition the potentials on the image; DeepLab refines network outputs with a fully connected one.