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 [18]. 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 [1016]. 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 [1719]. New hybrid models have been proposed with a view to integrate the main advantages of different methods that minimize noise, such as [2024]. 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, where represents the Gaussian type noise with mean 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 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 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.

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, 58]. Another good alternative is to use a simultaneous cartoon-texture decomposition model, which can be found in [914, 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: 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 , when and . 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]: with where 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, 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 . Then, 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 and , 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, 68]) 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 of the initial image , which resulted in an estimate for optimal time to stop the evolutionary process (3.1), given by 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 . For synthetic images, it is possible to calculate , but for real images it is usually not. In the latter case we try to estimate 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.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 , 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 , , , where and are two positive constants. According to [18], the elements of a frame of wave packets are called wave atoms when

To construct wave atoms for our problem, we first considered the case of a family of wave packets , , , in the scope of above same conditions. Let be a continuous real function given in , such that for , , and

Considering , with , the Fourier transform of is given by where , . In this case, must be such that .

Thus, we can write functions that make up the base as , whose coefficients can be obtained by . Here, represents the observed signal.

According [18], the extension to two-dimensional case ( version) can be computed by , where and a second equation is based on Hilbert transform relative to these wavelet packets. Therefore, the combination and 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 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: where denotes the wave atom transform (we use the version ), 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 . Next, we remove some insignificant coefficients from the wave atoms through thresholding (3.10). Finally, we apply the inverse transform 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.

3.4. Fuzzy Representation of Edges and Textures

In this step of the process, the method is to produce a fuzzy representation (in ) 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 . 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 of an image is defined by while the close Top-hat transform is given by 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 . The output image obtained in this process is here denoted by . The 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 and , that is, the component defined by intrinsic contours and parts of texture and the component , 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 and , so that absent characteristics of each component could be compensated by other components. More precisely, the texture contained in 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: 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 , 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 belongs to initial range .

Because of having applied the closing top-hat transform (3.13), pixels that represent texture in 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 . Therefore, the main advantages of each component can be used: oriented texture of 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, 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, 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 were used. All the tested images are constituted by matrices of dimension . In the case of synthetic images, it is possible calculating the standard deviation of the noise . However, for real images, must be estimated in some other way. Thus, in step 2 we choose based on the visual quality of . On the other hand, 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 and the temporal step in all examples considered. The number of iterations and the standard noise deviation 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 . Here, the number of iterations is determined by (3.7), while and (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, , while in Figure 4, . In Figure 5, we take and in Figure 6 a disk with radius .

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 () 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 ( and ). Taking , in the second step, adopting for the third step, and choosing a ball with radius and height 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.

In the second experiment we take a synthetic image of fingerprint having a considerable noise level (). 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—, ; step 2—; step 3—; step 4—B = ball with radius and height ), 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.

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 (, ) from Barbara's image. Figures 12(b) and 12(c) were obtained using curvelet transform [29, 30] () and wave atoms transform [18, 19] (), 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 , ) while Figure 12(e) shows the version recovered by the model described in [20] ( step size, iterations and ), 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) () 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.

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.


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.