# Preservation of Fine Structures in PDE-Based Image Denoising

Authors: Hakran Kim; Velinda R. Calvert; Seongjai Kim
Journal: Advances in Numerical Analysis
 (2012)

---

## Abstract

Image denoising processes often lead to significant loss of fine structures such as edges and textures. This paper studies various innovative mathematical and numerical methods applicable for conventional PDE-based denoising models. The method of diffusion modulation is considered to effectively minimize regions of undesired excessive dissipation. Then we introduce a novel numerical technique for residual-driven constraint parameterization, in order for the resulting algorithm to produce clear images whose corresponding residual is as free of image textures as possible. A linearized Crank-Nicolson alternating direction implicit time-stepping procedure is adopted to simulate the resulting model efficiently. Various examples are presented to show efficiency and reliability of the suggested methods in image denoising.


---

## Body

## 1. Introduction

During the last two decades, as the field of image processing requires higher reliability and efficiency, mathematical techniques have become important components of image processing. Various mathematical frameworks employing powerful tools of partial differential equations (PDEs) and functional analysis have emerged and successfully applied to various image processing tasks, particularly for image denoising and restoration [1–9], see also [10, 11]. Those PDE-based denoising techniques have allowed researchers and practitioners not only to introduce effective new models but also to improve traditional algorithms.However, most PDE-based models tend to either converge to a piecewise constant image or introduce image blur (undesired dissipation), partially because the models are derived by minimizing a functional of the image gradient. As a consequence, the conventional PDE-based models may lose interesting image fine structures during the denoising. In order to reduce the artifact, researchers have studied various mathematical and numerical techniques which either incorporate more effective constraint terms and iterative refinement [5, 12–14] or minimize a functional of second derivatives of the image [15–18]. These new mathematical models may preserve fine structures better than conventional ones; however, more advanced models and appropriate numerical procedures are yet to be developed.Most image denoising models incorporate parameters which are closely related to the noise level. Since it is often the case that the noise level is unknown, the problem of choosing parameters occasionally becomes a difficult task and, as a result, the resulting algorithm may produce unsatisfactory images.This paper suggests the method of diffusion modulation and a residual-driven constraint (RDC) parameterization to obtain clearer images, reduce excessive dissipation effectively, and minimize texture components in the residual. These new methods are modifications or improvements of those studied by one of the authors [19]. In this paper, the diffusion modulation model is derived in a different and clearer way (Section 3.1), the RDC parameterization is considered correspondingly (Section 3.3), an effective strategy is introduced for the evaluation of locally averaged diffusion terms (Section 3.4), and combined effects of these techniques are experimentally analyzed for various natural and medical images (Section 5). The new methods would be applicable for most conventional PDE-based denoising models.An outline of the paper is as follows. In the next section, we briefly review variational denoising models and their nonvariational variants. Section3 suggests a novel denoising model. The model incorporates methods of equalized net diffusion (END) and RDC parameterization, in order for the resulting algorithm to be able to restore clear images with the corresponding residuals left behind being texture-free. In Section 4, we present an efficient time-stepping procedure and anisotropic diffusion difference schemes for the model. Its stability is analyzed in the same section. Section 5 gives certain numerical examples to show effectiveness of the resulting model. It has been numerically verified that our new model is robust and works similarly well for wide ranges of algorithm parameters, producing clear images satisfactorily. It is successfully applicable for denoising of natural images and various medical images as well. Our ultimate goal is to produce clear images of texture-free residual.


## 2. Preliminaries

This section presents a brief review of variational approaches for image denoising, followed by nonvariational variants.


### 2.1. Variational Approaches

Given an observed (noisy) imageu0:Ω→ℝ, where Ω is the image domain which is an open subset in ℝ2, we consider a simple noise model of the following form:
(2.1)u0=u+v,
where u is the desired image and v denotes the noise of a zero mean. Then a common denoising technique is to minimize a functional of the following gradient:
(2.2)u=argminu{∫Ω‍ρ(|∇u|)dx+λ2∫Ω‍(u0-u)2dx},
where ρ is an increasing function (often, convex) and λ≥0 denotes the constraint parameter. It is often convenient to transform the minimization problem (2.2) into a differential equation, called the Euler-Lagrange equation, by applying the variational calculus [20] as follows:
(2.3)-∇·(ρ′(|∇u|)∇u|∇u|)=λ(u0-u).
For an edge-adaptive image denoising, it is required to hold ρ'(x)/x→0 as x→∞.For a convenient numerical simulation of (2.3), the energy descent direction may be parameterized by an artificial time t, that is, u can be considered as an evolutionary function and the corresponding evolutionary equation can be obtained by adding ∂u/∂t on the left side of (2.3).Whenρ(x)=x, the model (2.3) in its evolutionary form becomes the total variation (TV) model [8] as follows:
(2.4)∂u∂t-κ(u)=λ(u0-u),(TV),
where κ(u) is the mean curvature defined as
(2.5)κ(u)=∇·(∇u|∇u|).
It is often the case that the constraint parameter λ is set as a constant, as suggested by Rudin et al. [8]. In order to find the parameter, the authors merely multiplied  (2.4) by (u0-u) and averaged the resulting equation over the whole image domain Ω. Then, for its steady state, we have
(2.6)λ=-1σ21|Ω|∫Ω‍(u0-u)κ(u)dx,
where σ2 is the noise variance. (In [8], λ was evaluated after applying integration by parts, which could avoid approximations of second derivatives.)As another example of (2.3), the Perona-Malik (PM) model [7] can be obtained by setting ρ(x)=(1/2)K2ln(1+x2/K2), for some K>0, and λ=0 as follows:
(2.7)∂u∂t-∇·(c(|∇u|)∇u)=0,(PM),
where c(x)=ρ'(x)/x=(1+x2/K2)-1. Note that for the PM model, the function ρ is strictly convex for x<K and strictly concave for x>K. (K is a threshold.) Thus the model can flatten regions of slow transitions (where c(|∇u|) is big); however, it may enhance image content of large gradient magnitudes such as edges and speckles (where c(|∇u|) is small).


### 2.2. Nonvariational Reformulations

