Bayesian Image Denoising with Multiple Noisy Images

In this paper, we propose a fast image denoising method based on discrete Markov random fields and the fast Fourier transform. The purpose of the image denoising is to infer the original noiseless image from a noise corrupted image. We consider the case where several noisy images are available for inferring the original image and the Bayesian approach is adopted to create the posterior probability distribution of the denoised image. In the proposed method, the estimation of the denoised image is achieved using belief propagation and an expectation–maximization algorithm. We numerically verified the performance of the proposed method using several standard images.


Introduction
Bayesian image processing based on Markov random fields (MRFs) is an important framework in the field of image processing [1,2]. An MRF is a undirected graph representation of probability distribution, and many applications of MRFs exist in the image processing and computer vision fields [3][4][5]. MRFs have also been applied to other research fields, including traffic engineering [6,7] and earth science [8,9]. In Bayesian image processing, the objective image can be inferred based on the posterior probability distribution.
Recently, we proposed a fast image denoising method for the case where multiple noisy images are available for inferring the original noiseless image that is based on 1 3 Gaussian MRFs [10]. However, in the study in [10], for ease of mathematical treatment, we made an unnatural assumption that pixel values are continuous. In general, a pixel takes a discrete value from 0 to 255, and an additional framework is required to treat pixel values as discrete instead of continuous values. Therefore, in this paper, we focus on the Bayesian image denoising problem of inferring the original noiseless image from multiple noisy images when the pixel values are treated as discrete values. We created a probability model for image denoising based on the discrete MRF and Bayesian perspective. A major disadvantage of an image processing model based on discrete MRFs is the computational complexity. In fact, the inference problem from discrete MRFs belongs to the NP-hard class. Therefore, an approximate inference technique is required to infer the objective image from a discrete MRF. Belief propagation [11] is known as one such effective technique. In this paper, we propose an effective image denoising algorithm for multiple noisy images that applies belief propagation. The main contributions of this paper are that an MRF model for image denoising with multiple noisy images is defined and a fast effective denoising algorithm based on our discrete MRF model and the fast Fourier Transform (FFT) is proposed.
The remainder of this paper is organized as follows. In Sect. 2, we define a probability model for image denoising with multiple noisy images based on the discrete MRF and Baysian perspective. In Sect. 3, we derive an image denoising algorithm based on the posterior probability distribution defined in Sect. 2. In Sect. 4, we describe the framework for estimating the parameters in the posterior probability distribution. We explain the implementation of our denoising method using the FFT in Sect. 5. In Sect. 6, we describe the numerical verification of the performance of the proposed method. Finally, in Sect. 7, we present our concluding remarks.

Framework of Bayesian Image Denoising Method
In this section, we briefly explain the framework of the Bayesian image denoising method for the case where multiple noisy images are available. Suppose that we have K degraded images that are independently obtained by adding additive white Gaussian noise (AWGN) to the original image. We assume that the images are composed of sional vectors corresponding to the original image and the k-th noisy image, respectively. Vectors x and y (k) can be easily obtained by raster scanning the images. We assume that The purpose of the image denoising is to infer the original noiseless image x from K noisy images Y = y (1) , y (2) , … , y (K) . In the Bayesian framework, the original image x can be inferred using the posterior probability distribution P(x|Y) that is expressed as where ∑ x denotes the multiple summations over all the possible L N states of x . The framework of the proposed Bayesian image denoising method is illustrated in Fig. 1. From the definition of y (k) , the probability density function P(Y|x) is expressed as

3
The Review of Socionetwork Strategies (2019) 13:267-280 where V = {1, 2, … , N} and 2 is the variance of the AWGN. We express the parameters of the probability model by its arguments after the semicolon as Eq. (2). We define the prior probability distribution as where E is a set of edges of the h × w lattice graph and (x) is a downward convex even function taking its minimum at x = 0 . In this study, we assumed the periodic boundary condition on the graph structure E, as demonstrated in Fig. 2. > 0 is the parameter of the prior probability distribution; if is set to a large value, neighboring x i and x j tend to take close values. Z prior ( ) is a normalization constant defined as By substituting Eqs. (2) and (3) into Eq. (1), the posterior probability distribution P(x|Y) is expressed as

Inference Algorithm Based on Belief Propagation
The image denoising is achieved by finding the image x that maximizes the posterior probability distribution in Eq. (5): However, the problem of determining such an image is intractable, because this maximization problem belongs to the NP-hard class. Therefore, we need an approximate inference method to find x . In this section, we explain an effective approximate inference method called belief propagation for inferring the denoised image x. Belief propagation is a method of computing the approximate marginal distributions b i x i and b ij x i , x j for each i ∈ V and ij ∈ E . In the belief propagation framework, the approximate marginal distributions b i x i and b ij x i , x j are given by where Z j→i is a normalization constant to ensure algorithmic stability. The estimation of the denoised image x = x i |i ∈ V is achieved by finding for i ∈ V in the belief propagation framework.

Parameter Estimation Using Expectation-Maximization Algorithm
In the preceding section, we explained the method for inferring the denoised image x based on belief propagation. In our framework, the denoised image is inferred from the posterior probability distribution in Eq. (5), which has two parameters, and 2 . It is obvious that the inferred denoised image x depends on these parameters. In this section, we explain the method for determining these parameters from degraded images Y based on the expectation-maximization (EM) algorithm [12].
The EM algorithm is a statistical inference method to infer the maximum likelihood estimates by an iterated method. In the EM algorithm framework, the parameters and 2 are estimated by iterative maximization of the Q function defined as where t and 2 t are estimates of the parameters at the t-th iteration and Using belief propagation, we can approximate the expectations in Eq. (17) as and respectively, where b (t) i x i and b (t) ij x i , x j are the approximate marginal distributions of the posterior probability distribution P x|Y; t , 2 t computed using Eqs. (9) and (11). The parameter update at iteration t is given by and the maximum likelihood estimates in Eq. (16) are given as the convergence point of the above iterative estimation. By differentiating the Q function with respect to and 2 and considering the conditions for the extremal value, the updated parameter ; this expectation can also be computed approximately using belief propagation similarly to Eq. (20). Using the bisection method, we can easily find t+1 that satisfies Eq. (22).

