Abstract

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.

1. 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 [14]. Over the last decades, most of the literature deal with the additive noise and multiplicative noise model [511]. Given a blurred and noisy image or , where is the original image, 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 [1219] 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) [2024]. In a Poisson noise model, the recorded image is a realization of a random variable . A statistical model for is given by where 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 that these random variables are independent, the probability distribution of , for given , , and , can be written as 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 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 , we obtain Therefore, the maximum likelihood estimator of the original image is the minimizer with respect to of the negative-log Poisson likelihood functional .

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 is a regularization function and is the regularization parameter that controls the balance between and . 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 [2628].

It is shown in [29] that the choice of the regularization function 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 , 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 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 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 so-called 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 second-order 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.

2. 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 . It is clear that 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 , which has the effect of penalizing reconstructions with large norm. For 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 .

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 [2124] for more details. In the other hand, some Newton-related gradient methods were proposed for deblurring Poisson-corrupted 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 [3133] 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 . It leads to present the following functional for restoring blurred images corrupted with Poisson noise: where the parameter 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 such that ; in particular, must be positive almost everywhere. Some basic notations and properties concerning the space and 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 with Lipschitz boundary. Assume that is a positive bounded function and is positive definite. Then for such that has a unique minimizer for the model (13).

Proof. First, a straightforward computation can show that is bounded below by . Since the penalty function , we obtain that is bounded below. Thus, we can choose a sequence such that . Hence, and are bounded. Using Jensen’s inequality, we have which indicates that is bounded. This and the boundedness of yields that the sequence is bounded in . Following the compactness of in , we deduce that there exists such that a subsequence converges to pointwise almost everywhere in . By the lower semicontinuity of the space, we obtain that . Since is bounded below, we get that by using Fatou’s lemma. Therefore, and minimizes for all such that .
Since is positive definite and is positive, it immediately follows that is strictly convex. Obviously, 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.

3. 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 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 , we employ the following discretization used in the implementations; see Table 1.

In the discretization, the notation and 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 first-order and second-order variation approach, it is desirable that the weighting parameter along edges and in flat regions, emphasizing the restoration properties for the first-order total variation. To emphasize the restoration properties of the second-order total variation in smooth regions, we want . Specially, the resulting algorithm of (18) is just the TV regularization method for Poisson noise removal problem when , while the resulting algorithm is the high-order TV regularization method when . 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: where , and it means that only the absolute largest jumps are unaffected of the high-order regularization. As reported in [40], for large and small values of the parameter is closer to , and for intermediate values of the parameter approaches 0, which means that the high-order filter dominates the computation and the staircase effect is suppressed. Since only small jumps should be suppressed with the high-order regularization, it is a good choice for .

We are now in a position to describe the time-marching gradient descent algorithm (Algorithm 1) for restoring blurred images corrupted by Poisson noise.

Input , , , MaxIter, and . Initialize
For k = 1: MaxIter
 compute using the explicit scheme (18).
 update the parameter according to (19).
end for

4. Numerical Experiments

In this section, we present some numerical results to illustrate the performance of the proposed model for Poisson noise removal problem. We compare our model with the one proposed in [21] (TV method), the one proposed in [31] (Tikhonov method), and the one proposed in [47] (HTV method). All computations of the present paper were carried out in Matlab 7.10. The results were obtained by running the Matlab codes on an Intel(R) Core(TM) i3 CPU (2.27 GHz, 2.27 GHz) computer with RAM of 2048 M.

The quality of the restoration results with different methods is compared quantitatively by using the signal-to-noise ratio (SNR), the relative error (ReErr), and structural SIMilarity index (SSIM). They are defined as follows: where and are the ideal image and the restored image, respectively, and where and are averages of and , respectively, and are the variance of and , respectively, and is the covariance of and . The positive constants and can be thought of as stabilizing constants for near-zero denominator values. The SSIM is a well-known quality metric used to measure the similarity between two images. The SSIM method was developed by Wang et al. [48] and is based on three specific statistical measures that are much closer to how the human eye perceives differences between images. In the following experiments, we will use SSIM map to reveal areas 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 time-marching gradient descent algorithm for the four models. In the time-marching gradient descent algorithm, the parameter 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 satisfy 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 . The degraded image in Figure 1(b) is obtained by performing the blurring operation psfGauss proposed in [49] on the original image in Figure 1(a) with the background 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 “Lena” image shown in Figure 3(a) is degraded by the Moffat blur psfMoffat proposed in [50] with 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 original astronomical object shown in Figure 5(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.

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.

5. 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 first-order 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.

Conflict of Interests

All of the coauthors do not have a direct financial relation with the trademarks mentioned in our paper that might lead to a conflict of interests for any of the coauthors.

Acknowledgment

L. Jiang and J. Huang are supported by NSFC (10871034), X.-G. Lv and J. Liu are supported by NSFC (60973015 and 61170311), Sichuan Province Sci. & Tech. Research Project (2011JY0002 and 12ZC1802), and Chinese Universities Specialized Research Fund for the Doctoral Program (20110185110020).