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

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 u denote the image intensity function, t the time. The anisotropic diffusion as formulated by Perona and Malik [1] can be presented as


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][2][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 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 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 where  is the so-called constant parameter.Then the Perona-Malik's model is equivalent to minimizing where Ω ⊂  2 is the image domain.The energy functional (4) is minimized when |∇| 2 is minimum, which leads to piecewise constant approximation of .Therefore, formation of staircase on the ramp edges is unavoidable.
The solutions of the minimization problem of (5) after using Euler-Lagrange equation followed by gradient descent procedure is given by 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 fourthorder 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][21][22][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.

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 where Ω ⊂  2 ,  1 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 On the other hand, we replace (5) by where Ω ⊂  2 and  2 is a fixed positive constant with the contribution as  1 .Then, the minimizing functional (10) yields the associated Euler-Lagrange equation In this section, we have treated  1 () and  2 () 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.

Convex Combination of the Two Minimization Problems
In this section, we denote the solutions ( 9), ( 11) by  and V, respectively.It follows from the Euler-Lagrange variation principle that the minimizer of  and the minimizer of V can be interpreted as the steady-state solution of the nonlinear diffusion process with initial data (, , 0) =  0 (, ), and with the same initial data V(, , 0) =  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 V 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 − )V 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.

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, . . ., , where  ×  is the size of the image, and then a 4-nearestneighbors discretization of the Laplacian operator can be used: where Secondly, ( 13) still can be discredited on a square lattice as described above.We calculate the Laplacian of the image intensity function as with symmetric boundary conditions , ), which can be discredited as with symmetric boundary conditions Thus, the numerical approximation to the differential equation ( 13) is given as Thirdly, we deal with the convex combination Noticing that  and V 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,  1 = 0, and  = 0,  2 = 0 in (22), respectively.

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     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 = Variance of image Variance of noise , PSNR = 10 log 10 ( 255 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.

Conclusions
This paper proposes a new model for noise removal.The new model is based on a convex combination of the secondorder 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.
(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,  1 = 0.02,  2 = 0.002,  = 10,  1 = 25, and  2 = 100.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 Figure2.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.

Figure 1 :
Figure 1: Recovered results via our proposed model and compared with the Perona-Malik's model and the You-Kaveh model.

Figure 2 :
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.

Table 1 :
The comparison of the SNRs and PSNRs for experiments.