The TV and PM models tend to converge to a piecewise constant image. Such a phenomenon is called thestaircasing effect. In order to suppress it, Marquina and Osher [5] suggested to multiply the stationary TV model by a factor of |∇u| as follows:
(2.8)∂u∂t-|∇u|κ(u)=λ|∇u|(u0-u).(ITV)
Since |∇u| vanishes only on flat regions, its steady state is analytically the same as that of the TV model  (2.4). We will call (2.8) theimproved TV (ITV) model, as called in [21]. Such a nonvariational reformulation may reduce the staircasing effect successfully; however, it is yet to be improved for a better preservation of fine structures.As another variant, we may setρ(x)=x2-q,  0≤q<2, in  (2.3) and multiply the resulting equation by |∇u|q [12] as follows:
(2.9)∂u∂t-|∇u|q∇·(∇u|∇u|q)=β(u0-u),(CCAD),
where
(2.10)β=λ|∇u|q2-q.
The second-order differential operator in (2.9) is closely related to that of the PM model (2.7), in particular when q→2, and therefore we will call  (2.9) the convex-concave anisotropic diffusion (CCAD) model. The CCAD model can be implemented as a stable numerical algorithm for all q∈[0,2). It has been numerically verified [12] that for 1<q<2, the CCAD model is superior to the ITV model. Note that the ITV model is a special case (q=1) of the CCAD model.These reformulations of variational models can restore better images and their properties are now well understood. However, they may still leave interesting image fine structures to the residual. Here, we close the section by writing variational denoising models and their nonvariational reformulations in the following general form:(2.11)∂u∂t+𝒮(u)=𝒞(u)(u0-u),
where 𝒮(u) and 𝒞(u) denote the diffusion term and the constraint coefficient, respectively.


## 2.1. Variational Approaches

Given an observed (noisy) imageu0:Ω→ℝ, where Ω is the image domain which is an open subset in ℝ2, we consider a simple noise model of the following form:
(2.1)u0=u+v,
where u is the desired image and v denotes the noise of a zero mean. Then a common denoising technique is to minimize a functional of the following gradient:
(2.2)u=argminu{∫Ω‍ρ(|∇u|)dx+λ2∫Ω‍(u0-u)2dx},
where ρ is an increasing function (often, convex) and λ≥0 denotes the constraint parameter. It is often convenient to transform the minimization problem (2.2) into a differential equation, called the Euler-Lagrange equation, by applying the variational calculus [20] as follows:
(2.3)-∇·(ρ′(|∇u|)∇u|∇u|)=λ(u0-u).
For an edge-adaptive image denoising, it is required to hold ρ'(x)/x→0 as x→∞.For a convenient numerical simulation of (2.3), the energy descent direction may be parameterized by an artificial time t, that is, u can be considered as an evolutionary function and the corresponding evolutionary equation can be obtained by adding ∂u/∂t on the left side of (2.3).Whenρ(x)=x, the model (2.3) in its evolutionary form becomes the total variation (TV) model [8] as follows:
(2.4)∂u∂t-κ(u)=λ(u0-u),(TV),
where κ(u) is the mean curvature defined as
(2.5)κ(u)=∇·(∇u|∇u|).
It is often the case that the constraint parameter λ is set as a constant, as suggested by Rudin et al. [8]. In order to find the parameter, the authors merely multiplied  (2.4) by (u0-u) and averaged the resulting equation over the whole image domain Ω. Then, for its steady state, we have
(2.6)λ=-1σ21|Ω|∫Ω‍(u0-u)κ(u)dx,
where σ2 is the noise variance. (In [8], λ was evaluated after applying integration by parts, which could avoid approximations of second derivatives.)As another example of (2.3), the Perona-Malik (PM) model [7] can be obtained by setting ρ(x)=(1/2)K2ln(1+x2/K2), for some K>0, and λ=0 as follows:
(2.7)∂u∂t-∇·(c(|∇u|)∇u)=0,(PM),
where c(x)=ρ'(x)/x=(1+x2/K2)-1. Note that for the PM model, the function ρ is strictly convex for x<K and strictly concave for x>K. (K is a threshold.) Thus the model can flatten regions of slow transitions (where c(|∇u|) is big); however, it may enhance image content of large gradient magnitudes such as edges and speckles (where c(|∇u|) is small).


## 2.2. Nonvariational Reformulations

The TV and PM models tend to converge to a piecewise constant image. Such a phenomenon is called thestaircasing effect. In order to suppress it, Marquina and Osher [5] suggested to multiply the stationary TV model by a factor of |∇u| as follows:
(2.8)∂u∂t-|∇u|κ(u)=λ|∇u|(u0-u).(ITV)
Since |∇u| vanishes only on flat regions, its steady state is analytically the same as that of the TV model  (2.4). We will call (2.8) theimproved TV (ITV) model, as called in [21]. Such a nonvariational reformulation may reduce the staircasing effect successfully; however, it is yet to be improved for a better preservation of fine structures.As another variant, we may setρ(x)=x2-q,  0≤q<2, in  (2.3) and multiply the resulting equation by |∇u|q [12] as follows:
(2.9)∂u∂t-|∇u|q∇·(∇u|∇u|q)=β(u0-u),(CCAD),
where
(2.10)β=λ|∇u|q2-q.
The second-order differential operator in (2.9) is closely related to that of the PM model (2.7), in particular when q→2, and therefore we will call  (2.9) the convex-concave anisotropic diffusion (CCAD) model. The CCAD model can be implemented as a stable numerical algorithm for all q∈[0,2). It has been numerically verified [12] that for 1<q<2, the CCAD model is superior to the ITV model. Note that the ITV model is a special case (q=1) of the CCAD model.These reformulations of variational models can restore better images and their properties are now well understood. However, they may still leave interesting image fine structures to the residual. Here, we close the section by writing variational denoising models and their nonvariational reformulations in the following general form:(2.11)∂u∂t+𝒮(u)=𝒞(u)(u0-u),
where 𝒮(u) and 𝒞(u) denote the diffusion term and the constraint coefficient, respectively.


## 3. The New Denoising Model

In this section, we will build a new denoising model, incorporating the method of diffusion modulation and a variable constraint parameterization; the resulting model can restore clear images with the corresponding residuals being as free of image fine structures as possible.We will first try to find the source of undesired diffusion, exemplifying the TV model  (2.4) for simplicity. Let u and r be, respectively, the desired image and the residual, r=u0-u. Then, the associated residual equation reads
(3.1)∂r∂t+λr=-κ(u).
Although the given image u0 is piecewise smooth and is the same as the desired image u at t=0 (i.e., r|t=0≡0), the residual at t>0 becomes positive or negative at pixels where the image is concave or convex, respectively. Thus the solution of the TV model at t>0, u(t)=u0-r(t) must involve undesired dissipation wherever its curvature is nonzero; the larger the curvature is in modulus, the more undesired dissipation occurs. The above observation exemplified with the TV model is clearly applicable for more general PDE-based denoising models of the form (2.11).


### 3.1. The Method of Diffusion Modulation

