# Poisson Noise Removal Scheme Based on Fourth-Order PDE by Alternating Minimization Algorithm

Authors: Weifeng Zhou; Qingguo Li
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

To overcome the staircasing effects introduced by the TV regularization in image
restoration, this paper investigates a fourth-order partial differential equation (PDE) filter for
removing Poisson noise. In consideration of the slow convergence property of the classical
gradient descent method, we adopt the alternating minimization algorithm for realizing this
scheme. Compared with the corresponding total variation based one, numerical simulations
distinctly indicate the superiority of our proposed strategy in handling smooth regions of
Poissonian images and improving the computational speed.


---

## Body

## 1. Introduction

Noise removal occupies an important position in image processing. The classical total variation (TV) regularization model proposed by Rudin et al. [1] can be characterized as(1.1)min⁡u⁡∫Ω|Du|+λ2∫Ω(u-f)2dx,
where Ω stands for a spatial domain with Lipschitz boundary, u and f are the clean image and the initial noisy image, respectively, and Du represents the distributional derivative of u. This strategy can do a better job in preserving the sharp edges while removing additive Gaussian noise. However, in the real world, many important images, for instance, medical imaging [2, 3], astronomy [4], were contaminated by Poisson noise more frequently than additive Gaussian noise. Although the ROF model (1.1) is effective for removing Gaussian noise, it is inferior in processing Poissonian images, which contain noise that is signal dependent. To remove Poisson noise effectively, Le et al. [5] introduced the following TV-regularization model(1.2)min⁡u⁡∫Ω|Du|+λ∫Ω(u-flog⁡u)dx,
where u must be positive almost everywhere. Numerical experiments demonstrate the better performance of scheme (1.2) for Poisson noise removal than model (1.1).Although model (1.1) and model (1.2) have been indicated to be capable of achieving a good tradeoff between noise removal and edge preservation, unfortunately, the unpleasant blocky effect inevitably emerges in the homogeneous flat regions of restoration images owing to the use of second-order scheme. One of the most well-known solutions to overcome this drawback is to use high-order partial differential equations instead of second-order ones. Up till the present, high-order diffusion strategies have made great progress in effectively suppressing Gaussian noise such as schemes in [6–15]. Thereinto, the model proposed by Lysaker et al. [10] for Gaussian noise removal has been paid more attention, which can be formulated as(1.3)min⁡u⁡∫Ω|D2u|+λ2∫Ω(u-f)2dx.
Numerical experiments obviously exhibit the superiority of the model mentioned above over second-order version in relieving the staircase artifacts.Motivated by the above models (1.2) and (1.3), we propose the following fourth-order PDE regularization scheme to suppress Poisson noise:(1.4)min⁡uE(u)=min⁡u⁡∫Ω|D2u|+λ∫Ω(u-flog⁡u)dx.
In consideration of the slow convergence speed of the traditional implicit steepest descent method, we put forward the alternating minimization algorithm for obtaining its numerical solution.In this paper, we have two principal contributions. For one thing, a fourth-order regularization strategy for restoring Poissonian image is proposed, which substantially improves the quality of the reconstructed images, especially in relieving staircase artifacts. For another, we introduce the alternating minimization (AM) algorithm for solving our proposed fourth-order scheme (1.4).The remainder of this paper is organized as follows. Section2 shows the existence and uniqueness of the solution of our proposed model (1.4). Section 3 focuses on the numerical algorithm for solving the proposed strategy. The convergence analysis is arranged in Section 4. Results of numerical experiments to indicate the advantage of our approach in avoiding blocky effects and enhancing convergence speed, compared with the classical TV regularization based one, are shown in Section 5. Finally, we make the conclusion in Section 6.


## 2. Existence and Uniqueness

