118 American Academic Scientific Research Journal for Engineering, Technology, and Sciences ISSN (Print) 2313-4410, ISSN (Online) 2313-4402 http://asrjetsjournal.org/ A PDE-based Mathematical Method in Image Processing: Digital-Discrete Method for Perona-Malik Equation Ahmet Yıldırıma*, İsmet Karacab a,bEge University, Faculty of Science, Department of Mathematics, 35100, Bornova, İzmir, Turkey aEmail: yahmet49ege@gmail.com bEmail: ismet.karaca@ege.edu.tr Abstract In this study, we propose a new and effective algorithm for image processing. The method based on the combination of digital topology, partial differential equations and finite difference scheme is called the digital- discrete method. We try to solve the Perona-Malik equation using the digital-discrete method. We use the MATLAB package program when analyzing images. The analyzes we make on the images show how the algorithm is useful, effective and open to development. Keywords: Perona-Malik equation; digital-discrete method; digital topology; finite difference method; image processing. 1. Introduction Image processing problem is one of the most interested subjects in the modern world. There have been great advances in imaging processes, especially in the last 20 years. The widespread use of mobile phones and computers, the continuous development of TV screens, and advances in biomedical imaging have increased the interest in imaging processes. The main aim of our study is to contribute to this development by using mathematical tools. Our work is an interdisciplinary study. In addition to mathematical structures such as Digital Topology, Partial Differential Equations, and finite differences, image processing systems and MATLAB package program formed the basis of our work. ------------------------------------------------------------------------ * Corresponding author. http://asrjetsjournal.org/ American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 84, No1, pp 118-129 119 Especially in the last 20 years, the use of PDEs in image processing has started to increase. Mikula and Ramarosy implemented semi-implicit finite volume scheme in image processing [1]. Sarti and his colleagues used evolutionary PDEs in medical imaging [2]. Weickert applied additive operator splitting method to nonlinear diffusion equations [3]. Mikula performed image analysis using different partial differential equations [4]. Gibou and his colleagues suggested automatic algorithms for the segmentation phase of radiotherapy treatment planning [5]. Angenent and his colleagues used PDEs in medical imaging [6]. Kuijper used geometric and variational PDEs for image processing [7]. Nadernejad and his colleagues implemented PDEs-based model in image enhancement [8]. Kim and Lim applied fourth order PDEs-based methods in image denoising [9]. Lin and his colleagues designed PDEs and differential invariants for image problems [10]. Belaid used topological gradient approach in image processing [11]. Niang and his colleagues used PDEs-based method for image problems [12]. Shen and his colleagues implemented an adaptive partial differential equation in image processing [13]. Ghanbari and his colleagues used a restarted iterative HAM for TV models [14]. Liu and his colleagues applied an adaptive relaxation method for a fourth order PDE for image problems [15]. Nnolim proposed partial differential equation-based method for underwater image processing [16]. Yu and his colleagues implemented quasi-interpolation operators for bivariate quintic spline spaces in image processing [17]. Huang and his colleagues implemented hybrid analog-digital approach for image problems [18]. Benseghir and his colleagues used a partial differential equation that is based on a nonlinear structure tensor in image processing [19]. In this paper, we will use digital-discrete approach for Perona-Malik equation in image processing. 2. Digital Continuous Function (DCF) and Gradually Varied Function (GVF) 2.1. Digital Function (DF) [20] DCFs and GVFs were developed in 1980s. Rosenfeld suggested digital DCFs in digital imaging [21]. Chen proposed GVFs for interpolating a digital surface [22]. Khalimsky [23], Kong [24], Boxer [25], Rosenfeld [26] developed methods for digital deformations . Agnarsson and Chen established association between graph homomorphism and GVFs and [27]. Chen found systematic digital-discrete method [28] and Chen and Luo suggested harmonic functions in image analysis [29]. The DF is essentially a function defined in digital spaces. We can see that digital space is a subspace of Euclidean space. Digital space contains all grid integer points. We generally represent m-dimensional digital space with ∑𝑚. Definition 2.1.1. [20,21]: A DF is a function from digital space I defined to integers. Usually we use {1,2,...,n} as the display set of the function. For example, 𝑓: ∑2 → {1,2, . . . , 𝑛} is a digital function. Also we can see the DF as a function in general discrete space. Definition 2.1.2. [20,21]: A DCF can be defined as a function with integer values of digital points are the same as its neighbor or differ by at most 1. American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 84, No1, pp 118-129 120 Let 𝑓: ∑𝑘 → {… ,1,2,3, … } be a digital function and x, y be two close points in ∑2. If | 𝑓(𝑥) − 𝑓(𝑦)| ≤ 1, then we say that f is a DCF. Proposition 2.1.1. [20,21]: Any DCF is a Lipschitz function. Because | 𝑓(𝑥) − 𝑓(𝑦)| ≤ 𝑑(𝑥, 𝑦) is satisfied. 2.2. Gradually Varied Function (GVF) [20] GVF was proposed by Chen in discrete and digital spaces [20,22,30]. 𝐴1, 𝐴2,….,𝐴𝑚 ∈ ℝ, 𝐴1 < 𝐴2 < ⋯ < 𝐴𝑚 and 𝑝, 𝑞 ∈ ∑2 are given. Let 𝑓: ∑2(𝑜𝑟 𝑎𝑛𝑜𝑡ℎ𝑒𝑟 𝑑𝑖𝑠𝑐𝑟𝑒𝑡𝑒 𝑠𝑝𝑎𝑐𝑒) → {𝐴1, 𝐴2,….,𝐴𝑚} be a function. If 𝑓(𝑝) = 𝐴𝑖 and 𝑓(𝑞) = 𝐴𝑗 , then level difference is | i-j | between f(p) and f(q). The GVF can be defined as follows. Definition 2.2.1. [20,22,30]: For any p, q close pair in ∑2, if 𝑓(𝑝) = 𝐴𝑖 and 𝑓(𝑞) = 𝐴𝑖−1, 𝐴𝑖 𝑜𝑟 𝐴𝑖+1 , then f is GVF on p and q. Definition 2.2.2. [20,22,30 ]: If f is GVF on any p, q close pair in ∑2, then f is GVF. Theorem 2.2.1. [20,22,30]: Let p, q be points in J, if a gradually varied interpolation exists, then length of the shortest path in D between p and q is not less than level difference between f(p) and f(q). Gradually Varied Function (GVF) Algorithm [20,31]: J is a non-empty subset of D. fJ is a function defined on J. Step 1: We test all p and p/ points in J. If d(p, p/ ) ≥ LD(p, p/) is not satisfied, then there is no GVF. Take D0 ← J. Step 2: We take x from D − D0 where x has an adjacent vertex r in D0. Suppose fD (r) = Ai. Step 3: We take fD(x) = fD(r) = Ai. We test x against every vertex p in D0: If there is a p∈ D0 when d(x, p) < LD(x, p), change fD(x) to Ai−1 when fD(p) < Ai or change fD(x) to Ai+1 when fD(p) > Ai. Step 4: Let D0 ← D0 ∪ {x}. Step 5: We repeat 2–4 until D0 = D. American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 84, No1, pp 118-129 121 3. Digital-Discrete method for Perona-Malik equation 3.1. Perona-Malik equation [32,33] The main purpose of nonlinear diffusion models is to create a measure for the preservation and enhancement of edges with images. The most important tools for this are nonlinear partial differential equations. The first nonlinear diffusion model used in image processing is called anisotropic diffusion. The anisotropic diffusion model was first proposed by Perona and Malik [34]. The Perona-Malik equation is as follows [32,33,34]: 𝜕𝑢 𝜕𝑡 = 𝑑𝑖𝑣(𝑔‖∇𝑢‖∇𝑢), 𝑡 ≥ 0, 𝑥 ∈ Ω 𝑢0(𝑥) = 𝑓(𝑥) (3.1.1) The Perona-Malik equation has been used for multiple scaling, enhancement and splitting of images. In equation (3.1.1), div is the divergence operator and u is the smooth image at time t. ‖∇𝑢‖ is the gradient length of u and 𝑔‖∇𝑢‖ is the diffusivity function. f(x) is initial image. The diffusivity function g is non-negative and monotonically decreasing. It has following properties: 𝑔(0) = 1, 𝑔(𝑠) ≥ 0 𝑎𝑛𝑑 lim 𝑠→∞ 𝑔(𝑠) = 0. Perona and Malik proposed two different choices for the diffusivity function [32,33,34]: 𝑔(𝑠) = 1 1+ 𝑠2 𝜆2 (3.1.2) 𝑔(𝑠) = 𝑒 − 𝑠2 𝜆2 (3.1.3) 3.2. A numerical approach for the Perona-Malik equation with the finite difference method [32,33,35] Let's make the Perona-Malik equation given in (3.1.1) discrete. We can divide the ranges in the form 𝑥𝑖 = 𝑖ℎ1, 𝑖 = 1,2,3, … , 𝑁 𝑦𝑗 = 𝑗ℎ2, 𝑗 = 1,2,3, … , 𝑀 (3.2.1) 𝑡𝑘 = 𝑘. ∆𝑡, 𝑘 = 1,2,3, … , 𝑛 Where, h1 and h2 indicate pixel placements in the x and y directions, respectively. Pixels are generally American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 84, No1, pp 118-129 122 considered to have a unit length, h1=h2 =1. We obtain a numerical approximation for the gradient length ‖∇𝑢𝑖,𝑗‖ by using the central finite difference method and we use this approximation to calculate the diffusivity function 𝑔‖∇𝑢𝑖,𝑗‖ [32,33]: ‖∇𝑢𝑖,𝑗‖ = √( 𝑑𝑢𝑖,𝑗 𝑑𝑥 ) 2 + ( 𝑑𝑢𝑖,𝑗 𝑑𝑦 ) 2 ≈ √( 𝑢𝑖+1,𝑗 − 𝑢𝑖−1,𝑗 2ℎ1 ) 2 + ( 𝑢𝑖,𝑗+1 − 𝑢𝑖,𝑗−1 2ℎ2 ) 2 (3.2.2) The left side of equation (3.1.1) can be discretized as [32,33]: 𝜕𝑢 𝜕𝑡 = 𝑢𝑖,𝑗 𝑘+1−𝑢𝑖,𝑗 𝑘 ∆𝑡 (3.2.3) The right side of equation (3.1.1) can be discretized as [32,33]: 𝜕 𝜕𝑥 (𝑔‖∇𝑢‖ )𝑢𝑥 + 𝜕 𝜕𝑦 (𝑔‖∇𝑢‖ )𝑢𝑦 = 𝑔 𝑖+ 1 2,𝑗 (𝑢𝑖+1, 𝑗 𝑘 − 𝑢𝑖, 𝑗 𝑘 )−𝑔 𝑖− 1 2,𝑗 (𝑢𝑖,𝑗 𝑘 − 𝑢𝑖−1,𝑗 𝑘 ) ℎ1 2 + 𝑔 𝑖,𝑗+ 1 2 (𝑢𝑖,𝑗+1 𝑘 − 𝑢𝑖,𝑗 𝑘 )−𝑔 𝑖,𝑗− 1 2 (𝑢𝑖,𝑗 𝑘 − 𝑢𝑖,𝑗−1 𝑘 ) ℎ2 2 (3.2.4) The lattice we use for discretization in equations (3.2.2), (3.2.3),(3.2.4) is as follows [32,33]: 𝑢𝑖−1,𝑗−1 𝑢𝑖,𝑗−1 𝑢𝑖+1,𝑗−1 𝑔 𝑖,𝑗− 1 2 𝑢𝑖−1,𝑗 𝑔 𝑖− 1 2 ,𝑗 𝑢𝑖,𝑗 𝑔 𝑖+ 1 2 ,𝑗 𝑢𝑖+1,𝑗 𝑔 𝑖,𝑗+ 1 2 𝑢𝑖−1,𝑗+1 𝑢𝑖,𝑗+1 𝑢𝑖+1,𝑗+1 Figure 3.2.1 In the discretized state, mid-pixel points need calculated diffusivities. Therefore, we can simply calculate these values by taking advantage of the diffusivity in neighboring pixels. We can do this as follows [32,33]: American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 84, No1, pp 118-129 123 𝑔 𝑖+ 1 2 ,𝑗 = 𝑔𝑖+1,𝑗 + 𝑔𝑖,𝑗 2 𝑔 𝑖− 1 2 ,𝑗 = 𝑔𝑖−1,𝑗 + 𝑔𝑖,𝑗 2 𝑔 𝑖,𝑗+ 1 2 = 𝑔𝑖,𝑗+1 + 𝑔𝑖,𝑗 2 𝑔 𝑖,𝑗− 1 2 = 𝑔𝑖,𝑗−1 + 𝑔𝑖,𝑗 2 (3.2.5) As a result of the operations, the discretized version of the Perona-Malik equation is as follows [32,33]: 𝑢𝑖,𝑗 𝑘+1 − 𝑢𝑖,𝑗 𝑘 ∆𝑡 = (𝑔𝑖+1,𝑗 𝑘 + 𝑔𝑖,𝑗 𝑘 )[𝑢𝑖+1,𝑗 𝑘 − 𝑢𝑖,𝑗 𝑘 ] − (𝑔𝑖,𝑗 𝑘 + 𝑔𝑖−1,𝑗 𝑘 )[𝑢𝑖,𝑗 𝑘 − 𝑢𝑖−1,𝑗 𝑘 ] +(𝑔𝑖,𝑗+1 𝑘 + 𝑔𝑖,𝑗 𝑘 )[𝑢𝑖,𝑗+1 𝑘 − 𝑢𝑖,𝑗 𝑘 ] − (𝑔𝑖,𝑗 𝑘 + 𝑔𝑖,𝑗−1 𝑘 )[𝑢𝑖,𝑗 𝑘 − 𝑢𝑖,𝑗−1 𝑘 ] (3.2.6) ‖∇𝑢𝑖,𝑗‖ = √( 𝑢𝑖+1,𝑗 − 𝑢𝑖−1,𝑗 2 ) 2 + ( 𝑢𝑖,𝑗+1 − 𝑢𝑖,𝑗−1 2 ) 2 (3.2.7) 3.3. Image processing analysis by using digital-discrete method and Perona-Malik equation We discretized the Perona-Malik equation using the finite difference method in Section 3.2. We do digital- discrete adaptation in (k+1) so that 𝑢𝑖,𝑗 𝑘+1 ← ( 𝑢𝑖,𝑗 𝑘+1 + 𝑔𝑣𝑓(𝑖, 𝑗)/2) in this chapter. We continue to adapt in this way. While creating our algorithm, we take advantage of the gradually varied function structure [20]. We apply the finite difference scheme while constructing the numerical structure of the Perona-Malik equation [32,33,35]. Our algorithm, which we have established with the help of the Perona-Malik equation and the digital-discrete method, is given below [20,32,33,35]: Algorithm (Digital-Discrete method for Perona-Malik equation in image processing) Step 1: Upload the main points. Data points with observation values are loaded. Step 2: Identify the solution. Arrange the points in the lattice space. Step 3: Expand the function according to Theorem 2.2.1 using the local Lipschitz condition. GVFs are obtained. (We use Gradually Varied Function (GVF) Algorithm) Step 4: Start getting the lattice points with the iteration of the Perona-Malik equation that we obtained by the finite difference method: 𝑢𝑖,𝑗 𝑘+1 − 𝑢𝑖,𝑗 𝑘 ∆𝑡 = (𝑔𝑖+1,𝑗 𝑘 + 𝑔𝑖,𝑗 𝑘 )[𝑢𝑖+1,𝑗 𝑘 − 𝑢𝑖,𝑗 𝑘 ] − (𝑔𝑖,𝑗 𝑘 + 𝑔𝑖−1,𝑗 𝑘 )[𝑢𝑖,𝑗 𝑘 − 𝑢𝑖−1,𝑗 𝑘 ] American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 84, No1, pp 118-129 124 +(𝑔𝑖,𝑗+1 𝑘 + 𝑔𝑖,𝑗 𝑘 )[𝑢𝑖,𝑗+1 𝑘 − 𝑢𝑖,𝑗 𝑘 ] − (𝑔𝑖,𝑗 𝑘 + 𝑔𝑖,𝑗−1 𝑘 )[𝑢𝑖,𝑗 𝑘 − 𝑢𝑖,𝑗−1 𝑘 ] ‖∇𝑢𝑖,𝑗‖ = √( 𝑢𝑖+1,𝑗 − 𝑢𝑖−1,𝑗 2 ) 2 + ( 𝑢𝑖,𝑗+1 − 𝑢𝑖,𝑗−1 2 ) 2 Step 5: Restore the current values obtained in step 4 using a gradually varied function: We continue to use digital-discrete adaptation in (k+1) so that 𝑢𝑖,𝑗 𝑘+1 ← (𝑢𝑖,𝑗 𝑘+1 + 𝑔𝑣𝑓(𝑖, 𝑗)/2) Step 6: Perform image processing analysis with the help of the MATLAB-R2020b package program using our algorithm. We obtained images for different t values in Figure 3.3.1 and Figure 3.3.2 by using our algorithm. Original Image American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 84, No1, pp 118-129 125 Figure 3.3.1 Original Image American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 84, No1, pp 118-129 126 Figure 3.3.2 American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 84, No1, pp 118-129 127 4. Conclusion In this study, we applied a new and effective algorithm for image processing. We used the MATLAB-R2020b package program when analyzing images. The greatest progress we have achieved in image analysis is the preservation of image quality and the preservation of images as the t values increase. The solution algorithm we use is simpler than other algorithms and it is easier to apply to package programs such as MATLAB. Since the algorithm is based on topology, graph theory and continuity, the values we obtained in the revised finite difference approach with digital topology are efficient and stable. We are considering using the digital-discrete method and the methods produced from this method in imaging systems in other fields such as underwater imaging systems and biomedical imaging systems. Acknowledgement This study was produced from the first author's doctoral thesis. References [1]. K.Mikula, N.Ramarosy, Semi-implicit finite volume scheme for solving nonlinear diffusion equations in image processing, Numerische Mathematik 89 (2001) 561–590 [2]. A.Sarti, K.Mikula, F.Sgallari, C. Lamberti, Evolutionary partial differential equations for biomedical image processing, Journal of Biomedical Informatics 35 (2002) 77–91 [3]. J.Weickert, Efficient image segmentation using partial differential equations and morphology, Pattern Recognition 34 (2001) 1813-1824 [4]. K.Mikula, Image processing with partial differential equations, Modern Methods in Scientific Computing and Applications 75 (2002) 283-321 [5]. F.Gibou, Partial differential equations-based segmentation for radiotherapy treatment planning, Mathematical Biosciences and Engineering 2(2) (2005) 209-226 [6]. S.Angenent, E.Pichon, A.Tannenbaum, Mathematical Methods in Medical Image Processing, Bulletin (New Series) of the American Mathematical Society 43(3) (2006) 365–396 [7]. A.Kuijper, Image Processing with Geometrical and Variational PDEs, 31st annual workshop of the Austrian Association for Pattern Recognition (OAGM/AAPR), OEAGM07 (Schloss Krumbach, Austria, May 3-4, (2007) 89-96 [8]. E.Nadernejad , H.Koohi , H.Hassanpour, PDEs-Based Method for Image Enhancement, Applied Mathematical Sciences, 2(20) (2008) 981-993 [9]. S.Kim. H.Lim, Fourth-order partial differential equations for effective image denoising, Seventh American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 84, No1, pp 118-129 128 Mississippi State - UAB Conference on Differential Equations and Computational Simulations, Electronic Journal of Differential Equations, Conf. 17 (2009) 107–121 [10]. Z.Lin, W.Zhang, X.Tang, Designing Partial Differential Equations for Image Processing by Combining Differential Invariants, Microsoft Research Newsletter, MSR-TR-2009-192 [11]. L.J.Belaid, An Overview of the Topological Gradient Approach in Image Processing: Advantages and Inconveniences, Journal of Applied Mathematics Volume 2010, Article ID 761783, 19 pages [12]. O.Niang, A.Thioune, M.C.E.Gueirea, E.Delechelle, J.Lemoine, Partial Differential Equation-Based Approach for Empirical Mode Decomposition: Application on Image Analysis, IEEE Transactions on Image Processing 21(9) (2012) 3991-4001 [13]. C.Shen, C. Li, P.Wang, An Adaptive Partial Differential Equation for Noise Removal, Proceedings of the 2nd International Symposium on Computer, Communication, Control and Automation 4 (2013) 55- 58 [14]. B.Ghanbari, L.Rada, K.Chen, A restarted iterative homotopy analysis method for two nonlinear models from image processing, International Journal of Computer Mathematics 91(3) (2014) 661–687 [15]. X.Y. Liu, C.H. Lai, K.A. Pericleous, A fourth-order partial differential equation denoising model with an adaptive relaxation method, International Journal of Computer Mathematics, 92(3) (2015) 608-622 [16]. U.A.Nnolim, Analysis of proposed PDE-based underwater image enhancement algorithms, arXiv:1612.04447v1 [cs.CV] [17]. R.Yu, C.Zhu, X.Hou, L.Yin, Quasi-Interpolation Operators for Bivariate Quintic Spline Spaces and Their Applications, Mathematical and Computational Applications 22(1) (2017) 10: 15 pages [18]. Y.Huang, N.Guo, M.Seok, Y.Tsividis, K.Mandli, S.Sethumadhavan, Hybrid Analog-Digital Solution of Nonlinear Partial Differential Equations, In Proceedings of MICRO-50, Cambridge, MA, USA, October 14–18 (2017) 14 pages [19]. M.Benseghir, F.Z.Nouri, P.C.Tauber, A New Partial Differential Equation for Image Inpainting, Boletim da Sociedade Paranaense de Matemática 39(3) (2021) 137-155 [20]. L.M. Chen, Digital Functions and Data Reconstruction, Digital-Discrete Methods, Springer-Verlag New York (2013) [21]. A.Rosenfeld, ‘Continuous’ functions on digital pictures, Pattern Recognition Letters 4(3) (1986) 177- 184 [22]. L. Chen, The necessary and sufficient condition and the efficient algorithms for gradually varied fill, American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 84, No1, pp 118-129 129 Chinese Science Bulletin 35 (10) (1990) 870–873 [23]. E.Khalimsky, Motion, deformation, and homotopy in finite spaces, In Proceedings IEEE international conference on systems, man, and cybernetics, Chicago (1987) 227–234 [24]. T.Y.Kong, A digital fundamental group, Computers and Graphics 13 (1989) 159–166 [25]. L.Boxer, Digitally continuous functions, Pattern Recognition Letters 15(8) (1994) 833–839 [26]. A.Rosenfeld, , Contraction of digital curves, University of Maryland’s Technical Report (1996) [27]. G.Agnarsson, L.Chen, On the extension of vertex maps to graph homomorphisms, Discrete Mathematics 306(17) (2006) 2021–2030 [28]. L.Chen, A digital-discrete method for smooth-continuous data reconstruction, Journal of the Washington Academy of Sciences 96(2) (2010) 47–65 [29]. L.Chen, F.Luo, Harmonic Functions for Data Reconstruction on 3D Manifolds, arXiv:1102.0200v3 [math.NA] [30]. L.Chen, The properties and the algorithms for gradually varied fill, Chinese Journal of Computers 14(3) (1991) [31]. T.H.Cormen, C.E.Leiserson, R.L.Rivest, Introduction to algorithms, MIT, New York (1993) [32]. E.Erdem, Nonlinear Diffusion, Hacettepe University, Department of Computer Engineering Lecture Notes (2013) 1-15 [33]. O.Demirkaya, M.H. Asyali, P.K. Sahoo, Image Processing with MATLAB: Applications in Medicine and Biology, CRC Press: Taylor& Francis Group, Boca Raton, USA (2009) [34]. P.Perona, J.Malik, Scale-Space and Edge Detection Using Anisotropic Diffusion, IEEE Transactions on Pattern Analysis and Machine Intelligence 12(7) (1990) 629-639 [35]. G. D. Smith, Numerical Solution of Partial Differential Equations: Finite Difference Methods, Oxford Applied Mathematics and Computing Science Series (1985)