In order to eliminate the undesired dissipation effectively, we may consider a variable constraint parameterλ(x), which has motivated the method of diffusion modulation.An alternative to (2.6) is to get a variable parameter λ=λ(x) by averaging locally(3.2)λ(x)=-1σloc2(x)1|Ωx|∫Ωx‍(u0-u)κ(u)dx,
where Ωx is a neighborhood of x and σloc2(x) denotes the local noise variance measured over Ωx.Note that the right side of (3.2) can be considered as an inner product of u0-u and -κ(u) defined on Ωx, scaled by 1/σloc2(x), that is,
(3.3)λ(x)=1σloc2(x)〈u0-u,-κ(u)〉Ωx≤1σloc2(x)∥u0-u∥loc,x·∥κ(u)∥loc,x,
where ∥·∥loc,x denotes a local average over Ωx and the Cauchy-Schwarz inequality has been applied. It should also be noticed that the constraint parameter λ in (2.6) is nonnegative for an effective denoising; we will require its local evaluation λ(x) be nonnegative at all image points x. The corresponding modification of the constraint parameter λ(x) in (3.2) reads
(3.4)λ~(x)=max(λ(x),0)=θloc(x)σloc2(x)∥u0-u∥loc,x·∥κ(u)∥loc,x,
for 0≤θloc(x)≤1. Thus, the stationary TV model incorporated with λ~(x) and scaled by locally averaged curvature ∥κ(u)∥loc can be formulated as
(3.5)-1∥κ(u)∥locκ(u)=βloc(u0-u),
where
(3.6)βloc(x)=θloc(x)σloc2(x)∥u0-u∥loc,x.The model in (3.5) is a variant of the stationary TV model, of which the diffusion is modulated by a locally averaged curvature. Thus its solution must be essentially the same as that of the stationary TV model at pixels where ∥κ(u)∥loc,x is nonzero. When the average curvature term is regularized, the diffusion modulation of the evolutionary TV model (2.4) can be written as
(3.7)∂u∂t-1∥κ(u)∥loc+ε0κ(u)=βloc(u0-u),
where ε0>0 is incorporated to prevent the denominator approaching zero.The steady state of (3.7) must be similar to that of the TV model (2.4) incorporating  (3.4), when ε0 is small. However, their numerical solutions, which are in practice obtained terminating much earlier than reaching the steady state, may differ much from each other. The model  (3.7) requires evaluating βloc(x) at each point x in the image domain. In Section 3.3, we will consider an effective alternative to the constraint coefficient βloc. Here, our finding is that (a)
the diffusion operator can be modulated by multiplying the reciprocal of a local average of the diffusion term itself, and(b)
the constraint coefficient can be determined as a function of the residual magnitude|u0-u|.Thus, we may consider the following reformulation of (2.11), of which the diffusion term is modulated by a function of diffusion operator itself:
(3.8)∂u∂t+ℳ(𝒮(u))𝒮(u)=𝒞(u)(u0-u),
where ℳ(s) is a positive function called a modulator. The modulator will play an important role for suppressing undesired excessive dissipation on the regions of large diffusion magnitude |𝒮(u)|, while the constraint coefficient 𝒞(u) must become larger at pixels where the residual reveals structural components of the given image. Such a combination of diffusion modulation and residual-driven variable constraint may (a) minimize undesired dissipation in the first place and (b) return important image features (if any in the residual) back to the restored image. See examples in Section 5.Now, let us find an effective modulator. Let(3.9)𝒩(s)=ℳ(s)·s
be the net diffusion function. Then, an effective modulator can be defined to impose the net diffusion, 𝒩(𝒮(u)), approximately equal over a wide range of |𝒮(u)|≥s0, for some s0>0. However, the net diffusion function must be strictly increasing and origin-symmetric. Note that the model  (3.8) converges in the direction in which the net diffusion decreases in modulus. The convergence must introduce denoising, that is, 𝒮(u) becomes smaller in modulus, which requires 𝒩 to be strictly increasing. The origin symmetry of 𝒩 implies that 𝒩(-s)=-𝒩(s), with which 𝒩 becomes equally diffusive for both concavities (up and down). Such a net diffusion function can be defined, for example, as
(3.10)𝒩(s)=ℳ(s)s=γ1+η|s|s,
for some constants η,γ>0. See Figure 1, where |𝒩(s)| takes virtually the same values except on smooth regions (where |s| is small) and therefore the function 𝒩 may introduce an equalized net diffusion in practice.Figure 1
The modulatorℳ(s) and the net diffusion function 𝒩(s) for γ=1 and η=2.Incorporating (3.10), the model (3.8) can be rewritten as
(3.11)∂u∂t+γ1+η∥𝒮(u)∥loc𝒮(u)=𝒞(u)(u0-u),(END)
which can be viewed as a generalization of the diffusion-modulated TV model  (3.7). We will call  (3.11) the equalized net diffusion (END) formulation of  (2.11). In a time-stepping procedure of  (3.11), the locally averaged curvature term can be evaluated utilizing the last iterate of the solution.In the following, we will present effective strategies for the determination ofγ and η, a variable constraint coefficient 𝒞(u), and ∥𝒮(u)∥loc.


### 3.2. Determination ofη and γ

Letun-1 be the solution in the last time level tn-1. Then, for the computation of un, we will try to find the constants η and γ to satisfy
(3.12)𝒩(T)=χγη=T,0<χ<1,
for some threshold T>0 and a user-specified parameter χ. The first identity in  (3.12) determines the sharpness of 𝒩 near the origin; it becomes sharper, as χ→1.Note that the net diffusion magnitude|𝒩(𝒮(un-1))| satisfying  (3.12) is smaller than the magnitude of the original diffusion, |𝒮(un-1)|, at pixels where the image content changes rapidly (|𝒮(un-1)|>T), while it becomes larger in smooth regions. This would make the END model suppress excessive dissipation in regions of fast transitions and denoise more efficiently in smooth regions.The equations in (3.12) can be easily solved for η and γ, as follows:
(3.13)η=χ1-χ·1T,γ=1+ηT=11-χ.
Then, since 𝒩(s) is an increasing function, we have
(3.14)T=𝒩(T)≤|𝒩(s)|=|γs1+η|s||≤γη=T1χ,for|s|≥T.
Thus the net diffusion on oscillatory regions (|𝒮(un-1)|≥T) can differ only by a factor of 1/χ at most. The parameter χ must be large enough to equalize the net diffusion on oscillatory regions; however, it should not be too large, because otherwise the (almost) flat net diffusion will hardly be effective in denoising.The thresholdT must be small enough to equalize the net diffusion on every interesting oscillatory region including edges and textures. It has been numerically verified that T can be chosen as an average of |𝒮(un-1)|,
(3.15)T=S0:=(1|Ω|∫Ω‍|𝒮(un-1)|2dx)1/2.
Since the diffusion magnitude |𝒮(un-1)| evaluated from oscillatory regions is typically larger than its L2-average S0, the threshold T in  (3.15) suffices to equalize the net diffusion for regions of fine structures.The above arguments for the choice ofη and γ can be summarized as follows.(1)
Select a constantχ, 0<χ<1, and compute the L2-average of |𝒮(un-1)| as
(3.16)S0=(1|Ω|∫Ω‍|𝒮(un-1)|2dx)1/2.(2)
Determine the parametersη and γ for the computation of un as
(3.17)η=χ1-χ·1S0,γ=11-χ.Thus END requires the user to select only a single parameter,χ, which determines the sharpness of the net diffusion function 𝒩 near the origin. It has been numerically verified that it works well with χ=0.5~0.8. (We set χ=0.6 for all experiments presented in this paper.) Note that when χ=0, we have ℳ(s)≡1 and therefore the END model  (3.11) is reduced to the conventional model (2.11).


