Sixth Mississippi State Conference on Differential Equations and Computational Simulations, Electronic Journal of Differential Equations, Conference 15 (2007), pp. 175–192. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu ftp ejde.math.txstate.edu (login: ftp)

A NON-CONVEX DIFFUSION MODEL FOR SIMULTANEOUS IMAGE DENOISING AND EDGE ENHANCEMENT

SEONGJAI KIM, HYEONA LIM

Abstract. Mathematical restoration models, in particular, total variationbased models can easily lose fine structures during image denoising. In order to overcome the drawback, this article introduces two strategies: the nonconvex (NC) diffusion and the texture-free residual (TFR) parameterization. A non-standard numerical procedure is suggested and its stability is analyzed to effectively solve the NC diffusion model which is mathematically unstable. It has been numerically verified that the resulting algorithm incorporating the NC diffusion and TFR parameterization is able to not only reduce the noise satisfactorily but also enhance edges effectively, at the same time. Various numerical examples are shown to confirm the claim.

1. Introduction

Image restoration is an important image processing (IP) step for various realworld image-related applications and is often necessary as a pre-processing for other imaging techniques such as segmentation and compression. Thus image restoration methods have occupied a peculiar position in IP, computer graphics, and their applications [14, 18, 19, 20, 26].

There have been various partial differential equation (PDE)-based models in image restoration such as the Perona-Malik model [24], the motion by mean curvature [15], the total variation (TV) minimization [17, 25], and color restoration models [4, 11, 15, 16, 27]. These PDE-based models have been extensively studied to answer fundamental questions in image restoration and have allowed researchers and practitioners not only to introduce new mathematical models but also to improve traditional algorithms [1, 5, 10, 21, 29]. Good references to work on them are e.g. Aubert-Kornprobst [2], Osher-Fedkiw [23], and Sapiro [26].

However, most of these PDE-based restoration models and their numerical realizations show a common drawback: loss of fine structures. In particular, they often introduce an unnecessary numerical dissipation on regions where the image content changes rapidly such as on edges and textures. Therefore it is very important to develop mathematical models and/or numerical techniques which can effectively preserve fine structures during the restoration.

2000 Mathematics Subject Classification. 35K55, 65M06, 65M12.

Key words and phrases. Fine structures; denoising; edge enhancement;

nonphysical dissipation; total variation (TV) model; non-convex (NC) diffusion model;

texture-free residual (TFR) parameterization.

\begin{matrix}c \\ ◯ \end{matrix}2007 Texas State University - San Marcos.

Published February 28, 2007.

Partially supported by grants DMS-0312223 and DMS-0630798 from the NSF.

The main objective in this article is to develop a model for simultaneous image denoising and edge enhancement. The objective is quite challenging, because image denoising filters out high-frequency components of the image while the edges in principle show high frequencies. To overcome difficulties in conventional PDEmodels, we will consider a non-convex model incorporating an effective constraint parameter which can both suppress noise and enhance edges at the same time.

An outline of the paper is as follows. In the next section, we briefly review TV-based restoration models and recent studies for the reduction of nonphysical dissipation. Section 3 analyzes sources of the nonphysical dissipation. Based on the analysis, we suggest a non-convex (NC) diffusion model which is unstable mathematically. A non-standard numerical procedure is introduced for the new model and analyzed for stability. In Section 4, we study an effective strategy for the choice of the constraint parameter, called the texture-free residual (TFR) parameterization. Section 5 presents various numerical examples to show effectiveness of the NC diffusion and the TFR parameterization. It has been numerically verified that the resulting algorithm can satisfactorily reduce the noise and enhance edges simultaneously.

2. Preliminaries

2.1. TV-based models. Let f be the observed image of the form

f = K ∗ u + v,  (2.1)

where K is a blurring operator, ∗ represents the convolution, v denotes the noise, and u is the desired image to find. In this article, our main concern is to remove noise. Thus we assume that K = I, the identity operator. A common denoising technique is to minimize a functional of gradient, given as

min_{u} F_{p}(u), F_{p}(u) = ∫_{Ω} |∇u|^{p} dx + \frac{λ}{2} ∥f − u∥^{2},  (2.2)

where λ ≥ 0 and ∥ · ∥ denotes the L^{2}-norm. When p = 1, the first term in F_{1}(u) is called the total variation (TV).

It is often convenient to transform the minimization problem (2.2) into a differential equation, called the Euler-Lagrange equation. Recall that for the minimization problem in 2D of the form

min_{u} ∫_{Ω} G(x, u, ∇u)dx,

the Euler-Lagrange equation [28, §3.3] reads

( − \frac{∂}{∂x} \frac{∂}{∂u_{x}} − \frac{∂}{∂y} \frac{∂}{∂u_{y}} + \frac{∂}{∂u})G = 0.

By applying the variational calculus, one can transform the minimization problem (2.2) into an equivalent differential equation:

−p ∇ · (\frac{∇u}{|∇u|^{2−p}} ) − λ (f − u) = 0.  (2.3)

For a convenient numerical simulation of (2.3), we may parameterize the energy descent direction by an artificial time t; the resulting evolutionary Euler-Lagrange equation can be formulated as:

\frac{∂u}{∂t} − ∇ · (\frac{∇u}{|∇u|^{2−p}} ) = \frac{λ}{p}(f − u).  (2.4)

The restored image becomes closer to f as λ grows. For the above equation, the given image is set for the initial value, i.e., u(t = 0) = f. For the boundary condition, the no-flux condition has often been adopted for simplicity.

Note that to prevent the denominator |∇u| approaching zero, it must be regularized as

|∇u| ≈ |∇^{ε}u| := (u_{x}^{2} + u_{y}^{2} + ε^{2})^{1/2},  (2.5)