In this section, the existence and uniqueness of the minimizer for the model (1.4) are proved. Some basic notations and properties concerning with the BV2 space can be seen in [8, 13]. Motivated by [5], we have the following existence and uniqueness result for the optimization problem (1.4).Theorem 2.1.
AssumeΩ is an open, bounded Lipschitz domain in ℝn and f is a positive bounded function. Let u∈BV2(Ω) such that logu∈L1(Ω). Then there exists a unique minimizer for model (1.4) in BV2(Ω).Proof.
Obviously, the objective functionalE(u) is lower bounded. Then there exists a minimizing sequence {ui}i=1∞. Let F(u)=∫Ω(u-flogu)dx. Then, using Jensen inequality, we conclude that
(2.1)‖ui‖1-‖f‖∞log⁡⁡‖ui‖1≤F(ui),
which indicates that ∥ui∥1 is bounded. This, together with the boundedness of {∫Ω|D2ui|}, yields that {ui}i=1∞ is bounded in BV2(Ω). Then, following the compactness theorem in BV2(Ω) space, we deduce that there exists a function u*∈BV2(Ω) such that a subsequence of {ui} denoted also by {ui} converges to u* a.e. in L1(Ω). In addition, by the lower semicontinuity of the BV2(Ω) space, we have liminfi→∞∥D2ui∥1≥∥D2u*∥1. Combining this with Fatou's Lemma, we obtain
(2.2)∫Ω|D2u*|+λ∫Ω(u*-flog⁡u*)dx≤liminf⁡i→∞{∫Ω|D2ui|+λ∫Ω(ui-flog⁡ui)dx},
which ensures that u* is the minimizer of (1.4). Moreover, notice that the optimization problem (1.4) is strictly convex, then u* is unique by convex analysis [16]. This completes the proof.


## 3. Computational Method

This section focuses on minimizing the energy functionalE(u) of the model (1.4). Generally speaking, minimizing (1.4) can be achieved by the traditional gradient descent method [1, 17], which amounts to solve its formal Euler-Lagrange equation in W2,1(Ω):(3.1)(uxx|∇2u|)xx+(uxy|∇2u|)yx+(uyx|∇2u|)xy+(uyy|∇2u|)yy+λu(u-f)=0,
where |∇2u|=(uxx2+uxy2+uyx2+uyy2)1/2. Introducing the auxiliary time variable t, we derive the associated heat flow(3.2)∂u∂t=-(uxx|∇2u|)xx-(uxy|∇2u|)yx-(uyx|∇2u|)xy-(uyy|∇2u|)yy+λu(f-u).
Here, we omit the boundary conditions, but it can be deduced easily. Finite difference approaches [1, 18–20] can be chosen to resolve problem (3.2). As everyone knows, |∇2u| is usually replaced by |∇2u|2+ɛ because of the nondifferentiability of |∇2u|, where ɛ is a small positive parameter. The tradeoff in choosing the smoothing parameter ɛ is the reconstruction error versus the convergence speed. Therefore, exploring fast numerical computational methods is an important assignment. In the last few years, many effective algorithms have sprung up for handling problems like this. However, a great many current algorithms cannot be directly applied to problem (1.4) because of its nonquadratic and nonsmooth higher order properties, for instance, dual algorithm [20–22] and split Bregman iteration [23].Motivated by [24–26], we adopt the alternating minimization (AM) algorithm to acquire the optimal solution of the scheme (1.4). Take the numeric span of the variable u of (1.4) into consideration, we substitute the variable z for logu. So instead of the model (1.4), we consider the following minimization problem:(3.3)min⁡z∫Ω|D2z|+λ∫Ω(ez-fz)dx.
Owing to the calculation difficulties brought by the non-differentiability term ∫Ω|D2z|, we introduce an auxiliary variable w to add a quadratic fitting term ∥w-z∥22 to the model (3.3). Then the following approximation scheme for problem (3.3)(3.4)min⁡z,wE(z,w)=min⁡z,w∫Ω(ez-fz)dx+α‖w-z‖22+β∫Ω|D2w|
is achieved. Here, α and β are two positive parameters. Usually, α is larger than β, since relatively small α will result in inferior restoration quality. While for many problems, too large value for α makes (3.4) extremely difficult to solve numerically. Notice that the objective function of problem (3.4) is strictly convex and differential with respect to the first variable z, then the approach (3.4) can be decomposed into two decoupled z,w subproblems, namely,(3.5)zk=argmin⁡z∫Ω(ez-gz)dx+α‖z-wk-1‖22,wk=argmin⁡wα‖zk-w‖22+β∫Ω|D2w|.For thez subproblem, the objective function is strictly convex, then the solution is unique. In view of the differentiability of the objective function with regard to z, to determine the solution of the first equation of the system (3.5) amounts to solve(3.6)ez-g+2α(z-wk-1)=0.
Obviously, Newton iteration method is a good choice for efficiently solving these kinds of problems.As for the second subproblem, we make use of dual algorithm, which was first introduced by Chen et al. [20] for TV regularization and then was extended to the fourth-order situation for solving model (1.3) [21]. The dual formulas for the w subproblem in the discrete setting can be briefly depicted as(3.7)pi,jk+1=pi,jk-τ(∇2(div2pn-(β/2α)zk))i,j1+τ|∇2(div2pn-(β/2α)zk)i,j|.
Subsequently, for all time steps  τ with 0<τ≤1/64, as n→∞, we get the minimizer of the w subproblem(3.8)wk=zk-2αβdiv2p.Then, given the initial valuew0, based on the above-mentioned AM approach, we obtain the sequence given below(3.9)z1,w1,z2,w2,…,zk,wk,….
Thanks to the computationally inexpensive of both the subproblems w and z, the AM approach occupies advantage in convergence speed so that the calculation time cost is cut down to a great extent.


