A Decomposition and Noise Removal Method Combining Diffusion Equation and Wave Atoms for Textured Images

Abstract

We propose a new method that is aimed at denoising images having textures. The method combines a balanced nonlinear partial differential equation driven by optimal parameters, mathematical morphology operators, weighting techniques, and some recent works in harmonic analysis. Furthermore, the new scheme decomposes the observed image into three components that are well defined as structure/cartoon, texture, and noise-background. Experimental results are provided to show the improved performance of our method for the texture-preserving denoising problem.

1. Introduction

A major topic in the image processing community is concerning the extraction of features. Some of the current research has been aimed at decomposing an image into various auxiliary images, each representing a specific set of characteristics such as edges, contours, structure, noise, texture, among others. In this context, an important application is the denoising problem, where many models assume that a noisy image 𝑓 is the sum of two components, 𝑓=𝑢+𝑣, so that 𝑢 contains the structure and the objects of the image (cartoon) while 𝑣 contains the oscillatory characteristics; that is, texture and noise.

In most cases, when an image has only few regions defined by texture, the term 𝑣 is discarded, being characterized by the equation 𝑣=𝑓−𝑢; that is, the models in this category assume that component 𝑣 is only defined by noise, such as is the case of nonlinear models based on adaptive smoothing, anisotropic diffusion, variational methods, and inverse scale space flow [1–8]. On the other hand, there are models that simultaneously explicit both terms: structure 𝑢 and oscillatory 𝑣, so that 𝑓=𝑢+𝑣. This simultaneous decomposition was initially proposed by Meyer [9] from the study of a generalized space to model oscillatory patterns having zero mean, such as noise and texture, which has been recently studied and refined as in [10–16]. In this way, in both described cases, if the observed image 𝑓 has noise and a high texture concentration, then the term 𝑢 will not faithfully represent the recovered image (without noise), since a good part of the fine details of the image (texture) is oscillatory; so consequently they will be annexed to the noise component 𝑣. In all those models, texture and noise are treated equally, which makes it difficult to identify the oscillatory component 𝑣. For example, for real images, it is practically impossible to obtain satisfactory results using some 𝑓=𝑢+𝑣 decomposition model, due to the complex structure and irregular detailing in this category of images.

In order to overcome these problems there is even a third strand of studies that considers a combination of the mentioned techniques allied with recent approaches of harmonic analysis, as in [17–19]. New hybrid models have been proposed with a view to integrate the main advantages of different methods that minimize noise, such as [20–24]. Nevertheless, as most of these models are directly based on curvelets or wave atoms, they tend to make a reconstructed image a bit opaque, besides reproducing oscillating (Gibbs) phenomena, which is an oscillation problem in the frequency domain.

To solve this problem, in this paper we propose a methodology capable of restoring a noisy image having a high concentration of textures and fine details. The proposed model not only keeps the oscillatory well-placed characteristics in the image but also maintains more sensitive textures like intrinsic contours, and edges. Moreover, motivated by [9, 13], we propose a representation of a given image 𝑓 into three terms, 𝑓=𝑢+𝑣+𝑤, with 𝑢 representing the structure or cartoon of the image, 𝑣 only the texture and intrinsic contours and 𝑤 the noise and background. In this sense, another great advantage of our method is that it was constructed so as to satisfy this decomposition.

The proposed scheme seeks to combine the ideas described in [2, 25], that is, a nonlinear PDE (Partial Differential Equation) balanced with an automatic selector of best parameters, based on recent works by [18, 19, 26] on wave atoms, with mathematic morphology operators, such as the top-hat transform one described in [27] and also with the ideas introduced here to synthesize each component of the new decomposition in three terms.

The remainder of the paper is outlined as follows: in Section 2 we briefly describe our method since we use a combination of different approaches in order to treat textured images with noise. We motivate and present our denoising framework in Section 3. Our method is then validated experimentally and statistically in Section 4. We finally summarize and conclude our work in Section 5.

2. The Proposed Method

2.1. Description of the Problem

Let 𝑓 be the observed image with noise and ℎ the original image (without noise), represented by the functions 𝑓∶Ω→ℝ and ℎ∶Ω→ℝ, where Ω is a rectangular region of ℝ². Here, we assumed that both 𝑓 and ℎ can be periodically extended to the ℝ², so as to have 𝑓,ℎ∈𝐿²(ℝ²). We supposed that the noise is additive, that is, 𝑓(𝑥)=ℎ(𝑥)+𝑛(𝑥),𝑥∈Ω, (2.1) where 𝑛 represents the Gaussian type noise with mean 0 e variance 𝜎²_𝑛.

Furthermore, we assumed that the original image ℎ is composed of a structure and an oscillatory pattern component such as texture and irregular details.

The objective here is to minimize the noise level of the input image 𝑓; that is, the impact of noise 𝑛(𝑥) should be minimum in the output image, making it visibly closer to the original image ℎ. The texture and intrinsic contours built-in the image must be maintained and highlighted so that output images 𝑢, 𝑣, and 𝑤 satisfy ℎ≅𝑢+𝑣,𝑛≅𝑤 (2.2) and represent the following features of 𝑓:(i)𝑢: structure, skeleton, or cartoon,(ii)𝑣: texture, intrinsic contours, and irregular details, (iii)𝑤: noise and background.

