Second-order TGV model for Poisson noise image restoration

Restoring Poissonian noise images have drawn a lot of attention in recent years. There are many regularization methods to solve this problem and one of the most famous methods is the total variation model. In this paper, by adding a quadratic regularization on TGV regularization part, a new image restoration model is proposed based on second-order total generalized variation regularization. Then the split Bregman iteration algorithm was used to solve this new model. The experimental results show that the proposed model and algorithm can deal with Poisson image restoration problem well. What’s more, the restoration model performance is significantly improved both in visual effect and objective evaluation indexes.

where the first term is the total variation regularization term, the second term is fidelity term and λ is a positive fidelity parameter. Although the TV model is effective in image processing, specially for Gaussian white noise. However, this model is not effective in restoring Poisson noise images such as astronomical (Starck and Murtagh 2007), biomedical (Dey et al. 2006;Sarder and Nehorai 2006;Willett and Nowak 2003), and photographic imaging (Foi et al. 2005). For this reason, scholars have proposed some new methods to deal with Poisson noise image problems. Le et al. (2007) proposed a new total variation model to recover image corrupted by Poisson noise, the new total variation model with fidelity term is suitable for Poisson noise. The new model can be written as The authors used the gradient descent method to obtain its optimum solution. However, this method can not obtain optimal approximation when the image is both high intensity noise and low intensity features. In Sawatzky et al. (2013), an efficient EM-TV algorithm is presented to speed the computation of the optimization problem (2). In addition, alternating split Bregman iterative algorithm (Setzer et al. 2010) is also used to solve the question (2), since this algorithm does not contain iterations and also not produce negative values. In addition, Figueiredo and Bioucas-Dias (2010) proposed an approach based on alternating direction optimization method for deconvolving Poissonian images.
Recently, based on Chavent and Kunisch (1997) and Liu and Huang (2012) proposed another new total bounded variation-based Poissonian images restoration model Although the above model is better than the total variation model and has a competitive superiority, there also exist some shortcomings, for example, sometimes it will cause undesired oil painting artifacts. In order to avoid the staircase effect, many methods have been proposed. A well-known method to eliminate staircase effect is the TGV (Bredies et al. 2010;Bredies and Valkonen 2011;Bredies et al. 2013) regularization. The TGV regularizer can effectively eliminate the staircase effect but there are still some shortcomings, it tends to introduce some blurring on image edges and texture regions as the existence of high-order derivative term. More seriously, some small details will be lost during the denoising. For this reason, we consider combine the TGV and u 2 as one regularization term to solve the poisson noise image restoration problem.
The rest of this article is organized as follows. In "Total generalized variation (TGV)" section, we briefly review the total generalized variation (TGV). The proposed model and algorithm are presented in "The proposed Poisson noise recovering model and algorithm" section. In "Experimental results and discussions" section, experimental results are illustrated to show the consistent performance of the proposed method. Finally, conclusions are given in "Conclusions" section.
(2) min Bredies et al. (2010) proposed the concept of total generalized variation (TGV), which is considered to be the generalization of TV. For convenience, some concepts of TGV are given as follows.

Total generalized variation (TGV)
Definition 1 (Bredies et al. 2010) Let ⊂ R d be a domain, k ≥ 1 and α = (α 0 , . . . , α k−1 ) > 0.Then, the total generalized variation of order k with weight α for u ∈ L 1 loc (�) is defined as the value of the function where Sym k (R d ) denotes the space of symmetric tensors of order k with arguments in R d , and α l are fixed positive parameters.
Definition 2 (Bredies et al. 2010) The space of bounded generalized variation is defined as Here BGV k α (�) is a Banach space independent of the weight vector α.
Definition 3 (Bredies et al. 2010) The "dualization" in the definition of the functional TGV k α can also be informally interpreted in terms of iterated Fenchel duality.
Note that the tensor field u l are in different spaces for varying l. Moreover, the operator ε u l−1 denotes the symmetrized gradient operator In this paper, we use k = 2 in the proposed model. Thus, the second-order TGV can be written as where S d×d denotes the space of symmetric d × d matrices. And the first and second divergences are defined as (4) In addition, according to Bredies and Valkonen (2011), the energy term TGV 2 α can be formulated as where ε(p) can be separately expressed as

The proposed model for Poisson noise image
We assume that u ∈ R N + is the original image, f ∈ R N is the observed image, K ∈ R N ×N is a linear blurring operator related with the spread point function (PSF). Then the degradation model can be described as where P denotes the Poisson distribution function.
Based on Le et al. (2007), we define the Bayes Law as follows.
According to (14), for each u ∈ , we have Next, we assume that the prior distribution P(u) is TGV and u 2 , which can be written as where λ is the regularization parameter. Thus, we obtain a model for restoring the Poissonian noise image as follows.
where l s is the indicator function of set S (10) By reformulating TGV as a minimization in the discrete setting, the proposed model can be written as

