A Convex Adaptive Total Variation Model Based on the Gray Level Indicator for Multiplicative Noise Removal

Abstract

This paper focuses on the problem of multiplicative noise removal. Using a gray level indicator, we derive a new functional which consists of the adaptive total variation term and the global convex fidelity term. We prove the existence, uniqueness, and comparison principle of the minimizer for the variational problem. The existence, uniqueness, and long-time behavior of the associated evolution equation are established. Finally, experimental results illustrate the effectiveness of the model in multiplicative noise reduction. Different from the other methods, the parameters in the proposed algorithms are found dynamically.

1. Introduction

Multiplicative noise occurs while one deals with active imaging system, such as laser images, microscope images, and SAR images. Given a noisy image 𝑓∶Ω→𝑅, where Ω is a bounded open subset of ℝ², we assume 𝑓=𝑢𝜂, (1) where 𝑢 is the true image and 𝜂 is the noise. In what follows we will always assume that 𝑓>0 and 𝑢>0. Some prior information about the mean and variance of the multiplicative noise is as follows: 1|Ω|_Ω1𝜂𝑑𝑥=1,|Ω|_Ω𝜂−1²𝑑𝑥=𝜎², (2)(3) where |Ω|=∫_Ω1𝑑𝑥. Our purpose is to remove noise in images while preserving the maximum information about 𝑢.

The goal of this paper is to propose a globally strictly convex functional well-adapted to removing multiplicative noise, which is as follows: 𝐸(𝑢)=_Ω𝛼(𝑥)|𝐷𝑢|+_Ω1𝑢+𝑓log𝑢𝑑𝑥, (4) where 𝑓 is noisy image and ∫_Ω𝛼(𝑥)|𝐷𝑢| stands for the adaptive total variation of 𝑢. In the new model, we introduce a control factor, 𝛼(𝑥), which controls the speed of diffusion at different region: at the low gray level (𝛼(𝑥)→0), the speed is slow; at the high gray level (𝛼(𝑥)→1), the speed is fast. The second term in the functional, that is, the fidelity term is global convex, which implies the constraint (2).

Various adaptive filters for multiplicative noise removal have been proposed. In the beginning, variational methods for multiplicative noise reduction deal with the Gaussian multiplicative noise [1]. In actual life, the speckle noise [2] is more widespread, such as synthetic aperture radar (SAR) imagery [3, 4]. Then, Aubert and Aujol [5] propose a new variational method for Gamma multiplicative noise reduction. By the logarithm transform, a large variety of methods rely on the conversion of the multiplicative noise into additive noise.

1.1. Gaussian Noise

In the additive noise case, the most classical assumption is that the noise is a white Gaussian noise. So one case in which dealing with multiplicative noise is the white Gaussian noise. Using the framework in [6], Rudin et al. [1] consider the following optimization problem for Gaussian multiplicative noise reduction: 𝑢=argmin_{𝑢∈BV(Ω)}𝐽(𝑢)+𝜆𝐻𝑢,𝑓, (5) where 𝐽(𝑢)=∫_Ω|∇𝑢| stands for the total variation of 𝑢, and 𝐻(𝑢,𝑓) is a fidelity term which consists of two integrals with two Lagrange multipliers: 𝐻𝑢,𝑓=𝜆_1Ω𝑓𝑢𝑑𝑥+𝜆_2Ω𝑓𝑢−1²𝑑𝑥. (6)

In order to make sure the two constraints (2) and (3) are always satisfied during the evolution, by the gradient projection method, the authors evolve the following evolution equation: 𝜕𝑢𝜕𝑡=div∇𝑢|∇𝑢|+𝜆₁𝑓𝑢²+𝜆₂𝑓²𝑢³. (7)

If by the gradient projection method the values of 𝜆₁ and 𝜆₂ are found dynamically, the method is not always convex; while if 𝜆₁,𝜆₂>0 are fixed, the corresponding minimization problem will lead to a sequence of constant function 𝑢 approaching +∞.

1.2. Gamma Noise

Generally, the speckle noise is treated as Gamma noise with mean equal to one. The probability distribution function 𝑔(𝜂) of the noise 𝜂 takes the following form: 𝑔𝜂=𝐿^𝐿Γ(𝐿)𝜂^𝐿−1exp−𝐿𝜂1_{𝜂≥0}. (8)

Gamma noise is more complex than Gaussian noise [2]. Based on maximum a posteriori (MAP) on 𝑝(𝑢∣𝑓), Aubert and Aujol [5] assume that the noise 𝜂 follows a Gamma probability distribution with mean equal to one and 𝑝(𝑢) following a Gibbs prior and then derive a functional formulating the following minimization problem (called AA model): 𝑢=argmin_{𝑢∈𝑆(Ω)Ω}|𝐷𝑢|+𝜆_Ω𝑓log𝑢+𝑢𝑑𝑥, (9) where 𝑆(Ω)={𝑢>0,𝑢∈BV(Ω)}. The new fidelity term 𝐻(𝑢,𝑓)=∫_Ω(log𝑢+(𝑓/𝑢))𝑑𝑥 is strictly convex for 𝑢∈(0,2𝑓). The authors prove the existence of a minimizer 𝑢∈𝑆(Ω) to the minimization problem and derive existence and uniqueness results of the solution to the associated evolution equation: 𝜕𝑢𝜕𝑡=div∇𝑢|∇𝑢|+𝜆𝑓−𝑢𝑢². (10)

1.3. The Model Based on the Logarithm Transform

The simplest idea is to take the log of both sides of (1), log𝑓=log𝑢+log𝜂=𝑤+𝑛, (11) which essentially converts the multiplicative problem into an additive one. If the distribution of the noise 𝜂 takes the form (8), the expectation and the variance of the log-noise 𝑛 are 𝔼[𝑛]=Ψ₀(𝐿)−Ψ_𝑘(𝐿),Var[𝑛]=Ψ₁(𝐿), (12) where Ψ_𝑘(𝑧)=𝑑𝑑𝑧^𝑘+1logΓ(𝑧) (13) is the polygamma function [7]. In [8], Shi and Osher use relaxed inverse scale space (RISS) flows to deal with various noises and provide iterative TV regularization. First, using the log-data log𝑓, they propose the following model: 𝑢=argmin_{𝑤∈BV(Ω)Ω}||+𝜆𝐷𝑤2𝑤−log𝑓²_𝐿2(Ω)𝑑𝑥, (14) where 𝑤=log𝑢. The corresponding RISS flow reads 𝜕_𝑡𝑤=div∇𝑤||,𝜕∇𝑤+𝜆log𝑓−𝑤+𝑣_𝑡𝑣=𝛼log𝑓−𝑤, (15) with 𝑣(0)=0, 𝑤(0)=𝑐₀, for 𝑐₀=∫_Ωlog𝑓𝑑𝑥/|Ω|. Then, they generalize several multiplicative noise models [5, 9, 10], and for convergence reasons, use iterative TV regularization on exp(𝑤) to obtain the RISS flow. And they investigate the properties of the flow and study the dependence on flow parameters.

1.4. The Model Based on the Exponent Transform

Recently, Huang et al. [11] utilize an exponential transformation 𝑢→𝑒^𝑢 in the fidelity term of AA model to propose the following denoising model: min_𝑢,𝑤Ω𝑢+𝑓𝑒^−𝑢𝑑𝑥+𝛼_1Ω|𝑢−𝑤|²𝑑𝑥+𝛼_2Ω|𝐷𝑤|, (16) where 𝛼₁>0 and 𝛼₂>0 are regularization parameters, and 𝑤 is an auxiliary variable. Next, they further develop an alternating minimization algorithm for the model (16) by incorporating another way of modified TV regularization in [12]. However, the mathematical analysis to the variational problem (16) is not given in [11]. By the exponential transformation 𝑢→𝑒^𝑢, Jin and Yang [13] change the fidelity term log𝑢+𝑓𝑢 of AA model (2) into 𝑢+𝑓𝑒^𝑢 and then study the following denoising model: min_𝑢Ω|𝐷𝑢|+𝜆_Ω𝑢+𝑓𝑒^−𝑢𝑑𝑥. (17)

Notice that the fidelity term 𝐻(𝑢,𝑓)=∫_Ω(𝑢+𝑓𝑒^−𝑢)𝑑𝑥 is globally strictly convex. Based on this, they prove the uniqueness of the solutions to the variational problem (17) and show the existence and uniqueness of the weak solution for the following evolution equation corresponding to (17): 𝜕𝑢𝜕𝑡=div∇𝑢|∇𝑢|+𝜆𝑓𝑒^−𝑢−1. (18)

The paper is organized as follows. In Section 2, inspired from [1, 5, 14], based on the gray level indicator 𝛼(𝑥), we derive a convex adaptive total variation model (4) for multiplicative noise removal. In Section 3, we prove the existence and uniqueness of the solution of the minimization problem (4). In Section 4, we study the associated evolution problem to the minimization problem (4). Specifically, we define the weak solution of the evolution problem, derive estimates for the solution of an approximating problem, prove existence and uniqueness of the weak solution, and discuss the behavior of the weak solution as 𝑡→∞. In Section 5, we provide two numerical algorithms and experimental results to illustrate the effectiveness of our algorithms in image denoising. We also compare new algorithms with other ones.

2. A New Variational Model for Multiplicative Noise Removal

The goal of this section is to propose a new variational model for multiplicative noise removal. First, we proposed a new fidelity with global convexity, which can always satisfy the constraint (2) during the evolution. Then, by analyzing the properties of the noise, we propose the adaptive total variation based on a gray level indicator 𝛼(𝑥).

2.1. A Global Convex Fidelity Term

