bayesImageS v0.5-3

0

Monthly downloads

0th

Percentile

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) <doi:10.1016/j.csda.2014.12.001>. 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 Bayesian indirect likelihood (BIL). Refer to <doi:10.1007/s11222-014-9525-6> and <arXiv:1503.08066> for further details.

Readme

bayesImageS

cran version rstudio mirror downloads

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

Installation Instructions

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

install.packages("bayesImageS")

The current development version can be installed from Bitbucket:

devtools::install_git("https://bitbucket.org/Azeari/bayesimages/")

Example Usage

To generate synthetic data for a known value of β:

set.seed(1234)
library(bayesImageS)
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$mu.sd <- rep(0.5,k)
priors$sigma <- rep(0.25,k)
priors$sigma.nu <- rep(3, k)
priors$beta <- c(0,1.3*log(1 + sqrt(k)))

m0 <- sort(rnorm(priors$k,priors$mu,priors$mu.sd))
SS0 <- priors$sigma.nu*priors$sigma^2
s0 <- 1/sqrt(rgamma(priors$k,priors$sigma.nu/2,SS0/2))
l <- as.vector(z)
y <- m0[l] + rnorm(nrow(neigh),0,s0[l])
library(lattice)
levelplot(matrix(y, nrow=nrow(mask)))

Image segmentation using ABC-SMC:

res.smc <- smcPotts(y, neigh, blocks, priors=priors)
#> Initialization took 6sec
#> Iteration 1
#> previous epsilon 7 and ESS 10000 (target: 9500)
#> Took 1sec to update epsilon=2.625 (ESS=9505.29)
#> Took 5sec for 8918 RWMH updates (bw=0.497509)
#> Took 1sec for 10000 iterations to calculate S(z)=7
#> Iteration 2
#> previous epsilon 2.625 and ESS 9505.29 (target: 9030.02)
#> Took 9sec to update epsilon=1 (ESS=7970.86)
#> Took 6sec for 7671 RWMH updates (bw=0.466951)
#> Took 1sec for 10000 iterations to calculate S(z)=6
#> Iteration 3
#> previous epsilon 1 and ESS 7970.86 (target: 7572.32)
#> Took 9sec to update epsilon=4.66632e-302 (ESS=7949.67)
#> Took 6sec 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)))
plot(c(0.5,3.5),c(0.5,3.5),type='n',xaxt='n',yaxt='n',xlab="",ylab="",asp=1)
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[is.na(w)] <- 0
plot(density(res.smc$beta, weights=w),main=expression(paste("Posterior for ",beta)))
abline(h=0,lty=3)
abline(v=beta,lty=2,col=4)
abline(v=log(1 + sqrt(k)),lty=3,col=2) # critical point


library(coda)
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")
HPDinterval(res.coda)
#>            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
m0
#> [1] -2.279614  0.156580  2.208122
s0
#> [1] 0.2785175 0.6092555 0.3153176
beta
#> [1] 0.7

Functions in bayesImageS

Name Description
testResample Test the residual resampling algorithm.
mcmcPotts Fit the hidden Potts model using a Markov chain Monte Carlo algorithm.
mcmcPottsNoData Simulate pixel labels using chequerboard Gibbs sampling.
res2 Simulation from the Potts model using single-site Gibbs updates.
res Simulation from the Potts model using single-site Gibbs updates.
bayesImageS Package bayesImageS
exactPotts Calculate the distribution of the Potts model using a brute force algorithm.
sufficientStat Calculate the sufficient statistic of the Potts model for the given labels.
swNoData Simulate pixel labels using the Swendsen-Wang algorithm.
getNeighbors Get Neighbours of All Vertices of a Graph
gibbsGMM Fit a mixture of Gaussians to the observed data.
res5 Simulation from the Potts model using single-site Gibbs updates.
smcPotts Fit the hidden Potts model using approximate Bayesian computation with sequential Monte Carlo (ABC-SMC).
gibbsNorm Fit a univariate normal (Gaussian) distribution to the observed data.
initSedki Initialize the ABC algorithm using the method of Sedki et al. (2013)
res3 Simulation from the Potts model using single-site Gibbs updates.
res4 Simulation from the Potts model using single-site Gibbs updates.
getBlocks Get Blocks of a Graph
getEdges Get Edges of a Graph
synth Simulation from the Potts model using Swendsen-Wang.
No Results!

Vignettes of bayesImageS

Name
Background.Rmd
InverseTemperature.bib
bcrit2d-eps-converted-to.pdf
bcrit2d.eps
bcrit2d_sd-eps-converted-to.pdf
bcrit2d_sd.eps
bcrit3d-eps-converted-to.pdf
bcrit3d.eps
bcrit3d_sd-eps-converted-to.pdf
bcrit3d_sd.eps
beta_ct-eps-converted-to.pdf
beta_ct.eps
ct_hist-eps-converted-to.pdf
ct_hist.eps
elapsed2D-eps-converted-to.pdf
elapsed2D.eps
elapsed3D-eps-converted-to.pdf
elapsed3D.eps
elapsed_ct-eps-converted-to.pdf
elapsed_ct.eps
exact_exp_k-eps-converted-to.pdf
exact_exp_k.eps
exact_exp_n-eps-converted-to.pdf
exact_exp_n.eps
exact_var_k-eps-converted-to.pdf
exact_var_k.eps
exact_var_n-eps-converted-to.pdf
exact_var_n.eps
fanbeam.jpeg
introduction.Rnw
jsslogo.jpg
mcmcPotts.Rmd
mcmcPottsNoData.Rmd
mcmc_err2D_ABC-eps-converted-to.pdf
mcmc_err2D_ABC.eps
mcmc_err2D_MAVM-eps-converted-to.pdf
mcmc_err2D_MAVM.eps
mcmc_err2D_PL-eps-converted-to.pdf
mcmc_err2D_PL.eps
mcmc_err2D_TI-eps-converted-to.pdf
mcmc_err2D_TI.eps
ndvi_hist-eps-converted-to.pdf
ndvi_hist.eps
ndvi_image.jpeg
path2d-eps-converted-to.pdf
path2d.eps
path3d-eps-converted-to.pdf
path3d.eps
pl_exp_n12k3-eps-converted-to.pdf
pl_exp_n12k3.eps
pl_sd_n12k3-eps-converted-to.pdf
pl_sd_n12k3.eps
refs.bib
swNoData.Rmd
No Results!

Last month downloads

Details

Type Package
Date 2018-08-29
License GPL (>= 2) | file LICENSE
URL https://bitbucket.org/Azeari/bayesimages, https://mooresm.github.io/bayesImageS/
BugReports https://bitbucket.org/Azeari/bayesimages/issues
LinkingTo Rcpp, RcppArmadillo
VignetteBuilder knitr
RoxygenNote 6.1.0
NeedsCompilation yes
Packaged 2018-08-30 01:39:34 UTC; mmoores
Repository CRAN
Date/Publication 2018-08-30 05:35:00 UTC

Include our badge in your README

[![Rdoc](http://www.rdocumentation.org/badges/version/bayesImageS)](http://www.rdocumentation.org/packages/bayesImageS)