2.2. The Proposed Scheme

Let 𝑓 be the observed image contaminated with noise as in (2.1). The proposed algorithm can be described as an integrated system having six essential stages, being the first (step 1), the second-to-last (step 5), and the last (step 6) stages of the algorithm where output images are produced, while remaining stages (steps 2, 3, and 4) are where support images are generated that help output components which were synthesized in the previously mentioned steps. All these steps are presented according to the following description.(1)Classical image decomposition. The image decomposition 𝑓 is calculated into two components, 𝑢 and 𝑣, so that 𝑓=𝑢+𝑣, where 𝑢 contains the skeleton or cartoon of 𝑓 and 𝑣 contains the oscillating elements of the image such as noise, irregular details, and texture.(2)Denoising texture-noise component. An iterative procedure is applied to remove noise present in component 𝑣, thereby producing an auxiliary component 𝑣₁ that contains intrinsic contours, parts of the texture, and parts of the edges of 𝑣. (3)Oriented texture support component. An auxiliary image 𝑣₂ is produced containing the warped oscillatory patterns of 𝑣 such as regions characterized by the texture or oscillating details more subtle.(4)Fuzzy representation of edges and textures. A specific morphological filter is applied in 𝑣₂ to obtain a fuzzy edge-texture representation 𝜔_𝑣2∶Ω→[0,1] of the image 𝑣₂. (5)Output of component having only texture. A combination of component 𝑣₁, obtained in step 2, with the component 𝑣₂, obtained in 4, is done, producing a compound image of only texture, intrinsic contours, fine and irregular details, but not the noise. The mentioned image is here identified as 𝑣.(6)Output of restored image and residual component. A composition of the component 𝑢, obtained in 1, is done with the component 𝑣, generated in the previous step. In addition, the image 𝑤 characterized by noise and background (residue) can be obtained, thus satisfying a decomposition of three terms.

Figure 1 shows all the steps previously described while the Figures 2, 3, 5, 6, 7 and 8 show details of each of those steps, respectively. In the next section we will describe each of these steps.

Figure 1 

Illustrative diagram of the proposed algorithm.

(a)

(b)

(c)

Figure 2 

Details of step 1. (a) Input image $f$ with noise, (b) structure/cartoon $u$, and (c) oscillatory patterns $v$ (noise/texture).

(a)

(b)

Figure 3 

Details of step 2. (a) Input image $v$ and (b) intrinsic contours and parts of the texture, here represented by $v_1$. Note in $v_1$ that although there are texture losses close to the neck and background, face intrinsic contours are preserved, such as eyes, mouth, and hair.

3. Description of the Proposed Method

3.1. Classical Image Decomposition

In the first step of the proposed method, the idea is to decompose the initial image 𝑓 into two components 𝑢 and 𝑣, with 𝑓=𝑢+𝑣, such as previously described in this paper. For this purpose, we used the nonlinear anisotropic PDE, proposed by Barcelos et al. [2], with selector of best parameters present in [25]. It is a nonlinear parabolic PDE formulated as from a combination of the classical model from Perona and Malik [6], with the model from Alvarez et al. [1], and with the model proposed by Nordström [8]. Its differential is to eliminate noise through a smoothing process that preserves image edges and contours using automatic detectors. This nonlinear model is presented and detailed in the next section.

In this step, the algorithm aims to smooth the observed image 𝑓 by applying PDE [2], obtaining as a result the component cartoon 𝑢. Then, we do a simple operation 𝑓−𝑢 to determine the component having texture and noise, that is, 𝑣∶=𝑓−𝑢. As the objective is to obtain 𝑢 and 𝑣 well characterized as to structure and texture/noise, respectively, the smoothing process is intensified, being possible to control diffusion velocity through the inclusion of PDE base parameters.

An alternative to implement this step of the algorithm is to use any model that is capable of a good image smoothing, such as [1, 3, 5–8]. Another good alternative is to use a simultaneous cartoon-texture decomposition model, which can be found in [9–14, 16].

The justification to use the nonlinear PDE proposed in [2] is that, from the computational point of view, it is more practical, because besides results being similar to others in literature, it is necessary for use in another step of the our method.

The numerical algorithm used to implement this step of our scheme follows from [2].

3.2. Denoising Texture-Noise Component

