Fast Bayesian Restoration of Poisson Corrupted Images with INLA

TL;DR

This paper introduces a fast Bayesian image restoration technique for Poisson noise using INLA, significantly improving computational speed while maintaining accuracy, especially for images modeled as ICAR.

Contribution

The paper presents a novel application of INLA for Poisson noise removal in images, offering faster computation compared to existing methods like MCMC and belief propagation.

Findings

INLA-based method is faster than MCMC and belief propagation.

The approach maintains high accuracy in restoring Poisson-corrupted images.

Effective for images modeled as ICAR.

Abstract

Photon-limited images are often seen in fields such as medical imaging. Although the number of collected photons on an image sensor statistically follows Poisson distribution, this type of noise is intractable, unlike Gaussian noise. In this study, we propose a Bayesian restoration method of Poisson corrupted image using Integrated Nested Laplace Approximation (INLA), which is a computational method to evaluate marginalized posterior distributions of latent Gaussian models (LGMs). When the original image can be regarded as ICAR (intrinsic conditional auto-regressive) model reasonably, our method performs very faster than well-known ones such as loopy belief propagation-based method and Markov chain Monte Carlo (MCMC) without decreasing the accuracy.

Figures

Tables

Table 1. TABLE I: Comparison of parameter estimating methods with PSNR, SSIM, and CPU time.

Image	Method	PSNR	SSIM	Time(s)
Lena	INLA	22.24	0.966	226.8
	MCMC	24.97	0.968	2083.5
	LBP	21.77	0.928	18794.7
Boat	INLA	21.76	0.947	209.5
	MCMC	24.03	0.960	1585.3
	LBP	21.50	0.915	18891.5
Cameraman	INLA	10.88	0.916	113.5
	MCMC	21.82	0.956	1360.7
	LBP	19.61	0.923	18809.8

Equations

$$p(\mathbf{x})=\mathcal{N}(\mathbf{x}|\mathbf{\mu }({\mathbf{\theta }}_{\mathrm{lat}}),{\mathbf{\Sigma }}^{-1}({\mathbf{\theta }}_{\mathrm{lat}})),$$

$$p(\mathbf{y})=\prod _{i=1}^{n}p({\mathbf{y}}_i|\mathbf{x},{\mathbf{\theta }}_{\mathrm{obs}}),$$

$$\mathbf{\theta }=\{{\mathbf{\theta }}_{\mathrm{lat}},{\mathbf{\theta }}_{\mathrm{obs}}\}.$$

$$p(\mathbf{\theta }|\mathbf{y})\propto \frac{p(\mathbf{\theta })p(\mathbf{x}|\mathbf{\theta })p(\mathbf{y}|\mathbf{x},\mathbf{\theta })}{p(\mathbf{x}|\mathbf{\theta },\mathbf{y})}$$

$$p({\mathbf{x}}_i|\mathbf{y})=\int p({\mathbf{x}}_i|\mathbf{y},\mathbf{\theta })p(\mathbf{\theta }|\mathbf{y})\mathrm{d}\mathbf{\theta }.$$

$$\tilde{p}(\mathbf{\theta }|\mathbf{y})\propto \frac{p(\mathbf{\theta })p(\mathbf{x}|\mathbf{\theta })p(\mathbf{y}|\mathbf{x},\mathbf{\theta })}{\tilde{p}(\mathbf{x}|\mathbf{\theta },\mathbf{y})},$$

$$p({\mathbf{x}}_i|\mathbf{y})\approx \sum _{h=1}^{H}p({\mathbf{x}}_i|\mathbf{y},{\mathbf{\theta }}_h)p({\mathbf{\theta }}_h|\mathbf{y}){\Delta }_h,$$

$$p(x_i)=\mathcal{N}\left(x_i|\frac{1}{|C_i|}\sum _{j\in C_i}{x}_j,\frac{{\sigma }^2}{|C_i|}\right),$$

$$p(x_i)=\mathcal{N}\left(x_i|\frac{1}{|C_i|+d}\sum _{j\in C_i}{x}_j,\frac{{\sigma }^2}{|C_i|+d}\right),$$

$$p(x_l)=\mathcal{N}\left(x_l|{\mu }_l,\frac{{\sigma }^2}{4+d}\right),$$

$${\mu }_l=\frac{1}{4+d}(x_{l-1}+x_{l+1}+x_{l-L}+x_{l+L}).$$

$$p(y_i)=\mathrm{Poisson}(y_i|x_i)$$

$$p(\mathbf{y})=\prod _{i=1}^n\mathrm{Poisson}(y_i|x_i).$$