## 4. Convergence Analysis

Letzk=R(wk-1), wk=S(zk), and P=S∘R. Then we have wk=P(wk-1). By Opial theorem [27], firstly, we should show that the operator P is nonexpansive and asymptotically regular. Recall that for a proper, convex, lower semicontinuous function φ, the operator proxφ:Rn2→Rn2:(4.1)proxφ(x)=argmin⁡y12‖x-y‖22+φ(y)
is called the proximal operator of φ. According to the definition of proximal operator, it is easy to check that the operators R and S are proximal operators of (1/2α)∫Ω(ez-gz)dx and (β/2α)∫Ω|D2w|, respectively. Combining this with Lemma 2.4 [28], we obtain that both R and S are nonexpansive. Then P is nonexpansive.Lemma 4.1.
The operatorP is asymptotically regular.Proof.
For an initial guess,w0∈Rn2. To show, P is asymptotically regular, by the definition, we have to show that the sequence {Pk+1(w0)-Pk(w0)} tends to zero as k→∞, that is, limk→∞∥wk+1-wk∥2=0. Let H1(z,w)=∥z-w∥22 and H2(w)=∫Ω|D2w|. Then,
(4.2)E(zk+1,wk)-E(zk+1,wk+1)=α[H1(zk+1,wk)-H1(zk+1,wk+1)]+β[H2(wk)-H2(wk+1)].
By the Taylor expansion of H1(z,w) in the second variable and the fact that ∂2H1/∂w2=2I, we have
(4.3)H1(zk+1,wk)-H1(zk+1,wk+1)=(wk-wk+1)t∂H1∂w(zk+1,wk+1)+‖wk+1-wk‖22.
Further, considering the convexity property of H2, we obtain
(4.4)H2(wk)-H2(wk+1)≥(wk-wk+1)t∂H2∂w(wk+1).
(4.2), (4.3), together with (4.4), deduce
(4.5)E(zk+1,wk)-E(zk+1,wk+1)≥α‖wk+1-wk‖22.
Notice that α(∂H1/∂w)(zk+1,wk+1)+β(∂H2/∂w)(wk+1)=0, as wk+1 is the minimizer of E(zk+1,w). In addition, combining (4.5) with the fact that E(zk+1,wk+1)≥E(zk+2,wk+1), we have
(4.6)E(zk+1,wk)-E(zk+2,wk+1)≥E(zk+1,wk)-E(zk+1,wk+1)≥α‖wk+1-wk‖22.
Therefore,
(4.7)∑k=0∞‖wk+1-wk‖22≤1αE(z1,w0),
which guarantees that ∥wk+1-wk∥2 tends to zero as k→∞.Lemma 4.2.
The set of fixed points ofP is nonempty.Proof.
First of all, similar to [25], we can derive that the function E is coercive. Moreover, as E is closed and proper, by the Weierstrass's Theorem [16], we know that the set of minimizers of E is nonempty. Suppose (z,w) is a minimizer of E(z,w), then we have
(4.8)∂E∂z(z,w)=0,∂E∂w(z,w)=0.
It follows that
(4.9)S(z)=w=argmin⁡wE(z,⋅),(w)=z=argmin⁡zE(⋅,w).
Therefore, w=S(R(w))=P(w), which indicates that w is a fixed point of P.Then, by Opial theorem [27], we conclude that {wk} converges to a fixed point of P. As E is differentiable with regard to the second variable z and E is strictly convex, then the unique minimizer of E is just the fixed point of P. Hence, {wk} converges to the unique minimizer of E.