for some ε > 0 small. Such a regularization can be a significant source of nonphysical dissipation to be discussed later in this article. As most other PDE-based models, the model (2.4), 1 ≤ p ≤ 2, can easily lose fine structures (due to nonphysical dissipation).

An interesting case in (2.4) is when p = 1, the TV model:

\frac{∂u}{∂t} − ∇ · (\frac{∇u}{|∇u|}) = λ(f − u).  (2.6)

It is known that the TV model tends to transform the image into a collection of locally constant portions, which is called the staircasing. The removal of staircasing has been an interesting research topic [3, 4, 9]. As an anti-staircasing approach, Marquina and Osher [17] introduced the improved TV (ITV) model

\frac{∂u}{∂t} − |∇u| ∇ · (\frac{∇u}{|∇u|}) = λ |∇u| (f − u).  (2.7)

In order to obtain the above model, the authors have first scaled the stationary TV model (2.3) by a factor of |∇u| and then introduced the time parameterization. The scaling can be beneficial in some aspects, as indicated in [23, §11.3]. However, as Cha and Kim have analyzed in their recent paper [6], such a scaling is hardly advantageous for the reduction of nonphysical dissipation. Note that since |∇u| vanishes only in flat regions, the ITV model must have the same steady states as the TV model (2.6).

2.2. Recent studies for the reduction of nonphysical dissipation. For the TV-based image restoration, the staircasing effect is now well understood and relatively easy to handle compared with nonphysical dissipation. Here we review briefly three of recent studies for the reduction of nonphysical dissipation: employment of Besov norm [18], iterative refinement [22], and the method of nonflat time evolution (MONTE) [6].

Employment of Besov norm: As Meyer [18] analyzed, the L^{2}-norm applied to the residual (f − u) in the TV minimization (2.2) is not sensitive enough to distinguish the noise from textures. As a consequence, the residual can easily contain not only the noise but also fine structures and therefore the restored image u turns out to be erroneous and blurry. To reduce the blur associated with the TV model, Meyer suggested the following modified variational problem:

u = arg min_{u∈BV (Ω)}(|u|_{BV} + λ∥f − u∥_{∗}),  (2.8)

where | · |_{BV} is the bounded variation (BV) seminorm defined as

|u|_{BV} = ∫_{Ω} |∇u| dx

and ∥·∥_{∗} denotes the norm in the Besov space B_{∞}^{−1,∞}; see [18, §1.13-1.15] for details.

The Besov norm has shown a better ability in the distinction of different textures and therefore it can preserve more of fine structures. However, the model (2.8) is difficult to minimize, in particular, utilizing the Euler-Lagrange equation approach.

Iterative refinement: In order to effectively suppress nonphysical dissipation, Osher et al. [22] recently suggested an iterative refinement procedure:

• Initialize: u_{0} = e_{0} = 0.

• For k = 1, 2, ⋯: compute u_{k} as the minimizer of the following model

u_{k} = arg min_{u∈BV (Ω)} (|u|_{BV} + \frac{λ}{2}∥f + e_{k−1} − u∥^{2})  (2.9)

and update

e_{k} = e_{k−1} + f − u_{k}.  (2.10)

The authors have proved that u_{k} converges monotonically in L^{2} to f, the noisy image, as k → ∞. However, it should be noticed that such a mathematical property is not necessarily advantageous to a denoising algorithm, because the iterates u_{k} can get noisier as the iteration goes. The iterative refinement recovers not only fine structures but also the noise [22].

The method of nonflat time evolution (MONTE) [6]: The MONTE is a numerical approach for an effective reduction of nonphysical dissipation. Consider a model of the form

\frac{∂u}{∂t} + Lu = R(f − u),  (2.11)

where L is a nonlinear diffusion operator and R denotes a nonnegative constraint term. The above equation is a generalized restoration model which fits with most of known PDE-based models including the TV model (2.6), the ITV model (2.7), the Meyer model (2.8), the Perona-Malik model [24], and the motion by mean curvature.

The MONTE begins with a keen observation: Plugging v = f −u into the model (2.11), we obtain the associated residual equation

\frac{∂v}{∂t} + Rv = Lu.  (2.12)

Then one can see from (2.12) that the residual v increases as the diffusion term Lu grows. Note that the diffusion term is often larger in modulus on fine structures than on smooth regions. Thus the residual becomes extensive on fine structures; as a result, the restored image u must involve a lavish nonphysical dissipation wherever the diffusion term is large in modulus. In order to reduce such nonphysical dissipation effectively, Cha and Kim [6] has introduced the MONTE in which the timestep size (in the numerical time integration) is determined as

Δt(x, t^{n}) = τ F (Lu(x, t^{n−1})),  (2.13)

where τ > 0 is a constant and, for some γ, η > 0,

F (s) = \frac{γ}{1 + η |s|}.

Thus the solution is computed on a nonflat surface of time evolution and the corresponding residual becomes a standard numerical realization of

\frac{∂v}{∂t} + F (Lu) Rv = F (Lu) Lu.  (2.14)

Note that when the ratio γ/η is fixed, the function G(s) (:= F (s)s) becomes flatter on {s : |s| ≥ s_{0} > 0}, approaching γ/η, as η increases. Thus with appropriately chosen parameters (γ and η), the MONTE allows the residual to involve equalized diffusion except for very smooth regions. As a consequence, the restored image can preserve fine structures more effectively, which has also been numerically verified [6]. Since the MONTE can be viewed as a time-stepping numerical procedure, it is applicable to various restoration models. The MONTE also improves computational efficiency for all aforementioned PDE-based models; see [6] for details.

However, to further improve edge-preservation capability, the operator F (Lu) R must be well understood, because otherwise the diffusion equalization can be degraded in (2.14).

The main objectives in this article are (1) to analyze sources of nonphysical dissipation of TV-based models, (2) to introduce a non-convex diffusion model as a remedy for nonphysical dissipation and as a method of edge enhancement, and (3) to study the operator R; in order to remove noise, preserving fine structures more effectively, and at the same time, to enhance edges.