$$x_i=\frac{{\lambda }_{\max }-{\lambda }_{\min }}{I_{\max }-I_{\min }}(I_i-I_{\min })+{\lambda }_{\min },$$

$$PSNR=10{\log }_{10}\frac{(\max \{\mathbf{G},\mathbf{H}\}-\min \{\mathbf{G},\mathbf{H}\}{)}^2}{\mathrm{MSE}(\mathbf{G},\mathbf{H})}$$

$$MSE(\mathbf{G},\mathbf{H})=\frac{1}{n}\sum _i^n(G_i-H_i{)}^2,$$

$$SSIM=\frac{(2{\mu }_G{\mu }_H+c_1)(2{\sigma }_{G,H}+c_2)}{({\mu }_G^2+{\mu }_H^2+c_1)({\sigma }_G^2+{\sigma }_H^2+c_2)},$$

Full text

I Introduction

Estimating the original image from a noisy observation is one of the representative issues in the field of image processing. So far, a lot of image denoising methods have been proposed, but many of them such as Wiener filter require the assumption of additive white Gaussian noise. However, for example, in medical X-ray imaging systems, the number of detected photons stochastically fluctuates following Poisson distribution[1].

From a Bayesian viewpoint, it is natural to adopt Poisson observation modeling explicitly, and then solve the inverse problem to deal with such images. In fact, some Bayesian methods to restore Poisson corrupted images have been already proposed. Le et al. introduced a variational model with total-variation regularization[2]. Lefkimmiatis et al. applied quadtree decomposition and then estimated parameters using expectation-maximization (EM) algorithm to handle Poisson noise[3]. Shouno derived the fixed point equations of loopy belief propagation (LBP) by approximating Poisson distribution with binomial distribution[4]. Furthermore, Tachella et al. compared the performances of multiple MCMC methods to restore Poisson corrupted images [5].

These methods are frequently used to evaluate the posterior (or its point estimate) of each pixel. In more general applications, Markov chain Monte Carlo (MCMC) methods succeeded due to their usability and accuracy. On the other hand, recently, integrated nested laplace approximation (INLA) is proposed by Rue et al. [6], then the validity has been reported mainly in the fields of spatial statistics and epidemiology [7]. For example, INLA is used disease mapping[8], analysing spatio-temporal evolutions of Ebola virus[9], and Bayesian outbreak detection[10]. Also there is a similar research to ours, INLA-based super-resolution method has been proposed [11].

INLA is a powerful framework to evaluate posterior marginals on latent Gaussian models (LGMs). LGM is a class of statistical models which includes a lot of standards, for example, auto-regressive (AR), conditional auto-regressive (CAR), and generalized linear mixture (GLM) models[7, 12]. INLA consists of two layers of units: latent variables following Gaussian Markov random field (GMRF) and their observations. Latent variables requires normality on their interactions, but observation processes do not. INLA is based on Laplace approximations and numerical integrations as the name suggests. As the interface for using INLA, R-INLA is provided by R-INLA project 111http://www.r-inla.org/.

In this study, we try applying INLA to image restoring as the application to the new field, through dealing with Poisson corrupted images. When the original image seems to be reasonable to assume intrinsic CAR (ICAR) model, which means almost of adjacent pixels do not vary sharply, the proposed method obtained equal to the result by an MCMC simulation rapidly.

This paper organized as following: first, we give a brief description of INLA in Section 2. Then, we explain our model to infer the original images in Section 3, and to evaluate proposed method we show configurations and results of the computational simulations in Section 4. At last, in Section 5, we summarize and conclude about our method.

II Methodology

II-A Latent Gaussian Models

A latent Gaussian Models (LGM) is a class of statistical models, which include many commonly used statistical models such as auto-regressive (AR) model, conditional auto-regressive (CAR) model, and generalized linear mixture (GLM) model.

Fig.1 shows the outline of an LGM. In LGMs, latent variables $\boldsymbol{x}=(\boldsymbol{x}_{1},\boldsymbol{x}_{2},\dots,\boldsymbol{x}_{n})^{\top}$ follow Gaussian Markov random field (GMRF)

[TABLE]

where $\boldsymbol{\theta}_{\mathrm{lat}}$ is a set of hyperparameters of latent GMRF. In addition, mean vector $\boldsymbol{\mu}$ and precision matrix $\boldsymbol{\Sigma}^{-1}$ are parameterized by $\boldsymbol{\theta}_{\mathrm{lat}}$.

Observations $\boldsymbol{y}=(\boldsymbol{y}_{1},\boldsymbol{y}_{2},\dots,\boldsymbol{y}_{n})^{\top}$ should be conditionally independent and identically distributed on a parameter set $\boldsymbol{\theta}_{\mathrm{obs}}$

