{{ message }}

Instantly share code, notes, and snippets.

# msilvestro/!image_restore.R

Last active Apr 26, 2021
Coupling from the past (CFTP) for image restoration
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 # Created by Matteo Silvestro for an essay about CFTP algorithm # 16/06/2016 # CC BY library(pixmap) # to extract a matrix of zeros and ones from a binary image library(grid) # to display a matrix of zeros and ones in R # I define as one-line-matrix a square matrix with all elements in a single vector. # It is easy to get a virtual matrix from such vector by setting number of cols/rows = sqrt(length of vector) # draw the one-line-matrix of -1 and +1 as an image in R drawImage <- function(image) { N <- sqrt(length(image)) # compute the size of the squared image grid.newpage() # delete all previous stuff on the drawing surface imagen <- matrix(1-(image+1)/2, ncol=N) # change of values from -1 -> 0, +1 -> 1 grid.raster(imagen, interp=F) # with a matrix of 0 and 1, the grid package can draw the image correctly } # return a vector with the neighbour elements of the k-th element of the one-line-matrix x neighbour <- function(x, k) { N <- sqrt(length(x)) # compute the size of the virtual matrix neigh <- NULL # this will be the vector of all neighbours # here there is a bit of algebra to get the right elements, remember that x is a one-line-matrix # (one-line-matrix vector index) k -> (floor(k/N), k % N) (virtual matrix index) # soppose (i, j) are the coordinate of x[k], then the neighbours are # (i, j-1), (i, j+1), (i-1, j), (i+1, j) if 1 <= i, j <= N if (k-N >= 1) neigh <- c(neigh, x[k-N]) if (k+N <= N*N) neigh <- c(neigh, x[k+N]) if (k %% N != 1) neigh <- c(neigh, x[k-1]) if (k %% N != 0) neigh <- c(neigh, x[k+1]) return(neigh) } fullcond <- function(x, y, k, beta, p) { return(1 / ( 1 + exp(-2*beta*sum(neighbour(x, k)) - log((1-p)/p)*y[k] ) )) } # apply the Propp-Wilson algorithm rimage <- function(y, beta, p, N, seed) { Time <- -0.5 # actual starting time is -1, since it is doubled immediately xmax <- rep(1, N*N) xmin <- rep(-1, N*N) while (!all(xmax == xmin)) { # continue until convergence does occur xmax <- rep(1, N*N) # image with all black pixels xmin <- rep(-1, N*N) # image with all white pixels Time <- Time*2 # double the time, as suggested by Propp-Wilson set.seed(seed) # seed must be fixes, since we use random mappings utot <- runif(N*N*(-Time)) # create the needed random mappings random numbers for (t in Time:-1) { u <- utot[ (N*N*(-t-1)+1) : (N*N*(-t)) ] # extract the random numbers to be used is this iteration for (k in 1:(N*N)) { # apply the Gibbs sampler, updating a pixel at a time xmax[k] = (u[k] < fullcond(xmax, y, k, beta, p))*2-1 xmin[k] = (u[k] < fullcond(xmin, y, k, beta, p))*2-1 } } } print(c(seed, Time)) return(xmax) } restoreImage <- function(file, p, beta, simn) { image <- read.pnm(file, cellres=1) # cellres is to avoid the warning image <- 1-as.vector(image@grey) N <- sqrt(length(image)) # make an array with the binary original binary image, with +1 as black and -1 as white x <- image*2-1 drawImage(x) # degrade the image with a Bernoulli noise with p as probability to degrade y <- x * (rbinom(N*N, 1, 1-p)*2-1) # apply Bernoulli noise drawImage(y) # apply the restoring algorithm, based on Gibbs sampler with CFTP sample <- NULL for (i in 1:simn) { sample <- c(sample, rimage(y, beta, p, N, i)) # make one long vector with all samples } sample <- matrix(sample, nrow=simn, ncol=N*N, byrow=T) # and then split it in a matrix with as each row a sample # MPM (marginal posterior mode) estimator mode <- rep(0, N*N) for (i in 1:(N*N)) { mode[i] <- (sum(sample[,i]) > 0)*2-1 # make the mode of each pixel, i.e. for each column of the matrix take the most frequent value (+1 or -1) } drawImage(mode) # compute the misclassification rate mis <- sum(mode != x)/(N*N) # return a list with all informations return(list(original=x, degraded=y, restored=mode, error=mis)) } # ptm <- proc.time() # result <- restoreImage("img/penguin.pbm", p=0.3, beta=0.45, simn=1) # time <- proc.time() - ptm

# Coupling from the past (CFTP) for image restoration

Markov chain Monte Carlo methods allow to compute difficult integrals and sample from complex distributions, opening new possibilities. Gibbs sampler gained a lot of popularity, dealing with multidimensional problems. The main issue is to be able to sample directly from the stationary distribution of the Markov chain, that is to get perfect samples.

An efficient algorithm to accomplish that is the coupling from the past (CFTP) algorithm. In this essay, we will describe and prove the functionality of CFTP algorithm. Furthermore, a practical application in image restoration is given, applying CFTP to a Gibbs sampler.