3. Control of Nonphysical Dissipation

This section begins with an analysis for sources of nonphysical dissipation: regularization and convexity. Then we will introduce a new non-convex diffusion model and its non-standard numerical procedure which can control dissipation satisfactorily and enhance edges effectively.

3.1. Regularization and convexity. Due to the regularization in (2.5), the minimization problem in (2.2) must be modified as

min_{u} F_{ε, p}(u), F_{ε, p}(u) = ∫_{Ω} |∇^{ε}u|^{p} dx + \frac{λ}{2} ∥f − u∥^{2},  (3.1)

and its corresponding Euler-Lagrange equation becomes

−p ∇ · (\frac{∇u}{|∇^{ε}u|^{2−p}} ) − λ (f − u) = 0.  (3.2)

In the following, we will see that the regularization becomes a source of numerical dissipation. To understand the combined effects of p and ε in (3.1) (and therefore in (3.2)), we will first measure the functional F_{ε, p}, with λ = 0, for the cartoon images as in Figure 1.

Table 1 presents the values of F_{ε, p} for various selections of ε and p (with λ = 0). Consider the first case F_{0, 2}; the value is smaller for the blurry image in Figure 1(b), which implies that the minimization with (ε, p) = (0, 2) invokes blur. Applying the same logic, one can see that F_{0, 1} does not introduce blur during the minimization. But, when we set (ε, p) = (0.1, 1) as to prevent the denominator approaching zero in the Euler-Lagrange equation, the minimization algorithm will make the image blurry, because the blurry image reduces the value of the functional. Here an interesting case is when p < 1. As one can see from the table that when p = 0.9, the functional has a smaller value for the sharper image. Thus the functional F_{ε, p}, (a)

(b)

Figure 1. A cartoon image in one-dimensional space (a) and its blurry version (b). Gradient magnitudes are indicated in the bottom of each figure.

Table 1. The values of the functional F_{ε, p}, with λ = 0, for images in Figure 1.

F_{ε, p}  Figure 1(a)  Figure 1(b)
F_{0, 2} 0 + 0 + 1^{2} + 0 = 1  0 + 0.5^{2} + 0.5^{2} + 0 = 0.5
F_{0, 1} 0 + 0 + 1 + 0 = 1  0 + 0.5 + 0.5 + 0 = 1
F_{0.1, 1} 0.1 + 0.1 + \sqrt{1.01} + 0.1 ≈ 1.31 0.1 + \sqrt{0.26} + \sqrt{0.26} + 0.1 ≈ 1.22
F_{0, 0.9} 0 + 0 + 1^{0.9} + 0 = 1  0 + 0.5^{0.9} + 0.5^{0.9} + 0 ≈ 1.07

with appropriate p < 1 and ε > 0, can be utilized for edge enhancement, at least for the cartoon image in Figure 1.

We summarize what we have seen with F_{ε, p}, as follows:

• The strictly convex minimization (p > 1) makes images blurrier.

• The TV model itself (p = 1 and ε = 0) may not introduce “blur”. However, its regularization (p = 1 and ε > 0) does.

• When p < 1, the model can make the image sharper. For the non-convex functional, there may not exist the minimizer; however, this does not mean that the discrete version of the equation is not interesting by itself. This case has motivated one of the authors to develop an image denoising model which can preserve fine structures satisfactorily for color images [15].

Based on the above observation, we introduce the following model:

\frac{∂u}{∂t} − |∇^{ε}u|^{1+ω}∇ · (\frac{∇u}{|∇^{ε}u|^{1+ω}} ) = β (f − u), ω ∈ [−1, 1),  (3.3)

where β = β(x, t) ≥ 0. Following the idea of the ITV model (2.7), the parameter β can be chosen as

β = \frac{λ}{1 − ω} |∇^{ε}u|^{1+ω},

which is related to the minimization problem (3.1) with p = 1 − ω. There is no mathematical criterion for the choice of ω; we will set it heuristically between 0 and 1 (non-convex diffusion).

The constraint parameter can play a crucial role in the preservation of fine structures. In §4, we will study an effective numerical procedure for the determination of β(x, t).

In order to understand characteristics of the new model (3.3), we consider its one-dimensional formulation:

\frac{∂u}{∂t} − (u_{x}^{2} + ε^{2})^{(1+ω)/2}(\frac{u_{x}}{(u_{x}^{2} + ε^{2})^{(1+ω)/2}} )_{x} = β (f − u),

or equivalently,

\frac{∂u}{∂t} − \frac{ε^{2} − ωu_{x}^{2}}{u_{x}^{2} + ε^{2}} u_{xx} = β (f − u).  (3.4)

Thus, when ω > 0, the equation (3.4) (and therefore (3.3)) acts as a diffusion model for small derivatives (|u_{x}| ≤ ε/\sqrt{ω}) and as a reverse-time heat equation on regions where the slope is relatively large (|u_{x}| > ε/\sqrt{ω}). Such a property will play an important role for a simultaneous denoising and edge enhancement, although the mathematical model is unstable. In the following subsection, we will introduce an efficient linearized time-stepping procedure along with a non-standard spatial numerical scheme, which makes the resulting algorithm stable for reasonable choices of simulation parameters.

3.2. Numerical procedures. Here we present numerical procedures for (3.3). The spatial scheme can be viewed as a weighted average and the temporal approximation is based on a linearized θ-method and the alternating direction implicit (ADI) perturbation; the resulting algorithm is unconditionally stable for θ = 1 (fully implicit) and stable for a reasonably large timestep size when θ = 1/2 (Crank-Nicolson).

Denote the timestep size by Δt. Set t^{n} = nΔt and u^{n} = u(·, t^{n}) for n ≥ 0. Then, the problem (3.3) can be linearized by evaluating the nonlinear parts from the previous time level. Consider the linearized θ-method for (3.3) of the form:

\frac{u^{n} − u^{n−1}}{Δt} + (A^{n−1} + β^{n}I) [θu^{n} + (1 − θ)u^{n−1}] = β^{n} f, 0 ≤ θ ≤ 1,  (3.5)

where A^{n−1} is a spatial approximation of a linearized diffusion operator, i.e.,

A^{n−1}u^{m} ≈ −|∇^{ε}u^{n−1}|^{1+ω}∇ · (\frac{∇u^{m}}{|∇^{ε}u^{n−1}|^{1+ω}} ), m = n − 1, n.

For an efficient computation of (3.5), the matrix A^{n−1} can be integrated to be separable, i.e., A^{n−1} = A_{1}^{n−1} + A_{2}^{n−1}, with

A_{ℓ}^{n−1}u^{m} ≈ −|∇^{ε}u^{n−1}|^{1+ω} ∂_{x_{ℓ}}(\frac{∂_{x_{ℓ}}u^{m}}{|∇^{ε}u^{n−1}|^{1+ω}} ), ℓ = 1, 2,  (3.6)

where ∇ = (∂_{x_{1}}, ∂_{x_{2}})^{T} = (∂_{x}, ∂_{y})^{T}. Let

B_{ℓ}^{n−1} = A_{ℓ}^{n−1} + \frac{1}{2}β^{n}I, ℓ = 1, 2.

Then, the alternating direction implicit (ADI) method [12, 13] is a perturbation of (3.5), with a splitting error of O(Δt^{2}):

\begin{aligned}[1 + θΔtB_{1}^{n−1}] u^{∗} & = [1 − (1 − θ)ΔtB_{1}^{n−1} − ΔtB_{2}^{n−1}] u^{n−1} + Δtβ^{n} f, \\ [1 + θΔtB_{2}^{n−1}] u^{n} & = u^{∗} + θΔtB_{2}^{n−1}u^{n−1}, \end{aligned}  (3.7)

where u^{∗} is an intermediate solution. When the matrices A_{ℓ}^{n−1} are composed with a 3-point stencil, each sweep in (3.7) can be carried out by inverting a series of tri-diagonal matrices; the ADI is efficient.

Now, we will construct the matrix A_{1}^{n−1}; one can obtain A_{2}^{n−1} similarly. Let D u_{i−1/2,j}^{n−1} be a finite difference approximation of |∇u^{n−1}| evaluated at x_{i−1/2,j}, the mid point of x_{i−1,j} and x_{i,j}. For example, a second-order scheme reads

Du_{i−1/2,j}^{n−1} = ((u_{i,j}^{n−1} − u_{i−1,j}^{n−1})^{2} + [\frac{1}{2}(\frac{u_{i−1,j+1}^{n−1} + u_{i,j+1}^{n−1}}{2} − \frac{u_{i−1,j−1}^{n−1} + u_{i,j−1}^{n−1}}{2} )]^{2})^{1/2}.

(3.8)

Define

d_{ij,W}^{n−1} = [(D u_{i−1/2,j}^{n−1})^{2} + ε^{2}]^{(1+ω)/2}, d_{ij,E}^{n−1} = d_{i+1,j,W}^{n−1} .  (3.9)

Then the differential operators in (3.6) can be approximated as

\begin{matrix}−∂_{x_{1}}(\frac{∂_{x_{1}}u^{m}}{|∇^{ε}u^{n−1}|^{1+ω}} ) ≈ −\frac{1}{d_{ij,W}^{n−1}} u_{i−1,j}^{m} + (\frac{1}{d_{ij,W}^{n−1}} + \frac{1}{d_{ij,E}^{n−1}} )u_{i,j}^{m} − \frac{1}{d_{ij,E}^{n−1}} u_{i+1,j}^{m}, \\ |∇^{ε}u^{n−1}|^{1+ω} ≈ 2\frac{d_{ij,W}^{n−1} · d_{ij,E}^{n−1}}{d_{ij,W}^{n−1} + d_{ij,E}^{n−1}} . \end{matrix}

(3.10)

By combining the right sides of (3.10), we can obtain the three consecutive non-zero elements of the matrix A_{1}^{n−1} corresponding to the pixel x_{ij}:

[A_{1}^{n−1}]_{ij} = (−a_{ij,W}^{n−1} , 2, −a_{ij,E}^{n−1}),  (3.11)

where

a_{ij,W}^{n−1} = \frac{2 d_{ij,E}^{n−1}}{d_{ij,W}^{n−1} + d_{ij,E}^{n−1}} , a_{ij,E}^{n−1} = \frac{2 d_{ij,W}^{n−1}}{d_{ij,W}^{n−1} + d_{ij,E}^{n−1}} .  (3.12)

Note that a_{ij,W}^{n−1} + a_{ij,E}^{n−1} = 2. The above non-standard numerical scheme has been successfully applied as an edge-forming formula for image zooming of arbitrary magnification factors [7, 8].

The following theorem analyzes stability for the θ-method (3.5).

Theorem 3.1. Define β_{0} = min_{i,j,n} β_{ij}^{n} and β_{1} = max_{i,j,n} β_{ij}^{n}. Suppose that the θ-method (3.5) incorporate the numerical scheme in (3.8)-(3.12) and let

(4 + β_{1}) (1 − θ)Δt ≤ 1.  (3.13)

Then (3.5) is stable and holds the maximum principle and its solution satisfies

∥u^{n} − f∥_{∞} ≤ \frac{4}{4 + β_{0}} ∥f∥_{∞}, n ≥ 0.  (3.14)

Proof. We will first prove the following inequality:

min_{i,j} f_{ij} ≤ u_{ij}^{n} ≤ ∥f∥_{∞}, n ≥ 0,  (3.15)

