Abstract

Mimetic finite difference methods in image processing

C. Bazan^I; M. Abouali^I; J. Castillo^I; P. Blomgren^II

^IComputational Science Research Center, San Diego State University. 5500 Campanile Drive, San Diego, CA 92182-1245, U.S.A. E-mails: carlos.bazan@sdsu.edu / maboualiedu@gmail.com / castillo@myth.sdsu.edu

^IIDepartment of Mathematics and Statistics, San Diego State University. 5500 Campanile Drive, San Diego, CA 92182-7720, U.S.A. E-mail: blomgren.peter@gmail.com

ABSTRACT

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.

1 Introduction

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 u₀(x), u₀→ , defined in a bounded domain Ω ⊂ ², and with Lipschitz continuous boundary ∂ Ω. Usually, the observed image u₀(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 [30]. 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 [16]. Most of today's nonlinear diffusion models stem from the work of Perona and Malik [35], 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|²) → 0 when | ∇u| → ∞ and g(| ∇u|²) → 1 when | ∇u| → 0. One of the diffusivities that Perona and Malik proposed is g(| ∇u|²) = (1 + | ∇u|² / λ²) ^–1, where λ > 0 is a threshold (contrast) parameter that separates forward and backward diffusion [38]. 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 [15]. They proposed replacing the diffusivity g(| ∇u|²) of the Perona-Malik model by a slight variation g(| ∇u _σ|²) with u _σ = G _σ * u, where G _σ is a smooth kernel (Gaussian of variance σ²). We should note that this spatial regularization model belongs to a class of well-posed problems (existence and uniqueness were proven in [15]), and that its successful implementation is contingent on the choosing of an appropriate value for the additional regularization parameter σ. Whitaker and Pizer [43] and Li and Chen [28] suggested making the parameters σ and λ time-dependent, while Benhamouda [6] 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 [15] 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 A_h 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 [1].

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) = u₀(x) in Ω. Then, we can look for a function uⁿ 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 uⁿ 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 W^1,2( Ω) of L₂( Ω)-functions with doubly integrable weak derivatives, and integrate (using Green's theorem and homogeneous Neumann boundary conditions) to obtainthe weak form [30]

Then, for each scale step n, we look for a continuous piecewise linear function u_hⁿ ∈ V_h that satisfies

for all v_h ∈ V_h. Considering the standard Lagrangian base functions ϕ_q∈ V_h, q = 1, ..., L, given by ϕ_q(x_p) = δ_qp (Kronecker delta) for all nodes, the function u_hⁿ is given by

Substituting Eq. (8) into Eq. (7) and considering as test functions v_h = ϕ_q for q = 1, ..., L, we get the Ritz-Galerkin equation for the nodal values u_pⁿ, of the piecewise linear function u_hⁿ:

Then, in each scale step we need to assemble and solve a linear system of the form

for the vector of unknowns (DOF) uⁿ.

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 [8] 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 [10]. 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., e^TD = (-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 [33] 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 [35] suggested using the idea presented by Canny [9] 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 [33]. A time step of δt = 10^-2 was chosen to update all the models. Weickert, Romeny, and Viergever [41] 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, ⁶ , and ⁸ ). The quality of the image restoration after applying the three solution methods to the noisy benchmark images is illustrated in Figs. 3, ⁵ , ⁷ , and ⁹ .

5 Conclusion

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.

Received: 01/XII/10.

Accepted: 05/IV/11.

#CAM-303/10.

Publication Dates

History