The split Bregman algorithm for Poisson noise removal
The split Bregman algorithm (Goldstein and Osher 2009;Wang et al. 2008) has been widely used in image processing, which is easy to be realized and has fast convergence (Cai et al. 2009;Jia et al. 2009). Therefore, we use split Bregman algorithm to solve our minimization problem (19). Firstly, by introducing new auxiliary variables w, x, y and z, the problem (19) can be reformulated as the following constrained optimization problem For the above constrained problem (20), we transform it into the corresponding unconstrained problem where µ i (i = 1, . . . , 4) are positive penalty parameters. Thus the split Bregma iterative algorithm for solving the question (20) can be described as where the updates of the multipliers b 1 , b 2 , b 3 , b 4 is described as follows Since the updates of b k 1 , b k 2 , b k 3 , b k 4 are merely simple calculations, then the minimization question (22) can be divided into the following several subproblems: Given initial value the split Bregman algorithm can be written as For w-subproblem, note that it is separable with respect to each component. It is easy to solve and the solution of w may be written as As for solving x, y-subproblem, we can directly obtain the solutions by using shrinkage operator: The x-subproblem can be solved by The solution of the y-subproblem is similarly obtained The (u, p)-subproblem is a saddle-point problem, which can be divided into the following two subproblems: 1. For u, we have which can be solved by considering the following normal equation Finally, u is solved by 2. For the sub-problem p, it can be written as the following minimization problem where p = (p 1 , p 2 ) T is a 2 × 1 vector, ε(p) is a 2 × 2 matrix. For p 1 , it can be solved by considering the following linear system Therefore, Similarly, we can obtain the solution of p 2 as

Experimental results and discussions
In this section, we illustrate some numerical results of the proposed model for the Poisson noise removal problem. We compare our method with the one proposed in Figue (29) � . (30) We terminate the iterations for these methods by the following stopping criterion The quality of the restoration results is compared quantitatively by using the Signal-to-Noise Ratio (SNR), the Peak Signal-to-Noise (PSNR), the relative error (RelErr) and the Structural SIMilarity index (SSIM). They are defined as follows where u and u are the ideal image and the restored image, respectively.
where µ u and µ u are averages of u and u, respectively. σ u and σ u are the variance of u and u, respectively. σ u u is the covariance of u and u. The positive constants C 1 and C 2 can be thought of as stabilizing constants for near-zero denominator values. Generally speaking, the more bigger value of SNR, PSNR or the smaller value of RelErr is, the better quality of the reconstructed image is. The Poissonian images used for our experiments are generated as follows: the original images are convoluted with a blur kernel and additionally contaminated by Poisson noise, here we use the poissrnd function in MATLAB's Statistics Toolbox after blurring the true images with the given point spread functions to generate the blurred and noise images.
The selection of the regularization parameters highly affects the image restoration results, and related to make the fair comparison with different denoising models. The penalty parameters µ which relies on unknown noise level highly influences the speed of the algorithms. In experiments, we set µ = [0.01, 0.001] in the PIDAL algorithm. In the PID-Split algorithms, we choose µ = [0.0004; 0.1, 0.0001]. In the TGV model, we set µ = [0.1; 10, 5, 3]. The penalty parameter in the proposed method is empirically set µ = [0.1; 0.6; 0.1; 0.02]. Thus, we may have a good restoration results.
In the first experiment, we used the image "Woman" (512 × 512) in Fig. 1a. We perform the blurring operation psfGauss(5, 2) proposed in Nagy et al. (2004) on the original image and add the Poisson noise to the blurred data to generate the degraded image in Fig. 1b. The parameter of this test, we set β = 120 in PIDAL algorithm, β = 6, = 0.01 for PID-Split algorithms,due to the TGV model we set β = 450, α = [8, 10], set β = 54, = 0.001, α = [16,9] for the proposed model. The pictures of Fig. 1c-f are the restoration images, which represent the difference between the three methods. From these pictures, we can see the proposed model have more advantages. In order to more effectively reflect the experiment result, Fig. 1g-j present the residual images refer to the difference of the original image and the restoration image. From these pictures, we can see that the proposed model can preserve more details than other methods. In the Table 1, the SNR, PSNR, RelErr and SSIM values of the restored images by the proposed model are better than other methods.
In second experiment, we use the image "House" with size (512 × 512) in Fig. 2a, which contains a lot of edge details. We perform the blurring operation with radius 5 (psf = ones(5, 5)/25) on the original image and add the Poisson noise to the blurred data to generate the degraded image in Fig. 2b. As for parameter selection, we choose β = 200 for the PIDAL algorithm, β = 20, λ = 0.00001 for the PID-Split   Fig. 2c-f, we can see that the proposed model compared to the PIDAL method and PID-Split algorithms have better restoration results. In Fig. 2g-j, we have enlarged some details of the images, which can be clearly see the advantages of the proposed model for the recovery of edge details. The SNR, RelRrr and SSIM values in Table 1 showed that the proposed model have a better restoration result. In order to further verify the validity of the new model, in third experiment, we use the image which contains a lot of edge details. We perform the blurring operation by a line motion blur. The point spread function for the linear motion blur is returns a filter to approximate, once convolved with an image, the linear motion of a camera by r pixels, with an angle of θ degrees in a counter-clockwise direction. In this example, r = 2 and θ = 45, then add the Poisson noise to the blurred data to generate the degraded image in Fig. 3b. The parameters choose as the same as those in second experiment and also may be adjusted.
Finally, let us choose a human brain MR image of size 240 × 240 as the test image. We zoom in a marked close-up region which is abundant in texture-like features to better visual comparison. We can clearly see that the produced textures by our proposed method are better quality than the other methods from Fig. 4.

Conclusions
In this paper, we investigate the second-order total generalized variation with a quadratic regularization to deal with the Poissonian images restoration problem. The proposed model is solved efficiently by split Bregman iterative algorithm in this way the calculation speed is fast. Numerical results show that our proposed method is particularly advantageous for restoration the Poisson images in terms of SNR, SSIM and RelErr quality compared to other methods. In the model, the parameters selection is a difficult problem which needs further study.