The purpose of this step is to remove noise from component 𝑣 (obtained in previous step) in order to minimize the loss of edges, intrinsic contours, textures, and fine details. Again, we used the nonlinear model [2], but now the applied equation is helped by the best parameter selector proposed in [25]. The adopted model is based on the following nonlinear parabolic equation:𝜕𝑣^(𝑡)𝜕𝑡=𝑔|∇𝑣^(𝑡)|div(∇𝑣^(𝑡)|∇𝑣^(𝑡)|)−𝜆(1−𝑔)(𝑣^(𝑡)−𝑣),𝑣(𝑥)⁽⁰⁾=𝑣(𝑥),𝜕𝑣(𝑥)^(𝑡)𝜕𝑛|_𝜕Ω×ℝ+=0,𝑥∈Ω,𝑡∈ℝ₊, (3.1)(3.2) where 𝑣 represents the initial image having texture and noise, 𝑣^(𝑡) is its version on scale 𝑡, 𝑔=𝑔(|∇𝐺_𝜎∗𝑣^(𝑡)|) is a nonnegative and nonincreasing function called diffusivity term, 𝐺_𝜎∗𝑣^(𝑡) determines the convolution of signal 𝑣^(𝑡) with the Gaussian function 𝐺_𝜎, and 𝜆 is a weighting parameter. Here, the symbol |⋅| denotes the Euclidean norm while the constant 𝜎 denotes the standard noise deviation of the image 𝑣.

Generally, the diffusivity term 𝑔, besides being nonincreasing and nonnegative, is such that 𝑔(0)=1, 𝑔(𝑠)→0 when 𝑠→∞ and 𝑔(|𝑠|)∈[0,1]. In this paper we choose 𝑔 based on the Perona and Malik [6] diffusivity term, together with the ideas of Alvarez et al. [1] and Barcelos et al. [25]:𝑔(𝑠)=𝑔(|∇𝐺_𝜎∗𝑣^(𝑡)1|)=1+𝑘|∇𝐺_𝜎∗𝑣^(𝑡)|² (3.3) with 𝐺_𝜎=𝐺_𝜎̂1(𝑥,𝑡)=2𝜎𝜋𝑡̂exp(−|𝑥|²2𝜎𝑡̂), (3.4) where 𝑘=𝑘(𝜎)≥0 is a 𝜎-dependent constant and 𝑡̂ is a scalar variable related to the space of the scale produced by the Gaussian (3.4). The best choice for the scale 𝑡̂ will be better described in following pages.

With the intent to automate the computation of (3.3) and to avoid antiquated choices for 𝑘, we adopt it according the ideas in [25], that is,𝑎𝑘(𝜎)=¹exp(𝑎₂𝑏𝜎),𝜎∈[0,201],₁exp(𝑏₂𝜎),𝜎∈(201,350], (3.5) where 𝑎₁, 𝑎₂, 𝑏₁, and 𝑏₂ are constants such as in [25]. Such a choice brings great advantages such as eliminating an entry parameter in the application of the model and obtaining a good choice for the diffusivity term (3.3). In (3.3), 𝑘 works as an edge selection parameter: for a fixed image with a high 𝑘 value, false edges may be identified while with a small 𝑘 value, only prominent edges will be selected.

The equation (3.1) can be seen as a balancing between smoothing and “keeping close to component 𝑣”. This balancing is managed by the diffusivity term 𝑔, which is used as an edge detector and also to control diffusion velocity. It can be observed that in homogeneous regions of the image we have |∇𝐺_𝜎∗𝑣^(𝑡)| small, which implies in 𝑔∼1. Then, (1−𝑔)∼0 and also the reaction term (𝑣^(𝑡)−𝑣) act in a practically insignificant way in the application of (3.1). Consequently, the diffusion process done by the first parcel of (3.1) is intense; that is, the smoothing will be incisive in these regions. In contrast, for contour regions where |∇𝐺_𝜎∗𝑣^(𝑡)| is big we have 𝑔∼0 and (1−𝑔)∼1, implying that the reaction term (𝑣^(𝑡)−𝑣) proposed by Nordström [8] will retain strongly the initial features of the image 𝑣 under analysis.

The first great advantage of using the nonlinear model (3.1) instead of the classical models in literature (see [1, 3, 6–8]) is that it applies a balanced diffusion controlled by sensitive detector of contours, and since the image 𝑣 can be simultaneously characterized for texture, irregular details, and noise, only regions where there are no contours will be subject to equation diffusion (3.1). Thus, a large part of the texture, edges, and intrinsic contours will be kept in this process. In contrast, it is true that part of those irregular details (warped texture and fine details) will be smoothed in the process. Nevertheless, this deficiency is bypassed by our algorithm in the oriented-texture support component synthesizing step, described in the next section.

The second great advantage is that there are only two parameters in the numerical solution of the model given in [2] to be determined: the constant 𝑘 and the best scale for 𝑡̂, which can be automatically computed.

In [25], the authors linked the 𝑡̂ scale of the Gaussian kernel 𝐺_𝜎 with the noise standard deviation 𝜎_noi of the initial image 𝑣=𝑣(𝑥)^(𝑡=0), which resulted in an estimate for optimal time 𝑇 to stop the evolutionary process (3.1), given by𝜎𝑇=^2noi𝑎, (3.6) where 𝜎 is the standard deviation of image 𝑣 with noise and 𝑎=𝜎, with 𝜎 being the constant present in the Gaussian kernel (3.4).