Proposed Algorithm: Fast Implementation Based on the Fast Fourier Transform
The image denoising algorithm based on our probabilistic model in Eq. (5)  It should be noted that the computation time of the message update can be reduced to O(L log L) using the FFT [13]. It has been confirmed that this FFT-based method in fact accelerates the message computation for the probabilistic image denoising model in Eq. (5) for the case where K = 1 [14]. However, there exist 1 3 additional O(L 2 ) computation terms in Algorithm 1: Steps 11 and 12. In this section, we show that these O(L 2 ) computation terms can also be computed in O(L log L) using the FFT. Therefore, we can reduce the worst computation time of Algorithm 1 to O T EM T BP NL log L .
The key idea for accelerating the message updates in Eq. (13)    In Eq. (22), we need to calculate the messages and expectation of prior probability distribution P x; t to find t+1 that satisfies this equation using the bisection method. It should be noted that, because we assume the periodic boundary condition for the graph structure E, we can calculate the messages and expectation of prior probability distribution faster than those of posterior probability distribution by considering the translational symmetry assumption. If we assume both a periodic boundary condition and translational symmetry, the messages and expectation of prior probability distribution become not dependent on the position of the edges ij ∈ E . Therefore, the message update rule and expectation calculation for prior probability distribution are expressed as and respectively, where M prior x i is a message of the prior probability distribution. The calculation of the message and the expectation of the prior probability distribution in Eqs.

3 6 Numerical Experiments
In this section, we describe the numerical verification of the proposed method. We used the standard images in Fig. 3, which are widely used in the image processing research field. The pixel values of these images take L = 256 different values. All the experiments were implemented using C++ and were run single-threaded on an Ubuntu 18.04.1 LTS (64 bit) machine with an Intel Core i7-6850K CPU running at 3.60 GHz and 128 GB RAM. In this experiment, we defined the function (x) as and set the parameters of the proposed method as follows. The initial parameter 0 was set at 0.005 and 2 0 was set at the sample variance calculated from Y. The maximum numbers of iterations T EM and T BP were set at 100 and 1000, respectively. We considered that the messages converged if the absolute value of the average change in the messages was smaller than 10 −4 ; the same applied to the parameters and 2 . We set the search interval of the bisection method for computing t+1 as [0, 2 t ].
First, we compared the computation time of the proposed method with that of the previous methods (belief propagation FFT (BP-FFT) and Naive). The Naive method is the naive implementation version of Algorithm 1; the worst computation time is O T EM T BP NL 2 . BP-FFT is the method used in the study presented in [14], where only the message updates in Eqs. (13) and (14) were speeded-up by using the FFT. Tables 1  and 2 show the average computation times over 10 trials for each method where K noisy images Y were generated by adding an AWGN of = 15 to the original noiseless images. According to the results, the proposed method was faster than the other methods. It should be noted that the difference between the three methods is in whether Algorithm 1 is implemented using the FFT. Therefore, the image denoising results of these method are all the same.     Proposed ( Figure 4 shows the average computation time versus image size for K = 1 and K = 5 over 10 trials, where the noise level of AWGN was = 15 . The computation time of our denoising methods grows approximately linearly with the increase in the image size (it dose not grow strictly linearly, because we break off the message and parameter updates according to the convergence condition). Figure 5 shows the denoising performance of the proposed method versus various values of K for two levels of AWGN ( = 15 and = 30 ). We evaluated the performance of the method according to the average peak signal-to-noise ratio (PSNR) over 10 trials. The PSNR is defined as (40) PSNR = 10 log 10 255 2 MSE ,  Fig. 3a. b Fig. 3b. c Fig. 3c. d Fig. 3d 1 3 where MSE is the mean squared error between the original noiseless image and the inferred denoised image x . Figure 5 conforms that the image denoising results improve as the value of K is increased. Example of image denoising results for Fig. 3a are shown in Fig. 6 for K = 1 and K = 5 , respectively.

Concluding Remarks
In this paper, we defined a discrete MRF model for the Bayesian image denoising problem with multiple noisy images. We proposed a fast denoising algorithm for inferring a denoised image in O T EM T BP NL log L -time by using belief propagation and an EM algorithm based on our MRF model and FFT. We numerically verified the proposed denoising method using standard images. The results show that the proposed algorithm inferred the denoised image faster than previous implementation methods that use belief propagation. We believe that the proposed method is the most fastest implementation of an image denoising algorithm based on a discrete MRF model that uses belief propagation and an EM algorithm. However, the method cannot yet be used for real-time processing. Therefore, we need to seek a further effective fast approximate method that preserves the restoration quality for the discrete MRF model. In our experiment, we adopted the quadratic function as the form of function (x) . However, the proposed method is not restricted to the quadratic function: we can apply it to other types of the function (x) , such as the Huber prior [15] and generalized sparse prior [16]. Moreover, because it can be used in any discrete MRF with (x i − x j ) potential for ij ∈ E interaction, it is expected that the proposed method is applicable to not only image denoising but also other inference problems such as sparse modeling [17,18]. We intend to develop the method in these directions.