A NALYSIS AND S IMULATION OF A C OUPLED D IFFUSION BASED I MAGE D ENOISING M ODEL arXiv:1907.04526v2 [math.NA] 7 Aug 2019 A P REPRINT Subit K. Jain School of Basic Sciences Indian Institute of Technology Mandi PIN 175005, INDIA jain.subit@gmail.com Sudeb Majee School of Basic Sciences Indian Institute of Technology Mandi PIN 175005, INDIA sudebmajee@gmail.com Rajendra K. Ray School of Basic Sciences Indian Institute of Technology Mandi PIN 175005, INDIA rajendra@iitmandi.ac.in Ananta K. Majee Department of Mathematics Indian Institute of Technology Delhi PIN 110016, INDIA majee@maths.iitd.ac.in August 8, 2019 A BSTRACT In this study, a new coupled Partial Differential Equation (CPDE) based image denoising model incorporating space-time regularization into non-linear diffusion is proposed. This proposed model is fitted with additive Gaussian noise which performs efficient image smoothing along with the preservation of edges and fine structures. For this purpose, we propose a new functional minimization framework to remove the image noise, which results in solving a system of three partial differential equations (PDEs). Our proposed model is dissimilar from the existing CPDE models as it includes two additional evolution equations to handle edge strength function and data fidelity term. These two evolution equations control the smoothing process and force the resultant denoised solution to be close to the initial solution. To the best of our knowledge, the proposed model is the only work, which deciphers the combined effect of both the terms using separate PDEs. Furthermore, we establish the existence and uniqueness of a weak solution of the proposed system using the time discretization method with H 1 initial data. Finally, we used a generalized weighted average finite difference scheme to efficiently solve the coupled system and experiment results show the effectiveness of the proposed CPDE model. Keywords Image denoising · Space-Time regularization · Non-linear diffusion · Semi discretization · Weak solution 1 Introduction Reduction of generated noise in the restored images without loss of spatial resolution is highly desirable in many image processing applications. The behavior of noise reduction methods is entirely different for large areas with welldefined borders as well as for poorly defined borders. It has been observed that the goal of noise reduction with feature preservation and applicability to the multiple acquisition systems is difficult to achieve through a single method [1, 2]. Therefore, the accurate reconstruction of key features in an image from the noisy measurement is the main objective of our study. The noise caused by image acquisition devices is often modeled by Gaussian random distribution. Generally, the degraded image can be modeled as J = I + η, (1) A PREPRINT - AUGUST 8, 2019 where J is the noisy image, degraded by Gaussian noise with zero mean, defined on an image space Ω, and I depicts the clean image, and η is the additive noise function. Our main task is to obtain the reconstructed image from the observed noisy image. In recent years, PDE based methods have become widely applicable for noise removal and signal reconstruction, due to their well-established mathematical properties and their general wellposedness [1, 3, 4, 5, 6, 7, 8]. Obtaining the steady-state solution is the basic objective of PDE based image restoration problem. This goal can be achieved through two approaches. In the first approach, we directly find a steady state solution of the diffusion model. Whereas in the second approach, the steady-state solution is obtained from the evolution equation of the Euler-Lagrange equation associated with the energy functional. Several techniques are employed to handle the challenges posed by restoration problems. The availability of vast literature does not allow us to review all of them. However, we briefly report some of the robust techniques such that anisotropic diffusion [3, 6, 7, 9, 10, 11], fourth order diffusion [12, 13, 14], complex diffusion [15, 16] and variational based methods [17, 18, 19, 20, 21, 22]. In all of the above-mentioned methods, spatial regularization is used for diffusion function, which is not able to inject the past information into the iterative process of diffusion. Therefore, the scheme may not be effective for solving image denoising problems which contain huge textures in their image domain. To handle this, applications of time-delay regularization with spatial regularization are also available in the literature [23, 24, 25, 26, 27, 28]. The space-time regularization is more suitable approach than Gaussian filtering as it enables to incorporate the information obtained along with the scales into the diffusion process. This type of regularization is used to recover fine structures, especially for images degraded with higher noise levels. Hence apart from a single PDE based approaches, several studies are utilizing coupled PDE to perform image restoration [29]. In this regard, Nitzberg and Shiota [27] introduce a delay in time to calculate diffusivity term as,   1 It = ∇ ∇I , in Ω × (0, T ), (2) 1 + u2 where u in diffusion function is calculated using the relation: ut = ϕ(Gξ ∗ |∇I|2 − u) and T is a fixed time. Here, ξ is a positive smoothing parameter required for the analytical well-posedness of the model. So, the term u updates the edge information from the gradient of the updated image and past information of itself. Hence, instead of spatial regularization, u represents a temporal regularization in the formation of the diffusion term. This coupled PDE system is quite similar to the Perona-Malik model [10] as there is no spatial smoothing term. In [27], the authors show that the coupled PDE system admits a unique classical solution (I, u) in any dimension and satisfy the maximum principle. They also explain that the system does not produce spurious features. A similar idea of coupled PDE system using time-delay regularization is proposed by Luo et al. [26], for image denoising. They propose the estimation of a better edge map by substituting the isotropic diffusion by a nonlinear diffusion equation, to produce a regularized version of the diffusion function for image smoothing. Subsequently, Guo et al. [25] propose and analyze the well-posedness of a reaction-diffusion system for image denoising. They deploy the H −1 norm instead of L1/L2 norm with the total variation model [30], to preserve the oscillatory and texture patterns, as follows, Z (3) |∇I| dΩ + λ||J − I||2H −1 (Ω) . EH −1 = Ω R Here, = Ω |∇∆ (.)| dΩ. This model is efficient for denoising and decomposing the image into cartoon plus texture. Now evolution equation of the Euler-Lagrange equation of the above functional, in which two PDEs are interacting with each other to calculate noise fidelity term, can be given as,    ∇I It = ∇ − 2λv, in Ω × (0, T ), |∇I| (4)  vt = ∆v − (J − I), in Ω × (0, T ). ||.||2H −1 (Ω) −1 2 In the above model, another PDE is used to calculate the fidelity (v) term efficiently. The main difficulty with this 1 model is that the term |∇I| is singular when |∇I| = 0. However, the numerical simulations involve the regularization of |∇I| with |∇I| + ε. A good survey on the existing CPDE approaches can be found in [29]. The objective of the present work is to systematically develop a new framework of coupled nonlinear diffusion model. The present approach differs from previous existing CPDE based works as it utilizes separate evolution equations to calculate edge strength function as well as the data fidelity term. The originality of the proposed coupled non-linear PDEs algorithm lies in the combined effect of: i. diffusion coefficient with space-time regularization instead of space regularization only; ii. H −1 norm based time-varying fidelity term. 2 A PREPRINT - AUGUST 8, 2019 Moreover, due to the interest in theoretical research, we establish the well-posedness of the proposed CPDE model (10)-(13). To prove the well-posedness, we use a time discretization method along with essential a-priori estimates and classical results of compact inclusion in Sobolev spaces [31]; by utilizing some important inequalities, lemmas and theorems [32] and then tends the time discretization parameter to zero. This imparts a major role of this study as it is very important for the numerical computation. Further, to simulate the image denoising, we apply a higher order accurate generalized weighted average finite difference scheme with an advanced iterative solver (Hybrid Bi-Conjugate Gradient Stabilized method) to solve the algebraic system of equations generated from the discretization process [33]. The image quality of denoised images using the proposed CPDE model is compared against several existing nonlinear diffusion models, variational based models and CPDE based models. In practice, our calculations indicate that the proposed approach allows effective image smoothing on fine scales even when the images are degraded with a higher level of noise. The rest of the paper is arranged in the following sequence. Section 2 describes the mathematical formulation of the proposed coupled nonlinear diffusion model. In section 3 we establish the existence and uniqueness of a weak solution to the proposed PDE model. Section 4 shows an appropriate numerical realization in terms of higher order accurate implicit finite difference scheme. Numerical validation of the proposed approach using experimental study is carried out in section 5. At the end, we summarize our observations in section 6. 2 Proposed Coupled Diffusion Model The main aim of this work is to present a novel approach for image smoothing by the evolution of coupled non-linear PDEs. Among the existing CPDE based models mentioned in the previous section, the choices of diffusion function and fidelity term include the space and time regularization. Numerical simulations and experimental results depict that the space-time regularization plays a important role in the quality of recovered images degraded by noise. We note that all of these existing coupled PDEs described in the last section utilizes the time-delay regularization to calculate either diffusion coefficient or fidelity term, and sometimes fail to preserve significant contents of the image. Thus, one issue of these smoothing algorithms is how to design a non-linear diffusion algorithm, which can remove the noise and preserve the small features, simultaneously. For this purpose, inspired by the impressive performance of time-delay regularization in image smoothing, we proposed the following energy functional: Z 1 minEu (I) = g(u)|∇I|2 dΩ + λ||I0 − I||2H −1 (Ω) , (5) 2 Ω where I0 represents the high pass filtered version of image J and λ is weight parameter. Generally, λ is chosen to be inversely proportional to the variance of the noise in the given image. The first term represents a regularizing term producing a smooth varying variable function. The term g(u) serves the purpose of detecting the edges in the image, and u is the edge variable. The second term forces the denoised solution I to be a close approximation to I0 and referred to as data fidelity term. In the data fidelity term, we have replaced L2 norm of (I0 − I) by H −1 norm to preserve the oscillatory and texture information more appropriately in a denoised image. The solution of the above minimization problem (5) can be given by the steady-state of the Euler-Lagrange equation for Eu (I): 1 ∇(g(u)∇I) = ∆−1 (I0 − I). 2λ (6) This is equivalent to following equations, ∇(g(u)∇I) − 2λv = 0, ∇(∇v) − (I0 − I) = 0.  (7) To this end, we consider the function u as the following, ut = α(|∇Iξ |2 − u + β2 ∆u), 2 (8) where α > 0 and β > 0 are parameters to be specified; average time delay (α−1 ) and amount of the spatial smoothing (β) control the deformation and smoothness of the edge information and chosen as described in [27]. And, Iξ = Gξ (X) ∗ I where two dimensional Gaussian function Gξ (.) at each pixel, X ∈ Ω is adopted as,    √ −1 −|X|2 Gξ (X) = ξ 2π , (9) exp 2ξ 2 3 A PREPRINT - AUGUST 8, 2019 where “ ∗ ” denote the convolution operator and ξ > 0 is the standard deviation of the function. Finally, our proposed model in which a set of PDEs interacting with each other, can be expressed as follows, ∂I = ∇(g(u)∇I) − 2λv, in ΩT := (0, T ) × Ω, ∂t β2 ∂u = α(h(|∇Iξ |2 ) − u + △ u), in ΩT , ∂t 2 ∂v = ∇(∇v) − (I0 − I), in ΩT , ∂t I(0, x) = I0 (x), v(0, x) = 0, u(0, x) = Gξ ∗ |∇I0 |2 , in Ω. (10) (11) (12) (13) In the above model, diffusion coefficient g(u) is chosen as g(u) = 1 1+ |Gξ ∗u| k2 , (14) where k > 0 is a threshold parameter, and u represents edge strength at each iteration. Note that we replace |∇Iξ |2 in (8) by h(|∇Iξ |2 ), where h being a sort of truncation, see [24]. Our proposed model is different from the models reported earlier in [29], as it yields two separate evolution equations (PDEs) to handle the diffusion coefficient and data fidelity term. These two evolution equations are responsible for stopping the image smoothing at the edges as well as textures and forces the resulting solution to be a close approximation to the given initial image. The proposed model is in a spirit similar to the reaction-diffusion model equation in [25]; however, the diffusion function is very different in our case. In (10)-(13), the smoothing equation (10) has different edge variable u and fidelity term v, obtained from two different PDEs. The data fidelity term between I and I0 can be handled by function v, which can be obtained from equation (12). Whereas the function u in the diffusion coefficient can be calculated from equation (11). In summary, the proposed model achieves a suitable edge map and fidelity between noisy image and restored image at each iteration, which ultimately leads to quality denoising results. Hence, the proposed model provides a potential approach for image noise removal and enhancement. 3 Existence and uniqueness of weak solution In this section, we establish the well-posedness of (10)-(13) with Dirichlet’s boundary condition for the fidelity variable and Neumann’s boundary condition for other two variables. Furthermore, for simplicity we choose all the constants involved in the equations (10)-(13) equals to 1. Note that, h : R+ → R+ is a bounded, Lipschitz continuous function with Lipschitz constant ch such that δ ≤ h(ũ) ≤ 1 ∀ ũ ∈ R+ , (15) for some δ > 0. Also from (14) we observe that, g : R → R+ is a bounded, decreasing and Lipschitz continuous C function with Lipschitz constant k2ξ . Moreover, g(0) = 1 and lim g(u) = 0. u→+∞ 3.1 Technical framework We denote by H k (Ω), k is a positive integer, the set of all functions I : Ω → R such that I and its distributional P2 ∂mI 2 k derivatives ∂x m of order |m| = j=1 mj ≤ k all belongs to L (Ω). H (Ω) is a Hilbert space endowed with the norm  X Z ∂ m I 2 1/2 kIkH k := ||I||H k (Ω) = dx . m Ω ∂x |m|≤k We denote by Lp (0, T ; H k (Ω)), p > 1 and k is a positive integer, the set of all measurable functions I : [0, T ] → H k (Ω) such that Z T 1/p kIkLp(0,T ;H k (Ω)) := ||I(t)||pH k (Ω) dt < ∞. 0 We denote by H (Ω) the dual of H (Ω), and H (Ω) the dual of H01 (Ω). For any f ∈ H 1 (Ω)′ , we define a norm as n o ||f ||H 1 (Ω)′ = sup hf, ui : u ∈ H 1 (Ω) , ||u||H 1 (Ω) ≤ 1 . 1 ′ 1 −1 4 A PREPRINT - AUGUST 8, 2019 We introduce the solution space W (0, T ) := W1 (0, T ) × W1 (0, T ) × W2 (0, T ) for the problem (10)-(13), where n o ∂w W1 (0, T ) = w : w ∈ L2 (0, T ; H 1(Ω)) , ∈ L2 (0, T ; (H 1 (Ω))′ ) , ∂t o n ∂w 2 1 ∈ L2 (0, T ; H −1 (Ω)) . W2 (0, T ) = w : w ∈ L (0, T ; H0 (Ω)) , ∂t Note that, the space Wi (0, T ) (i = 1, 2) is a Hilbert space for the graph norm, see [34]. Definition 3.1 (Weak solution) Let I0 ∈ H 1 (Ω). We say that a triple (I, u, v) is a weak solution of (10)-(13) if (i) I, u ∈ W1 (0, T ), and v ∈ W2 (0, T ). (ii) For all Ψ ∈ H 1 (Ω) and Φ ∈ H01 (Ω), there holds Z T Z Z ∂I vΨ dx dt = 0 , g(u) ∇I · ∇Ψ dx dt + 2 , Ψ dt + ∂t 0 ΩT ΩT Z T Z Z ∂u uΨ dx dt h(|∇Iξ |2 )Ψ dx dt + , Ψ dt − ∂t 0 ΩT ΩT Z ∇u · ∇Ψ dx dt = 0 , + Z 0 T ∂v , Φ dt + ∂t Z ΩT ∇v · ∇Ψ dx dt + Z (16) (17) ΩT ΩT (I0 − I)Ψ dx dt = 0 . (18) (iii) (13) holds. 3.2 Discretized system and existence of its weak solution To prove existence of a weak solution of (10)-(13), we use semi discretization in time with a semi-implicit Euler T method. Let 0 = t0 < t1 < t2 · · · < tN = T be a uniform partition of [0, T ] with time-step size τ = N for some + 2 1 1 N ∈ N , i.e, tn = τ n for 0 ≤ n ≤ N . For simplicity, in the following we write L , H , H0 instead of L2 (Ω), H 1 (Ω), H01 (Ω) respectively. We may then consider time discretization of (10)-(13) as follows: for 0 ≤ n ≤ N − 1 , ∀ Ψ ∈ H 1 , and ∀ Φ ∈ H01 iterate: (i) compute un+1 ∈ H 1 such that