which implies the maximum principle and stability. The equation (3.5) at a point x_{ij} can be written as

\begin{matrix}[1 + (4 + β_{ij}^{n})θΔt] u_{ij}^{n} \\ = θΔt[a_{ij,W}^{n−1} u_{i−1,j}^{n} + a_{ij,E}^{n−1} u_{i+1,j}^{n} + a_{ij,S}^{n−1} u_{i,j−1}^{n} + a_{ij,N}^{n−1} u_{i,j+1}^{n}] \\ + (1 − θ)Δt[a_{ij,W}^{n−1} u_{i−1,j}^{n−1} + a_{ij,E}^{n−1} u_{i+1,j}^{n−1} + a_{ij,S}^{n−1} u_{i,j−1}^{n−1} + a_{ij,N}^{n−1} u_{i,j+1}^{n−1}] \\ + [1 − (4 + β_{ij}^{n})(1 − θ)Δt] u_{ij}^{n−1} + Δtβ_{ij}^{n} f_{ij}. \end{matrix} (3.16)

Let u_{ij}^{n} be a local minimum. Then, it follows from (3.13) and the identity

a_{ij,W}^{n−1} + a_{ij,E}^{n−1} + a_{ij,S}^{n−1} + a_{ij,N}^{n−1} = 4  (3.17)

that each of coefficients in the right side of (3.16), including the term of f_{ij}, is nonnegative and their sum becomes 1 + (4 + β_{ij}^{n})θΔt. Thus, since u_{ij}^{n} is smaller than or equal to the neighboring values, we must have

f_{ij} ≤ u_{ij}^{n}.

The inequality holds for all local minima, which proves the first inequality in (3.15). The same argument can be applied for local maxima to verify the other inequality.

Now, to prove (3.14), let δ_{ij}^{n} = u_{ij}^{n} − f_{ij}, n ≥ 0. Then it follows from (3.16) that

\begin{matrix}[1 + (4 + β_{ij}^{n}) θΔt]δ_{ij}^{n} \\ = θΔt[a_{ij,W}^{n−1} u_{i−1,j}^{n} + a_{ij,E}^{n−1} u_{i+1,j}^{n} + a_{ij,S}^{n−1} u_{i,j−1}^{n} + a_{ij,N}^{n−1} u_{i,j+1}^{n}] \\ + (1 − θ)Δt[a_{ij,W}^{n−1} u_{i−1,j}^{n−1} + a_{ij,E}^{n−1} u_{i+1,j}^{n−1} + a_{ij,S}^{n−1} u_{i,j−1}^{n−1} + a_{ij,N}^{n−1} u_{i,j+1}^{n−1}] \\ − 4Δtf_{ij} + [1 − (4 + β_{ij}^{n})(1 − θ)Δt] δ_{ij}^{n−1}. \end{matrix} (3.18)

Thus, utilizing (3.15) and (3.17), we have

|δ_{ij}^{n}| ≤ \frac{4Δt}{1 + (4 + β_{0}) θΔt} ∥f∥_{∞} + γ_{0} ∥δ^{n−1}∥_{∞},  (3.19)

where

γ_{0} = \frac{1 − (4 + β_{0})(1 − θ)Δt}{1 + (4 + β_{0}) θΔt} = 1 − \frac{(4 + β_{0})Δt}{1 + (4 + β_{0}) θΔt}.

The inequality (3.13) also holds when β_{1} is replaced by β_{0}, i.e.,

(4 + β_{0})Δt ≤ 1 + (4 + β_{0}) θΔt,

which implies 0 ≤ γ_{0} < 1. Since ∥δ^{0}∥_{∞} = 0, we can have

\begin{aligned}|δ_{ij}^{n}| & ≤ \frac{4Δt}{1 + (4 + β_{0}) θΔt} · ∑_{k=0}^{n−1} γ_{0}^{k} · ∥f∥_{∞} \\ & ≤ \frac{4Δt}{1 + (4 + β_{0}) θΔt} · \frac{1}{1 − γ_{0}} · ∥f∥_{∞} \\ & = \frac{4}{4 + β_{0}} ∥f∥_{∞}, \end{aligned}  (3.20)

which completes the proof. □

The above analysis deserves a few remarks.

(1) The stability condition in Theorem 3.1 holds independently on ω ∈ [−1, 1).

(2) The θ-method is unconditionally stable for θ = 1. When θ = 1/2 (Crank-Nicolson), the stability condition reads

Δt ≤ \frac{2}{(4 + β_{1})}.

However, in practice, it is stable for reasonable choices of the timestep size, say, Δt ≤ 2. The Crank-Nicolson scheme is preferred due to a smaller numerical error.

(3) The main diagonal of the matrix A^{n−1} is 4 for all n ≥ 1, independently on ω and image contents as well. The numerical scheme in (3.8)-(3.12), producing such a diffusion matrix, plays important roles in mathematical analysis and practical computation.

(4) As one can see from (3.14), the residual (f −u) becomes smaller in modulus as β increases. A question may occur whether the estimate can be sharper or not. To answer the question, we consider an example. Let f be zero at all pixels except one point, say x_{0}. Then the desired image to recover must be zero: u ≡ 0, which is possible when β = 0 in a neighborhood of x_{0}. In the case, ∥u − f∥_{∞} = ∥f∥_{∞} = f(x_{0}). It seems to the authors that the estimate (3.14) may not be sharper, as long as the maximum norm is the measure.

We will close the section, figuring out the stationary solution of the θ-method (3.5) incorporating the numerical scheme in (3.8)-(3.12). By removing the superscripts n and n − 1 from (3.16), one can obtain

(4 + β_{ij})u_{ij} = u_{i−1,j} + u_{i+1,j} + u_{i,j−1} + u_{i,j+1} + β_{ij}f_{ij},

which is also satisfied by a second-order finite difference solution of the following elliptic equation

−∇ · ∇u = β(f − u).  (3.21)