Encouraged by the authors of [25], we take the optimal smoothing time 𝑡̂ as in (3.6), that is, 𝑡̂=𝑇. Having the optimal stopping time 𝑇, it is possible to automatically obtain (unless the temporal step Δ𝑡) the number of iterations of the model and also the best scale 𝑡̂ for the Gaussian function (3.4). The advantage of being able to automatically obtain the incoming parameters of the model (3.1) makes it pretty efficient and practical, considerably minimizing user intervention at this step of the algorithm. On the other hand, it is true that optimal time 𝑇 is directly related to the standard noise deviation 𝜎_noi. For synthetic images, it is possible to calculate 𝜎_noi, but for real images it is usually not. In the latter case we try to estimate 𝜎_noi based on the visual quality of the image.

To implement the numerical equation (3.1), the idea is to construct an iterative process whose stopping criteria is based on optimal time 𝑇 (3.6), as presented in [2]. In this case, 𝑇 is used in the Gaussian function (3.3), as was previously described. Also, the 𝑁 number of iterations is calculated based on 𝑇, which is given by 𝑇𝑁=.Δ𝑡 (3.7)

3.3. Oriented Texture Support Component

That is one of the most important steps of the proposed method, because it is in it that we extract the oriented texture and most of the oscillatory details of the image. To do this, we use a recent study of wavelet variants presented in [18, 19] for texture analysis, which became known as wave atoms.

Wave atoms are a variant obtained through a 2D wavelet packet obeying the important parabolic scaling relation wavelenght∼(diameter)², which improves the sparse representation of certain oscillatory patterns when compared to more traditional expansions, such as wavelets, gabor atoms, or curvelets. To be more precise, it means that the warped oscillatory functions (oriented textures) have a significant sparse expansion in wave atoms than in other representations of the literature.

Compared to other transforms, wave atoms have two great advantages: the ability to arbitrarily adapt in localities defined by a certain pattern and the ability to sparsely represent anisotropic patterns aligned with the axes. Wave atoms composition elements have a high direction sensibility and anisotropy, which makes them ideal to apply wherever the intention is to identify regions characterized by oscillatory patterns such as texture, as is the case here presented.

In the following, based on [18], we will give a brief explanation about wave atoms mathematical precedents. For more details, also see [19, 28].

Consider wave atoms given by 𝜑_𝜇, with subscript 𝜇=(𝑗,𝐦,𝐧)=(𝑗,𝑚₁,𝑚₂,𝑛₁,𝑛₂). The five quantities above are integer values and index a point (𝑥_𝜇,𝜔_𝜇) in a space-phase such as 𝑥_𝜇=2^−𝑗𝐧, 𝜔_𝜇=𝜋2^𝑗𝐦, 𝑐₁2^𝑗≤max_𝑖=1,2|𝑚_𝑖|≤𝑐₂2^𝑗, where 𝑐₁ and 𝑐₂ are two positive constants. According to [18], the elements of a frame of wave packets 𝜑_𝜇 are called wave atoms when|𝜑_𝜇(𝜔)|≤𝐶_𝑀⋅2^−𝑗(1+2^−𝑗|𝜔−𝜔_𝜇|)^−𝑀+𝐶_𝑀⋅2^−𝑗(1+2^−𝑗|𝜔+𝜔_𝜇|)^−𝑀|𝜑_𝜇(𝑥)|≤𝐶_𝑀⋅2^𝑗(1+2^𝑗|𝑥−𝑥_𝜇|)^−𝑀,∀𝑀>0. (3.8)

To construct wave atoms for our problem, we first considered the case of a family 1𝐷 of wave packets 𝜓^𝑗_𝑚,𝑛(𝑥)≥0, 𝑚≥0, 𝑛∈ℕ, in the scope of above same conditions. Let 𝜙 be a continuous real function given in [−7𝜋/6,5𝜋/6], such that for |𝜔|≤𝜋/3, 𝜙((𝜋/2)−𝜔)²+𝜙((𝜋/2)+𝜔)²=1, and 𝜙(−(𝜋/2)−2𝜔)=𝜙((𝜋/2)+𝜔).

Considering 𝜓⁰_𝑚(𝑡)=2Re{exp(𝑖𝜋(𝑚+1/2)𝑡)𝜈((−1)^𝑛(𝑡−1/2))}, with 𝜈(𝑡)=(1/2𝜋)∫𝜙(𝜔)exp(𝑖𝜔𝑡)𝑑𝜔, the Fourier transform of 𝜓⁰_𝑚 is given by𝜓⁰_𝑚(𝜔)=exp(−𝑖𝜔2)[exp(𝑖𝛼_𝑚)𝜙(𝜖_𝑚1(𝜔−𝜋(𝑚+2)))+exp(−𝑖𝛼_𝑚)𝜙(𝜖_𝑚+11(𝜔+𝜋(𝑚+2)))], (3.9) where 𝜖_𝑚=(−1)^𝑚, 𝛼_𝑚=(𝜋/2)(𝑚+(1/2)). In this case, 𝜙 must be such that ∑^∞_𝑚=0|𝜓⁰_𝑚|²=1.

