Abstract: In this work, we present a new and highly efficient Markov chain Monte Carlo (MCMC) methodology to perform Bayesian inference in low-photon imaging problems, with particular attention given to situations involving observation noise processes that deviate significantly from Gaussian noise, such as binomial, geometric, and low-intensity Poisson noise. These problems are challenging for many reasons. From an inferential viewpoint, low-photon numbers lead to severe identifiability issues, poor stability, and high uncertainty about the solution. Moreover, low-photon models often exhibit poor regularity properties that make efficient Bayesian computation difficult, e.g., hard nonnegativity constraints, nonsmooth priors, and log-likelihood terms with exploding gradients.

We address these difficulties by proposing an MCMC methodology based on a reflected and regularized Langevin SDE, which is shown to be well-posed and exponentially ergodic under mild and easily verifiable conditions. For more information, you can find our published paper here https://epubs.siam.org/doi/10.1137/22M1502240.

In this repository, you can find Pytorch codes of our suggested reflected proximal Langevin MCMC algorithms to perform Bayesian computation in low-photon imaging problems. The proposed approach is demonstrated with a range of experiments related to image deblurring, denoising, and inpainting under binomial, geometric, and Poisson noise.

Our repository is devided according to the experiments conducted in the paper under different settings:

├── Binomial cases # all experiments with Binomial noise
│   ├── Binomial denoising_MIV_0.1_t_100 
│   ├── Binomial_denoising_MIV_1_t_10
├── Geometric_cases # all experiments with Geometric noise
│   ├── Geometric_inpainting_MIV_0.01
│   ├── Geometric_inpainting_MIV_0.1
├── Poisson_cases # all experiments with Poisson noise
│   ├── Poisson_deblurring_MIV_1

* MIV stands for the Mean Intensity Value. Before we generate the synthetic data, we scale the true image (signal) to this value. That's because the aforementioned noise processes are signal-dependant. This means that so to generate the desired level of noise we need to scale the true signal accordingly.

In each of the cases, you will find codes for 4 different sampling algorithms which can be run for example as :

python3 Poisson_deblurring_SKROCK.py
python3 Poisson_deblurring_MYULA.py
python3 Poisson_deblurring_MYUULA.py
python3 Poisson_deblurring_MYMALA.py

It should be noted that by far the faster algorithm is R-SKROCK.

Below, we present an experiment under severe Poisson noise which leads to low photon count in the data. The first reconstruction is the Minimum Mean Squared Error (MMSE) estimator calculated by our suggested algorithm R-SKROCK. The second reconstuction is the Maximum-a-Posteriori (MAP) estimator calculated by PIDAL [1], an ADMM algorithm adapted for Poisson noise. We used TV convex regularization for these recontructions (however data-driven convex regularizations can also be considered).

Ground truth	Noisy image (Poisson)	MMSE (PSNR 20.53 dB)	MAP (PSNR 19.21 dB)

The use of MCMC methods allows for uncertainty quantification (not usually possible when optimization methods or end-to-end networks are used). Below, we present uncertainty visualisation plots presenting the marginal standard deviation of pixels at different resolutions. See our paper for more information.

Below, we illustrate the convergence speed of the different MCMC methods considered in our paper. More precisely, we plot the evolution of the NRMSE estimation for the MMSE solutions as a function of the number of gradient evaluations (in log-scale). We see that R-SKROCK has the highest convergence speed followed by SPA [2].

To reproduce the results of R-SKROCK, after cloning the repository, just navigate to the Poisson_cases/Poisson_deblurring_MIV_1 folder and run:

python3 Poisson_deblurring_SKROCK.py

Below, we present an experiment under the very difficult case of Geometric noise in the data under a masking operator. The first reconstruction is the MMSE calculated by our suggested algorithm R-SKROCK. To the best of our knowledge, there is no other sampler in the literature that can deal with such noise process. The second reconstuction is the MAP. To calculate the MAP, we modified PIDAL [1] to work under the Geometric likelihood. We used TV convex regularization for these recontructions (however data-driven convex regularizations can also be considered).

Ground truth	Noisy image (Poisson)	MMSE (PSNR 19.6 dB)	MAP (PSNR 19.7 dB)

Below, we present uncertainty visualisation plots presenting the marginal standard deviation of pixels at different resolutions. See our paper for more information.

To reproduce the results of R-SKROCK, after cloning the repository, just navigate to the Geometric_cases/Geometric_inpainting_MIV_0.01 folder and run:

python3 Geometric_inpainting_SKROCK.py

[1] M. A. T. Figueiredo and J. M. Bioucas-Dias, Restoration of Poissonian Images using Alternating Direction Optimization, IEEE Transactions on Image Processing, 19 (2010), pp. 3133–3145.

[2] M. Vono, N. Dobigeon, and P. Chainais, Bayesian image restoration under Poisson noise and log-concave prior, in ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2019, pp. 1712–1716.