1 (un+1 − un ), Ψ L2 + ∇un+1 , ∇Ψ L2 + un+1 , Ψ L2 τ
= h(|∇Gξ ∗ In |2 ), Ψ L2 , (ii) compute vn+1 ∈ H01 (iii) compute In+1 ∈ H 1 such that


1 1 vn+1 , Φ L2 + ∇vn+1 , ∇Φ L2 = vn − (I0 − In ), Φ L2 τ τ


1 1 In+1 , Ψ L2 + g(un )∇In+1 , ∇Ψ L2 = In − 2vn , Ψ L2 . τ τ (19) (20) (21) Thanks to (15), the existence of un+1 resp. vn+1 in step (i) resp. in step (ii) easily follows from Lax-Milgram lemma. Moreover, these solutions are unique. To show the well-posedness of (21), we use standard Lax-Milgram lemma. Define a continuous bilinear form C : H 1 × H 1 → R, and a linear bounded functional ℓ : H 1 → R as


1 1 C(Θ, Ψ) = Θ, Ψ L2 + g(un )∇Θ, ∇Ψ L2 ; ℓ(Ψ) = In − 2vn , Ψ L2 . (22) τ τ Then, by Lax-Milgram lemma, there exists a unique In+1 ∈ H 1 satisfying (21) provided C is H 1 -coercive. Indeed, thanks to Gagliardo-Nirenberg inequality,there exists a constant C > 0, depending only on Ω and Gξ , such that 1 ≤ g(un ) , (23) C||un ||H 1 1+ k2 o n 1 1 2 2 e e := min 1 , yielding ||∇Ψ|| ≥ C||Ψ|| where C and hence C(Ψ, Ψ) ≥ τ1 ||Ψ||2L2 + 2 1 L H τ C||un ||H 1 C||un ||H 1 1+ 1+ k2 k2 the H 1 -coercivity of C. 5 A PREPRINT - AUGUST 8, 2019 3.3 A-priori estimates Let 0 ≤ n ≤ N − 1. Choose Ψ = un+1 as test function in (19) and use the algebraic identity ha, a − bi =
1 |a|2 − |b|2 + |a − b|2 , 2 ∀ a, b ∈ Rp (p ≥ 1) , (24) Cauchy-Schwarz and Young’s inequalities, (15), and then sum over n = 0, 1, 2, . . . , j − 1 with 1 ≤ j ≤ N . The result is ||uj ||2L2 + j−1 X n=0 ||un+1 − un ||2L2 j−1 X ≤ ||u0 ||2L2 + |Ω| n=0 + 2τ j−1 X n=0 ||∇un+1 ||2L2 τ ≤ ||u0 ||2L2 + T |Ω| . (25) Hence, there exists a constant C1 > 0, independent of τ > 0, such that max ||uj ||2L2 + 1≤j≤N N −1 X n=0 ||un+1 − un ||2L2 + 2τ N −1 X n=0 ||∇un+1 ||2L2 ≤ C1 . (26) Similarly, by choosing (formally) Ψ = −∆un+1 in (19) and using the integration by parts formula, (24), CauchySchwarz and Young’s inequalities, (15), and then summing over n = 0, 1, 2, . . . , j − 1 with 1 ≤ j ≤ N , we obtain max 1≤j≤N ||∇uj ||2L2 + N −1 X n=0 ||∇(un+1 − un )||2L2 +τ N −1 X n=0 ||∆un+1 ||2L2 ≤ C2 , (27) where C2 > 0 is a constant, independent of τ > 0. In view of (26) and (27), there exists a constant C3 > 0, independent of τ , such that (28) max kuj kH 1 (Ω) ≤ C3 . 0≤j≤N Therefore, from (23), the positive lower bound for g(un ) is given by ν := 1 1+ CC3 k2 (29) ≤ g(un ) (0 ≤ n ≤ N ) . Again, one can use the test function Φ = vn+1 and Φ = −∆vn+1 (formally) in (20), and proceed as above (under the cosmetic changes) to obtain (1 ≤ j ≤ N ): ||vj ||2L2 + j−1 X n=0 ||vn+1 − vn ||2L2 + τ ≤ ||∇vj ||2L2 + j−1 X n=0 j−1 X n=0 kv0 k2L2 ||∇vn+1 ||2L2 + ||∇(vn+1 − vn )||2L2 + τ CP2 T kI0 k2L2 j−1 X n=0 + τ CP2 j−1 X n=0 kIn k2L2 , (30) ||∆vn+1 ||2L2 ≤ ||∇v0 ||2L2 + T kI0 k2L2 + τ j−1 X n=0 ||In ||2L2 , (31) where CP is the Poincaré constant. Note that, because of coupled system, we are not able to find the bound of vn ’s from (30) and (31). May be we need to combine with the estimates coming from (21) by choosing appropriate test 6 A PREPRINT - AUGUST 8, 2019 functions. Taking Ψ = In+1 in (21) and using (24), (29), and (30), we have ||Ij ||2L2 + j−1 X n=0 ||In+1 − In ||2L2 + 2τ ν ≤ ||I0 ||2L2 + 2τ j−1 X n=0 j−1  X n=0 ≤ C4 + C5 τ k=0 n=0 ||∇In+1 ||2L2 ||In+1 ||2L2 + 2τ j−1 X j−1 X kv0 k2L2 + CP2 T kI0 k2L2 + τ CP2 n X k=0 kIk+1 k2L2  kIk+1 k2L2 , (32) for some constants C4 , C5 > 0, only depend on T, CP , v0 , I0 . We now apply discrete Gronwall’s lemma (implicit form) and obtain the following: there exists τ1 > 0 such that for all time step sizes 0 < τ < τ1 , max ||Ij ||2L2 + 0≤j≤N N −1 X n=0 ||In+1 − In ||2L2 + 2τ ν N −1 X n=0 ||∇In+1 ||2L2 ≤ C6 , (33) where C6 > 0 is independent of τ > 0. We now use (33) in (30) and (31), and get max ||vj ||2H 1 + 0≤j≤N N −1 X n=0 ||vn+1 − vn ||2H 1 + τ  ||∇vn+1 ||2L2 + k∆vn+1 k2L2 ≤ C7 , N −1  X n=0 (34) for some constant C7 > 0, independent of τ with 0 < τ < τ1 . Again one may choose Ψ = −∆In+1 (formally) in (21) and use the integration by parts
formula, (24), Cauchy-Schwarz and Young’s inequalities, (29), (28), (34) and the fact that k∇g(un )kL∞ ≤ C ξ, Ω, k kun kH 1 to get, after the application of discrete Gronwall’s lemma (implicit form), max ||∇Ij ||2L2 + 1≤j≤N N −1 X n=0 ||∇(In+1 − In )||2L2 + τ j−1 X n=0 ||∆In+1 ||2L2 ≤ C8 , (35) with C8 > 0, independent of τ with 0 < τ < τ2 for some τ2 > 0. Note that, for τ0 = min{τ1 , τ2 }, (33), (34) and (35) all hold true. Putting things together, we arrive at the following theorem. Theorem 3.1 Let I0 ∈ H 1 . Then there exist constants τ0 > 0 and C > 0 such that for all time step sizes τ > 0 with 0 < τ < τ0 , there holds   max kuj k2H 1 + kvj k2H 1 + kIj k2H 1 1≤j≤N + N −1  X n=0 +τ ||un+1 − un ||2H 1 + ||vn+1 − vn ||2H 1 + ||In+1 − In ||2H 1  k∆un+1 k2L2 + k∆vn+1 k2L2 + k∆In+1 k2L2 ≤ C , N −1  X n=0  (36) where uj , vj and Ij solves (19), (20) and (21) respectively. 3.4 Continufication and existence of weak solution Let {tn }N n=0 be a uniform partition of [0, T ] with time-step size τ as described in Subsection 3.2. For any sequence xn+1 − xn {xn } ⊂ X, where X is a Banach space, we define the difference quotient dt xn+1 = for 0 ≤ n ≤ N − 1. τ The globally time interpolant Xτ ∈ C(X) of {xn } is defined via Xτ (t) := t − tn tn+1 − t xn+1 + xn , τ τ 7 ∀ t ∈ (tn , tn+1 ] . A PREPRINT - AUGUST 8, 2019 − Moreover, we define the piecewise constant in time interpolants X+ τ (t) and Xτ (t) as follows: X+ τ (t) := xn+1 , X− τ (t) := xn , ∀ t ∈ (tn , tn+1 ] , with x−1 = x0 and xN +1 = 0. With the above notation, we now show the boundedness of the sequence {∂t Uτ }τ >0 in L2 (H 1 (Ω)′ ). To do so, let Ψ ∈ H 1 (Ω)′ \ {0}, and t ∈ (tn , tn+1 ] where 0 ≤ n ≤ N − 1. Then from (19), we have, by using Cauchy-Schwarz inequality and (15),  1  ∂t Uτ , Ψ ≤ |Ω| 2 + max kun+1 kH 1 kΨkH 1 , (37) 0≤n≤N −1 and hence, thanks to Theorem 3.1, −1  NX |Ω| + T ||∂t Uτ ||2L2 (H 1 (Ω)′ ) ≤C τ n=0 max 0≤n≤N −1  kun+1 k2H 1 ≤ C . (38) Similarly, the sequences {∂t Vτ }τ >0 and {∂t Iτ }τ >0 are bounded in L2 (H −1 (Ω)) and L2 (H 1 (Ω)′ ) respectively. Moreover, there exists a constant C > 0 such that for all 0 < τ < τ0 , ||∂t Vτ ||2L2 (H −1 (Ω)) + ||∂t Iτ ||2L2 (H 1 (Ω)′ ) ≤ C . (39) Thanks to (38), (39) and Theorem 3.1, by using classical results of compact inclusion in Sobolev spaces [31], there exists (I, u, v) ∈ W1 (0, T ) × W1 (0, T ) × W2 (0, T ) such that along a subsequence (still we denote it by the same index) the following hold:  + − ∗ Iτ , Iτ ⇀ I in L∞ (H 1 ) ; ∂t Iτ ⇀ ∂t I in L2 (H 1 (Ω)′ ) ,     Iτ , Iτ± ⇀ I in L2 (H 1 ) ; Iτ , Iτ± → I in L2 (L2 ) ,      |∇Gξ ∗ Iτ− |2 → |∇Gξ ∗ I|2 in L2 (L2 ) and a.e. (0, T ) × Ω ,    2 2 2 2  h(|∇Gξ ∗ Iτ | ) → h(|∇Gξ ∗ I| ) in L (L ) , + − ∗ ∞ 1 (40) Uτ , Uτ ⇀ u in L (H ) ; ∂t Uτ ⇀ ∂t u in L2 (H 1 (Ω)′ ) , (τ → 0) ,   ± 2 1 ± 2 2  U , U ⇀ u in L (H ) ; U , U → u in L (L ) , τ τ  τ τ    g(Uτ+ ) → g(u) in L2 (L2 ) ,     V + , V − ⇀∗ v in L∞ (H01 ) ; ∂t Vτ ⇀ ∂t v in L2 (H −1 (Ω)) ,   τ ±τ Vτ , Vτ ⇀ v in L2 (H01 ) ; Vτ , Vτ± → v in L2 (L2 ) . For any Ψ ∈ H 1 and Φ ∈ H01 , we rewrite (19), (20) and (21) in terms of Uτ , Uτ+ , Uτ− , Vτ , Vτ+ , Vτ− , Iτ , Iτ+ , and Iτ− , integrate with respect to the time variable and make use of (40) to pass to the limit as τ → 0 in the resulting variational formulation to arrive at (16), (17) and (18). In other words, the triple (I, u, v) is a weak solution of (10)-(13). 3.5 Uniqueness of weak solution We use standard methodology [32] to prove the uniqueness of weak solution of (10)-(13). Let (I1 , u1 , v1 ) and (I2 , u2 , v2 ) be two sets of solution for the system (10)-(13) with I1 6= I2 , u1 6= u2 and v1 6= v2 . Let e := u1 − u2 , Ve := v1 − v2 and gi = g(ui ) (i = 1, 2). Then, the following equations hold in I˜ := I1 − I2 , U the sense of distribution: ∂ I˜ ˜ = ∇((g1 − g2 )∇I2 ) − 2Ve , − ∇(g1 ∇I) in ΩT , (41) ∂t   e ∂U e = h(|∇Gξ ∗ I1 |2 ) − h(|∇Gξ ∗ I2 |2 ) − U e, − ∆U in ΩT , (42) ∂t ∂ Ve − ∆Ve = I˜ , in ΩT , (43) ∂t ˜ x) = 0 = U e (0, x) = Ve (0, x) , I(0, in Ω . (44) Note that 1 ν1 := 1+ Cku1 kL∞ (H 1 ) k2 ≤ g(u1 ) ; ||(g1 − g2 )||L∞ ≤ Cξ ||(u1 − u2 )||L2 . 8 A PREPRINT - AUGUST 8, 2019 Similar to a-priori estimates, we have for a.e. t ∈ (0, T ) with C̃I := k∇Ik2L∞ (L2 ) , 1 d ˜ 2 ˜ 2 2 + kVe k2 2 kIkL2 ≤ C̃I C(ν1 )||(g1 − g2 )||2L∞ + kIk L L 2 dt ˜ 2 2 + ||Ve ||2 2 . e U e ||2 2 + ||I|| ≤ C|| L L L e and integrating over Ω, we get Similarly, by multiplying (42) by U d e 2 e ||2 2 . e ||2 2 ≤ C||h(|∇Gξ ∗ I1 |2 ) − h(|∇Gξ ∗ I2 |2 )||2 2 + ||U ||U ||L2 + 2||∇U L L L dt By using Lipschitz continuity of h along with Gagliardo-Nirenberg inequality, we have and hence ˜ 22 , ||h(|∇Gξ ∗ I1 |2 ) − h(|∇Gξ ∗ I2 |2 )||2L2 ≤ C(ch , ξ, I1 , I2 )kIk L In a similar way, one can easily deduce that d e 2 ˜ 22 . ||U ||L2 ≤ C||I|| L dt (45) (46) (47) (48)   d e 2 ˜ 2 2 + ||Ve ||2 2 . (49) ||V ||L2 ≤ C ||I|| L L dt e ≡ 0 and Ve ≡ 0. In other words, We now add (45)-(49) and then apply Gronwall’s inequality to infer that I˜ ≡ 0, U weak solution of the proposed model (10)-(13) is unique. 4 Numerical Approximation In this section, we discuss the numerical implementation of the proposed model (10)-(13). To solve the model numerically, we have applied the generalized weighted average finite difference scheme which has a combined nature of the forward Euler method and the backward Euler method, at nth and (n + 1)th step, respectively. The image variable I, edge variable u and fidelity term v are calculated at each iteration. In the derivation of finite difference formulas, h and τ are considered as spatial step size and time step size, respectively. n Let I(xi , yj , tn ) express the gray level of the image plane Ii,j , where xi = ih, yj = jh and tn = nτ . Also, we consider 0 0 2 vi,j = 0 and u0i,j = Gξ ∗ |∇Ii,j | . Derivatives are approximated by central difference formula as, n n n n Ii+1,j − Ii−1,j Ii,j+1 − Ii,j−1 ∂I ∂I ≈ , ≈ , ∂x 2h ∂y 2h n n n n n n Ii+1,j − 2Ii,j + Ii−1,j Ii,j+1 − 2Ii,j + Ii,j−1 ∂2I ∂2I ≈ , ≈ . ∂2x h2 ∂2y h2 The discrete form of (10)-(13) can be written as;

τ τ n+1 n+1 n n Ii,j − , (∇(g(u)∇I)n+1 = Ii,j + (∇(g(u)∇I)ni,j − 2λvi,j i,j − 2λvi,j 2 2   n un+1 ψ2 i,j − ui,j n+1 , + = ϕ h(|∇Iξni,j |2 ) − un+1 ∆u i,j i,j τ 2 n+1 n
vi,j − vi,j n+1 n 0 , = ∆vi,j − Ii,j − Ii,j τ where, n n n n n n wi+1,j − 2wi,j + wi−1,j wi,j+1 − 2wi,j + wi,j−1 + , h2 h2 for w = u, v. Moreover the diffusion term in (10)-(13) approximated using the central difference scheme as, n ∆wi,j = ∇(g(u)∇I)i,j = 1 (g 1 (Ii+1,j − Ii,j ) − gi− 21 ,j (Ii,j − Ii−1,j )+ h2 i+ 2 ,j gi,j+ 12 (Ii,j+1 − Ii,j ) − gi,j− 12 (Ii,j − Ii,j−1 )) , 9 (50) (51) (52) (53) (54) A PREPRINT - AUGUST 8, 2019 with the Neumann boundary conditions: Ii,0 − Ii,1 = 0, Ii,N −1 − Ii,N = 0, for 0 ≤ i ≤ M, I0,j − I1,j = 0, IM−1,j − IM,j = 0, for 0 ≤ j ≤ N. Further, to solve the algebraic system of the the form AI n+1 = B, generated from numerical discretization, hybrid Bi-Conjugate Gradient Stabilized solver has been used [33, 37]. Apart from the numerical discretization of (10)-(13), a stopping criterion is needed to terminate the diffusion process. To achieve this goal, we start with an initial image I0 and apply the system (10)-(13) repeatedly. This results in a family of smoother images I(t, x); t > 0, which depicts refined versions of I0 . After sufficient iterations, changes between two consecutive iterations become redundant. At this point, the convergence of the iterative procedure has been achieved. Now, we use the following measure as a stopping criterion, 2 ||I k+1 − I k ||2 2 ||I k ||2 (55) ≤ ε. th Here, we use ε = 10−4 as a fixed threshold. I k and I (k+1) depict the image planes at the k th and (k + 1) iteration, respectively. Parameter values used for the numerical experiments are mentioned in the caption of each figure. For Algorithm 1 Algorithm of CPDE model for Image Denoising 1: Inputs: I (0) = I0 as initial solution or initial noisy image, v (0) = 0 as initial edge variable and u(0) = Gξ ∗ |∇I (0) |2 , τ as time step, k = 0. 2: Result: Denoised image (I) 3: Calculate I (k+1) , u(k+1) and v (k+1) from     k+1   k τ 0 τλ 1 − τ2 A(I k+1 ) I + 2 (A(I k ) − 2λv k ) I 2 k+1   =  uk + τ ϕh(||∇Iξk ||2 )  0 1 + τ (1 − ψ2 ∆) 0  u v k − τ (I k − I 0 ) v k+1 0 0 1 − τ∆ . 4: Solve the algebraic system (obtained using step 3), using Hybrid BiCGStab solver . 5: Check if 2 ||I k+1 − I k ||2 ≤ ε = 10−4 ||I k ||22 then stop or go to Step 3. 6: end the TV and NLM model, we have reproduced the results by utilizing the fact as mention in [19] and [39] respectively. Apart from the parameters displayed in the captions, we have chosen ϕ = 1 and ξ = 1 for the present model and a uniform time step size τ = 0.1 for the current as well as the other discussed approaches. 5 Results and Discussion The proposed CPDE model was employed using the finite difference scheme given in the previous section. Image denoising using (10)-(13) was compared with the results of other state-of-art models available in the literature. Especially, TV model [19], NS model [27], Luo model [26], RD model [25], PV model [28], SYS model [38] are considered for the comparison. Apart from the diffusion-based approaches, non-iterative approach (Non-Local Means method) [39] is also used for comparison with the present method. Since the space-time regularization based proposed model is claimed to be an improvement over the existing CPDE models, our main aim was to compare the edge detection and denoising results with these models (NS Model, and Luo model). In this process, the considered non-linear diffusion models are solved by the existing numerical schemes. And also, to stop the iterative process of each smoothing algorithm discussed stopping criteria equation (55) is employed. The effectiveness of results was evaluated through several standard gray level test images which are degraded with additive Gaussian noise of zero mean and different levels of standard deviations. We have artificially added additive Gaussian noise of different standard deviations ranging from 20 to 50 by using our MATLAB program. Also, to evaluate the efficiency of our model, quantitative comparisons in terms of Peak signal to noise ratio (PSNR)[2], mean structural similarity index measure (MSSIM)[40], Gradient of peak signal to noise ratio (PSNRGrad )[41], Improvement in signal to noise ratio (ISNR)[2] have been shown with existing models. A higher value of quantity metrics suggests that the filtered image is closer to the noise-free image. 10 A PREPRINT - AUGUST 8, 2019 ˆ (a). Peak signal to noise ratio (PSNR)[2] can measure the match between the clean (I) and denoised image (I), ˆ = 10log PSNR(I, I) 10 MN|max(I) − min(I)|2 . 2 ||Iˆ − I|| 2 L (b). In addition to PSNR, gradient of peak signal to noise ratio (PSNRGrad ) is also used to measure the match between the derivatives of reconstructed and true image, defined as PSNRGrad = 1 (PSNR(Ix , Iˆx ) + PSNR(Iy , Iˆy )) , 2 ˆ where (Ix , Iy ) and (Iˆx , Iˆy ) are the derivatives of ground truth image (I) and denoised image (I). (c). Structural similarity index (SSIM) [40], is used to calculate the similarity between structure of clean and reconstructed images, and can be given as, SSIM(X, Y ) = (2µx µy + c1 )(2σxy + c2 ) , (µ2x + µ2y + c1 )(σx2 + σy2 + c2 ) here µx , µy , σx2 , σy2 , σxy are the average, variance and covariance of X and Y , respectively. The variables c1 and c2 are used to stabilize the division with weak denominator. To find overall performance of the image, we used a mean SSIM index (MSSIM), N MSSIM(X, Y ) = 1X SSIM(Xi , Yi ) , N i=1 where X and Y are used to represent the original and the denoised images, respectively. N is the number of local windows in the image, and Xi and Yi are the image contents at the ith local window. (d). Improvement in signal to noise ratio (ISNR) [2], is simply the difference between the improved and the original signal to noise ratio,  M N  PP 2 (I − I ) 0    i=1j=1  ISNR = 10log10  M N . P P  2 (I − It ) i=1 j=1 To test the effectiveness of the proposed model, figure[1] presents the denoising results for Walkbridge test image that contain additive Gaussian noise of σ = 40, which make the features hard to visualize. The restored outputs obtained from the state-of-the-art diffusion models depict that the denoised images are not satisfactory and the texture information of the image is ruined. From the figure[1e], it is easy to perceive that the proposed model can effectively eliminate all noise particles and preserve the original structure of the image. Further, to showcase the efficiency of our model, in figure 2, we compare the visual quality of the mosaic image computed by the present model along with the other existing models. Here we present the denoised images and 3D surface plots of the denoised images side by side for better comparison. From the figure, it can be observed that our model works better in terms of noise removal as compared to other models. From the 3D surface plots of the denoised images, one can see that our model leave fewer fluctuations with compare to the other models, it indicates that the proposed model not only removes noise efficiently but also preserves the fine structure as compared to the other models. To further confirm the ability of the proposed model, in figure[3], we present the comparison of the denoised Pirate image which was initially degraded with additive Gaussian noise with very low SNR value 3.70. From figure[3], it is clear that the proposed algorithm retains more texture information in addition to noise reduction, which makes our filtered image better in comparison to other models considered here. Along with the qualitative analysis using full image surface, we also study the quality of the resultant images considering a slice of the Boat image and the Livingroom image with different noise levels. Figure 4 shows the signals of the original, noisy and restored images obtained by the proposed CPDE model and other discussed CPDE models. From these figures, it is easy to conclude that the restored signals computed by the proposed model are more closure to the clean signals in comparison to other discussed models. In addition to qualitative comparisons, the quantitative results, in terms of PSNR and MSSIM values for different test images as well as noise levels (σ = 20, 30, 40, 50), are shown in tables [1]-[4]. To ease the comparison, the highest 11 A (a) Clear PREPRINT - AUGUST 8, 2019 (b) Noisy (c) PV (d) SYS (e) Proposed Figure 1: (a) Clear image; (b) Noisy image with Gaussian noise of mean 0.0 and σ = 40; Denoised image using (c) PV Model, λ = 0.01, K = 8; (d) SYS Model, λ = 4, K = 8; (e) Proposed model, ψ = 0.1, k = 4.45. Table 1: Comparison of MSSIM and PSNR values for various approaches and proposed model. Clean image is degraded by additive Gaussian noise of mean 0.0 and σ = 20. Images Measure TV[19] NS[27] Luo[26] RD[25] NLM[39] PV[28] SYS[38] Proposed Boat MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSI PSNR MSSIM PSNR MSSIM PSNR 0.8232 28.10 0.8539 28.37 0.8790 28.51 0.8958 27.51 0.8615 28.54 0.9077 26.79 0.8007 30.35 0.9699 34.28 0.8883 28.81 0.8996 28.07 0.8698 28.62 0.9001 27.71 0.9053 29.40 0.9141 26.83 0.9191 33.83 0.9091 32.71 0.8774 29.40 0.8775 28.65 0.8855 28.68 0.8909 27.77 0.8764 29.20 0.9070 27.05 0.8257 31.31 0.9270 33.66 0.8886 29.61 0.8788 28.61 0.8849 28.80 0.8997 27.95 0.8714 29.07 0.9101 26.99 0.8206 31.29 0.9600 34.44 0.9105 30.21 0.9140 29.27 0.8870 29.11 0.8684 27.62 0.8908 29.69 0.8747 26.72 0.9243 33.64 0.9588 34.12 0.8900 29.69 0.9060 28.87 0.8903 28.87 0.9015 28.24 0.8936 29.51 0.9093 26.99 0.9262 34.47 0.9675 33.38 0.8943 29.84 0.9074 29.08 0.8912 28.93 0.9029 28.37 0.8943 29.53 0.9089 27.02 0.9225 34.59 0.9290 32.54 0.9219 30.24 0.9268 29.57 0.9019 29.22 0.9034 28.46 0.9055 29.90 0.9143 27.24 0.9396 35.03 0.9780 34.74 Lake Livingroom Mandril Pirate Walkbridge Woman Mosaic Table 2: Comparison of MSSIM and PSNR values for various approaches and proposed model. Clean image is degraded by additive Gaussian noise of mean 0.0 and σ = 30. Images Measure TV[19] NS[27] Luo[26] RD[25] NLM[39] PV[28] SYS[38] Proposed Boat MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR 0.7157 25.28 0.7644 25.27 0.7859 25.91 0.8188 24.94 0.7745 26.18 0.8391 24.60 0.6834 27.58 0.9558 31.63 0.6923 25.71 0.7676 25.41 0.8510 27.28 0.7633 25.46 0.8348 27.05 0.8194 24.07 0.8559 31.16 0.8827 30.38 0.7324 25.83 0.7751 25.72 0.7922 26.14 0.8262 25.55 0.7765 26.26 0.8431 24.91 0.6805 27.47 0.8800 30.68 0.7327 26.01 0.7820 26.01 0.7969 26.40 0.8268 25.52 0.7905 26.86 0.8453 24.99 0.7108 28.58 0.9467 31.65 0.8568 28.02 0.8705 27.38 0.8229 26.90 0.7873 25.43 0.8346 27.72 0.8015 24.96 0.8861 31.38 0.9477 31.69 0.8360 27.79 0.8576 27.12 0.8273 26.97 0.8309 26.13 0.8392 27.86 0.8462 25.19 0.8963 32.58 0.9520 31.10 0.8400 27.96 0.8637 27.33 0.8321 27.02 0.8395 26.27 0.8400 27.81 0.8508 25.22 0.8922 32.64 0.8764 30.24 0.8841 28.47 0.8956 27.87 0.8515 27.39 0.8454 26.47 0.8575 28.20 0.8566 25.47 0.9200 33.25 0.9659 31.97 Lake Livingroom Mandril Pirate Walkbridge Woman Mosaic 12 A (a) Clear (c) PV PREPRINT - AUGUST 8, 2019 (b) Noisy (d) SYS (e) Proposed (f) Noisy (g) PV (h) SYS (i) Proposed Figure 2: (a) Clear image (b) Noisy image with Gaussian noise of mean 0.0 and s.d. 40 (c) PV Model, λ = 0.01, K = 1; (d) SYS Model, λ = 0.1, K = 4; (e) Proposed Model, ψ = 1, k = 3.75. (f-i) 3D surfaces of images 13 A (a) Original PREPRINT - AUGUST 8, 2019 (b) Noisy (c) PV Model (d) SYS Model (e) Proposed Model Figure 3: (a) Clear image (b)Heavily noised image with low SNR=3.70; Improved SNR values by different models: (c) SNR=15.41; λ = 0.01, K = 11 (d) SNR=15.37; λ = 1, K = 10 (e) SNR=15.53; ψ = 1, k = 14. Table 3: Comparison of MSSIM and PSNR values for various approaches and proposed model. Clean image is degraded by additive Gaussian noise of mean 0.0 and σ = 40. Images Measure TV[19] NS[27] Luo[26] RD[25] NLM[39] PV[28] SYS[38] Proposed Boat MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR 0.6305 23.43 0.6964 23.67 0.6941 23.85 0.7409 23.14 0.6882 24.32 0.7728 23.05 0.5976 25.75 0.9384 29.43 0.6146 23.41 0.7052 23.70 0.6989 23.36 0.7456 23.24 0.6913 24.23 0.7883 23.48 0.7317 27.59 0.8892 28.51 0.6290 23.35 0.6927 23.62 0.6855 23.56 0.7462 23.41 0.6800 24.04 0.7722 23.08 0.5671 25.93 0.8315 28.38 0.6481 24.19 0.7117 24.33 0.7095 24.50 0.7525 23.82 0.7067 25.06 0.7821 23.54 0.6235 26.61 0.9328 29.46 0.8022 26.40 0.8245 25.93 0.7600 25.37 0.7164 24.01 0.7849 26.35 0.7371 23.76 0.8513 29.75 0.9388 29.66 0.7862 26.51 0.8154 25.94 0.7699 25.68 0.7692 24.75 0.7869 26.64 0.7922 24.01 0.8690 31.00 0.9303 29.04 0.7798 26.60 0.8173 26.05 0.7731 25.70 0.7777 24.90 0.7870 26.58 0.7944 24.05 0.8641 30.96 0.8373 28.15 0.8470 27.21 0.8644 26.62 0.7991 26.11 0.7844 25.10 0.8132 27.06 0.7985 24.24 0.9037 31.65 0.9508 29.77 Lake Livingroom Mandril Pirate Walkbridge Woman Mosaic Table 4: Comparison of MSSIM and PSNR values for various approaches and proposed model. Clean image is degraded by additive Gaussian noise of mean 0.0 and σ = 50. Images Measure TV [19] NS [27] Luo [26] RD [25] NLM[39] PV[28] SYS[38] Proposed Boat MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR 0.5648 22.02 0.6404 22.47 0.6265 22.51 0.6707 21.77 0.6275 23.11 0.7150 22.00 0.5318 24.30 0.9203 27.61 0.5793 21.70 0.6657 22.09 0.6424 21.57 0.6177 21.48 0.6486 21.66 0.7046 21.75 0.5702 22.02 0.8404 26.93 0.5519 21.50 0.6232 21.98 0.5981 21.54 0.6765 22.06 0.5969 22.15 0.7025 21.52 0.5626 22.56 0.8662 26.55 0.5813 22.75 0.6540 23.06 0.6429 23.17 0.6846 22.49 0.6450 23.78 0.7243 22.44 0.6143 25.02 0.9189 27.66 0.7523 25.14 0.7828 24.75 0.7090 24.28 0.6505 22.99 0.7423 25.24 0.6838 22.84 0.8150 28.23 0.9144 27.71 0.7425 25.48 0.7747 24.89 0.7245 24.75 0.7102 23.69 0.7483 25.68 0.7414 23.10 0.8430 29.31 0.9078 27.34 0.7333 25.46 0.7710 24.94 0.7229 24.72 0.7173 23.87 0.7454 25.60 0.7428 23.16 0.8354 29.25 0.7730 26.51 0.8159 26.16 0.8334 25.53 0.7562 25.17 0.7261 24.06 0.7762 26.10 0.7483 23.34 0.8830 29.81 0.9333 27.87 Lake Livingroom Mandril Pirate Walkbridge Woman Mosaic 14 A 300 PREPRINT 300 Noisy Signal 250 250 200 200 Intensity Intensity Noisy Signal 150 150 100 100 50 50 0 0 0 50 100 150 200 250 300 0 50 100 Pixel Coordinate 250 300 300 Clean signal PV Model 250 Clean signal PV Model 250 200 Intensity Intensity 200 (b) 300 150 200 150 100 100 50 50 0 0 0 50 100 150 200 250 300 0 50 100 Pixel Coordinate 150 200 250 300 Pixel Coordinate (c) λ = 0.1, K = 5 (d) λ = 0.1, K = 15 300 300 Clean signal SYS Model 250 Clean signal SYS Model 250 200 Intensity Intensity 150 Pixel Coordinate (a) 150 100 200 150 100 50 50 0 0 0 50 100 150 200 250 300 0 50 100 Pixel Coordinate 150 200 250 300 Pixel Coordinate (e) λ = 2, K = 6 (f) λ = 4, K = 9 300 300 Clean Signal Proposed Model 250 Clean Signal Proposed Model 250 200 Intensity Intensity - AUGUST 8, 2019 150 100 200 150 100 50 50 0 0 0 50 100 150 200 250 300 Pixel Coordinate 0 50 100 150 200 250 300 Pixel Coordinate (g) ψ = 1, k = 4 (h) ψ = 1, k = 4.45 Figure 4: (a) The 219th slice of Boat image corrupted with Gaussian noise of mean 0.0 and σ = 40; (b) The 95th slice of Livingroom image corrupted with additive Gaussian noise of mean 0.0 and σ = 50; denoised using (c-d) PV model; (e-f) SYS model, and (g-h) Proposed model. 15 Table 5: Comparison of ISNR and PSNRgrad values for various approaches and proposed model. Clean image is degraded by additive Gaussian noise of mean 0.0 and σ = 40. Images Measure TV[19] NS[27] Luo[26] RD[25] NLM[39] PV[28] SYS[38] Proposed Boat ISNR PSNRgrad ISNR PSNRgrad ISNR PSNRgrad ISNR PSNRgrad ISNR PSNRgrad ISNR PSNRgrad ISNR PSNRgrad ISNR PSNRgrad 7.02 27.04 7.05 27.23 7.51 27.68 7.12 26.59 7.92 28.11 6.59 26.55 9.11 30.05 12.57 37.34 7.01 27.01 7.12 27.30 7.02 27.33 7.19 26.70 7.83 27.57 6.88 26.65 10.95 32.36 11.64 35.68 6.93 26.86 7.00 27.12 7.22 27.51 7.31 26.83 7.64 27.41 6.62 26.59 9.29 30.72 11.51 35.48 7.78 28.07 7.71 28.05 8.16 28.62 7.58 27.44 8.66 29.14 7.07 27.14 9.97 31.46 12.60 37.37 9.90 30.84 9.12 30.11 9.03 29.84 7.77 27.91 9.96 30.59 7.30 27.68 13.11 34.77 12.74 37.52 9.99 30.90 9.19 30.19 9.20 30.04 8.51 28.48 10.14 30.62 7.55 27.71 13.94 37.37 12.17 36.01 10.19 31.28 9.43 30.21 9.36 30.25 8.66 28.58 9.99 30.60 7.59 27.73 13.87 37.30 11.29 34.92 10.80 31.64 10.00 30.62 9.77 30.62 8.86 28.88 10.66 31.26 7.78 27.98 15.01 38.80 12.86 37.56 Lake Livingroom Mandril Pirate Walkbridge Woman Mosaic values of both MSSIM and PSNR measures are highlighted in each table. The highest values of PSNR and MSSIM for each image, clearly shows that the proposed model is better than all other models considered in the table. We have a significant observation that proposed model performs better for higher noise densities, as can be seen from tables[3]-[4]. For low noise densities, such as σ = 20, our results are still batter or very closed to the best values. Hence, analysis of the results presented in the tables, reveals that the performance of the present model increases with increasing noise density. In table [5], apart from the MSSIM and PSNR, we show the quantitative comparison with gradient PSNR and ISNR values of the proposed model and the alternative approaches. Finally, we note that the different quantitative measures (in terms of PSNR, MSSIM, gradient PSNR, and ISNR values) confirm the superiority of the proposed approach in comparison with other models. Our comparison study shows that the proposed CPDE model preserves image structure efficiently in comparison to other earlier reported methods when signal-to-noise ratio is low. On the other hand, when signal-to-noise ratio is high, the present model restores image structure better than the other existing models. From the above discussion, it is confirmed that the proposed model is robust and more efficient than the available techniques considered here. 6 Conclusion In this work, we propose a new CPDE based denoising framework. The proposed space-time regularization based image denoising approach highlights the choice of diffusion function and data fidelity term. For this purpose, we use the evolution equations to determine both the terms in the proposed model. It may be treated as the non-linear diffusion model which yield two separate PDEs to remove the noise with preservation of significant edges and fine structures. First, we establish the well-posedness of the proposed model using a time discretization method. Then, to solve the proposed model numerically, an implicit finite difference scheme along with advanced iterative solver has been used. We validate our model with different standard test images. Experimental studies confirm that the proposed model is more efficient than the existing models, in terms of quantitative measures and visual quality. Especially, the present model produces best results for high noise levels. We believe that present model has multidimensional applications which will ultimately benefit the society. Competing interests The authors declare that they have no competing interests. Acknowledgments Authors sincerely acknowledge Dr. Arnav Bhavsar, School of Computing and Electrical Engineering, Indian Institute of Technology Mandi, Himachal Pradesh, India for some valuable discussions and suggestions during the preparation of this manuscript. A PREPRINT - AUGUST 8, 2019 References [1] L. Alvarez, F. Guichard, P.-L. Lions, J.-M. Morel, Axioms and fundamental equations of image processing, Archive for rational mechanics and analysis 123 (3) (1993) 199-257. [2] R. C. Gonzalez, R. E. Woods, Digital image processing (2002). [3] G. Aubert, P. Kornprobst, Mathematical problems in image processing: partial differential equations and the calculus of variations, Vol. 147, Springer, 2006. [4] K. Mikula, Image processing with partial differential equations, in: Modern methods in scientifc computing and applications, Springer, 2002, pp. 283-321. [5] J. Nolen, Partial differential equations and difusion processes, Tech. rep., Technical report, Stanford University. Department of Mathematics (2009). [6] J. Weickert, Anisotropic diffusion in image processing, Vol. 1, Teubner Stuttgart, 1998. [7] J. Weickert, Applications of nonlinear diffusion in image processing and computer vision, Acta Math. Univ. Comenianae 70 (1) (2001) 33-50. [8] A. P. Witkin, Scale-space filtering: A new approach to multi-scale description,in: Acoustics, Speech, and Signal Processing, IEEE International Conference on ICASSP’84., Vol. 9, IEEE, 1984, pp. 150-153. [9] F. Catté, P.-L. Lions, J.-M. Morel, T. Coll, Image selective smoothing and edge detection by nonlinear diffusion, SIAM Journal on Numerical analysis 29 (1) (1992) 182-193. [10] P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, Pattern Analysis and Machine Intelligence, IEEE Transactions on 12 (7) (1990) 629-639. [11] J. Weickert, A review of nonlinear diffusion filtering, in: Scale-space theory in computer vision, Springer, 1997, pp. 1-28. [12] X. Liu, L. Huang, Z. Guo, Adaptive fourth-order partial differential equation filter for image denoising, Applied Mathematics Letters 24 (8) (2011) 1282-1288. [13] A. Siddig, Z. Guo, Z. Zhou, B. Wu, An image denoising model based on a fourth-order nonlinear partial differential equation, Computers & Mathematics with Applications 76 (5) (2018) 1056-1074. [14] Y.-L. You, M. Kaveh, Fourth-order partial differential equations for noise removal, IEEE Transactions on Image Processing 9 (10) (2000) 1723-1730. [15] A. Araújo, S. Barbeiro, P. Serranho, Stability of finite difference schemes for complex diffusion processes, SIAM Journal on Numerical Analysis 50 (3) (2012) 1284-1296. [16] A. Araújo, S. Barbeiro, P. Serranho, Stability of finite difference schemes for nonlinear complex reactiondiffusion processes, IMA Journal of Numerical Analysis 35 (3) (2015) 1381-1401. [17] A. Chambolle, P.-L. Lions, Image recovery via total variation minimization and related problems, Numerische Mathematik 76 (2) (1997) 167-188. [18] C. Elliott, S. Smitheman, Numerical analysis of the tv regularization and H −1 fidelity model for decomposing an image into cartoon plus texture, IMA Journal of Numerical Analysis 29 (3) (2009) 651-689. [19] L. I. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D: Nonlinear Phenomena 60 (1) (1992) 259-268. [20] V. S. Prasath, A. Singh, A hybrid convex variational model for image restoration, Applied Mathematics and Computation 215 (10) (2010) 3655-3664. [21] Y.-H. R. Tsai, S. Osher, Total variation and level set methods in image science, Acta Numerica 14 (2005) 509573. [22] R. Zanella, F. Porta, V. Ruggiero, M. Zanetti, Serial and parallel approaches for image segmentation by numerical minimization of a second-order functional, Applied Mathematics and Computation 318 (2018) 153-175. [23] H. Amann, Time-delayed perona-malik type problems, Acta Math. Univ. Comenianae 76 (1) (2007) 15-38. [24] A. Belahmidi, A. Chambolle, Time-delay regularization of anisotropic diffusion and image processing, ESAIM: Mathematical Modelling and Numerical Analysis 39 (2) (2005) 231-251. [25] Z. Guo, J. Yin, Q. Liu, On a reaction-diffusion system applied to image decomposition and restoration, Mathematical and Computer Modelling 53 (5) (2011) 1336-1350. [26] H. Luo, L. Zhu, H. Ding, Coupled anisotropic diffusion for image selective smoothing, Signal Processing 86 (7) (2006) 1728-1736. 17 A PREPRINT - AUGUST 8, 2019 [27] M. Nitzberg, T. Shiota, Nonlinear image filtering with edge and corner enhancement, IEEE transactions on pattern analysis and machine intelligence 14 (8) (1992) 826-833. [28] V. S. Prasath, D. Vorotnikov, On a system of adaptive coupled pdes for image restoration, Journal of mathematical imaging and vision 48 (1) (2014) 35-52. [29] P. Guidotti, Anisotropic diffusions of image processing from perona-malik on, Advanced Studies in Pure Mathematics 99 (2015) 20XX. [30] S. Osher, A. Solé, L. Vese, Image decomposition and restoration using total variation minimization and the H −1 norm, Multiscale Modeling and Simulation 1 (3) (2003) 349-370. [31] R. Adam, Sobolev spaces, in: Pure and Applied Mathematics Series of Monographs and Textbooks, Vol. 65, Academic Press, Inc., New York, San Francisco, London„ 1975. [32] L. Evans, Partial Differential Equations, in: Graduate Studies in Mathematics, Vol. 19, American Mathematical Society, Providence, Rhode Island, 1998. [33] S. K. Jain, R. K. Ray, A. Bhavsar, Iterative solvers for image denoising with diffusion models: A comparative study, Elsevier, 2015. [34] J. L. Lions, Contrôle Optimal de Systèmes Gouvernés par des Équations aux Dérivées Partielles, Dunod, Paris, 1968. [35] M. Chipot, I. Shafrir, V. Valente, G. V. Caffarelli, On a hyperbolic- parabolic system arising in magnetoelasticity, Journal of Mathematical Analysis and Applications 352 (1) (2009) 120-131. [36] S. Zheng, Nonlinear parabolic equations and hyperbolic-parabolic coupled systems, CRC Press, 1995. [37] S. K. Jain, R. K. Ray, A. Bhavsar, A comparative study of iterative solvers for image denoising, in: Proceedings of the 3rd International Conference on Frontiers of Intelligent Computing: Theory and Applications (FICTA) 2014, Springer, 2015, pp. 307-314. [38] J. Sun, J. Yang, L. Sun, A class of hyperbolic-parabolic coupled systems applied to image restoration, Boundary Value Problems 2016 (1) (2016) 187. [39] A. Buades, B. Coll, J.-M. Morel, A non-local algorithm for image denoising, in: Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on, Vol. 2, IEEE, 2005, pp. 60-65. [40] Z. Wang, A. C. Bovik, H. R. Sheikh, E. P. Simoncelli, Image quality assessment: from error visibility to structural similarity, Image Processing, IEEE Transactions on 13 (4) (2004) 600-612. [41] X.-C. Tai, K.-A. Lie, T. F. Chan, S. Osher, Image Processing Based on Partial Differential Equations: Proceedings of the International Conference on PDE-Based Image Processing and Related Inverse Problems, CMA,Oslo, August 8-12, 2005, Springer Science & Business Media, 2006. 18