Based on the idea in [15], we can obtain the following fidelity: 𝐻1𝑢,𝑓=𝑢+𝑓log𝑢. (19)

Note that 𝐻^′(𝑢,𝑓)=1−(𝑓/𝑢). Let us consider the following Euler-Lagrange equation for some functional-based problem: 0=ℱ∇𝑢,∇²𝑢𝑓−𝜆1−𝑢,inΩ,𝜕𝑢𝜕𝑛=0,on𝜕Ω, (20)(21) where ∇𝑢 and ∇²𝑢 stand, respectively, for the gradient and the Hessian matrix of 𝑢 with respect to the space variable 𝑥, and ℱ(∇𝑢,∇²𝑢) is divergence operator from the functional. By integrating (20) in space, using integration by parts and the boundary condition (21) in the sense of distributions, we have _Ω𝑓1−𝑢𝑑𝑥=0, (22) which implies the constraint (2). Moreover, using the idea in [6], the parameter 𝜆 can be calculated as 1𝜆=𝜎²|Ω|_Ωℱ∇𝑢,∇²𝑢𝑓1−𝑢𝑑𝑥. (23)

Remark 1. (1) It is easy to check that the function 𝑢→𝑢+𝑓log(1/𝑢) reaches its minimum value 𝑓+𝑓log(1/𝑓) over ℝ⁺ for 𝑢=𝑓.
(2) Let us denote by ℎ(𝑢)1=𝑢+𝑓log𝑢. (24)
We have ℎ^′(𝑢)=1−(𝑓/𝑢)=(𝑢−𝑓)/𝑢, and ℎ^′(𝑢)=(𝑓/𝑢²)>0, which deduce that the new fidelity term 𝐻(𝑢,𝑓) is globally strictly convex.

2.2. The Adaptive Total Variation Model

Assume 𝑢 is a piecewise function, that is, 𝑢=∑^𝑁_𝑖=1𝑔_𝑖1_Ω𝑖, where Ω_𝑖∩Ω_𝑗=∅, for 𝑖≠𝑗, ⋃^𝑁_𝑖=1Ω_𝑖=Ω, and 𝑔_𝑖 is the gray level. Moreover, we assume that the samples of the noise on each pixel 𝑥 are mutually independent and identically distributed (i.i.d.) with the probability density function 𝜂(𝑥). For 𝑥∈Ω, 𝑓(𝑥)=𝑔_𝑖𝜂(𝑥), and therefore, Var[𝑓]=𝑔²_𝑖Var[𝜂]=𝑔²_𝑖𝜎², where Var[𝑓] and Var[𝜂] are the variance of the noise image 𝑓 and the noise at the pixel 𝑥, respectively. It is noticed that the variance of the noise is the constant 𝜎² on each pixel, but the variance of the noise image 𝑓 is influenced by the gray level. The higher the gray level is, the more remarkable the influence of the noise is. Especially, 𝑓=𝑢 when 𝑢=0, and therefore 𝑓 is noise free in this case. The fact is illustrious in Figure 1. It can be seen that in despite of the independent identically distribution of noise (see Figure 1(b)), the noise image shows different features on the pixel, where the gray levels are different (see Figure 1(d)).

(a) Noise-free signal

(b) Speckled noise

(c) Noise signal

(d) Original/noise signal

Figure 1 

The relation between the influence of the noise and the gray levels. (a) 1D signal 𝑓; (b) 1D speckled noise with mean equal to 1 and 𝐿=5; (c) 𝑓 degraded by the speckled noise; (d) the compare between 𝑓 and the noise signal 𝑓.

In [14], Chan and Strong propose the following adaptive total variation model: min_Ω𝑔(𝑥)|∇𝑢|𝑑𝑥, (25) where the weight function 𝑔(𝑥) controls the speed of diffusion at different region. Utilizing this idea, we proposed a gray level indicator 𝛼(𝑥). The indicator 𝛼(𝑥) has such properties as follows: 𝛼(𝑠) is monotonically increasing, 𝑔(0)=0, 𝛼(𝑠)≥0, and 𝛼(𝑠)→1, as 𝑠→sup_𝑥∈Ω𝑢. Therefore, we propose the following gray indicators 𝛼(𝑥): 𝛼(𝑥)=11−1+𝑘𝐺_𝜎∗𝑓²1+𝑘𝑀²𝑘𝑀², (26) or 𝛼(𝑥)=𝐺_𝜎∗𝑓𝑀, (27) where 𝑀=sup_𝑥∈Ω(𝐺_𝜎∗𝑓)(𝑥), 𝐺_𝜎(𝑥)=(1/4𝜋𝜎)exp{−|𝑥|²/4𝜎²}, 𝜎>0, and 𝑘>0 are parameters. With this choice, 𝛼(𝑥) is a positive-valued continuous function; 𝛼(𝑥) is much smaller value at low gray levels (𝛼(𝑥)→0) than at high gray levels (𝛼(𝑥)→1) so that some small features at low gray levels are much less smoothed and therefore are preserved. As previously stated, at high gray levels, the region of the noise image is degraded more, while the region at low gray levels is degraded less (see Figure 1(d)). Then, from (26)/(27), at the high gray levels, 𝛼(𝑥)→1, the new model is more smooth at the region; at low gray levels, 𝛼(𝑥)→0, the new model is less smooth at the region.

The previous analysis leads us to propose a convex adaptive total variation model for multiplicative noise removal, min_{𝑢∈BV(Ω)Ω}𝛼(𝑥)|𝐷𝑢|+𝜆_Ω1𝑢+𝑓log𝑢𝑑𝑥. (28)

The evolution of the Euler-Lagrange equation for (28) is as follows: 𝜕𝑢𝛼(𝑥)𝜕𝑡=div∇𝑢||𝑓∇𝑢−𝜆1−𝑢,(𝑥,𝑡)∈Ω×(0,𝑇),𝜕𝑢𝜕𝑛()[],𝑢()=0,𝑥,𝑡∈𝜕Ω×0,𝑇0,𝑥=𝑓,𝑥∈Ω. (29)(30)(31)

3. The Minimization Problem (28)

In this section, we study the existence and uniqueness of the solution to the minimization problem (28), and then we consider the comparison principle for the problem (28).

3.1. Preliminaries

If 𝑓∈𝐿^∞(Ω) with 0<inf_𝑥∈Ω𝑓≤𝑓≤sup_𝑥∈Ω𝑓, then 0<𝛼inf_𝑥∈Ω𝑓≤𝛼(𝑥)≤1,𝑥∈Ω. (32)

We always assume that Ω is a bounded open subset of ℝ^𝑛 with Lipschitz boundary. As in [16], we introduce the following definition of 𝛼-total variation.

Definition 2. A function 𝑢∈𝐿¹(Ω) has bounded 𝛼-total variation in Ω, if sup_{𝜙∈Φ𝛼Ω}𝑓div𝜙𝑑𝑥<∞, (33) where Φ_𝛼=∶𝜙∈𝐶¹₀Ω,ℝ^𝑛,𝜙≤𝛼. (34)

Remark 3 (see [16]). (1) If 𝑢∈𝐿¹(Ω) has bounded 𝛼-total variation in Ω, there is a Radon vector measure 𝐷𝑢 on Ω such that _Ω𝛼|𝐷𝑢|=∶sup_{𝜙∈Φ𝛼Ω}𝑓div𝜙𝑑𝑥. (35)
(2) From (32), 𝑢∈𝐿¹(Ω) having bounded 𝛼-total variation in Ω implies that 𝑢∈BV(Ω).
Now, we directly show some lemmas on 𝛼-total variation from [16].

Lemma 4. Assume that {𝑢_𝑘}^∞_𝑘=1⊂BV(Ω) and 𝑢_𝑘→𝑢 in 𝐿¹(Ω). Then, 𝑢∈BV(Ω), and _Ω𝛼𝐷𝑢_𝑘≤liminf_𝑘→∞Ω𝛼𝐷𝑢_𝑘. (36)

Lemma 5. Assume that 𝑢∈BV(Ω). Then, there exists a sequence {𝑢_𝑗}^∞_𝑗=1∈𝐻¹(Ω) such that 𝑢_𝑗⟶𝑢in𝐿¹(Ω),_Ω𝛼𝐷𝑢_𝑗⟶_Ω𝛼|𝐷𝑢|. (37)

By a minor modification of the proof of Lemma 1 in Section 4.3 of [17], we can have the following lemma.

Lemma 6. Let 𝑢∈BV(Ω), and let 𝜑_𝑎,𝑏 be the cut-off function defined by 𝜑_𝑎,𝑏(𝑥)=𝑎,if𝑥≤𝑎,𝑥,if𝑎≤𝑥≤𝑏,𝑏,if𝑥≥𝑏. (38) Then, 𝐷𝜑_𝑎,𝑏(𝑢)∈BV(Ω), and _Ω𝛼𝐷𝜑_𝑎,𝑏(𝑢)≤_Ω𝛼|∇𝑢|. (39)

3.2. Existence and Uniqueness of the Problem (28)

In this subsection, we show that problem (28) has at least one solution in BV(Ω).

Theorem 7. Let 𝑓∈𝐿^∞(Ω) with inf_𝑥∈Ω𝑓>0. Then, the problem (28) admits a unique solution 𝑢∈𝐵𝑉(Ω) such that 0<inf_𝑥∈Ω𝑓≤𝑢≤sup_𝑥∈Ω𝑓. (40)