### 3.3. Residual-Driven Variable Constraint Coefficients

The determination of the constraint parameter has been an interesting problem for PDE-based denoising models, of which the basic mechanism 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 so that the resulting models can either smear out fine structures excessively or maintain an objectionable amount of noise into the restored image.In order to overcome the difficulty, the parameter must be set variable, more precisely,edge-adaptive. Our strategy toward the objective is to (a)
initialize the parameter to be small, and(b)
allow the parameter grow wherever undesired dissipation is excessive, keeping it small elsewhere.Note that the parameter would better be initialized small so that in the early stage of computation, the END model (3.11) can remove the noise effectively and equally everywhere. Then, by letting the parameter grow, the model can return structural components (lost in the residual) back to the image.An automatic and effective numerical method for the determination of the constraint coefficient𝒞, as a function of (x,t), can be formulated as follows. (1)
Select a desirable intervalIc=[c0,c1] for which 𝒞(x,t)∈Ic, where c0≥0 is sufficiently small.(2)
Initialize𝒞 as a constant as
(3.18)𝒞0=𝒞(x,t=0)=c0.(3)
Set𝒞1=𝒞0 and for n=2,3,….
(3a)
Compute the absolute residualℛn-1 and the correction vector ℋn-1 as follows:
(3.19)ℛn-1=|u0-un-1|,ℋn-1=max(0,Gm(ℛn-1)-Av(ℛn-1)),whereGm is a localized Gaussian smoothing of radius m and Av(ℛn-1) denotes the L2-average of ℛn-1.(3b)
Update(3.20)𝒞n=𝒞n-1+ξnℋn-1,whereξn is a scaling factor. For example, when the constraint coefficient is to be limited in a prescribed interval [c0,c1], that is, 𝒞(x,t)∈[c0,c1] for all (x,t), the scaling factor ξn can be chosen as
(3.21)ξn=12n-1·c1-c0∥ℋn-1∥∞,n=2,3,….Remark 3.1.
TheL2-average of ℛn-1 is the standard deviation (SD) of the the residual, that is,
(3.22)Av(ℛn-1)=(1|Ω|∫Ω‍(u0-un-1)2dx)1/2=:σn-1.
It must be larger than the true SD of the noise, when the intermediate residual (u0-un-1) involves structural components of the image. Thus the correction vector ℋn-1 in  (3.19) is nonzero mainly on regions where undesired dissipation is excessive. See Figure 4 below.The above procedure has been motivated also from the observation presented in the beginning of the current section: PDE-based denoising models tend to introduce a large numerical dissipation near fine structures such as edges and textures and the tendency in turn makes the residual have structural components there. Such structural components can be viewed as an indicator for an undesired dissipation. By adding the components to the constraint coefficient𝒞, we may reduce the undesired dissipation from the resulting image. We call the procedure the residual-driven constraint (RDC) parameterization.


### 3.4. Evaluation of∥𝒮(u)∥loc

The incorporation of∥𝒮(u)∥loc in the END model (3.11) has been initiated from a local evaluation of the constraint parameter. At pixels where noise is involved, the magnitude of the diffusion operator 𝒮(u) is large, which in turn makes the diffusion modulator, ℳ(𝒮(u)), small. Thus, the denoising may become slow unless the diffusion operator included in the modulator is appropriately evaluated from its local average.For an effective evaluation of∥𝒮(u)∥loc, it is quite natural to choose a bigger averaging window when the given image involves a larger noise level. Here we present an effective strategy for the evaluation of ∥𝒮(u)∥loc.Consider an averaging algorithm of the following form:(3.23)Sk=ℒSk-1,
where S0=𝒮(un-1) and ℒ represents a local averaging scheme of which the weights for each point are given as
(3.24)116[121242121].
Then we can evaluate ∥𝒮(un-1)∥loc at the pixel point x as
(3.25)∥𝒮(un-1)∥loc(x)=|Sk(x)|,
where the iteration count k is set large in the beginning of the simulation (n=1) and becomes smaller as n grows. It has been numerically verified that the resulting END model is efficient when k is chosen to be at least 4 and not larger than 10. For all examples presented in Section 5 below, the smoothing iteration is chosen as
(3.26)k=max(4,11-n).
Thus for the evaluation of ∥𝒮(u0)∥loc(x), for example, the algorithm first makes an average of 𝒮(u0) over a 21×21 window centered at x, with weights being diminished exponentially in radial directions, and then measures the absolute value of the averaged diffusion at the same point.


## 3.1. The Method of Diffusion Modulation

