Bayesian Methods for Image Segmentation using a Potts Model

Various algorithms for segmentation of 2D and 3D images, such as computed tomography and satellite remote sensing. This package implements Bayesian image analysis using the hidden Potts model with external field prior of Moores et al. (2015) . Latent labels are sampled using chequerboard updating or Swendsen-Wang. Algorithms for the smoothing parameter include pseudolikelihood, path sampling, the exchange algorithm, approximate Bayesian computation (ABC-MCMC and ABC-SMC), and the parametric functional approximate Bayesian (PFAB) algorithm. Refer to and for further details.

cranversion rstudio mirrordownloads

bayesImageS implements algorithms for segmentation of 2D and 3D images, such as computed tomography (CT) and satellite remote sensing. This R package provides functions for Bayesian image analysis using a hidden Potts/Ising model with external field prior. Latent labels are updated using chequerboard Gibbs sampling or Swendsen-Wang. Algorithms for the smoothing parameter include:

  • pseudolikelihood
  • path sampling (thermodynamic integration)
  • approximate exchange algorithm (AEA)
  • approximate Bayesian computation (ABC-MCMC and ABC-SMC)
  • Bayesian indirect likelihood (BIL), including the parametric functional approximate Bayesian (PFAB) algorithm

Stable releases, including binary packages for Windows & Mac OS, are available from CRAN:


The current development version can be installed from Bitbucket:


Example Usage

To generate synthetic data for a known value of (\beta):

mask <- matrix(1,3,3)
neigh <- getNeighbors(mask, c(2,2,0,0))
blocks <- getBlocks(mask, 2)
k <- 3
beta <- 0.7
res.sw <- swNoData(beta, k, neigh, blocks, niter=200)
z <- matrix(max.col(res.sw$z)[1:nrow(neigh)], nrow=nrow(mask))
image(z, xaxt = 'n', yaxt='n', col=rainbow(k), asp=1)

Now add some Gaussian noise to the labels, according to the prior:

