Image denoising by MCMC

Using Markov Chain Monte Carlo in conjunction with an Ising model to clean up a noisy binary image

In this article, I will be demonstrating the use of Markov Chain Monte Carlo to denoise a binary image.

Markov Chain Monte Carlo or MCMC for short refers to a class of techniques used for estimating a probability distribution by sampling from it. The various techniques comprising MCMC are differentiated from each other based on the method used for drawing samples. Some of the more well known MCMC techniques are Metropolis — Hastings, Gibbs sampling and Hamiltonian Monte Carlo. The technique that I will be using is Gibbs sampling.

Gibbs sampling is a method for sampling from a multivariate distribution given all other variable are kept the same.

For example, if there is distribution with only two variables x1 and x2, the sampling is as follows:

Example distribution

It starts at point 1 and then when the second point is to be sampled, it goes as:

P(x2|x1), it looks for the next sample on the same line as x1 (keeping it constant)

Following which the third point is sampled as follows:

P(x1|x2) , it looks for the next sample on the same line as x2 (keeping it constant)

And so on, the sampling process continues for a set number of points allowing it to traverse the entire space.

Ising Model

The Ising model is a mathematical model corresponding to a square lattice used for modelling phase transitions. Each element in the lattice can exist in two discrete states and can be represented by +1 and -1. Each element exerts an influence on all its neighbouring elements and tries to reach a state of equilibrium where all elements exist in the same state.

Application:

Conceptualized model of image

The binary image can be thought of in the form of a lattice where each pixel represents one element. The pixel’s state can be represented as 1 or -1 depending on the colour of the pixel. The image can be imagined consisting of two layers, the underlying layer which represents the true image and the above layer representing the noise. The Gaussian noise is said to be superimposed upon the image with it matching in some places to the actual image and in some places taking the opposite values.

The Ising model is applied to the above layer consisting of noise. The noise is influenced by each of its neighbours depending on how closely they are tied to each other, represented by the edge potential. The closer they are tied to each other, the more they try to be in the same state. The formula for the edge potential is given by:

exp(J, X ₐ , X ₙ )

Here J is the coupling constant which denotes how closely neighbours are tied to each other. Xₐ represents the pixel under consideration and Xₙ represents the observed values of its neighbours.

A Gaussian observation model is used to model the relation between the noise and the pre-existing pixels. The noise is said to be a function of the actual underlying pixel value and of the standard deviation from it. It can be represented by:

N(Yₐ | Xₐ, σ ²)

Here Yₐ represents the observed value of the pixel, Xₐ the pixel under consideration and σ² the standard deviation.

It is the interplay between these two forces that decide the eventual value of the pixel. The coupling constant tries to keep all the neighbouring pixels(noise) in the same state or colour. The gaussian model tries to match the actual colour of the pixel to the observed colour(noise). The influence of these two forces can be manipulated by changing the values of J and that of σ².

Gibbs sampling applied would sample each pixel conditioned on all its neighbours and the underlying true value. It would then fix this value and proceed to the next element and repeat the same operation. One iteration would be completed when it has finished traversing through the entire lattice. Based on the number of iterations, the final image quality can differ.

Xₐ the pixel under consideration can take either the values +1 or -1. The above expression gives the likelihood of Xₐ taking the values +1 or -1.

The numerator gives the likelihood of Xₐ to be, for instance, +1, and checking its relationship with the observed value of the pixel Yₐ, and with the observed values of its neighbours Xₙ.

It is divided by the sum of the likelihood of it taking the values +1 and -1 each to give the probability of the actual value of Xₐ to be +1.