In order to eliminate the undesired dissipation effectively, we may consider a variable constraint parameterλ(x), which has motivated the method of diffusion modulation.An alternative to (2.6) is to get a variable parameter λ=λ(x) by averaging locally(3.2)λ(x)=-1σloc2(x)1|Ωx|∫Ωx‍(u0-u)κ(u)dx,
where Ωx is a neighborhood of x and σloc2(x) denotes the local noise variance measured over Ωx.Note that the right side of (3.2) can be considered as an inner product of u0-u and -κ(u) defined on Ωx, scaled by 1/σloc2(x), that is,
(3.3)λ(x)=1σloc2(x)〈u0-u,-κ(u)〉Ωx≤1σloc2(x)∥u0-u∥loc,x·∥κ(u)∥loc,x,
where ∥·∥loc,x denotes a local average over Ωx and the Cauchy-Schwarz inequality has been applied. It should also be noticed that the constraint parameter λ in (2.6) is nonnegative for an effective denoising; we will require its local evaluation λ(x) be nonnegative at all image points x. The corresponding modification of the constraint parameter λ(x) in (3.2) reads
(3.4)λ~(x)=max(λ(x),0)=θloc(x)σloc2(x)∥u0-u∥loc,x·∥κ(u)∥loc,x,
for 0≤θloc(x)≤1. Thus, the stationary TV model incorporated with λ~(x) and scaled by locally averaged curvature ∥κ(u)∥loc can be formulated as
(3.5)-1∥κ(u)∥locκ(u)=βloc(u0-u),
where
(3.6)βloc(x)=θloc(x)σloc2(x)∥u0-u∥loc,x.The model in (3.5) is a variant of the stationary TV model, of which the diffusion is modulated by a locally averaged curvature. Thus its solution must be essentially the same as that of the stationary TV model at pixels where ∥κ(u)∥loc,x is nonzero. When the average curvature term is regularized, the diffusion modulation of the evolutionary TV model (2.4) can be written as
(3.7)∂u∂t-1∥κ(u)∥loc+ε0κ(u)=βloc(u0-u),
where ε0>0 is incorporated to prevent the denominator approaching zero.The steady state of (3.7) must be similar to that of the TV model (2.4) incorporating  (3.4), when ε0 is small. However, their numerical solutions, which are in practice obtained terminating much earlier than reaching the steady state, may differ much from each other. The model  (3.7) requires evaluating βloc(x) at each point x in the image domain. In Section 3.3, we will consider an effective alternative to the constraint coefficient βloc. Here, our finding is that (a)
the diffusion operator can be modulated by multiplying the reciprocal of a local average of the diffusion term itself, and(b)
the constraint coefficient can be determined as a function of the residual magnitude|u0-u|.Thus, we may consider the following reformulation of (2.11), of which the diffusion term is modulated by a function of diffusion operator itself:
(3.8)∂u∂t+ℳ(𝒮(u))𝒮(u)=𝒞(u)(u0-u),
where ℳ(s) is a positive function called a modulator. The modulator will play an important role for suppressing undesired excessive dissipation on the regions of large diffusion magnitude |𝒮(u)|, while the constraint coefficient 𝒞(u) must become larger at pixels where the residual reveals structural components of the given image. Such a combination of diffusion modulation and residual-driven variable constraint may (a) minimize undesired dissipation in the first place and (b) return important image features (if any in the residual) back to the restored image. See examples in Section 5.Now, let us find an effective modulator. Let(3.9)𝒩(s)=ℳ(s)·s
be the net diffusion function. Then, an effective modulator can be defined to impose the net diffusion, 𝒩(𝒮(u)), approximately equal over a wide range of |𝒮(u)|≥s0, for some s0>0. However, the net diffusion function must be strictly increasing and origin-symmetric. Note that the model  (3.8) converges in the direction in which the net diffusion decreases in modulus. The convergence must introduce denoising, that is, 𝒮(u) becomes smaller in modulus, which requires 𝒩 to be strictly increasing. The origin symmetry of 𝒩 implies that 𝒩(-s)=-𝒩(s), with which 𝒩 becomes equally diffusive for both concavities (up and down). Such a net diffusion function can be defined, for example, as
(3.10)𝒩(s)=ℳ(s)s=γ1+η|s|s,
for some constants η,γ>0. See Figure 1, where |𝒩(s)| takes virtually the same values except on smooth regions (where |s| is small) and therefore the function 𝒩 may introduce an equalized net diffusion in practice.Figure 1
The modulatorℳ(s) and the net diffusion function 𝒩(s) for γ=1 and η=2.Incorporating (3.10), the model (3.8) can be rewritten as
(3.11)∂u∂t+γ1+η∥𝒮(u)∥loc𝒮(u)=𝒞(u)(u0-u),(END)
which can be viewed as a generalization of the diffusion-modulated TV model  (3.7). We will call  (3.11) the equalized net diffusion (END) formulation of  (2.11). In a time-stepping procedure of  (3.11), the locally averaged curvature term can be evaluated utilizing the last iterate of the solution.In the following, we will present effective strategies for the determination ofγ and η, a variable constraint coefficient 𝒞(u), and ∥𝒮(u)∥loc.


## 3.2. Determination ofη and γ

Letun-1 be the solution in the last time level tn-1. Then, for the computation of un, we will try to find the constants η and γ to satisfy
(3.12)𝒩(T)=χγη=T,0<χ<1,
for some threshold T>0 and a user-specified parameter χ. The first identity in  (3.12) determines the sharpness of 𝒩 near the origin; it becomes sharper, as χ→1.Note that the net diffusion magnitude|𝒩(𝒮(un-1))| satisfying  (3.12) is smaller than the magnitude of the original diffusion, |𝒮(un-1)|, at pixels where the image content changes rapidly (|𝒮(un-1)|>T), while it becomes larger in smooth regions. This would make the END model suppress excessive dissipation in regions of fast transitions and denoise more efficiently in smooth regions.The equations in (3.12) can be easily solved for η and γ, as follows:
(3.13)η=χ1-χ·1T,γ=1+ηT=11-χ.
Then, since 𝒩(s) is an increasing function, we have
(3.14)T=𝒩(T)≤|𝒩(s)|=|γs1+η|s||≤γη=T1χ,for|s|≥T.
Thus the net diffusion on oscillatory regions (|𝒮(un-1)|≥T) can differ only by a factor of 1/χ at most. The parameter χ must be large enough to equalize the net diffusion on oscillatory regions; however, it should not be too large, because otherwise the (almost) flat net diffusion will hardly be effective in denoising.The thresholdT must be small enough to equalize the net diffusion on every interesting oscillatory region including edges and textures. It has been numerically verified that T can be chosen as an average of |𝒮(un-1)|,
(3.15)T=S0:=(1|Ω|∫Ω‍|𝒮(un-1)|2dx)1/2.
Since the diffusion magnitude |𝒮(un-1)| evaluated from oscillatory regions is typically larger than its L2-average S0, the threshold T in  (3.15) suffices to equalize the net diffusion for regions of fine structures.The above arguments for the choice ofη and γ can be summarized as follows.(1)
Select a constantχ, 0<χ<1, and compute the L2-average of |𝒮(un-1)| as
(3.16)S0=(1|Ω|∫Ω‍|𝒮(un-1)|2dx)1/2.(2)
Determine the parametersη and γ for the computation of un as
(3.17)η=χ1-χ·1S0,γ=11-χ.Thus END requires the user to select only a single parameter,χ, which determines the sharpness of the net diffusion function 𝒩 near the origin. It has been numerically verified that it works well with χ=0.5~0.8. (We set χ=0.6 for all experiments presented in this paper.) Note that when χ=0, we have ℳ(s)≡1 and therefore the END model  (3.11) is reduced to the conventional model (2.11).


## 3.3. Residual-Driven Variable Constraint Coefficients

