An anisotropic PDE model for image inpainting

Computers and Mathematics with Applications xxx (xxxx) xxx

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

An anisotropic PDE model for image inpainting Abdul Halim a,b , B.V. Rathish Kumar a , a b

∗

Department of Mathematics and Statistics, IIT Kanpur, 208016, India Department of Mathematics, Hari Singh College, Munger, Bihar 81213, India

article

info

Article history: Received 9 July 2019 Received in revised form 19 November 2019 Accepted 1 December 2019 Available online xxxx Keywords: Inpainting 4th-order PDE Cahn-Hilliard equation Fidelity term Convexity splitting

a b s t r a c t A PDE model based on an image specific auto generated multi-well potential function and a 4th order anisotropic diffusion with good performance with respect to loss due to smoothing effects has been proposed for grayscale image inpainting. An unconditionally stable numerical scheme using the notion of convexity splitting in time and Fourier spectral method in space has been derived. The stable scheme is both consistent and convergent. Numerical computations are done for some standard test images and results are compared with the results of other existing models in the literature using image analysis specs such as PSNR, SNR, and SSIM. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Image inpainting is an important area of research in the field of image processing. It is a process of filling in the missing part of an image from the information of its surroundings. It has lot of applications in real life such as restoration of ancient frescoes, removal of scratches of photos and it is used in medical fields to reduce artifacts in MRI, CT, PET images. There are several methods for inpainting such as exemplar-based [1], texture synthesis [2], stochastic methods [3], waveletbased [4,5] and PDE methods [6–9]. We will focus on PDE-based methods. PDE-based methods are both theoretically and computationally sound as they are totally based and backed by the well established theory of PDEs and numerical analysis. Many a time they provide a revealing insight into the problem which helps one to come up with a new model and or robust and cost-effective method of solution. The first work in this direction was proposed by Bertalmio et al. [10] in 2000. They followed the technique of museum artists and have proposed a 3rd order nonlinear PDE model for inpainting. In the subsequent work [11] authors bring out that the nonlinear PDE proposed in the earlier work [10] is related to two dimensional Navier–Stokes (NS) equation. They have shown that the problem proposed in [10] is analogous to solving the inviscid Euler equation for incompressible flow where the image intensity acts like stream function in the fluid problem and the boundary condition and the isophote lines imposed on the boundary are nothing but the no slip boundary condition that is required for NS equation to introduce a diffusion term. Chan and Shen [7] introduced another class of inpainting model called variational inpainting model. They have used the idea of Rudin, Osher, and Fatemi model [12] for image denoising. The equation is of 2nd order. The problem with the 2nd order models is that they are not able to fill the large gaps and are unable to connect the level lines in the missing domain. So the authors Chan and Shen proposed [8] a third order model based on total variation (TV) which is called CDD method. Their work is motivated by the work of Perona and Malik [13]. To overcome the problem of 2nd order model ∗ Corresponding author. E-mail addresses: [email protected] (A. Halim), [email protected], [email protected] (B.V. Rathish Kumar). https://doi.org/10.1016/j.camwa.2019.12.002 0898-1221/© 2019 Elsevier Ltd. All rights reserved.

Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

2

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

the authors Chan et al. [14] have proposed a 4th-order model called Euler Elastica model. This is motivated by work of Nitzberg, Mumford, and Shiota [15] model. Esedoglu and Shen [16] utilizes the Mumford and Shah image segmentation model to inpainting problem. The corresponding Euler–Lagrange equation results in a 2nd order PDE and hence faces the connectivity problem. So to improve this model they have used the idea Euler-Elastica [14] model and proposed Mumford–Shah–Euler image model in the same paper [16]. This model results in a 4th-order nonlinear parabolic PDE containing a small parameter. Then Bertozzi, Esedoglu, and Gillette [9] proposed a new inpainting model based on Cahn–Hilliard equation and the model is analyzed in [17]. This model is faster but it smoothens out the edges of images during the inpainting of binary images and cannot handle gray scale images due to the presence of the double-well potential. Schönlieb et al. [6] have generalized the inpainting model proposed in [9] for gray scale image, known as TV − H −1 model. But the TV − H −1 still has the problem of smoothening the edges overlapping with inpainting region and it is also costly. Hence there is a need for a cost effective and a robust PDE based model for grayscale image inpainting free from edge loss. To overcome the restriction of binary images a multi-well potential based on histogram information is introduced in place of double-well potential of Cahn–Hilliard (CH) model [9]. Next, to avoid the edge smoothing effect the biharmonic term in the CH model is replaced by the 4th-order anisotropic diffusion term analogous to the one in TV − H −1 model. The presence of multi-well potential forces a pixel intensity towards the nearest minima and thereby accelerates the solution convergence process. Absence of such a term in TV − H −1 model makes it slower w.r.t. solution convergence. Since our image specific multi-well potential reduces to double-well potential for binary images by default our model is valid for binary images. It may also be noted that CH model is executed in two steps. In the first step it repairs the image and the 2nd step improves the image. We also the take advantage of the two stage processes by using two different ϵ . But the TV − H −1 model has no ϵ and it smoothens the whole image without any scope for recovery of some of the lost information due to smoothing operation. But for our case, in the 2nd stage it recovers the brightness of the image except for the domain of inpainting in other words our model does not make any changes outside the domain of inpainting. In a way our idea is a technical amalgamation of these two notions to enable us to get a faster and edge-preserving model and also to get output as close as to the original one. For validating this closeness we also calculated the SSIM [18] which calculates the similarity of two images. The proposed model is performing better than the model of Schönlieb [6]. The paper is organized as follows: In Section 2, we will present convexity splitting idea and related topics. In Section 3, we discuss the proposed model and its numerical scheme. In Section 4, we have discussed the consistency, stability and convergence of the time discretized scheme. Numerical results are presented in Section 5. Finally, in Section 6 some conclusions are provided. 2. Convexity splitting In this section, we will discuss the convexity splitting method. Since our model is a 4th order model if we solve it by explicit method then it may impose restriction on time step size of form (∆x)4 where ∆x denotes the step size of the spatial domain. Hence we need to take a very small time step size which results in a slow convergent method. If fully implicit method is used then due to nonlinear terms it will involve inner iterations and it turns out to be time consuming. This kind of methods are widely used to solve optimization problem but it was originally developed for solving gradient system. Authors in [19] have shown that this method is also applicable for evolution equations not following variational principle. Firstly, we will discuss the convexity splitting method for gradient system. Consider the following gradient system: ut = −∇ E(u) u(·, x) = u0

in

in Ω

(1)

Ω

(2) n

n

where E is a twice continuously differentiable function from R → R and u0 ∈ R . For strictly convex functional E, unconditionally stable schemes are available but for non-convex functionals E the convexity splitting methods are used to get stable methods. In convexity splitting, we split the functional E as a sum of a convex and a concave functional and in the semi-implicit scheme, the convex part is treated implicitly and the concave part explicitly in time. Thus we write E as E(u) = Ec (u) − Ee (u) where Ec and Ee are strictly convex functional from Rn → R for all u ∈ Rn . The semi-implicit scheme for the gradient system (1) is as follows:

(

Uk+1 − Uk = −∆t ∇ Ec (Uk+1 ) − ∇ Ee (Uk )

)

(3)

with initial iterate U0 = u0 and time step ∆t. Definition ([19,20]). A one-step numerical integration scheme is said to be unconditionally gradient stable if a function F : Rn → R exists which satisfies the following: Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

3

F (u) ≥ 0 ∀u ∈ Rn F (u) → ∞ as u → ∞ F (uk+1 ) ≤ F (uk ) for all u ∈ Rn and if F (uk ) = F (u0 ) ∀k ≥ 0, implies u0 is a zero of ∇ F .