Hence we may conclude that the iterates of the θ-method (3.5) approach the solution of (3.21), provided that the numerical scheme (3.8)-(3.12) is incorporated and the simulation parameters convince (3.13). Note that the solution of (3.21) can be controlled in some degrees by selecting β appropriately.

In practice, the computation in image restoration stops much earlier than reaching the stationary state; however, the restored image must show some characteristics of the stationary solution. In this point of view, the choice of the constraint parameter β can be crucial for an effective image restoration.

4. The Constraint Parameter

The selection of the constraint parameter has been an interesting problem for PDE-based models, in particular, for the TV model (2.6) and its relatives. The basic mechanism of the PDE-based denoising is diffusion. Thus the parameter β cannot be too large; it must be small enough to introduce a sufficient amount of diffusion. On the other hand, it should be large enough to keep the details in the image. However, in the literature, the parameter has been chosen constant for most cases; the resulting models either smear out fine structures more excessively than desired or leave a certain amount of noise in the restored image.

In order to overcome the difficulty, the parameter must be variable; it can become larger wherever dissipation is excessive, being small else where. In the following, we will consider an automatic and effective numerical method for the determination of the constraint function β(x, t):

(1) Set β as a constant:

β(x, 0) = β_{0}.  (4.1)

(2) For n = 2, 3, ⋯ (2a) Compute the absolute residual and a quantity G_{Res}^{n−1}:

\begin{aligned}R^{n−1} & = |f − u^{n−1}|, \\ G_{Res}^{n−1} & = max (0, S_{m}(R^{n−1}) − \overline{R^{n−1}}), \end{aligned}  (4.2)

where S_{m} is a smoother and \overline{R^{n−1}} denotes the L^{2}-average of R^{n−1}. (2b) Update:

β^{n} = β^{n−1} + γ^{n} G_{Res}^{n−1},  (4.3)

where γ^{n} is a scaling factor having the property: γ^{n} → 0 as n → ∞.

The above procedure has been motivated from the following observation. The PDE-based denoising algorithms tend to have a larger numerical dissipation near fine structures. The tendency in turn makes the residual have structural components on regions where dissipation is excessive. It is more apparent when the absolute residual is smoothed and adjusted by an average. Such structural components in the residual can be viewed as an indicator for nonphysical dissipation. By adding the components to the constraint parameter β, and by allowing the parameter to grow at pixels wherever diffusion is inordinate, we may return fine structures in the residual back to the restored image. The objective is to restore images whose corresponding residual is free of structures and textures. We will call the above procedure, (4.1)-(4.3), the texture-free residual (TFR) parameterization.

In practice, one may wish to limit the parameter in a prescribed interval, i.e., β(x, t) ∈ [β_{0}, β_{1}]. Also one may want to update the parameter a few times only. Note that the θ-method (3.7) converges fast, not larger than 10 iterations, for most images having reasonable noise levels. Thus we can update the TFR parameter β in the first few iterations, say, four. In the case, the scaling factor γ^{n} can be selected as follows:

γ^{n} ∥G_{res}^{n−1}∥_{∞} = η (β_{1} − β_{0}), η = \begin{cases}0.4, & n = 2, \\ 0.3, & n = 3, \\ 0.2, & n = 4, \\ 0.1, & n = 5. \end{cases}  (4.4)

Remark. One may compute the quantity R^{n−1} in a different way:

R^{n−1} = |∇(f − u^{n−1})|.  (4.5)

The above and the one in (4.2) have made no significant differences in practice.

5. Numerical Experiments

In this section, we present some numerical examples to show effectiveness of the non-convex model, (3.3) with ω > 0, and the TFR parameterization. For comparison purposes, we consider four different models: the ITV model with constant β(= β_{c}) (ITV), the non-convex model with constant β(= β_{c}) (NC), the ITV model with variable β (ITV-TFR), and the non-convex model with variable β (NC-TFR).

(a) (b)

(c) (d)

Figure 2. Lenna: (a) The original image g, (b) a noisy image f with a Gaussian noise (PSNR=29.6), and restored images u by using (c) ITV and (d) ITV-TFR.

For all experiments, images are first scaled to have values in [0, 1] and after denoising, it is scaled back for the 8-bit display. For ITV-TFR and NC-TFR, the variable β is chosen as in (4.1)-(4.4). We set the parameters as ε = 0.05, β_{c} = 0.6, β_{0} = 0.5, β_{1} = 5.0, and Δt = 1. For non-convex models, ω = 0.9. (Exceptionally, for the last example in Figure 5, we set ω = 1.5.) The residual is measured in the peak signal-to-noise ratio (PSNR) defined as

PSNR ≡ 10 log_{10} (\frac{∑_{ij} 255^{2}}{∑_{ij}(g_{ij} − u_{ij})^{2}} )dB,

where g denotes the original image and u is the restored image from a noisy image of g.

In Figure 2, we verify effectiveness of the TFR parameterization presented in Section 4. The Lenna image, Figure 2(a), is contaminated by a Gaussian noise of PSNR=29.6 as in Figure 2(b) and restored by the ITV and ITV-TFR models. As one can see from Figure 2(c), the ITV model has introduced a great deal of nonphysical dissipation which in turn makes the restored image blurry (PSNR=33.3). On the other hand, the ITV-TFR model has restored an image as in Figure 2(d),

Table 2. PSNR analysis.

