Introduction

Image segmentation partitions an image into regions with homogeneous appearance and is fundamental to applications in medical imaging , satellite remote sensing , and materials science. The hidden Potts model provides a principled Bayesian framework for this task, jointly inferring pixel labels, Gaussian mixture parameters, and a spatial smoothing parameter \(\beta\) from the observed image data.

The package is available from CRAN at . It provides reference implementations of the principal algorithms for Bayesian inference in the hidden Potts model, with particular emphasis on handling the intractable normalising constant in the posterior for \(\beta\). A theoretical overview of these algorithms is given in ; the scalable PFAB algorithm is developed and compared with alternatives in ; and the external field prior was introduced in . The present paper focuses on the software: how the package is structured, how to call its functions, and how to interpret the output.

Related packages include and , which provide simulation from the Potts model without data, and , which uses a hidden Markov random field for MRI tissue classification. extends these approaches by providing multiple algorithms for estimating \(\beta\), supporting an external field prior , and handling 2D and 3D irregular masks.

The remainder of this paper is organised as follows. Section \(\ref{s:model}\) describes the hidden Potts model and the priors implemented in the package. Section \(\ref{s:lattice}\) explains how to construct the required lattice structures. Section \(\ref{s:algorithms}\) describes each algorithm for estimating \(\beta\). Section \(\ref{s:results}\) illustrates the package with synthetic data. Section \(\ref{s:summary}\) concludes.

The Hidden Potts Model

Model for the Pixel Labels

The Potts model is a type of Markov random field. It specifies a probability distribution over a vector of \(n\) discrete labels \(\mathbf{z} = (z_1, \ldots, z_n)\), each taking a value in \(\{1, \ldots, k\}\): \[ p(\mathbf{z} \mid \beta) \;=\; \frac{1}{C(\beta)} \exp\bigl\{ \beta \, S(\mathbf{z}) \bigr\}. \] The sufficient statistic \[ S(\mathbf{z}) = \sum_{\{i,j\} \in \mathcal{E}} \delta(z_i, z_j) \] counts matching label pairs across the edges \(\mathcal{E}\) of the lattice graph, where \(\delta(\cdot,\cdot)\) is the Kronecker delta function. Larger values of \(\beta \ge 0\) favour smooth, piecewise-constant labellings; \(\beta = 0\) gives independent pixels. The normalising constant \(C(\beta) = \sum_{\mathbf{z}} \exp\{\beta\, S(\mathbf{z})\}\) sums over all \(k^n\) label configurations and is intractable for all but the smallest lattices . Full derivations and motivation are given in and .

Observation Model and Priors

Observed pixel intensities \(y_i\) are modelled as a Gaussian mixture conditioned on the latent labels: \[ y_i \mid z_i = j \;\sim\; \mathcal{N}(\mu_j,\, \sigma_j^2), \qquad i = 1,\ldots,n, \quad j = 1,\ldots,k. \] Conjugate priors are placed on the mixture parameters: \(\mu_j \sim \mathcal{N}(\mu_0, \sigma_{\mu}^2)\) and \(\sigma_j^{-2} \sim \mathrm{Gamma}(\nu/2,\, s/2)\). The argument accepted by and is a named list with the following elements:

: number of labels \(k\)
, : prior mean and standard deviation for each Gaussian component mean \(\mu_j\)
, : rate and shape for the inverse-Gamma prior on each component variance \(\sigma_j^2\)
: two-element vector giving the bounds of a uniform prior on \(\beta\)

External Field Prior

The package also supports an external field prior , in which pixel-specific log-weights modulate label assignment using auxiliary information, such as a prior scan or calibration data.

Lattice Construction

Before calling any MCMC or SMC function, three lattice structures must be pre-computed from a binary mask matrix. This pre-computation is done once and the structures are reused across algorithm calls.

library(bayesImageS)

# Define a 100x100 pixel image domain
mask <- matrix(1, nrow = 100, ncol = 100)

# Four-connected neighbourhood (N/S/E/W)
neigh <- getNeighbors(mask, c(2, 2, 0, 0))

# Chequerboard partition for parallel Gibbs updates
block <- getBlocks(mask, 2)

# Undirected edge list for the sufficient statistic
edges <- getEdges(mask, c(2, 2, 0, 0))
maxS  <- nrow(edges)   # maximum possible value of S(z)

getNeighbors() returns an \(n \times m\) matrix where row \(i\) lists the indices of the neighbors of pixel \(i\); missing neighbors (boundary pixels) are coded as \(n+1\). The argument specifies connectivity: gives four-connectivity; gives eight-connectivity. Three-dimensional arrays are also accepted, making the package suitable for volumetric CT data.

getBlocks() partitions pixels into two independent sets using a chequerboard pattern, such that all neighbors of any pixel in block~1 belong to block~2 and vice versa. This enables parallel conditional updates in the Gibbs sampler for the labels .

getEdges() enumerates the unique undirected edges of the lattice. The number of edges equals the maximum value of \(S(\mathbf{z})\), which is useful for scaling convergence diagnostics.

Algorithms for the Inverse Temperature

The central computational challenge is Bayesian inference for \(\beta\), since the likelihood \(p(\mathbf{z} \mid \beta)\) involves the intractable normalising constant \(C(\beta)\). provides five algorithms for this problem. See for a unified theoretical treatment.

Pseudolikelihood