## 5. Numerical Experiments

In this section, we pay attention to illustrate the efficiency and superiority of our proposed scheme. We compare the TV-regularization-based model (1.2) by gradient descent method with the fourth-order (1.4) by gradient descent method and the fourth-order scheme (3.4) by alternating minimization method, that is, Algorithm 1, in measuring the convergence speed and image restoration quality. Let GD-(1.2), GD-(1.4), and AM-(3.4) denote the above three schemes, respectively. The performance of the three schemes will be tested on three images: “lenna”, “pepper”, and “panda”. All of our experiments are processed under Windows XP and MATLAB R2009a. The Poisson noise is added by using the MATLAB library function imnoise (I, “poisson”). It is worth to mention that the poisson noise magnitude in an image increases with the intensity of the region of the image, which is different from Gaussian noise that is signal independent. We choose time step τ=1/64 for all the following experiments.For well assessing the quality of recovered images, the relative error (ReErr) [25] is adopted as follows:(5.1)ReErr=‖u-f‖2‖f‖2,
where u and f represent reconstructed image and the original image, respectively. Generally speaking, the smaller the ReErr value, the better the restored image quality.Example 5.1.
The first test image “lenna” (Figure1(a)) is a gray image of size 256×256 pixels. The corresponding degraded image corrupted by Poisson noise is demonstrated in Figure 1(b). Figures 1(c), 1(d), and 1(e) denote the restoration images via 460 iterations by GD-(1.2), 280 iterations by GD-(1.4), and 45 iterations by AM-(3.4), respectively. Their corresponding enlarged partial images are displayed in Figures 2(a)–2(e). Table 1 lists the number of iterations, computational time, and the value of ReErr of the three methods about this experiment. From Table 1, we can see that the AM-(3.4) scheme can achieve lower value of ReErr with less time than GD-(1.2) and GD-(1.4). Comparing GD-(1.2) with GD-(1.4), although GD-(1.2) costs less time per iteration, GD-(1.4) is slightly fast as a whole. In terms of the restoration result, the results images processed by the fourth-order approaches GD-(1.4) (see Figure 2(d)) and AM-(3.4) (see Figure 2(e)) are similar, but the one (see Figure 2(c)) processed by second-order approach GD-(1.2) suffers from disgusting staircasing effects just as the ROF model for Gaussian noise removal. It is turned out that fourth-order approaches successfully overcome the piecewise constant of the solution and can obtain relative natural recovered results than the second-order one. The parameters with respect to the above test are λ=0.01,α=100,andβ=1.6.Table 1
Number of iterations, computational time, and ReErr of the three methods for the results in Figures1 and 2.
LennaIterationsRun time (second)ReErrGD-(1.2)451.840.0785GD-(1.4)452.460.0725GD-(1.2)46017.500.0562GD-(1.4)28013.230.0543AM-(3.4)451.870.0505Recovered results of different models by different numerical algorithms.
(a)
Original image(b)
Degenerated image(c)
460 iterations by GD-(1.2)(d)
280 iterations by GD-(1.4)(e)
45 iterations by AM-(3.4)Enlarged partial recovered results by different approaches.
(a)
Original image(b)
Degenerated image(c)
460 iterations by GD-(1.2)(d)
280 iterations by GD-(1.4)(e)
45 iterations by AM-(3.4)Example 5.2.
To further evaluate the performance of our proposed fourth-order scheme by AM method, we test the three methods on the512×512 image “pepper” (Figure 3(b)). The corresponding parameters used in this experiment are set as follows: λ=0.01,α=100,andβ=1.24. From Figure 5(c), we can see that the staircasing effects still appears in the homogeneous flat regions of restoration image by second-order scheme (1.2). Figures 5(d) and 5(e) denote the results processed by fourth-order PDE, which are more close to the original one. The number of iterations, computational time, and value of ReErr of the three methods about this experiment are listed in Table 2, and the relations between ReErr and the number of iterations are illustrated in Figure 4. Figure 4 indicates that the AM algorithm for (3.4) converges so fast that 50 iterations are enough to obtain relative ideal result with ReErr=0.0410, but the numbers of iterations needed by GD-(1.2) (see Figure 5(a)) and GD-(1.4) (see Figure 5(b)) are about 700 and 350, respectively. As a result, our proposed scheme, that is, AM-(3.4), can obtain higher quality restoration images with less time, especially for some larger images.Table 2
Number of iterations, computational time, and ReErr of the three methods for the results in Figures3 and 5.
PepperIterationsRun time (second)ReErrGD-(1.2)700150.150.0420GD-(1.4)350113.480.0426AM-(3.4)5515.960.0410Recovered results of different models by different numerical algorithms.
(a)
Original image(b)
Degenerated(c)
700 iterations by GD-(1.2)(d)
350 iterations by GD-(1.4)(e)
55 iterations by AM-(3.4)Relations between number of iterations (abscissa axis) and value ofReErr (vertical axis) by the three methods for “pepper”.
(a)
GD-(1.2)(b)
GD-(1.4)(c)
AM-(3.4)Recovered results by the three methods.
(a)
Original image(b)
Degenerated image(c)
700 iterations by GD-(1.2)(d)
350 iterations by GD-(1.4)(e)
55 iterations by AM-(3.4)Example 5.3.
We use the256×256 image “panda” (Figure 6(a)) as our third test image. Unlike “lenna” and “pepper”, the image “panda” has more small details. This experiment will show that our proposed scheme also performs well on the images with small details. The degenerated “panda” corrupted by Poisson noise is displayed in Figure 6(b). Parameters are chosen to be λ=0.01,α=110,andβ=1.9. Table 3 lists the number of iterations, computational time, and ReErr of the three methods about this experiment. The relation between value of ReErr and the number of iterations can be seen in Figure 7. From Table 3 and Figure 7, we can see that fourth-order schemes can achieve higher quality images, evaluated by lower value of ReErr, than the second one. Moreover, our proposed AM-(3.4) outperform the other both in reducing numerical computational time and image restoration quality.Table 3
Number of iterations, computational time, and ReErr of the three methods for the results in Figure6.
PandaIterationsRun time (second)ReErrGD-(1.2)45017.210.0541GD-(1.4)25010.040.0530AM-(3.4)602.680.0496Recovered results by the three methods.
(a)
Original image(b)
Degenerated image(c)
450 iterations by GD-(1.2)(d)
250 iterations by GD-(1.4)(e)
60 iterations by AM-(3.4)Relations between number of iterations (abscissa axis) and value ofReErr (vertical axis) by the three methods for “panda”.
(a)
GD-(1.2)(b)
GD-(1.4)(c)
AM-(3.4)


## 6. Conclusion

This paper investigates a fourth-order regularization based scheme for restoring corrupted images by Poisson noise. The existence and uniqueness of the minimizer for the proposed model are presented in detail. To quickly obtain the optimum solution of our proposed scheme, the alternating minimization method is employed. Compared with the recovered results by the TV-regularization-based one, simulation results evidently demonstrate the competitive superiority of our proposed approach especially in overcoming staircase effects and improving the computational speed while removing Poisson noise.Algorithm 1:Alternating minimization algorithm for solving the model (3.4).
Initialization: w0=0 and p0=0;Whilek=1,2,3,…N.Computezk according to the first equation of (3.5) for fixed wk-1.Computewk according to the second equation of (3.5) for fixed zk.EndComputeu=exp(w).



---

*Source: 965281-2012-02-27.xml*