This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

∎

11institutetext:
Department of Applied Mathematics
King Abdullah University of Science and Technology, Saudi Arabia22institutetext: 22email: [email protected] 33institutetext: Department of Mathematics, Raiganj University, West Bengal, India
33email: [email protected]
44institutetext: Department of Mathematics and Statistics , IIT Kanpur, Uttar Pradesh, India
44email: [email protected] 55institutetext: Department of Mathematics, Raganj University, India

A Fractional Image Inpainting Model Using a Variant of Mumford-Shah Model

Abdul Halim    Abdur Rohim    B.V. Rathish Kumar    Ripan Saha
Abstract

In this paper, we propose a fourth order PDE model for image inpainting based on a variant of the famous Mumford-Shah (MS) image segmentation model. Convexity splitting is used to discrtised the time and we replace the Laplacian by its fractional counterpart in the time discretised scheme. Fourier spectral method is used for space discretization. Consistency, stability and convergence of the time discretised model has been carried out. The model is tested on some standard test images and compared them with the result of some models existed in the literature.

Keywords:
Image Processing Inpainting Fourth Order Partial differential equation (PDE)Convexity splitting.

1 Introduction

Image inpainting is a process of filling in the missing part of an image from the information available outside the missing part. In mathematical term it is an extrapolation of the image. Although there are several inpainting methods are available in the literature, we consider the Partial differential equation (PDE) based methods because of its sound established theory.

Since image inpainting is relatively new topic in the field of image processing, a lot of inpainting models are based on the models of denoising and segmentation. Like in CS1 Chan and Shen have used a modified version of Rudin-Osher-Fatemi (ROF) denoising model rof for image inpainting and the resulting inpainting model is known as T​V−L2TV-L^{2} inpainting model. Following the T​V−L2TV-L^{2} model a series of models sch ; CS2 ; CS3 have been proposed which fall in the category of variational inpainting. In variational inpainting, the PDE is obtained by minimizing an energy functional of the following form:

E​(u)=R​(u)+λ​F​(f−u)E(u)=R(u)+\lambda F(f-u) (1)

where ff is the given damaged image and uu is the inpainted image and