[TABLE]

and conditional distribution $p(\boldsymbol{y}_{i}|\boldsymbol{x},\boldsymbol{\theta}_{\mathrm{obs}})$ must be belong to exponential family. Hereafter, for simplicity, we will put hyperparameter sets into $\boldsymbol{\theta}$,

[TABLE]

II-B Integrated Nested Laplace Approximation

Integrated Nested Laplace Approximation (INLA) is a fast and accurate method to approximate posterior distributions of LGMs[6]. Especially when the number of hyperparamters is small (empirically $|\boldsymbol{\theta}|\lesssim 5$[7]), INLA performs very faster than Markov chain Monte Carlo (MCMC).

In LGMs, our aim is to get marginalized posteriors of hyperparameters

[TABLE]

and latent variables

[TABLE]

In equation (4), each term in the numerator can be calculated by forward computation, but the denominator can not. Hence, first we apply Laplace (or simplified Laplace, which can fit with a skew normal distribution[13, 6]) approximation to the denominator

[TABLE]

where $\tilde{p}(\cdot)$ denotes an approximation of $p(\cdot)$.

For latent variables (5), applying numerical integration as following:

Exploring the mode of $\mathrm{log}~{}\tilde{p}(\boldsymbol{\theta}|\boldsymbol{y})$ using an optimization algorithm (e.g., quasi-Newton method) with respect to $\boldsymbol{\theta}$. 2. 2.

Arranging integration points in accordance with grid or central composite design (CCD)[14] strategy (see Fig.2). 3. 3.

Calculating $\mathrm{log}~{}\tilde{p}(\boldsymbol{\theta}_{h}|\boldsymbol{y})$ at arranged points, where $h\in\{1,\ldots,H\}$ is the index of each point. 4. 4.

Approximating equation (5) with weighted summation on $\theta_{h}$:

[TABLE]

where $\Delta_{h}$ is the weight of $\boldsymbol{\theta}_{h}$.

III Our Model

III-A ICAR Models

To use INLA, we should choose a latent structure carefully. This is because the computaional performance of INLA strongly depends on the number of hyperparameters of the LGM.

ICAR model is one of the commonly used representations of spatial interaction[15]. In ICAR model, $x_{i}$ has interaction with $x_{j}~{}(j\in C_{i})$, where $C_{i}$ indicates the set of indices of nearest neighbor nodes of $x_{i}$ (for simplicity, each variable is assumed to be a scalar value). Here it assumed that $x_{i}$ follows

[TABLE]

where $\sigma^{2}$ is the hyperparameter indicates variance between nodes. This model means some assumptions on $\boldsymbol{x}$:

• •

Each node hardly differs from adjacent ones.

• •

The similarity between adjacent nodes is controlled by the fixed effect $\sigma^{2}$.

• •

As $|C_{i}|$ is increasing, the variance of $x_{i}$ decreases.

At a glance, this model seems not to be inappropriate. However, the precision matrix of $\boldsymbol{x}$ does not satisfy the condition of positive definiteness, then ICAR models become improper. Therefore, we need to regularize them to guarantee its computational stabiliy.

III-B Proper ICAR Models

To guarantee properness of ICAR models, two reguralization ways are frequently implemented: scaling $\boldsymbol{x}$ by $\sum_{i}x_{i}=0$ or adding extra terms to the parameters of equation (8). We adopt the latter strategy, since the number of detected photons must be non-negative. We use following proper ICAR model for the latent GMRF of adopted model,

[TABLE]

where $d>0$ is the extra hyperparameter for stability. In this study, we configure interaction in a grid without periodic boundary condition for the latent GMRF to deal with images. For instance, we consider the probability of a non-corner pixel $p(x_{l})$ of an $L\times L$ image. Here, we can deform equation (9) into

[TABLE]

where

[TABLE]

III-C Observation Model

In photon-limited images, the number of observated photons $y_{i}$ fluctuates as following Poisson distribution with the true value $x_{i}$. Hence, we use

[TABLE]

as the observation model straightforwardly, and $y_{i}$ is assumed to be independent from other observations with given $\boldsymbol{x}$. Therefore, the joint distribution of $\boldsymbol{y}$ becomes

[TABLE]

III-D Hyperpriors

In our model described above, we have two hyperparameters $\sigma^{2}$ and $d$. As their hyperpriors, we use uninformative prior on $[0,\infty)$.

IV Experiments

IV-A Preprocessing

To handle general images with our model and control their contrast, we firstly transform the pixel value $I_{i}$ in the original image into positive value $x_{i}$, where the original image $\boldsymbol{I}=(I_{1},\ldots,I_{n})^{\top}$ and its transformation $\boldsymbol{x}$ consist of $n$ pixels. In particular, we give linear transformation

