Image Restoration Combining the Second-Order and Fourth-Order PDEs

Abstract

A noise removal technique using partial differential equations (PDEs) is proposed. It combines a second-order filter with a fourth-order filter. The combined method takes the advantage of both filters since it can preserve edges and at the same time avoid the blocky effects in smooth regions. The experimental results illustrate the effectiveness of the model in image restoration.

1. Introduction

In the last two decades, the second-order partial differential equations have been well studied by many scholars as one of the useful tools for the image restoration problem. For instance, the anisotropic diffusion model [1–3], the total variation models [4], and the curve evolution equations [5], have been demonstrated to be effective for removing noise and edge preservation. However, the images resulting from these second-order models are often piecewise constant, and therefore, the processed image suffers from the so-called blocky effects, which make it be visually uncomfortable.

To be precise, we first give a brief description about the blocky effects associated with anisotropic diffusion. Let 𝑢 denote the image intensity function, 𝑡 the time. The anisotropic diffusion as formulated by Perona and Malik [1] can be presented as 𝜕𝑢𝑔,𝜕𝑡=∇⋅(|∇𝑢|)∇𝑢 (1) where 𝑔 is the diffusion coefficient and ∇⋅ and ∇ denote the divergence and the gradient, respectively. You et al. [6] carried out a detailed analysis to show that the solution of (1) is equal to the minimization of energy functional 𝐸(𝑢)=_Ω𝑓(|∇𝑢|)𝑑𝑥𝑑𝑦. (2) From energy functional, it is obvious that level images are the global minima of the energy functional. The analysis in [6] indicates that when there is no backward diffusion, a level image is the only minimum of the energy functional, so Perona-Malik’s model will evolve toward the formation of a level image function. Since Perona-Malik’s model is designed such that smooth areas are diffused faster than less smooth ones, blocky effects will appear in the early stage of the diffusion and will develop as time evolves.

In particular, one of the classical diffusivity functions defined in [1] is given by 𝑔(𝑥)=11+(𝑥/𝑘)², (3) where 𝑘 is the so-called constant parameter. Then the Perona-Malik’s model is equivalent to minimizing 𝐸(𝑢)=_Ω𝑘²2ln𝑘²+|∇𝑢|²𝑑𝑥𝑑𝑦, (4) where Ω⊂𝑅² is the image domain. The energy functional (4) is minimized when |∇𝑢|² is minimum, which leads to piecewise constant approximation of 𝑢. Therefore, formation of staircase on the ramp edges is unavoidable.

To reduce the blocky effect, high-order PDEs (typically, fourth-order PDEs) have been introduced into image restoration [7–18]. In 2000, You and Kaveh [7] proposed a family of fourth-order partial differential equations (the You-Kaveh model). They considered the second-order functional 𝐸(𝑢)=_Ω𝑓∇²𝑢𝑑𝑥𝑑𝑦, (5) where 𝑓^′(𝑠)=𝑠𝑔(𝑠) and |∇²𝑢| is simply an absolute value of Laplacian of 𝑢 approximated by |𝑢_𝑥𝑥+𝑢_𝑦𝑦|. For the diffusivity function in (3), the energy functional (5) is in the form of 𝐸(𝑢)=_Ω𝑘²2ln𝑘²+∇²𝑢𝑑𝑥𝑑𝑦, (6) meaning that (6) minimized when |∇²𝑢| is minimum. Therefore, the ramp region of 𝑢 (i.e., the regions where |∇²𝑢|=0) are fit in the solution of the associate fourth-order PDE. The solutions of the minimization problem of (5) after using Euler-Lagrange equation followed by gradient descent procedure is given by 𝜕𝑢𝜕𝑡=−∇²𝑔∇²𝑢∇²𝑢. (7) The You-Kaveh model replaces the gradient operator in the Perona-Malik’s model with a Laplacian operator. Due to the fact that the Laplacian of an image at a pixel is zero only if the image is planar in its neighborhood, the You-Kaveh fourth-order PDE attempts to remove noise and preserve edge by approximating an observed image with a piece planar image. It is well known that piecewise smooth images look more natural.

The further theoretical analysis in [10, 19] shows that fourth-order equations have advantages over second-order equations in some aspects. First, fourth-order linear diffusion dampens oscillations at high frequencies (i.e., noise) much faster than second order diffusion. Second, there is the possibility of having schemes that include effects of curvature (i.e., the second derivatives of the image) in the dynamics, thus creating a richer set of functional behaviors [19]. Therefore, the blocky effect will be reduced and image will look more natural. However, the fourth-order equation of the type You-Kaveh model tends to leave images with speckle artifacts.