f  ITV  NC ITV-TFR NC-TFR
Lenna 29.6 33.3 (3) 34.2 (3) 35.0 (6) 35.7 (8)
24.7 32.1 (3) 32.6 (5) 32.8 (5) 33.0 (4)
Clock 29.6 32.9 (3) 33.8 (3) 35.3 (7) 35.8 (8)
24.7 31.6 (4) 32.1 (4) 32.5 (5) 32.7 (4)
Park 28.5 28.3 (3) 28.7 (3) 32.1 (26) 32.2 (26)
24.8 27.9 (3) 28.5 (3) 30.1 (17) 30.1 (17)
Milk 28.5 33.8 (5) 34.4 (4) 35.9 (20) 36.2 (18)
22.9 31.9 (6) 32.3 (5) 32.5 (5) 32.8 (3)
Tank 28.5 31.5 (4) 31.8 (4) 33.1 (23) 33.2 (22)
23.9 30.7 (5) 30.8 (4) 31.1 (6) 31.1 (5)

which shows much less nonphysical dissipation (PSNR=35.0). It is apparent from the example that the TFR parameterization is effective in preserving fine structures.

It has been observed from various numerical examples that the NC model improves the restoration quality over the ITV model. However, for most cases, ITV-TFR outperforms the NC model and NC-TFR is the best. In the following example, we will compare performances of the four models.

Figure 3 presents restored images for the Clock image carried out by four different models. The noisy image in Figure 3(b) contains a Gaussian noise of PSNR=29.6. Again for this example, the ITV model introduces a large amount of nonphysical dissipation to lose fine structures as in Figure 3(c). One can see from Figures 3(d) and 3(e) that the ITV-TFR outperforms the NC model, which has been observed for most cases. On the other hand, the NC-TFR model produces the best restored image as shown in Figure 3(f). It preserves fine structures of the image effectively and therefore most details are restored satisfactorily.

In the following, we present a PSNR analysis for a more systematic comparison of the four models. We select three additional images in Figure 4, in addition to the ones dealt in Figures 2 and 3.

Table 2 shows PSNRs of noised images f and restored images u by applying the four different models. For each image, Gaussian noise is applied in two different levels. The integers in parentheses represent the number of CN-ADI iterations which can restore the best image, measured in PSNR and visual verification. As shown in the table, the NC and ITV-TFR models have restored the images better than the ITV model; the NC-TFR model has obtained the best PSNR values for all cases. The TFR parameterization (ITV-TFR and NC-TFR) takes in general more iterations to reach its best restored images. However, this observation does not imply that the TFR parameterization makes the denoising procedure slow down. Instead it is related to the capability that can continue improving image quality for a longer time. For example, for the Park image with PSNR=28.5 in the noisy image f, three iterations of the ITV-TFR and NC-TFR models produce restored images of PSNR=29.5 and PSNR=29.8, respectively. These intermediate images have better PSNRs (and better in visual verification) than those of the constant parameterization. It is safe to say that the TFR parameterization does not slow down the restoration procedure but runs more steps to continuously improve the image quality in higher levels.

(a) (b)

(c) (d)

(e) (f)

Figure 3. Clock: (a) The original image, (b) a noisy image with a Gaussian noise (PSNR=29.6), and restored images by using (c) ITV, (d) NC, (e) ITV-TFR, and (f) NC-TFR.

The inclusion of the NC diffusion (ω > 0) improves the performance only slightly when compared by the PSNR. For example, PSNRs of the NC model are slightly higher than those of the ITV model and PSNRs of the NC-TFR model are improved by only a small amount over the ITV-TFR model. However, the idea of the NC Figure 4. Additional images utilized for a PSNR analysis: Park, Milk, and Tank.

(a) (b)

(c) (d)

Figure 5. Human Brain: (a) The original image and denoised images by (b) ITV, (c) ITV-TFR, and (d) NC-TFR.

diffusion should not be underestimated. In the following example, we will explore characteristics of the NC diffusion.

Figure 6. Horizontal line cuts of the Brain image and its restored images in Figure 5. Each of the curves depicts the first 48 values of the 200th row of the corresponding image.

Figure 5 shows an MRI image for a human brain and the restored images by ITV, ITV-TFR, and NC-TFR. (No noise is applied.) For this real-world example, we have reset ω = 1.5 for NC-TFR to perform its best. The other parameters (ε, β_{c}, β_{0}, β_{1}, and Δt) are kept the same as in previous examples, because they have been verified little sensitive to the image content. The ITV model has reduced the noise reasonably well, although its restored image looks somewhat blurry as in Figure 5(b). However, the ITV-TFR and NC-TFR models produce better images as shown in Figures 5(c) and 5(d), respectively in six and seven iterations. These TFRparameterized models can reduce the noise effectively, preserving fine structures satisfactorily.

In order to show details of the restoration, Figure 6 presents horizontal line segments of the original brain image and the restored images. Each of the curves depicts the first 48 values of the 200th row of the corresponding image. The solid, dot-dashed, dot-solid, and dashed curves are respectively for the original Brain image and three restored images by the ITV, ITV-TFR, and NC-TFR models in Figure 5. The ITV model dissipates the image in some degree, while the TFRparameterized models represent the original image quite accurately. It should be noticed that the restored image by the NC-TFR model shows edges sharper than the original image, particularly on pixels 13-20 and 32-40. It is clear to see that the NC-TFR model has an ability not only to reduce the noise but also to enhance edges. Such simultaneous image denoising and edge enhancement (without an apparent torture of image content) are unique characteristics of the NC-TFR model.

Conclusions

PDE-based models in image restoration show a common drawback, loss of fine structures, due to a severe nonphysical dissipation. After presenting preliminaries, we have analyzed sources of the nonphysical dissipation focusing on the total variation-based models. Based on the analysis, a non-convex (NC) diffusion model has been introduced in order to enhance edges effectively by compensating nonphysical dissipation. Then a non-standard numerical procedure has been suggested and analyzed for stability. In order to further improve the edge-preserving capability, this article has studied a strategy for the choice of an effective constraint parameter, called the texture-free residual (TFR) parameterization. It has been numerically verified that (1) the TFR parameterization improves the image quality continuously and in higher levels, (2) the NC diffusion shows an ability to make the edges sharper than the original real-world images of interests, and (3) the resulting algorithm incorporating both the NC diffusion and the TFR parameterization can simultaneously suppress the noise and enhance edges.