The determination of the constraint parameter has been an interesting problem for PDE-based denoising models, of which the basic mechanism 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 so that the resulting models can either smear out fine structures excessively or maintain an objectionable amount of noise into the restored image.In order to overcome the difficulty, the parameter must be set variable, more precisely,edge-adaptive. Our strategy toward the objective is to (a)
initialize the parameter to be small, and(b)
allow the parameter grow wherever undesired dissipation is excessive, keeping it small elsewhere.Note that the parameter would better be initialized small so that in the early stage of computation, the END model (3.11) can remove the noise effectively and equally everywhere. Then, by letting the parameter grow, the model can return structural components (lost in the residual) back to the image.An automatic and effective numerical method for the determination of the constraint coefficient𝒞, as a function of (x,t), can be formulated as follows. (1)
Select a desirable intervalIc=[c0,c1] for which 𝒞(x,t)∈Ic, where c0≥0 is sufficiently small.(2)
Initialize𝒞 as a constant as
(3.18)𝒞0=𝒞(x,t=0)=c0.(3)
Set𝒞1=𝒞0 and for n=2,3,….
(3a)
Compute the absolute residualℛn-1 and the correction vector ℋn-1 as follows:
(3.19)ℛn-1=|u0-un-1|,ℋn-1=max(0,Gm(ℛn-1)-Av(ℛn-1)),whereGm is a localized Gaussian smoothing of radius m and Av(ℛn-1) denotes the L2-average of ℛn-1.(3b)
Update(3.20)𝒞n=𝒞n-1+ξnℋn-1,whereξn is a scaling factor. For example, when the constraint coefficient is to be limited in a prescribed interval [c0,c1], that is, 𝒞(x,t)∈[c0,c1] for all (x,t), the scaling factor ξn can be chosen as
(3.21)ξn=12n-1·c1-c0∥ℋn-1∥∞,n=2,3,….Remark 3.1.
TheL2-average of ℛn-1 is the standard deviation (SD) of the the residual, that is,
(3.22)Av(ℛn-1)=(1|Ω|∫Ω‍(u0-un-1)2dx)1/2=:σn-1.
It must be larger than the true SD of the noise, when the intermediate residual (u0-un-1) involves structural components of the image. Thus the correction vector ℋn-1 in  (3.19) is nonzero mainly on regions where undesired dissipation is excessive. See Figure 4 below.The above procedure has been motivated also from the observation presented in the beginning of the current section: PDE-based denoising models tend to introduce a large numerical dissipation near fine structures such as edges and textures and the tendency in turn makes the residual have structural components there. Such structural components can be viewed as an indicator for an undesired dissipation. By adding the components to the constraint coefficient𝒞, we may reduce the undesired dissipation from the resulting image. We call the procedure the residual-driven constraint (RDC) parameterization.


## 3.4. Evaluation of∥𝒮(u)∥loc

The incorporation of∥𝒮(u)∥loc in the END model (3.11) has been initiated from a local evaluation of the constraint parameter. At pixels where noise is involved, the magnitude of the diffusion operator 𝒮(u) is large, which in turn makes the diffusion modulator, ℳ(𝒮(u)), small. Thus, the denoising may become slow unless the diffusion operator included in the modulator is appropriately evaluated from its local average.For an effective evaluation of∥𝒮(u)∥loc, it is quite natural to choose a bigger averaging window when the given image involves a larger noise level. Here we present an effective strategy for the evaluation of ∥𝒮(u)∥loc.Consider an averaging algorithm of the following form:(3.23)Sk=ℒSk-1,
where S0=𝒮(un-1) and ℒ represents a local averaging scheme of which the weights for each point are given as
(3.24)116[121242121].
Then we can evaluate ∥𝒮(un-1)∥loc at the pixel point x as
(3.25)∥𝒮(un-1)∥loc(x)=|Sk(x)|,
where the iteration count k is set large in the beginning of the simulation (n=1) and becomes smaller as n grows. It has been numerically verified that the resulting END model is efficient when k is chosen to be at least 4 and not larger than 10. For all examples presented in Section 5 below, the smoothing iteration is chosen as
(3.26)k=max(4,11-n).
Thus for the evaluation of ∥𝒮(u0)∥loc(x), for example, the algorithm first makes an average of 𝒮(u0) over a 21×21 window centered at x, with weights being diminished exponentially in radial directions, and then measures the absolute value of the averaged diffusion at the same point.


## 4. Numerical Schemes

This section presents an efficient time-stepping procedure which incorporates anisotropic diffusion difference schemes for the END model (3.11), followed by its stability analysis.


### 4.1. A Linearized Time-Stepping Procedure

LetΔt be the timestep size and tn=nΔt, n≥0. Define un=u(·,tn). For the diffusion operator, we will exemplify the CCAD model, that is,
(4.1)𝒮(u)=-|∇u|q∇·(∇u|∇u|q),0≤q<2.
For ℓ=1,2 and m=n,n-1, let 𝒮ℓn-1 be a diffusion matrix approximating the directional diffusion operator as follows:
(4.2)𝒮ℓn-1um≈-|∇un-1|q∂xℓ(∂xℓum|∇un-1|q).
(See Section 4.2 below for the spatial numerical scheme.) Let 𝒮n-1=𝒮1n-1+𝒮2n-1 and
(4.3)ℱn-1=γ1+η∥𝒮n-1un-1∥loc.
Define 𝒜n-1=𝒜1n-1+𝒜2n-1, where
(4.4)𝒜ℓn-1:=ℱn-1𝒮ℓn-1+12𝒞n.
Here 𝒞n is the variable constraint coefficient, RDC, presented in Section 3.3. Then, a linearized Crank-Nicolson time-stepping procedure for the END model (3.11) can be formulated as
(4.5)un-un-1Δt+12𝒜n-1(un+un-1)=𝒞nu0,
or equivalently
(4.6)(1+Δt2𝒜n-1)un=(1-Δt2𝒜n-1)un-1+Δt𝒞nu0.One may solve the above system of equations by applying an iterative algebraic solver. However, for an efficiency reason, we in this paper will employ thealternating direction implicit (ADI) method [22, 23] for (4.6) as follows:
(4.7)(1+Δt2𝒜1n-1)u*=(1-Δt2𝒜1n-1-Δt𝒜2n-1)un-1+Δt𝒞nu0,(1+Δt2𝒜2n-1)un=u*+Δt2𝒜2n-1un-1,
where u* is an intermediate solution. Since 𝒜1n-1 and 𝒜2n-1 are tridiagonal matrices, each step of the Crank-Nicolson ADI (CN-ADI) procedure  (4.7) can be carried out by inverting a series of tri-diagonal matrices.


### 4.2. Spatial Difference Schemes