Therefore, both the Perona-Malik’s model and the You- Kaveh model have their strengths and weaknesses depending on the characteristics of the image of interest. Motivated by [1, 6, 7, 12, 17], the aim of this paper is to generate a new solution by taking the best from each of the two methods by a convex combination. For other recent studies on the noise removal by using the second- or fourth-order diffusion PDEs, we refer to [20–23].

The outline of this paper is as follows. Section 2 gives a detailed description of two minimization problems. A fourth-order PDE together with a second-order is the basic ingredients in our proposed model. The way these two PDEs interfere with each other is discussed in Section 3. Section 4 elaborates on the numerical method for our proposed model. And experimental results are provided in Section 5, followed by some conclusions in Section 6.

2. Description of Two Minimization Problems

We use functionals 𝐸_𝑖, 𝑖=1,2 to measure the quality of the restoration process. Smaller values of 𝐸_𝑖 correspond to a result that reflects features (flat, smooth, and jumps) in a better way than larger values do. Instead of (2), we consider 𝐸₁(𝑢)=_Ω𝑓+𝜆(|∇𝑢|)¹2𝑢−𝑢₀²𝑑𝑥𝑑𝑦, (8) where Ω⊂𝑅², 𝜆₁ is a fixed positive constant that balances the regularity of the solution and the fidelity. The minimizing functional (8) yields the associated Euler-Lagrange equation ∇⋅𝑔(|∇𝑢|)∇𝑢−𝜆₁𝑢−𝑢₀=0. (9) On the other hand, we replace (5) by 𝐸₂(𝑢)=_Ω𝑓∇²𝑢+𝜆₂2𝑢−𝑢₀²𝑑𝑥𝑑𝑦, (10) where Ω⊂𝑅² and 𝜆₂ is a fixed positive constant with the contribution as 𝜆₁. Then, the minimizing functional (10) yields the associated Euler-Lagrange equation −∇²𝑔∇²𝑢∇²𝑢−𝜆₂𝑢−𝑢₀=0. (11)

In this section, we have treated 𝐸₁(𝑢) and 𝐸₂(𝑢) and their associated Euler-Lagrange equations separately. However, we want to establish a positive interaction between these equations, and that is the topic for the next section.

3. Convex Combination of the Two Minimization Problems

In this section, we denote the solutions (9), (11) by 𝑢 and 𝑣, respectively. It follows from the Euler-Lagrange variation principle that the minimizer of 𝑢 and the minimizer of 𝑣 can be interpreted as the steady-state solution of the nonlinear diffusion process 𝑢_𝑡=∇⋅𝑔(|∇𝑢|)∇𝑢−𝜆₁𝑢−𝑢₀, (12) with initial data 𝑢(𝑥,𝑦,0)=𝑢₀(𝑥,𝑦), and 𝑣_𝑡=−∇²𝑔∇²𝑣∇²𝑣−𝜆₂𝑣−𝑣₀, (13) with the same initial data 𝑣(𝑥,𝑦,0)=𝑢₀(𝑥,𝑦), respectively. As mentioned in the last section, each of the above PDEs substantially suppress noise, but (12) is designed such that smooth areas are diffused faster than less smooth ones and thus the blocky effects will appear, while (13) attempts to preserve edges by approximating an observed image with a piecewise planar image at the cost of leaving images with speckle artifacts. Then, we do not expect their solutions 𝑢 and 𝑣 to be equal all over the image domain Ω.

Considering that the methods in (12) and (13) have their strengths and weakness, we try to generate a new model by a convex combination 𝑤=𝛼𝑢+(1−𝛼)𝑣 with 𝛼∈[0,1] to fully take advantage of the strengths of (12) and (13). We prefer that the weighting constant 𝛼 can be found adaptively. Through several different approaches to calculate the weighting constant, we have found that the assumption 𝛼≤1/2 could give good results. Indeed, we will take 𝛼=0.315. The details of the algorithm we have used are given in the next section. We remark that the theoretical analysis of the best constant 𝛼 for the convex combination is out of the scope of this paper.

4. Discredited Numerical Scheme

In this section, we use a simple numerical scheme that discrete (12), (13) and then combine them. For this purpose, we divide it into three steps.

