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

Authors: Tianhua Liu; Zhaoyin Xiang
Journal: Mathematical Problems in Engineering
 (2013)

---

## 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.


---

## Body

## 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. Letu denote the image intensity function, t the time. The anisotropic diffusion as formulated by Perona and Malik [1] can be presented as
(1)∂u∂t=∇·(g(|∇u|)∇u),
where g 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)E(u)=∫Ωf(|∇u|)dxdy.
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
(3)g(x)=11+(x/k)2,
where k is the so-called constant parameter. Then the Perona-Malik’s model is equivalent to minimizing
(4)E(u)=∫Ωk22ln(k2+|∇u|2)dxdy,
where Ω⊂R2 is the image domain. The energy functional (4) is minimized when |∇u|2 is minimum, which leads to piecewise constant approximation of u. 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)E(u)=∫Ωf(|∇2u|)dxdy,
where f′(s)=sg(s) and |∇2u| is simply an absolute value of Laplacian of u approximated by |uxx+uyy|. For the diffusivity function in (3), the energy functional (5) is in the form of
(6)E(u)=∫Ωk22ln(k2+|∇2u|)dxdy,
meaning that (6) minimized when |∇2u| is minimum. Therefore, the ramp region of u (i.e., the regions where |∇2u|=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)∂u∂t=-∇2(g(|∇2u|)∇2u).
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. Section2 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 functionalsEi, i=1,2 to measure the quality of the restoration process. Smaller values of  Ei correspond to a result that reflects features (flat, smooth, and jumps) in a better way than larger values do. Instead of (2), we consider
(8)E1(u)=∫Ωf(|∇u|)+λ12(u-u0)2dxdy,
where Ω⊂R2, λ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
(9)∇·(g(|∇u|)∇u)-λ1(u-u0)=0.
On the other hand, we replace (5) by
(10)E2(u)=∫Ωf(|∇2u|)+λ22(u-u0)2dxdy,
where Ω⊂R2 and λ2 is a fixed positive constant with the contribution as λ1. Then, the minimizing functional (10) yields the associated Euler-Lagrange equation
(11)-∇2(g(|∇2u|)∇2u)-λ2(u-u0)=0.In this section, we have treatedE1(u) and E2(u) 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 u and v, respectively. It follows from the Euler-Lagrange variation principle that the minimizer of u and the minimizer of v can be interpreted as the steady-state solution of the nonlinear diffusion process
(12)ut=∇·(g(|∇u|)∇u)-λ1(u-u0),
with initial data u(x,y,0)=u0(x,y), and
(13)vt=-∇2(g(|∇2v|)∇2v)-λ2(v-v0),
with the same initial data v(x,y,0)=u0(x,y), 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 u 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 w=αu+(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.


## 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 h denotes the spatial mesh size and Δt the temporal step length. We quantize the space and time coordinates as follows:
(14)t=ki*Δt,ki=0,1,2,…,(i=1,2),x=i*h,i=0,1,2,…,M,y=j*h,j=0,1,2,…,N,
where M×N is the size of the image, and then a 4-nearest-neighbors discretization of the Laplacian operator can be used:
(15)u~i,jk1+1=ui,jk1+Δt*((ui,jk1-u0k1)(cN·∇Nu+cS·∇Su00000000000+cE·∇Eu+cW·∇Wu)i,jk10000000000-λ1(ui,jk1-u0k1)),
where
(16)∇Nui,jk1=ui-1,jk1-ui,jk1,∇Sui,jk1=ui+1,jk1-ui,jk1,∇Eui,jk1=ui,j+1k1-ui,jk1,∇Wui,jk1=ui,j-1k1-ui,jk1,cNk1=g(|∇Nui,jk1|),cSk1=g(|∇Sui,jk1|),cEk1=g(|∇Eui,jk1|),cWk1=g(|∇Wui,jk1|).Secondly, (13) still can be discredited on a square lattice as described above. We calculate the Laplacian of the image intensity function as
(17)∇2vi,jk2=vi+1,jk2+vi-1,jk2+vi,j+1k2+vi,j-1k2-4vi,jk2
with symmetric boundary conditions
(18)v-1,jk2=v0,jk2,vM+1,jk2=vM,jk2,j=0,1,…,N,vi,-1k2=vi,0k2,vi,N+1k2=vi,Nk2,i=0,1,…,M.
Let ci,jk2=g(|∇2vi,jk2|∇2vi,jk2), which can be discredited as
(19)∇2ci,jk2=ci+1,jk2+ci-1,jk2+ci,j+1k2+ci,j-1k2-4ci,jk2h2,
with symmetric boundary conditions
(20)c-1,jk2=c0,jk2,cM+1,jk2=cM,jk2,j=0,1,…,N,ci,-1k2=ci,0k2,ci,N+1k2=ci,Nk2,i=0,1,…,M.
Thus, the numerical approximation to the differential equation (13) is given as
(21)v~i,jk2+1=vi,jk2-Δt(∇2ci,jk2-λ2(vi,jk2-v0k2)).Thirdly, we deal with the convex combination(22)u=αu~k1+1+(1-α)v~k2+1.
Noticing that u 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 Δt 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.


## 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 a256×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 Δt=0.2,  k=10 for iteration 25, Figure 1(d) is obtained with Δt=0.2, k=5 for iteration 300, while our result is carried out by setting α=0.315, Δt=0.185, λ1=0.02, λ2=0.002, k=10, k1=25, and  k2=100.Recovered results via our proposed model and compared with the Perona-Malik’s model and the You-Kaveh model.
(a)
Original image
(b)
Noisy image
(c)
Perona-Malik’s model
(d)
You-Kaveh’s model
(e)
Our proposed modelIn 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.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.
(a)
Original image
(b)
Noisy image
(c)
Perona-Malik’s model
(d)
You-Kaveh’s model
(e)
Our proposed modelThe 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(23)SNR=VarianceofimageVarianceofnoise,PSNR=10log10(2552∑ij(gij-hij)2),
and, respectively, where g is the original image, h denotes the compared image, and the unit of SNR(PSNR) is decibel (dB).In Table1, 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.
Image
Peronal-Malik's model
You-Kaveh's model
Our proposed model
SNR (dB)
23.14387
21.56607
23.29253
PSNR (dB)
28.0408
26.46297
28.18947


## 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.



---

*Source: 743891-2013-06-27.xml*