Proof. Let us rewrite 𝐸(𝑢)=_Ω𝛼(𝑥)|𝐷𝑢|+𝜆_Ω1𝑢+𝑓log𝑢𝑑𝑥. (41) From Remark 1(1), we have 𝐸(𝑢)≥𝜆_Ω1𝑓+𝑓log𝑓𝑑𝑥, (42) for 𝑢∈BV(Ω) with 𝑢>0. This implies that 𝐸(𝑢) has a lower bound for all 𝑢∈BV(Ω) with 𝑢>0. Hence, there exists a minimizing sequence {𝑢_𝑛}⊂BV(Ω) for the problem (28).
Step 1. Let us denote 𝑎=inf_𝑥∈Ω𝑓 and 𝑏=sup_𝑥∈Ω𝑓. We first claim that 0<𝑎≤𝑢_𝑛≤𝑏.
In fact, we remark that ℎ(𝑠) is decreasing for 𝑠∈(0,𝑓) and increasing for 𝑠∈(𝑓,+∞). Therefore, if 𝑀≥𝑓, one always has ((+𝑓logmin𝑠,𝑀))()𝑓min𝑠,𝑀≤log𝑠+𝑠. (43) Hence, if 𝑀=𝑏=sup_𝑥∈Ω𝑓, we have _Ω((+𝑓loginf𝑢,𝑀))inf(𝑢,𝑀)𝑑𝑥≤_Ω𝑓log𝑢+𝑢𝑑𝑥. (44) Moreover, utilizing Lemma 6, we can have _Ω𝛼|𝐷(inf(𝑢,𝑏))|≤_Ω𝛼|𝐷𝑢|. (45) Combining (44) and (45), we deduce that 𝐸(inf(𝑢,𝑏))≤𝐸(𝑢). (46) On the other hand, in the same way, we get that 𝐸(sup(𝑢,𝑎))≤𝐸(𝑢). Therefore, we can assume that 𝑎≤𝑢_𝑛≤𝑏 without restriction.
Step 2. Let us prove that there exists 𝑢∈BV(Ω) such that 𝐸(𝑢)=min_{𝑣∈BV(Ω)}𝐸(𝑣). (47) The proof in the Step 1 implies in particular that 𝑢_𝑛 is bounded in 𝐿¹(Ω). Since 𝑎≤𝑢_𝑛≤𝑏 and ℎ(𝑠)∈𝐶[𝑎,𝑏], ℎ(𝑢_𝑛) is bounded. Moreover, by the definition of {𝑢_𝑛}, there exists a constant 𝐶 such that _Ω𝛼𝐷𝑢_𝑛+_Ωℎ𝑢_𝑛𝑑𝑥≤𝐶. (48) Then, _Ω𝛼𝐷𝑢_𝑛≤𝐶, (49) which implies that the sequence {𝑢_𝑛}^∞_𝑛=1 is bounded in BV(Ω). Consequently, there exists a function 𝑢∈BV(Ω) and a subsequence {𝑢_𝑛}^∞_𝑛=1, denoted by itself, such that as 𝑛→∞, 𝑢_𝑛⇀𝑢inweak^∗(Ω),𝑢inBV_𝑛⟶𝑢stronglyin𝐿¹(Ω),𝑎≤𝑢≤𝑏. (50) Utilizing Lemma 4 and Fatou’s lemma, we get that 𝑢 is a solution of the problem (28).
Finally, from Remark 1(2), ℎ is strictly convex as 𝑓>0, and then the uniqueness of the minimizer follows from the strict convexity of the energy functional in (28).

3.3. Comparison Principle

In this subsection, we state a comparison principle for problem (28).

Proposition 8. Let 𝑓₁,  𝑓₂∈𝐿^∞(Ω) with inf_𝑥∈Ω𝑓₁,  inf_𝑥∈Ω𝑓₂>0, and 𝑓₁<𝑓₂. Assume that 𝑢₁ (resp., 𝑢₂) is a solution of the problem (28) for 𝑓=𝑓₁ (resp., 𝑓=𝑓₂). Then, 𝑢₁≤𝑢₂ a.e. in Ω.

Proof. Let us denote 𝑢∨𝑣=sup(𝑢,𝑣), and 𝑢∧𝑣=inf(𝑢,𝑣).
From Theorem 7, It is noticed that there exist the solutions 𝑢₁ and 𝑢₂ for 𝑓₁ and 𝑓₂, respectively. Since 𝑢_𝑖 is a minimizer with data 𝑓_𝑖, for 𝑖=1,2, _Ω𝛼𝐷𝑢₁∧𝑢₂+𝜆_Ω𝑢₁∧𝑢₂+𝑓₁1log𝑢₁∧𝑢₂≥_Ω𝛼(𝑥)𝐷𝑢₁+𝜆_Ω𝑢₁+𝑓₁1log𝑢₁𝑑𝑥,_Ω𝛼𝐷𝑢₁∨𝑢₂+𝜆_Ω𝑢₁∨𝑢₂+𝑓₁1log𝑢₁∨𝑢₂≥_Ω𝛼(𝑥)𝐷𝑢₂+𝜆_Ω𝑢₂+𝑓₂1log𝑢₂𝑑𝑥. (51) Adding these two inequalities and by a minor modification of the facts in [18, 19], which is as follows: _Ω𝛼𝐷𝑢₁∧𝑢₂+_Ω𝛼𝐷𝑢₁∨𝑢₂≤_Ω𝛼𝐷𝑢₁+_Ω𝛼𝐷𝑢₂, (52) then we can deduce that _Ω𝑢₁∧𝑢₂−𝑢₁𝑑𝑥+_Ω𝑓₁1log𝑢₁∧𝑢₂1−log𝑢₁+𝑑𝑥_Ω𝑢₁∨𝑢₂−𝑢₂+𝑑𝑥_Ω𝑓₂1log𝑢₁∨𝑢₂1−log𝑢₂𝑑𝑥≥0. (53) Writing Ω={𝑢₁>𝑢₂}∪{𝑢₁≤𝑢₂}, we easily deduce that _{{𝑢1>𝑢2}}log𝑢₁−log𝑢₂𝑓₁−𝑓₂𝑑𝑥≥0. (54) Since 𝑓₁<𝑓₂, the set {𝑢₁>𝑢₂} is a zero Lebesgue measure, that is, 𝑢₁≤𝑢₂ a.e. in Ω.

4. The Associated Evolution Equations (29)–(31)

In this section, by an approach from the theory used in both [16, 20], we define the weak solution of the evolution problems (29)–(31), derive estimates for the solution of an approximating problem, prove existence and uniqueness of the weak solution, and discuss the behavior of the weak solution as 𝑡→∞.

4.1. Definition of a Pseudosolution to the Problems (29)–(31)

Denote 𝑄_𝑇=Ω×[0,𝑇), 0<𝑇≤∞. Suppose that 𝑣∈𝐿²(0,𝑇;𝐻¹(Ω))∩𝐿^∞(𝑄_𝑇) with 𝑣>0, 𝜕𝑣/𝜕𝑛=0 and 𝑢 is a classical solution of (29)–(31) with 𝑢>0. Multiplying (29) by (𝑣−𝑢), integrating over Ω, we have _Ω𝜕_𝑡𝑢(𝑣−𝑢)𝑑𝑥+_Ω𝛼(𝑥)∇𝑢||∇()∇𝑢𝑣−𝑢𝑑𝑥+𝜆_Ω𝑢−𝑓𝑢(𝑣−𝑢)𝑑𝑥=0. (55) Since _Ω∇𝑢||∇()=∇𝑢𝑣−𝑢𝑑𝑥_Ω𝛼(𝑥)∇𝑢|∇𝑢|∇𝑣𝑑𝑥−_Ω𝛼(𝑥|≤)|∇𝑢𝑑𝑥_Ω𝛼(𝑥)|∇𝑣|𝑑𝑥−_Ω𝛼(𝑥)|∇𝑢|𝑑𝑥, (56) and utilizing the Lagrange mean value theorem, _Ω1𝑣+𝑓log𝑣𝑑𝑥−_Ω1𝑢+𝑓log𝑢𝑑𝑥=_Ω(𝑣−𝑢)𝑑𝑥−_Ω𝑓𝜉()≥𝑣−𝑢𝑑𝑥_Ω(𝑣−𝑢)𝑑𝑥−_Ω𝑓𝑢(𝑣−𝑢)𝑑𝑥, (57) for either 𝑣≥𝜉≥𝑢>0 or 𝑢≥𝜉≥𝑣>0, we have _Ω𝜕_𝑡(𝑣−𝑢)𝑑𝑥+𝐸(𝑣)≥𝐸(𝑢). (58) Integrating over [0,𝑡] for any 𝑡∈[0,𝑇] then yields ^𝑡_0Ω𝜕_𝑠𝑢(𝑣−𝑢)𝑑𝑥𝑑𝑠+^𝑡₀𝐸(𝑣)𝑑𝑠≥^𝑡₀𝐸(𝑢)𝑑𝑠. (59)

On the other hand, let 𝑣=𝑢+𝜖𝜙 in (59) with 𝜙∈𝐶^∞(Ω), 𝜕𝜙/𝜕𝑛=0, and 𝜙>0. Then left-hand side of the inequality (59) has a minimum at 𝜖=0. Hence, if 𝑢 satisfies (59) and 𝑢∈𝐿²(0,𝑇;BV∩𝐿²(Ω))∩𝐿^∞(𝑄_𝑇) with 𝑢_𝑡∈𝐿²(𝑄_𝑇), 𝑢>0, and 𝜕𝑢/𝜕𝑛=0, 𝑢 is also a solution of the problems (29)–(31) in the sense of distribution. Based on this fact which is similar to the ideas in [13, 16, 21], we give the following definition.

