POISSON IMAGE DENOISING BASED ON FRACTIONAL-ORDER TOTAL VARIATION

,

1. Introduction.In photon-counting devices such as positron emission tomography [54] and astronomical imaging [23], noise is inevitably dependent on the number of photons.Since the number of photons follows a Poisson distribution, this type of noise is often referred to as Poisson noise.In terms of image denoising, the smaller the peak value of an original image is, the noisier the image becomes when corrupted by Poisson noise.Based on the statistical property of the Poisson distribution and maximum likelihood estimation [17], Poisson denoising and restoration models often adopt the following data fitting term (1) min where f is an observed (noisy) image, u is an image to be recovered, and Ω is a bounded open domain in R 2 (typically rectangular for images).
Recently, many regularization techniques have been investigated for Poisson denoising, including Tikhonov regularization [51], total variation (TV) [44], and highorder TV [8].In particular, the TV-based Poisson denoising model [24] can be expressed as (2) min where β is a positive parameter.There are a number of efficient algorithms for solving the TV-regularized Poisson denoising model (2) such as an alternating extragradient method [2], an inexact iteratively reweighted approach [13], a primal-dual algorithm [58], a dual algorithm [47], an augmented Lagrangian method [18], and the split Bregman method [30,59,68].In addition to denoising, other applications involving TV and Poisson noise include image decomposition [36], image deblurring [31,48,56,66], and phase retrieval [9].As one of the most popular regularizations, TV performs well on piecewise constant images for preserving edges while removing noises.However, it causes staircasing artifacts and the loss of image contrasts, latter of which is particularly obvious for piecewise smooth images [32].To address these issues, some high-order regularizations for the Gaussian denoising problem have been introduced, including non-local total variation [21], TV combined with a fourth-order diffusive term [8], Euler's elastic models [15,49,50,61], a mean curvature model [69], a total generalized variation (TGV) model [4,55], a second-order TGV [26], and TGV together with sparsity and/or shearlets [19,27,40].Specifically for the Poisson denoising, some high-order methods include PDE based models [62,67], a hybrid regularization combining TV with a fourth-order variation [20], and TGV-based Poisson denoising model [29].In particular, the Lysaker-Lundervold-Tai (LLT) model [32] was proposed to minimize an objective functional involving the Laplacian, thus leading to a fourth-order partial differential equation (PDE).Despite of its effective suppression of staircasing artifacts, the LLT method introduces speckle artifacts.
Different from other types of high-order TV, the fractional-order total variation (FOTV) uses the derivatives with order greater than or equal to one, leading to a compact discrete form and thereby bringing computational convenience.It has been empirically shown that FOTV is able to suppress both staircasing and speckle artifacts [1].By taking the neighboring image intensities into account, FOTV can preserve local geometric characteristics.For example, FOTV can preserve textures [37,38,39] and image contrasts [60].By extending the definition from a 2D domain to a graph, the FOTV regularization has been successfully applied to the electroencephalography (EEG) source imaging [28,41,42].Other applications involving the FOTV regularization include image restoration from Gaussian noise [10,14,25,33,60], multiplicative noise removal [52,53,64], impulse noise removal [7], image decomposition [34], image reconstruction [65], and fluorescence microscopy image deconvolution [43].
Motivated by the success of FOTV in preserving image features and contrasts, we propose a Poisson denoising model based on FOTV to preserve high-order smoothness of an image.To the best of our knowledge, this is the first time that FOTV is considered for Poisson noise.The major contributions of the paper are three-fold: (1) We provide a rigorous discussion of the existence and uniqueness of the solution to the proposed model, based on the properties of the corresponding function space.(2) We derive three algorithms, which can shed lights on computational efficiency for such a nonlinear problem.Convergence guarantee of each proposed algorithm is also analyzed, which offers some guidance on parameter selection.(3) We conduct extensive experiments including a toy example and several natural images to demonstrate that our approach has superior performance over the state-of-the-art methods in Poisson denoising.
The rest of the paper is organized as follows.Section 2 provides both continuous and discrete formulations of the FOTV regularization together with analysis on existence and uniqueness.Three numerical schemes are detailed in Section 3 including the Chambolle-Pock's primal-dual (PD) method, a forward-backward splitting (FBS) scheme, and an alternating direction method of multiplier (ADMM).Experimental results on a synthetic image and five standard images are presented in Section 4 to illustrate the effectiveness and efficiency of our proposed methods.Finally, Section 5 concludes the paper.
2. Our proposed model.In this section, we first introduce the definition of FOTV and its related function space.Then we propose a variational model of FOTV-regularized Poisson denoising.Using the properties of FOTV, we can show that there uniquely exists a solution to our model.To ease the exposition of numerics, we provide discrete gradients and divergence for the implementations of FOTV.
2.1.Fractional-order total variation.Let C α 0 (Ω) with α > 0 be the space of α-order continuously differentiable functions defined on Ω with compact support.The α-order total variation is then defined as (3) Here the α-order divergence is given by where ∂ α φi ∂z α with z ∈ {x, y} is the α-order derivative of φ i along the z-direction for i = 1, 2. Due to the convenience in numerical implementation, we adopt the Grünwald-Letnikov (GL) fractional-order derivative [63].Based on the seminorm defined in (3), we define the space of functions with bounded α-order total variation, denoted by BV α (Ω) := {u ∈ L 1 (Ω) : T V α (u) < ∞}.Similar to the space of functions with bounded variations, it can be shown that BV α (Ω) is a Banach space and any functional J : BV α (Ω) → R + is lower semi-continuous [60].
where T V α (u) is defined in (3) and β > 0 is a weighting parameter to balance the regularization term and the data fitting.Theorem 2.1 establishes the existence and uniqueness of solutions to the proposed model (4).
By using the fact that FOTV is convex, the proof is similar to [24,29,30], thus omitted.
2.3.Discrete fractional-order gradient and divergence.For u ∈ BV α (Ω), we discretize the image domain Ω as a rectangular grid {(x i , y j ) : 1 ≤ i ≤ N, 1 ≤ j ≤ M }.Then an image can be represented as a matrix in the Euclidean space R N ×M , denoted as u i,j = u(x i , y j ).Based on the GL fractional-order derivatives [63], the discrete fractional-order gradient is defined as where the discrete gradients D α 1 u, D α 2 u ∈ R N ×M along the x-axis and the y-axis are given by (6) Here K is the number of neighboring pixels that are used to compute the fractional-order derivative at each pixel.The coefficients {C α k } K−1 k=0 are given by Γ(k+1)Γ(α+1−k) with the Gamma function Γ(x).Then the discrete fractionalorder total variation of u ∈ Ω is defined as ( 7) Using the relation that (∇ α ) * = (−1) α div α , the discrete fractional-order divergence div α p ∈ R N ×M for p = (p (1) , p (2) i+k,j + p i,j+k ).