(i) (ii) (iii) (iv)

for all ∆t ≥ 0 and for any given data. The above definition is valid for a gradient system. But Cahn–Hilliard model and TV − H −1 models are not a gradient flow, so the authors of [19] define the notion of unconditional stability for these inpainting models as follows: Definition ([19]). Let Ω ⊂ R2 be open and bounded and u belongs to some appropriate function space H defined on Ω × [0, T ], with T ≥ 0. Let ut = G(u, Dα u) be a partial differential equation where Dα u are derivatives in space for α = 1, . . . , 4 and Uk+1 = Uk + ∆tGk (Uk , Uk+1 , Dα Uk , Dα Uk+1 ),

(4)

be a discrete time stepping method where Gk is an approximation of G in Uk and Uk+1 . Then the scheme (4) is said to be (a) unconditionally stable, if solutions of (4) are bounded for all ∆t ≥ 0 and all k satisfying k∆t ≤ T i.e. the stability of the scheme is not dependent on ∆t, (b) consistent if lim∆t →0 τk (∆t) = 0 i.e the discretized PDE is consistent with the continuous PDE as ∆t → 0 where τk is the local truncation error of the scheme having the form:

τk (∆t) =

uk+1 − uk

∆t

− Gk (uk , uk+1 , Dα uk , Dα uk+1 )

(5)

where uk = u(k∆t) denote the exact solution at time k∆t. The global truncation error is defined as

τ (∆t) = max ∥τk (∆t)∥H k

Moreover, if τ (∆t) = O(∆t p ) for all ∆t → 0, then p is called the order of the numerical in time. 3. Proposed model and numerical scheme 3.1. Proposed model Let Ω be the image domain and D be the part of the image where information is missing called as the inpainting domain. Let u(x) be the inpainted image of the given image f (x). To obtain the solution u(x) Bertozzi et al. [9] have used the following modified version of Cahn–Hilliard equation: ut = −ϵ ∆2 u +



1

ϵ

∆W ′ (u) + λ(f − u)     fidelity part

Cahn−Hilliard

(6)

part

where

λ(x) =

{

0 x∈D

λ0 x ∈ Ω \ D

(7)

and W (u) = u2 (u − 1)2 , the double-well potential function. But due to the isotropic nature of the biharmonic term this model smoothen out the edges of the image. So to avoid this we replace the biharmonic term by the 4th order nonlinear diffusion term which is known to be edge preserving smoothing term. Hence the modified binary model takes the form:

(

) ( ∇u ) 1 + ∆W ′ (u) + λ(f − u) ut = −ϵ ∆ ∇ · |∇ u| ϵ

(8)

But for the theoretical and numerical purpose, we replace |∇ u| by its regularized version

|∇ u|δ =

√ |∇ u|2 + δ 2

(9)

where δ > 0 is a constant. Then Eq. (8) reduces to

(

) ( ∇u ) 1 ut = −ϵ ∆ ∇ · + ∆W ′ (u) + λ(f − u) |∇ u|δ ϵ

(10)

Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

4

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

3.2. Multi-well potential based grayscale inpainting model Now we will discuss briefly, the construction of multi-well potential from the histogram which was proposed in [21]. First authors have generalized the histogram by using Dirac delta symbol and further for theoretical and numerical soundness they have used a regularized version of Dirac delta to define the generalized histogram. The generalized histogram of a given image u0 has the following form:

∫ ∫ Hu0 ,α (k) := where δα (s) =

Ω

1 d

δα (|u0 (x) − k| + |y|)dydx

(

1 + tan−1

π ds

∫ ∫ Hu0 ,α (k) :=

(11)

R

Ω

R

( s )) α

. After simplifying the generalized histogram has the following form:

α 1 dydx. π (|u0 (x) − k| + |y|)2 + α 2

(12)

We will construct the multi-well potential from this generalized histogram. For simplicity let Hu0 be the histogram of the initial image u0 . Calculate the normalized histogram ˜ H(k) := |Ω1 | Hu0 (k). Then the required multi-well potential is

)2

