bayesImageS: an R package for Bayesian image analysis Matt Moores @MooresMt https://mooresm.github.io/bayesImageS NIASRA, University of Wollongong, NSW, Australia useR! 2018 1 / 22 

Image Segmentation Classifying the observed pixel intensities y := {yi }n i=1 ∈ Rn according to a set of latent labels z := {zi }n i=1 ∈ [1, . . . , k]n Assuming additive Gaussian noise: yi | zi = j ∼ N µj, σ2 j The resulting augmented likelihood is a mixture of Gaussians: p(y, z | µ, σ2, β) = n i=1 p yi | µzi , σ2 zi p(zi | z\i , β) 2 / 22 Rydén & Titterington (1998) J. Comput. Graph. Stat. 7(2): 194–211. 

Synthetic Data z <- sample.int(3, 49, TRUE) mu <- c(-1, 0, 1) sigma <- rep(0.3, 3) y <- rnorm(49, mu[z], sigma[z]) 2 4 6 Y −1 0 1 Z 0.0 0.2 0.4 0.6 density Component Means −1 0 1 3 / 22 

Mixture Model library(bayesImageS) lPriors <- list(k=3, mu=c(-1,0,1), mu.sd=rep(0.25,3), sigma=rep(0.5,3), sigma.nu=rep(5,3), lambda=rep(1,3)) mix.fit <- gibbsGMM(y, priors=lPriors) mu[1] mu[2] mu[3] 0 200 400 600 800 1000 0 200 400 600 800 1000 0 200 400 600 800 1000 0.4 0.8 1.2 1.6 −0.8 −0.4 0.0 0.4 −1.50 −1.25 −1.00 −0.75 −0.50 mu[1] mu[2] mu[3] −1.4 −1.2 −1.0 −0.8 −0.6 −0.4 −0.8 −0.4 0.0 0.4 0.4 0.8 1.2 1.6 Estimated True 4 / 22 

Pixel Classification Ground Truth Classification 5 / 22 

Potts model In many images, neighbouring pixels tend to share the same label. This tendency can be represented using a hidden Markov random field: p(zi | z\i , β) = exp{β i∼ δ(zi , z )} k j=1 exp{β i∼ δ(j, z )} where: β is the inverse temperature parameter, i ∼ are the neighbours of pixel i, and δ(x, y) is the Kronecker delta function. 6 / 22 Potts (1952) Proceedings of the Cambridge Philosophical Society 48(1) 

Inverse Temperature 7 / 22 

Synthetic Data Using the Swendsen-Wang algorithm to simulate from the Potts model with known β: mask <- matrix(1, nrow=7, ncol=7) neigh <- getNeighbors(mask, c(2,2,0,0)) blocks <- getBlocks(mask, 2) k <- 3 beta <- 0.88 res.sw <- swNoData(beta, k, neigh, blocks, niter=200) z <- matrix(max.col(res.sw$z)[1:nrow(neigh)], nrow=nrow(mask)) y <- rnorm(49, mu[z], sigma[z]) 8 / 22 

Fitting the Model lPriors <- list(k=3, mu=c(-1,0,1), mu.sd=rep(0.25,3), sigma=rep(0.5,3), sigma.nu=rep(5,3), beta=c(0,2)) mh <- list(algorithm="ex",bandwidth=0.001,auxiliary=200) mrf.fit <- mcmcPotts(y, neigh, blocks, lPriors, mh, niter=500, nburn=250) 0.0 0.5 1.0 0 100 200 300 400 500 beta beta 0.0 0.5 1.0 Estimated True 9 / 22 

Pixel Classification Ground Truth Classification 10 / 22 

Intractable Likelihood p(β|z) = C−1(β)eβS(z)π(β) β C−1(β)eβS(z)π(dβ) (1) The normalising constant has computational complexity O(nkn), since it involves a sum over all possible combinations of the labels z ∈ Z: C(β) = z∈Z eβS(z) (2) S(z) is the sufficient statistic of the Potts model: S(z) = i∼ ∈E δ(zi , z ) (3) where E is the set of all unique neighbour pairs. 11 / 22 

Algorithms for β The argument mh$alg selects which algorithm to use for sampling from p(β|z): pseudolikelihood path sampling (thermodynamic integration) approximate exchange algorithm (AEA) approximate Bayesian computation (ABC-MCMC and ABC-SMC) Bayesian indirect likelihood (BIL) 12 / 22 Moores, Drovandi, Mengersen & Robert (2015) Stat. Comput. 25(1): 23–33. Moores, Pettitt & Mengersen (2015; v2 2018) arXiv:1503.08066 [stat.CO] 

Satellite Remote Sensing −3050000 −3045000 −3040000 −3035000 −3030000 490000 495000 500000 505000 510000 515000 −1.0 −0.5 0.0 0.5 1.0 NDVI Frequency −1.0 −0.5 0.0 0.5 1.0 0 50000 100000 150000 13 / 22 

Image segmentation with k = 5 labels 14 / 22 

Posterior Distributions Pixel Intensity Density −1.0 −0.5 0.0 0.5 1.0 0.0 0.5 1.0 1.5 2.0 1.193 1.194 1.195 1.196 0 200 400 600 800 N = 450 Bandwidth = 0.0001203 Density AEA BIL ABC 15 / 22 

Segmentation of Anatomical Structures 16 / 22 Radiography courtesy of Cathy Hargrave, Radiation Oncology Mater Centre, Queensland Health 

External Field Prior p(zi |zi∼ , β, µ, σ2, yi ) = exp {αi,zi + π(αi,zi )} k j=1 exp {αi,j + π(αi,j)} π(zi |zi∼ , β) (4) Isotropic translation: π(αi,j) = log    1 nj h∈j φ ∆(h, i)|µ∆ = 1.2, σ2 ∆ = 7.32    (5) where nj is the number of voxels in object j h ∈ j are the voxels in object j ∆(u, v) is the Euclidean distance between the coordinates of pixel u and pixel v µ∆, σ2 ∆ are parameters that describe the level of spatial variability of the object j 17 / 22 

Patient- and Organ-Specific Prior αi (prostate) ∼ MVN       0.1 −0.5 0.2    ,    4.12 0 0 0 2.92 0 0 0 0.92       18 / 22 

Seminal Vesicles αi (SV) ∼ MVN       1.2 −0.7 −0.9    ,    7.32 0 0 0 4.52 0 0 0 1.92       19 / 22 

Gaussian Random Field 20 / 22 

Results −300 −250 −200 −150 150 200 250 300 right−left (mm) posterior−anterior (mm) −300 −250 −200 −150 150 200 250 300 right−left (mm) posterior−anterior (mm) 21 / 22 Moores, Hargrave, Deegan, Poulsen, Harden & Mengersen (2015) CSDA 86: 27–41. 

Summary bayesImageS can be used for segmentation of 2D and 3D images: Update labels using chequerboard Gibbs sampling or Swendsen-Wang Posterior for β using the exchange algorithm, ABC, pseudolikelihood, path sampling, or BIL Fast implementation using RcppArmadillo Parallelism using OpenMP PkgDown documentation: https://mooresm.github.io/bayesImageS Version 0.5-2 currently available on CRAN: https://CRAN.R-project.org/package=bayesImageS 22 / 22