Thus, we can write functions that make up the base as 𝜓^𝑗_𝑚,𝑛(𝑥)=𝜓^𝑗_𝑚(𝑥−2^−𝑗𝑛)=2^𝑗/2𝜓⁰_𝑚(2^𝑗𝑥−𝑛), whose coefficients can be obtained by 𝑐_{𝑗,𝑚,𝑛}=∫𝜓^𝑗_𝑚,𝑛𝑣(𝑥)𝑑𝑥=(1/2𝜋)∫exp(𝑖2^−𝑗−𝑛𝜔)𝜓^𝑗_𝑚(𝜔)𝑣(𝜔)𝑑𝜔. Here, 𝑣 represents the observed signal.

According [18], the extension to two-dimensional case (2𝐷 version) can be computed by 𝜑⁺_𝜇(𝑥₁,𝑥₂)∶=𝜓^𝑗_𝑚1,𝑛1(𝑥₁)𝜓^𝑗_𝑚2,𝑛2(𝑥₂),𝜑⁻_𝜇(𝑥₁,𝑥₂)∶=𝐻𝜓^𝑗_𝑚1,𝑛1(𝑥₁)𝐻𝜓^𝑗_𝑚2,𝑛2(𝑥₂), where 𝜇=(𝑗,𝐦,𝐧) and a second equation is based on Hilbert transform relative to these wavelet packets. Therefore, the combination 𝜑_𝜇⁽¹⁾=(𝜑⁺_𝜇+𝜑⁻_𝜇)/2 and 𝜑_𝜇⁽²⁾=(𝜑⁺_𝜇−𝜑⁻_𝜇)/2 make up a wave atom frame.

For the denoising problem, it is recommendable to use wave atom shrinkage, which is formulated in most cases by 𝑣_𝑐=∑_𝜇𝜃(𝑐_𝜇⁽¹⁾(𝑣))𝜑_𝜇⁽¹⁾+𝜃(𝑐_𝜇⁽²⁾(𝑣))𝜑_𝜇⁽²⁾, where 𝜃=𝜃_𝛾 is a thresholding function, here adopted by𝜃_𝛾=𝑥,|𝑥|≥𝛾,0,|𝑥|<𝛾, (3.10) where 𝛾 begin the threshold value.

In this work we use wave atom shrinkage to extract oriented texture from image component 𝑣. We take a transform based on wave atoms as follows: 𝑣₂∶=𝑇_𝛾𝑣=(WA)⁻¹𝜃_𝛾(WA)(𝑣), (3.11) where WA denotes the wave atom transform (we use the version 2𝐷), (WA)⁻¹ the inverse transform (see [18, 19]), 𝜃_𝛾 is the threshold function (3.10) mentioned earlier, and 𝑣₂ is the output component. Here, the nonlinear operator 𝑇 will keep and highlight important characteristics of the examined image 𝑣 such as warped texture, oscillating details, and irregular patterns.

The implementation of our wave-atoms shrinkage transform consists of the execution of three steps. First, we apply wave atom transform WA. Next, we remove some insignificant coefficients from the wave atoms through thresholding (3.10). Finally, we apply the inverse transform WA⁻¹ in order to reconstruct the signal containing remaining wave atoms coefficients. For a computational discretization of wave atoms present in the proposed model, we used the WaveatomLab packet, which can be found on the site http://www.waveatom.org/.

3.3.1. Wave Atoms × Other Systems

The main advantage of methods based on wavelet variants is space-frequency localization and multiscale view of the features of surfaces. However, it is known that traditional wavelets are not good to analyze surfaces with “scratches” or textures, due to wavelets ignoring properties defined by geometric features of edges and textures, which leads to strong oscillation along these “scratches”.

In contrast, curvelet transforms, such as [29, 30], are multiscale geometric transforms, which constitute an optimal sparse representation of objects characterized by singularities 𝐶². Nevertheless, they do not work so efficiently when the objective is to represent oscillating textures; that is, they are not efficient to characterize surfaces having warped textures such as fingerprints, photographs, among other types of images.

Curvelets are good for representing edges while wave atoms are good for representing oscillatory patterns and textures. Wave atom texture-shape elements not only capture the coherence along the oscillations like curvelets but also take into consideration patterns across the oscillations (see Figures 4(a) and 4(b)). Since the objective is to characterize surfaces with oriented textures, we have a great advantage in applying wave atom transforms for this purpose.

(a)

(b)

Figure 4 

Elements of the curvelets and wave atoms. (a) A digital curvelet [29] and (b) a digital wave atom [18].

(a)

(b)

Figure 5 

Details of step 3. (a) Input image $v$ and (b) oriented texture and irregular details $v_2$. Note that contrary to what happened in Figure 3(b), all oriented texture is highlighted in $v_2$; however, intrinsic contours (details of face) are lost or left with no definition.

(a)

(b)

Figure 6 

Details of step 4. (a) Input image $v_2$ and (b) highlight of oriented texture of $v_2$ with background correction (in the scale $[\mathrm{0,1}]$), here denoted by ${\omega }_{v_2}$.

(a)

(b)

(c)

Figure 7 

Details of step 5. (a) Input image $v_1$, which is characterized by intrinsic contours, (b) fuzzy representation ${\omega }_{v_2}$ of warped texture, and (c) combination of features of $v_1$ with ${\omega }_{v_2}$, resulting in the final component $\tilde{v}$ defined by oscillatory patterns of observed image $f$, except noise.