Definition 9. A function 𝑢∈𝐿²(0,𝑇;BV(Ω)∩𝐿²(Ω))∩𝐿^∞(𝑄_𝑇) is called a pseudosolution of (29)–(31), if 𝜕𝑢/𝜕𝑡∈𝐿²(𝑄_𝑇), 𝜕𝑢/𝜕𝑛=0, 𝑢>0, and 𝑢 satisfies (59) for all 𝑡∈[0,𝑇], 𝑣∈𝐿²(0,𝑇;BV(Ω)∩𝐿²(Ω))∩𝐿^∞(𝑄_𝑇) with 𝑣>0 and 𝜕𝑣/𝜕𝑛=0. Moreover, lim_𝑡→0+Ω𝑢(𝑥,𝑡)−𝑓(𝑥)²𝑑𝑥=0. (60)

4.2. The Approximating Problem for Problems (29)–(31)

In this subsection, we consider the approximating problem 𝜕𝑢_𝑝,𝛿𝛼𝜕𝑡−div∇𝑢_𝑝,𝛿^𝑝−2∇𝑢_𝑝,𝛿𝑢+𝜆^𝑝,𝛿−𝑓_𝛿𝑢_𝑝,𝛿()=0,in𝑥,𝑡∈Ω×ℝ⁺,𝜕𝑢_𝑝,𝛿𝜕𝑛=0,on(𝑥,𝑡)∈𝜕Ω×ℝ⁺,𝑢_𝑝,𝛿(𝑥,0)=𝑓_𝛿,on𝑥∈Ω, (61)(62)(63) where 1<𝑝≤2, 𝑓_𝛿∈𝐻¹(Ω)∩𝐿^∞(Ω) such that as 𝛿→0, 𝑓_𝛿⟶𝑓in𝐿²(Ω),_Ω𝛼𝐷𝑓_𝛿⟶_Ω𝛼𝐷𝑓,0<inf_𝑥∈Ω𝑓≤𝑓_𝛿≤sup_𝑥∈Ω𝑓. (64) From Lemma 5, we can have the existence of 𝑓_𝛿.

Let us denote 𝐸_𝑝(𝑢)=1𝑝_Ω𝛼(𝑥)|∇𝑢|^𝑝+𝜆_Ω1𝑢+𝑓log𝑢𝑑𝑥, (65) which is associate with the problems (61)–(63). Combining the proof of Theorem  3.4 in [16] with the proof of Theorem 7, we also have the following lemma.

Lemma 10. Let 𝑓_𝛿∈𝐿^∞(Ω) with inf_𝑥∈Ω𝑓_𝛿>0. Then, there is a unique solution to the problem min_{𝑢∈𝑊+1,𝑝(Ω)∩𝐿∞(Ω)}𝐸_𝑝(𝑢), (66) which is satisfying 0<inf_𝑥∈Ω𝑓_𝛿≤𝑢≤sup_𝑥∈Ω𝑓_𝛿, (67) where 𝑊₊^1,𝑝(Ω)={𝑢>0,𝑢∈𝑊^1,𝑝(Ω)}.

Based on this fact, we have the following existence and uniqueness result for the problems (61)–(63).

Theorem 11. Let 𝑓_𝛿∈𝐿^∞(Ω) with inf_𝑥∈Ω𝑓_𝛿>0 and sup_𝑥∈Ω𝑓≤1. Then, the problems (61)–(63) admit a unique pseudosolution 𝑢_𝑝,𝛿∈𝐿^∞(0,∞;𝑊^1,𝑝(Ω)∩𝐿²(Ω)) with 𝜕_𝑡𝑢_𝑝,𝛿∈𝐿²(𝑄_∞), that is, ^𝑡_0Ω𝜕_𝑠𝑢𝑣−𝑢_𝑝,𝛿𝑑𝑥𝑑𝑠+^𝑡₀𝐸_𝑝(𝑣)𝑑𝑠≥^𝑡₀𝐸_𝑝𝑢_𝑝,𝛿𝑑𝑠, (68) for any 𝑡∈[0,𝑇] and 𝑣∈𝐿²(0,𝑇;𝐻¹(Ω)) with 𝑣>0 and 𝜕𝑣/𝜕𝑛=0. Moreover, 0<inf_𝑥∈Ω𝑓_𝛿≤𝑢≤sup_𝑥∈Ω𝑓_𝛿, (69) and for any 𝑇>0, ^𝑇_0Ω𝜕_𝑡𝑢_𝑝,𝛿²𝑑𝑥𝑑𝑡+sup_{𝑡∈[0,𝑇]}1𝑝_Ω𝛼∇𝑢_𝑝,𝛿^𝑝𝑑𝑥+𝜆_Ω𝑢_𝑝,𝛿+𝑓_𝛿1log𝑢_𝑝,𝛿≤1𝑑𝑥𝑝_Ω𝛼∇𝑓_𝛿^𝑝𝑑𝑥+𝜆_Ω𝑓_𝛿+𝑓_𝛿1log𝑓_𝛿𝑑𝑥. (70)

Proof. Let us fix 𝑘=inf_𝑥∈Ω𝑓_𝛿>0 and introduce the following functions: [𝑢]_𝑘=max{𝑢,𝑘}=∶𝑢,if𝑢≥𝑘,𝑘,if𝑢<𝑘. (71) We consider the auxiliary problem as follows: 𝜕𝑢_𝑝,𝛿𝛼𝜕𝑡−div∇𝑢_𝑝,𝛿^𝑝−2∇𝑢_𝑝,𝛿𝑢+𝜆^𝑝,𝛿−𝑓_𝛿[𝑢]_𝑘()=0,in𝑥,𝑡∈Ω×ℝ⁺,𝜕𝑢_𝑝,𝛿𝜕𝑛=0,on(𝑥,𝑡)∈𝜕Ω×ℝ⁺,𝑢_𝑝,𝛿(𝑥,0)=𝑓_𝛿,on𝑥∈Ω. (72)(73)(74) Since 𝑝-Laplacian operator is a maximal monotone operator, using Galerkin method and Lebesgue convergence theorem, from standard results for parabolic equations [22], we get that the problems (72)–(74) admit a unique weak solution 𝑢_𝑝,𝛿 such that 𝑢_𝑝,𝛿∈𝐿^∞0,∞;𝑊^1,𝑝(Ω)∩𝐿²(Ω,𝜕_𝑡)𝑢_𝑝,𝛿∈𝐿²𝑄_∞. (75)
Next, let us verify that the truncated function [⋅]_𝑘 in the problems (72)–(74) can be omitted. Multiply (72) by (𝑢_𝑝,𝛿−𝑘)₋, where 𝑢_𝑝,𝛿−𝑘₋=𝑢_𝑝,𝛿−𝑘,if𝑢_𝑝,𝛿−𝑘≤0,0,otherwise, (76) and integrate over Ω to get 12𝑑𝑑𝑡_Ω𝑢_𝑝,𝛿−𝑘²₋𝑑𝑥+_Ω𝛼∇𝑢_𝑝,𝛿−𝑘₋^𝑝+𝑑𝑥_Ω𝑢_𝑝,𝛿−𝑓_𝛿𝑢_𝑝,𝛿−𝑘₋[𝑢]_𝑘𝑑𝑥=0. (77) Then, 12𝑑𝑑𝑡_Ω𝑢_𝑝,𝛿−𝑘²₋𝑑𝑥≤0. (78) Therefore, (1/2)(𝑑/𝑑𝑡)∫_Ω(𝑢_𝑝,𝛿−𝑘)²₋𝑑𝑥 is decreasing in 𝑡 and since 12_Ω𝑢_𝑝,𝛿−𝑘²₋1𝑑𝑥≥0,2_Ω𝑢_𝑝,𝛿−𝑘²₋𝑑𝑥|_𝑡=0=0, (79) we have that 12_Ω𝑢_𝑝,𝛿−𝑘²₋𝑑𝑥=0, (80) for all 𝑡∈[0,𝑇], and so 𝑢(𝑡)≥𝑘=inf_𝑥∈Ω𝑓_𝛿>0ℒ−a.e.onΩ∀𝑡>0. (81) Then, 𝑢_𝑝,𝛿 is also the solution to the problems (61)–(63). Let us fix 𝐾=sup_𝑥∈Ω𝑓_𝛿. Multiplying (61) by (𝑢−𝐾)₊, where 𝑢_𝑝,𝛿−𝐾₊=𝑢_𝑝,𝛿−𝐾,if𝑢_𝑝,𝛿−𝐾≥0,0,otherwise, (82) a similar argument yields that 𝑢(𝑡)≤𝐾=sup_𝑥∈Ω𝑓_𝛿 for all 𝑡, and then (69) follows.
Moreover, multiplying (61) by 𝜕_𝑡𝑢_𝑝,𝛿 and integrating it over Ω_𝑡 yield ^𝑡_0Ω𝜕_𝑠𝑢_𝑝,𝛿²1𝑑𝑥𝑑𝑠+𝑝^𝑡₀𝜕_𝑠Ω𝛼∇𝑢_𝑝,𝛿^𝑝𝑑𝑥𝑑𝑠+𝜆^𝑡₀𝜕_𝑠𝑢_𝑝,𝛿1+𝑓log𝑢_𝑝,𝛿𝑑𝑥𝑑𝑠=0, (83) for 𝑡∈[0,𝑇]. Then (70) follows.
Finally, we will verify that the above weak solution will be the pseudosolution for the problems (61)–(63). Suppose that 𝑣∈𝐿²(0,𝑇;𝐻¹(Ω)) with 𝑣>0 and 𝜕𝑣/𝜕𝑛=0. Multiplying (61) by 𝑣−𝑢_𝑝,𝛿, then integrating by parts over Ω_𝑇, we obtain ^𝑇_0Ω𝜕_𝑡𝑢_𝑝,𝛿𝑣−𝑢_𝑝,𝛿+𝑑𝑥𝑑𝑡_𝑇0Ω𝛼∇𝑢_𝑝,𝛿^𝑝−2∇𝑣−∇𝑢_𝑝,𝛿+𝑑𝑥𝑑𝑡_𝑇0Ω𝑢_𝑝,𝛿−𝑓_𝛿𝑣−𝑢_𝑝,𝛿𝑢_𝑝,𝛿𝑑𝑥𝑑𝑡=0. (84) Utilizing Young’s inequality with 𝜖=|∇𝑢_𝑝,𝛿|^2−𝑝 yields for the second term in the left-hand side of the above inequality (84), we obtain ∇𝑢_𝑝,𝛿^𝑝−2∇𝑢_𝑝,𝛿∇𝑣−∇𝑢_𝑝,𝛿^𝑝≤|∇𝑣|^𝑝𝑝+∇𝑢_𝑝,𝛿^𝑝𝑝. (85) Utilizing the Lagrange mean value theorem yields _Ω𝑣+𝑓_𝛿1log𝑣𝑑𝑥−_Ω𝑢_𝑝,𝛿+𝑓_𝛿1log𝑢_𝑝,𝛿=𝑑𝑥_Ω𝑣−𝑢_𝑝,𝛿𝑑𝑥−_Ω𝑓_𝛿𝜉𝑣−𝑢_𝑝,𝛿≥𝑑𝑥_Ω𝑣−𝑢_𝑝,𝛿𝑑𝑥−_Ω𝑓_𝛿𝑢_𝑝,𝛿𝑣−𝑢_𝑝,𝛿𝑑𝑥, (86) for either 𝑣≥𝜉≥𝑢_𝑝,𝛿>0 or 𝑢_𝑝,𝛿≥𝜉≥𝑣>0. Returning to (84), we therefore conclude ^𝑇_0Ω𝜕_𝑡𝑢_𝑝,𝛿𝑣−𝑢_𝑝,𝛿1𝑑𝑥𝑑𝑡+𝑝^𝑇_0Ω𝛼|∇𝑣|^𝑝𝑑𝑥+𝜆^𝑇_0Ω𝑣+𝑓_𝛿1log𝑣≥1𝑑𝑥𝑑𝑡𝑝^𝑇_0Ω𝛼∇𝑢_𝑝,𝛿^𝑝𝑑𝑥+𝜆^𝑇_0Ω𝑢_𝑝,𝛿+𝑓_𝛿1log𝑢_𝑝,𝛿𝑑𝑥𝑑𝑡, (87) that is, ^𝑡_0Ω𝜕_𝑠𝑢𝑣−𝑢_𝑝,𝛿𝑑𝑥𝑑𝑠+^𝑡₀𝐸_𝑝(𝑣)≥𝑑𝑠_𝑡0𝐸_𝑝𝑢_𝑝,𝛿𝑑𝑠, (88) for any 𝑡∈[0,𝑇]. This fact implies that 𝑢_𝑝,𝛿 is a pseudosolution for the problems (61)–(63). The uniqueness follows by the uniqueness of the weak solution for the problems (61)–(63).

