Restoring Poissonian Images by a Combined First-Order and Second-Order Variation Approach

The restoration of blurred images corrupted by Poisson noise is an important topic in imaging science. The problem has recently received considerable attention in recent years. In this paper, we propose a combined first-order and second-order variation model to restore blurred images corrupted by Poisson noise. Our model can substantially reduce the staircase effect, while preserving edges in the restored images, since it combines advantages of the first-order and second-order total variation. We study the issues of existence and uniqueness of a minimizer for this variational model. Moreover, we employ a gradient descent method to solve the associated Euler-Lagrange equation. Numerical results demonstrate the validity and efficiency of the proposed method for Poisson noise removal problem.


Introduction
Image restoration problem has been widely studied in the areas of image processing.The goal of image restoration is to reconstruct an approximation of an original image from a blurred and noisy one [1][2][3][4].Over the last decades, most of the literature deal with the additive noise and multiplicative noise model [5][6][7][8][9][10][11].Given a blurred and noisy image  =  + V or  = V, where  is the original image, V represents the noise, and the blurring operator  is a point spread function (PSF); the restoration problem involving the additive noise or multiplicative noise is to recover  from the observed image .The additive noise and multiplicative noise models have been extensively studied.We refer the reader to [12][13][14][15][16][17][18][19] and references therein for a review of various methods.In fact, there are other types of noise such as Poisson noise.In this paper, we consider the problem of seeking the approximations of original images from blurred images corrupted by Poisson noise.The restoration of blurred images corrupted by Poisson noise is an important task in various applications such as astronomical imaging, electronic microscopy, single-particle emission-computed tomography (SPECT), and positron emission tomography (PET) [20][21][22][23][24].
In a Poisson noise model, the recorded image  is a realization of a random variable ĝ.A statistical model for ĝ is given by ĝ ∼ Poiss ( + ) , where  > 0 is the expected value of the background which is assumed to be constant.The given data of the problem are the imaging operator , the expected value of the background , and the detected image .It is known that the major difficulty for finding the original image  of (1) is the models and algorithms developed for deconvolution of additive noise, or multiplicative noise images cannot be directly applied to the restoration of Poissonian images.Therefore, it is important to devise efficient and reliable algorithms for restoring blurred images corrupted by Poisson noise.The problem of restoration of Poissonian images has received considerable attention in recent years.We note that, in a discrete setting, the intensity   in the pixel  is a random variable that follows a Poisson law of mean ( + )  where the structured matrix  related to the boundary conditions is called a blurring matrix.Considering The use of Stirlings formula leads to the well-known Csiszár [25] -divergence measure between  +  and .Indeed, by taking the negative logarithm of the likelihood function, we obtain T 0 () = − log [ ()] = ∑  ( + )  −   +   log   ( + )  . ( This measure generalizes the Kullback-Leibler divergence or cross-entropy measure to accommodate functions whose integrals are not constant, as they would be if they were probability distributions.Dropping the terms independent of  in L(), we obtain Therefore, the maximum likelihood estimator of the original image is the minimizer with respect to  of the negative-log Poisson likelihood functional J 0 ().
As is well known, the structured matrix  has many singular values of different orders of magnitude close to the origin.This makes the minimizer of (4) very sensitive to the noise in the right-hand side .Thus a regularization method may be employed to compute the approximate solution that is less sensitive to noise than the naive solution.In general, regularization methods formulate the image restoration problem as a minimization problem of the form: where J  () is a regularization function and  > 0 is the regularization parameter that controls the balance between J 0 and J  .We note that the regularization parameter  is an important quantity which controls the properties of the regularized solution, and  should therefore be chosen with care.Throughout the years a variety of parameter choice strategies such as the discrepancy principle, the -curve, and generalized cross-validation (GCV) have been proposed for Poisson noise removal problem [26][27][28].
It is shown in [29] that the choice of the regularization function J  () is related to the features of the images to be restored.The use of different regularization functions is an active research area.Probably one of the most popular regularization functions is the Tikhonov regularization function J  () = (1/2)‖‖ 2  2 , where  is an auxiliary operator chosen among the identity and low-order differential operators; see [30] for more details.In [31], Landi and Piccolomini proposed a quasi-Newton projection method for deblurring Poisson-corrupted images after formulating the image restoration problem as a nonnegatively constrained minimization problem where the objective function is the sum of the Kullback-Leibler divergence, used to express fidelity to the data in the presence of Poisson noise, and of a Tikhonov regularization term.Their numerical results show the potential of the method both in terms of relative error reduction and computational efficiency.In [32,33], the authors developed a cost functional which incorporates the statistics of the noise in the image data and Tikhonov regularization to induce stability and employed an efficient hybrid gradient projection-reduced Newton (GPRNCG) method for the problem of restoration of Poissonian images.Finally, they gave a comparison between the Richardson-Lucy [34,35] iterative regularization method and GPRNCG algorithm for the restoration of blurred images corrupted by Poisson noise.By the measures of performance used in this comparison, GPRNCG was more efficient than Richardson-Lucy iteration.
It is known that Tikhonov regularization estimate is similar to low-pass filtering; therefore, it produces a smoothing effect on the restored images; that is, it penalizes edges, which is not a good approximation of the original image if it contains edges.To overcome this shortcoming, Rudin et al. [5] proposed a total variation-(TV-) based regularization technique, which preserves the edge information in the restored image.In the case of TV regularization, the estimated solution is obtained by minimizing the object function: where ‖‖ TV = ‖∇‖ 1 is the discrete TV regularization term.It is shown in [36] that the TV penalty function acts more like a model selection device by identifying a small number of critical points at which  is allowed to jump.The number of these selected jump points is controlled by .
Especially, if  is piecewise constant with a finite number of jump discontinuities, then ‖‖ TV gives the sum of magnitudes of the jumps.
In [21], Bardsley and Luttman considered the TV regularization for poissonian images and showed that the problem of computing the minimizer of the resulting TV-penalized Poisson likelihood functional is well posed.They then proved that as the errors in the data and in the blurring operator tend to zero, the resulting minimizer converges to the minimizer of the exact likelihood function.Finally, the practical effectiveness of the approach is demonstrated on synthetically generated data, and a nonnegatively constrained, projected quasi-Newton method was introduced.In [22], Figueiredo and Bioucas-Dias proposed an approach to deconvolving Poissonian images with the TV function, which is based upon an alternating direction optimization method.Using standard convex analysis tools, they presented sufficient conditions for existence and uniqueness of solutions of these optimization problems.In [23], Setzer et al. considered the restoration of blurred images corrupted by Poisson noise by minimizing an energy functional consisting of the Kullback-Leibler divergence as similarity term and the TV regularization term.Their minimizing algorithm uses alternating split Bregman techniques which can be reinterpreted as Douglas-Rachford splitting applied to the dual problem.In contrast to recently developed iterative algorithms, their algorithm contains no inner iterations and produces nonnegative images.The high efficiency of their algorithm in comparison to other recently developed algorithms is demonstrated by artificial and real-world numerical examples.We refer the reader to [24,37] and references therein for other efficient methods solving the restoration of blurred images corrupted by Poisson noise with the TV function.
Although the total variation regularization is extremely popular in a variety of applications, it also gives rise to some undesired effects.Both from a theoretical and experimental point of view, it has been shown that the TV norm transforms smooth signal into piecewise constants, the socalled staircase effect.Staircase solutions fail to satisfy the evaluation of visual quality, and they can develop false edges that do not exist in the true image.To attenuate the staircase effect, there is a growing interest in the literature for replacing the TV norm by a high-order TV norm.The motivation behind such attempt is to restore potentially a wider range of images, which comprise more than merely piecewise constant regions.The majority of the high-order norms involve second-order differential operators because piecewise-vanishing second-order derivatives lead to piecewise-linear solutions that better fit smooth intensity changes; see [38] for more details.
Second-order regularization schemes have been considered so far in the literature mainly for dealing with the staircase effect while preserving the edge information in the restored image.There are two classes of second-order regularization methods for image restoration problems.The first class combines a second-order regularizer with the TV norm.For example, a technique in [39,40] combining the TV filter with a fourth-order PDE filter was proposed to preserve edges and at the same time avoid the staircase effect in smooth regions for noise removal.In [41], Papafitsoros and Schönlieb considered a high-order model involving convex functions of the TV and the TV of the first derivatives for image restoration problems and used the split Bregman method to solve numerically the corresponding discretized problem.Chan et al. proposed a second-order model to substantially reduce the staircase effect, while preserving sharp jump discontinuities for image restoration in [42].The second class employs a second-order regularizer in a standalone way.The secondorder regularizer was first proposed in [43] for additive noise removal in medical magnetic resonance images.The LLT model presented in [43] can overcome the staircase effect that occurs with the TV norm filter and better preserve the fine details in nonblocky images.In order to test its practical potential the authors have applied their method to a wide range of real images, including structural and functional MRI data sets.The main strength of their method was the ability to process signals with a smooth change in the intensity value.Numerical results have shown that, compared with some related partial differential equation (PDE) models, the LLT model is rather robust in removing noise and handling edges.In [44,45], the authors reconsidered the fourth-order PDE model for additive noise removal and employed the dual algorithm of Chambolle for solving the high-order problems.
In order to get restored images with edge preserved and simultaneously smooth region reconstructed for Poisson noise removal, it is natural to utilize a combined first-order and second-order total variation technique.In this paper, we consider to modify the total variation model by adding a high-order functional for restoring blurred images corrupted by Poisson noise.Our model can substantially reduce the staircase effect, while preserving edges in the restored images, since it combines advantages of the first-order and second-order total variation.We study the issues of existence and uniqueness of a minimizer for this variational model.Moreover, we employ a gradient descent method to solve the associated Euler-Lagrange equation.Some numerical results demonstrate the validity and efficiency of the proposed method for Poisson noise removal problem.
The organization of this paper is outlined as follows.In the next section, we present a combined first-order and second-order variation model to restore blurred images corrupted by Poisson noise and study the issues of existence and uniqueness of a minimizer for this variational model.A gradient descent method to solve the associated Euler-Lagrange equation is presented in Section 3. Some numerical experiments are given to illustrate the performance of the proposed algorithm in Section 4. We give concluding remarks in the last section.

A Combined First-Order and Second-Order Variation Model
In the Bayesian approach, a prior probability density   () for  is also specified, and the posterior density: given by Bayes's law, is maximized with respect to .In (7),   () defines the prior probability distribution of the measurements , which are given by measurements, and the likelihood function   ( | ) is the conditional probability of observing a fixed  for a variable .The maximizer of   ( | ) is called the maximum a posteriori (MAP) estimate.
Since   () does not depend on , it is not needed in the computation of the MAP estimate and thus can be ignored.
Taking the natural logarithm of (7) and dropping the terms independent of , we obtain that maximizing ( 7) is equivalent to minimizing with respect to , where J 0 () = ∫ Ω (( + ) −  log( + )).It is clear that − log   () corresponds to the regularization term in the classical penalized likelihood approach to regularization.However in the Bayesian setting,   () is the probability density known as the prior from which the unknown  is assumed to arise.Thus the prior knowledge regarding the characteristics of  can be formulated in the form of a probability density   (), and this yields to a natural, and statistically rigorous, motivation for the regularization method [27].Standard Tikhonov regularization corresponds to the following choice of the prior: This corresponds to the assumption that the prior for  is a zero-mean Gaussian random variable with covariance matrix  −1 , which has the effect of penalizing reconstructions with large  2 norm.For  2 norm of the gradient regularization, the penalty has the similar form where ∇ is the gradient operator.The use of this regularization function has the effect of penalizing reconstructions that are not smooth.As an alternative to Tikhonov regularization including standard Tikhonov regularization and gradient regularization for Poisson noise removal problem, total variation regularization is another regularization technique that allows for the presence of sharp edges in the resulting reconstruction.In the case of total variation regularization, we have where ‖∇‖ 1 = ∫ Ω √ 2  +  2  .Optimization techniques such as alternation direction method of multipliers (ADMM) and alternating split Bregman algorithm were proposed for the restoration of blurred images corrupted by Poisson noise by minimizing an energy functional consisting of the Kullback-Leibler divergence as the similarity term and the TV regularization term; see [21][22][23][24] for more details.In the other hand, some Newton-related gradient methods were proposed for deblurring Poissoncorrupted images by formulating the image restoration problem as a nonnegatively constrained minimization problem where the objective function is the sum of the Kullback-Leibler divergence, used to express fidelity to the data in the presence of Poisson noise, and of a Tikhonov regularization term; see [31][32][33] for more details.
To attenuate the staircase effect in the TV regularization, we consider a combined first-order and second-order variation model to restore blurred images corrupted by Poisson noise in this work.In the proposed model, the prior probability density is of the form: where where the parameter  ∈ [0,1] is used to control the balance between the edges and the smooth surface and  is the regularization parameter which measures the tradeoff between the fidelity term and the regularized term.The parameter  can be computed by using the information of the edges and smooth regions of the resulting image obtained by smoothing the observed image  with the low-pass filters such as the median filter and the Gauss filter.The functional () in ( 13) is defined on the set of  ∈ (Ω) ∩  2 (Ω) such that log  ∈  1 (Ω); in particular,  must be positive almost everywhere.Some basic notations and properties concerning the  space and  2 space can be found in [39,40].Motivated by [46], we have the following result to show the existence and uniqueness of the minimizer for the model (13).
Theorem 1.Let Ω be a bounded, open subset of R 2 with Lipschitz boundary.Assume that  is a positive bounded function and  is positive definite.Then () for  ∈ (Ω) ∩  2 (Ω) such that log  ∈  1 (Ω) has a unique minimizer for the model (13).
Since  is positive definite and  is positive, it immediately follows that J 0 () is strictly convex.Obviously, J  () is a convex function.We conclude that () is strictly convex, as the sum of a convex function and of a strictly convex function.Therefore, the minimizer  * is unique.

Computational Method
In this section, our aim is to propose a time-marching gradient descent method for computing the minimizer of the model (13).From [5,39,40] we know that minimizing (13) for a given constant  yields the associate Euler-Lagrange equation: with the boundary conditions where  = ( 1 ,  2 ) denotes the unit outernormal vector of Ω.
Using the gradient descent method, we are able to derive the associated heat flow for the model ( 13): In this work, we use the finite difference scheme to discretize (17); see [5,39,40] for more details.Denoting the step space by ℎ = 1, we employ the following discretization used in the implementations; see Table 1.
In the discretization, the notation [, ] = ((sgn  + sgn )/2) ⋅ min(||, ||) and  > 0 is near 0. Denoting the time step by , we get the following explicit computational scheme: Now we discuss the choice of the weighting parameter  in (18).Due to the strengths and weaknesses of the firstorder and second-order variation approach, it is desirable that the weighting parameter  = 1 along edges and in flat regions, emphasizing the restoration properties for the firstorder total variation.To emphasize the restoration properties of the second-order total variation in smooth regions, we want 0 ≤  < 1. Specially, the resulting algorithm of ( 18) is just the TV regularization method for Poisson noise removal problem when  = 1, while the resulting algorithm is the high-order TV regularization method when  = 0. Usually, we may compute the parameter  by using the information of the edges and smooth regions of the resulting image obtained by smoothing the observed image  with the low-pass filters such as the median filter and the Gauss filter.We have carried out some numerical experiments.We observe that the fixed  obtained by computing the observed image can give very good results.In this work, in order to detect edges and smoothing regions much better, we adopt the method for updating  as proposed in [40].The results obtained by carrying out various numerical examples show that the updating procedure behaves better for our model than the fixed .It is because, as the iteration proceeds, the edges and smoothing regions of recovered image are closer to the original image, then the parameter  computed by the updating scheme can be better suitable for restoration.So we employ the method in [40] for updating  in our numerical experiments.More precisely, suppose that   is the th iterative solution, and we update the parameter  as follows:  of high or low similarity between two images, the whiter SSIM map, and the closer between the two images.
In order to give a good comparison, we employ the timemarching gradient descent algorithm for the four models.In the time-marching gradient descent algorithm, the parameter  = 10 −4 is introduced to avoid divisions by zero.The initial guess is chosen to be the degraded image in all tests.We set the step size  as 0.1 in order to obtain a stable iterative procedure.We terminate the iterations for the methods and accept   as the computed approximation of the ideal image as soon as the maximum number of allowed outer iterations has been carried out or the relative differences between consecutive iterates  1 , In addition, we carry out many experiments with different regularization parameters in the models, andthe one with the best results is presented in this work.In this way, we have a fair comparison since we compare the restorations for four different methods.
In the first experiment, we consider the "Cameraman" image of size 256 × 256.The degraded image in Figure 1(b) is obtained by performing the blurring operation psfGauss (5,1) proposed in [49] on the original image in Figure 1(a) with the background  = 10 and adding the Poisson noise.Note that there is no parameter associated with the Poisson noise, but the noise magnitude depends on the absolute image intensities.The amount of noise in a region of the image increases with the intensity of the image there.The restored images by all four methods are shown in Figures 1(c)-1(f).From these figures, compared with the tikhonov regularization method, the TV regularization method, and the high TV regularization method, the proposed approach yields better results in image restoration since it avoids the staircase effect of the general TV methods, and the edges are preserved as well as the general TV methods.In Figures 2(c)-2(f), we have enlarged some details of the four restored images.As it is seen in the zoomed parts, the proposed method outperforms the other three methods.In Figures 2(c)-2(f), we also present the SSIM maps of the restored images in Figures 2(g)-2(j).Note that the SSIM maps of the restored image by the proposed method are whiter than the  other methods; that is, our method can get better results.In Table 2, we compare their restoration results in SNRs and ReErrs.We observe from Table 2 that both the SNR and ReErr values of the restored images by the proposed method are better than those by the other three methods.Moffat blur is considered in the second example.It is known that the PSF of an astronomical telescope is often modeled by the Moffat function.The 256 × 256 "Lena" image shown in Figure 3(a) is degraded by the Moffat blur psfMoffat (5,2,5) proposed in [50] with  = 10 and the Poisson noise.The information of restored images by the four methods is displayed in Figures 3(c)-3(f).From the visual quality of restored images, the proposed regularization method is better than the other three methods.In Figures 4(c)-4(f), we display the zoom parts of the restored images (the part shown as the white rectangle in Figure 3(a)).From Figures 4(c)-4(f), we see that the restored image obtained by our method has more details than those by the other methods.The comparison of SSIM maps shown in Figures 4(g)-4(j) also proves that our method can get better results.We report the SNR and RelErr values in Table 2. From the table, we know that our method behaves much better.
In the third example, the 256 × 256 original astronomical object shown in Figure 5 In Figures 5(e)-5(h), we show the restored images by the four different methods.From these figures, we observe that the restored image by our model is more better than the other models in terms of the staircase effect and edge preservation.In addition, we plot the choice of  in Figure 5(d).From the figure, we see that only small jumps are suppressed with the second-order regularization.The SNRs and relative errors are reported in Table 2.The SNR by our method is higher than those by the other methods, and the relative error by our method is smaller than those by the other methods.From these figures and Table 2, we observe that the proposed method outperforms the tikhonov regularization method, the TV regularization method, and the high TV regularization method.

Conclusions
In this paper, we propose a new variational model to restore blurred images corrupted by Poisson noise.Based on the good feature of high-order functional, we propose a model by adding an extra high-order functional term in the total variation model.Our model combines advantages of the firstorder and second-order total variation.It can substantially reduce the staircase effect, while preserving edges in the restored images.The issues of existence and uniqueness of a minimizer for this variational model is discussed.At last, we employ a gradient descent method to solve the associated Euler-Lagrange equation.The numerical experiments show that the proposed method outperforms some existing restoration methods in terms of the SNR, relative error, and SSIM map for Poisson noise removal problem.The comparisons between the reconstructed images obtained by the four methods show that the proposed one can alleviate the staircase effect significantly while preserve edges.

Figure 2 :
Figure 2: Comparison of a small portion of the restored images and SSIM maps for Example 1.

Figure 3 :
Figure 3: The "Barbara" image is degraded by the Moffat blur and Poisson noise.
(a) is degraded by the given PSF in Figure 5(b) and the Poisson noise.The degraded image is displayed in Figure 5(c).The relative error between the noisy image and the original image is 0.3150.

Figure 4 :
Figure 4: Comparison of a small portion of the restored images and SSIM maps for Example 2.

Table 2 :
Restoration results with different methods for Example 1, Example 2, and Example 3.