(a)

(b)

(c)

(d)

Figure 8 

Details of final step 6. (a) Input image $u$, (b) texture $\tilde{v}$, (c) recovered image $\tilde{h}$, obtained by the addition of $u$ to $\tilde{v}$, and (d) residual $w$ containing noise and small details of observed image $f$. Here, the reconstruction process of each component validates the decomposition of three terms $f=\tilde{h}+w=u\hspace{0.17em}+\hspace{0.17em}\tilde{v}\hspace{0.17em}+\hspace{0.17em}w$.

3.4. Fuzzy Representation of Edges and Textures

In this step of the process, the method is to produce a fuzzy representation (in [0,1]) of the features, contours, and principally of oriented texture (nonintrinsic) of the support image 𝑣₂ generated in the previous step. For this purpose, we first apply a morphological filter to simultaneously remove background and heterogeneous regions of 𝑣₂, seeking to highlight the oscillatory characteristics of that image. Next, the algorithm normalizes the image; that is, it translates the preprocessed image to the interval [0,1]. Here, the main idea was to use a morphological filter as is presented here in after.

3.4.1. Mathematical Morphology

For the treatment of the image 𝑣₂, we used an approach based on morphological filters, which has been showing to be a powerful tool for analyzing the structure of an image as well as for investigating the geometry of objects that constitute an image. As the objective here is to maintain oscillatory features of component 𝑣₂, a good alternative to doing this is to use a combination between the analyzed image and its version obtained by an opening or closing operator. Those combinations create a class of transformations called Top-hats.

According to [27], the open Top-hat transform WTH of an image 𝐼 is defined by WTH(𝐼)=𝐼−𝛾(𝐼), (3.12) while the close Top-hat transform BTH is given by BTH(𝐼)=𝜌(𝐼)−𝐼, (3.13) where 𝛾 and 𝜌 are opening and closing operators, respectively. For details, see [27].

In this step of the algorithm, the objective is to emphasize the texture and simultaneously remove heterogeneous parts of the image. In such case, we opted to use a Top-hat transform. Precisely, we applied transformation (3.13) to the image 𝑣₂, obtained in the previous step, to highlight its oscillatory features at the same time as correcting its background. Resulting image will have oscillatory details characterized by shades of gray very close to the black while the interior of the objects and background will be characterized by shades of gray close to the white.

To finish this step, we convert the preprocessed image to the interval of shades of gray [0,1]. The output image obtained in this process is here denoted by 𝜔_𝑣2. The 𝜔_𝑣2 component is a fuzzy representation of contours, and warped texture, which serve as a guide to generate final component 𝑣 defined by oriented texture, intrinsic contours and irregular details of the initial image.

3.5. Output of Component Having Only Texture

This step is aimed to synthesize final component 𝑣 representing all oscillatory features of the image, excepting noise. This component must be composed of the oriented texture (nonintrinsic), intrinsic contours, edges, and irregular details.

To do this, the idea is to combine auxiliary components generated in steps 2 and 4, that is, the component defined by intrinsic contours and parts of texture 𝑣₁ and the component 𝜔_𝑣2, characterized by the fuzzy representation of oriented texture from 𝑣, respectively. Motivated by the ideas of [1, 2, 6], we introduced an efficient weighted technique between 𝜔_𝑣2 and 𝑣₁, so that absent characteristics of each component could be compensated by other components. More precisely, the texture contained in 𝜔_𝑣2 will be superimposed on the regions where there is an absence of this information in support image 𝑣₁. Therefore, the second component of our decomposition in three terms 𝑣 is given as follows: 𝑣(𝑥)∶=𝜔_𝑣2(𝑥)⋅𝑣₁(𝑥),𝑥∈Ω, (3.14) where the above product is computed pixel by pixel.

Here, the proposed idea is very similar to that used by the diffusivity term 𝑔 studied in step 2, which balances and attributes weights to each pixel according to its classification in the image.

This highlighting among pixels does not contribute for a variation in the range of the input image 𝑣₁, since 𝜔_𝑣2 belongs to initial range [0,1].

Because of having applied the closing top-hat transform (3.13), pixels that represent texture in 𝜔_𝑣2 will be closer to zero while those that represent background and homogeneous regions will be close to one. On the other hand, in component 𝑣₁ there are no pixels having high variations, since noise was previously eliminated. Moreover, the component 𝑣₁ preserves the edges and intrinsic contours of the observed image 𝑓, which does not happen with term 𝜔_𝑣2. Therefore, the main advantages of each component can be used: oriented texture of 𝜔_𝑣2 and intrinsic contours of 𝑣₁.

3.6. Output of Restored Image and Residual Component

The last step of the propose scheme consists of obtaining the recovered image ℎ and the image composed of noise and background 𝑤. Moreover, this step finishes the decomposition process of initial image 𝑓 in terms of the three components previously described: 𝑢 (structure), 𝑣 (texture), and 𝑤 (noise and background).