4.3. Existence and Uniqueness of the Problems (29)–(31)

In this subsection, we will prove the main theorem for the existence and uniqueness for the solution to the problems (29)–(31).

Theorem 12. Suppose 𝑓∈BV(Ω)∩𝐿^∞(Ω) with inf_𝑥∈Ω𝑓>0 and sup_𝑥∈Ω𝑓≤1. Then, there exists a unique pseudosolution 𝑢∈𝐿^∞(0,∞,BV(Ω)∩𝐿^∞(Ω)) to the problems (29)–(31). Moverover, 0<inf_𝑥∈Ω𝑓≤𝑢≤sup_𝑥∈Ω𝑓,_∞0Ω𝜕_𝑡𝑢²𝑑𝑥𝑑𝑡+sup_{[)𝑡∈0,∞Ω}𝛼|𝐷𝑢|+𝜆_Ω1𝑢+𝑓log𝑢≤𝑑𝑥_Ω𝛼𝐷𝑓+𝜆_Ω1𝑓+𝑓log𝑓𝑑𝑥. (89)(90)

Proof
Step 1. First we fix 𝛿>0 and pass to the limit 𝑝→1.
Let 𝑢_𝑝,𝛿 be the pseudosolution solution of (61)–(63). From (69)-(70), we know that, for fixed 𝛿>0, 𝑢_𝑝,𝛿 is uniformly bounded in 𝐿^∞(0,∞,BV(Ω)∩𝐿^∞(Ω)) and 𝜕_𝑡𝑢_𝑝,𝛿 is uniformly bounded in 𝐿²(𝑄_∞), that is, 𝑢_{𝑝,𝛿𝐿∞(𝑄∞)}+𝑢_{𝑝,𝛿𝐿∞(0,∞,BV(Ω))}𝜕≤𝐶,_𝑡𝑢_{𝑝,𝛿𝐿2(𝑄∞)}≤𝐶, (91) where 𝐶>0 is constant. Now we claim that there exists a sequence of functions 𝑢_{𝑝𝑗,𝛿} and a function 𝑢_𝛿∈𝐿^∞(0,∞,BV(Ω)∩𝐿^∞(Ω)) such that, as 𝑝_𝑗→1, 𝜕_𝑡𝑢_{𝑝𝑗,𝛿}⟶𝜕_𝑡𝑢_𝛿weaklyin𝐿²𝑄_∞,𝑢_{𝑝𝑗,𝛿}⟶𝑢_𝛿weakly^∗in𝐿^∞𝑄_∞with0<inf_𝑥∈Ω𝑓_𝛿≤𝑢_𝛿≤sup_𝑥∈Ω𝑓_𝛿,𝑢_{𝑝𝑗,𝛿}⟶𝑢_𝛿stronglyin𝐿¹(Ω)[),𝑢foreach𝑡∈0,∞_{𝑝𝑗,𝛿}⟶𝑢_𝛿in𝐿²(Ω)uniformlyin𝑡,lim_𝑡→0+Ω𝑢_{𝑝𝑗,𝛿}(𝑥,𝑡)−𝑓_𝛿(𝑥)²𝑑𝑥=0. (92)(93)(94)(95)(96)
In fact, from (91), there is a sequence 𝑢_{𝑝𝑗,𝛿} and a function 𝑢_𝛿∈𝐿^∞(𝑄_∞) with 𝜕_𝑡𝑢_𝛿∈𝐿²(𝑄_∞) such that (92) and (93) hold.
Note that, for any 𝜓∈𝐿²(Ω), as 𝑗→∞, _Ω𝑢_{𝑝𝑗,𝛿}(𝑥,𝑡)−𝑓_𝛿(𝑥)𝜓(𝑥)=𝑑𝑥_𝑄𝑡𝜕_𝑠𝑢_{𝑝𝑗,𝛿}(𝑥,𝑠)1_[0,𝑡](𝑠)𝜓(𝑥)⟶𝑑𝑥𝑑𝑠_𝑄𝑡𝜕_𝑠𝑢_𝛿(𝑥,𝑠)1_[0,𝑡](𝑠)𝜓(𝑥)=𝑑𝑥𝑑𝑠_Ω𝑢_𝛿(𝑥,𝑡)−𝑓_𝛿(𝑥)𝜓(𝑥)𝑑𝑥, (97) which shows that, for each 𝑡, 𝑢_{𝑝𝑗,𝛿}⟶𝑢_𝛿weaklyin𝐿²(Ω). (98) By (69) and (70), for each 𝑡∈[0,∞), 𝑢_𝑝𝑗(⋅,𝑡) is a bounded sequence in 𝑊^1,1(Ω). Hence, combining this with (98), we get that, for each 𝑡 as 𝑝_𝑗→1, 𝑢_{𝑝𝑗,𝛿}⟶𝑢_𝛿in𝐿¹(Ω), (99) which implies (94).
Since 𝑢_{𝑝𝑗,𝛿}(⋅,𝑡)−𝑢_{𝑝𝑗,𝛿}(⋅,𝑡^′)²_𝐿2(Ω)≤𝑡−𝑡^′_𝑄𝑇𝜕_𝑡𝑢_{𝑝𝑗,𝛿}²𝑑𝑥𝑑𝑡, (100) (96) follows.
To see (95), noting that from (100), 𝑡→𝑢_{𝑝𝑗,𝛿}(⋅,𝑡)∈𝐿²(Ω) is equicontinuous, and from (93) and (99), we have, for each 𝑡∈[0,∞), 𝑢_{𝑝𝑗,𝛿}⟶𝑢_𝛿stronglyin𝐿²(Ω). (101) Then, a standard argument yields (95).
From (70) and (92), we also obtain that 𝑢_𝛿∈𝐿^∞(0,∞,BV(Ω)∩𝐿^∞(Ω)) with 𝜕_𝑡𝑢_𝛿∈𝐿²(𝑄_∞). The claim is proved.
Next we show that, for all 𝑣∈𝐿²(0,𝑇;𝐻¹(Ω)) with 𝑣>0 and 𝜕𝑣/𝜕𝑛=0, and for each 𝑡∈[0,∞), ^𝑡_0Ω𝜕_𝑠𝑢_𝛿𝑣−𝑢_𝛿𝑑𝑥𝑑𝑠+^𝑡_0Ω𝛼|∇𝑣|𝑑𝑥𝑑𝑠+𝜆^𝑡_0Ω𝑣+𝑓_𝛿1log𝑣≥𝑑𝑥𝑑𝑠_𝑡0Ω𝛼𝐷𝑢_𝛿𝑑𝑠+𝜆^𝑡_0Ω𝑢_𝛿+𝑓_𝛿1log𝑢_𝛿𝑑𝑥𝑑𝑠. (102) By Theorem 11, we obtain that ^𝑇_0Ω𝜕_𝑡𝑢_{𝑝𝑗,𝛿}𝑣−𝑢_{𝑝𝑗,𝛿}1𝑑𝑥𝑑𝑡+𝑝_𝑗^𝑇_0Ω𝛼|∇𝑣|^𝑝𝑗𝑑𝑥𝑑𝑡+𝜆^𝑇_0Ω𝑣+𝑓_𝛿1log𝑣≥1𝑑𝑥𝑑𝑡𝑝_𝑗^𝑇_0Ω𝛼∇𝑢_{𝑝𝑗,𝛿}^𝑝𝑗𝑑𝑥𝑑𝑡+𝜆^𝑇_0Ω𝑢_{𝑝𝑗,𝛿}+𝑓_𝛿1log𝑢_{𝑝𝑗,𝛿}𝑑𝑥𝑑𝑡, (103) Notice that, using Lemma 4, we have _Ω𝛼𝐷𝑢_𝛿≤liminf_𝑗→∞Ω𝛼𝐷𝑢_{𝑝𝑗,𝛿}≤liminf_𝑗→∞Ω𝛼∇𝑢_{𝑝𝑗,𝛿}^𝑝𝑗𝑑𝑥^1/𝑝𝑗_Ω𝛼(𝑥)𝑑𝑥^{1−1/𝑝𝑗}=liminf_𝑗→∞1𝑝_𝑗Ω𝛼∇𝑢_{𝑝𝑗,𝛿}^𝑝𝑗𝑑𝑥^1/𝑝𝑗. (104) Letting 𝑗→∞, 𝑝_𝑗→1 in (103) and using (92), (94), and (104), we obtain (102) for all 𝑣∈𝐿²(0,𝑇;𝐻¹(Ω)) with 𝑣>0 and 𝜕𝑣/𝜕𝑛=0.
Step 2. Now it only remains to pass to the limit as 𝛿→0 in (102) to complete the existence of the solution to (29)–(31).
Replacing 𝑢_𝑝 by 𝑢_𝑝𝑗 in (70), letting 𝑗→∞ (𝑝_𝑗→1), and using (92)–(95), (104), and (64), we obtain ^∞_0Ω𝜕_𝑡𝑢_𝛿²𝑑𝑥𝑑𝑡+sup_{[)𝑡∈0,∞Ω}𝛼𝐷𝑢_𝛿+𝜆_Ω𝑢_𝛿+𝑓_𝛿1log𝑢_𝛿≤𝑑𝑥_Ω𝛼∇𝑓_𝛿𝑑𝑥+𝜆_Ω𝑓_𝛿+𝑓_𝛿1log𝑓_𝛿𝑑𝑥. (105) Combining this with (69), we have 𝑢_𝛿isuniformlyboundedin𝐿^∞0,∞,BV(Ω)∩𝐿^∞(Ω,𝜕_𝑡)𝑢_𝛿isuniformlyboundedin𝐿²𝑄_∞. (106)
Then, by a processing similar to the one for getting (92)–(96), we can obtain a sequence of functions 𝑢_𝛿𝑗 and a function 𝑢∈𝐿^∞(0,∞;BV(Ω)∩𝐿^∞(Ω) such that, as 𝑗→∞, 𝛿_𝑗→0, 𝜕_𝑡𝑢_𝛿𝑗⟶𝜕_𝑡𝑢weaklyin𝐿²𝑄_∞,𝑢_𝛿𝑗⟶𝑢weakly^∗in𝐿^∞𝑄_∞with0<inf_𝑥∈Ω𝑓≤𝑢≤sup_𝑥∈Ω𝑢𝑓,_𝛿𝑗⟶𝑢stronglyin𝐿¹(Ω)[),𝑢foreach𝑡∈0,∞_𝛿𝑗⟶𝑢in𝐿²(Ω)uniformlyin𝑡,lim_𝑡→0+Ω𝑢_𝛿𝑗(𝑥,𝑡)−𝑓²𝑑𝑥=0. (107)(108)(109)(110)(111)
Replacing 𝑢_𝛿 by 𝑢_𝛿𝑗 in (102), letting 𝑗→∞, 𝛿_𝑗→0, and using Lemma 4, we can have from (107)–(110) ^𝑡_0Ω𝜕_𝑠𝑢(𝑣−𝑢)𝑑𝑥𝑑𝑠+^𝑡_0Ω𝛼|∇𝑣|𝑑𝑥𝑑𝑠+𝜆^𝑡_0Ω1𝑣+𝑓log𝑣≥𝑑𝑥𝑑𝑠_𝑡0Ω𝛼|𝐷𝑢|𝑑𝑠+𝜆^𝑡_0Ω1𝑢+𝑓log𝑢𝑑𝑥𝑑𝑠, (112) for all 𝑣∈𝐿²(0,𝑇;𝐻¹(Ω)) with 𝑣>0 and 𝜕𝑣/𝜕𝑛=0. Then, the existence of the pseudosolution to (29)–(31) is completed.
Moreover, replacing 𝑢_𝛿 by 𝑢_𝛿𝑗, letting 𝑗→∞ in (105), and using (107)–(110) and (64), we have ^∞_0Ω𝜕_𝑡𝑢²𝑑𝑥𝑑𝑡+sup_{𝑡∈[0,∞)Ω}𝛼|𝐷𝑢|+𝜆_Ω1𝑢+𝑓log𝑢≤𝑑𝑥_Ω𝛼𝐷𝑓𝑑𝑥+𝜆_Ω1𝑓+𝑓log𝑓𝑑𝑥. (113)
Finally, the uniqueness of pseudosolutions to the problems (29)–(31) follows as in [21, 23]. Let 𝑢₁, 𝑢₂ be two pseudosolutions to (29)–(31) with 𝑢₁(𝑥,0)=𝑢₂(𝑥,0)=𝑓. Then, we have ^𝑡_0Ω𝜕_𝑠𝑢₁𝑢₂−𝑢₁𝑑𝑥𝑑𝑠+^𝑡_0Ω𝛼∇𝑢₂𝑑𝑥𝑑𝑠+𝜆^𝑡_0Ω𝑢₂1+𝑓log𝑢₂≥𝑑𝑥𝑑𝑠_𝑡0Ω𝛼𝐷𝑢₁𝑑𝑠+𝜆^𝑡_0Ω𝑢₁1+𝑓log𝑢₁𝑑𝑥𝑑𝑠,_𝑡0Ω𝜕_𝑠𝑢₂𝑢₁−𝑢₂𝑑𝑥𝑑𝑠+^𝑡_0Ω𝛼∇𝑢₁𝑑𝑥𝑑𝑠+𝜆^𝑡_0Ω𝑢₁1+𝑓log𝑢₁≥𝑑𝑥𝑑𝑠_𝑡0Ω𝛼𝐷𝑢₂𝑑𝑠+𝜆^𝑡_0Ω𝑢₂1+𝑓log𝑢₂𝑑𝑥𝑑𝑠. (114) Adding the above two inequalities, we get ^𝑡_0Ω𝜕_𝑠𝑢₁−𝑢₂²𝑑𝑥𝑑𝑠≤0, (115) for all 𝑡>0. This implies that 𝑢₁=𝑢₂ a.e. in (𝑥,𝑡)∈𝑄_∞.

4.4. Long-Time Behavior

At last, we will show the asymptotic limit of the solution 𝑢(⋅,𝑡) as 𝑡→∞.

Theorem 13. As 𝑡→∞, the pseudosolutions to the problems (29)–(31), 𝑢(𝑥,𝑡), converges strongly in 𝐿²(Ω) to the minimizer 𝑢 of the functional 𝐸(𝑢), that is, the solution of the problem (28).

Proof. Take a function 𝑣∈BV(Ω)∩𝐿^∞(Ω) with 𝑣>0 in (59), then _Ω𝑢()(𝑥𝑣(𝑥)1𝑥,𝑡−𝑓)𝑑𝑥−2_Ω𝑢²(𝑥,𝑡)−𝑓²(𝑥)𝑑𝑥+𝑡_Ω𝛼|∇𝑣|𝑑𝑥𝑑𝑠+𝜆𝑡_Ω1𝑣+𝑓log𝑣≥𝑑𝑥𝑑𝑠_𝑡0Ω𝛼|𝐷𝑢|𝑑𝑠+𝜆^𝑡_0Ω1𝑢+𝑓log𝑢𝑑𝑥𝑑𝑠. (116) As in [16], let 𝑤()=1𝑥,𝑡𝑡^𝑡₀𝑢(𝑥,𝑠)𝑑𝑠. (117) Since 𝑢∈𝐿^∞(0,∞;BV(Ω)∩𝐿^∞(Ω)) with 𝑢>0, for each 𝑡>0, we have that 𝑤(⋅,𝑡)∈BV(Ω)∩𝐿^∞(Ω) with {𝑤(⋅,𝑡)} uniformly bounded in BV(Ω) and 𝐿^∞(Ω). Then, there exists a subsequence {𝑤(⋅,𝑡_𝑖)} of {𝑤(⋅,𝑡)} and a function 𝑢∈BV(Ω)∩𝐿^∞(Ω), such that as 𝑡_𝑖→∞, 𝑤⋅,𝑡_𝑖⟶𝑢stronglyin𝐿¹(Ω),𝑤⋅,𝑡_𝑖⟶𝑢weakly^∗inBV(Ω). (118) Since {𝑤(⋅,𝑡)} is uniformly bounded in 𝐿^∞(Ω), we have 𝑤⋅,𝑡_𝑖⟶𝑢stronglyin𝐿²(Ω). (119)
By dividing 𝑡 in (116) and then taking the limit along 𝑠_𝑖→∞, we get that, for any 𝑣∈BV(Ω)∩𝐿^∞ with 𝑣>0, _Ω𝛼|𝐷𝑣|𝑑𝑥+𝜆_Ω1𝑣+𝑓log𝑣≥𝑑𝑥_Ω𝛼|𝐷𝑢|𝑑𝑠+𝜆_Ω1𝑢+𝑓log𝑢𝑑𝑥, (120) which implies that 𝑢 is the minimizer of the problem (28).

5. Numerical Methods and Experimental Results

We present in this section some numerical examples illustrating the capability of our model. We also compare it with the known model (AA). In the next two subsections, two numerical discrete schemes, the 𝛼-total variation (𝛼-TV) scheme and the 𝑝-Laplace approximate (𝑝-LA) scheme, will be proposed.

5.1. 𝛼-Total Variation Scheme

Numerically we get a solution to the problem (28) by computing the associated equation (29) to a steady state. To discretize (29), the finite difference scheme in [6] is used. Denote the space step by ℎ=1 and the time step by 𝜏. Thus, we have 𝐷^±_𝑥𝑢_𝑖,𝑗=±𝑢_𝑖±1,𝑗−𝑢_𝑖,𝑗,𝐷^±_𝑦𝑢_𝑖,𝑗=±𝑢_𝑖,𝑗±1−𝑢_𝑖,𝑗,𝐷_𝑥𝑢_𝑖,𝑗=𝐷⁺_𝑥𝑢_𝑖,𝑗²+𝑚𝐷⁺_𝑦𝑢_𝑖,𝑗,𝐷⁻_𝑦𝑢_𝑖,𝑗²𝐷+𝜖,_𝑥𝑢_𝑖,𝑗=𝐷⁺_𝑦𝑢_𝑖,𝑗²+𝑚𝐷⁺_𝑥𝑢_𝑖,𝑗,𝐷⁻_𝑥𝑢_𝑖,𝑗²+𝜖, (121) where 𝑚[𝑎,𝑏]=(sign𝑎+sign𝑏)min(|𝑎|,|𝑏|)/2 and 𝜖>0 is the regularized parameter chosen near 0.

The numerical algorithms for the problems (29)–(31) are given as follows: 𝛼_𝑖,𝑗=11−1+𝑘𝐺_𝜎∗𝑓²_𝑖,𝑗1+𝑘𝑀²𝑘𝑀²,𝑍^𝑛_𝑖,𝑗=𝐷⁻_𝑥𝛼_𝑖,𝑗𝐷⁺_𝑥𝑢^𝑛_𝑖,𝑗𝐷⁺_𝑥𝑢^𝑛_𝑖,𝑗+𝐷⁻_𝑦𝛼_𝑖,𝑗𝐷⁺_𝑦𝑢^𝑛_𝑖,𝑗𝐷⁺_𝑦𝑢^𝑛_𝑖,𝑗,𝜆^𝑛=∑_𝑖,𝑗𝑍^𝑛_𝑖,𝑗1−𝑓/(𝑢+𝜖)𝐿,𝑢𝐼𝐽_{𝑛+1𝑖,𝑗}=𝑢^𝑛_𝑖,𝑗+𝜏𝑍^𝑛_𝑖,𝑗−𝜏𝜆^𝑛𝑓1−,𝑢𝑢+𝜖_0𝑖,𝑗=𝑓_𝑖,𝑗𝑢=𝑓𝑖ℎ,𝑗ℎ,0≤𝑖≤𝐼,0≤𝑗≤𝐽,_𝑛𝑖,0=𝑢^𝑛_𝑖,1,𝑢^𝑛_0,𝑗=𝑢^𝑛_1,𝑗𝑢^𝑛_𝐼,𝑖=𝑢^𝑛_{𝐼−1,𝑖},𝑢^𝑛_𝑖,𝐽=𝑢^𝑛_{𝑖,𝐽−1}. (122)(123)(124)(125)(126)(127)

Here the MATLAB function “conv2" is used to represent the two-dimensional discrete convolution transform of the matrix 𝑢_𝑖,𝑗, that is, 𝐺∗𝑢. Through the above lines, we can obtain 𝑢¹_𝑖,𝑗 by 𝑢⁰_𝑖,𝑗. The program will stop when it achieves our desire.

5.2. 𝑝-Laplace Approximate Scheme

From the proof of Theorem 12, we can know that the term ∫_Ω𝛼|𝐷𝑢| can be approximated by the term (1/𝑝)∫_Ω𝛼|∇𝑢|^𝑝𝑑𝑥. Based on this, we can use the numerical algorithms for the problems (61)–(63) to get a solution to the problem (28), as 𝑝→1. As in [24], the numerical discrete scheme for the problems (61)–(63) is given as follows: 𝛼_𝑖,𝑗=11−1+𝑘𝐺_𝜎∗𝑓²_𝑖,𝑗1+𝑘𝑀²𝑘𝑀²,𝑍^𝑛_1𝑖,𝑗=𝛼_𝑖,𝑗𝐷⁺_𝑥𝑢^𝑛_𝑖,𝑗²+𝐷⁺_𝑦𝑢^𝑛_𝑖,𝑗^{2(𝑝𝑛−2)/2}𝐷⁺_𝑥𝑢^𝑛_𝑖,𝑗,𝑍^𝑛_2𝑖,𝑗=𝛼_𝑖,𝑗𝐷⁺_𝑥𝑢^𝑛_𝑖,𝑗²+𝐷⁺_𝑦𝑢^𝑛_𝑖,𝑗^{2(𝑝𝑛−2)/2}𝐷⁺_𝑦𝑢^𝑛_𝑖,𝑗,𝜆^𝑛=∑_𝑖,𝑗𝐷⁻_𝑥𝑍^𝑛_1𝑖,𝑗+𝐷⁻_𝑦𝑍^𝑛_2𝑖,𝑗1−𝑓/(𝑢+𝜖)𝐿,𝑢𝐼𝐽_{𝑛+1𝑖,𝑗}=𝑢^𝑛_𝑖,𝑗+𝜏𝐷⁻_𝑥𝑍^𝑛_1𝑖,𝑗+𝐷⁻_𝑦𝑍^𝑛_2𝑖,𝑗−𝜏𝜆^𝑛𝑓1−,𝑝𝑢+𝜖_𝑛+1=𝑝^𝑛+𝑛𝑝_𝐿−𝑝_𝑈𝑁,𝑝⁰=𝑝_𝑈,𝑢⁰_𝑖,𝑗=𝑓_𝑖,𝑗,𝑢=𝑓𝑖ℎ,𝑗ℎ0≤𝑖≤𝐼,0≤𝑗≤𝐽,_𝑛𝑖,0=𝑢^𝑛_𝑖,1,𝑢^𝑛_0,𝑗=𝑢^𝑛_1,𝑗,𝑢^𝑛_𝐼,𝑖=𝑢^𝑛_{𝐼−1,𝑖},𝑢^𝑛_𝑖,𝐽=𝑢^𝑛_{𝑖,𝐽−1}, (128) where 𝑝_𝐿=1, 1<𝑝_𝑈≤2 is upper limit of the exponent 𝑝, 𝑛=1,2,…,𝑁, and 𝑁 is the iteration time.

5.3. Comparison with Other Methods

In this section, we used the similar way for numerical experiments as [25]. For comparison purposes, some very recent multiplicative noise removal algorithms from the literature are considered, such as the SO Algorithm [8] (see (16)-(17)) and the AA Algorithm [5] (see (10)). As recommended in [8], the stopping rule about SO Algorithm is to reach 𝐾 such that 𝐾=max{𝑘∈ℕ∶Var[𝑤_𝑘−𝑤_𝑂]≥Var[𝜂]=Ψ₁(𝐿)}, where 𝑤_𝑂 is the underlying log-image, and 𝜂 is the relevant noise, and the variance is defined by (16).

The denoising algorithms were tested on three images: a synthetic image (300×300 pixels), an aerial image (256×256 pixels), and a cameraman image (256×256 pixels). In all numerical experiments by our algorithms, the images do not need to be normalized in the range [1,256]. This is different from the other algorithms. For each image, a noisy observation is generated by multiplying the original image by the speckle noise according to the model in (8) with the choice 𝐿∈{1,4,10}.

For a noise-free image 𝑢_𝑂 and its denoised version by any algorithm 𝑢, the denoising performance is measured in terms of peak signal to noise ratio (PSNR) [25], PSNR=10log₁₀𝑀𝑁max𝑢_𝑂−min𝑢_𝑂²𝑢−𝑢_𝑂²_𝐿2dB, (129) and mean absolute-deviation error (MAE), MAE=𝑢−𝑢_𝑂𝐿1,𝑀𝑁 (130) where |max𝑢_𝑂−min𝑢_𝑂| gives the gray-scale range of the original image, 𝑀,𝑁 is the size of the image.

For fair comparison, the parameters of SO and AA were tweaked manually to reach their best performance level. Their values are summarized in Table 1. In the 𝛼-total variation (𝛼-TV) scheme, there are four parameters: the influencing factor 𝑘, the scale of convolution 𝜎, the time step 𝜏, and the turning factor 𝜆. However, the parameter 𝜆 is found dynamically by (124). To ensure stability as well as optimal results, we choose 𝜏=0.05. In 𝑝-Laplace approximate (𝑝-LA) scheme, besides the same four parameters with 𝛼-TV scheme, there is a new parameter 𝑝_𝑈 with 1<𝑝_𝑈≤2. We need not the the exact value of 𝑝_𝑈 but an approximate estimate (e.g., 𝑝_𝑈=1.2). Notice that the parameters of our method are very stable with respect to the image.

Table 1 

Parameters used in the comparison study.

The results are depicted in Figures 2, 3, and 4 for the synthetic image, Figures 5, 6, and 7 for the aerial image, and Figures 8, 9, and 10 for the cameraman image. Our methods do a good job at restoring faint geometrical structures of the images even for low values of 𝐿, see for instance the results on the aerial image for 𝐿=1 and 𝐿=4. Our algorithm performs among the best and even outperforms its competitors most of the time both visually and quantitatively as revealed by the PSNR and MAE values. For SO method, the number of iterations necessary to satisfy the stopping rule rapidly increases when 𝐿 decreases [25]. For AA method, the appropriate parameter 𝜆 is necessary.

(a) Noisy: 𝐿=1, PSNR = −3.46

(b) Original

(c) 𝛼-TV: PSNR = 19.67, MAE = 2.45

(d) 𝑝-LA: PSNR = 19.76, MAE = 2.51

(e) SO: PSNR = 4.14, MAE = 23.94

(f) AA: PSNR = 16.33, MAE = 4.50

Figure 2 

Synthetic image (300×300). (a) Noisy image corrupted by speckle noise for 𝐿=1 in (8). (b) Original image. (c) Our algorithm by 𝛼-TV, 𝜎=2, 𝑘=0.05, and 𝜏=0.05. (d) Our algorithm by 𝑝-LA, 𝑝𝑈=1.3. (e) SO algorithm, 𝜆=0.1, 𝛼=0.25. (f) AA algorithm, 𝜆=50.

(a) Noisy: 𝐿=4, PSNR = 2.52

(b) Original

(c) 𝛼-TV: PSNR = 23.86, MAE = 1.27

(d) 𝑝-LA: PSNR = 23.14, MAE = 1.43

(e) SO: PSNR = 15.18, MAE = 6.15

(f) AA: PSNR = 20.14, MAE = 2.94

Figure 3 

Synthetic image (300×300) (a) Noisy image corrupted by speckle noise for 𝐿=4 in (8). (b) Original image. (c) Our algorithm by 𝛼-TV, 𝜎=2, 𝑘=0.05, 𝜏=0.05. (d) Our algorithm by 𝑝-LA, 𝑝𝑈=1.3. (e) SO algorithm, 𝜆=0.7, 𝛼=0.25. (f) AA algorithm, 𝜆=150.

(a) Noisy: 𝐿=10, PSNR = 6.52

(b) Original

(c) 𝛼-TV: PSNR = 26.37, MAE = 0.92

(d) 𝑝-LA: PSNR = 25.68, MAE = 0.96

(e) SO: PSNR = 21.00, MAE = 2.84

(f) AA: PSNR = 23.94, MAE = 1.49

Figure 4 

Synthetic image (300×300) (a) Noisy image corrupted by speckle noise for 𝐿=10 in (8). (b) Original image. (c) Our algorithm by 𝛼-TV, 𝜎=2, 𝑘=0.05, 𝜏=0.05. (d) Our algorithm by 𝑝-LA, 𝑝𝑈=1.2. (e) SO algorithm, 𝜆=0.7, 𝛼=0.25. (f) AA algorithm, 𝜆=240.

(a) Noisy: 𝐿=1, PSNR = 12.39

(b) Original

(c) 𝛼-TV: PSNR = 23.25, MAE = 12.00

(d) 𝑝-LA: PSNR = 23.47, MAE = 11.64

(e) SO: PSNR = 17.82, MAE = 25.44

(f) AA: PSNR = 22.34, MAE = 13.35

Figure 5 

Aerial image (512×512). (a) Noisy image corrupted by speckle noise for 𝐿=1 in (8). (b) Original image. (c) Our algorithm by 𝛼-TV, 𝜎=2, 𝑘=0.03, and 𝜏=0.05. (d) Our algorithm by 𝑝-LA, 𝑝𝑈=1.4. (e) SO algorithm, 𝜆=0.1, 𝛼=0.25. (f) AA algorithm, 𝜆=25.

(a) Noisy: 𝐿=4, PSNR = 18.39

(b) Original

(c) 𝛼-TV: PSNR = 25.82, MAE = 8.95

(d) 𝑝-LA: PSNR = 26.16, MAE = 8.81

(e) SO: PSNR = 24.00, MAE = 11.05

(f) AA: PSNR = 24.20, MAE = 10.30

Figure 6 

Aerial image (512×512). (a) Noisy image corrupted by speckle noise for 𝐿=4 in (8). (b) Original image. (c) Our algorithm by 𝛼-TV, 𝜎=2, 𝑘=0.015, and 𝜏=0.05. (d) Our algorithm by 𝑝-LA, 𝑝𝑈=1.4. (e) SO algorithm, 𝜆=0.3, 𝛼=0.25. (f) AA algorithm, 𝜆=120.

(a) Noisy: 𝐿=10, PSNR = 22.45

(b) Original

(c) 𝛼-TV: PSNR = 27.92, MAE = 7.12

(d) 𝑝-LA: PSNR = 28.18, MAE = 6.98

(e) SO: PSNR = 27.27, MAE = 7.42

(f) AA: PSNR = 26.04, MAE = 8.13

Figure 7 

Aerial image (512×512). (a) Noisy image corrupted by speckle noise for 𝐿=10 in (8). (b) Original image. (c) Our algorithm by 𝛼-TV, 𝜎=2, 𝑘=0.015, and 𝜏=0.05. (d) Our algorithm by 𝑝-LA, 𝑝𝑈=1.4. (e) SO algorithm, 𝜆=1.2, 𝛼=0.25. (f) AA algorithm, 𝜆=130.

(a) Noisy: 𝐿=1, PSNR = 5.26

(b) Original

(c) 𝛼-TV: PSNR = 20.81, MAE = 14.65

(d) 𝑝-LA: PSNR = 21.03, MAE = 13.94

(e) SO: PSNR = 14.15, MAE = 38.92

(f) AA: PSNR = 18.81, MAE = 22.14

Figure 8 

Cameraman image (512×512). (a) Noisy image corrupted by speckle noise for 𝐿=1 in (8). (b) Original image. (c) Our algorithm by 𝛼-TV, 𝜎=2, 𝑘=0.01, and 𝜏=0.05. (d) Our algorithm by 𝑝-LA, 𝑝𝑈=1.3. (e) SO algorithm, 𝜆=0.04, 𝛼=0.25. (f) AA algorithm, 𝜆=125.

(a) Noisy: 𝐿=4, PSNR = 11.23

(b) Original

(c) 𝛼-TV: PSNR = 24.22, MAE = 8.88

(d) 𝑝-LA: PSNR = 24.10, MAE = 9.18

(e) SO: PSNR = 21.61, MAE = 14.64

(f) AA: PSNR = 21.91, MAE = 14.25

Figure 9 

Cameraman image (512×512). (a) Noisy image corrupted by speckle noise for 𝐿=4 in (8). (b) Original image. (c) Our algorithm by 𝛼-TV, 𝜎=2, 𝑘=0.005, and 𝜏=0.05. (d) Our algorithm by 𝑝-LA, 𝑝𝑈=1.2. (e) SO algorithm, 𝜆=0.2, 𝛼=0.25. (f) AA algorithm, 𝜆=240.

(a) Noisy: 𝐿=10, PSNR = 15.27

(b) Original

(c) 𝛼-TV: PSNR = 26.29, MAE = 7.14

(d) 𝑝-LA: PSNR = 26.03, MAE = 7.43

(e) SO: PSNR = 24.96, MAE = 8.36

(f) AA: PSNR = 24.25, MAE = 9.97

Figure 10 

Cameraman image (512×512). (a) Noisy image corrupted by speckle noise for 𝐿=10 in (8). (b) Original image. (c) Our algorithm by 𝛼-TV, 𝜎=8/6, 𝑘=0.005, and 𝜏=0.05. (d) Our algorithm by 𝑝-LA, 𝑝𝑈=1.2. (e) SO algorithm, 𝜆=1, 𝛼=0.25. (f) AA algorithm, 𝜆=390.

In the numerical experiments, we can see that for the nontexture image, our methods and AA method work well (see Figure 4); for the texture image, our methods and SO method work well (see Figure 7). The denoising performance results are tabulated in Table 2, where the best PSNR and MAE values are shown in boldface. The PSNR improvement brought by our approach can be quite high particularly for 𝐿=1 (see e.g., Figures 2–4), and the visual resolution is quite respectable. This is an important achievement since in practice 𝐿 has a small value, usually 1, rarely above 4. This improvement becomes less salient as 𝐿 increases which is intuitively expected. But even for 𝐿=10, the PSNR of our algorithm can be higher than AA and SO methods (see Table 2).

Table 2 

PSNR and MAE.

Acknowledgments

The authors would like to express their sincere thanks to the referees for their valuable suggestions for the revision of the paper which contributed greatly to this work. The authors would also like to thank Jalal Fadili for providing them the MATLAB code of his algorithm. They would like to express their deep thanks to the referees for their suggestions while revising the paper, yet again. This work was partially supported by the Fundamental Research Funds for the Central Universities (Grant nos. HIT.NSRIF.2011003 and HIT.NSRIF.2012065), the National Science Foundation of China (Grant no. 11271100), the Aerospace Supported Fund, China, under Contract no. 2011-HT-HGD-06, China Postdoctoral Science Foundation funded project, Grant no. 2012M510933, and also the 985 project of Harbin Institute of Technology.