2.4.
Connections to the existing models.For α = 1 in (6), it is straightforward to have , which are the first-order backward differences along the coordinate axes.In this case, FOTV is reduced to the standard TV when backward difference is adopted as a discretization scheme.
Similarly, (D 2 2 u) i,j = u i,j −2u i,j−1 +u i,j−2 .As a consequence, FOTV with α = 2 can be expressed by , where u xx and u yy are discretized by (11) ( This is equivalent to the LLT model [32].In summary, FOTV is a general regularization framework including TV and LLT as its special case for α = 1 and α = 2, respectively. 3. Algorithms.For the purpose of implementation, we consider a discretized version of the Poisson denoising model (4) as follows, (12) min where 1 Ω denotes the indicator function on a discrete grid of the image domain Ω.We propose three efficient numerical algorithms based on three popular optimization methods, i.e., PD, FBS and ADMM.By discussing various optimization schemes, we can select the one with the optimal computational efficiency.More importantly, we can provide some rule of thumb regarding the optimal scheme for such nonlinear optimization problem.
3.1.Primal-dual algorithm.We rewrite the proposed model ( 12) as a minimax problem in order to apply the Chambolle-Pock primal-dual (PD) method [6].It leads to two subproblems, each of which can be solved efficiently.The convergence of the scheme can be guaranteed as discussed in [6].By letting we express the proposed model ( 12) in the following form, It follows from the property of convex conjugate F * that the convex conjugate function of F is given by ( 14) We can then rewrite (13) as a minimax form Directly applying the Chambolle-Pock algorithm [6] for the minimax problem (15) yields the following iterative scheme, ( 16) where σ, τ and θ are positive tuning parameters.Note that the u-subproblem is derived due to the relationship between the two operators ∇ α and div α , i.e., In what follows, we describe the respective closed-form solutions for the two subproblems in (16).Specifically, the p-subproblem is equivalent to ( 17) It has a closed-form solution, given by ( 18) , otherwise, for p = p k + σ∇ α u k .As for the u-subproblem, the solution can be obtained by solving the corresponding Euler-Lagrange equation of the form where ũ = u k − τ (−1) α div α p k+1 .By using the quadratic formula, we get a nonnegative solution of u, i.e., We summarize the FOTV regularized Poisson denoising via PD in Algorithm 1.The convergence of this algorithm is guaranteed by the following theorem, whose proof is presented in [6].
Remark.We seek an upper bound of L that enables the convergence of the proposed PD algorithm as well as provides guidance on choosing appropriate values of the parameters τ and σ.For this purpose, we estimate For the special case of TV, i.e., α = 1, K = 2, we have L 2 ≤ 8, which is consistent in [5].Following discussion in Section 2.4, FOTV with α = 2 has fewer terms than anisotropic LLT and hence the upper bound of L with α = 2, K = 3 (L 2 ≤ 36) is smaller than that of LLT (L 2 ≤ 64 in [11]).
3.2.Forward-backward splitting algorithm.Motivated by the works [46,47], we consider the first-order optimality condition of (12), i.e., (20) v where v ∈ ∂J(u) is a subderivative of the functional J at u.At the k-th iteration, we adopt the following discretization scheme where the constant term 1 Ω is approximated by u k+1 u k and this method is referred to as the forward-backward splitting (FBS) method.It follows from the optimality condition (21) that u k+1 is an optimal solution of the problem, i.e., (22) u k+1 = arg min We want to point out that the convergence of the FBS scheme ( 21) can be proved in a similar way as in [46].
We then describe how to solve for the subproblem of ( 22) by introducing an auxiliary variable z in such a way that ( 22) is equivalent to (23) min To efficiently handle the L 1 term of this model, we choose a dual formulation [5].Specifically, we construct the corresponding Lagrangian functional Since L is separable with respect to u and z, the dual functional is given by ( 24) The minimizer u of ( 24) can be expressed as By the dual norm of L 1 -norm, we can get a closed-form solution of z using a characteristic function inf where S = {λ | λ ∞ ≤ 1}.Therefore, the dual problem of ( 22) becomes ( 26) whose (i, j)-th entry is denoted by θ(λ; u k ) i,j .Then the Karush-Kuhn-Tucker conditions yield the existence of a Lagrange multiplier γ i,j ≥ 0, associated to each component λ i,j .For each index pair (i, j), we have In either case, we have Therefore, we consider to solve for λ i,j from the following equation via the semi-implicit gradient descent algorithm.More specifically, given a step size τ > 0, the l-th iteration goes as follows, (28) Once λ is obtained, the update for u is given by ( 25).The entire FBS scheme (22) with a dual algorithm for the subproblem is summarized in Algorithm 2.
3.3.Augmented Lagrangian method.We adopt the alternating direction method of multipliers (ADMM) [3] to solve the proposed model (12).In particular, we introduce an auxiliary variable z to reformulate (12) as, (30) min The augmented Lagrangian functional corresponding to ( 30) is where λ is Lagrangian multiplier and µ is a positive parameter.Then ADMM yields the following algorithm, To solve the u-subproblem in (32), we consider the Euler-Lagrange equation of (31) with respect to u, which has the form In order to efficiently solve the nonlinear equation ( 33), we replace u in the denominator by the previous iteration u k and solve for u k+1 by using the fast Fourier transform (FFT) with periodic boundary conditions.More specifically, u k+1 is given by ( 34) where F is the Fourier transform and • is the component-wise multiplication.Note that the update in ( 34) is an approximation to the solution of the nonlinear problem (33).We can iterate (34) for a few times in order to get a better solution, but empirical evidence shows that one iteration is sufficient to give a reasonable solution.In addition, there is a closed-form solution for the z-subproblem given by ( 35) where shrink(s, γ) := sgn(s) max{|s| − γ, 0}.The ADMM algorithm for solving the FOTV-based Poisson denoising model is summarized in Algorithm 3. Convergence of ADMM for a general convex problem was proven by Eckstein and Bertsekas [16], which is described as follows.
Theorem 3.2 (Eckstein and Bertsekas [16]).Consider a convex problem of the form where F, G are proper convex functionals.Suppose we are given µ > 0 and for all k.Then if (36) has a solution u * , the sequence {u k } converges to u * .If (36) has no solution, then at least one of the sequences {z k } or {λ k } is unbounded.
Remark.It is straightforward to see that the functionals F, G in our problem (12) are proper and convex.The existence of solutions is guaranteed by Theorem 2.1.
If u and z subproblems can be solved with sufficient accuracy, the convergence of Algorithm 3 is guaranteed by Theorem 3.2.

Numerical experiments.
In this section, we present extensive numerical results of the proposed Poisson denoising method, in comparison with some state-ofthe-art methods.Since Poisson noise depends on the pixel intensity, the noise level can be controlled by the peak intensity of the original image.In order to examine .By the theory of kernel functions of FOTV [60], the optimal fractional order for denoising this synthetic image is α = a + 1 = 1.8.In Figure .1, we show this ground-truth image together with two noisy images of peak values at 55 and 255, respectively.One can see that the image with peak 55 is more noisy compared the one with peak 255.We present the denoising results at peak 255 using different fractional orders α = 1, 1.8, and 2.4.The bottom row of Figure . 1 illustrates staircasing artifacts when α = 1 in FOTV (i.e., TV) and the over-smoothing phenomenon when α = 2.4.
Next we compare three minimizing algorithms for the proposed FOTV model, followed by comparison with some popular Poisson denoising methods, including NPTool [22], NL-PCA [45], and BM3D [35].NPTool [22] is a MATLAB toolbox for the nonnegative image restoration with Newton projection methods, among which we use the total variation model.NL-PCA [45] is based on the idea of nonlocal patches together with principal component analysis (PCA) for Poisson noise reduction.BM3D [35] for Poisson denoising involves the Anscombe transformation that converts Poisson noise into Gaussian noise, followed by the Gaussian denoising using BM3D [12].
We use the peak-signal-to-noise-ratio (PSNR) and structural-similarity (SSIM) index [57] for quantitative comparison.The PSNR is defined as (38) PSNR(X, where R is the maximum peak value of the original image X and Y is the restored image.As for SSIM, we define the local similarity index computed on two small patches x and y, , where µ x /σ 2 x , µ y /σ 2 y are the average/variance of x, y, σ xy is the covariance of x, y, and two positive parameters c 1 , c 2 are set to avoid the denominator being zero.The overall SSIM is the mean of local similarity indexes, i.e., (40) SSIM(X, Y ) = where x i , y i are the corresponding patches indexed by i and P is the number of patches.All numerical experiments are performed under Windows 7 and MATLAB R2017b running on a desktop with Intel(R) Core(TM) i7-4790 CPU @ 3.60GHz.

4.1.
Comparison of proposed algorithms.We compare the efficiency of the three proposed algorithms, i.e., Algorithms 1-3, by minimizing the same objective function (12).In all experiments, the stop criterion is For the synthetic image, we fix α = 1.8 (ground-truth value), K = 20, and β = 7.We optimize other algorithmic parameters for each algorithm when plotting the energy and PSNR values with respect to time (in seconds).The results are provided in Figure 2, which shows that Algorithm 2 is the fastest algorithm to minimize the energy, but Algorithm 1 achieves the largest PSNR.Another example using Barbara image at peak 55 is also included in Figure 2, where we set α = 1.6,K = 20 and β = 12.We observe that Algorithm 3 achieves the best results in terms of PNSR.As ADMM involves the least number of parameters compared to the other two algorithms, we will use it (Algorithm 3) for the rest of experiments.All the plots in Figure 2 numerically demonstrate the convergence of each algorithm.
In Figure 3, we illustrate how the fractional order α affects the image restoration.For each order in FOTV, we choose the best combination of parameters from the parameter sets of β ∈ {0.1, 1, 7, 10} and µ ∈ {10 −5 , 10 −4 , 10 −3 , 2 × 10 −3 }.We then plot the largest PSNR as a function of fractional order for two peak levels 55 and 255.For the synthetic image, when peak is 255, the optimal order indeed occurs at the theoretical value α = 1.8.However, it is not the case for peak value 55.This phenomenon occurs due to the presence of the truncation error when we approximate the fractional order derivative by only taking K adjacent pixels into consideration.In other words, when the noise level is high, we could increase the accuracy in discretization of the fractional order derivatives.We also include the PSNR value versus fractional order for the Barbara image in Figure 3  shows advantages of fractional order derivatives over integer order ones.However, with complicated image contents, natural images usually have less obvious denoising improvements by using FOTV compared to synthetic ones, e.g., 5 db PSNR gain for the synthetic image versus 0.5 db for the Barbara image when comparing the optimal fractional order and the standard TV (α = 1).is a TV-regularized denoising method and NL-PCA/BM3D are designed for texture images.Although BM3D works well on the Barbara image, it yields the worst restoration performance on the Mandril image.Since both the Barbara and Mandril images have rich textures, it indicates that BM3D is not robust to textures.The computational times for these methods are provided in Table 2.We show some visual results of Poisson denoising in Figures 4-7, each with a zoom-in patch and the corresponding PSNR/SSIM values.In particular, Figures 4-5 illustrate staircasing artifacts of NPTool and oversmoothing artifacts of BM3D.Oversmoothing is very obvious when BM3D is applied on the Mandril image in Figure 6, which is consistent with Table 1 that BM3D performs consistently poorly on the Mandril image.Figure 7 shows the denoising results of the Cameraman and Penguin images with peak values of 255 and 155, respectively.5. Conclusions.In this paper, we proposed a novel Poisson denoising approach based on the fractional-order TV regularization.Using the properties of the function space of bounded FOTV, we proved the existence and uniqueness of the solution to the proposed model.Based on the primal-dual, forward-backward splitting, and ADMM optimization methods, we proposed three numerical algorithms, each with guaranteed convergence.A variety of numerical experiments were performed  to justify the effectiveness of the proposed approaches in Poisson denoising.In particular, we used a synthetic image to investigate the optimal order of FOTV both theoretically and numerically.Furthermore, we compared our method with the state-of-the-art in Poisson denoising by testing five natural images with various levels of Poisson noise.Experimental results demonstrated the superior performance of the proposed approach over other competing ones.In the future work, we will explore acceleration techniques to further reduce the cost of computing fractionalorder derivatives while improving the restoration accuracy.

InverseFigure 1 .
Figure 1.A synthetic example.Top: the synthetic ground-truth image and two noisy images with peak values 55 and 255, respectively.Bottom: FOTV Poisson denoising results ( via Algorithm 3) with the respective fractional order 1, 1.8, and 2.4 for the case with peak value 255.

Figure 2 .
Figure 2. Algorithmic comparison in terms of energy (left) and PSNR (right) by denoising the noisy synthetic image with peak value of 255 (top) and the Barbara image with peak value of 55 (bottom).