proposed replacing the joint Gibbs distribution with a product of full conditional distributions, each of which depends only on a pixel’s immediate neighbors and can be evaluated in closed form. The pseudolikelihood approximation is fast but can be biased, particularly near the phase transition of the Potts model . In , pseudolikelihood is available via .

Exchange Algorithm

The exchange algorithm of , generalised to doubly-intractable distributions by , introduces an auxiliary draw from the Potts model at the proposed \(\beta\) value such that the intractable normalising constants cancel in the Metropolis–Hastings ratio. Each iteration requires one complete simulation of the Potts model, making the exchange algorithm exact at the cost of being computationally intensive for large lattices. This is the default method in :

q      <- 2
priors <- list(k        = q,
               mu       = c(-1, 1),
               mu.sd    = rep(0.5, q),
               sigma    = rep(2, q),
               sigma.nu = rep(1.5, q),
               beta     = c(0, 2))

mh  <- list(algorithm = "ex", bandwidth = 0.5,
            adaptive = 500, auxiliary = 100)
res <- mcmcPotts(y, neigh, block, priors, mh,
                 niter = 10000, nburn = 2000)

The element specifies the number of Gibbs sweeps used to generate the auxiliary lattice. The element enables the Robbins–Monro adaptive scaling of the random-walk proposal .

Thermodynamic Integration

Thermodynamic integration , also called path sampling, estimates \(\log C(\beta)\) by numerically integrating \(\mathrm{E}[S(\mathbf{z}) \mid \beta]\) over a grid of \(\beta\) values. The resulting log-likelihood surface is used within a Metropolis–Hastings step for \(\beta\). Select this method with .

ABC-SMC

Approximate Bayesian computation avoids evaluating the likelihood by simulating pseudo-data and accepting parameter values whose simulated sufficient statistic falls within a tolerance of the observed value. The ABC-SMC implementation in propagates a population of \(N\) particles through a sequence of shrinking tolerances, with the schedule adapted automatically :

param <- list(npart = 10000, nstat = 50)
smc   <- smcPotts(y, neigh, block, param, priors)

A pre-processing step based on is available via initSedki(), which initialises the particle population efficiently using summary statistics derived from the image histogram.

PFAB

The parametric functional approximate Bayesian (PFAB) algorithm fits a parametric curve to the relationship \(\beta \mapsto \mathrm{E}[S(\mathbf{z}) \mid \beta]\) using a small number of pilot simulations. This surrogate function replaces the costly Potts simulations during the main SMC run, reducing the total number of model evaluations by orders of magnitude at negligible loss of accuracy. PFAB is the recommended method for images with more than \(10^5\) pixels, and is accessed via .

Examples

We demonstrate the package using synthetic data generated from a \(500 \times 500\) Ising lattice (\(q = 2\)) at four values of the inverse temperature: \(\beta \in \{0.22, 0.44, 0.88, 1.32\}\). These values span a range from near-independence (\(\beta = 0.22\)) through to strong spatial dependence above the critical temperature (\(\beta_c \approx 0.88\) for the Ising model). The data are included in the package as the dataset:

library(bayesImageS)
data("synth", package = "bayesImageS")

Each element of corresponds to one value of \(\beta\) and contains: true labels , noisy Gaussian observations , mixing parameters and , and the true sufficient statistic . Lattice structures for the \(500 \times 500\) grid are constructed as in Section~\(\ref{s:lattice}\).

We fit the model using the exchange algorithm with \(\beta\) fixed at the true value for the first dataset, to assess convergence of the label sampler:

mask  <- matrix(1, nrow = 500, ncol = 500)
neigh <- getNeighbors(mask, c(2, 2, 0, 0))
block <- getBlocks(mask, 2)

q      <- 2
priors <- list(k        = q,
               mu       = c(-1, 1),
               mu.sd    = rep(0.5, q),
               sigma    = rep(2, q),
               sigma.nu = rep(1.5, q),
               beta     = rep(synth[[1]]$beta, 2))

mh  <- list(algorithm = "ex", bandwidth = 1,
            adaptive = NA, auxiliary = 1)
res <- mcmcPotts(synth[[1]]$y, neigh, block,
                 priors, mh, niter = 100, nburn = 50)

data("res", package = "bayesImageS")
print(res$tm)

##    user  system elapsed 
##  23.299   3.729  13.109

mean(res$sum[51:100])

## [1] 277581.5

print(synth[[1]]$sum)

## [1] 277664

ts.plot(res$sum,
        xlab = "MCMC iteration",
        ylab = expression(S(z)),
        main = expression(paste(beta, " = 0.22")))
abline(h = synth[[1]]$sum, col = 4, lty = 2)

The trace plot of \(S(\mathbf{z})\) shows rapid mixing at \(\beta = 0.22\): within a dozen iterations the chain is oscillating around the true value (dashed line). Pre-computed results for all four values of \(\beta\) are stored as , , , and . The vignette provides a systematic comparison of convergence across these values and contrasts with the no-data sampler .

Summary

The package provides a unified interface to five algorithms for Bayesian inference in the hidden Potts model, supporting 2D and 3D images, irregular masks, external field priors, and adaptive MCMC. The computational core, written in via and , enables efficient processing of images with hundreds of thousands of pixels. For large-scale applications, the PFAB algorithm is recommended; the exchange algorithm provides exact inference at greater computational cost; and pseudolikelihood gives a fast approximation suitable for exploratory analysis. Theoretical background for all five algorithms is given in .

true