[TABLE]

where $\lambda_{\mathrm{max}}$ and $\lambda_{\mathrm{min}}$ are given parameters indicate maximum and minimum values of $\boldsymbol{x}$ respectively, then $I_{\mathrm{max}}$ and $I_{\mathrm{min}}$ are the maximum and minimum values of $\boldsymbol{I}$. We can control the range of intensity by varying $\lambda_{\mathrm{max}}$ and $\lambda_{\mathrm{min}}$.

IV-B Evaluation Criteria

In our simulations, we use peak signal-to-noise ratio (PSNR) and structural similarity (SSIM)[16] to measure similarities between original and denoised images. PSNR is defined by

[TABLE]

where $\boldsymbol{G}=(G_{1},\ldots,G_{n})^{\top}$ and $\boldsymbol{H}=(H_{1},\ldots,H_{n})^{\top}$ denote images have $n$ pixels. Here, PSNR just focuses on the each corresponding pixel value. To consider the perceptual similarity between two images, we also introduce SSIM for the evaluation criterion. SSIM is defined as follows:

[TABLE]

where $\mu_{G}$ and $\mu_{H}$ are the means of $\boldsymbol{G}$ and $\boldsymbol{H}$ respectively, and $\sigma^{2}_{G}$ and $\sigma^{2}_{H}$ are the variances of $\boldsymbol{G}$ and $\boldsymbol{H}$ respectively. Also, $\sigma_{G,H}$ indicates the covariance between $\boldsymbol{G}$ and $\boldsymbol{H}$. In addition, $c_{1},c_{2}$ are small positive real numbers for the stability.

Here, PSNR can take from [math] to infinity, and SSIM is in range of $-1$ to $1$. In each criterion, the larger value is, the closer two images are.

IV-C Simulations

In this study, we compared denoising simulation results with INLA, MCMC (Hamiltonian Monte Carlo with No U-Turn Sampler[17, 18]), and LBP using the described model in Section 3. We simulated MCMC for 2000 steps and set first 1000 steps as burn-in. As the original images, we use “Lena”, “Boat”, and “Cameraman” for the evaluation (in Fig.3). Here, we fixed the contrast parameters as $\lambda_{\mathrm{min}}=2$ and $\lambda_{\mathrm{max}}=25$, and ran the simulations on the computer has AMD Ryzen Threadripper 2990WX (32 cores, 3.0 GHz) as its CPU and 128 GB memory. To compute point estimates of $\boldsymbol{x}$, we used expected a posteriori (EAP) estimation in INLA and MCMC.

Table I shows the comparison of results of INLA, MCMC, and LBP. In Table I, we can find that INLA is very faster than other algorithms with little decreasing PSNR and SSIM from MCMC, expect for “Cameraman”. LBP is too slow because it requires spectral decomposition of an $n\times n$ matrix in its algorithm. Fig.4 is the denoised images. Qualitatively, the images restored by INLA look clearer and unhazier than by other methods on a whole. On the other hand, we could not estimate the original “Cameraman” appropriately using INLA (Fig.4(c)). Specifically, the result of INLA is apparently different from others in its brightness. Then, to give additional consideration, we present histograms of intensities in Fig.5. Seeing Fig.5(b), we can find that brighter pixels have been hardly restored in “Cameraman”. In “Lena”, the shape of true histogram is slightly flat. In contrast, the true histogram of “Cameraman” has two sharp peaks. Therefore, it is conceivable that the error of INLA mainly caused by this illness. In this case, INLA might fail to estimate the hyperparameters due to the illness of the original image, then also failed to get appropriate posteriors with respect to the latent variables.

Additionaly we indicate evaluated posteriors or its point estimates of some pixels of each image in Fig.6. The differences of posteriors between INLA and MCMC are unnoticeable in “Lena” (Fig.6(a)). In “Boat”, the pixel which is near a sharp edge is slightly different the histogram by MCMC. Moreover, we can also find that brighter pixels could not be restored in “Cameraman” by seeing Fig.6(c).

V Conclusion

In this study, we proposed and evaluated a denoising methods for Poisson corrupted images based on INLA. INLA requires an assumption of GMRF to latent variables. Hence, its accuracy decreases when the assumption is inappropriate. Whereas when the assumption is suitable for the original image, we can reduce much time compared to other Bayesian computational algorithms, such as LBP or MCMC with enough accuracy.

As a future work, we can also consider to use other models for latent models. For example, to apply segmentation to noisy image firstly, and then evaluate posteriors of each segmentation independently. In addition, we should adopt our method to real X-ray imaging data, and then evaluate its effectivity.