λ​(x)={λ0​i​n​Ω∖D0​i​n​D\lambda(x)=\begin{cases}\lambda_{0}\hskip 2.84544ptin\hskip 5.69046pt\Omega\setminus D\\ 0\hskip 5.69046ptin\hskip 5.69046ptD\end{cases} (2)

with λ0>>1\lambda_{0}>>1 and DD is the damaged part of the image. The term R​(u)R(u) is called regularization term, F​(⋅)F(\cdot) is called fidelity term and λ\lambda is called fidelity parameter which forces the inpainting image to remain closer to the original image outside the missing/ degraded region DD.

Another variational model for inpainting has been proposed by Esedoglu and Shen mse . They have brought the Mumford-Shah model of segmentation and used it for image inpainting. The Mumford-Shah(MS) functional is non-convex and contains an unknown set of lower dimension, so it is difficult to solve the minimization problem. Several approximations have been proposed in the literature surveyMS . Esedoglu and Shen mse have used the Ambrosio and Tortorelli’s approximation of the MS energy functional for the inpainting model. Applying the gradient descent to the approximations give raise to a parabolic problem of second order with a parameter ϵ\epsilon. Since the second order models are not good for inpainting tvh , as second order models are unable to fill large gaps, they have modified their second order model using Euler-Elastica energy to get a fourth order model in the same paper mse , known as Mumford-Shah-Euler (MSE) model.

The parameter ϵ\epsilon of the MSE model reminds the authors of mCH about the Cahn-Hilliard equation used to model phase-separation phenomena of material science. They have modified the Cahn-Hilliard equation and proposed a model for image inpainting in mCH . Here after we will call this model as mCH model. The mCH model has few drawbacks, like it is applicable for binary images only. Then several authors have come up with inapinting models cCH ; fCH ; anis ; lininp ; nonlininp by modifying the Cahn-Hilliard model to overcome the drawbacks of mCH model.

More recently, Cai et.al. twophase have proposed a segmentation model which is also derived from the Mumford-Shah model. This segmentation model is also used for inpainting in the subsequent paper of the same authors in  threephase . But this model will leads to second order equation. As we know that the second order model are not able to fill large gaps tvh , so we will propose a fourth order model in this article using a variant of the MS functional. We will use convexity splitting tvh in time and Fourier spectral methods in space for the proposed model.

The paper is organised is as follows: In section 2 we have discussed some relevant existing models and thereby proposed our model. Section 3, talks about convexity splitting and the numerical scheme for our proposed model. Also we have established the consistency, stability and convergence of the time discretised model. Further, we have introduced the fractional time discretised model. In section 4, we have presented the fractional version of the time discretized scheme and present the scheme with complete discretization. Then the numerical results of our model have been presented and compared them with the same of models like T​V−L2TV-L^{2} and T​V−H−1TV-H^{-1} and the second order model of Cai et.al. twophase . Finally, in section 5 we draw some conclusions.

2 Proposed model and some existing models

In this section, we will discuss few inpainting models from the literature which are relevant to us and then proposed our model.

2.1 T​V−L2TV-L^{2} and T​V−H−1TV-H^{-1} model

The inpainting model T​V−L2TV-L^{2} and T​V−H−1TV-H^{-1} are of variational type that is they are obtained by minimizing an energy functionl of the form (1). For both the models the regularization term is the total variation of uu that is R​(u)=∫Ω|∇u|​𝑑xR(u)=\int_{\Omega}|\nabla u|dx. The fidelity term for the T​V−L2TV-L^{2} model is ‖f−u‖2\|f-u\|_{2} and T​V−H−1TV-H^{-1} is ‖f−u‖H−1​(Ω)\|f-u\|_{H^{-1}(\Omega)}.

With the above choices, the minimizing energy of the T​V−L2TV-L^{2} model is :

E​(u)=∫Ω|∇u|​𝑑x+λ2​∫Ω(f−u)2​𝑑x.E(u)=\int_{\Omega}|\nabla u|dx+\frac{\lambda}{2}\int_{\Omega}(f-u)^{2}dx. (3)

The gradient descent of the corresponding Euler-Lagrange equation is

ut=∇⋅∇u|∇u|+λ​(f−u).u_{t}=\nabla\cdot\frac{\nabla u}{|\nabla u|}+\lambda(f-u). (4)

To avoid division by zero the above equation is modified CS1 as:

ut=∇⋅∇u|∇u|2+δ2+λ​(f−u),u_{t}=\nabla\cdot\frac{\nabla u}{\sqrt{|\nabla u|^{2}+\delta^{2}}}+\lambda(f-u), (5)

where δ<<1\delta<<1 is a parameter.

Similarly, we get the PDE for the T​V−H−1TV-H^{-1} model as

ut=−Δ​(∇⋅∇u|∇u|2+δ2)+λ​(f−u),u_{t}=-\Delta\Big{(}\nabla\cdot\frac{\nabla u}{\sqrt{|\nabla u|^{2}+\delta^{2}}}\Big{)}+\lambda(f-u), (6)

2.2 Mumford-Shah model

As we have mentioned in the introduction that our model will be derived from MS-model. In MS- model on has to minimize the following energy functional

EM​S​(u,Γ)=η2​∫Ω(f−u)2​𝑑x+μ2​∫Ω∖Γ|∇u|2​𝑑x+L​e​n​g​t​h​(Γ),E_{MS}(u,\Gamma)=\frac{\eta}{2}\int_{\Omega}(f-u)^{2}dx+\frac{\mu}{2}\int_{\Omega\setminus\Gamma}|\nabla u|^{2}dx+Length(\Gamma), (7)

where η,μ\eta,\mu are parameters and Γ\Gamma is a curve inside the image domain Ω.\Omega. But the functional EM​SE_{MS} is non-convex, so it is not solvable directly. So there are several simplified version has been proposed like proposed in twophase ; mse . A detail survey on Mumford-Shah model and its variant applied to image processing can be found in surveyMS .

In twophase authors have proved that minimizing the functional (7) is equivallent to minimizing the following functional:

E​(u)=η2​∫Ω(f−u)2​𝑑x+μ2​∫Ω|∇u|2​𝑑x+∫Ω|∇u|​𝑑xE(u)=\frac{\eta}{2}\int_{\Omega}(f-u)^{2}dx+\frac{\mu}{2}\int_{\Omega}|\nabla u|^{2}dx+\int_{\Omega}|\nabla u|dx (8)

For segmentation η\eta is constant throughout the domain but for inpainting η\eta is some constant outside the damaged part DD and zero in DD. Changing the fidelity parameter η\eta of the model (8), the same model has been used in threephase in context of colour image which is of the form:

E​(u)=λ2​∫Ω(f−u)2​𝑑x+μ2​∫Ω|∇u|2​𝑑x+∫Ω|∇u|​𝑑xE(u)=\frac{\lambda}{2}\int_{\Omega}(f-u)^{2}dx+\frac{\mu}{2}\int_{\Omega}|\nabla u|^{2}dx+\int_{\Omega}|\nabla u|dx (9)

where λ\lambda is defined in (2). The corresponding Euler-Lagrange equation will be

−λ​(f−u)−μ​Δ​u−∇⋅∇u|∇u|=0.-\lambda(f-u)-\mu\Delta u-\nabla\cdot\frac{\nabla u}{|\nabla u|}=0. (10)

which is equivalent to solve the steady state solution of the equation

ut=μ​Δ​u+∇⋅∇u|∇u|+λ​(f−u).u_{t}=\mu\Delta u+\nabla\cdot\frac{\nabla u}{|\nabla u|}+\lambda(f-u). (11)

Here after we will call this model as convex variant of Mumford-Shah (CVMS) model.

2.3 Proposed Model

But this is a second order model so it will not be able to fill large gap. So we will propose a higher order model. To do this we modify the MS-functional as follows:

EM​S​(u,Γ)=η2​∫Ω(f−u)2​𝑑x+μ2​∫Ω∖Γ(ux1​x12+2​ux1​x22+ux2​x22)​𝑑x+L​e​n​g​t​h​(Γ),E_{MS}(u,\Gamma)=\frac{\eta}{2}\int_{\Omega}(f-u)^{2}dx+\frac{\mu}{2}\int_{\Omega\setminus\Gamma}(u_{x_{1}x_{1}}^{2}+2u_{x_{1}x_{2}}^{2}+u_{x_{2}x_{2}}^{2})dx+Length(\Gamma), (12)

then following the steps of twophase we can prove that this is equivalent to minimising

E​(u)=λ2​∫Ω(f−u)2​𝑑x+μ2​∫Ω|Δ​u|2​𝑑x+∫Ω|∇u|​𝑑x.E(u)=\frac{\lambda}{2}\int_{\Omega}(f-u)^{2}dx+\frac{\mu}{2}\int_{\Omega}|\Delta u|^{2}dx+\int_{\Omega}|\nabla u|dx. (13)

Now using the regularized version of the total-variation the functional reduced to:

E​(u)=λ2​∫Ω(f−u)2​𝑑x+μ2​∫Ω|Δ​u|2​𝑑x+∫Ω|∇u|2+δ2.E(u)=\frac{\lambda}{2}\int_{\Omega}(f-u)^{2}dx+\frac{\mu}{2}\int_{\Omega}|\Delta u|^{2}dx+\int_{\Omega}{\sqrt{|\nabla u|^{2}+\delta^{2}}}. (14)

where δ\delta is a small parameter. For all numerical results presented in this work we choose δ=0.01\delta=0.01. The corresponding descent model will be

ut=−μ​Δ2​u+∇⋅∇u|∇u|2+δ2+λ​(f−u).u_{t}=-\mu\Delta^{2}u+\nabla\cdot\frac{\nabla u}{\sqrt{|\nabla u|^{2}+\delta^{2}}}+\lambda(f-u). (15)
Image
α\alpha
PSNR
SNR
SSIM
Iteration
Time taken
Dog 1 37.50 30.46 0.9412 2999 17.29
1.2 38.78 31.73 0.9704 1580 9.13
1.4 39.14 32.09 0.9767 950 5.41
1.6 39.13 32.09 0.9764 830 4.80
1.8 39.02 31.97 0.9758 810 4.74
2 38.88 31.83 0.9748 720 4.21
GrayShade 1 35.11 31.41 0.7853 1180 9.10
1.2 42.39 38.69 0.9627 1200 9.37
1.4 44.73 41.03 0.9872 1210 9.16
1.6 45.39 41.69 0.9925 1220 9.10
1.8 45.40 41.70 0.9925 1170 9.15
2 45.38 41.67 0.9925 1200 9.42
Table 1: Comparison of results with different fractional power α\alpha, in terms of PSNR, SNR and SSIM.
Image
Model
PSNR
SNR
SSIM
Iteration
parameters
CPU time
GrayShade T​V−L2TV-L^{2} 41.64 37.94 0.9822 33965 λ\lambda=10 18.81
T​V−H−1TV-H^{-1} 45.23 41.54 0.9909 6576 λ\lambda=10 38.84
Our 45.40 41.70 0.9925 1210 α\alpha=1.8 8.90
Kaleidoscope T​V−L2TV-L^{2} 26.24 18.51 0.8827 66365 λ\lambda=50 52.35
T​V−H−1TV-H^{-1} 26.90 19.17 0.8770 2280 λ\lambda=50 22.87
Our 27.89 20.16 0.9130 2020 α\alpha=1.8 23.07
Elephant 1 T​V−L2TV-L^{2} 35.02 28.20 0.9694 2720 λ\lambda=200 12.47
T​V−H−1TV-H^{-1} 33.93 27.12 0.9396 2030 λ\lambda=200 11.41
Our 36.16 29.35 0.9694 490 α\alpha=1.4 3.96
Elephant 2 T​V−L2TV-L^{2} 27.53 20.71 0.9009 2400 λ\lambda=500 10.64
T​V−H−1TV-H^{-1} 28.29 21.47 0.8997 756 λ\lambda=500 4.78
Our 28.51 21.70 0.8961 260 α\alpha=1.8 2.27
Dog T​V−L2TV-L^{2} 38.68 31.63 0.9867 40265 λ\lambda=500 14.37
T​V−H−1TV-H^{-1} 38.13 31.09 0.9758 3946 λ\lambda=500 18.01
Our 39.14 32.09 0.9764 830 α\alpha=1.6 4.80
Table 2: Comparison of Inpainting results with T​V−L2TV-L^{2} and T​V−H−1TV-H^{-1} model on few Real world Images.
Image
Model
PSNR
SNR
SSIM
Iteration
parameters
CPU time
Lena 1 T​V−L2TV-L^{2} 31.58 25.92 0.9167 1091 λ\lambda=100 14.45
CVMS 31.99 26.33 0.9174 645 λ\lambda=100 8.63
Our 32.92 27.27 0.9321 450 α\alpha=1.6 9.57
Lena2 T​V−L2TV-L^{2} 34.86 29.20 0.9456 1990 λ\lambda=100 26.91
CVMS 34.97 29.31 0.9451 835 λ\lambda= 100 11.30
Our 35.62 29.96 0.9597 945 α\alpha=1.4 19.55
Lena3 T​V−L2TV-L^{2} 26.08 20.42 0.7972 11550 λ\lambda =10 159.31
CVMS 27.67 22.01 0.8808 1520 λ\lambda=100 21.32
Our 28.34 22.68 0.8963 1795 α\alpha=1.4 37.60
Barbara 1 T​V−L2TV-L^{2} 26.93 21.04 0.8851 1265 λ\lambda=100 16.74
CVMS 27.70 21.81 0.8867 970 λ\lambda=100 13.10
Our 28.33 22.45 0.9055 485 α\alpha= 1.2 10.40
Barbara 2 T​V−L2TV-L^{2} 29.87 23.98 0.9346 1890 λ\lambda=100 25.72
CVMS 30.56 24.67 0.9340 1020 λ\lambda=100 13.74
OUR 31.40 25.52 0.9505 1055 α\alpha=1.2 22.09
Barbara3 T​V−L2TV-L^{2} 23.10 17.21 0.6797 7660 λ\lambda=10 105.13
CVMS 25.77 19.88 0.8547 1330 λ\lambda=100 17.96
OUR 26.42 20.53 0.8720 1610 α\alpha=1.2 35.52
Table 3: Comparison of Inpainting results with T​V−L2TV-L^{2} and T​V−H−1TV-H^{-1} model on Lena and Barbara Image.

3 Convexity splitting and numerical scheme

We will solve our proposed model (15) using the convexity splitting in time and Fourier spectral method in space. We replaced the Laplacian in the time discretized scheme by its fraction version and get the fractional time discretized model of the proposed model. Convexity splitting methods are used to solve gradient system. Split the energy functional into a convex part and a concave part and in the time discretized scheme the convex part is considered implicitly and the other one explicitly.

The proposed model (15) can be obtained as the gradient flow of the energies E1E_{1} and E2E_{2} in L2L^{2} norm where

E1​(u)=∫Ω(μ2​|Δ​u|2+|∇u|2+δ2)​𝑑x,E2​(u)=∫Ωλ2​(f−u)2​𝑑x.E_{1}(u)=\int_{\Omega}\Big{(}\frac{\mu}{2}|\Delta u|^{2}+\sqrt{|\nabla u|^{2}+\delta^{2}}\Big{)}dx,\quad E_{2}(u)=\int_{\Omega}\frac{\lambda}{2}(f-u)^{2}dx. (16)

To apply convexity splitting idea we split E1E_{1} as E11−E12E_{11}-E_{12} and E2E_{2} as E21−E22E_{21}-E_{22} where

E11=∫Ω(μ2​|Δ​u|2+C12​|∇u|2)​𝑑x and E_{11}=\int_{\Omega}\Big{(}\frac{\mu}{2}|\Delta u|^{2}+\frac{C_{1}}{2}|\nabla u|^{2}\Big{)}dx\quad\text{ and }
E12=∫Ω(C12​|∇u|2−|∇u|2+δ2)​𝑑x,E_{12}=\int_{\Omega}\Big{(}\frac{C_{1}}{2}|\nabla u|^{2}-\sqrt{|\nabla u|^{2}+\delta^{2}}\Big{)}dx,

and

E21=∫ΩC22​|u|2​𝑑x, and E22=∫Ω(−λ2​(f−u)2+C22​|u|2)​𝑑x,E_{21}=\int_{\Omega}\frac{C_{2}}{2}|u|^{2}dx,\quad\text{ and }\quad E_{22}=\int_{\Omega}\Big{(}-\frac{\lambda}{2}(f-u)^{2}+\frac{C_{2}}{2}|u|^{2}\Big{)}dx,

where C1,C2C_{1},C_{2} are positive and should be large enough so that Ei​jE_{ij} for i,j=1,2i,j=1,2 becomes convex.

So the semi-discrete time-stepping scheme for our model is

Uk+1−UkΔ​t=−∇L2(E11​(Uk+1)−E12​(Uk))−∇L2(E21​(Uk+1)−E22​(Uk)).\frac{U_{k+1}-U_{k}}{\Delta t}=-\nabla_{L^{2}}(E_{11}(U_{k+1})-E_{12}(U_{k}))-\nabla_{L^{2}}(E_{21}(U_{k+1})-E_{22}(U_{k})). (17)

Simplifying we get the time discretized model as:

Uk+1−UkΔ​t\displaystyle\frac{U_{k+1}-U_{k}}{\Delta t} +μ​Δ2​Uk+1−C1​Δ​Uk+1+C2​Uk+1\displaystyle+\mu\Delta^{2}U_{k+1}-C_{1}\Delta U_{k+1}+C_{2}U_{k+1}
=∇⋅(∇Uk|∇Uk|2+δ2)−C1​Δ​Uk+λ​(f−Uk)+C2​Uk\displaystyle=\nabla\cdot\Big{(}\frac{\nabla U_{k}}{\sqrt{|\nabla U_{k}|^{2}+\delta^{2}}}\Big{)}-C_{1}\Delta U_{k}+\lambda(f-U_{k})+C_{2}U_{k} (18)

with the boundary condition

∇U⋅n=∇(Δ​U)⋅n=0on∂Ω.\nabla U\cdot n=\nabla(\Delta U)\cdot n=0\quad\text{on}\quad\partial\Omega. (19)
Refer to caption
(a) Damaged Dog image
Refer to caption
(b) α=1\alpha=1 (37.50)
Refer to caption
(c) α=1.2\alpha=1.2 (38.78)
Refer to caption
(d) α=1.6\alpha=1.6 (39.13)
Refer to caption
(e) α=1.8\alpha=1.8 (39.02)
Refer to caption
(f) α\alpha=2 (38.88)
Figure 1: Inpainted results of our model with different fractional power α\alpha, tested on Dog image. In the bracket we have reported the PSNR.
Refer to caption
(a) Damaged Grayshade
Refer to caption
(b) α=1\alpha=1 (35.11)
Refer to caption
(c) α=1.2\alpha=1.2 (42.39)
Refer to caption
(d) α=1.6\alpha=1.6 (45.39)
Refer to caption
(e) α=1.8\alpha=1.8 (45.40)
Refer to caption
(f) α\alpha=2 (45.38)
Figure 2: Inpainted results of our model with different fractional power of Laplacian α\alpha, tested on Grayshade image. In the bracket we have reported the PSNR.
Refer to caption
(a) Damaged Grayshade
Refer to caption
(b) Damaged Kaleidoscope
Refer to caption
(c) T​V−L2TV-L^{2}(41.64)
Refer to caption
(d) T​V−H−1TV-H^{-1}(45.23)
Refer to caption
(e) Our (45.40)
Refer to caption
(f) T​V−L2TV-L^{2}(26.24)
Refer to caption
(g) T​V−H−1TV-H^{-1}(26.90)
Refer to caption
(h) Our (27.89)
Figure 3: Comparison of our inpainted results with the results of T​VL2TV_{L}^{2} and T​V−H−1TV-H^{-1} model on Grayshade and Kaleidoscope image. In the bracket we have reported the PSNR.
Refer to caption
(a) Damaged Elephant 1
Refer to caption
(b) Damaged Elephant 2
Refer to caption
(c) Damaged Dog image
Refer to caption
(d) T​V−L2TV-L^{2}(35.02)
Refer to caption
(e) T​V−H−1TV-H^{-1}(33.93)
Refer to caption
(f) Our (36.16)
Refer to caption
(g) T​V−L2TV-L^{2}(27.53)
Refer to caption
(h) T​V−H−1TV-H^{-1}(28.29)
Refer to caption
(i) Our (28.51)
Refer to caption
(j) T​V−L2TV-L^{2}(38.68)
Refer to caption
(k) T​V−H−1TV-H^{-1}(38.13)
Refer to caption
(l) Our (39.14)
Figure 4: Comparison of our inpainted results with the results of T​VL2TV_{L}^{2} and T​V−H−1TV-H^{-1} model on elephant and dog image. In the bracket we have reported the PSNR.
Refer to caption
(a) Damaged Lena 1
Refer to caption
(b) Damaged Lena 2
Refer to caption
(c) Damaged Lena 3
Refer to caption
(d) T​V−L2TV-L^{2}(31.58)
Refer to caption
(e) CVMS (31.99)
Refer to caption
(f) Our (32.92)
Refer to caption
(g) T​V−L2TV-L^{2}(34.86)
Refer to caption
(h) CVMS (34.97)
Refer to caption
(i) Our (35.62)
Refer to caption
(j) T​V−L2TV-L^{2}(26.08)
Refer to caption
(k) CVMS (27.67)
Refer to caption
(l) Our (28.34)
Figure 5: Comparison of our inpainted results with the results of T​VL2TV_{L}^{2} and T​V−H−1TV-H^{-1} model on Lena image. In the bracket we have reported the PSNR.
Refer to caption
(a) Damaged Barbara 1
Refer to caption
(b) Damaged Barbara 2
Refer to caption
(c) Damaged Barbara 3
Refer to caption
(d) T​V−L2TV-L^{2} (26.93)
Refer to caption
(e) CVMS (27.70)
Refer to caption
(f) Our (28.33)
Refer to caption
(g) T​V−L2TV-L^{2}(29.87)
Refer to caption
(h) CVMS (30.56)
Refer to caption
(i) Our (31.40)
Refer to caption
(j) T​V−L2TV-L^{2}(23.10)
Refer to caption
(k) CVMS (25.77)
Refer to caption
(l) Our (26.42)
Figure 6: Comparison of inpainted results of our model with the results of T​VL2TV_{L}^{2} and T​V−H−1TV-H^{-1} model on Barbara image. In the bracket we have reported the PSNR.

Now, we replace the Laplacian Δ\Delta in (3) by it’s fractional counterpart −(−Δ)α2-(-\Delta)^{\frac{\alpha}{2}} in the spirit of fCH and obtained a fractional model. In general, fractional models are more effective and flexible for image processing problems fCH . So the equation (3) reduces to:

Uk+1−UkΔ​t\displaystyle\frac{U_{k+1}-U_{k}}{\Delta t} +μ​(−Δ)α​Uk+1+C1​(−Δ)α2​Uk+1+C2​Uk+1\displaystyle+\mu(-\Delta)^{\alpha}U_{k+1}+C_{1}(-\Delta)^{\frac{\alpha}{2}}U_{k+1}+C_{2}U_{k+1}
=∇⋅(∇Uk|∇Uk|2+δ2)+C1​(−Δ)α2​Uk+λ​(f−Uk)+C2​Uk\displaystyle=\nabla\cdot\Big{(}\frac{\nabla U_{k}}{\sqrt{|\nabla U_{k}|^{2}+\delta^{2}}}\Big{)}+C_{1}(-\Delta)^{\frac{\alpha}{2}}U_{k}+\lambda(f-U_{k})+C_{2}U_{k} (20)

with the boundary condition

∇U⋅n=∇(Δ​U)⋅n=0on∂Ω.\nabla U\cdot n=\nabla(\Delta U)\cdot n=0\quad\text{on}\quad\partial\Omega. (21)

Here α\alpha is fractional power for 0<α≤2.0<\alpha\leq 2. We know that the Laplacian −Δ-\Delta has a complete orthonormal eigenvectors {φm,n}\{\varphi_{m,n}\} satisfying the Neumann boundary condition corresponding to eigenvalues λm,n\lambda_{m,n}. Consider the eigenvalue problem in Ω=[0,a]×[0,b]\Omega=[0,a]\times[0,b]

−Δ​φm,n=λm,n​φm,n\displaystyle-\Delta\varphi_{m,n}=\lambda_{m,n}\varphi_{m,n}
∇φm,n⋅n=0on∂Ω,\displaystyle\nabla\varphi_{m,n}\cdot n=0\quad\text{on}\quad\partial\Omega,

for m,n=1,2,…m,n=1,2,\dots where

λm,n=π2​((m−1)2a2+(n−1)2b2)\lambda_{m,n}=\pi^{2}\Big{(}\frac{(m-1)^{2}}{a^{2}}+\frac{(n-1)^{2}}{b^{2}}\Big{)}
φm,n=1a​b​c​o​s​((m−1)​π​xa)​c​o​s​((n−1)​π​ya).\varphi_{m,n}=\frac{1}{\sqrt{ab}}cos\Big{(}\frac{(m-1)\pi x}{a}\Big{)}cos\Big{(}\frac{(n-1)\pi y}{a}\Big{)}.

Let us define the space fCH ,

𝒰α:={u=∑m,n=1∞u^m,n​φm,n;u^m,n=⟨u,φm,n⟩L2:∑m,n=1∞|u^m,n|2​|λm,n|α2<∞,0<α≤2}\mathscr{U}_{\alpha}:=\Big{\{}u=\sum\limits_{m,n=1}^{\infty}\hat{u}_{m,n}\varphi_{m,n};\hat{u}_{m,n}=\langle u,\varphi_{m,n}\rangle_{L^{2}}:\\ \sum\limits_{m,n=1}^{\infty}|\hat{u}_{m,n}|^{2}|\lambda_{m,n}|^{\frac{\alpha}{2}}<\infty,0<\alpha\leq 2\Big{\}} (22)

Then for any u∈𝒰αu\in\mathscr{U}_{\alpha}, the fractional Laplace operator can be defined by

(−Δ)α2​u=∑m,n=1∞u^m,n​λm,nα2​φm,n.(-\Delta)^{\frac{\alpha}{2}}u=\sum\limits_{m,n=1}^{\infty}\hat{u}_{m,n}\lambda_{m,n}^{\frac{\alpha}{2}}\varphi_{m,n}. (23)

Hence (−Δ)α2(-\Delta)^{\frac{\alpha}{2}} has the same interpretation as −Δ-\Delta in terms of its spectral decomposition. Fourier spectral methods represent the truncated series expansion when a finite number of orthonormal eigenfunction {φm,n}\{\varphi_{m,n}\} is considered. Consistency, stabilit and convergence of the time discretized model for the case α=2\alpha=2 has been carried out in the Appendix B and the analysis for the fractional one is a matter of future research as it need fractional Green’s type theorem to prove the results. Taking Fourier transform of (3) and using the spectral decomposition of Laplacian (23) we get

U^k+1​(i,j)=1Δ​t+(C1​λi,jα2+C2)​U^k​(i,j)+λ​(f−Uk)^​(i,j)+κ^​(i,j)1Δ​t+μ​λi,jα+C1​λi,jα2+C2\widehat{U}_{k+1}(i,j)=\frac{\frac{1}{\Delta t}+(C_{1}\lambda_{i,j}^{\frac{\alpha}{2}}+C_{2})\hat{U}_{k}(i,j)+\widehat{\lambda(f-U_{k})}(i,j)+\hat{\kappa}(i,j)}{\frac{1}{\Delta t}+\mu\lambda_{i,j}^{\alpha}+C_{1}\lambda_{i,j}^{\frac{\alpha}{2}}+C_{2}} (24)

where κ=∇⋅(∇Uk|∇Uk|2+δ2)\kappa=\nabla\cdot\Big{(}\frac{\nabla U_{k}}{\sqrt{|\nabla U_{k}|^{2}+\delta^{2}}}\Big{)}. In every time step we will calculate the Uk^\widehat{U_{k}} and the real part of the inverse Fourier transform will be the solution UkU_{k}.

4 Numerical results and discussion

In this section, we present some numerical results on few standard test images and compared them with the results of T​V−H−1TV-H^{-1} and T​V−L2TV-L^{2} model. For our model, we set the parameters as Δ​t=1\Delta t=1 and δ=0.01,μ=0.9\delta=0.01,\mu=0.9, the fidelity parameter λ0=250\lambda_{0}=250 and the constant occurred in convexity splitting is set as C1=1δ,C2=50.C_{1}=\frac{1}{\delta},C_{2}=50. The fractional parameter α\alpha is different for different images which is mentioned in the corresponding table. For comparison we have used the code of T​V−L2TV-L^{2} and T​V−H−1TV-H^{-1} provided in codeinp . To compare the quality of the results we calculate three quality metrics namely peak signal to noise ratio (PSNR), siganl to noise ratio (SNR) and structural similarity index measure (SSIM) ssim . PSNR is the ratio between the maximum possible value of the original image and the mean squared error between the original and the resulting image in log scale and SNR is the ratio between the variance of the original image and the variance of the noise. SSIM represents the geometric similarity of the resulting image with the original image. For all the three quality metrics, higher the value implies better the result.

To investigate the performance with different fractional power α\alpha, we have tested our model with different values of α\alpha on few images. In Fig. 1, we have shown the inpainting results of dog image withe different value of α\alpha namely 1,1.2,1.6,1.8,21,1.2,1.6,1.8,2. The results corresponding to α=1,2\alpha=1,2 are not good enough. We have calculated the quality metrics PSNR, SNR, and SSIM for all the results and reported them in Table 1. From the table one can see that for α=1.4\alpha=1.4 and 1.61.6 the results are better than the others specially better than the results of the cases when α\alpha is an integer. Thus the model with fractional power is better than the the usual model.

Similarly, we have tested our model on grayshade image for different values of α\alpha. The inpainted results are shown in Fig. 2. To compare the results we have calculated PSNR, SNR and SSIM for all the cases and reported in Table 1. From the table one can see that in this case α=1.8\alpha=1.8 gives the best result in this case.

Now we will compare results of our model with that of T​V−H−1TV-H^{-1} and T​V−L2TV-L^{2} model. In Fig. 3, we have taken two images namely Grayshade image and Kaleidoscope image and corresponding damaged image has been shown in the first row. Inpainted results of all the three models for the grayshade image are presented in second row and results of Kaleidoscope image in third row. From the figure one can see that the results of T​V−L2TV-L^{2} is not good and results of T​V−H−1TV-H^{-1} and our model is better. For both the images the fractional power α\alpha is chosen as 1.8. PSNR, SNR and SSIM for these images are reported in Table 2. From the table one can see that for Grayshade image our model gives little better result. For the Kaleidoscope image our model gives much better result than the results of the other two model. For example, SSIM for T​V−H−1TV-H^{-1} model is 0.87 whereas SSIM for our model is 0.91.

In Fig. 4, we have presented the results of two images namely Elephant and dog. We have taken two different damaged region for the elephant image and one damage region for the dog image and all the damaged images are shown in the first row of Fig. 4. Inpainted results are shown in row two, three and four. The first column contains results of T​V−L2TV-L^{2} model, second column contains the result of T​V−H−1TV-H^{-1} and results of our model is displayed in third column. For these image the fractional parameter α\alpha of our model is chosen as 1.4,1.8 and 1.6 respectively. From the table one can see that for all the case our model perform better than other two models in terms of PSNR, SNR and SSIM. For the Dog and first elephant image there is significant improvement in the performance in our model.

In the next two experiments, we have compared our results with the results obtained by T​V−L2TV-L^{2} and CVMS twophase model. For the CVMS model we have discretized the equation (11) by convexity splitting in time and Fourier spectral method in space which is presented in Appendix A.

First experiment is on Lena image of size 512×512512\times 512. We select three different inpainting domain for Lena image and the corresponding damaged images are shown in the first row of Fig. 5. The results of first damaged image have been shown in second row and the results of second damaged image have been shown in third row and results of third damaged image is shown in fourth row. From the fourth row of Fig. 5 one can see that our model is giving better result than the other two. For all the models we have reported image quality metrics PSNR, SNR, and SSIM in Table 3. Also the best fractional power α\alpha is mentioned for each image in the same table. From the Table 3, it is evident that our model perform better than the other two models in terms of PSNR, SNR, and SSIM.

Finally, we did the experiment on the Barbara image of size 512×512512\times 512. Three different inpainting domain are chosen and the resulting degraded images are shown in the first row of Fig. 6. The results of T​V−L2TV-L^{2} model have been shown in first column, second column contains the result of CVMS model twophase and our results are in third column. In this case, also we can observe from the Table 3 our model perform much better than the other two in terms of PSNR, SNR and SSIM. Note that the Barbara image contains some structures and preserving them is a challenge and in this case also our model gives better result.

5 Conclusion

Here we have presented a new 4th order PDE model for grayscale image inpainting using a variant of Mumford-Shah model. Convexity splitting in time has been used for time discretization and further the time discretized model has been modified by replacing the Laplace term by its fractional version. Fourier spectral method in space has been to get the complete discretization. Consistency, stability and convergence have been established for the time discretized scheme. Numerical results show the superiority of our model over T​V−L2TV-L^{2} and T​V−H−1TV-H^{-1} and convex variant of Mumford-Shah model.

Appendix A

The equation (11) can be written as gradient descent of two functional E1,E2E_{1},E_{2} in L2L^{2} where

E1​(u)=∫Ω(μ2​|∇u|2+|∇u|2+δ2)​𝑑x,E2​(u)=∫Ωλ2​(f−u)2​𝑑x.E_{1}(u)=\int_{\Omega}\Big{(}\frac{\mu}{2}|\nabla u|^{2}+\sqrt{|\nabla u|^{2}+\delta^{2}}\Big{)}dx,\quad E_{2}(u)=\int_{\Omega}\frac{\lambda}{2}(f-u)^{2}dx. (25)

To apply convexity splitting idea we split E1E_{1} as E11−E12E_{11}-E_{12} and E2E_{2} as E21−E22E_{21}-E_{22} where

E11=∫Ω(μ2​|∇u|2+C12​|∇u|2)​𝑑x and E_{11}=\int_{\Omega}\Big{(}\frac{\mu}{2}|\nabla u|^{2}+\frac{C_{1}}{2}|\nabla u|^{2}\Big{)}dx\quad\text{ and }
E12=∫Ω(C12​|∇u|2−|∇u|2+δ2)​𝑑xE_{12}=\int_{\Omega}\Big{(}\frac{C_{1}}{2}|\nabla u|^{2}-\sqrt{|\nabla u|^{2}+\delta^{2}}\Big{)}dx

and

E21=∫ΩC22​|u|2​𝑑x, and E22=∫Ω(−λ2​(f−u)2+C22​|u|2)​𝑑x,E_{21}=\int_{\Omega}\frac{C_{2}}{2}|u|^{2}dx,\quad\text{ and }\quad E_{22}=\int_{\Omega}\Big{(}-\frac{\lambda}{2}(f-u)^{2}+\frac{C_{2}}{2}|u|^{2}\Big{)}dx,

where C1,C2C_{1},C_{2} are positive and should be large enough so that Ei​jE_{ij} for i,j=1,2i,j=1,2 becomes convex.

So the semi-discrete time-stepping scheme for our model is

Uk+1−UkΔ​t=−∇L2(E11​(Uk+1)−E12​(Uk))−∇L2(E21​(Uk+1)−E22​(Uk)).\frac{U_{k+1}-U_{k}}{\Delta t}=-\nabla_{L^{2}}(E_{11}(U_{k+1})-E_{12}(U_{k}))-\nabla_{L^{2}}(E_{21}(U_{k+1})-E_{22}(U_{k})). (26)
Uk+1−UkΔ​t\displaystyle\frac{U_{k+1}-U_{k}}{\Delta t} =−{−μ​Δ​Uk+1−C1​Δ​Uk+1+C1​Δ​Uk−∇⋅(∇Uk|∇Uk|2+δ2)}\displaystyle=-\{-\mu\Delta U_{k+1}-C_{1}\Delta U_{k+1}+C_{1}\Delta U_{k}-\nabla\cdot\Big{(}\frac{\nabla U_{k}}{\sqrt{|\nabla U_{k}|^{2}+\delta^{2}}}\Big{)}\}
−{C2​Uk+1−C2​Uk−λ​(f−Uk)}\displaystyle\qquad-\{C_{2}U_{k+1}-C_{2}U_{k}-\lambda(f-U_{k})\}

Simplifying we get the time discretized model as:

Uk+1−UkΔ​t\displaystyle\frac{U_{k+1}-U_{k}}{\Delta t} −μ​Δ​Uk+1−C1​Δ​Uk+1+C2​Uk+1\displaystyle-\mu\Delta U_{k+1}-C_{1}\Delta U_{k+1}+C_{2}U_{k+1}
=∇⋅(∇Uk|∇Uk|2+δ2)−C1​Δ​Uk+λ​(f−Uk)+C2​Uk\displaystyle=\nabla\cdot\Big{(}\frac{\nabla U_{k}}{\sqrt{|\nabla U_{k}|^{2}+\delta^{2}}}\Big{)}-C_{1}\Delta U_{k}+\lambda(f-U_{k})+C_{2}U_{k} (27)

To get the complete discretized scheme we will use Fourier spectral method in space. Let Uk^\widehat{U_{k}} be the discrete Fourier trans form of UkU_{k}. Taking Fourier transform of (Appendix A) and using the relation Δ​u^=L​u^\widehat{\Delta u}=L\hat{u} and rearranging we get the final numerical scheme for the first model as:

U^k+1​(i,j)=(1Δ​t−C1​Li,j+C2)​Uk^​(i,j)+κk^​(i,j)+λ​(f−Uk)^​(i,j)1Δ​t−(μ+C1)​Li,j+C2\widehat{U}_{k+1}(i,j)=\frac{(\frac{1}{\Delta t}-C_{1}L_{i,j}+C_{2})\widehat{U_{k}}(i,j)+\widehat{\kappa_{k}}(i,j)+\widehat{\lambda(f-U_{k})}(i,j)}{\frac{1}{\Delta t}-(\mu+C_{1})L_{i,j}+C_{2}} (28)

where κk=∇⋅∇Uk|∇Uk|2+δ2\kappa_{k}=\nabla\cdot\frac{\nabla U_{k}}{\sqrt{|\nabla U_{k}|^{2}+\delta^{2}}} and L=(Li,j)L=(L_{i,j}) is matrix with size m×nm\times n with

Li,j=2(Δ​x)2​(c​o​s​(2​π​im)−1)+2(Δ​y)2​(c​o​s​(2​π​jn)−1).L_{i,j}=\frac{2}{(\Delta x)^{2}}\Big{(}cos(\frac{2\pi i}{m})-1\Big{)}+\frac{2}{(\Delta y)^{2}}\Big{(}cos(\frac{2\pi j}{n})-1\Big{)}.

Appendix B

Now, we will discuss the consistency, stability and convergence of the time discretized model (3) that is for the fractional model with α=2.\alpha=2.

Theorem 5.1 (Consistency, Stability and Convergence)

Let uu be the exact solution of (15) and uk=u​(k​Δ​t)u_{k}=u(k\Delta t) be the exact solution at time k​Δ​tk\Delta t for a time step Δ​t>0\Delta t>0 and k∈ℕ.k\in\mathbb{N}. Let UkU_{k} be the kk-th iterate of(3) with constant C​1>1δ,C2>λ0C1>\frac{1}{\delta},C_{2}>\lambda_{0}. Then the following statements are true:

  1. (a)

    Under the assumption that ‖ut​t‖2\|u_{tt}\|_{2} and ‖ut‖2\|u_{t}\|_{2} are bounded, the numerical scheme is consistent with the continuous equation and of order one in time.

  2. (b)

    the solution UkU_{k} is bounded on a finite time interval [0,T][0,T], for all Δ​t>0\Delta t>0. In particular for k​Δ​t≤Tk\Delta t\leq T, T>0T>0 fixed, we have for every Δ​t>0\Delta t>0

    ‖∇Uk‖2+Δ​t​K1​‖Δ​Uk+1‖2+Δ​t​K2​‖∇Δ​Uk‖2\displaystyle\|\nabla U_{k}\|^{2}+\Delta tK_{1}\|\Delta U_{k+1}\|^{2}+\Delta tK_{2}\|\nabla\Delta U_{k}\|^{2} (29)
    ≤eK​T​(‖∇U0‖2+Δ​t​K1​‖Δ​U0‖2+Δ​t​K2​‖∇Δ​U0‖2+Δ​t​T​C​(Ω,D,λ0,f)),\displaystyle\leq e^{KT}\Big{(}\|\nabla U_{0}\|^{2}+\Delta tK_{1}\|\Delta U_{0}\|^{2}+\Delta tK_{2}\|\nabla\Delta U_{0}\|^{2}+\Delta tTC(\Omega,D,\lambda_{0},f)\Big{)},
  3. (c)

    The discretization error ek=uk−Uke_{k}=u_{k}-U_{k}. For smooth solution uku_{k} and UkU_{k}, the error eke_{k} converges to 0 as Δ​t→0\Delta t\to 0. In particular we have for k​Δ​t≤Tk\Delta t\leq T, T>0T>0 fixed, that

    ‖∇ek‖2+Δ​t​M1​‖Δ​Uk+1‖2+Δ​t​M2​‖∇Δ​ek‖2≤TM3​eM4​T​(Δ​t)2\|\nabla e_{k}\|^{2}+\Delta tM_{1}\|\Delta U_{k+1}\|^{2}+\Delta tM_{2}\|\nabla\Delta e_{k}\|^{2}\leq\frac{T}{M_{3}}e^{M_{4}T}(\Delta t)^{2}

    for suitable constants M​1,M2,M3,M​4.M1,M_{2},M_{3},M4.

Proof

(a) Putting u=uku=u_{k} in equation (15) and Uk=ukU_{k}=u_{k} in (3) and subtracting the later we get the error τk\tau_{k} at kt​hk^{th} time step as

τk=uk+1−ukΔ​t−ut​(k​Δ​t)+μ​Δ2​(uk+1−uk)−C1​Δ​(uk+1−uk)+C2​(uk+1−uk)\tau_{k}=\frac{u_{k+1}-u_{k}}{\Delta t}-u_{t}(k\Delta t)+\mu\Delta^{2}(u_{k+1}-u_{k})-C_{1}\Delta(u_{k+1}-u_{k})+C_{2}(u_{k+1}-u_{k}) (30)

Then applying Taylors theorem we get,

τk=(ut​(k​Δ​t)+Δ​t​ut​t​(ξ1))−ut​(k​Δ​t)+Δ​t​(μ​Δ2​ut​(ξ2)−C1​Δ​ut​(ξ3)+C2​ut​(ξ4))\tau_{k}=(u_{t}(k\Delta t)+\Delta tu_{tt}(\xi_{1}))-u_{t}(k\Delta t)+\Delta t(\mu\Delta^{2}u_{t}(\xi_{2})-C_{1}\Delta u_{t}(\xi_{3})+C_{2}u_{t}(\xi_{4}))

Now using the boundedness assumptions we get ‖τk‖=O​(Δ​t)\|\tau_{k}\|=O(\Delta t). Hence the global truncation error τ=maxk⁡‖τk‖=O​(Δ​t)\tau=\max\limits_{k}\|\tau_{k}\|=O(\Delta t).

Proof

(b) Consider the time discretized scheme

Uk+1−UkΔ​t\displaystyle\frac{U_{k+1}-U_{k}}{\Delta t} +μ​Δ2​Uk+1−C1​Δ​Uk+1+C2​Uk+1\displaystyle+\mu\Delta^{2}U_{k+1}-C_{1}\Delta U_{k+1}+C_{2}U_{k+1}
=∇⋅(∇Uk|∇Uk|2+δ2)−C1​Δ​Uk+λ​(f−Uk)+C2​Uk\displaystyle=\nabla\cdot\Big{(}\frac{\nabla U_{k}}{\sqrt{|\nabla U_{k}|^{2}+\delta^{2}}}\Big{)}-C_{1}\Delta U_{k}+\lambda(f-U_{k})+C_{2}U_{k} (31)

Multiplying the time discretised scheme (Proof) by −Δ​Uk+1-\Delta U_{k+1} and integrating over Ω\Omega we get

1Δ​t(∥∇Uk+1∥2−\displaystyle\frac{1}{\Delta t}\Big{(}\|\nabla U_{k+1}\|^{2}- ⟨∇Uk,Uk+1⟩)+μ∥∇ΔUk+1∥2+C1∥ΔUk+1∥2+C2∥∇Uk+1∥2\displaystyle\langle\nabla U_{k},U_{k+1}\rangle\Big{)}+\mu\|\nabla\Delta U_{k+1}\|^{2}+C_{1}\|\Delta U_{k+1}\|^{2}+C_{2}\|\nabla U_{k+1}\|^{2}
=⟨∇∇⋅∇Uk|∇Uk|2+δ2,∇Uk+1⟩+C1​⟨Δ​Uk,Δ​Uk+1⟩\displaystyle=\Big{\langle}\nabla\nabla\cdot\frac{\nabla U_{k}}{\sqrt{|\nabla U_{k}|^{2}+\delta^{2}}},\nabla U_{k+1}\Big{\rangle}+C_{1}\langle\Delta U_{k},\Delta U_{k+1}\rangle
+⟨∇λ​(f−Uk),∇Uk+1⟩+C2​⟨∇Uk,∇Uk+1⟩\displaystyle\qquad+\langle\nabla\lambda(f-U_{k}),\nabla U_{k+1}\rangle+C_{2}\langle\nabla U_{k},\nabla U_{k+1}\rangle

Using Young’s inequality on the product terms we get,

12​Δ​t​(‖∇Uk+1‖2−‖∇Uk‖2)+μ​‖∇Δ​Uk+1‖2+C1​‖Δ​Uk+1‖2+C2​‖∇Uk+1‖2\displaystyle\frac{1}{2\Delta t}\Big{(}\|\nabla U_{k+1}\|^{2}-\|\nabla U_{k}\|^{2}\Big{)}+\mu\|\nabla\Delta U_{k+1}\|^{2}+C_{1}\|\Delta U_{k+1}\|^{2}+C_{2}\|\nabla U_{k+1}\|^{2}
≤δ1​‖∇Uk+1‖2+14​δ1​‖∇∇⋅∇Uk|∇Uk|2+δ2‖2+C1​δ2​‖Δ​Uk+1‖2+C14​δ2​‖Δ​Uk‖2\displaystyle\quad\leq\delta_{1}\|\nabla U_{k+1}\|^{2}+\frac{1}{4\delta_{1}}\|\nabla\nabla\cdot\frac{\nabla U_{k}}{\sqrt{|\nabla U_{k}|^{2}+\delta^{2}}}\|^{2}+C_{1}\delta_{2}\|\Delta U_{k+1}\|^{2}+\frac{C_{1}}{4\delta_{2}}\|\Delta U_{k}\|^{2}
+δ3​‖∇Uk+1‖2+14​δ3​‖λ​(f−Uk)‖2+C2​δ4​‖∇Uk+1‖2+C24​δ4​‖∇Uk‖2\displaystyle\quad+\delta_{3}\|\nabla U_{k+1}\|^{2}+\frac{1}{4\delta_{3}}\|\lambda(f-U_{k})\|^{2}+C_{2}\delta_{4}\|\nabla U_{k+1}\|^{2}+\frac{C_{2}}{4\delta_{4}}\|\nabla U_{k}\|^{2} (32)

Using the following two estimates nondenoising in (Proof)

‖λ​(f−Uk)‖2≤2​λ02​‖∇Uk‖2+C​(Ω,D,λ0,f)\|\lambda(f-U_{k})\|^{2}\leq 2\lambda_{0}^{2}\|\nabla U_{k}\|^{2}+C(\Omega,D,\lambda_{0},f) (33)

and

‖∇∇⋅∇Uk|∇Uk|2+δ2‖2≤C¯​(δ,Ω)​(‖∇Uk‖2+‖∇Δ​Uk‖2)\|\nabla\nabla\cdot\frac{\nabla U_{k}}{\sqrt{|\nabla U_{k}|^{2}+\delta^{2}}}\|^{2}\leq\bar{C}(\delta,\Omega)(\|\nabla U_{k}\|^{2}+\|\nabla\Delta U_{k}\|^{2}) (34)

we get,

(12​Δ​t+C2​(1−δ4)−(δ1+δ3))​‖∇Uk+1‖2+C1​(1−δ2)​‖Δ​Uk+1‖2+μ​‖∇Δ​Uk+1‖2\displaystyle\Big{(}\frac{1}{2\Delta t}+C_{2}(1-\delta_{4})-(\delta_{1}+\delta_{3})\Big{)}\|\nabla U_{k+1}\|^{2}+C_{1}(1-\delta_{2})\|\Delta U_{k+1}\|^{2}+\mu\|\nabla\Delta U_{k+1}\|^{2}
≤(12​Δ​t+C¯4​δ1+C24​δ4+2​λ024​δ3)​‖∇Uk‖2+C14​δ2​‖Δ​Uk‖2+C¯4​δ1​‖∇Δ​Uk+1‖2+C​(Ω,D,λ0,f)\displaystyle\quad\leq\Big{(}\frac{1}{2\Delta t}+\frac{\bar{C}}{4\delta_{1}}+\frac{C_{2}}{4\delta_{4}}+\frac{2\lambda_{0}^{2}}{4\delta_{3}}\Big{)}\|\nabla U_{k}\|^{2}+\frac{C_{1}}{4\delta_{2}}\|\Delta U_{k}\|^{2}+\frac{\bar{C}}{4\delta_{1}}\|\nabla\Delta U_{k+1}\|^{2}+C(\Omega,D,\lambda_{0},f)

Let us choose δ1=δ3=14\delta_{1}=\delta_{3}=\frac{1}{4} and δ2=δ4=12\delta_{2}=\delta_{4}=\frac{1}{2}, we get,

(12​Δ​t+C1−12)​‖∇Uk+1‖2+C12​‖Δ​Uk+1‖2+μ​‖∇Δ​Uk+1‖2\displaystyle\Big{(}\frac{1}{2\Delta t}+\frac{C_{1}-1}{2}\Big{)}\|\nabla U_{k+1}\|^{2}+\frac{C_{1}}{2}\|\Delta U_{k+1}\|^{2}+\mu\|\nabla\Delta U_{k+1}\|^{2}
≤(12​Δ​t+C¯+C22+2​λ02)​‖∇Uk‖2+C12​‖Δ​Uk‖2+C¯​‖∇Δ​Uk+1‖2+C​(Ω,D,λ0,f)\displaystyle\quad\leq\Big{(}\frac{1}{2\Delta t}+\bar{C}+\frac{C_{2}}{2}+2\lambda_{0}^{2}\Big{)}\|\nabla U_{k}\|^{2}+\frac{C_{1}}{2}\|\Delta U_{k}\|^{2}+\bar{C}\|\nabla\Delta U_{k+1}\|^{2}+C(\Omega,D,\lambda_{0},f)

Since C2>1C_{2}>1, all the coefficients are positive. Now multiplying both sides by 2​Δ​t2\Delta t we get,

A1​‖∇Uk+1‖2+Δ​t​C1​‖Δ​Uk+1‖2+2​Δ​t​μ​‖∇Δ​Uk+1‖2\displaystyle A_{1}\|\nabla U_{k+1}\|^{2}+\Delta tC_{1}\|\Delta U_{k+1}\|^{2}+2\Delta t\mu\|\nabla\Delta U_{k+1}\|^{2}
≤A2​‖∇Uk‖2+Δ​t​C1​‖Δ​Uk‖2+2​Δ​t​C¯​‖∇Δ​Uk‖2+2​Δ​t​C​(Ω,D,λ0,f),\displaystyle\quad\leq A_{2}\|\nabla U_{k}\|^{2}+\Delta tC_{1}\|\Delta U_{k}\|^{2}+2\Delta t\bar{C}\|\nabla\Delta U_{k}\|^{2}+2\Delta tC(\Omega,D,\lambda_{0},f), (35)

where A1=1+Δ​t​(C2−1),A2=1+Δ​t​(C2+2​C¯+4​λ02)A_{1}=1+\Delta t(C_{2}-1),A_{2}=1+\Delta t(C_{2}+2\bar{C}+4\lambda_{0}^{2}). Dividing both sides by A1A_{1} we get,

‖∇Uk+1‖2+Δ​t​C1A1​‖Δ​Uk+1‖2+Δ​t​2​μA1​‖∇Δ​Uk+1‖2\displaystyle\|\nabla U_{k+1}\|^{2}+\Delta t\frac{C_{1}}{A_{1}}\|\Delta U_{k+1}\|^{2}+\Delta t\frac{2\mu}{A_{1}}\|\nabla\Delta U_{k+1}\|^{2}
≤A2A1​(‖∇Uk‖2+Δ​t​C1A2​‖Δ​Uk‖2+Δ​t​2​C¯A2​‖∇Δ​Uk‖2)+2​Δ​tA1​C​(Ω,D,λ0,f),\displaystyle\quad\leq\frac{A_{2}}{A_{1}}\Big{(}\|\nabla U_{k}\|^{2}+\Delta t\frac{C_{1}}{A_{2}}\|\Delta U_{k}\|^{2}+\Delta t\frac{2\bar{C}}{A_{2}}\|\nabla\Delta U_{k}\|^{2}\Big{)}+\frac{2\Delta t}{A_{1}}C(\Omega,D,\lambda_{0},f), (36)

Choose A1,A2A_{1},A_{2} in such a way that A2A1>1\frac{A_{2}}{A_{1}}>1 and μ​A2C¯​A1>1\frac{\mu A_{2}}{\bar{C}A_{1}}>1 (this is possible if we choose big enough λ0\lambda_{0}), then above equation reduces to

‖∇Uk+1‖2+Δ​t​C1A1​‖Δ​Uk+1‖2+Δ​t​2​μA1​‖∇Δ​Uk+1‖2\displaystyle\|\nabla U_{k+1}\|^{2}+\Delta t\frac{C_{1}}{A_{1}}\|\Delta U_{k+1}\|^{2}+\Delta t\frac{2\mu}{A_{1}}\|\nabla\Delta U_{k+1}\|^{2}
≤A2A1​(‖∇Uk‖2+Δ​t​C1A1​‖Δ​Uk‖2+Δ​t​2​μA1​‖∇Δ​Uk‖2)+Δ​t​C​(Ω,D,λ0,f).\displaystyle\quad\leq\frac{A_{2}}{A_{1}}\Big{(}\|\nabla U_{k}\|^{2}+\Delta t\frac{C_{1}}{A_{1}}\|\Delta U_{k}\|^{2}+\Delta t\frac{2\mu}{A_{1}}\|\nabla\Delta U_{k}\|^{2}\Big{)}+\Delta tC(\Omega,D,\lambda_{0},f). (37)

So by induction we will get,

‖∇Uk+1‖2+Δ​t​C1A1​‖Δ​Uk+1‖2+Δ​t​2​μA1​‖∇Δ​Uk+1‖2\displaystyle\|\nabla U_{k+1}\|^{2}+\Delta t\frac{C_{1}}{A_{1}}\|\Delta U_{k+1}\|^{2}+\Delta t\frac{2\mu}{A_{1}}\|\nabla\Delta U_{k+1}\|^{2}
≤(A2A1)k​(‖∇U0‖2+Δ​t​C1A1​‖Δ​U0‖2+Δ​t​2​μA1​‖∇Δ​U0‖2)+Δ​t​∑i=0k−1(A2A1)i​C​(Ω,D,λ0,f).\displaystyle\quad\leq\Big{(}\frac{A_{2}}{A_{1}}\Big{)}^{k}\Big{(}\|\nabla U_{0}\|^{2}+\Delta t\frac{C_{1}}{A_{1}}\|\Delta U_{0}\|^{2}+\Delta t\frac{2\mu}{A_{1}}\|\nabla\Delta U_{0}\|^{2}\Big{)}+\Delta t\sum\limits_{i=0}^{k-1}\Big{(}\frac{A_{2}}{A_{1}}\Big{)}^{i}C(\Omega,D,\lambda_{0},f). (38)

Expressing A2A1\frac{A_{2}}{A_{1}} as =1+K3​Δ​t1+K_{3}\Delta t and using the inequality (1+n​x)k≤ek​x(1+nx)^{k}\leq e^{kx}, we get for k​Δ​t<Tk\Delta t<T,

‖∇Uk+1‖2+Δ​t​C1A1​‖Δ​Uk+1‖2+Δ​t​2​μA1​‖∇Δ​Uk+1‖2\displaystyle\|\nabla U_{k+1}\|^{2}+\Delta t\frac{C_{1}}{A_{1}}\|\Delta U_{k+1}\|^{2}+\Delta t\frac{2\mu}{A_{1}}\|\nabla\Delta U_{k+1}\|^{2}
≤eK3​T​(‖∇U0‖2+Δ​t​C1A1​‖Δ​U0‖2+Δ​t​2​μA1​‖∇Δ​U0‖2+T​C​(Ω,D,λ0,f)).\displaystyle\quad\leq e^{K_{3}T}\Big{(}\|\nabla U_{0}\|^{2}+\Delta t\frac{C_{1}}{A_{1}}\|\Delta U_{0}\|^{2}+\Delta t\frac{2\mu}{A_{1}}\|\nabla\Delta U_{0}\|^{2}+TC(\Omega,D,\lambda_{0},f)\Big{)}. (39)
Proof

(c) Putting the value of ut​(k​Δ​t)u_{t}(k\Delta t) from (15) in equation (30) and simplifying we will have,

τk=uk+1−ukΔ​t\displaystyle\tau_{k}=\frac{u_{k+1}-u_{k}}{\Delta t} −∇⋅∇uk|∇uk|2+δ2−λ​(f−uk)+μ​Δ2​uk+1\displaystyle-\nabla\cdot\frac{\nabla u_{k}}{\sqrt{|\nabla u_{k}|^{2}+\delta^{2}}}-\lambda(f-u_{k})+\mu\Delta^{2}u_{k+1}
−C1​Δ​(uk+1−uk)+C2​(uk+1−uk)\displaystyle-C_{1}\Delta(u_{k+1}-u_{k})+C_{2}(u_{k+1}-u_{k}) (40)

Now subtracting (Proof) from (Proof) and putting ek=uk−Uke_{k}=u_{k}-U_{k} we get,

ek+1−ekΔ​t+μ​Δ2​ek+1−C1​Δ​ek+1+C2​ek+1\displaystyle\frac{e_{k+1}-e_{k}}{\Delta t}+\mu\Delta^{2}e_{k+1}-C_{1}\Delta e_{k+1}+C_{2}e_{k+1}
=(∇⋅∇uk|∇uk|2+δ2−∇⋅∇Uk|∇Uk|2+δ2)−C1​Δ​ek+C2​ek−λ​ek+τk\displaystyle=\Big{(}\nabla\cdot\frac{\nabla u_{k}}{\sqrt{|\nabla u_{k}|^{2}+\delta^{2}}}-\nabla\cdot\frac{\nabla U_{k}}{\sqrt{|\nabla U_{k}|^{2}+\delta^{2}}}\Big{)}-C_{1}\Delta e_{k}+C_{2}e_{k}-\lambda e_{k}+\tau_{k} (41)

Multiplying equation (Proof) by ek+1e_{k+1} and integrating over Ω\Omega and arguing as proof (b) and following tvh we will get the required result.

References

  • (1) M. Burger, L. He and C. Schonlieb, Cahn-Hilliard inpainting and a generalization for grayvalue images, SIAM J. Imaging Sci., 2 (4), 1129-1167, 2009.
  • (2) C.B. Schonlieb, and A. Bertozzi, Unconditionally stable schemes for higher order inpainting.Commun. Math. Sci. 9 (2) (2011) 413-457.
  • (3) X. Cai, R. Chan, and T. Zeng, A Two-Stage Image Segmentation Method Using a Convex Variant of the Mumford–Shah Model and Thresholding, SIAM Journal of Imaging Sciences, vol. 6(1), pp. 368-390, 2013.
  • (4) X. Cai, R. Chan, M. Nikolova, and T. Zeng, A Three-Stage Approach for Segmenting Degraded Color Images: Smoothing, Lifting and Thresholding (SLaT). J Sci Comput 72, 1313–1332 (2017).
  • (5) T.F. Chan and J. Shen, Mathematical models for local non-texture inpaintings, SIAM J. Appl. Math., 62(3), 1019-1043, 2001.
  • (6) T.F. Chan and J. Shen, Non-texture inpainting by curvature driven diffusions (CDD), J. Visual Commun. Image Rep., 12(4), 436-449, 2001.
  • (7) T.F. Chan, S.H. Kang and J. Shen, Euler’s elastica and curvature-based inpainting, SIAM J. Appl. Math., 63(2), 564–592, 2002.
  • (8) L. I. Rudin, S. Osher and E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D, 60, 259-268, 1992.
  • (9) S. Esedoglu and J.H. Shen, Digital inpainting based on the Mumford-Shah-Euler image model, Eur. J. Appl. Math., 13(4), 353-370, 2002.
  • (10) L. Bar, T. F. Chan, G. Chung, M. Jung, N. Kiryati, N. Sochen, and L. A. Vese, Mumford and Shah Model and its Applications to Image Segmentation and Image Restoration. In: Handbook of Mathematical Methods in Imaging. Springer, New York, NY, 2011.
  • (11) A. Bertozzi, S. Esedoglu and A. Gillette, Inpainting of binary images using the Cahn-Hilliard equation, IEEE Trans. Image Proc., 16(1), 285-291, 2007.
  • (12) A. Bertozzi, S. Esedoglu and A. Gillette, Analysis of a two-scale Cahn-Hilliard model for image inpainting, Multiscale Modeling and Simulation, 6(3), 913-936, 2007.
  • (13) T.F. Chan, J.H. Shen and H.M. Zhou, Total variation wavelet inpainting, J. Math. Imaging Vis., 25(1), 107-125, 2006.
  • (14) L. Cherfils, H. Fakih and A. Miranville, A complex version of the Cahn-Hilliard equation for grayscale image inpainting, Multiscale Model. Simul. 15 (2017), 575-605.
  • (15) J. Bosch and M. Stoll, A Fractional Inpainting Model Based on the Vector-Valued Cahn-Hilliard Equation, SIAM Journal on Imaging Sciences, vol. 8(4), pp. 2352-2382, 2015.
  • (16) Anis Theljani, Zakaria Belhachmi, Moez Kallel, Maher Moakher. Multiscale Fourth- Order Model for Image Inpainting and Low-Dimensional Sets Recovery. 2016.
  • (17) B.V.R. Kumar, A. Halim, A Linear Fourth Order PDE Based Grayscale Image Inpainting Model, Computational and Applied Mathematics, vol. 38(1), pp. 1-21,2019.
  • (18) A. Halim, B.V.R. Kumar, An anisotropic PDE model for image inpainting, Computers and Mathematics with Applications, vol.79(9), pp. 2701–2721, 2020.
  • (19) Z. Wang, A.C. Bovik, H.R. Sheikh, and E.P. Simoncelli, Image quality assessment: from error visibility to structural similarity, IEEE Transactions on Image Processing, vol. 13 (4), pp. 600–612, 2004.
  • (20) A. Halim, B.V. R. Kumar, A T​V−L2−H−1TV-L^{2}-H^{-1} PDE model for effective denoising, Computers and Mathematics with Applications, vol. 80(9) , pp. 2176–2193, 2020.
  • (21) C.B. Schnlieb, Higher-order total variation inpainting - File Exchange - MATLAB Central.