versión On-line ISSN 1807-0302
Comput. Appl. Math. vol.30 no.3 São Carlos 2011
Mimetic finite difference methods in image processing
C. BazanI; M. AboualiI; J. CastilloI; P. BlomgrenII
IComputational Science Research Center, San Diego State University. 5500 Campanile Drive, San Diego, CA 92182-1245, U.S.A. E-mails: firstname.lastname@example.org / email@example.com / firstname.lastname@example.org
IIDepartment of Mathematics and Statistics, San Diego State University. 5500 Campanile Drive, San Diego, CA 92182-7720, U.S.A. E-mail: email@example.com
We introduce the use of mimetic methods to the imaging community, for the solution of the initial-value problems ubiquitous in the machine vision and image processing and analysis fields. PDE-based image processing and analysis techniques comprise a host of applications such as noise removal and restoration, deblurring and enhancement, segmentation, edge detection, inpainting, registration, motion analysis, etc. Because of their favorable stability and efficiency properties, semi-implicit finite difference and finite element schemes have been the methods of choice (in that order of preference). We propose a new approach for the numerical solution of these problems based on mimetic methods. The mimetic discretization scheme preserves the continuum properties of the mathematical operators often encountered in the image processing and analysis equations. This is the main contributing factor to the improved performance of the mimetic method approach, as compared to both of the aforementioned popular numerical solution techniques. To assess the performance of the proposed approach, we employ the Catté-Lions-Morel-Coll model to restore noisy images, by solving the PDE with the three numerical solution schemes. For all of the benchmark images employed in our experiments, and for every level of noise applied, we observe that the best image restored by using the mimetic method is closer to the noise-free image than the best images restored by the other two methods tested. These results motivate further studies of the application of the mimetic methods to other imaging problems.
Mathematical subject classification: Primary: 68U10; Secondary: 65L12.
Key words: mimetic methods, image processing, discrete operators, conservative methods.
The aim of this paper is to introduce the mimetic methods to the imaging community, for the solution of the initial-value problem ubiquitous in the machine vision and image processing and analysis fields. PDE-based image processing and analysis techniques comprise a host of applications such as noise removal and restoration, deblurring and enhancement, segmentation, edge detection, inpainting, registration, motion analysis, etc. In this context, a gray-scale image is modeled as a real-valued function u0(x), u0 → , defined in a bounded domain Ω ⊂ 2, and with Lipschitz continuous boundary ∂ Ω. Usually, the observed image u0(x) = u(x, 0) is associated with a sequence of images u(x, t), where the evolution depends on the abstract parameter t > 0, called the scale. Hence the name image multiscale analysis given to this approach . The numerical solution of this problem is normally based on semi-discretization in scale and on finite difference or finite element discretization in space.
This paper is organized as follows: In Section 2, we briefly describe one of the most popular nonlinear diffusion models applied in image processing for the reduction of noise and the detection of edges. This serves as background for the readers who might not be familiar with PDE-based image processing techniques. In Section 3, we describe the mimetic discretization formulation for the solution of the initial-value problem. We also present the two most popular numerical solutions to this problem, namely finite difference and finite element methods. In Section 4, we present some computational examples of the performance of the proposed method as compared to the other two methods described in the previous section. We conclude the paper in Section 5 with a discussion of the reasons for the improved results obtained by applying the mimetic method. We also outline some possible future improvements to the this approach, and other areas within the imaging field where the method can be applied successfully.
2 Nonlinear diffusion models in image processing
During the last two decades PDE-based models in the field of image processing and analysis have become very popular . Most of today's nonlinear diffusion models stem from the work of Perona and Malik , who introduced an approach that allows for the reduction of the noise in images while retaining and enhancing edges. In order to achieve these goals, the Perona-Malik's nonlinear diffusion model aimed at avoiding the blurring of edges, and other localization problems presented by linear diffusion models [9, 26, 29, 44]. Their model accomplishes this by applying a process that reduces the diffusivity in places having higher likelihood of belonging to edges. This likelihood is measured by a function of the local gradient, | ∇u|. The Perona-Malik model can be written as
where ‹g · ∇u, n › = 0 denotes homogeneous Neumann boundary conditions. In this model the diffusivity has to be such that g(| ∇u|2) → 0 when | ∇u| → ∞ and g(| ∇u|2) → 1 when | ∇u| → 0. One of the diffusivities that Perona and Malik proposed is g(| ∇u|2) = (1 + | ∇u|2 / λ2) 1, where λ > 0 is a threshold (contrast) parameter that separates forward and backward diffusion . The model accomplishes the long sought effect of blurring small fluctuations (possible noise) while enhancing edges. The results obtained by Perona and Malik were visually very impressive.
Notwithstanding the practical success of the Perona-Malik model, it presents some serious theoretical problems: (i) none of the classical well-posedness frameworks is applicable to the Perona-Malik model, i.e. we can not ensure well-posedness results [34, 42]; (ii) uniqueness and stability with respect to the initial image should not be expected, i.e. solvability is a difficult problem, in general [15, 21, 22, 25, 36]; (iii) the regularizing effect of the discretizationplays too much of an important role in the solution [6, 17]. The latter is perhaps the key element in the success or failure of the model. Most practical applications work very well provided that the numerical schemes stabilize the process through some implicit (or explicit) regularization.
This observation motivated much research towards the introduction of the regularization directly into the PDE to avoid the dependence on the numerical schemes [15, 34]. A variety of spatial, spatio-temporal, and temporal regularization procedures have been proposed over the years [5, 15, 28, 38, 40, 43]. The one that has attracted much attention is the mathematically sound formulation due to Catté, Lions, Morel and Coll . They proposed replacing the diffusivity g(| ∇u|2) of the Perona-Malik model by a slight variation g(| ∇u σ|2) with u σ = G σ * u, where G σ is a smooth kernel (Gaussian of variance σ2). We should note that this spatial regularization model belongs to a class of well-posed problems (existence and uniqueness were proven in ), and that its successful implementation is contingent on the choosing of an appropriate value for the additional regularization parameter σ. Whitaker and Pizer  and Li and Chen  suggested making the parameters σ and λ time-dependent, while Benhamouda  performed a systematic study of the influence of these parameters for the one-dimensional case.
3 Numerical solution to the nonlinear diffusion models
Digital images are given on discrete (regular) grids. This lends itself for discretizing the PDEs to obtain numerical schemes that can be solved on a computer. Because of their favorable stability and efficiency properties, semi-implicit schemes have been the methods of choice for the scale discretization [3, 4, 15, 16, 18, 19, 20, 24, 27, 31, 32, 37, 39, 41]. As for the space discretization, the most popular choices are finite difference [15, 39, 41] and finiteelement [3, 4, 16, 24, 37, 39, 41] methods (in that order of preference). We propose a new approach in image processing based on mimetic discretization.
3.1 Finite difference implementation
The numerical solution to the Catté-Lions-Morel-Coll model proposed in  is as follows. Given an N × M image we introduce the coordinates lattice (ih, jh, n Δt) where h is the pixel size, and 0 < i < N + 1, 0 < j < M + 1. We consider as an approximation of u(ih, jh, n Δt), and as an approximation of g(| ∇u σ|)(ih, jh, n Δt). Then we discretise g(| ∇u σ|) ∂u / ∂x ∂x by ∂u / ∂x (ih, jh, (n + 1) Δt) and ∂/ ∂x [g(| ∇u σ|) ∂u / ∂x] by
and similarly for ∂/ ∂y [g(| ∇u σ|) ∂u / ∂y], by exchanging the roles of parameters i and j,
The finite difference scheme will be given by
Then, the discrete problem can be written as a system
where the matrix of coefficients Ah is positive definite and block-tridiagonal.
3.2 Finite element implementation
The starting point for the finite element method is to partition the geometry (domain) into small units (elements or cells) of simple shape joined together at the vertices (nodes). This will constitute our finite element space (mesh or grid). Once we have our mesh (see Fig. 1), the idea is to approximate the dependent variables with functions that we can describe with a finite number of parameters (degrees of freedom, DOF). Inserting this approximation into the weak form of the equation for the Catté-Lions-Morel-Coll model generates a system of equations for the degrees of freedom .
As mentioned above, we need to perform discretizations in scale and space. We perform the semi-discretization in scale by letting Q ∈ , and Δt = T / Q be fixed numbers (here, T represents the last scale state we want to reach), and letting u( x, 0) = u0(x) in Ω. Then, we can look for a function un for every n = 1, ..., Q, such that it is a solution to the equation
It is shown in [4, 24] that there exist unique variational solutions un of Eq. (4) at every discrete scale step, for which the following stability estimates hold:
where C is a general (large) constant (here, h represents a typical element size). To discretize the problem in space we can take advantage of the pixel structure of the image. For this case, the finite element method assumes that the approximation of the solution to the PDE is continuous piecewise linear. This means that the discrete intensity values are regarded as approximations of the continuous intensity function in the center of the pixels (see Fig. 1). We can multiply Eq. (4) by an arbitrary test function v ∈ V, where V is the Sobolev space W1,2( Ω) of L2( Ω)-functions with doubly integrable weak derivatives, and integrate (using Green's theorem and homogeneous Neumann boundary conditions) to obtainthe weak form 
Then, for each scale step n, we look for a continuous piecewise linear function uhn ∈ Vh that satisfies
for all vh ∈ Vh. Considering the standard Lagrangian base functions ϕq ∈ Vh, q = 1, ..., L, given by ϕq(xp) = δqp (Kronecker delta) for all nodes, the function uhn is given by
Substituting Eq. (8) into Eq. (7) and considering as test functions vh = ϕq for q = 1, ..., L, we get the Ritz-Galerkin equation for the nodal values upn, of the piecewise linear function uhn:
Then, in each scale step we need to assemble and solve a linear system of the form
for the vector of unknowns (DOF) un.
3.3 Mimetic discretization implementation
In the mimetic discretization approach, instead of discretizing the actual equation, we use the equivalent discrete version of its mathematical operators. The majority of the PDE-based models in image processing can be written in terms of gradient, divergence, and curl operators. The mimetic method, as its name implies, mimics the properties of these operators in their continuous form and thus preserve their continuum properties exactly, not approximately [2, 8, 7, 11, 13, 14]. In fact, Bohner and Castillo  have stated that "... a discretization of the first derivative and the one-dimensional integral are mimetic if they are analogous of the fundamental theorem and the integration by parts formula that are exact." For this reason, the mimetic discretization method has shown to be more stable and accurate than other numerical discretization methods [11, 12, 8, 23]. Here, we present only the second-order accuracy mimetic gradient and divergence operators developed by Castillo and Grone [10, 14] that were used in our numerical experiments. For a detailed explanation of how these operators were developed we refer the reader to [10, 14] and the references therein.
The Catté-Lions-Morel-Coll variant to Eq. (1) can be written using mimetic operators as
where D is the mimetic divergence operator, ∇·, and G is the mimetic gradient operator, ∇. In one dimension and on a uniform staggered grid, the mimetic gradient operator with second order accuracy is defined as
where h is the grid spacing. Likewise, the mimetic divergence is defined as
The divergence matrix D satisfies some desired properties listed in . These properties are:
D has zero row sums, i.e., De = 0, where e = ( 1, 1, ..., 1)T.
D has column sums -1, 0, ..., 0, 1, i.e., eT D = (-1, 0, 0, ..., 0, 1).
D is banded.
D has a Toeplitz-type structure on the interior rows and is defined independently of the number of grid points.
One of the mimetic method's interesting properties is that the matrix D and G are centro-skew-symmetric. It has to be noted that the gradient and the divergence are not calculated at the same position. The gradient is calculated at the edges of the cell, whereas the divergence is calculated at the center of the cell (see Fig. 2).
In two dimensions, the positions where the gradient is calculated is shown in Fig. 3. Notice that the x-component and y-component of the gradient are calculated separately and in different positions. This is important when designing a parallel code. Similarly, the position where the divergence is calculated isshown in Fig. 4. Notice that the mimetic gradient and mimetic divergence are designed in such way that the output of the mimetic gradient is positionedexactly where the mimetic divergence reads its input.
One of the advantages of the mimetic gradient and mimetic divergence developed by Castillo and Grone is the way in which they represent the boundary conditions. The Neumann boundary conditions imposed in Eq. (1) dictate that the flux through the boundaries is zero. This is equivalent to saying that thenormal component of the gradient to the boundary is zero. By referring to Fig. 3, we can understand that, to impose this boundary condition, all we have to do is set the gradient at the boundary to be zero. This has a great advantage. Usually the boundary conditions and their implementation is a major source of errors in numerical modeling and simulation. To implement the boundary conditions, normally some type of interpolation, extrapolation, or the use of ghost/dummy nodes is needed. Here, in this article, the boundary conditions dictate that the normal component of the gradient at the edges of the image must be zero,Eq. (1). To calculate the normal component of the gradient at the outer edges of the image using the mimetic method, we need the value of u at the center of the two cells adjacent to the edge and also its value on the edge itself, see Fig. 3 and Eq. (12). But the u here refers to the value in the image, i.e., we have the value of u just at the cell centers and not the outer edges. To get the value of u at the outer edges we have to use the extrapolation techniques. The value of u at the outer edges are needed only to calculate the gradient.The value of the gradient at the outer edges is given by the boundary condition exactly at the same position and for the same component that our mimetic method calculates the gradient.Hence, there is no need for any extrapolation or use of dummy nodes to implement the boundary conditions. Therefore, using Castillo Grone mimetic method, one big source of error can be eliminated.
4 Numerical experiments and results
In order to compare the performance of the mimetic discretization we implemented the Perona-Malik variant by Catté, Lions, Morel and Coll,
It has been shown  that σ = 1 is sufficient for a large interval of noise variances, provided that the noise in neighboring pixels is uncorrelated and that the grid size is one. There are several ways to set the parameter λ > 0. Perona and Malik  suggested using the idea presented by Canny  and set λ as a percentile, p, of the image gradient magnitudes at each iteration. (The recommended value is commonly p = 90%.) A by-product of this approach is a decreasing λ, which has an stabilizing effect on the diffusion process . A time step of δt = 10-2 was chosen to update all the models. Weickert, Romeny, and Viergever  have shown that, for explicit discretization schemes, the stability condition, assuming δx = 1 and ∀s : g(s) < 1, is δt < 1 / (2d), with d being the number of dimensions of the data (which for a 2D image d = 2).
The experiment consisted in trying to restore the noise-free image f (x), that has been perturbed by additive Gaussian noise of zero mean and variance 0.001 > ν > 0.02. Figure 1 shows the images of the Cameraman, the Baboon, the Boats, and Barbara. These are the benchmark images that will be used in the comparative experiments. The three solution methods, finite difference, finite element, and mimetic discretization, were run for 50 iterations and the correlation coefficient between the noise-free image and each of the filtered images was computed at each iteration. This measure indicates how similar is the filtered image to the noise-free image after restoration. For every benchmark image and for every level of noise, we observe that the best image restored by the mimetic discretization is closer to the noise-free image than the best images restored by the other two methods tested (see Figs. 2, 4, 6, and 8). The quality of the image restoration after applying the three solution methods to the noisy benchmark images is illustrated in Figs. 3, 5, 7, and 9.
In this paper we introduced the mimetic discretization method for the numerical solution of the PDE-based image processing and analysis models. Themimetic discretization scheme preserves the continuum properties of the mathematical operators often encountered in the image processing and analysis equations. This contributes to the stability and accuracy of the numerical solution which allows for the improved performance of the approach, as compared to the two very popular numerical solution techniques employed in our experiments. In these experiments, we applied a wide range of noise levels to benchmark images commonly used in the imaging field. These images were restored by applying the well stablished Catté-Lions-Morel-Coll model. Our results show that, for each noise level, the best image that has been restored by solving the PDE with the mimetic discretization scheme, outperforms the best images that have been restored by solving the PDE with the other two methods.
Acknowledgments. This work has been supported in part by NIH RoadMap Initiative award R90 DK07015, and the Computational Science Research Center.
 COMSOL AB, COMSOL Multiphysics Modeling Guide v3.2. COMSOL AB, Burlington, Massachusetts (2005). [ Links ]
 B. Aulbach and S. Hilger, Linear dynamic processes with inhomogeneous time scale. Nonlinear Dynamics and Quantum Dynamical Systems (1990). [ Links ]
 E. Bänsch and K. Mikula, A coarsening finite element strategy in image selective smoothing. Computing and Visualization in Science, 1(1) (1997), 53-61. [ Links ]
 E. Bänsch and K. Mikula, Adaptivity in 3D image processing. Computing and Visualization in Science, 4(1) (2001), 21-30. [ Links ]
 G.I. Barenblatt, M. Bertsch, R. Dal Passo and M. Ughi, A degenerate pseudo-parabolic regularization of a nonlinear forward-backward heat equation arising in the theory of heat and mass exchange in stably stratified turbulent shear flow. SIAM Journal on Mathematical Analysis, 24(6) (1993), 1414-1439. [ Links ]
 B. Benhamouda, Parameter adaptation for nonlinear diffusion in image processing. Master's thesis, University of Kaiserslautern, Kaiserslautern, Germany(1994). [ Links ]
 P.B. Bochev and J.M. Hyman, Principles of mimetic discretizations of differential operators. The IMA Volumes in Mathematics and its Applications, 142 (2006), 89-119. [ Links ]
 M. Bohner and J.E. Castillo, Mimetic methods on measure chains. Computers and Mathematics with Applications, 42 (2001), 705-710. [ Links ]
 A. Canny, A computational approach to edge detection. IEEE Transactions on Pattern Analysis and Machine Intelligence, 8(6) (1986), 679-698. [ Links ]
 J.E. Castillo and R.D. Grone, A matrix analysis approach to higher-order approximations for divergence and gradients satisfying a global conservation law. SIAM Journal on Matrix Analysis and Applications, 25(1) (2003), 128-142. [ Links ]
 J.E. Castillo and J.M. Hyman, Fourth and sixth-order conservative finite difference approximations of the divergence and gradient. Applied Numerical Mathematics, 37(1-2) (2001), 171-187. [ Links ]
 J.E. Castillo, J.M. Hyman, M.J. Shashkov and S. Steinberg, The sensitivity and accuracy of forth order finite-difference schemes on non uniform grids in one dimension. Computers and Mathematics with Applications, 30(8) (1995), 41-55. [ Links ]
 J.E. Castillo, J.M. Hyman, M.J. Shashkov and S. Steinberg, High-order mimetic finite difference methods on non-uniforms grids. Houston J. Math. ICOSAHOM 95, (1996). [ Links ]
 J.E. Castillo and M. Yasuda, Linear systems arising for second-order mimetic divergence and gradient discretizations. Journal of Mathematical Modeling and Algorithms, 4 (2005), 67-82. [ Links ]
 F. Catté, P.-L. Lions, J.-M. Morel and T. Coll, Image selective smoothing and edge detection by nonlinear diffusion. SIAM Journal of Numerical Analysis, 29(1) (1992), 182-193. [ Links ]
 U. Diewald, T. Preusser, M. Rumpf and R. Strzodka, Diffusion models and their accelerated solution in image and surface processing. Acta Mathematica Universitatis Comenianae, 70(1) (2001), 15-34. [ Links ]
 J. Fröhlich and J. Weickert, Image processing using a wavelet algorithm for nonlinear diffusion. Report 104, Laboratory of Technomathematics, University of Kaiserslautern, Kaiserslautern, Germany (1994). [ Links ]
 A. Handlovicová, K. Mikula and A. Sarti, Numerical solution of parabolic equations related to level set formulation of mean curvature flow. Computing and visualization in Science, 1(2) (1999), 179-182. [ Links ]
 A. Handlovicová, K. Mikula and F. Sgallari, Variational numerical methods for solving nonlinear diffusion equations arising in image processing. Journal for Visual Communication and Image Representation, 13(1-2) (2002), 217-237. [ Links ]
 A. Handlovicová, K. Mikula and F. Sgallari, Semi-implicit complementary volume scheme for solving level set like equations in image processing and curve evolution. Numerische Mathematik, 93(4) (2003), 675-669. [ Links ]
 K. Höllig, Existence of infinitely many solutions for a forward-backward heat equation. Transactions of the American Mathematical Society, 278(1) (1983), 299-319. [ Links ]
 K. Höllig and J.A. Nohel, A diffusion equation with a non-monotone constitutive function. In: J.M. Ball (ed.), Proceedings of NATO/London Mathematical Society Conference on Systems of Partial Differential Equation, pages 409-422, (1983). [ Links ]
 J.M. Hyman and S. Steinberg, The convergence of mimetic discretization for rough grids. Computers and Mathematics with Applications, 47(10-11) (2004), 1565-1610. [ Links ]
 J. Kacur and K. Mikula, Solution of nonlinear diffusion appearing in image smoothing and edge detection. Applied Numerical Mathematics, 17(1) (1995), 47-59. [ Links ]
 S. Kichenassamy, The Perona-Malik paradox. SIAM Journal of Applied Mathematics, 57(5) (1997), 1328-1342. [ Links ]
 J.J. Koenderink, The structure of images. Biological Cybernetics, 50(5) (1984), 363-370. [ Links ]
 Z. Krivá and K. Mikula, An adaptive finite volume scheme for solving nonlinear diffusion in image processing. Journal for Visual Communication and Image Representation, 13(1-2) (2002), 22-35. [ Links ]
 X. Li and T. Chen, Nonlinear diffusion with multiple edginess thresholds. Pattern Recognition, 27(8) (1994), 1029-1037. [ Links ]
 D. Marr and E. Hildreth, Theory of edge detection. In: Proceedings of Royal Society of London. Series B, Biological Sciences, Royal Society Publishing, 207 (1980), 187-217. [ Links ]
 K. Mikula, Image processing with partial differential equations, volume 75 of Modern Methods in Scientific Computing and Applications, NATO Science Series II, pages 283-322. Kluwer Academic Publishers, Dodrecht, Netherlands (2002). [ Links ]
 K. Mikula and N. Ramarosy, Semi-implicit finite volume schemes for solving nonlinear diffusion equations in image processing. Numerische Mathematik, 89(3) (2001), 561-590. [ Links ]
 K. Mikula, A. Sarti and C. Lamberti, Geometrical diffusion in 3-d-echocardiography. In: J. Kacur and K. Mikula (eds.), Proceedings of ALGORITMY '97, Conference on Scientific Computing, volume 67, pages 167-181. Acta Mathematica Universitatis Comenianae (1998). [ Links ]
 P. Mrázek, Nonlinear Diffusion for Image Filtering and Monotonicity Enhancement. PhD thesis, Czech Technical University, Prague, Czech Republic (2001). [ Links ]
 M. Nitzberg and T. Shiota, Nonlinear image filtering with edge and corner enhancement. IEEE Transactions on Pattern Analysis and Machine Intelligence, 14(8) (1992), 826-833. [ Links ]
 P. Perona and J. Malik, Scale space and edge detection using anisotropic diffusion. IEEE Transactions on Pattern Analysis and Machine Intelligence, 12(7) (1990), 629-639. [ Links ]
 P. Perona, T. Shiota and J. Malik, Anisotropic diffusion. In: B.M. ter Haar Romeny (ed.), Geometry-Driven Diffusion in Computer Vision, volume 1 of Computational Imaging and Vision, pages 72-92. Springer, Kluwer (1994). [ Links ]
 T. Preusser and M. Rumpf, An adaptive finite element method for large scale image processing. In Proceedings of Scale-Space '99, pages 223-234, (1999). [ Links ]
 J. Weickert, Anisotropic Diffusion in Image Processing. PhD thesis, Universität Kaiserslautern, Kaiserslautern, Germany (1996). [ Links ]
 J. Weickert, Coherence-enhancing diffusion filtering. International Journal of Computer Vision, 31(2/3) (1999), 111-127. [ Links ]
 J. Weickert, Efficient image segmentation using partial differential equations and morphology. Pattern Recognition, 34(9) (2001), 1813-1824. [ Links ]
 J. Weickert, B.M.t.H. Romeny and M.A. Viergever, Efficient and reliable schemes for nonlinear diffusion filtering. IEEE Transactions on Image Processing, 7(3) (1998), 398-410. [ Links ]
 J. Weickert and C. Schnörr, PDE-based preprocessing of medical images. Kunstliche Intelligenz, 3 (2000), 5-10. [ Links ]
 R.T. Whitaker and S.M. Pizer, A multi-scale approach to non-uniform diffusion. Computer Vision, Graphics, and Image Processing: Image Understanding, 57(1) (1993), 99-110. [ Links ]
 A.P. Witkin, Space-scale filtering. In: A. Bundy (ed.), Proceedings of Eighth International Joint Conference on Artificial Intelligence, pages 1019-1022, San Francisco, California, (1983). Morgan Kaufmann. [ Links ]