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.

(ProQuest: ... denotes non-US-ASCII text omitted.)

Recommended by Panos Liatsis

DCCE/IBILCE, UNESP-São Paulo State University, Rua Cristóvão Colombo 2265, 15054-000 São José do Rio Preto, SP, Brazil

Received 29 June 2009; Revised 18 December 2009; Accepted 6 February 2010

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 f is the sum of two components, f=u+v , so that u contains the structure and the objects of the image (cartoon) while v contains the oscillatory characteristics; that is, texture and noise.

In most cases, when an image has only few regions defined by texture, the term v is discarded, being characterized by the equation v=f-u ; that is, the models in this category assume that component v 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 u and oscillatory v , so that f=u+v . 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 f has noise and a high texture concentration, then the term u 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 v . In all those models, texture and noise are treated equally, which makes it difficult to identify the oscillatory component v . For example, for real images, it is practically impossible to obtain satisfactory results using some f=u+v 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 f into three terms, f=u+v...+w , with u representing the structure or cartoon of the image, v... only the texture and intrinsic contours and w 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 f be the observed image with noise and h the original image (without noise), represented by the functions f:Ω[arrow right]... and h:Ω[arrow right]... , where Ω is a rectangular region of ^...2 . Here, we assumed that both f and h can be periodically extended to the ^...2 , so as to have f,h∈^L2 (^...2 ) . We supposed that the noise is additive, that is, [figure omitted; refer to PDF] where n represents the Gaussian type noise with mean 0 e variance ^σn2 .

Furthermore, we assumed that the original image h 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 f ; that is, the impact of noise n(x) should be minimum in the output image, making it visibly closer to the original image h . The texture and intrinsic contours built-in the image must be maintained and highlighted so that output images u , v... , and w satisfy [figure omitted; refer to PDF] and represent the following features of f :

(i) u : structure, skeleton, or cartoon,

(ii) v... : texture, intrinsic contours, and irregular details,

(iii): w : noise and background.

2.2. The Proposed Scheme

Let f 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 f is calculated into two components, u and v , so that f=u+v , where u contains the skeleton or cartoon of f and v 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 v , thereby producing an auxiliary component _v1 that contains intrinsic contours, parts of the texture, and parts of the edges of v .

(3) Oriented texture support component . An auxiliary image _v2 is produced containing the warped oscillatory patterns of v 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 _v2 to obtain a fuzzy edge-texture representation _ωv2 :Ω[arrow right][0,1] of the image _v2 .

(5) Output of component having only texture . A combination of component _v1 , obtained in step 2, with the component _v2 , 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 v... .

(6) Output of restored image and residual component . A composition of the component u , obtained in 1, is done with the component v... , generated in the previous step. In addition, the image w 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.

[figure omitted; refer to PDF]

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

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

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

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

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 f into two components u and v , with f=u+v , 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 f by applying PDE [2], obtaining as a result the component cartoon u . Then, we do a simple operation f-u to determine the component having texture and noise, that is, v:=f-u . As the objective is to obtain u and v 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 v (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: [figure omitted; refer to PDF] [figure omitted; refer to PDF] where v represents the initial image having texture and noise, ^v(t) is its version on scale t , g=g(|∇_Gσ *^v(t) |) is a nonnegative and nonincreasing function called diffusivity term, _Gσ *^v(t) determines the convolution of signal ^v(t) with the Gaussian function _Gσ , and λ is a weighting parameter. Here, the symbol |·| denotes the Euclidean norm while the constant σ denotes the standard noise deviation of the image v .

Generally, the diffusivity term g , besides being nonincreasing and nonnegative, is such that g(0)=1 , g(s)[arrow right]0 when s[arrow right]∞ and g(|s|)∈[0,1] . In this paper we choose g based on the Perona and Malik [6] diffusivity term, together with the ideas of Alvarez et al. [1] and Barcelos et al. [25]: [figure omitted; refer to PDF] with [figure omitted; refer to PDF] where k=k(σ)≥0 is a σ -dependent constant and t... is a scalar variable related to the space of the scale produced by the Gaussian (3.4). The best choice for the scale t... will be better described in following pages.

With the intent to automate the computation of (3.3) and to avoid antiquated choices for k , we adopt it according the ideas in [25], that is, [figure omitted; refer to PDF] where _a1 , _a2 , _b1 , and _b2 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), k works as an edge selection parameter: for a fixed image with a high k value, false edges may be identified while with a small k value, only prominent edges will be selected.