(

W1 (u) := 1 − ˜ H(u)

. Further, for theoretical purpose we project the multi-well potential to the space of polynomial of

degree 2n. Suppose that W1 (u) is the multi-well potential then it has the form, W1 (u) =

2n ∑

ai ui

with

a2n > 0

(13)

bi ui

(14)

i=1

and let W1′ (u) =

2n−1

d dx

W1 (u) =

∑ i=1

So for gray image inpainting we use the equation

( ut = −ϵ ∆ ∇ ·

∇u

( √

)

)

|∇ u| + δ 2 2

1

+ ∆W1′ (u) + λ(f − u) ϵ

(15)

Since the multi-well potential is image specific by construction it naturally takes the form of double-well potential for binary images i.e (15) reduces to (10) for binary images. 3.3. Numerical scheme For simplicity we will discuss the numerical scheme of the proposed model (10). It will also hold for model (15) by replacing W (u) with W1 (u). For time discretization we will use convexity splitting method and for the space discretization Fourier spectral method is used. Eq. (10) can also be obtained as a gradient descent of two energy functional. The first two terms can be obtained as a gradient of the energy E1 (u) =

∫ ( Ω

) 1 ϵ|∇ u|δ + W (u) dx ϵ

(16)

in H −1 norm and the fidelity part can be obtained as gradient flow in L2 norm of the energy

∫ E2 (u) =

Ω

λ(f − u)2 dx.

(17)

For convexity splitting, we split both the energy E1 and E2 as a sum of a convex and a concave part i.e. E1 = E1c − E1e and E2 = E2c − E2e where

∫ E1c (u) = E1e (u) =

C1 ( Ω

2

∫ ( C1 Ω

2

) |∇ u|2 + |u|2 dx

(18)

(|∇ u|2 + |u|2 ) − ϵ|∇ u|δ −

1

ϵ

)

W (u) dx

(19)

and

∫ E2c =

Ω

C2 2

|u| dx, 2

E2e =

∫ ( C2 Ω

2

) |u|2 − λ(f − u)2 dx

(20)

Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

5

C1 and C2 are positive constant and are chosen in such a way that the functionals E1c , E1e and E2c , E2e become strictly convex. Here we have to choose C1 > max{ 1ϵ , 1δ } and C2 > λ0 . The resulting time stepping scheme for an initial condition U0 = u0 is given by Uk+1 − Uk

∆t

( ) ( ) = −∇H −1 E1c (Uk+1 ) − E1e (Uk ) − ∇L2 E2c (Uk+1 ) − E2e (Uk )

where ∇H −1 and ∇L2 denote gradient with respect to the inner product H −1 and L2 respectively. This leads to the scheme of the form Uk+1 − Uk

∆t

(

(

+ C1 ∆ Uk+1 − C1 ∆Uk+1 + C2 Uk+1 =C1 ∆ Uk − C1 ∆Uk − ϵ ∆ ∇ · √ 2

2

+ λ(f − Uk )

∇ Uk |∇ Uk | + δ 2

) 2

)

1

+ ∆W ′ (Uk ) ϵ (21)

We impose the following Neumann boundary conditions

∇ Uk+1 · η = ∇ ∆Uk+1 · η = 0 on ∂ Ω

(22)

where η is the outward normal vector to the boundary. For space discretization, we will apply discrete Fourier transformation (DFT) to the time-discretized scheme (21). We ˆ will compute Uk+1 in the spectral domain. Let Uˆ be the DFT of U. Taking DFT of Eq. (21) and expressing ∆ U in terms of Uˆ we get the scheme of the form

( Uˆ k+1 (i, j) = Uˆ k (i, j) +

with M(i, j) =

)

M(i, j) + λ(fˆ − Uk )(i, j) ∆t (23)

1 + C2 ∆t + C1 ∆tL2i,j − C1 ∆tLi,j

( ˆ ) ′ (U )(i, j) − ϵ L ˆ W ∇ · ( |∇∇UUk| ) (i, j). L k i , j i , j ϵ

1

k δ

ˆ Here we have used the fact that ∆ U = LUˆ derived in [22] where L = (Li,j ) is a matrix whose size is m × n which is same as the size of U with Li,j =

2 (∆x)2

( cos

( 2π i ) m

) −1 +

2 (∆y)2

(

( 2π j )

cos

n

) −1 .

4. Consistency, stability and convergence Theorem 1. Let u be the exact solution of (10) and uk = u(k∆t) be the exact solution at time k∆t for a time step ∆t > 0 and k ∈ N. Let C1 > max{ 1ϵ , 1δ } and C2 > λ0 and Uk be the kth iterate of (21). Then the following statements are true: (a) (Consistency) Assuming ∥utt ∥−1 and ∥∇ ∆ut ∥2 are bounded, the numerical scheme (21) is consistent with Eq. (10) and of order 1 in time. (b) (Unconditional stability) Under the assumption W ′′ (Uk−1 ) ≤ K

(24)

the approximate solution Uk is bounded on a time interval [0, T ] where T is finite, for all ∆t > 0. In particular, for k∆t ≤ T , we have for every ∆t > 0,

∥∇ Uk ∥2 + ∆tK1 ∥∆Uk ∥2 + ∆tK2 ∥∇ ∆Uk ∥2 ( ) ≤ eK3 T ∥∇ U0 ∥2 + ∆tK1 ∥∆U0 ∥2 + ∆tK2 ∥∇ ∆U0 ∥2 + ∆tC (Ω , D, λ0 , f )

(25)

(c) (Convergence) The discretization error ek = uk − Uk . For smooth solution uk and Uk , the error ek converges to 0 as ∆t → 0. In particular, we have for k∆t ≤ T that

∥∇ ek ∥2 + ∆tM1 ∥∆ek ∥2 + ∆tM2 ∥∇ ∆ek ∥2 ≤

T M3

eM4 T (∆t)2

(26)

for some positive constants M1 , M2 , M3 and M4 . Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

6

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

Proof. (a) Putting u = uk in (10) and Uk = uk in (21) and subtracting former from the later one we get the truncation error as

) ( ∇u ) k τk = + C1 ∆ uk+1 − C1 ∆uk+1 + C2 uk+1 − C1 ∆ uk + C1 ∆uk + ϵ ∆ ∇ · ∆t |∇ uk |δ ) } { ( } ( ∇u ) 1 1 k ′ ′ − ∆W (uk ) − λ(f − uk ) − ∆W (uk ) − λ(f − uk ) − ut (k∆t) + ϵ ∆ ∇ · ϵ |∇ uk |δ ϵ uk+1 − uk − ut (k∆t) + C1 ∆2 (uk+1 − uk ) − C1 ∆(uk+1 − uk ) + C2 (uk+1 − uk ) = ∆t {

(

uk+1 − uk

2

2

(27)

Applying standard Taylor series arguments

τk =

{(

)

}

(

)

ut (k∆t) + ∆tutt (ξ1 ) − ut (k∆t) + C1 ∆t ∆2 ut (ξ2 ) − C1 ∆ut (ξ3 ) + C2 ut (ξ4 )

(28)

Using the bounded assumption we get, lim∆t →0 τk = 0 i.e the time discretized equation is consistent with the continuous equation. To get the order of numerical scheme we need to know the global truncation error τ . Notice that ∥τk ∥−1 = O(∆t), hence the global truncation error τ given by

τ = max ∥τk ∥−1 = O(∆t) as ∆t → 0. k

Hence the order of the numerical scheme is 1 in time.

□

Proof. (b) Let us consider the time discretized model Uk+1 − Uk

∆t

(

(

+ C1 ∆ Uk+1 − C1 ∆Uk+1 + C2 Uk+1 =C1 ∆ Uk − C1 ∆Uk − ϵ ∆ ∇ · √ 2

2

∇ Uk |∇ Uk | + δ 2

) 2

)

1

+ ∆W ′ (Uk ) ϵ

+ C2 Uk + λ(f − uk )

(29)

Multiplying the above equation with −∆Uk+1 and integrate over Ω we get 1 (

) ∥∇ Uk+1 ∥2 − ⟨∇ Uk , ∇ Uk+1 ⟩ + C1 ∥∇ ∆Uk+1 ∥2 + C1 ∥∆Uk+1 ∥2 + C2 ∥∇ Uk+1 ∥2 ∆t ⟨ ( ) ⟩ 1 ∇ Uk = C1 ⟨∇ ∆Uk , ∇ ∆Uk+1 ⟩ + C1 ⟨∆Uk , ∆Uk+1 ⟩ − ϵ ∇∇ · √ , ∇ ∆Uk+1 + ⟨W ′′ (Uk )∇ Uk , ∇ ∆Uk+1 ⟩ ϵ |∇ Uk |2 + δ 2 + C2 ⟨∇ Uk , ∇ Uk+1 ⟩ + ⟨∇λ(f − Uk ), ∇ Uk+1 ⟩ Using the Young’s inequality we obtain 1 ( 2∆t C1

≤

δ1

∥∇ Uk+1 ∥2 − ∥∇ Uk ∥

2

)

+ C1 ∥∇ ∆Uk+1 ∥2 + C1 ∥∆Uk+1 ∥2 + C2 ∥∇ Uk+1 ∥2

∥∇ ∆Uk+1 ∥2 + C1 δ1 ∥∇ ∆Uk+1 ∥2 +

C1 2

∥∆Uk ∥2 +

C1 2

∥∆Uk+1 ∥2 + δ2 ∥∇ ∆Uk+1 ∥2 +

( ) 2 ϵ ∇ Uk ∥∇∇ · √ ∥ δ2 |∇ Uk |2 + δ 2

δ3 C2 1 + ∥W ′′ (Uk )∇ Uk ∥2 + ∥∇ ∆Uk+1 ∥2 + ∥∇ Uk ∥2 + C2 δ4 ∥∇ Uk+1 ∥2 + ∥∇λ(f − Uk )∥2 + δ5 ∥∇ Uk+1 ∥2 (30) 2ϵδ3 2ϵ δ4 δ5 1

Now we use the following estimates

∥∇λ(f − Uk )∥2 ≤ 2λ20 ∥∇ Uk ∥2 + C (Ω , D, λ0 , f ) and ∥W ′′ (Uk )∇ Uk ∥2 ≤ K 2 ∥∇ Uk ∥2 (by using (24))

(31)

Applying Poincaré’s and Cauchy’s inequality we get

) ( ) 2 C( ∇ Uk ∥∇∇ · √ ∥ ≤ ∥∇ Uk ∥2 + ∥∆Uk ∥2 + ∥∇ ∆Uk ∥2 δ |∇ Uk |2 + δ 2 Expressing the L2 norm of ∆u in terms of the L2 norms of ∇ u and ∇ ∆u we get,

( ) 2 ∇ Uk ∥∇∇ · √ ∥ ≤ C¯ (∥∇ Uk ∥2 + ∥∇ ∆Uk ∥2 ) 2 |∇ Uk | + δ 2

(32)

Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

7

where C¯ = C ( 1δ , Ω ). Using (31) and (32) in (30) and rearranging the terms we get

( 1 ) ( C1 δ3 ) + C2 (1 − δ4 ) − δ5 ∥∇ Uk+1 ∥2 + ∥∆Uk+1 ∥2 + C1 (1 − δ1 ) − δ2 − ∥∇ ∆Uk+1 ∥2 2∆t 2 2ϵ ( 1 (C 2λ2 ) C¯ ϵ K2 C2 C1 C¯ ϵ ) 1 ≤ + + + + 0 ∥∇ Uk ∥2 + ∥∆Uk ∥2 + + ∥∇ ∆Uk ∥2 + C (Ω , D, λ0 , f ). 2∆t δ2 2ϵδ3 δ4 δ5 2 δ1 δ2 Now choosing δi =

1 2

for i = 1, 4, 5 and δ2 =

1 4

, δ3 =

ϵ 2

we obtain,

( 1 C2 − 1 ) C1 C1 − 1 + ∥∇ Uk+1 ∥2 + ∥∆Uk+1 ∥2 + ∥∇ ∆Uk+1 ∥2 2∆t 2 2 2 ( 1 ) C1 K2 ≤ + 4ϵ C¯ + 2 + 2C2 + 4λ20 ∥∇ Uk ∥2 ∥∆Uk ∥2 + (2C1 + 4ϵ C¯ )∥∇ ∆Uk ∥2 + C (Ω , D, λ0 , f ). 2∆t ϵ 2 Since C 1 > max{ 1ϵ , 1δ } > 1 and C2 > λ0 > 1 all the coefficients in (33) are positive.

(33)

2

Multiply the above inequality with 2∆t and by taking Ca = 1 + ∆t(C2 − 1), Cb = C1 − 1, Cc = 1 + 2∆t(4ϵ C¯ Kϵ 2 + 2C2 + 4λ20 ) and Cd = 4(C1 + 2ϵ C¯ ) we get Ca ∥∇ Uk+1 ∥2 + ∆tC1 ∥∆Uk+1 ∥2 + ∆tCb ∥∇ ∆Uk+1 ∥2

≤ Cc ∥∇ Uk ∥2 + ∆tC1 ∥∆Uk ∥2 + ∆tCd ∥∇ ∆Uk ∥2 + 2∆tC (Ω , D, λ0 , f ). Dividing by Ca (> 0) we have

∥∇ Uk+1 ∥2 + ∆t

C1

∥∆Uk+1 ∥2 + ∆t

Cb

∥∇ ∆Uk+1 ∥2 Ca Cc C1 Cd 2 ≤ ∥∇ Uk ∥2 + ∆t ∥∆Uk ∥2 + ∆t ∥∇ ∆Uk ∥2 + ∆tC (Ω , D, λ0 , f ) Ca Ca Ca Ca Ca

Rewrite the right hand as C1

Cb

∥∇ Uk+1 ∥2 + ∆t ∥∆Uk+1 ∥2 + ∆t ∥∇ ∆Uk+1 ∥2 Ca Ca ) ∆ tC ∆t 2 Cc Cd ( 1 1 ≤ ∥∇ Uk ∥2 + ∥∆Uk ∥2 + ∥∇ ∆Uk ∥2 + ∆tC (Ω , D, λ0 , f ) Cc

Cd

Cc Cd

Cc

Ca

Now multiplying the first term by Cd , second term by hand side we get, C1

Cc Cd ( Ca

> 1) and third term by

Cc Cb ( Ca

> 1) within the bracket of right

Cb

∥∇ Uk+1 ∥2 + ∆t ∥∆Uk+1 ∥2 + ∆t ∥∇ ∆Uk+1 ∥2 Ca Ca ) Cc Cd ( C Cb 2 1 ≤ ∥∇ Uk ∥2 + ∆t ∥∆Uk ∥2 + ∆t ∥∇ ∆Uk ∥2 + ∆tC (Ω , D, λ0 , f ) Cc

Ca

Ca

Ca

By induction it follows that

∥∇ Uk+1 ∥2 + ∆t ≤

Using

C1 Ca

∥∆Uk+1 ∥2 + ∆t

Cb Ca

∥∇ ∆Uk+1 ∥2

k−1 ( ( C C )k ( ) ∑ Cb Cc Cd )i 2 C1 c d ∥∇ U0 ∥2 + ∆t ∥∆U0 ∥2 + ∆t ∥∇ ∆U0 ∥2 + ∆t ∆tC (Ω , D, λ0 , f )

Cc Cd Cc

Cc

Ca

Ca

i=0

Cc

Ca

= 1 + K3 ∆t and the inequality (1 + x)k ≤ exk for any k∆t ≤ T we have C1

Cb

∥∇ Uk ∥2 + ∆t ∥∆Uk ∥2 + ∆t ∥∇ ∆Uk ∥2 Ca Ca ( ) C1 Cb K3 T 2 ≤e ∥∇ U0 ∥ + ∆t ∥∆U0 ∥2 + ∆t ∥∇ ∆U0 ∥2 + ∆tTC (Ω , D, λ0 , f ) , Ca

Ca

which is of the form (25) with the choice K1 =

C1 Ca

, K2 =

Cb . Ca

□

Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

8

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

To prove (c) we need the following two lemmas. Lemma 1. The error ek = uk − Uk between the exact and approximate solution satisfies

∫ Ω

ek dx = O(∆t 2 )

Proof. Proof is exactly similar to the proof of Lemma 3.2 of [19].

□

Lemma 2. Under the assumption that the exact solution uk of (10) is smooth at time t = k∆t and let T > 0. Then there exists a constant C > 0 such that ∇ uk ≤ C for all t = k∆t < T . Proof. Proof is exactly similar to the proof of Lemma 3.3 of [19] . □ Proof. (c) Putting value of ut (k∆t) in (27) we get the truncation error as uk+1 − uk

τk =

∆t

(

+ ϵ∆ ∇ · √

) )

∇ uk

(

|∇ uk | + δ 2

1

− ∆W ′ (uk ) + λ(f − uk ) + C1 ∆2 (uk+1 − uk ) ϵ

2

− C1 ∆(uk+1 − uk ) + C2 (uk+1 − uk )

(34)

Subtracting (21) from (34) and rearranging the terms we get that the discretization error ek satisfies ek+1 − ek

∆t

1

+ C1 ∆2 ek+1 − C1 ∆ek+1 + C2 ek+1 =C1 ∆2 ek − C1 ∆ek + C2 ek + (W ′ (uk ) − W ′ (Uk )) − λek ϵ ( ) ∇ uk ∇ Uk + ϵ∆ ∇ · ( √ ) − ∇ · (√ ) + τk |∇ uk |2 + δ 2 |∇ Uk |2 + δ 2

Now multiplying by −∆ek+1 and integrating over Ω we have, 1 (

∆t

) ∥∇ ek+1 ∥2 − ⟨∇ ek , ∇ ek+1 ⟩ + C1 ∥∇ ∆ek+1 ∥2 + C1 ∥∆ek+1 ∥2 + C2 ∥∇ ek+1 ∥2

⟨ ( = ϵ ∆∇ · √

∇ uk |∇ uk |2 + δ 2

−√

∇ Uk

) ⟩ 1 , ∆ek+1 − ⟨∆(W ′ (uk ) − W ′ (Uk )), ∆ek+1 ⟩ + C1 ⟨∇ ∆ek , ∇ ∆ek+1 ⟩ ϵ |∇ Uk |2 + δ 2

+ C1 ⟨∆ek , ∆ek+1 ⟩ + C2 ⟨∇ ek , ∇ ek+1 ⟩ − ⟨∇λek , ∇ ek+1 ⟩ + ⟨∇ ∆−1 τk , ∇ ∆ek+1 ⟩ Applying Young’s inequality we get, 1 2∆t

(∥∇ ek+1 ∥2 − ∥∇ ek ∥2 ) + C1 ∥∇ ∆ek+1 ∥2 + C1 ∥∆ek+1 ∥2 + C2 ∥∇ ek+1 ∥2

⟨ ( ≤ −ϵ ∇∇ · √

∇ uk |∇ uk |2 + δ 2

+ C1 δ1 ∥∇ ∆ek+1 ∥2 + +

C1

δ1

−√

∇ Uk

)

|∇ Uk |2 + δ 2

∥∇ ∆ek ∥2 +

C1 2

⟩ 1 , ∇ ∆ek+1 + ⟨∆(W ′′ (uk )∇ uk − W ′′ (Uk )∇ Uk ), ∇ ∆ek+1 ⟩ ϵ

∥∆ek ∥2 +

C1 2

∥∆ek+1 ∥2 + C2 δ2 ∥∇ ek+1 ∥2 +

C2

δ2

∥∇ ek ∥2 + δ3 ∥∇ ek+1 ∥2

λ20 1 δ4 ∥∇ ek ∥2 + ∥τk ∥2−1 + ∥∇ ∆ek+1 ∥2 δ3 2δ4 2

(35)

Following the similar argument as given in [19] and using Lemmas 1 and 2 we get the following inequality for the first and second term of RHS of (35) i.e.

⟨∆(W ′′ (uk )∇ uk − W ′′ (Uk )∇ Uk ), ∇ ∆ek+1 ⟩ ≤

( C C ) δ5 + δ6 + ∥∇ ek ∥2 + ∥∇ ∆ek ∥2 + ∥O(∆t)2 ∥2 2 δ5 2δ6 2

(36)

Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

⟨ ( ∇∇ · √ ≤

C

∇ uk |∇ uk | + δ 2 2

−√

∇ Uk |∇ Uk |2 + δ 2

3

2δ¯

9

⟩ ) , ∇ ∆ek+1 5

L

¯ ∆ek+1 ∥2 + BC ( + )∥∇ ∆ek ∥2 ∥∇ ∆ek ∥2 + CL + ∥∇ ek ∥2 + 6δ∥∇ 2δ¯ δ¯ δ¯

(37)

Using (36) and (37) in Eq. (35) and re arranging the term we get

) ( ) ( 1 C1 δ4 δ5 + δ5 + C2 (1 − δ2 ) − δ3 ∥∇ ek+1 ∥2 + ∥∆ek+1 ∥2 + C1 (1 − δ1 ) − −C − 6δ¯ ∥∇ ∆ek+1 ∥2 2∆t 2 2 2ϵ ( 1 ) ( 2 λ C2 C C 3ϵ CL C1 Cϵ 5ϵ BC ϵ BCL ) C1 ≤ + + 0 + + + ∥∇ ek ∥2 + + + + ∥∇ ∆ek ∥2 + ∥∆ek ∥2 ¯δ ¯ ¯ ¯ 2∆t δ2 δ3 2ϵδ5 2ϵδ6 δ1 2 2δ 2δ δ +

1

δ4

∥τk ∥2−1 + C ∥O(∆t)2 ∥2

(38)

Multiplying (38) by 2∆t and setting Ca = 1 + 2∆t(C2 (1 − δ2 ) − δ3 ), cc = 1 + 2∆t

(C

2

δ2

+

) δ4 δ5 + δ5 C − 6δ¯ 2 2ϵ (C 5ϵ BC ϵ BCL ) Cϵ 1 + + Cd = + ¯ ¯ δ1 2δ 2δ δ¯

(

Cb = 2 C1 (1 − δ1 ) −

λ20 C C 3ϵ CL ) + + + , δ3 2ϵδ5 2ϵδ6 δ¯

we have Ca ∥∇ ek+1 ∥2 + ∆tC1 ∥∆ek+1 ∥2 + ∆tCb ∥∇ ∆ek+1 ∥2

≤ Cc ∥∇ ek ∥2 + ∆tC1 ∥∆ek ∥2 + ∆tCd ∥∇ ∆ek ∥2 + ∆t

(1 δ4

∥τk ∥2−1 + C ∥O(∆t)2 ∥2

)

Now choose all δ s such that all the coefficients of (38) the inequality are positive and Dividing both sides of (39) by Ca we get

∥∇ ek+1 ∥2 + ∆t

C1 Ca

∥∆ek+1 ∥2 + ∆t

Ca

Cc Cb Ca

, CCc Ca d greater than one.

∥∇ ∆ek+1 ∥2

) ∆t ( 1 ) ∥∇ ∆ek ∥2 + ∥τk ∥2−1 + C ∥O(∆t)2 ∥2 Ca Cd Cc Cd Ca C a δ4 ( ) ) Cc Cd C1 Cb ∆t ( 1 ≤ ∥∇ ek ∥2 + ∆t ∥∆ek ∥2 + ∆t ∥∇ ∆ek ∥2 + ∥τk ∥2−1 + C ∥O(∆t)2 ∥2 Ca Ca Ca C a δ4 ≤

Cc Cd ( 1

∥∇ ek ∥2 +

∆tC1

Cb

(39)

∥∆ek ∥2 +

∆t

So by induction on k we get

∥∇ ek+1 ∥2 + ∆t

≤

( Cc Cd )k (

≤

∆t

Ca

Ca

C1 Ca

∥∆ek+1 ∥2 + ∆t

∥∇ e0 ∥2 +

keK3 k∆t

(1 δ4

C1 Ca

Ca

Ca

∥∇ ∆ek+1 ∥2

∥∆e0 ∥2 +

∆tCb Ca

k ) ∆t ∑ ) ( Cc Cd )i ( 1 ∥∇ ∆e0 ∥2 + ∥τk ∥2−1 + C ∥O(∆t)2 ∥2 Ca Ca δ4 i=0

max ∥τi ∥2−1 + C ∥O(∆t)2 ∥2

)

i≤k

where we used the fact that

∥∇ ek+1 ∥2 + ∆t

∆tC1

Cb

Cc Cd Ca

= 1 + M4 ∆t and e0 = 0. Hence for k∆t < T , we have

∥∆ek+1 ∥2 + ∆t

Cb Ca

∥∇ ∆ek+1 ∥2 ≤

which is of the form (26) with the choice M1 =

C1 Ca

T Ca

, M2 =

eM4 T O(∆t)2 ,

Cb Ca

and M3 = Ca .

□

Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

10

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 1. Left to right: Image to be inpainted, Inpainting results of CH model, TV − H −1 and Our binary model.

Fig. 2. Left to right: Image to be inpainted, Inpainting results of CH model, TV − H −1 and Our binary model.

Fig. 3. Left to right: Image to be inpainted, Inpainting results of CH model, TV − H −1 and Our binary model.

Fig. 4. Left to right: Image to be inpainted, Inpainting results of CH model, TV − H −1 and Our binary model.

Fig. 5. Left to right: Image to be inpainted, Inpainting results of CH model, TV − H −1 and Our binary model.

5. Numerical results In this section firstly we will present the results of our model (10) for the binary images and compare results with the results of Cahn–Hilliard (CH) model [9] and TV − H −1 model [6]. After that, we present the results of our model (15) for grayscale image and compare with the results of TV − H −1 and TV − L2 model. We will compare our results with the results of other models in terms of quality measures PSNR and SNR and SSIM [18] of resulting images. Now we will introduce the above mentioned measures. Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

11

Fig. 6. Left to right: Image to be inpainted, Inpainting results of CH model, TV − H −1 and Our binary model.

Fig. 7. Left to right: Image to be inpainted, Inpainting results of CH model, TV − H −1 and Our binary model.

Fig. 8. First row contains the original image and the image to be inpainted and second row contains Inpainting results of CH model, TV − H −1 and Our binary model respectively.

(a) Peak Signal to noise ratio (PSNR) is the ratio of maximum possible value of the original image and the mean squared error between the original and the resulting image. The formula for PSNR is as follows:

∑row ∑col i=1

PSNR = 10log10 ∑row ∑col i=1

j=1

j=1

2552

[O(i, j) − I(i, j)]2

where O is the original image and I is the recovered image. Higher PSNR indicates the better the processed image. (b) Signal to noise ratio (SNR) is defined as the ratio of the variance of the original image and recovered image. Higher SNR implies the better inpainting. (c) Structural Similarity Index (SSIM) is defined to measure the similarity between two images [18]. SSIM is defined as follows: SSIM(O, I) =

(2µO µI + c1 )(2σOI + c2 ) (µ2O + µ2I + c1 )(σO2 + σI2 + c2 )

Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

12

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 9. First row contains the original logoiImage and the image to be inpainted. Second row contains inpainting results of CH model, TV − H −1 and Our binary model respectively.

Fig. 10. First row: Original Image and damaged image. Second row: Inpainting results of CH model, TV − H −1 and Our binary model respectively.

where µ, σ represents mean and variance of the corresponding image and c1 = (0.01L)2 , c2 = (0.03L)2 , here dynamic range L = 255. Higher the SSIM implies the better recovery of the image. The parameters are: the coefficient of 4th-order term ϵ , δ is used to regularize the curvature term in (9), λ is fidelity parameter defined in (7), ∆t , ∆x, ∆y are the step sizes in the corresponding direction. The constant C1 , C2 comes from the convexity splitting. For our case C1 > max{ 1ϵ , 1δ }, for CH model C1 > 1ϵ and for TV − H −1 model C1 > 1δ and C2 > λ0 . For the parameters of CH model we have followed the thesis [22] and for TV − H −1 model we have followed [23]. For CH and TV − H −1 model, λ is varying image to image and for our model, we vary λ and C1 . From the expression of the λ in (7) it is clear that varying λ is equivalent to varying λ0 . So in the rest of the paper we will mention only λ0 . These Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

13

Table 1 Table for binary image. Image

Model

PSNR

SNR

SSIM

Iteration

CPU time

parameters

Triangle

CH TV − H −1 Our

24.52 24.29 27.78

23.05 22.83 26.32

0.9550 0.91772 0.9841

2900 + 4700 15 900 12 350 + 350

10.69 49.41 45.56

λ0 = 1e6 λ0 = 10 λ0 = 1, C1 = 100

Strip

CH TV − H −1 Our

32.09 36.71 62.11

31.22 35.85 61.24

0.9848 0.9877 1.0000

2750 + 4150 4000 450 + 100

8.95 9.06 1.48

λ0 = 1e8 λ0 = 10 λ0 = 1, C1 = 10

Single strip

CH TV − H −1 Our

32.61 37.11 63.18

24.76 29.26 55.33

0.9749 0.9872 1.00

9850 + 6000 6850 500 + 100

22.25 20.26 2.19

λ0 = 1e8 λ0 = 10 λ0 = 1, C1 = 10

Double strip

CH TV − H −1 Our

27.15 35.58 60.17

22.2 30.64 55.22

0.9711 0.9669 1.00

2050 + 2250 12 950 3350 + 100

6.39 38.89 12.18

λ0 = 1e7 λ0 = 10 λ0 = 1, C1 = 10

Large strip Fig. 5

CH TV − H −1 Our

26.34 30.12 63.18

21.50 25.28 58.34

0.9648 0.9329 1.00

1550 + 1700 17 550 3450 + 200

4.68 52.52 13.18

λ0 = 1e6 λ0 = 10 λ0 = 1, C1 = 100

Cross

CH TV − H −1 Our

23.28 24.57 26.85

17.88 19.17 21.44

0.9375 0.9732 0.9843

50 + 250 3350 400 + 50

0.44 9.89 1.67

λ0 = 5e5 λ0 = 10 λ0 = 1, C1 = 50

Circle

CH TV − H −1 Our

24.46 26.84 32.04

22.96 25.34 30.54

0.9512 0.9625 0.9948

50 + 1350 5700 1900 + 100

2.15 17.52 7.38

λ0 = 1e5 λ0 = 10 λ0 = 1, C1 = 10

IIT Fig. 8

CH TV − H −1 Our

19.82 19.41 20.54

12.07 11.67 12.80

0.9007 0.9370 0.9390

200 + 150 9550 1750 + 100

0.48 28.77 6.43

λ0 = 1e7 λ0 = 100 λ0 = 100, C1 = 300

logo

CH TV − H −1 Our

25.90 27.79 28.09

16.89 18.79 19.09

0.9776 0.9843 0.9862

900 + 9700 9550 4150 + 100

48.39 81.33 42.97

λ0 = 1e8 λ0 = 200 λ0 = 200, C1 = 200

Football

CH TV − H −1 Our

27.58 29.34 30.50

25.85 27.61 28.77

0.9864 0.9912 0.9937

650 + 2750 4000 1300 + 300

14.46 33.00 15.75

λ0 = 1e8 λ0 = 100 λ0 = 100, C1 = 200

parameters which are not fixed are reported in the table in the column bearing the name parameters. The parameters which are fixed for all the images for all the three models are mentioned in the next paragraph. The parameters for our model are as follows: ϵ = 0.9 in 1st stage and 0.3 in the 2nd stage, δ = 0.1, ∆t = 100, ∆x = ∆y = 1 and the constants C2 = λ0 . The parameters used (from [22]) for CH inpainting are: ϵ = 0.8 in the 1st phase and .01 in the 2nd phase the constants C1 = 300 and C2 = 3λ0 and ∆t = 1, ∆x = ∆y = 0.01. For TV − H −1 model, δ = .01, ∆t = 1, ∆x = ∆y = 1 and the constants C1 = 1δ , C2 = λ0 as set in [23]. In the table, we have mentioned the number of iterations in the format: iterations of first phase + iterations second phase, for CH and our model. The convergence criterion for all the models are based on the convergence of PSNR. The stopping criterion for all the models ensures that the variation in PSNR is below the tolerance 10−4 . 5.1. Results for binary image In Fig. 1 we have shown the inpainted result of a triangle by our model and compared the result with the result obtained by applying CH inpainting model and TV − H −1 inpainting model. The proposed model effectively restores the original image including edges but the other methods get less effective w.r.t edges. While TV − H −1 model generates edge with a blur especially in the inpainted portion, the CH method leads to a blurred edge. In Table 1 we have reported the PSNR, SNR, and SSIM of each result. Also we have reported the variable parameter values in parameters column. From Table 1 it is clear that our model gives better result than other two models in terms of PSNR, SNR, and SSIM and also we can visually see the difference between the images, In Figs. 2–5 we have taken corrupted horizontal and vertical strips of different shape and sizes and compared the results of our model with other models. The parameter values together with the values of PSNR, SNR, and SSIM for each of the resulting images are reported in Table 1. From the table, it is clear that our model gives better result than other two models in terms of the calculated metrics. Also, our model takes less time than TV − H −1 model. Visually also, our model gives better result than other two models. In Fig. 6 we have considered a cross image and shown the inpainting results. From the figure, one can see that the result given by CH model has blurry edges and the result of TV − H −1 keeps some amount of blurry portion especially in the inpainted domain and both these results are inferior to our result. In Table 1 we have reported the PSNR, SNR, Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

14

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 11. First row: Original and damaged image. Second row — results of TV − L2 , TV − H −1 , and Our model (15). In third row, the upper left part of each result has been visualized.

and SSIM for this image. From the table it is clear that our model gives a better result than other two models in terms of PSNR, SNR, and SSIM. Fig. 7 we have considered an image of a circle. In the result of CH model, edges are blurred and TV − H −1 model keeps some amount of blur in the inpainted domain but our model has given a clear result. PSNR, SNR, SSIM for all the three models have been reported in Table 1. Form the table one can see that our model gives a better result than other two models in terms of PSNR, SNR, SSIM. In Fig. 8 we have considered an image which contains the writing ‘IIT’ and some rectangular box of different width is considered. For this case also our model gives a better result than other two models in terms of PSNR, SNR, and SSIM. The details are reported in Table 1. From Fig. 8 one can also note that TV − H −1 model truncates the original image vividly seen in alphabet ‘I’. In order to explore the anisotropic capability of the proposed model, when the images are corrupted in a non-cartesian manner a set of different standard corrupted meaningful images have been considered. Results from Figs. 9–10 clearly portray the overall superiority of the proposed model in comparison to TV − H −1 and CH model. Table 1 contains the PSNR, SNR and SSIM of these images and all the models. From the table, one can see that our model gives better result than other two models in terms of PSNR, SNR, and SSIM. In the logo image visually difference is clear from the racket. 5.2. Results for grayscale image with multi-well potential Now, we compare the results of our multi-well potential model given by Eq. (15) with the result of TV − H −1 and TV − L2 model. For this case, we choose all the parameters of our model same as mentioned in the beginning of the current section except ϵ = .1 in the second phase. The constant C1 = 10 for all images, and the parameter α used in Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

15

Fig. 12. First row contains the original image of an elephant and a damaged version of it and the second row contains the inpainted results of TV − L2 , TV − H −1 , and Our model (15) respectively. In third row the lower right portion of each result has been shown. Table 2 Table for gray image. Image

Model

PSNR

SNR

SSIM

Iteration

CPU time

parameters

Grayshade

TV − L2 TV − H −1 Our

41.64 45.23 59.52

37.94 41.53 55.82

0.9822 0.9909 0.9997

3550 6650 1050 + 250

43.16 87.59 22.17

λ0 = 10 λ0 = 10 λ0 = 10, C1 = 10

Kaleidoscope

TV − L2 TV − H −1 Our

26.40 27.51 28.27

18.67 19.78 20.54

0.8965 0.9019 0.9274

10 000 3700 3000 + 250

233.51 92.85 108.84

λ0 = 100 λ0 = 100 λ0 = 50

Einstein 1

TV − L2 TV − H −1 Our

24.50 28.01 28.47

18.05 21.57 22.02

0.8070 0.8519 0.8772

8750 7800 4250 + 100

177.59 165.62 134.25

λ0 = 100 λ0 = 100 λ0 = 50

Einstein 2

TV − L2 TV − H −1 Our

22.68 25.59 25.81

16.24 19.14 19.25

0.7395 0.7990 0.8242

6950 7100 4650 + 150

152.66 151.45 140.02

λ0 = 100 λ0 = 100 λ0 = 50

Elephant 1 Fig. 12

TV − L2 TV − H −1 Our

27.36 28.04 28.62

20.54 21.23 21.81

0.8727 0.8685 0.9097

1200 1300 250 + 100

15.53 17.27 5.89

λ0 = 100 λ0 = 100 λ0 = 50

Elephant 2 Fig. 16

TV − L2 TV − H −1 Our

33.82 33.57 36.24

27.01 26.75 29.43

0.9500 0.9394 0.9804

3400 4000 1950 + 150

41.05 51.50 35.16

λ0 = 100 λ0 = 100 λ0 = 50

Dog

TV − L2 TV − H −1 Our

38.83 38.66 41.49

31.78 31.62 34.45

0.9897 0.9835 0.9910

7900 8150 1750 + 1800

84.92 94.81 50.39

λ0 = 1000 λ0 = 1000 λ0 = 200

Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

16

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 13. First row contain the original image of a kaleidoscope and a damaged version of it and the second row contains the inpainted results of TV − L2 , TV − H −1 , and Our model (15) respectively. In third row the lower right portion of each result has been shown.

Table 3 Comparison of results of Lena and Barbara image. Image

Model

PSNR

SNR

SSIM

Iteration

CPU time

parameters

Lena 1

TV − L2 TV − H −1 Our

31.57 32.87 33.26

25.92 27.22 27.61

0.9165 0.9100 0.9488

1200 1900 550 + 100

48.75 82.93 36.99

λ0 = 100 λ0 = 100 λ0 = 50

Lena 2

TV − L2 TV − H −1 Our

34.84 34.60 35.86

29.19 28.94 30.20

0.9456 0.9271 0.9735

2500 4900 3150 + 100

107.93 220.62 199.34

λ0 = 100 λ0 = 100 λ0 = 50

Lena 3

TV − L2 TV − H −1 Our

26.79 29.28 29.38

21.11 23.59 23.70

0.8795 0.8818 0.9180

20 000 20 000 7900 + 1550

818.30 850.05 547.88

λ0 = 100 λ0 = 100 λ0 = 50

Barbara 1

TV − L2 TV − H −1 OUr

26.92 26.93 28.26

21.04 21.05 22.37

0.8851 0.8516 0.9160

1650 1450 600 + 100

67.38 64.80 42.80

λ0 = 100 λ0 = 100 λ0 = 50

Barbara 2

TV − L2 TV − H −1 Our

29.87 28.65 31.70

23.98 22.76 25.82

0.9346 0.8872 0.9613

2000 4350 3300 + 100

81.37 206.74 204.44

λ0 = 100 λ0 = 100 λ0 = 50

Barbara 3

TV − L2 TV − H −1 Our

23.12 23.76 26.45

17.24 17.88 20.56

0.6802 0.6871 0.8571

9650 88 000 1900 + 1900

375.27 389.07 214.83

λ0 = 10 λ0 = 10 λ0 = 10

Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

17

Fig. 14. First row contains the original image and two damaged image. Second and third row contain the inpainted result of the above damaged image obtained by TV − L2 , TV − H −1 , and Our model (15) respectively.

Fig. 15. First row contains original images of a dog and a damaged image of it, second row contains inpainting results of dog by TV − L2 , TV − H −1 , and Our model (15) respectively.

Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

18

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 16. First row contains original images of an elephant and a damaged image of it, second row contains inpainting results of dog by TV − L2 , TV − H −1 , and Our model (15) respectively.

the formula (12) of generalized histogram is set as α = .01 except the first image where we have used α = 0.001. For TV − H −1 model we have used the same parameters except λ0 which varies image to image. For TV − L2 model also the convexity splitting method is used and whose code is available in [24] and the parameters are chosen same as TV − H −1 model. In Fig. 11 we have shown the inpainted result of a gray shade image. The gray shade image has intensity variation with sharp edges and we have taken strips of different size as inpainting domain. In the last row we have shown the upper left part of each result. From the figure, we can see that our model recovers the edges with out any blur where as TV − H −1 has little blurry edge. In Table 2 we have reported PSNR, SNR, SSIM for both the images. From the table, it is clear that our model gives a better result than the TV − H −1 and TV − L2 model. Also, notice that SSIM gives the similarity between original image and resulting image and for both the cases our SSIM is highest i.e. our result is more closer to the original image. So our model gives better result visually as well as in terms of PSNR, SNR and SSIM. In Fig. 12 we have presented the inpainting result of an elephant whose original and damaged version are included in the first row. In the last row of Fig. 12 we have shown the lower right part of each result presented in second row. From the figure, one can see that the results of TV − H −1 and TV − L2 model are more blurry than our result i.e. visually our model gives better result than the other two. Also Table 2 contains the PSNR, SNR, and SSIM values of these two images for all the three models. From the table, one can see that our model gives the higher PSNR, SNR, and SSIM which implies that our model gives the better result than the rest. In Fig. 13 we have presented the inpainting result of a kaleidoscope image whose original and damaged version are included in the first row. In the last row of Fig. 13 we have shown the lower right part of each result presented in second row. From the figure and Table 2 one can see that our model gives better result than other two models. In Fig. 14 we have taken the image of Einstein and took two different inpainting domain. We have presented the results of our model and compared with the results of TV − H −1 and TV − L2 model. From the figure, one can see that TV − L2 is not giving a good result but TV − H −1 and our model gives a good result. Our model gives the better result than the TV − H −1 model which is clear from the image itself (see the left cheek, eye). PSNR, SNR, and SSIM are reported in Table 2 for all the models. From the table, one can see that our model results in highest PSNR, SNR, and SSIM i.e. our model is better than other two models in terms of PSNR, SNR and SSIM. Also notice that there is a little difference between the PSNR values of TV − H −1 model and Our model but the difference in SSIM is more than 2 which implies that our result is more close to the original image. In Figs. 15–16 we have presented inpainting results of two images which are similar to those presented in [23]. For these two images, we have taken the parameters same as taken by the author of [23]. From the figure, one can notice that our model gives better result than other two models. Our model keeps the same brightness and the image is effected only on the inpainting domain but other two models smooth the whole image and hence less bright. We get this advantage because of the two-step process. In the first step, the image is repaired and 2nd step helps to improve the brightness of the image outside the inpainting domain. In Table 2 we have reported the PSNR, SNR, and SSIM of all the three models. From the table, one can see that there is huge difference in PSNR and SNR of our model and the other two. Also, SSIM for our model is higher than the rest which implies that our model results in an image which is more closer to the original image. Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

19

Fig. 17. First column contains three different damaged image of Lena image and rest of the rows contains the results of the corresponding image obtained by TV − L2 , TV − H −1 , and Our model (15) respectively.

In Fig. 17 we have taken different types of inpainting domain to the Lena image of size 512 × 512 and have shown the results of our model. Also, compare our result with the result of TV − H −1 model and TV − L2 model. In Table 3 we have reported the PSNR, SNR, and SSIM of the inpainted image. Our model is found to be superior not only w.r.t both the image quantity metrics but also w.r.t total computational time. In Fig. 18 we have taken inpainting domain of different size of the Barbara image of size 512 × 512 and have shown the result of our model and also compare the result with the result of TV − H −1 and TV − L2 model. In Table 3 we summarize the output results in terms of PSNR, SNR, and SSIM and again find that our model is superior to both TV − H −1 and TV − L2 model. 6. Conclusion Here we have presented a new 4th order anisotropic PDE model for grayscale images based on a multi-well potential using the image specific histogram data. The proposed PDE based model for image inpainting is superior to both TV − L2 Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

20

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 18. First column contains three different damaged image of Barbara image and rest of the cows contains the inpainting results of the corresponding image obtained by TV − L2 , TV − H −1 , and Our model (15) respectively.

model and TV − H −1 model. The introduction of histogram-based multi-well potential enables the cost-effective approach to restore grayscale images and the anisotropic 4th order diffusion term helps in effective handling of edges. The convexity splitting approach in conjunction with discrete Fourier method for solving PDE gifts an unconditionally stable scheme which is accurate, economical and robust. References [1] A. Criminisi, P. Perez, K. Toyama, Object removal by exemplar-based inpainting, IEEE Int. Conf. Comp. Vis. Pattern Recognit. 2 (2003) 721–728. [2] A.A. Efros, T.K. Leung, Texture synthesis by non-parametric sampling, in: IEEE International Conference on Computer Vision, Corfu, Greece, September, 1999. [3] Xin Li, Image recovery via hybrid sparse representations: A deterministic annealing approach, IEEE J. Sel. Top. Sign. Proces. 5 (5) (2011). [4] T.F. Chan, J.H. Shen, H.M. Zhou, Total variation wavelet inpainting, J. Math. Imaging Vision 25 (1) (2006) 107–125. [5] J.A. Dobrosotskaya, A.L. Bertozzi, A wavelet-laplace variational technique for image de-convolution and inpainting, IEEE Trans. Image Process. 17 (5) (2008) 657–663.

Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.

A. Halim and B.V. Rathish Kumar / Computers and Mathematics with Applications xxx (xxxx) xxx [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]

21

M. Burger, L. He, C. Schönlieb, Cahn-Hilliard inpainting and a generalization for gray value images, SIAM J. Imaging Sci. 2 (4) (2009) 1129–1167. T.F. Chan, J. Shen, Mathematical models for local non-texture inpaintings, SIAM J. Appl. Math. 62 (3) (2001) 1019–1043. T.F. Chan, J. Shen, Non-texture inpainting by curvature driven diffusions (CDD), J. Vis. Commun. Image Rep. 12 (4) (2001) 436–449. A. Bertozzi, S. Esedoglu, A. Gillette, Inpainting of binary images using the Cahn-Hilliard equation, IEEE Trans. Image Process. 16 (1) (2007) 285–291. M. Bertalmio, G. Sapiro, V. Caselles, C. Ballester, Image inpainting, in: Siggraph 2000, Computer Graphics Proceedings, 2000, pp. 417–424. M. Bertalmio, A.L. Bertozzi, G. Sapiro, Navier–Stokes, fluid dynamics, and image and video inpainting, in: Proceedings of the 2001 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, Vol. 1, 2001, pp. 355–362. L.I. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D 60 (1992) 259–268. P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Trans. Pattern Anal. Mach. Intell. 12 (7) (1990) 629–639. T.F. Chan, S.H. Kang, J. Shen, Eulers elastica and curvature-based inpainting, SIAM J. Appl. Math. 63 (2) (2002) 564–592. M. Nitzberg, D. Mumford, T. Shiota, Filtering, Segmentation, and Depth, in: Lecture Notes in Computer Science, Springer-Verlag, 1993, p. 662. S. Esedoglu, J.H. Shen, Digital inpainting based on the Mumford-Shah-Euler image model, Eur. J. Appl. Math. 13 (4) (2002) 353–370. A. Bertozzi, S. Esedoglu, A. Gillette, Analysis of a two-scale Cahn-Hilliard model for image inpainting, Multiscale Model. Simul. 6 (3) (2007) 913–936. Zhou Wang, A.C. Bovik, H.R. Sheikh, E.P. Simoncelli, Image quality assessment: from error visibility to structural similarity, IEEE Trans. Image Process. 13 (4) (2004) 600–612. C.B. Schönlieb, A. Bertozzi, Unconditionally stable schemes for higher order inpainting, Commun. Math. Sci. 9 (2) (2011) 413–457. D. Eyre, An unconditionally stable one-step scheme for gradient systems, 1998, unpublished. Rowthu Vijayakrishna, A Unified Model of Cahn-Hilliard Greyscale Inpainting and Multiphase Classification (Ph.D. thesis), IIT Kanpur, India, 2015. Alan Gillette, Image Inpainting using a Modified Cahn-Hilliard Equation (Ph.D. thesis), University of California, Los Angeles, 2006. Carola-Bibiane Schönlieb, Modern PDE Techniques for Image Inpainting (Ph.D. thesis), DAMTP, Centre for Mathematical Sciences University of Cambridge, 2009. C.B. Schönlieb, Higher-order total variation inpainting, File Exchange - MATLAB Central.

Please cite this article as: A. Halim and B.V. Rathish Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.12.002.