As component 𝑢 represents the structure/cartoon of the observed image 𝑓 and 𝑣 the warped texture, edges, and intrinsic contours (but not noise), according to the classic Meyer model [9], it is sufficient to add both to obtain the restored image. The great advantage is that noise was removed by applying the previous steps, as much as for 𝑢 as for 𝑣. Then, in this case, ℎ(𝑥)∶=𝑢(𝑥)+𝑣(𝑥),𝑥∈Ω, (3.15) where the sum is done pixel by pixel.

The characterizing noise is done by calculating the residue between the restored image ℎ and observed image 𝑓, that is, 𝑤(𝑥)∶=𝑓(𝑥)−ℎ(𝑥),𝑥∈Ω, (3.16) where the sum is calculated term by term. In this case, this operation defines not only noise added to 𝑓 but also small fragments in the background of the image.

Finally, besides generating reconstructed image ℎ, the algorithm also satisfies the decomposition of three terms 𝑓=𝑢+𝑣+𝑤 or of an classical decomposition of two components 𝑓=ℎ+𝑤 (see Figure 8).

4. Experimental Results

Now we present some experiments obtained by our scheme, where images in a grayscale defined in the standard interval [0,255] were used. All the tested images are constituted by matrices of dimension 256×256. In the case of synthetic images, it is possible calculating the standard deviation of the noise 𝜎_noi. However, for real images, 𝜎_noi must be estimated in some other way. Thus, in step 2 we choose 𝜎_noi based on the visual quality of 𝑣. On the other hand, 𝜎_noi can be provided by the user in the step 1 since the goal here is smoothing the image 𝑓.

In order to validity our approach with respect to the tested methods, we used the statistical measure PSNR (peak signal-to-noise ratio), which is measured in dB.

In the first step of the algorithm, we use (3.1) supported by (3.3) and (3.5), where (motivated by [25]) we adopt 𝜆=1 and the temporal step Δ𝑡=0.1 in all examples considered. The number of iterations 𝑁 and the standard noise deviation 𝜎_noi must be given because they vary with each experiment. In the second step we again use the equation described in the previous step and, once again, we adopted 𝜆=1. Here, the number of iterations 𝑁 is determined by (3.7), while Δ𝑡 and 𝜎_noi (for real images) continue being input parameters. In the third step, the only input parameter is the threshold value 𝛾 while in the fourth step we adopt two types of structuring elements in the top-hat transform: disk or ball. In the fifth and sixth steps there are no parameters to be determined or provided. In the case of images in Figure 3, 𝑁=50, 𝜎_noi=25 while in Figure 4, Δ𝑡=0.1. In Figure 5, we take 𝛾=0.23 and in Figure 6 a disk with radius 4.

In Section 4.1 we emphasize the decomposition technique in three terms previously mentioned. In Section 4.2 we evaluate the good performance of the proposed scheme in comparison to other recent models in literature.

4.1. Restoring and Decomposition Using the Proposed Scheme

In the following, we show two experiments done on images having different levels of complexity: a highly-detailed real image and one of fingerprint.

Our first experiment mentions the real image of Barbara. Here the image contaminated with noise (SNR=9.1) contains important features to be preserved such as textures on the legs, in the region near the neck, in the background, and intrinsic contours of the face. Figure 9(a) presents a noisy image 𝑓. The image in Figure 9(c) shows component 𝑢 characterized by the structure/cartoon of 𝑓 obtained in the first step of our algorithm (𝜎_noi=15 and 𝑁=30). Taking Δ𝑡=0.15, 𝜎_noi=5 in the second step, adopting 𝛾=0.06 for the third step, and choosing a ball with radius 12 and height 10 in the fourth step, the algorithm generates Figure 9(d), 𝑣 containing restored oscillatory details, excepting noise. Both intrinsic contours and oriented texture of 𝑓 are present in 𝑣; that is, there was no significant loss of any type of oscillatory features. We can see that Figure 9(c) as in Figure 9(d), the images remain well defined in the visual perception sense of classical decomposition models such as [9, 13]. In Figure 9(b) we present the recovered image ℎ and in Figure 9(e) the component 𝑤 containing the noise and the background. We can note that our algorithm works efficiently with noise removal, preserving texture and contours, besides producing, from a visual point of view, a well-defined three-term composition.

(a)

(b)

(c)

(d)

(e)

Figure 9 

Decomposition into three components. (a) Observed image $f$, (b) image recovered using proposed method, (c) structure/cartoon $u$, (d) texture and intrinsic details $\tilde{v}$, and (e) noise and residual parts $w$. Here, $f=u+\tilde{v}+w$.

In the second experiment we take a synthetic image of fingerprint having a considerable noise level (SNR=4.5). Figure 10(a) represents a version with noise. Figures 10(b), 10(c), and 10(d) denote the three components of the evaluated decomposition: structure/cartoon 𝑢, texture 𝑣, and noise background 𝑤, respectively, while Figures 11(a) and 11(b) show the original image and the restored image obtained with our scheme (step 1—𝜎_noi=25, 𝑁=200; step 2—Δ𝑡=0.06; step 3—𝛾=0.16; step 4—B = ball with radius 10 and height 3), respectively. As to visual quality, the restored image Figure 11(b) is fairly close to the original image Figure 11(a). Furthermore, the residual image 𝑤=𝑓−(𝑢+𝑣) represented by Figure 10(d) detailed only noise and parts of the background of the observed image 𝑓, not maintaining any type of texture traces of the fingerprint.