This subsection considers nonstandard difference schemes for𝒮ℓn-1 in  (4.2). Here we will present the scheme for 𝒮1n-1 only; the same scheme is applicable for 𝒮2n-1. Let 𝒟ui-1/2,jn-1 be a finite difference approximation of |∇un-1| evaluated at xi-1/2,j, the mid point of xi-1,j and xi,j. For example, a second-order scheme reads
(4.8)𝒟ui-1/2,jn-1=((ui,jn-1-ui-1,jn-1)2+[12(ui-1,j+1n-1+ui,j+1n-12-ui-1,j-1n-1+ui,j-1n-12)]2)1/2.
Define
(4.9)dij,Wn-1=[(𝒟ui-1/2,jn-1)2+ε2]q/2,dij,En-1=di+1,j,Wn-1,
where ε is a small positive constant introduced to prevent dij,Wn-1 from approaching zero. Then the differential operator in  (4.2) for ℓ=1 can be approximated as
(4.10)-∂x1(∂x1um|∇un-1|q)≈-1dij,Wn-1ui-1,jm+(1dij,Wn-1+1dij,En-1)ui,jm-1dij,En-1ui+1,jm,|∇un-1|q≈2dij,Wn-1·dij,En-1dij,Wn-1+dij,En-1,
where the last term is the harmonic average of dij,Wn-1 and dij,En-1. It follows from  (4.2) and  (4.10) that the three consecutive nonzero elements of the matrix 𝒮1n-1 corresponding to the pixel xij become
(4.11)[𝒮1n-1]ij=(-sij,Wn-1,2,-sij,En-1),
where
(4.12)sij,Wn-1=2dij,En-1dij,Wn-1+dij,En-1,sij,En-1=2dij,Wn-1dij,Wn-1+dij,En-1.
It should be noticed that sij,Wn-1+sij,En-1=2. The above non-standard numerical scheme has been successfully applied for image zooming of arbitrary magnification factors [24–26].The following theorem presents a stability analysis for the Crank-Nicolson method  (4.6) of the END model, which can be proved straightforwardly using our previous results [19, 25].Theorem 4.1.
Let the Crank-Nicolson method (4.6) incorporate difference schemes in (4.8)–(4.12) and
(4.13)F0=mini,j,nℱijn-1,F1=maxi,j,nℱijn-1;C0=mini,j,n𝒞ijn,C1=maxi,j,n𝒞ijn.
Suppose that
(4.14)(4F1+C1)Δt2≤1.
Then the Crank-Nicolson method holds the maximum principle, independently of 0≤q<2, and its solution satisfies
(4.15)∥un-u0∥∞≤4F14F0+C0∥u0∥∞,n≥0.The stability condition  (4.14) reads
(4.16)Δt≤2(4F1+C1).
In practice, it has been verified that the Crank-Nicolson method is stable for reasonable choices of timestep size, say, Δt≤2.One should notice that the entries of the main diagonal of𝒮n-1 are all 4 for n≥1, independently of both q and the image. The numerical schemes in  (4.8)–(4.12), producing such a diffusion matrix, play important roles in mathematical analysis and practical computation.


## 4.1. A Linearized Time-Stepping Procedure

LetΔt be the timestep size and tn=nΔt, n≥0. Define un=u(·,tn). For the diffusion operator, we will exemplify the CCAD model, that is,
(4.1)𝒮(u)=-|∇u|q∇·(∇u|∇u|q),0≤q<2.
For ℓ=1,2 and m=n,n-1, let 𝒮ℓn-1 be a diffusion matrix approximating the directional diffusion operator as follows:
(4.2)𝒮ℓn-1um≈-|∇un-1|q∂xℓ(∂xℓum|∇un-1|q).
(See Section 4.2 below for the spatial numerical scheme.) Let 𝒮n-1=𝒮1n-1+𝒮2n-1 and
(4.3)ℱn-1=γ1+η∥𝒮n-1un-1∥loc.
Define 𝒜n-1=𝒜1n-1+𝒜2n-1, where
(4.4)𝒜ℓn-1:=ℱn-1𝒮ℓn-1+12𝒞n.
Here 𝒞n is the variable constraint coefficient, RDC, presented in Section 3.3. Then, a linearized Crank-Nicolson time-stepping procedure for the END model (3.11) can be formulated as
(4.5)un-un-1Δt+12𝒜n-1(un+un-1)=𝒞nu0,
or equivalently
(4.6)(1+Δt2𝒜n-1)un=(1-Δt2𝒜n-1)un-1+Δt𝒞nu0.One may solve the above system of equations by applying an iterative algebraic solver. However, for an efficiency reason, we in this paper will employ thealternating direction implicit (ADI) method [22, 23] for (4.6) as follows:
(4.7)(1+Δt2𝒜1n-1)u*=(1-Δt2𝒜1n-1-Δt𝒜2n-1)un-1+Δt𝒞nu0,(1+Δt2𝒜2n-1)un=u*+Δt2𝒜2n-1un-1,
where u* is an intermediate solution. Since 𝒜1n-1 and 𝒜2n-1 are tridiagonal matrices, each step of the Crank-Nicolson ADI (CN-ADI) procedure  (4.7) can be carried out by inverting a series of tri-diagonal matrices.


## 4.2. Spatial Difference Schemes

This subsection considers nonstandard difference schemes for𝒮ℓn-1 in  (4.2). Here we will present the scheme for 𝒮1n-1 only; the same scheme is applicable for 𝒮2n-1. Let 𝒟ui-1/2,jn-1 be a finite difference approximation of |∇un-1| evaluated at xi-1/2,j, the mid point of xi-1,j and xi,j. For example, a second-order scheme reads
(4.8)𝒟ui-1/2,jn-1=((ui,jn-1-ui-1,jn-1)2+[12(ui-1,j+1n-1+ui,j+1n-12-ui-1,j-1n-1+ui,j-1n-12)]2)1/2.
Define
(4.9)dij,Wn-1=[(𝒟ui-1/2,jn-1)2+ε2]q/2,dij,En-1=di+1,j,Wn-1,
where ε is a small positive constant introduced to prevent dij,Wn-1 from approaching zero. Then the differential operator in  (4.2) for ℓ=1 can be approximated as
(4.10)-∂x1(∂x1um|∇un-1|q)≈-1dij,Wn-1ui-1,jm+(1dij,Wn-1+1dij,En-1)ui,jm-1dij,En-1ui+1,jm,|∇un-1|q≈2dij,Wn-1·dij,En-1dij,Wn-1+dij,En-1,
where the last term is the harmonic average of dij,Wn-1 and dij,En-1. It follows from  (4.2) and  (4.10) that the three consecutive nonzero elements of the matrix 𝒮1n-1 corresponding to the pixel xij become
(4.11)[𝒮1n-1]ij=(-sij,Wn-1,2,-sij,En-1),
where
(4.12)sij,Wn-1=2dij,En-1dij,Wn-1+dij,En-1,sij,En-1=2dij,Wn-1dij,Wn-1+dij,En-1.
It should be noticed that sij,Wn-1+sij,En-1=2. The above non-standard numerical scheme has been successfully applied for image zooming of arbitrary magnification factors [24–26].The following theorem presents a stability analysis for the Crank-Nicolson method  (4.6) of the END model, which can be proved straightforwardly using our previous results [19, 25].Theorem 4.1.
Let the Crank-Nicolson method (4.6) incorporate difference schemes in (4.8)–(4.12) and
(4.13)F0=mini,j,nℱijn-1,F1=maxi,j,nℱijn-1;C0=mini,j,n𝒞ijn,C1=maxi,j,n𝒞ijn.
Suppose that
(4.14)(4F1+C1)Δt2≤1.
Then the Crank-Nicolson method holds the maximum principle, independently of 0≤q<2, and its solution satisfies
(4.15)∥un-u0∥∞≤4F14F0+C0∥u0∥∞,n≥0.The stability condition  (4.14) reads
(4.16)Δt≤2(4F1+C1).
In practice, it has been verified that the Crank-Nicolson method is stable for reasonable choices of timestep size, say, Δt≤2.One should notice that the entries of the main diagonal of𝒮n-1 are all 4 for n≥1, independently of both q and the image. The numerical schemes in  (4.8)–(4.12), producing such a diffusion matrix, play important roles in mathematical analysis and practical computation.