The equation (3.1) can be seen as a balancing between smoothing and "keeping close to component v ". This balancing is managed by the diffusivity term g , 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 |∇_Gσ *^v(t) | small, which implies in g~1 . Then, (1-g)~0 and also the reaction term (^v(t) -v) 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 |∇_Gσ *^v(t) | is big we have g~0 and (1-g)~1 , implying that the reaction term (^v(t) -v) proposed by Nordström [8] will retain strongly the initial features of the image v 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 v 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 k and the best scale for t... , which can be automatically computed.

In [25], the authors linked the t... scale of the Gaussian kernel _Gσ with the noise standard deviation _σnoi of the initial image v=v(x^)(t=0) , which resulted in an estimate for optimal time T to stop the evolutionary process (3.1), given by [figure omitted; refer to PDF] where σ is the standard deviation of image v with noise and a=σ , with σ being the constant present in the Gaussian kernel (3.4).

Encouraged by the authors of [25], we take the optimal smoothing time t... as in (3.6), that is, t...=T . Having the optimal stopping time T , it is possible to automatically obtain (unless the temporal step Δt ) the number of iterations of the model and also the best scale t... 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 T 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 T (3.6), as presented in [2]. In this case, T is used in the Gaussian function (3.3), as was previously described. Also, the N number of iterations is calculated based on T , which is given by [figure omitted; refer to PDF]

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 _{[straight phi]μ} , with subscript μ=(j,m,n)=(j,_m1 ,_m2 ,_n1 ,_n2 ) . The five quantities above are integer values and index a point (_xμ ,_ωμ ) in a space-phase such as _xμ =^2-j n , _ωμ =π^2j m , _c1^2j ≤_{max i=1,2} |_mi |≤_c2^2j , where _c1 and _c2 are two positive constants. According to [18], the elements of a frame of wave packets _{[straight phi]μ} are called wave atoms when [figure omitted; refer to PDF]

To construct wave atoms for our problem, we first considered the case of a family 1D of wave packets ^ψm,nj (x)≥0 , m≥0 , n∈... , in the scope of above same conditions. Let [varphi] be a continuous real function given in [-7π/6, 5π/6] , such that for |ω|≤π/3 , [varphi]^((π/2)-ω)2 +[varphi]^((π/2)+ω)2 =1 , and [varphi](-(π/2)-2ω)=[varphi]((π/2)+ω).

