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.
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. Section 3 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) image , where is the image domain which is an open subset in , we consider a simple noise model of the following form: where is the desired image and 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 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 as .
For a convenient numerical simulation of (2.3), the energy descent direction may be parameterized by an artificial time , that is, can be considered as an evolutionary function and the corresponding evolutionary equation can be obtained by adding on the left side of (2.3).
When , the model (2.3) in its evolutionary form becomes the total variation (TV) model [8] as follows: where 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 and averaged the resulting equation over the whole image domain . Then, for its steady state, we have where 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 , for some , and as follows: where . Note that for the PM model, the function is strictly convex for and strictly concave for . ( is a threshold.) Thus the model can flatten regions of slow transitions (where is big); however, it may enhance image content of large gradient magnitudes such as edges and speckles (where is small).
2.2. 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 as follows: Since 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 , , in (2.3) and multiply the resulting equation by [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 , 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 . It has been numerically verified [12] that for , the CCAD model is superior to the ITV model. Note that the ITV model is a special case () 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 and 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 and be, respectively, the desired image and the residual, . Then, the associated residual equation reads Although the given image is piecewise smooth and is the same as the desired image at (i.e., ), the residual at becomes positive or negative at pixels where the image is concave or convex, respectively. Thus the solution of the TV model at , 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 , which has motivated the method of diffusion modulation.
An alternative to (2.6) is to get a variable parameter by averaging locally where is a neighborhood of and denotes the local noise variance measured over .
Note that the right side of (3.2) can be considered as an inner product of and defined on , scaled by , that is, where denotes a local average over 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 be nonnegative at all image points . The corresponding modification of the constraint parameter in (3.2) reads for . Thus, the stationary TV model incorporated with and scaled by locally averaged curvature 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 is nonzero. When the average curvature term is regularized, the diffusion modulation of the evolutionary TV model (2.4) can be written as where 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 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 at each point in the image domain. In Section 3.3, we will consider an effective alternative to the constraint coefficient . 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 .
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 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 , while the constraint coefficient 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, , approximately equal over a wide range of , for some . 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, becomes smaller in modulus, which requires to be strictly increasing. The origin symmetry of implies that , with which becomes equally diffusive for both concavities (up and down). Such a net diffusion function can be defined, for example, as for some constants . See Figure 1, where takes virtually the same values except on smooth regions (where is small) and therefore the function 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 , and .
3.2. Determination of and
Let be the solution in the last time level . Then, for the computation of , we will try to find the constants and to satisfy for some threshold and a user-specified parameter . The first identity in (3.12) determines the sharpness of near the origin; it becomes sharper, as .
Note that the net diffusion magnitude satisfying (3.12) is smaller than the magnitude of the original diffusion, , at pixels where the image content changes rapidly (), 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 is an increasing function, we have Thus the net diffusion on oscillatory regions () can differ only by a factor of 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 must be small enough to equalize the net diffusion on every interesting oscillatory region including edges and textures. It has been numerically verified that can be chosen as an average of , Since the diffusion magnitude evaluated from oscillatory regions is typically larger than its -average , the threshold 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 , , and compute the -average of as (2)Determine the parameters and for the computation of as
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 . (We set for all experiments presented in this paper.) Note that when , we have 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 , can be formulated as follows. (1) Select a desirable interval for which , where is sufficiently small. (2) Initialize as a constant as (3) Set and for . (3a) Compute the absolute residual and the correction vector as follows: where is a localized Gaussian smoothing of radius and denotes the -average of . (3b) Update where is a scaling factor. For example, when the constraint coefficient is to be limited in a prescribed interval , that is, for all , the scaling factor can be chosen as
Remark 3.1. The -average of is the standard deviation (SD) of the the residual, that is, It must be larger than the true SD of the noise, when the intermediate residual involves structural components of the image. Thus the correction vector 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
The incorporation of 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 is large, which in turn makes the diffusion modulator, , 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 , 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 .
Consider an averaging algorithm of the following form: where and represents a local averaging scheme of which the weights for each point are given as Then we can evaluate at the pixel point as where the iteration count is set large in the beginning of the simulation () and becomes smaller as grows. It has been numerically verified that the resulting END model is efficient when 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 Thus for the evaluation of , for example, the algorithm first makes an average of over a window centered at , 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 be the timestep size and , . Define . For the diffusion operator, we will exemplify the CCAD model, that is, For and , let be a diffusion matrix approximating the directional diffusion operator as follows: (See Section 4.2 below for the spatial numerical scheme.) Let and Define , where Here 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 or equivalently
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 is an intermediate solution. Since and 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 in (4.2). Here we will present the scheme for only; the same scheme is applicable for . Let be a finite difference approximation of evaluated at , the mid point of and . For example, a second-order scheme reads Define where is a small positive constant introduced to prevent from approaching zero. Then the differential operator in (4.2) for can be approximated as where the last term is the harmonic average of and . It follows from (4.2) and (4.10) that the three consecutive nonzero elements of the matrix corresponding to the pixel become where It should be noticed that . 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 Suppose that Then the Crank-Nicolson method holds the maximum principle, independently of , 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, .
One should notice that the entries of the main diagonal of are all 4 for , independently of both 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 ,(iii)END: the END-incorporated model (3.11) of CCAD, (iv)RDC: the RDC parameterization built over CCAD, (v): 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 where denotes the original image and is the restored image from a noisy image, , which is a contamination of by Gaussian noise. For every incorporation of END, the sharpness constant in (3.17) is set 0.6. The RDC coefficient begins with and is upper limited by . The Gaussian smoothing for the absolute residual considered in (3.19) is carried out with six iterations of the nearest 4-point averaging algorithm ().
Public domain images are downloaded, as shown in Figure 2, and then deteriorated by Gaussian noise. For the numerical schemes, we choose in the CN-ADI procedure (4.7) and in (4.9). The ADI iteration is stopped when
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
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 , 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 . 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 by CCAD, END, RDC, and , 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.
(a)
(b)
(c)
(d)
Figure 4 presents the RDC coefficient , generated by restoring Lena. We set the initial value and the upper limit . Then the RDC parameterization in Section 3.3 lets grow, as 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 Figure 5, we apply our new model, , 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, , as before. The residual contains no interesting image features, as one can see.
(a)
(b)
(c)
One may try to choose better parameters (, , and ) for the model , 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.
Acknowledgment
The authors are grateful for the constructive and helpful suggestions from two anonymous referees. The work of S. Kim is supported in part by the NSF Grant DMS-1228337.