(a)

(b)

(c)

(d)

Figure 10 

Decomposition into three components. (a) Noisy image $f$, (b) structure/cartoon $u$, (c) texture and intrinsic details $\tilde{v}$, and (d) noise and background $w$. All three components satisfy $f=u+\tilde{v}+w$.

(a)

(b)

Figure 11 

Restoring the image $f$. (a) Original image $h$ (without noise) and (b) recovered image $\tilde{h}$.

4.2. Comparison to Some Existing Methods

To attest to the good performance of the proposed method, we compared it to recent models in literature. Parameters adopted in each of the models tested were chosen according to the best visual quality obtained from each one of those models, in addition to computation of PSNR between the original image and the compared image. Classical models that remove noise but do not cover treating texture were not considered.

Figure 12(a) shows a noisy part (SNR=9.97, PSNR=18.58) from Barbara's image. Figures 12(b) and 12(c) were obtained using curvelet transform [29, 30] (PSNR=19.04) and wave atoms transform [18, 19] (PSNR=20.89), respectively. Here, both techniques were supported by the hard threshold. In the image using curvelet, texture was not appropriately recovered. Moreover, pseudo-Gibbs phenomena was present. In contrast, the wave atom transform correctly restored texture but produced a blurred image. Figure 12(d) shows the image restored by the model based on adaptive fidelity term [22] (with 𝐶=0, PSNR=20.95) while Figure 12(e) shows the version recovered by the model described in [20] (0.0005 step size, 11 iterations and PSNR=21.06), which combines nonlinear anisotropic diffusion with curvelet shrinkage. Although the model based on adaptive fidelity term has recovered texture, some important image details were excessively smoothed as with face and hand. Furthermore, a diffusion-based curvelet shrinkage did not produce any type of intensive smoothing but retains part of the noise in the reconstructed image, besides the PSNR being higher. Figure 12(f) (PSNR=22.12) is the image restored by the proposed scheme using the same parameters mentioned at the beginning of this section. In this case, both texture and image details (intrinsic contours) are recovered, besides achieving a satisfactory minimizing of the noise level without excessive smoothing. Furthermore, the PSNR from the proposed method presented the biggest value in comparison with considering techniques. Figures 13(a), 13(b), 13(c), 13(d), and 13(e) show residues in relation to Figure 12(a) obtained by Figures 12(b), 12(c), 12(d), 12(e), and 12(f), respectively. As mentioned before, we can see that the curvelet transform (Figure 13(a)) removed part of the texture. Contrasting this, the wave atoms transform (Figure 13(b)) maintained the texture but extracted some features from the image. In residue Figure 13(c) generated by the model [22] there is no texture; nevertheless, it can be seen that all intrinsic contours of the image stay annexed to the residue. This does not happen in component Figure 13(d) extracted by model [20], since part of the texture and details of the image were retained in the same proportions of the residue, but, a considerable part of the texture was removed in the process. Finally, component Figure 12(f) extracted by our algorithm shows that only noise and few relevant details were left annexed to the residue, therefore attesting to the efficiency of the proposed method.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 12 

Comparison to models existing in literature. (a) Noisy image, (b) denoising by curvelets, (c) by wave atoms, (d) by adaptive fidelity term, (e) by nonlinear diffusion combined with curvelet, and (f) by proposed scheme.

(a)

(b)

(c)

(d)

(e)

Figure 13 

Components removed using each of the methods. (a) Residual by curvelets, (b) by wave atoms, (c) by adaptive fidelity term, (d) by nonlinear diffusion combined with curvelet, and (e) by proposed scheme.

5. Conclusion

In this work we gather important mathematical techniques for image processing and we combine these techniques to generate an efficient algorithm for decomposition and noise removal in image processing. The new method has as its aim treating images contaminated with noise, having a high texture concentration, intrinsic contours, and irregular patterns. The scheme combines elementary techniques, such as classic morphological operators, with more sophisticated harmonic analysis models, such as wave atoms. Moreover, the scheme has, in some of its processing steps, a best parameter automatic selector. Based on the proposed method, we propose an efficient decomposition standard to separate the observed image into three well-defined components, as was shown previously. One of the advantages of this type of decomposition is that there is possible, from the degraded image, the individual treatment of each component, which allows a range of image processing applications such as image segmentation and digital inpainting. Experimental tests show the efficiency of the new method, even when compared to recent harmonic analysis techniques and to models based on nonlinear diffusion for the processing of images with texture.

Acknowledgments

The authors thank the São Paulo State Research Foundation (FAPESP) and the Brazilian Commission for Higher Education (CAPES) for financial support. They also thank Alagacone Sri Ranga and the unknown referees for suggestions on the improvement of the paper.