Considering ^ψm0 (t)=2Re{exp (iπ(m+1/2)t)ν((-1⁾ⁿ (t-1/2))} , with ν(t)=(1/2π)∫[varphi](ω)exp (iωt)dω , the Fourier transform of ^ψm0 is given by [figure omitted; refer to PDF] where _...m =(-1^)m , _αm =(π/2)(m+(1/2)) . In this case, [varphi] must be such that ^∑m=0∞ |^ψ...m0|2 =1 .

Thus, we can write functions that make up the base as ^ψm,nj (x)=^ψmj (x-^2-j n)=^2j/2ψm0 (^2j x-n) , whose coefficients can be obtained by _cj,m,n =∫^ψm,nj v(x)dx=(1/2π)∫exp (i^2-j -nω)^ψ...mj (ω)¯v...(ω)dω . Here, v represents the observed signal.

According [18], the extension to two-dimensional case (2D version) can be computed by ^{[straight phi]μ+} (_x1 ,_x2 ):=^{ψ_m1 ,_n1 j} (_x1 )^{ψ_m2 ,_n2 j} (_x2 ),^{[straight phi]μ-} (_x1 ,_x2 ):=H^{ψ_m1 ,_n1 j} (_x1 )H^{ψ_m2 ,_n2 j} (_x2 ) , where μ=(j,m,n) and a second equation is based on Hilbert transform relative to these wavelet packets. Therefore, the combination ^{[straight phi]μ(1)} =(^{[straight phi]μ+} +^{[straight phi]μ-} )/2 and ^{[straight phi]μ(2)} =(^{[straight phi]μ+} -^{[straight phi]μ-} )/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 _vc =_∑μ θ(^cμ(1) (v))^{[straight phi]μ(1)} +θ(^cμ(2) (v))^{[straight phi]μ(2)} , where θ=_θγ is a thresholding function, here adopted by [figure omitted; refer to PDF] where γ begin the threshold value.

In this work we use wave atom shrinkage to extract oriented texture from image component v . We take a transform based on wave atoms as follows: [figure omitted; refer to PDF] where WA denotes the wave atom transform (we use the version 2D ), (WA^)-1 the inverse transform (see [18, 19]), _θγ is the threshold function (3.10) mentioned earlier, and _v2 is the output component. Here, the nonlinear operator T... will keep and highlight important characteristics of the examined image v 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-1 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 ^C2 . 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.

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

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

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

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Details of step 4. (a) Input image _v2 and (b) highlight of oriented texture of _v2 with background correction (in the scale [0,1] ), here denoted by _ωv2 .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Details of step 5. (a) Input image _v1 , which is characterized by intrinsic contours, (b) fuzzy representation _ωv2 of warped texture, and (c) combination of features of _v1 with _ωv2 , resulting in the final component v... defined by oscillatory patterns of observed image f , except noise.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

Details of final step 6. (a) Input image u , (b) texture v... , (c) recovered image h... , obtained by the addition of u to 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=h...+w=u + v... + w .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

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 _v2 generated in the previous step. For this purpose, we first apply a morphological filter to simultaneously remove background and heterogeneous regions of _v2 , 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 _v2 , 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 _v2 , 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 I is defined by [figure omitted; refer to PDF] while the close Top-hat transform BTH is given by [figure omitted; refer to PDF] 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 _v2 , 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 _ωv2 . The _ωv2 component is a fuzzy representation of contours, and warped texture, which serve as a guide to generate final component v... 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 v... 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 _v1 and the component _ωv2 , characterized by the fuzzy representation of oriented texture from v , respectively. Motivated by the ideas of [1, 2, 6], we introduced an efficient weighted technique between _ωv2 and _v1 , so that absent characteristics of each component could be compensated by other components. More precisely, the texture contained in _ωv2 will be superimposed on the regions where there is an absence of this information in support image _v1 . Therefore, the second component of our decomposition in three terms v... is given as follows: [figure omitted; refer to PDF] where the above product is computed pixel by pixel.

Here, the proposed idea is very similar to that used by the diffusivity term g 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 _v1 , since _ωv2 belongs to initial range [0,1] .

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

3.6. Output of Restored Image and Residual Component

The last step of the propose scheme consists of obtaining the recovered image h... and the image composed of noise and background w . Moreover, this step finishes the decomposition process of initial image f in terms of the three components previously described: u (structure), v... (texture), and w (noise and background).

As component u represents the structure/cartoon of the observed image f and v... 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 u as for v... . Then, in this case, [figure omitted; refer to PDF] where the sum is done pixel by pixel.

The characterizing noise is done by calculating the residue between the restored image h... and observed image f , that is, [figure omitted; refer to PDF] where the sum is calculated term by term. In this case, this operation defines not only noise added to f but also small fragments in the background of the image.

Finally, besides generating reconstructed image h... , the algorithm also satisfies the decomposition of three terms f=u + v... + w or of an classical decomposition of two components f=h...+w (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 v . On the other hand, _σnoi can be provided by the user in the step 1 since the goal here is smoothing the image f .

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 Δt=0.1 in all examples considered. The number of iterations N 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 N is determined by (3.7), while Δt 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, N=50 , _σnoi =25 while in Figure 4, Δt=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 f . The image in Figure 9(c) shows component u characterized by the structure/cartoon of f obtained in the first step of our algorithm (_σnoi =15 and N=30 ). Taking Δt=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), v... containing restored oscillatory details, excepting noise. Both intrinsic contours and oriented texture of f are present in v... ; 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 h... and in Figure 9(e) the component w 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.

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

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

(e) [figure omitted; refer to PDF]

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 u , texture v... , and noise background w , respectively, while Figures 11(a) and 11(b) show the original image and the restored image obtained with our scheme (step 1--_σnoi =25 , N=200 ; step 2--Δt=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 w=f-(u+v...) represented by Figure 10(d) detailed only noise and parts of the background of the observed image f , not maintaining any type of texture traces of the fingerprint.

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

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

Restoring the image f . (a) Original image h (without noise) and (b) recovered image h... .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

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 C=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.

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) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

(e) [figure omitted; refer to PDF]

(f) [figure omitted; refer to PDF]

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.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

(e) [figure omitted; refer to PDF]

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.

References

[1] L. Alvarez, P. L. Lions, J. M. Morel, "Image selective smoothing and edge detection by nonlinear diffusion. II," SIAM Journal on Numerical Analysis , vol. 29, no. 3, pp. 845-866, 1992.

[2] C. A. Z. Barcelos, M. Boaventura, E. Silva Jr., "A well-balanced flow equation for noise remove and edge detection," IEEE Transactions on Image Processing , vol. 12, pp. 751-763, 2003.

[3] M. Burger, G. Gilboa, S. Osher, J. Xu, "Nonlinear inverse scale space methods," Communications in Mathematical Sciences , vol. 4, no. 1, pp. 179-212, 2006.

[4] M. Burger, S. Osher, J. Xu, G. Gilboa, "Nonlinear inverse scale space methods for image restoration," Lecture Notes in Computer Science , vol. 3752, pp. 85-96, 2005.

[5] T. F. Chan, S. Osher, J. Shen, "The digital TV filter and nonlinear denoising," IEEE Transactions on Image Processing , vol. 10, no. 2, pp. 231-241, 2001.

[6] P. Perona, J. Malik, "Scale-space and edge detection using anisotropic diffusion," IEEE Transactions on Pattern Analysis and Machine Intelligence , vol. 12, no. 7, pp. 629-639, 1990.

[7] L. I. Rudin, S. Osher, E. Fatemi, "Nonlinear total variation based noise removal algorithms," Physica D , vol. 60, no. 1-4, pp. 259-268, 1992.

[8] K. N. Nordström, "Biased anisotropic diffusion: a unified regularization and diffusion approach to edge detection," Image and Vision Computing , vol. 8, no. 4, pp. 318-327, 1990.

[9] Y. Meyer Oscillating Patterns in Image Processing and Nonlinear Evolution Equations , vol. 22, of University Lecture Series, American Mathematical Society, Providence, RI, USA, 2002.

[10] J. B. Garnett, T. M. Le, Y. Meyer, L. A. Vese, "Image decompositions using bounded variation and generalized homogeneous Besov spaces," Applied and Computational Harmonic Analysis , vol. 23, no. 1, pp. 25-56, 2007.

[11] S. E. Levine, M. Ramsey, T. Misner, S. Schwab, "An adaptive variational model for image decomposition," in Proceedings of 5th International Workshop on Energy Minimization Methods in Computer Vision and Pattern Recognition (EMMCVPR '05), vol. 3757, of Lecture Notes in Computer Science, pp. 382-397, St. Augustine, Fla, USA, November 2005.

[12] L. Lieu Contribution to problems in image restoration, decomposition, and segmentation by variational methods and partial differential equations , Ph.D. thesis, UCLA, Los Angeles, Calif, USA, 2006.

[13] L. A. Vese, S. J. Osher, "Modeling textures with total variation minimization and oscillating patterns in image processing," Journal of Scientific Computing , vol. 19, no. 1-3, pp. 553-572, 2003.

[14] W. Yin, D. Goldfarb, S. Osher, "Image cartoon-texture decomposition and feature selection using the total variation regularized ^L1 functional," Variational, Geometric, and Level Set Methods in Computer Vision , vol. 3752, of Lecture Notes in Computer Science, pp. 73-84, Springer, New York, NY, USA, 2005.

[15] T. F. Chan, S. Esedoglu, F. E. Park, "Image decomposition combining staircase reduction and texture extraction," Journal of Visual Communication and Image Representation , vol. 18, no. 6, pp. 464-486, 2007.

[16] J. F. Aujol, A. Chambolle, "Dual norms and image decomposition models," International Journal of Computer Vision , vol. 63, no. 1, pp. 85-104, 2005., [email protected]; [email protected]

[17] J. L. Starck, E. J. Candès, D. L. Donoho, "The curvelet transform for image denoising," IEEE Transactions on Image Processing , vol. 11, no. 6, pp. 670-684, 2002.

[18] L. Demanet, L. Ying, "Wave atoms and sparsity of oscillatory patterns," Applied and Computational Harmonic Analysis , vol. 23, no. 3, pp. 368-387, 2007.

[19] L. Demanet, L. Ying, "Curvelets and wave atoms for mirror-extended images," in Wavelets XII, vol. 6701, of Proceedings of SPIE, San Diego, Calif, USA, August 2007.

[20] J. Ma, G. Plonka, "Combined curvelet shrinkage and nonlinear anisotropic diffusion," IEEE Transactions on Image Processing , vol. 16, no. 9, pp. 2198-2206, 2007.

[21] W. C. O. Casaca, M. Boaventura, "A regularized nonlinear diffusion approach for texture image denoising," in Proceedings of the 22nd Brazilian Symposium on Computer Graphics and Image Processing (SIBGRAPI '09), IEEE Computer Society, Rio de Janeiro, Brazil, October 2009.

[22] G. Gilboa, N. Sochen, Y. Y. Zeevi, "Variational denoising of partly textured images by spatially varying constraints," IEEE Transactions on Image Processing , vol. 15, no. 8, pp. 2281-2289, 2006.

[23] J. Ma, "Image assimilation by geometric wavelet based reaction-diffusion equation," in Wavelets XII, vol. 6701, of Proceedings of SPIE, San Diego, Calif, USA, [email protected], August 2007.

[24] G. Plonka, J. Ma, "Nonlinear regularized reaction-diffusion filters for denoising of images with textures," IEEE Transactions on Image Processing , vol. 17, no. 8, pp. 1283-1294, 2008.

[25] C. A. Z. Barcelos, M. Boaventura, E. Silva Jr., "Edge detection and noise removal by use of a partial differential equation with automatic selection of parameters," Computational & Applied Mathematics , vol. 24, no. 1, pp. 131-150, 2005.

[26] L. Demanet Curvelets, wave atoms and wave equations , Ph.D. thesis, California Institute of Technology, Pasadena, Calif, USA, 2006.

[27] P. Soille Morphological Image Analysis. Principles and Applications , pp. xvi+391, Springer, Berlin, Germany, 2003., 2nd.

[28] L. F. Villemoes, "Wavelet packets with uniform time-frequency localization," Comptes Rendus Mathématique , vol. 335, no. 10, pp. 793-796, 2002.

[29] E. Candés, L. Demanet, D. Donoho, L. Ying, "Fast discrete curvelet transforms," Multiscale Modeling & Simulation , vol. 5, no. 3, pp. 861-899, 2006.

[30] E. J. Candès, D. L. Donoho, "New tight frames of curvelets and optimal representations of objects with piecewise ^C2 singularities," Communications on Pure and Applied Mathematics , vol. 57, no. 2, pp. 219-266, 2004.