Firstly, (12) can be discredited on a square lattice with the horizontal and vertical directions having the same step of space. Suppose that ℎ denotes the spatial mesh size and Δ𝑡 the temporal step length. We quantize the space and time coordinates as follows: 𝑡=𝑘_𝑖∗Δ𝑡,𝑘_𝑖=0,1,2,…,(𝑖=1,2),𝑥=𝑖∗ℎ,𝑖=0,1,2,…,𝑀,𝑦=𝑗∗ℎ,𝑗=0,1,2,…,𝑁, (14) where 𝑀×𝑁 is the size of the image, and then a 4-nearest-neighbors discretization of the Laplacian operator can be used: 𝑢^𝑘1+1_𝑖,𝑗=𝑢^𝑘1_𝑖,𝑗+Δ𝑡∗𝑐_𝑁⋅∇_𝑁𝑢+𝑐_𝑆⋅∇_𝑆𝑢+𝑐_𝐸⋅∇_𝐸𝑢+𝑐_𝑊⋅∇_𝑊𝑢^𝑘1_𝑖,𝑗−𝜆₁𝑢^𝑘1_𝑖,𝑗−𝑢^𝑘1₀, (15) where ∇_𝑁𝑢^𝑘1_𝑖,𝑗=𝑢^𝑘1_{𝑖−1,𝑗}−𝑢^𝑘1_𝑖,𝑗,∇_𝑆𝑢^𝑘1_𝑖,𝑗=𝑢^𝑘1_𝑖+1,𝑗−𝑢^𝑘1_𝑖,𝑗,∇_𝐸𝑢^𝑘1_𝑖,𝑗=𝑢^𝑘1_𝑖,𝑗+1−𝑢^𝑘1_𝑖,𝑗,∇_𝑊𝑢^𝑘1_𝑖,𝑗=𝑢^𝑘1_{𝑖,𝑗−1}−𝑢^𝑘1_𝑖,𝑗,𝑐^𝑘1_𝑁=𝑔∇_𝑁𝑢^𝑘1_𝑖,𝑗,𝑐^𝑘1_𝑆=𝑔∇_𝑆𝑢^𝑘1_𝑖,𝑗,𝑐^𝑘1_𝐸=𝑔∇_𝐸𝑢^𝑘1_𝑖,𝑗,𝑐^𝑘1_𝑊=𝑔∇_𝑊𝑢^𝑘1_𝑖,𝑗. (16)

Secondly, (13) still can be discredited on a square lattice as described above. We calculate the Laplacian of the image intensity function as ∇²𝑣^𝑘2_𝑖,𝑗=𝑣^𝑘2_𝑖+1,𝑗+𝑣^𝑘2_{𝑖−1,𝑗}+𝑣^𝑘2_𝑖,𝑗+1+𝑣^𝑘2_{𝑖,𝑗−1}−4𝑣^𝑘2_𝑖,𝑗 (17) with symmetric boundary conditions 𝑣^𝑘2_−1,𝑗=𝑣^𝑘2_0,𝑗,𝑣^𝑘2_𝑀+1,𝑗=𝑣^𝑘2_𝑀,𝑗𝑣,𝑗=0,1,…,𝑁,_{𝑘2𝑖,−1}=𝑣^𝑘2_𝑖,0,𝑣^𝑘2_𝑖,𝑁+1=𝑣^𝑘2_𝑖,𝑁,𝑖=0,1,…,𝑀. (18) Let 𝑐^𝑘2_𝑖,𝑗=𝑔(|∇²𝑣^𝑘2_𝑖,𝑗|∇²𝑣^𝑘2_𝑖,𝑗), which can be discredited as ∇²𝑐^𝑘2_𝑖,𝑗=𝑐^𝑘2_𝑖+1,𝑗+𝑐^𝑘2_{𝑖−1,𝑗}+𝑐^𝑘2_𝑖,𝑗+1+𝑐^𝑘2_{𝑖,𝑗−1}−4𝑐^𝑘2_𝑖,𝑗ℎ², (19) with symmetric boundary conditions 𝑐^𝑘2_−1,𝑗=𝑐^𝑘2_0,𝑗,𝑐^𝑘2_𝑀+1,𝑗=𝑐^𝑘2_𝑀,𝑗𝑐,𝑗=0,1,…,𝑁,_{𝑘2𝑖,−1}=𝑐^𝑘2_𝑖,0,𝑐^𝑘2_𝑖,𝑁+1=𝑐^𝑘2_𝑖,𝑁,𝑖=0,1,…,𝑀. (20) Thus, the numerical approximation to the differential equation (13) is given as 𝑣̃^𝑘2+1_𝑖,𝑗=𝑣^𝑘2_𝑖,𝑗−Δ𝑡∇²𝑐^𝑘2_𝑖,𝑗−𝜆₂𝑣^𝑘2_𝑖,𝑗−𝑣^𝑘2₀. (21)

