Markov random fields

Why you need this

A Markov chain generates a sequence: the next state depends on the current one. A side-scroller column has a previous column, so that fits. A tile on a top-down map has neighbors to the left, right, above, and below, and none of them is "previous". A Markov random field (MRF) is the model you get when you keep the local-dependence idea but drop the ordering. It is the general theory behind biome smoothing, texture-like fills, and coupled regions, and it is also the direct answer to the question "why not just run an HMM over the whole grid?"

The idea

Take a grid of $W \times H$ cells, width then height, and treat it as a graph (see Graphs and grids): each cell is a node and each cell is joined to its von Neumann neighbors (the 4 cells sharing an edge). Each cell $c$ holds a state $X_c$ from $N$ possible states, numbered $1, 2, \dots, N$.

The local Markov property on a graph. A cell is independent of every other cell once you know its neighbors:

$$P(X_c = s \mid \text{all other cells}) = P(X_c = s \mid \text{neighbors of } c).$$

On a chain, "neighbors" means the one cell before and the one after, and this reduces to the ordinary Markov property. On a grid, it means four cells, and there is no direction.

Why there is no transition matrix. A chain has a natural order, so you can define the whole joint distribution by a start distribution $\pi$ and a transition matrix $A$ with $A_{ij} = P(X_{t+1} = j \mid X_t = i)$, and sample it by walking left to right. A grid has no start and no order. Instead, the joint distribution is defined by potentials: a compatibility score $\phi(a, b) \geq 0$ for every neighboring pair of states, and optionally a per-cell score $\psi_c(a)$ that says how much cell $c$ likes state $a$ on its own. The probability of a whole grid configuration $x$ is proportional to the product of all those scores:

$$P(x) = \frac{1}{Z} \prod_{\text{cells } c} \psi_c(x_c) \prod_{\text{neighbor pairs } (c, d)} \phi(x_c, x_d).$$

The $Z$ is the sum of that product over every possible grid, which is what makes the probabilities add to 1. For a 16x16 grid with 2 states that sum has $2^{256}$ terms. Nobody computes $Z$. So you never compute $P(x)$ either; you sample grids from the distribution without knowing the normalizer.

Gibbs sampling. The trick is that the conditional probability of one cell given its neighbors does not need $Z$. The scores that involve cell $c$ are its own $\psi_c$ and the four $\phi$ terms with its neighbors; everything else cancels. So:

$$P(X_c = s \mid \text{neighbors}) = \frac{\psi_c(s) \prod_{d \in \text{nbrs}(c)} \phi(s, x_d)}{\sum_{s'=1}^{N} \psi_c(s') \prod_{d \in \text{nbrs}(c)} \phi(s', x_d)}.$$

The denominator is a sum over $N$ states, not $2^{256}$ grids. Gibbs sampling is then:

1. Fill the grid with random states from the seeded generator.
2. Visit every cell in some fixed order (row by row is fine). For each cell, compute the score of every state against the cell's current neighbors, divide by their sum to get a categorical distribution, draw $u \sim \mathrm{Uniform}(0, 1)$ from the seeded generator, and set the cell to the state the draw lands on.
3. One full pass is a sweep. Repeat for $S$ sweeps. Early sweeps still remember the random start; later sweeps are draws from the distribution the potentials define. Ten to a hundred sweeps is typical for a map.

Nothing is ever forbidden outright unless a potential is exactly 0. A cell can flip to a state its neighbors dislike, just rarely. That is the difference between a soft preference and a hard constraint.

Worked example

Two states, land (L) and water (W), on a $3 \times 3$ grid with von Neumann neighbors. No per-cell score, so $\psi_c(a) = 1$ for every cell and state. One pair potential that rewards matching neighbors:

$\phi(a, b)$	$b$ = L	$b$ = W
$a$ = L	2	1
$a$ = W	1	2

This is the Ising model from statistical physics, written with scores instead of energies: a physicist would write $\phi(a, b) = e^{J}$ when the pair matches and $e^{0} = 1$ when it does not, and our table is that with $J = \log 2 \approx 0.69$. Larger $J$ means blobbier maps.

The grid after the random fill, with cells named $(x, y)$, $x$ to the right and $y$ down, both from 0:

	x = 0	x = 1	x = 2
y = 0	L	L	W
y = 1	W	L	L
y = 2	L	W	W

Updating the center cell $(1, 1)$. Its four neighbors are $(1, 0)$ = L, $(0, 1)$ = W, $(2, 1)$ = L, $(1, 2)$ = W. Score each candidate state:

Candidate	Against L	Against W	Against L	Against W	Product
L	2	1	2	1	4
W	1	2	1	2	4

Two matching neighbors each way, so $P(\text{L}) = 4 / 8 = 0.5$ and $P(\text{W}) = 0.5$. Draw $u = 0.63$; L covers $[0, 0.5)$ and W covers $[0.5, 1)$, so the cell becomes W.

Updating the corner $(2, 0)$, currently W, with only two neighbors: $(1, 0)$ = L and $(2, 1)$ = L.

Candidate	Against L	Against L	Product
L	2	2	4
W	1	1	1

$P(\text{L}) = 4 / 5 = 0.8$, $P(\text{W}) = 0.2$. Draw $u = 0.41$, so the corner flips to L. The isolated water pixel in the corner has been smoothed away, which is exactly what this potential is for. After a full sweep and then several more, the grid settles into a few large patches of L and W with ragged edges, and the seed decides where they fall.

Adding a per-cell score. If a noise field says cell $(0, 2)$ has height 0.2 and you set $\psi(\text{W}) = 3$, $\psi(\text{L}) = 1$ for low cells, the products for that cell are multiplied by 3 for W and 1 for L before normalizing. The field now follows the height map where the height is confident and smooths where it is not.

How this relates to the neighboring pages

• Cellular automata are the deterministic cousin. Same grid, same neighborhood, but the update is a rule ("become wall if 5 or more of 8 neighbors are wall") instead of a draw from a distribution. If you take the Gibbs update and always pick the highest-scoring state, you get a cellular automaton.
• Wave Function Collapse is the hard-constraint cousin. Its adjacency rules are a pair potential that is 0 or 1, and instead of sampling in a way that could violate a rule, it never places a tile whose potential against a neighbor is 0. WFC gives you legal tilings; an MRF gives you textures.
• Hidden Markov models are the one-dimensional special case where the graph is a chain. And here is the direct answer to "why not run an HMM over a 2D map": the forward algorithm and Viterbi work because on a chain, everything to the left of column $t$ reaches everything to the right only through column $t$. That lets you summarize the whole left side in one vector of $N$ numbers per column and fill a table one column at a time, which is dynamic programming. On a grid, a row is not a bottleneck: cells in row 3 are joined to row 4 at every column, so the state of a whole row is one variable with $N^W$ values. For 6 states and a 16-wide chunk that is $6^{16} \approx 2.8 \times 10^{12}$ per row. Exact inference is gone, and sampling is what is left.

In a map generator

An MRF sits in the fill stage and in the regime stage.

• Texture-like biome fill. Inside a chunk labelled "swamp", the states are the swamp tileset's ground variants (mud, reeds, shallow water) with a potential that likes matching neighbors. Gibbs sweeps produce patches rather than salt-and-pepper noise, and the per-cell score can be a noise field.
• Smoothing a noisy biome map. The layout stage assigned a biome per chunk with a Markov chain walk, and the result has lone forest chunks inside a desert. Run a few Gibbs sweeps over the chunk grid with a matching potential and a per-cell score that is the original assignment. Isolated chunks get absorbed; large regions survive.
• Coupling regime chains over time. Each region's regime (calm, contested, depleted) is a regime chain that advances per tick. Add a pair potential between neighboring regions and update with a Gibbs step, and contested regions spread to their neighbors instead of flickering independently.

The data is the same as everywhere else in the pipeline: an integer state per cell plus a small $N \times N$ table of scores, all driven by one seeded generator.

Common mistakes

• Reading the grid from the wrong buffer. Gibbs updates a cell using the current neighbors, including ones already updated this sweep. Copying into a second buffer and reading from the old one is a different algorithm (synchronous update) that oscillates in checkerboards with a matching potential.
• Too few sweeps. The output still looks like the random fill with a little clumping. Symptom: uniform speckle. Double the sweep count until the picture stops changing.
• Too strong a potential. With $\phi(\text{same}) = 20$ the whole grid becomes one state after a few sweeps. Symptom: a chunk that is all water. Keep the ratio near 2 to 4, or add per-cell scores.
• Using it for hard rules. A cliff edge tile next to open grass is wrong, not merely unlikely. An MRF will produce it sometimes. Use WFC for anything with a rulebook.
• Unseeded draws. The draw at every cell must come from the seeded generator in a fixed visit order, or the same seed gives a different map.
• Expecting a transition matrix to exist. There is no row of numbers that says "given water on the left, the next cell is sand with probability 0.3". The potentials are relative scores, not probabilities, and only the conditional in the Gibbs step is a probability.

Cost

With $S$ sweeps over a $W \times H$ grid, $N$ states per cell, and a neighborhood of 4, each cell update scores $N$ states against 4 neighbors, so the total is $O(S \cdot W \cdot H \cdot N)$ time. Memory is one integer per cell plus the $N \times N$ potential table, $O(W H + N^2)$. It hurts when $S$ has to be large: a strongly coupled field on a big grid needs hundreds of sweeps to lose the memory of its random start, and a 256x256 grid with 8 states at 200 sweeps is about $10^8$ score evaluations. Run it per chunk, or run it on the coarse chunk grid rather than on tiles.

Going further

• The Ising model and the Potts model (the same idea with more than two states), read as an intuition pump for how the coupling strength controls patch size.
• Simulated annealing: start with a weak potential and strengthen it sweep by sweep, which reaches large patches faster.
• Belief propagation, the approximate inference method that generalizes the forward-backward algorithm to graphs with loops.
• Conditional random fields, where the per-cell scores come from observed data, as in labelling a photo's pixels.
• Texture synthesis by example, which is the MRF idea with potentials learned from a sample image.