Preservation of Fine Structures in PDE-Based Image Denoising

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.


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 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: 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: 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: where κ u is the mean curvature defined as 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 u 0 − u and averaged the resulting equation over the whole image domain Ω.Then, for its steady state, we have 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 K 2 ln 1 x 2 /K 2 , for some K > 0, and λ 0 as follows: where c x ρ x /x 1 x 2 /K 2 −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 Advances in Numerical Analysis 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 .

Nonvariational Reformulations
The TV and PM models tend to converge to a piecewise constant image.Such a phenomenon is called the staircasing effect.In order to suppress it, Marquina and Osher 5 suggested to multiply the stationary TV model by a factor of |∇u| as follows: 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 the improved 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 x 2−q , 0 ≤ q < 2, in 2.3 and multiply the resulting equation by |∇u| q 12 as follows: where 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: where S u and C u denote the diffusion term and the constraint coefficient, respectively.

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 u 0 − u.Then, the associated residual equation reads λr −κ u .

3.1
Although the given image u 0 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 u 0 − 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 .

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 where Ω x is a neighborhood of x and σ 2 loc 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 u 0 − u and −κ u defined on Ω x , scaled by 1/σ 2 loc x , that is, 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 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 where 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 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 |u 0 − 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: where M 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 |S u |, while the constraint coefficient C 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 be the net diffusion function.Then, an effective modulator can be defined to impose the net diffusion, N S u , approximately equal over a wide range of |S u | ≥ s 0 , for some s 0 > 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, S u becomes smaller in modulus, which requires N to be strictly increasing.The origin symmetry of N implies that N −s −N s , with which N becomes equally diffusive for both concavities up and down .Such a net diffusion function can be defined, for example, as for some constants η, γ > 0. See Figure 1, where |N s | takes virtually the same values except on smooth regions where |s| is small and therefore the function N may introduce an equalized net diffusion in practice.Incorporating 3.10 , the model 3.8 can be rewritten as 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 C u , and S u loc .

Determination of η and γ
Let u n−1 be the solution in the last time level t n−1 .Then, for the computation of u n , we will try to find the constants η and γ to satisfy for some threshold T > 0 and a user-specified parameter χ.The first identity in 3.12 determines the sharpness of N near the origin; it becomes sharper, as χ → 1.
Note that the net diffusion magnitude |N S u n−1 | satisfying 3.12 is smaller than the magnitude of the original diffusion, |S u n−1 |, at pixels where the image content changes rapidly |S u n−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: Then, since N s is an increasing function, we have

3.14
Thus the net diffusion on oscillatory regions |S u n−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 threshold T 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 |S u n−1 |,

3.15
Since the diffusion magnitude |S u n−1 | evaluated from oscillatory regions is typically larger than its L 2 -average S 0 , 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.

3.16
Advances in Numerical Analysis 9 2 Determine the parameters η and γ for the computation of u n as 3.17 Thus END requires the user to select only a single parameter, χ, which determines the sharpness of the net diffusion function N 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 M s ≡ 1 and therefore the END model 3.11 is reduced to the conventional model 2.11 .

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 C 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 C, as a function of x, t , can be formulated as follows.
1 Select a desirable interval I c c 0 , c 1 for which C x, t ∈ I c , where c 0 ≥ 0 is sufficiently small.
2 Initialize C as a constant as 3 Set C 1 C 0 and for n 2, 3, . ... 3a Compute the absolute residual R n−1 and the correction vector H n−1 as follows:

3.19
where G m is a localized Gaussian smoothing of radius m and 3b Update where ξ n is a scaling factor.For example, when the constraint coefficient is to be limited in a prescribed interval c 0 , c 1 , that is, C x, t ∈ c 0 , c 1 for all x, t , the scaling factor ξ n can be chosen as , n 2, 3, . . . .

3.21
Remark 3.1.The L 2 -average of R n−1 is the standard deviation SD of the the residual, that is,

3.22
It must be larger than the true SD of the noise, when the intermediate residual u 0 − u n−1 involves structural components of the image.Thus the correction vector H 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 C, we may reduce the undesired dissipation from the resulting image.We call the procedure the residual-driven constraint RDC parameterization.

Evaluation of S u loc
The incorporation of S 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 S u is large, which in turn makes the diffusion modulator, M S 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 S 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 S u loc .
Consider an averaging algorithm of the following form: where S 0 S u n−1 and L represents a local averaging scheme of which the weights for each point are given as

3.24
Then we can evaluate S u n−1 loc at the pixel point x as 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 k max 4, 11 − n .

3.26
Thus for the evaluation of S u 0 loc x , for example, the algorithm first makes an average of S u 0 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.

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.

A Linearized Time-Stepping Procedure
Let Δt be the timestep size and t n nΔt, n ≥ 0. Define u n u •, t n .For the diffusion operator, we will exemplify the CCAD model, that is, For 1, 2 and m n, n − 1, let S n−1 be a diffusion matrix approximating the directional diffusion operator as follows: See Section 4.2 below for the spatial numerical scheme.Let where Here C 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 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 the alternating direction implicit ADI method 22, 23 for 4.6 as follows: where u * is an intermediate solution.Since A n−1 1 and A n−1 2 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.

Spatial Difference Schemes
This subsection considers nonstandard difference schemes for S n−1 in 4.2 .Here we will present the scheme for S n−1 1 only; the same scheme is applicable for S n−1 2 .Let Du n−1 i−1/2,j 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

4.8
Advances in Numerical Analysis 13 Define where ε is a small positive constant introduced to prevent d n−1 ij,W from approaching zero.Then the differential operator in 4.2 for 1 can be approximated as where the last term is the harmonic average of d n−1 ij,W and d n−1 ij,E .It follows from 4.2 and 4.10 that the three consecutive nonzero elements of the matrix S n−1 1 corresponding to the pixel x ij become where

4.12
It should be noticed that 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 .

4.13
Suppose that

Advances in Numerical Analysis
Then the Crank-Nicolson method holds the maximum principle, independently of 0 ≤ q < 2, and its solution satisfies The stability condition 4.14 reads 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 S 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.

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 PSNR ≡ 10 log 10 where g denotes the original image and u is the restored image from a noisy image, u 0 , 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 c 0 0.5 and is upper limited by c 1 3.5.The Gaussian smoothing G m 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 Figure 2, 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 Table 1 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.
Figure 3 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.
Figure 4 presents the RDC coefficient C, generated by END RDC restoring Lena.We set the initial value C 0 ≡ c 0 0.5 and the upper limit c 1 3.5.Then the RDC 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 texturefree residual.

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.

Advances in Numerical Analysis
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.

Figure 1 :
Figure 1: The modulator M s and the net diffusion function N s for γ 1 and η 2.

Figure 2 :
Figure 2: 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 of 256 × 256 pixels.