References

[1] L. Alvarez, P. Lions, and M. Morel, “Image selective smoothing and edge detection by nonlinear diffusion. II,” SIAM J. Numer. Anal., vol. 29, pp. 845–866, 1992.

[2] G. Aubert and P. Kornprobst, Mathematical Problems in Image Processing, ser. Applied Mathematics Sciences. New York: Springer-Verlag, 2002, no. 147.

[3] P. V. Blomgren, “Total Variation Methods for Restoration of Vector Valued Images, (Ph.D. thesis),” UCLA Dept. of Math. CAM 98-30, 1998.

[4] P. V. Blomgren and T. F. Chan, “Color TV: Total variation methods for restoration of vector valued images,” IEEE Trans. Image Processing, vol. 7, no. 3, pp. 304–309, 1998.

[5] F. Catte, P. Lions, M. Morel, and T. Coll, “Image selective smoothing and edge detection by nonlinear diffusion.” SIAM J. Numer. Anal., vol. 29, pp. 182–193, 1992.

[6] Y. Cha and S. Kim, “MONTE: The method of nonflat time evolution in PDE-based image restoration,” (submitted).

[7] ——, “Edge-forming methods for color image zooming,” IEEE Trans. Image Process., vo. 15, no. 8, pp. 2315–2323, 2006.

[8] ——, “Edge-forming methods for image zooming,” J. Mathematical Imaging and Vision, vo. 25, no. 3, pp. 353–364, 2006.

[9] A. Chambolle and P. L. Lions, “Image recovery via Total Variational minimization and related problems,” Numer. Math., vol. 76, pp. 167–188, 1997.

[10] T. Chan, S. Osher, and J. Shen, “The digital TV filter and nonlinear denoising,” Department of Mathematics, University of California, Los Angeles, CA 90095-1555, Technical Report #99-34, October 1999.

[11] T. F. Chan and J. Shen, “Variational restoration of non-flat image features: Models and algorithms,” SIAM Journal of Applied Mathematics, vol. 61, no. 4, pp. 1338–1361, 2000.

[12] J. Douglas, Jr. and J. Gunn, “A general formulation of alternating direction methods Part I. Parabolic and hyperbolic problems,” Numer. Math., vol. 6, pp. 428–453, 1964.

[13] J. Douglas, Jr. and S. Kim, “Improved accuracy for locally one-dimensional methods for parabolic equations,” Mathematical Models and Methods in Applied Sciences, vol. 11, pp. 1563–1579, 2001.

[14] R. Gonzalez and R. Woods, Digital Image Processing, 2nd Ed. Upper Saddle River, New Jersey: Prentice-Hall, Inc., 2002.

[15] S. Kim, “PDE-based image restoration: A hybrid model and color image denoising,” IEEE Trans. Image Processing, vol. 15, no. 5, pp. 1163–1170, 2006.

[16] R. Kimmel and N. Sochen, “Orientation diffusion or how to comb a porcupine?” Special issue on PDEs in Image Processing, Computer Vision, and Computer Graphics, J. Visual Comm. Image Representation, vol. 13, pp. 238–248, 2002.

[17] A. Marquina and S. Osher, “Explicit algorithms for a new time dependent model based on level set motion for nonlinear deblurring and noise removal,” SIAM J. Sci. Comput., vol. 22, pp. 387–405, 2000.

[18] Y. Meyer, Oscillating Patterns in Image Processing and Nonlinear Evolution Equations, ser. University Lecture Series. Providence, Rhode Island: American Mathematical Society, 2001, vol. 22.

[19] S. Mitra and G. Sicuranza, Nonlinear Image Processing. San Diego, San Francisco, New York, Boston, London, Sydney, ToKyo: Academic Press, 2001.

[20] J.-M. Morel and S. Solimini, Variational Methods in Image Segmentation, ser. Progress in Nonlinear Differential Equations and Their Applications. Birkhäuser, Boston, 1995, vol. 14.

[21] M. Nitzberg and T. Shiota, “Nonlinear image filtering with edge and corner enhancement,” IEEE Trans. on Pattern Anal. Mach. Intell., vol. 14, pp. 826–833, 1992.

[22] S. Osher, M. Burger, D. Goldfarb, J. Xu, and W. Yin, “Using geometry and iterated refinement for inverse problems (1): Total variation based image restoration,” Department of Mathematics, UCLA, LA, CA 90095, CAM Report #04-13, 2004.

[23] S. Osher and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces. New York: Springer-Verlag, 2003.

[24] P. Perona and J. Malik, “Scale-space and edge detection using anisotropic diffusion,” IEEE Trans. on Pattern Anal. Mach. Intell., vol. 12, pp. 629–639, 1990.

[25] L. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D, vol. 60, pp. 259–268, 1992.

[26] G. Sapiro, Geometric partial differential equations and image analysis. Cambridge: Cambridge University Press, 2001.

[27] N. Sochen, R. Kimmel, and R. Malladi, “A general framework for low level vision,” IEEE Trans. Image Processing, vol. 7, no. 3, pp. 310–318, 1998.

[28] R. Weinstock, Calculus of Variations. New York: Dover Pubilcations, Inc., 1974.

[29] Y. You, W. Xu, A. Tannenbaum, and M. Kaveh, “Behavioral analysis of anisotropic diffusion in image processing,” IEEE Trans. Image Process., vol. 5, pp. 1539–1553, 1996.

Seongjai Kim Department of Mathematics and Statistics, Mississippi State University, Mississippi State, MS 39762-5921, USA

   E-mail address: skim@math.msstate.edu

Hyeona Lim Department of Mathematics and Statistics, Mississippi State University, Mississippi State, MS 39762-5921, USA

   E-mail address: hlim@math.msstate.edu