## 5. Numerical Experiments

This section gives numerical examples to show effectiveness of the method of diffusion modulation (in particular, END (3.11)) and the RDC parameterization introduced in Section 3.3. For comparison purposes, five different models are considered and noted as follows:  (i)
ITV: the ITV model (2.8) with constant λ,(ii)
CCAD: the CCAD model (2.9) with a constant β and q=1.7,(iii)
END: the END-incorporated model (3.11) of CCAD,(iv)
RDC: the RDC parameterization built over CCAD,(v)
END+RDC: the model incorporating both END and RDC over CCAD.For ITV and CCAD, the constant parametersλ and β are chosen from many trials to produce best images compared from PSNR and visual content. By PSNR, we mean the peak signal-to-noise ratio (PSNR) defined as
(5.1)PSNR≡10log10(∑ij2552∑ij(gij-uij)2)dB,
where g denotes the original image and u is the restored image from a noisy image, u0, which is a contamination of g by Gaussian noise. For every incorporation of END, the sharpness constant χ in  (3.17) is set 0.6. The RDC coefficient begins with c0=0.5 and is upper limited by c1=3.5. The Gaussian smoothing Gm for the absolute residual considered in  (3.19) is carried out with six iterations of the nearest 4-point averaging algorithm (m=6).Public domain images are downloaded, as shown in Figure2, and then deteriorated by Gaussian noise. For the numerical schemes, we choose Δt=1 in the CN-ADI procedure  (4.7) and ε=0.05 in (4.9). The ADI iteration is stopped when
(5.2)∥un-un-1∥∞<0.01.Sample images downloaded from public domains: (a) aircraft, (b) balloons, (c) Elaine, (d) fruits, (e) house, (f) Lena, (g) stones, and (h) swan. All are gray-scaled and of256×256 pixels.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Table1 presents PSNRs for the restored images by the five different models. The CCAD model can restore better images than the ITV model for all tested images; it has been numerically verified that the best images can be obtained when q=1.5~2.2, for all tested images. As one can see from the table, END and RDC have improved the restoration quality similarly and more dramatically over CCAD than CCAD does over ITV, for most cases. When END and RDC are combined, they can produce even superior restored images to their individual applications; one of combined effects is to raise the PSNR values by 0.5~1. In practice, the END incorporation increases the computational cost by 30–40% per iteration. However, it is still efficient computationally. The CN-ADI procedure  (4.7) has converged in 3–9 iterations for all models and for all images we have tested.Table 1
PSNR analysis.
Noise
ITV
CCAD
END
RDC
END + RDC
Aircraft
22.83
25.11
26.66
27.37
27.81
28.42
Balloons
22.83
28.06
29.81
31.49
31.34
31.63
Elaine
24.78
29.20
29.80
30.91
30.94
31.28
Fruits
24.77
23.44
24.08
24.81
25.30
25.77
House
22.83
25.33
26.02
27.59
27.84
28.49
Lena
21.25
27.14
27.84
28.86
28.78
29.50
Stones
25.22
28.43
29.31
31.12
30.87
31.61
Swan
23.94
28.89
29.27
30.96
30.85
31.43Figure3 depicts the restored images u by CCAD, END, RDC, and END+RDC, beginning from the noisy Lena image of the PSNR 21.25. As one may have expected, CCAD has restored a somewhat blurry image of which the corresponding residual holds a lot of texture components. END and RDC have produced better images than CCAD when they work separately; their combination has resulted in even a better image whose residual barely shows texture components although twice-magnified. END and RDC, when they are combined, can satisfactorily restore clear images of which the residuals contain no or little interesting image features.Lena: the restored imagesu and the corresponding residuals by (a) CCAD, (b) END, (c) RDC, and (d) END+RDC. For a visualization purpose, the displayed residuals are twice-magnified and shifted, 2(g-u)+128, where g is the original clean image.
(a)
(b)
(c)
(d)Figure 4
Lena: the residual-driven constraint coefficient𝒞, generated by END+RDC.Figure4 presents the RDC coefficient 𝒞, generated by END+RDC restoring Lena. We set the initial value 𝒞0≡c0=0.5 and the upper limit c1=3.5. Then the RDC parameterization in Section 3.3 lets 𝒞n grow, as n increases, wherever excessive undesired dissipation is detected. As one can see from the figure, the variable constraint coefficient exceeds 1 on only regions near textures and edges.In Figure5, we apply our new model, END+RDC, to an ultrasound image, using the same algorithm parameters chosen for the previous example. Denoising must be performed with a special care particularly for such a medical image; every image feature may carry an important information, although it looks like noise. Thus a denoising algorithm for medical images must be successful only if it does not lose any of the image contents, that is, the residual is texture-free. Figure 5(b) shows the restored image by six iterations of the Crank-Nicolson ADI algorithm  (4.7), while Figure 5(c) depicts the corresponding residual. The displayed residual is twice-magnified and shifted, 2*(g-u)+128, as before. The residual contains no interesting image features, as one can see.Ultrasound image in311×382 pixels: (a) the original g, (b) the restored image u obtained by END + RDC, and (c) the corresponding residual, which measures a PSNR of 32.89. For a better visualization, the displayed residual is twice-magnified and shifted, 2*(g-u)+128.
(a)
(b)
(c)One may try to choose better parameters (q, χ, and ξn) for the model END+RDC, a combination of END and RDC built over the CCAD model, to get better images with higher PSNR values. However, the resulting model is robust and works similarly well for wide ranges of the parameters. It is successfully applicable for denoising of natural images and various medical images as well. Here the ultimate goal is to produce clear images of texture-free residual.


## 6. Conclusions

Partial differential equation PDE-based denoising models often lose important fine structures due to an excessive dissipation. In order to minimize such undesired dissipation, we have considered mathematical and numerical techniques applicable for conventional PDE-based denoising models. The method of diffusion modulation is first presented; as an example, the paper has introduced the equalized net diffusion (END) technique, which may suppress excessive dissipation on most regions of fine structures. Then, an effective residual-driven constraint (RDC) parameterization has been studied in order for the resulting algorithm to be able to return important image features in the residual (if any) back to the restored image. These methods are incorporated with a semi-implicit time-stepping procedure, the Crank-Nicolson ADI method. The resulting algorithm has been analyzed for stability and tested to prove its efficiency and reliability in denoising for various natural and medical images. It restores clear images satisfactorily, preserving fine structures successfully.



---

*Source: 750146-2012-10-21.xml*