priors <- list()
priors$k <- k
priors$mu <- c(-2, 0, 2)
priors$ <- rep(0.5,k)
priors$sigma <- rep(0.25,k)
priors$ <- rep(3, k)
priors$beta <- c(0,1.3*log(1 + sqrt(k)))
m0 <- sort(rnorm(priors$k,priors$mu,priors$
SS0 <- priors$*priors$sigma^2
s0 <- 1/sqrt(rgamma(priors$k,priors$,SS0/2))
l <- as.vector(z)
y <- m0[l] + rnorm(nrow(neigh),0,s0[l])
levelplot(matrix(y, nrow=nrow(mask)))

Image segmentation using ABC-SMC:

res.smc <- smcPotts(y, neigh, blocks, priors=priors)
#> Initialization took 5sec
#> Iteration 1
#> previous epsilon 7 and ESS 10000 (target: 9500)
#> Took 0sec to update epsilon=2.625 (ESS=9505.29)
#> Took 4sec for 8918 RWMH updates (bw=0.497509)
#> Took 0sec for 10000 iterations to calculate S(z)=7
#> Iteration 2
#> previous epsilon 2.625 and ESS 9505.29 (target: 9030.02)
#> Took 8sec to update epsilon=1 (ESS=7970.86)
#> Took 4sec for 7671 RWMH updates (bw=0.466951)
#> Took 0sec for 10000 iterations to calculate S(z)=6
#> Iteration 3
#> previous epsilon 1 and ESS 7970.86 (target: 7572.32)
#> Took 8sec to update epsilon=4.66632e-302 (ESS=7949.67)
#> Took 3sec for 7968 RWMH updates (bw=0.466673)
#> Took 1sec for 10000 iterations to calculate S(z)=7
# pixel classifications
pred <- res.smc$alloc/rowSums(res.smc$alloc)
predMx <- as.raster(array(pred, dim=c(nrow(mask),ncol(mask),3)))
rasterImage(t(predMx)[nrow(mask):1,], 0.5, 0.5, 3.5, 3.5, interpolate = FALSE)

Note that CODA ignores the particle weights, so we need to resample to obtain accurate HPD intervals. This step is not usually necessary and does introduce some noise due to duplication of particles. Depending on how many SMC iterations have been performed, one or more resampling steps might have already been done (but not in this specific example).

seg <- max.col(res.smc$alloc) # posterior mode (0-1 loss)
all.equal(seg, l)
#> [1] TRUE
# filter weights to remove Ninf, NaN
w <- res.smc$wt
w[] <- 0
plot(density(res.smc$beta, weights=w),main=expression(paste("Posterior for ",beta)))
abline(v=log(1 + sqrt(k)),lty=3,col=2) # critical point

res.res <- testResample(res.smc$beta, w, cbind(res.smc$mu, res.smc$sigma))
#> Took 0sec to resample 10000 particles
res.coda <- mcmc(cbind(res.res$pseudo, res.res$beta))
varnames(res.coda) <- c(paste("mu",1:k), paste("sd",1:k), "beta")
#>            lower      upper
#> mu 1 -2.56248597 -1.9038983
#> mu 2 -0.50953511  0.4881316
#> mu 3  1.37933533  2.7939284
#> sd 1  0.11321785  0.4984628
#> sd 2  0.22470066  0.9154286
#> sd 3  0.11011734  0.9105506
#> beta  0.09338037  1.2852256
#> attr(,"Probability")
#> [1] 0.95
#> [1] -2.279614  0.156580  2.208122
#> [1] 0.2785175 0.6092555 0.3153176
#> [1] 0.7


bayesImageS 0.6-0

  • updated citation info for Bayesian Analysis (2018) paper
  • added a hex logo to the README
  • removed duplicate vignette entry

bayesImageS 0.5-3

  • modified vignette swNoData to only use PottsUtils conditionally
  • deleted vignette PFAB due to problems with Stan and Boost on gannet/ripley/tests-ATLAS

bayesImageS 0.5-2

  • changed maintainer email address due to new academic affiliation
  • reorganised package vignettes
  • NAMESPACE is now generated by roxygen2
  • added LICENSE file for GPL 2

bayesImageS 0.5-1

  • new vignette introducing the parametric functional approximate Bayesian (PFAB) algorithm
  • removed all @seealso from Roxygen docs due to CRAN WARNings because PottsUtils has now been permanently removed from the CRAN repository

bayesImageS 0.5-0

New features

  • improvements to Bayesian indirect likelihood (BIL) algorithm (mh$alg="aux"):

    • MH hyperparameter Ecrit for 2nd order phase transition when k > 4
    • MH hyperparameter factor to inflate the variance of the surrogate model
    • separate MH hyperparameters Vmax1,Vmax2 for $\beta < \beta_{crit}$ and $\beta > \beta_{crit}$, respectively
  • option mh$sort to impose an ordering constraint on the component means

Bug fixes

  • removed OpenMP options from Makevars to fix problems with RcppArmadillo. Users should edit ~/.R/Makevars to enabel multithreading.

  • no longer wait for mh$auxiliary iterations before beginning to update $\beta$ using the exchange algorithm or ABC. For $\beta > \beta_{crit}$, this increased the propensity of the Gibbs sampler to become stuck in a suboptimal local mode, far from the true parameter value.

  • added jss.bst and jss.cls to .Rbuildignore

bayesImageS 0.4-1

  • removed the deprecated slices parameter from mcmcPotts()

  • some minor edits to the Roxygen documentation


  • removed the R package PottsUtils from Suggested packages; added mcmcse

  • moved README-chunkname.png images from root directory to inst/images

bayesImageS 0.4-0

New features

  • Functions getNeighbors(), getEdges(), getBlocks() contributed by Dai Feng. These are currently exact copies of the equivalent functions in the R package PottsUtils. The reason for duplicating this code is to avoid picking up dependencies on miscF, R2jags, BRugs, etc.

  • Added new option to swNoData and mcmcPottsNoData to initialise the labels using either random or deterministic values. This can be useful when comparing the convergence of the Swendsen-Wang and Gibbs sampling algorithms to known distributions (e.g. when $\beta$=0 or $\infty$)

  • Additional documentation in with an example of usage

Bug fixes

bayesImageS 0.3-3

  • First version released on CRAN

bayesImageS 0.1-21

Reference manual

It appears you don't have a PDF plugin for this browser. You can click here to download the reference manual.


0.6-0 by Matt Moores, 2 years ago,

Report a bug at

Browse source code at

Authors: Matt Moores [aut, cre] , Kerrie Mengersen [aut, ths] , Dai Feng [ctb]

Documentation:   PDF Manual  

Task views: Bayesian Inference, Medical Image Analysis

GPL (>= 2) | file LICENSE license

Imports Rcpp

Suggests mcmcse, coda, PottsUtils, rstan, knitr, rmarkdown, lattice

Linking to Rcpp, RcppArmadillo

See at CRAN