Thirdly, we deal with the convex combination 𝑢=𝛼𝑢^𝑘1+1+(1−𝛼)𝑣̃^𝑘2+1. (22) Noticing that 𝑢 and 𝑣 can be found independently each other, we can combine them when they are convergent. Numerical tests indicate that a combination at convergence is most effective and accurate. Each of the numerical schemes (12) and (13) is stable if they are solved separately, as long as Δ𝑡 fulfills the Courant-Friedrichs-Lewy (CFL) condition. Note that the corresponding algorithm for the Perona-Malik’s model and the You-Kaveh’s model can be given by setting 𝛼=1, 𝜆₁=0, and 𝛼=0, 𝜆₂=0 in (22), respectively.

5. Experimental Results

In this section, we present some of the results obtained by the proposed model and compare them with the corresponding ones for the Perona-Malik’s model given by solving PDE (1) and the You-Kaveh’s model given by solving PDE (7). From the experimental results, the new model presented in this paper can performance better than Perona-Malik’s model and You-Kaveh’s model. In particular, the new model can reduce the blocky effects appeared in Perona-Malik’s model and avoid leaving the speckle artifacts appeared in the You-Kaveh’s model.

Our example is a 256×256 sized gray-scale image Lena, which is displayed in Figure 1(a). Figure 1(b) is its degraded version corrupted by white random Gaussian noise with standard deviation 15. Then, Figure 1(c) is the recovered results by employing the Perona-Malik’s model, Figure 1(d) is the recovered results by employing the You-Kaveh’s model, and Figure 1(e) is the recovered results by employing the proposed model. Figure 1(c) is obtained with Δ𝑡=0.2, 𝑘=10 for iteration 25, Figure 1(d) is obtained with Δ𝑡=0.2, 𝑘=5 for iteration 300, while our result is carried out by setting 𝛼=0.315, Δ𝑡=0.185, 𝜆₁=0.02, 𝜆₂=0.002, 𝑘=10, 𝑘₁=25, and 𝑘₂=100.

(a) Original image

(b) Noisy image

(c) Perona-Malik’s model

(d) You-Kaveh’s model

(e) Our proposed model

Figure 1 

Recovered results via our proposed model and compared with the Perona-Malik’s model and the You-Kaveh model.

In order to better understand the behavior of the proposed model in local regions, especially in regions with smooth signals and regions with discontinuities, we present the following zoomed-in local results.

A small part of the Lena image is shown in Figure 2. It is clear that the Perona-Malik’s model appears obvious blocky effect and the You-Kaveh model leaves the speckle artifacts. Our proposed model can avoid the staircase and the speckle artifacts while removing the noise.

(a) Original image

(b) Noisy image

(c) Perona-Malik’s model

(d) You-Kaveh’s model

(e) Our proposed model

Figure 2 

Partially enlarged results are displayed to compare the denoising performance of the Perona-Malik’s model and the You-Kaveh’s model with our proposed model.

The restoration quality can be quantitatively measured by the signal-to-noise ratio (SNR) and the peak signal-to-noise ratio (PSNR), which are defined as SNR=Varianceofimage,VarianceofnoisePSNR=10log10255²∑_𝑖𝑗𝑔_𝑖𝑗−ℎ_𝑖𝑗², (23) and, respectively, where 𝑔 is the original image, ℎ denotes the compared image, and the unit of SNR(PSNR) is decibel (dB).

In Table 1, we give the comparison of the SNR and PSNRs for Figure 1, which shows that our model has the better SNR and PSNR than those of the Perona-Malik’s model and the You-Kaveh’s model.

Table 1 

The comparison of the SNRs and PSNRs for experiments.

6. Conclusions

This paper proposes a new model for noise removal. The new model is based on a convex combination of the second-order filter with the fourth-order filter. We have tested our algorithm on images consisting of edges and smooth regions. From these experimental results, we observed that the proposed method is able to preserve edges while at the same time avoiding the blocky effects in smooth regions. In a word, the combined model reaps benefits of both the Perona-Malik’s model and the You-Kaveh’s model, surpassing each individually in image restoration.

Acknowledgments

The authors are very grateful to the referee for his/her valuable comments, which greatly improved the paper. This work was supported in part by the NNSF of China under Grant 11101068, the Sichuan Youth Science & Technology Foundation under Grant 2011JQ0003, and the SRF for ROCS, SEM.